From fb50ca87a6036f76cb4b59847114aa0bdcacc4ad Mon Sep 17 00:00:00 2001
From: Oliver Sander <oliver.sander@tu-dresden.de>
Date: Tue, 4 Jun 2024 07:32:23 +0200
Subject: [PATCH] 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 | 15 ++++++++++++++
dune/gfe/densities/chiralskyrmiondensity.hh | 14 +++++++++++++
dune/gfe/densities/duneelasticitydensity.hh | 22 +++++++++++++++++++--
dune/gfe/densities/harmonicdensity.hh | 8 ++++++++
dune/gfe/densities/localdensity.hh | 7 +++++++
test/localintegralstiffnesstest.cc | 11 ++++++-----
6 files changed, 70 insertions(+), 7 deletions(-)
diff --git a/dune/gfe/densities/bulkcosseratdensity.hh b/dune/gfe/densities/bulkcosseratdensity.hh
index 099b1583..824f46f6 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 32119a92..4af14f8e 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 c62af569..b23b05b0 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 644743c5..916ba6db 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 4b9de66d..00a5f1c2 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 f6887d7b..8c5c4504 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,
--
GitLab