diff --git a/dune/gfe/assemblers/localintegralenergy.hh b/dune/gfe/assemblers/localintegralenergy.hh
index aa4512473bd33d045f5cfd6baa985436893f75ce..385cb10b95695fe2475764651e1de27ed9ad3957 100644
--- a/dune/gfe/assemblers/localintegralenergy.hh
+++ b/dune/gfe/assemblers/localintegralenergy.hh
@@ -7,13 +7,11 @@
#include <dune/geometry/quadraturerules.hh>
#include <dune/gfe/assemblers/localenergy.hh>
-#include <dune/gfe/cosseratstrain.hh>
+#include <dune/gfe/densities/duneelasticitydensity.hh>
#include <dune/gfe/densities/localdensity.hh>
#include <dune/gfe/spaces/realtuple.hh>
#include <dune/gfe/spaces/rotation.hh>
-#include <dune/elasticity/materials/localdensity.hh>
-
namespace Dune::GFE {
/** \brief An energy given as an integral over a density
@@ -35,12 +33,6 @@ namespace Dune::GFE {
public:
- /** \brief Constructor with a Dune::Elasticity::LocalDensity
- */
- LocalIntegralEnergy(const std::shared_ptr<Elasticity::LocalDensity<gridDim,RT,DT> >& ld)
- : localDensityElasticity_(ld)
- {}
-
/** \brief Constructor with a Dune::GFE::LocalDensity
*/
LocalIntegralEnergy(const std::shared_ptr<GFE::LocalDensity<FieldVector<DT,gridDim>,TargetSpace> >& ld)
@@ -149,9 +141,35 @@ namespace Dune::GFE {
for (size_t comp=0; comp<deformationReferenceDerivative.N(); comp++)
jacobianInverseTransposed.mv(deformationReferenceDerivative[comp], deformationDerivative[comp]);
+ typename LocalOrientationGFEFunctionType::DerivativeType orientationDerivative;
+
// Integrate the energy density
- if (localDensityElasticity_)
- energy += weightWithintegrationElement * (*localDensityElasticity_)(x, deformationDerivative);
+ // We do a special-casing for DuneElasticityDensity:
+ // Of that one we know that it does not actually depend on the rotation.
+ using DEDensity = DuneElasticityDensity<FieldVector<DT,gridDim>,TargetSpace,0>;
+ std::shared_ptr<DEDensity> duneElasticityDensity = std::dynamic_pointer_cast<DEDensity>(localDensityGFE_);
+ if (duneElasticityDensity)
+ {
+ // Dummy local rotation value
+ ProductManifold<RealTuple<RT,gridDim>,Rotation<RT,gridDim> > value;
+ value[_0] = deformationValue;
+ value[_1] = TargetSpaceRotation::identity();
+
+ // Copy the two derivatives into a joint matrix object
+ // TODO: I am not sure about this. May the densities should get the
+ // separate derivatives.
+ FieldMatrix<RT,deformationDerivative.rows+orientationDerivative.rows, deformationDerivative.cols> derivative;
+
+ for (int i=0; i<deformationDerivative.rows; i++)
+ derivative[i] = deformationDerivative[i];
+
+ for (int i=0; i<orientationDerivative.rows; i++)
+ derivative[i+deformationDerivative.rows] = 0.0;
+
+ energy += weightWithintegrationElement * (*duneElasticityDensity)(x,
+ value,
+ derivative);
+ }
else if (localDensityGFE_) {
// The value of the local rotation
Rotation<RT,gridDim> orientationValue = localOrientationGFEFunction.evaluate(quadPos);
@@ -164,7 +182,6 @@ namespace Dune::GFE {
typename LocalOrientationGFEFunctionType::DerivativeType orientationReferenceDerivative = localOrientationGFEFunction.evaluateDerivative(quadPos,orientationValue);
// The derivative of the rotation defined on the actual element
- typename LocalOrientationGFEFunctionType::DerivativeType orientationDerivative;
for (size_t comp=0; comp<orientationReferenceDerivative.N(); comp++)
jacobianInverseTransposed.mv(orientationReferenceDerivative[comp], orientationDerivative[comp]);
@@ -190,8 +207,7 @@ namespace Dune::GFE {
}
protected:
- const std::shared_ptr<Elasticity::LocalDensity<gridDim,RT,DT> > localDensityElasticity_ = nullptr;
- const std::shared_ptr<GFE::LocalDensity<FieldVector<DT,gridDim>,TargetSpace> > localDensityGFE_ = nullptr;
+ const std::shared_ptr<GFE::LocalDensity<FieldVector<DT,gridDim>,TargetSpace> > localDensityGFE_;
};
} // namespace Dune::GFE
diff --git a/dune/gfe/densities/CMakeLists.txt b/dune/gfe/densities/CMakeLists.txt
index cea45c0e9d897b47477a9d570a8062154dff5b12..4101377ac2f47cb8cdfbcd69c4f35f967f17a449 100644
--- a/dune/gfe/densities/CMakeLists.txt
+++ b/dune/gfe/densities/CMakeLists.txt
@@ -3,6 +3,7 @@ install(FILES
bulkcosseratdensity.hh
chiralskyrmiondensity.hh
cosseratshelldensity.hh
+ duneelasticitydensity.hh
harmonicdensity.hh
localdensity.hh
DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/dune/gfe/densities)
diff --git a/dune/gfe/densities/duneelasticitydensity.hh b/dune/gfe/densities/duneelasticitydensity.hh
new file mode 100644
index 0000000000000000000000000000000000000000..9eb1ff98936a39eb80625aaa42cb1287bbe120b6
--- /dev/null
+++ b/dune/gfe/densities/duneelasticitydensity.hh
@@ -0,0 +1,65 @@
+#ifndef DUNE_GFE_DENSITIES_DUNEELASTICITYDENSITY_HH
+#define DUNE_GFE_DENSITIES_DUNEELASTICITYDENSITY_HH
+
+#include <dune/common/fmatrix.hh>
+
+#include <dune/gfe/densities/localdensity.hh>
+
+#include <dune/elasticity/materials/localdensity.hh>
+
+namespace Dune::GFE
+{
+
+ /** \brief Adapter for densities from the dune-elasticity module
+ *
+ * \tparam index The TargetSpace is typically a product with one factor being
+ * RealTuple. The 'index' argument selects which factor to assign the
+ * dune-elasticity density to.
+ */
+ template<class Position, class TargetSpace, std::size_t index=0>
+ class DuneElasticityDensity final
+ : public LocalDensity<Position,TargetSpace>
+ {
+ using ctype = typename Position::value_type;
+ using field_type = typename TargetSpace::field_type;
+
+ constexpr static auto dim = Position::size();
+ constexpr static auto embeddedBlocksize = TargetSpace::EmbeddedTangentVector::dimension;
+
+ public:
+
+ /** \brief Constructor with a Dune::Elasticity::LocalDensity
+ */
+ DuneElasticityDensity(const std::shared_ptr<Elasticity::LocalDensity<dim,field_type,ctype> >& elasticityDensity)
+ : elasticityDensity_(elasticityDensity)
+ {}
+
+ /** \brief Evaluate the density
+ *
+ * \param x The current position
+ * \param value The deformation at the current position
+ * \param derivative The derivative at the current position
+ */
+ virtual field_type operator() (const Position& x,
+ const TargetSpace& value,
+ const FieldMatrix<field_type,embeddedBlocksize,dim>& derivative) const override
+ {
+ // 'derivative' is the derivative of the entire TargetSpace point x.
+ // But (assuming that TargetSpace is a ProductManifold) we only want
+ // the parts that belongs to the factor given by the 'index' template parameter.
+ FieldMatrix<field_type,dim,dim> factorDerivative;
+
+ static_assert(index==0, "index!=0 is not implemented");
+ for (std::size_t i=0; i<dim; i++)
+ factorDerivative[i] = derivative[i];
+
+ return (*elasticityDensity_)(x, factorDerivative);
+ }
+
+ private:
+ const std::shared_ptr<Elasticity::LocalDensity<dim,field_type,ctype> > elasticityDensity_;
+ };
+
+} // namespace Dune::GFE
+
+#endif // DUNE_GFE_DENSITIES_DUNEELASTICITYDENSITY_HH
diff --git a/src/film-on-substrate.cc b/src/film-on-substrate.cc
index 97fc8ef8d204ed73bf7a83099e9be19e6980a47b..c8e0ee99878e4cf655292a2445b07b7f88b78618 100644
--- a/src/film-on-substrate.cc
+++ b/src/film-on-substrate.cc
@@ -62,6 +62,7 @@
#include <dune/gfe/neumannenergy.hh>
#include <dune/gfe/assemblers/surfacecosseratenergy.hh>
#include <dune/gfe/assemblers/sumenergy.hh>
+#include <dune/gfe/densities/duneelasticitydensity.hh>
#if MIXED_SPACE
#include <dune/gfe/mixedriemannianpnsolver.hh>
@@ -508,13 +509,17 @@ int main (int argc, char *argv[]) try
using ActiveRigidBodyMotion = GFE::ProductManifold<RealTuple<adouble,dim>, Rotation<adouble,dim> >;
+ // Wrap dune-elasticity density as dune-gfe density
+ auto elasticDensityWrapped = std::make_shared<GFE::DuneElasticityDensity<FieldVector<double,dim>,ActiveRigidBodyMotion,0> >(elasticDensity);
+
+
// Select which type of geometric interpolation to use
using LocalDeformationInterpolationRule = LocalGeodesicFEFunction<dim, GridType::ctype, decltype(deformationFEBasis.localView().tree().finiteElement()), RealTuple<adouble,dim> >;
using LocalOrientationInterpolationRule = LocalGeodesicFEFunction<dim, GridType::ctype, decltype(orientationFEBasis.localView().tree().finiteElement()), Rotation<adouble,dim> >;
using LocalInterpolationRule = std::tuple<LocalDeformationInterpolationRule,LocalOrientationInterpolationRule>;
- auto elasticEnergy = std::make_shared<GFE::LocalIntegralEnergy<CompositeBasis, LocalInterpolationRule, ActiveRigidBodyMotion> >(elasticDensity);
+ auto elasticEnergy = std::make_shared<GFE::LocalIntegralEnergy<CompositeBasis, LocalInterpolationRule, ActiveRigidBodyMotion> >(elasticDensityWrapped);
auto neumannEnergy = std::make_shared<GFE::NeumannEnergy<CompositeBasis, RealTuple<ValueType,targetDim>, Rotation<ValueType,dim> > >(neumannBoundary,neumannFunctionPtr);
auto surfaceCosseratEnergy = std::make_shared<GFE::SurfaceCosseratEnergy<
decltype(stressFreeShellFunction), CompositeBasis, ActiveRigidBodyMotion> >(
diff --git a/test/filmonsubstratetest.cc b/test/filmonsubstratetest.cc
index c9b40dda93a72713d1e2cc7b0034606e64856fed..a79ec09818520b909909aac0da22d6fb11f0d9b3 100644
--- a/test/filmonsubstratetest.cc
+++ b/test/filmonsubstratetest.cc
@@ -45,6 +45,7 @@
#include <dune/gfe/neumannenergy.hh>
#include <dune/gfe/assemblers/surfacecosseratenergy.hh>
#include <dune/gfe/assemblers/sumenergy.hh>
+#include <dune/gfe/densities/duneelasticitydensity.hh>
#if MIXED_SPACE
#include <dune/gfe/mixedriemannianpnsolver.hh>
@@ -331,10 +332,12 @@ int main (int argc, char *argv[])
materialParameters["b2"] = "1.0";
materialParameters["b3"] = "1.0";
+ using ActiveRigidBodyMotion = GFE::ProductManifold<RealTuple<adouble,dim>, Rotation<adouble,dim> >;
+
std::shared_ptr<Elasticity::LocalDensity<dim,ValueType> > elasticDensity;
elasticDensity = std::make_shared<Elasticity::MooneyRivlinDensity<dim,ValueType> >(materialParameters);
- using ActiveRigidBodyMotion = GFE::ProductManifold<RealTuple<adouble,dim>, Rotation<adouble,dim> >;
+ auto elasticDensityWrapped = std::make_shared<GFE::DuneElasticityDensity<FieldVector<double,dim>,ActiveRigidBodyMotion,0> >(elasticDensity);
// Select which type of geometric interpolation to use
using LocalDeformationInterpolationRule = LocalGeodesicFEFunction<dim, GridType::ctype, decltype(deformationFEBasis.localView().tree().finiteElement()), RealTuple<adouble,dim> >;
@@ -342,7 +345,7 @@ int main (int argc, char *argv[])
using LocalInterpolationRule = std::tuple<LocalDeformationInterpolationRule,LocalOrientationInterpolationRule>;
- auto elasticEnergy = std::make_shared<GFE::LocalIntegralEnergy<CompositeBasis, LocalInterpolationRule, ActiveRigidBodyMotion> >(elasticDensity);
+ auto elasticEnergy = std::make_shared<GFE::LocalIntegralEnergy<CompositeBasis, LocalInterpolationRule, ActiveRigidBodyMotion> >(elasticDensityWrapped);
auto neumannEnergy = std::make_shared<GFE::NeumannEnergy<CompositeBasis, RealTuple<ValueType,targetDim>, Rotation<ValueType,dim> > >(neumannBoundary,neumannFunction);
auto surfaceCosseratEnergy = std::make_shared<GFE::SurfaceCosseratEnergy<