Allow LocalDensity objects to hand out adouble-versions of themselves

Currently the user code that plugs together the assembler for a
particular problems has to know that ADOL-C is used internally
to compute certain derivatives.  In shows because certain objects
need to be instantiated with adouble as the number type.  One such
example is all objects derived from LocalDensity.

However, if this is to be avoided, then assemblers need to have
a way to turn a (e.g.) 'double' version of a density into an
'adouble' one.  This commit provides this way, by adding a
'makeActiveDensity' method.

dune/gfe/densities/bulkcosseratdensity.hh

@@ -25,6 +25,9 @@ namespace Dune::GFE {
     static const int gridDim = Position::size();
     static const int embeddedDim = Rotation<field_type,gridDim>::embeddedDim;
 
+    // The target space with 'adouble' as the number type
+    using ATargetSpace = GFE::ProductManifold<RealTuple<adouble,3>,Rotation<adouble,3> >;
+
     static FieldMatrix<field_type,gridDim,gridDim> curl(const Tensor3<field_type,gridDim,gridDim,gridDim>& DR)
     {
       FieldMatrix<field_type,gridDim,gridDim> result;
@@ -63,6 +66,12 @@ namespace Dune::GFE {
       b3_ = parameters.template get<double>("b3");
     }
 
+    /** \brief Constructor with a set of material parameters
+     */
+    BulkCosseratDensity(double mu, double lambda, double mu_c, double L_c, double q, const std::array<double,3>& b)
+      : mu_(mu), lambda_(lambda), mu_c_(mu_c), L_c_(L_c), q_(q), b1_(b[0]), b2_(b[1]), b3_(b[2])
+    {}
+
     /** \brief The energy \f$ W_{mp}(\overline{U}) \f$, as written in
      * the first equation of (4.4) in Neff's paper from 2006: A geometrically exact planar Cosserat shell model with microstructure: Existence of minimizers for zero Cosserat couple modulus
      * OR: the equation (2.27) of 2019: Refined dimensional reduction for isotropic elastic Cosserat shells with initial curvature
@@ -165,6 +174,12 @@ namespace Dune::GFE {
       return strainEnergyDensity;
     }
 
+    // Construct a copy of this density but using 'adouble' as the number type
+    virtual std::unique_ptr<LocalDensity<Position,ATargetSpace> > makeActiveDensity() const
+    {
+      return std::make_unique<BulkCosseratDensity<Position,adouble> >(mu_, lambda_, mu_c_, L_c_, q_, std::array<double,3>{b1_, b2_, b3_});
+    }
+
     /** \brief The density depends on the microrotation value, but not on the deformation
      */
     virtual bool dependsOnValue(int factor=-1) const override

dune/gfe/densities/chiralskyrmiondensity.hh

 #ifndef DUNE_GFE_DENSITIES_CHIRALSKYRMIONDENSITY_HH
 #define DUNE_GFE_DENSITIES_CHIRALSKYRMIONDENSITY_HH
 
+#include <adolc/adouble.h>
+
 #include <dune/common/fmatrix.hh>
 #include <dune/common/parametertree.hh>
 
@@ -22,6 +24,8 @@ namespace Dune::GFE
     using TargetSpace = UnitVector<field_type,3>;
     using Derivative = FieldMatrix<field_type, TargetSpace::embeddedDim, Position::size()>;
 
+    using ATargetSpace = typename TargetSpace::template rebind<adouble>::other;
+
   public:
 
     ChiralSkyrmionDensity(const Dune::ParameterTree& parameters)
@@ -30,6 +34,10 @@ namespace Dune::GFE
       kappa_ = parameters.template get<double>("kappa");
     }
 
+    ChiralSkyrmionDensity(double h, double kappa)
+      : h_(h), kappa_(kappa)
+    {}
+
     //! Dimension of a tangent space
     constexpr static int blocksize = TargetSpace::TangentVector::dimension;
 
@@ -68,6 +76,12 @@ namespace Dune::GFE
       return density;
     }
 
