diff --git a/dune/gfe/densities/bulkcosseratdensity.hh b/dune/gfe/densities/bulkcosseratdensity.hh
index 099b1583a0792cd1cfb388ff5a0c5ba9177d2610..824f46f638d7387f70488b9c2ae0ba48f9d259cf 100644
--- a/dune/gfe/densities/bulkcosseratdensity.hh
+++ b/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: Reļ¬ned 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
diff --git a/dune/gfe/densities/chiralskyrmiondensity.hh b/dune/gfe/densities/chiralskyrmiondensity.hh
index 32119a92bb7b06279e116625bb89f24e8ddeaac3..4af14f8e8dcae30852e37453b9a21a82e3ff85fe 100644
--- a/dune/gfe/densities/chiralskyrmiondensity.hh
+++ b/dune/gfe/densities/chiralskyrmiondensity.hh
@@ -1,6 +1,8 @@
#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
{
diff --git a/dune/gfe/densities/duneelasticitydensity.hh b/dune/gfe/densities/duneelasticitydensity.hh
index c62af5698961e002ec24b7c5509e183555662d66..b23b05b05ea517e36bf7c5dc6d1cff5bb91b71e9 100644
--- a/dune/gfe/densities/duneelasticitydensity.hh
+++ b/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;
diff --git a/dune/gfe/densities/harmonicdensity.hh b/dune/gfe/densities/harmonicdensity.hh
index 644743c54cd6272a8cd897021807f7e704cffc16..916ba6db38e63c1fec093c042ec9cb44dac4b048 100644
--- a/dune/gfe/densities/harmonicdensity.hh
+++ b/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
{
diff --git a/dune/gfe/densities/localdensity.hh b/dune/gfe/densities/localdensity.hh
index 4b9de66daf6c4bd3d560d6f91484250a5bbc04f5..00a5f1c2c98bf6787a981cdfb317325999530d4f 100644
--- a/dune/gfe/densities/localdensity.hh
+++ b/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
diff --git a/test/localintegralstiffnesstest.cc b/test/localintegralstiffnesstest.cc
index f6887d7bcabb7dad13fbf8538da9c46bdd748350..8c5c450453f0850b326d5da04c9cd5255aec9bfc 100644
--- a/test/localintegralstiffnesstest.cc
+++ b/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,