+    // Construct a copy of this density but using 'adouble' as the number type
+    virtual std::unique_ptr<LocalDensity<Position,ATargetSpace> > makeActiveDensity() const
+    {
+      return std::make_unique<ChiralSkyrmionDensity<Position,adouble> >(h_, kappa_);
+    }
+
     /** \brief The density depends on the value */
     virtual bool dependsOnValue([[maybe_unused]] int factor=-1) const override
     {

dune/gfe/densities/duneelasticitydensity.hh

@@ -26,6 +26,8 @@ namespace Dune::GFE
     constexpr static auto dim = Position::size();
     constexpr static auto embeddedBlocksize = TargetSpace::EmbeddedTangentVector::dimension;
 
+    using ATargetSpace = typename TargetSpace::template rebind<adouble>::other;
+
   public:
 
     /** \brief Constructor with a Dune::Elasticity::LocalDensity
@@ -34,6 +36,12 @@ namespace Dune::GFE
       : elasticityDensity_(elasticityDensity)
     {}
 
+    /** \brief Constructor with a Dune::Elasticity::LocalDensity
+     */
+    DuneElasticityDensity(std::unique_ptr<Elasticity::LocalDensity<dim,field_type,ctype> >&& elasticityDensity)
+      : elasticityDensity_(std::shared_ptr<Elasticity::LocalDensity<dim,field_type,ctype> >(elasticityDensity.release()))
+    {}
+
     /** \brief Evaluate the density
      *
      * \param x The current position
@@ -56,13 +64,23 @@ namespace Dune::GFE
       return (*elasticityDensity_)(x, factorDerivative);
     }
 
-    /** \brief The density depends on the value */
+    // Construct a copy of this density but using 'adouble' as the number type
+    virtual std::unique_ptr<LocalDensity<Position,ATargetSpace> > makeActiveDensity() const
+    {
+      // The active dune-elasticity density
+      auto activeDensity = elasticityDensity_->makeActiveDensity();
+
+      // Wrap it as a dune-gfe density
+      return std::make_unique<DuneElasticityDensity<Position,ATargetSpace,index> >(std::move(activeDensity));
+    }
+
+    /** \brief The density does not depend on the value */
     virtual bool dependsOnValue([[maybe_unused]] int factor=-1) const override
     {
       return false;
     }
 
-    /** \brief The density does not depend on the derivative */
+    /** \brief The density does depend on the derivative */
     virtual bool dependsOnDerivative([[maybe_unused]] int factor=-1) const override
     {
       return factor==-1 || factor==index;

dune/gfe/densities/harmonicdensity.hh

@@ -16,6 +16,8 @@ namespace Dune::GFE
     constexpr static auto dim = Position::size();
     constexpr static auto embeddedBlocksize = TargetSpace::EmbeddedTangentVector::dimension;
 
+    using ATargetSpace = typename TargetSpace::template rebind<adouble>::other;
+
   public:
 
     /** \brief Evaluate the density
@@ -31,6 +33,12 @@ namespace Dune::GFE
       return 0.5 * derivative.frobenius_norm2();
     }
 
+    // Construct a copy of this density but using 'adouble' as the number type
+    virtual std::unique_ptr<LocalDensity<Position,ATargetSpace> > makeActiveDensity() const
+    {
+      return std::make_unique<HarmonicDensity<Position,ATargetSpace> >();
+    }
+
     /** \brief The density does not depend on the value */
     virtual bool dependsOnValue([[maybe_unused]] int factor=-1) const override
     {

dune/gfe/densities/localdensity.hh

@@ -3,6 +3,9 @@
 
 #include <dune/common/fmatrix.hh>
 
+#include <adolc/adalloc.h>
+#include <adolc/adolc.h>
+
 namespace Dune::GFE {
 
   /** \brief A base class for energy densities to be evaluated in an integral energy
@@ -14,6 +17,7 @@ namespace Dune::GFE {
   class LocalDensity
   {
     using field_type = typename TargetSpace::field_type;
+    using ATargetSpace = typename TargetSpace::template rebind<adouble>::other;
     using DerivativeType = FieldMatrix<field_type,TargetSpace::EmbeddedTangentVector::dimension,Position::size()>;
 
   public:
@@ -28,6 +32,9 @@ namespace Dune::GFE {
                                    const typename TargetSpace::CoordinateType& value,
                                    const DerivativeType& derivative) const = 0;
 
+    // Construct a copy of this density but using 'adouble' as the number type
+    virtual std::unique_ptr<LocalDensity<Position,ATargetSpace> > makeActiveDensity() const = 0;
+
     /** \brief Whether the density depends on the 'value' parameter
      *
      * If TargetSpace is a ProductManifold, then this method returns the information

test/localintegralstiffnesstest.cc

@@ -101,7 +101,8 @@ int testHarmonicMapIntoSphere(TestSuite& test, const GridView& gridView)
   std::shared_ptr<LocalGeodesicFEStiffness<FEBasis,TargetSpace> > localIntegralStiffness;
 
   using LocalCoordinate = typename GridView::template Codim<0>::Geometry::LocalCoordinate;
-  auto harmonicDensity = std::make_shared<GFE::HarmonicDensity<LocalCoordinate,ATargetSpace> >();
+  auto harmonicDensity = std::make_shared<GFE::HarmonicDensity<LocalCoordinate,TargetSpace> >();
+  auto harmonicDensityA = std::make_shared<GFE::HarmonicDensity<LocalCoordinate,ATargetSpace> >();
 
   if constexpr (interpolationType==Geodesic)
   {
@@ -117,7 +118,7 @@ int testHarmonicMapIntoSphere(TestSuite& test, const GridView& gridView)
         ATargetSpace>;
 
     // Assemble using the old assembler
-    auto energy = std::make_shared<GFE::LocalIntegralEnergy<FEBasis, ALocalInterpolationRule,ATargetSpace> >(harmonicDensity);
+    auto energy = std::make_shared<GFE::LocalIntegralEnergy<FEBasis, ALocalInterpolationRule,ATargetSpace> >(harmonicDensityA);
 
     auto localGFEADOLCStiffness = std::make_shared<LocalGeodesicFEADOLCStiffness<FEBasis,TargetSpace> >(energy);
     assemblerADOLC = std::make_shared<GeodesicFEAssembler<FEBasis,TargetSpace> >(feBasis, localGFEADOLCStiffness);
@@ -139,7 +140,7 @@ int testHarmonicMapIntoSphere(TestSuite& test, const GridView& gridView)
         ATargetSpace>;
 
     // Assemble using the old assembler
-    auto energy = std::make_shared<GFE::LocalIntegralEnergy<FEBasis, ALocalInterpolationRule,ATargetSpace> >(harmonicDensity);
+    auto energy = std::make_shared<GFE::LocalIntegralEnergy<FEBasis, ALocalInterpolationRule,ATargetSpace> >(harmonicDensityA);
 
     auto localGFEADOLCStiffness = std::make_shared<LocalGeodesicFEADOLCStiffness<FEBasis,TargetSpace> >(energy);
     assemblerADOLC = std::make_shared<GeodesicFEAssembler<FEBasis,TargetSpace> >(feBasis, localGFEADOLCStiffness);
@@ -333,7 +334,7 @@ int testCosseratBulkModel(TestSuite& test, const GridView& gridView)
     assemblerADOLC = std::make_shared<MixedGFEAssembler<CompositeBasis,RigidBodyMotion> >(compositeBasis, localGFEADOLCStiffness);
 
     // Assemble using the new assembler
-    localIntegralStiffness = std::make_shared<GFE::LocalIntegralStiffness<CompositeBasis,LocalInterpolationRule,RigidBodyMotion> >(aBulkCosseratDensity);
+    localIntegralStiffness = std::make_shared<GFE::LocalIntegralStiffness<CompositeBasis,LocalInterpolationRule,RigidBodyMotion> >(bulkCosseratDensity);
   }
   else
   {
@@ -357,7 +358,7 @@ int testCosseratBulkModel(TestSuite& test, const GridView& gridView)
     assemblerADOLC = std::make_shared<MixedGFEAssembler<CompositeBasis,RigidBodyMotion> >(compositeBasis, localGFEADOLCStiffness);
 
     // Assemble using the new assembler
-    localIntegralStiffness = std::make_shared<GFE::LocalIntegralStiffness<CompositeBasis,LocalInterpolationRule,RigidBodyMotion> >(aBulkCosseratDensity);
+    localIntegralStiffness = std::make_shared<GFE::LocalIntegralStiffness<CompositeBasis,LocalInterpolationRule,RigidBodyMotion> >(bulkCosseratDensity);
   }
 
   MixedGFEAssembler<CompositeBasis,RigidBodyMotion> mixedAssemblerSmart(compositeBasis,