dune/gfe/assemblers/localintegralstiffness.hh

@@ -11,7 +11,7 @@
 #include <dune/gfe/assemblers/localgeodesicfestiffness.hh>
 #include <dune/gfe/assemblers/localintegralenergy.hh>
 #include <dune/gfe/densities/localdensity.hh>
-#include <dune/gfe/interpolationderivatives.hh>
+#include <dune/gfe/functions/interpolationderivatives.hh>
 
 
 namespace Dune::GFE
@@ -32,6 +32,7 @@ namespace Dune::GFE
   class LocalIntegralStiffness
     : public LocalGeodesicFEStiffness<Basis,TargetSpace>
   {
+    using Element = typename Basis::GridView::template Codim<0>::Entity;
   public:
     constexpr static int gridDim = Basis::GridView::dimension;
 
@@ -54,8 +55,28 @@ namespace Dune::GFE
 
     using LocalCoordinate = typename GridView::template Codim<0>::Geometry::LocalCoordinate;
 
-    LocalIntegralStiffness(const std::shared_ptr<GFE::LocalDensity<LocalCoordinate,TargetSpace> >& ld)
-      : localDensity_(ld)
+    /** \brief Constructor from a GFE function as a shared pointer
+     *
+     * \param localGFEFunction The geometric finite element function that
+     * the density will be evaluated on
+     * \param density The density that will be integrated
+     */
+    LocalIntegralStiffness(std::shared_ptr<LocalInterpolationRule> localGFEFunction,
+                           const std::shared_ptr<GFE::LocalDensity<Element,TargetSpace> >& ld)
+      : localGFEFunction_(localGFEFunction)
+      , localDensity_(ld)
+    {}
+
+    /** \brief Constructor from a GFE function r-value reference
+     *
+     * \param localGFEFunction The geometric finite element function that
+     * the density will be evaluated on
+     * \param ld The density that will be integrated
+     */
+    LocalIntegralStiffness(LocalInterpolationRule&& localGFEFunction,
+                           const std::shared_ptr<GFE::LocalDensity<Element,TargetSpace> >& density)
+      : localGFEFunction_(std::make_shared<LocalInterpolationRule>(std::move(localGFEFunction)))
+      , localDensity_(density)
     {}
 
     virtual RT
@@ -66,16 +87,20 @@ namespace Dune::GFE
 
       if constexpr (not Impl::LocalEnergyTypes<TargetSpace>::isProductManifold)
       {
-        const auto& localFiniteElement = localView.tree().finiteElement();
-        LocalInterpolationRule localInterpolationRule(localFiniteElement,localCoefficients);
+        // Bind GFE function to the current element
+        const auto& element = localView.element();
+        localGFEFunction_->bind(element,localCoefficients);
 
-        int quadOrder = (localFiniteElement.type().isSimplex())
-           ? (localFiniteElement.localBasis().order()-1) * 2
-           : (localFiniteElement.localBasis().order() * gridDim - 1) * 2;
+        // Bind density to the element
+        localDensity_->bind(element);
 
-        const auto element = localView.element();
+        // Get a suitable quadrature rule
+        const auto localGFEOrder = localGFEFunction_->localFiniteElement().localBasis().order();
+        int quadOrder = (element.type().isSimplex())
+           ? (localGFEOrder-1) * 2
+           : (localGFEOrder * gridDim - 1) * 2;
 
-        const auto& quad = QuadratureRules<double, gridDim>::rule(localFiniteElement.type(), quadOrder);
+        const auto& quad = QuadratureRules<double, gridDim>::rule(element.type(), quadOrder);
 
         for (auto&& qp : quad)
         {
@@ -90,13 +115,13 @@ namespace Dune::GFE
           {
             if (localDensity_->dependsOnDerivative())
             {
-              auto [value, derivative] = localInterpolationRule.evaluateValueAndDerivative(quadPos);
+              auto [value, derivative] = localGFEFunction_->evaluateValueAndDerivative(quadPos);
               derivative = derivative * geometryJacobianInverse;
               energy += qp.weight() * integrationElement * (*localDensity_)(quadPos,value.globalCoordinates(),derivative);
             }
             else
             {
-              auto value = localInterpolationRule.evaluate(quadPos);
+              auto value = localGFEFunction_->evaluate(quadPos);
               typename LocalInterpolationRule::DerivativeType dummyDerivative;
               energy += qp.weight() * integrationElement * (*localDensity_)(quadPos,value.globalCoordinates(),dummyDerivative);
             }
@@ -106,7 +131,7 @@ namespace Dune::GFE
             if (localDensity_->dependsOnDerivative())
             {
               typename TargetSpace::CoordinateType dummyValue;
-              auto derivative = localInterpolationRule.evaluateDerivative(quadPos);
+              auto derivative = localGFEFunction_->evaluateDerivative(quadPos);
               derivative = derivative * geometryJacobianInverse;
               energy += qp.weight() * integrationElement * (*localDensity_)(quadPos,dummyValue,derivative);
             }
@@ -148,9 +173,14 @@ namespace Dune::GFE
                                             std::vector<double>& localGradient,
                                             typename GFE::Impl::LocalStiffnessTypes<TargetSpace>::Hessian& localHessian) const override
     {
-      LocalInterpolationRule localGFEFunction(localView.tree().finiteElement(),localCoefficients);
+      // Bind GFE function to the current element
+      const auto& element = localView.element();
+      localGFEFunction_->bind(element,localCoefficients);
+
+      // Bind density to the element
+      localDensity_->bind(element);
 
-      InterpolationDerivatives<LocalInterpolationRule> interpolationDerivatives(localGFEFunction,
+      InterpolationDerivatives<LocalInterpolationRule> interpolationDerivatives(*localGFEFunction_,
                                                                                 localDensity_->dependsOnValue(),
                                                                                 localDensity_->dependsOnDerivative());
 
@@ -175,14 +205,12 @@ namespace Dune::GFE
 
       const auto& localFiniteElement = localView.tree().finiteElement();
 
-      const auto& element = localView.element();
-
       // The range of input variables that the density depends on
       const size_t begin = (localDensity_->dependsOnValue()) ? 0 : TargetSpace::CoordinateType::dimension;
       const size_t end = (localDensity_->dependsOnDerivative()) ? m : TargetSpace::CoordinateType::dimension;
 
       // The quadrature rule
-      int quadOrder = (localFiniteElement.type().isSimplex())
+      int quadOrder = (element.type().isSimplex())
            ? (localFiniteElement.localBasis().order()-1) * 2
            : (localFiniteElement.localBasis().order() * gridDim - 1) * 2;
 
@@ -359,8 +387,12 @@ namespace Dune::GFE
                                  hessianDensity);
     }
 
+    // The value and derivative of this function are evaluated at the quadrature points,
+    // and given to the density.
+    const std::shared_ptr<LocalInterpolationRule> localGFEFunction_;
+
     // The density that is being integrated over
-    const std::shared_ptr<GFE::LocalDensity<LocalCoordinate,TargetSpace> > localDensity_ = nullptr;
+    const std::shared_ptr<GFE::LocalDensity<Element,TargetSpace> > localDensity_ = nullptr;
   };
 
 
@@ -372,6 +404,7 @@ namespace Dune::GFE
   class LocalIntegralStiffness<Basis,LocalInterpolationRule,ProductManifold<FactorSpaces...> >
     : public LocalGeodesicFEStiffness<Basis,ProductManifold<FactorSpaces...> >
   {
+    using Element = typename Basis::GridView::template Codim<0>::Entity;
   public:
 
     using TargetSpace = ProductManifold<FactorSpaces...>;
@@ -405,8 +438,28 @@ namespace Dune::GFE
 
     using LocalCoordinate = typename GridView::template Codim<0>::Geometry::LocalCoordinate;
 
-    LocalIntegralStiffness(const std::shared_ptr<GFE::LocalDensity<LocalCoordinate,TargetSpace> >& ld)
-      : localDensity_(ld)
+    /** \brief Constructor from a GFE function as a shared pointer
+     *
+     * \param localGFEFunction The geometric finite element function that
+     * the density will be evaluated on
+     * \param density The density that will be integrated
+     */
+    LocalIntegralStiffness(std::shared_ptr<LocalInterpolationRule> localGFEFunction,
+                           const std::shared_ptr<GFE::LocalDensity<Element,TargetSpace> >& ld)
+      : localGFEFunction_(localGFEFunction)
+      , localDensity_(ld)
+    {}
+
+    /** \brief Constructor from a GFE function r-value reference
+     *
+     * \param localGFEFunction The geometric finite element function that
+     * the density will be evaluated on
+     * \param ld The density that will be integrated
+     */
+    LocalIntegralStiffness(LocalInterpolationRule&& localGFEFunction,
+                           const std::shared_ptr<GFE::LocalDensity<Element,TargetSpace> >& density)
+      : localGFEFunction_(std::make_shared<LocalInterpolationRule>(std::move(localGFEFunction)))
+      , localDensity_(density)
     {}
 
     virtual ~LocalIntegralStiffness() {}
@@ -457,7 +510,12 @@ namespace Dune::GFE
                                             typename Dune::GFE::Impl::LocalStiffnessTypes<TargetSpace>::CompositeHessian& localHessian) const override;
 
   protected:
-    const std::shared_ptr<GFE::LocalDensity<LocalCoordinate,TargetSpace> > localDensity_ = nullptr;
+
+    // The value and derivative of this function are evaluated at the quadrature points,
+    // and given to the density.
+    const std::shared_ptr<LocalInterpolationRule> localGFEFunction_;
+
+    const std::shared_ptr<GFE::LocalDensity<Element,TargetSpace> > localDensity_ = nullptr;
 
   private:
 
@@ -522,16 +580,22 @@ namespace Dune::GFE
   {
     using namespace Dune::Indices;
 
+    // Bind GFE functions to the current element
+    const auto& element = localView.element();
+    std::get<0>(*localGFEFunction_).bind(element,localCoefficients[_0]);
+    std::get<1>(*localGFEFunction_).bind(element,localCoefficients[_1]);
+
+    // Bind density to the element
+    localDensity_->bind(element);
+
+    // Construct factory objects for the GFE derivatives
     using DeformationLocalInterpolationRule = typename std::tuple_element<0,LocalInterpolationRule>::type;
     using OrientationLocalInterpolationRule = typename std::tuple_element<1,LocalInterpolationRule>::type;
 
-    DeformationLocalInterpolationRule localDeformationGFEFunction(localView.tree().child(_0,0).finiteElement(),localCoefficients[_0]);
-    OrientationLocalInterpolationRule localOrientationGFEFunction(localView.tree().child(_1,0).finiteElement(),localCoefficients[_1]);
-
-    InterpolationDerivatives<DeformationLocalInterpolationRule> deformationInterpolationDerivatives(localDeformationGFEFunction,
+    InterpolationDerivatives<DeformationLocalInterpolationRule> deformationInterpolationDerivatives(std::get<0>(*localGFEFunction_),
                                                                                                     localDensity_->dependsOnValue(0),
                                                                                                     localDensity_->dependsOnDerivative(0));
-    InterpolationDerivatives<OrientationLocalInterpolationRule> orientationInterpolationDerivatives(localOrientationGFEFunction,
+    InterpolationDerivatives<OrientationLocalInterpolationRule> orientationInterpolationDerivatives(std::get<1>(*localGFEFunction_),
                                                                                                     localDensity_->dependsOnValue(1),
                                                                                                     localDensity_->dependsOnDerivative(1));
 
@@ -570,8 +634,10 @@ namespace Dune::GFE
 
     const auto& deformationLocalFiniteElement = localView.tree().child(_0,0).finiteElement();
 
-    const auto& element = localView.element();
+    // Bind density to the element
+    localDensity_->bind(element);
 
+    // Get a suitable quadrature rule
     int quadOrder = (element.type().isSimplex()) ? deformationLocalFiniteElement.localBasis().order()
                                                  : deformationLocalFiniteElement.localBasis().order() * gridDim;
 

dune/gfe/assemblers/mixedgfeassembler.hh

@@ -12,121 +12,126 @@
 #include <dune/gfe/assemblers/localgeodesicfestiffness.hh>
 
 
-/** \brief A global FE assembler for problems involving functions that map into non-Euclidean spaces
- */
-template <class Basis, class TargetSpace>
-class MixedGFEAssembler {
-
-  typedef typename Basis::GridView GridView;
-  typedef typename GridView::template Codim<0>::template Partition<Dune::Interior_Partition>::Iterator ElementIterator;
-
-  using TargetSpace0 = std::tuple_element_t<0,TargetSpace>;
-  using TargetSpace1 = std::tuple_element_t<1,TargetSpace>;
-  using LocalStiffness = LocalGeodesicFEStiffness<Basis, TargetSpace>;
-
-
-  //! Dimension of the grid.
-  constexpr static int gridDim = GridView::dimension;
-
-  //! Dimension of a tangent space
-  constexpr static int blocksize0 = TargetSpace0::TangentVector::dimension;
-  constexpr static int blocksize1 = TargetSpace1::TangentVector::dimension;
-
-  //!
-  typedef Dune::BCRSMatrix<Dune::FieldMatrix<double, blocksize0, blocksize0> > MatrixBlock00;
-  typedef Dune::BCRSMatrix<Dune::FieldMatrix<double, blocksize0, blocksize1> > MatrixBlock01;
-  typedef Dune::BCRSMatrix<Dune::FieldMatrix<double, blocksize1, blocksize0> > MatrixBlock10;
-  typedef Dune::BCRSMatrix<Dune::FieldMatrix<double, blocksize1, blocksize1> > MatrixBlock11;
-
-public:
-  typedef Dune::MultiTypeBlockMatrix<Dune::MultiTypeBlockVector<MatrixBlock00,MatrixBlock01>,
-      Dune::MultiTypeBlockVector<MatrixBlock10,MatrixBlock11> > MatrixType;
-  const Basis basis_;
-
-  std::shared_ptr<LocalStiffness> localStiffness_;
-
-public:
-
-  /** \brief Constructor for a given basis */
-  template <class LocalStiffnessT>
-  MixedGFEAssembler(const Basis& basis,
-                    LocalStiffnessT&& localStiffness)
-    : basis_(basis),
-    localStiffness_(Dune::Solvers::wrap_own_share<LocalStiffness>(std::forward<LocalStiffnessT>(localStiffness)))
-  {}
-
-  /** \brief Constructor for a given grid */
-  MixedGFEAssembler(const Basis& basis,
-                    LocalGeodesicFEStiffness<Basis, Dune::GFE::ProductManifold<TargetSpace0, TargetSpace1> >* localStiffness)
-    : basis_(basis),
-    localStiffness_(localStiffness)
-  {}
-
-  /** \brief Assemble the tangent stiffness matrix and the functional gradient together
-   *
-   * This is more efficient than computing them separately, because you need the gradient
-   * anyway to compute the Riemannian Hessian.
-   */
-  virtual void assembleGradientAndHessian(const std::vector<TargetSpace0>& configuration0,
-                                          const std::vector<TargetSpace1>& configuration1,
-                                          Dune::BlockVector<Dune::FieldVector<double, blocksize0> >& gradient0,
-                                          Dune::BlockVector<Dune::FieldVector<double, blocksize1> >& gradient1,
-                                          MatrixType& hessian,
-                                          bool computeOccupationPattern=true) const;
-
-  /** \brief Compute the energy of a deformation state */
-  virtual double computeEnergy(const std::vector<TargetSpace0>& configuration0,
-                               const std::vector<TargetSpace1>& configuration1) const;
-
-  //protected:
-  void getMatrixPattern(Dune::MatrixIndexSet& nb00,
-                        Dune::MatrixIndexSet& nb01,
-                        Dune::MatrixIndexSet& nb10,
-                        Dune::MatrixIndexSet& nb11) const;
-
-}; // end class
-
-
-
-template <class Basis, class TargetSpace>
-void MixedGFEAssembler<Basis,TargetSpace>::
-getMatrixPattern(Dune::MatrixIndexSet& nb00,
-                 Dune::MatrixIndexSet& nb01,
-                 Dune::MatrixIndexSet& nb10,
-                 Dune::MatrixIndexSet& nb11) const
+namespace Dune::GFE
 {
-  nb00.resize(basis_.size({0}), basis_.size({0}));
-  nb01.resize(basis_.size({0}), basis_.size({1}));
-  nb10.resize(basis_.size({1}), basis_.size({0}));
-  nb11.resize(basis_.size({1}), basis_.size({1}));
-
-  // A view on the FE basis on a single element
-  auto localView = basis_.localView();
 
-  // Loop over grid elements
-  for (const auto& element : elements(basis_.gridView(), Dune::Partitions::interior))
+  /** \brief A global FE assembler for problems involving functions that map into non-Euclidean spaces
+   */
+  template <class Basis, class TargetSpace>
+  class MixedGFEAssembler {
+
+    typedef typename Basis::GridView GridView;
+    typedef typename GridView::template Codim<0>::template Partition<Dune::Interior_Partition>::Iterator ElementIterator;
+
+    using TargetSpace0 = std::tuple_element_t<0,TargetSpace>;
+    using TargetSpace1 = std::tuple_element_t<1,TargetSpace>;
+    using LocalStiffness = LocalGeodesicFEStiffness<Basis, TargetSpace>;
+
+
+    //! Dimension of the grid.
+    constexpr static int gridDim = GridView::dimension;
+
+    //! Dimension of a tangent space
+    constexpr static int blocksize0 = TargetSpace0::TangentVector::dimension;
+    constexpr static int blocksize1 = TargetSpace1::TangentVector::dimension;
+
+    //!
+    typedef Dune::BCRSMatrix<Dune::FieldMatrix<double, blocksize0, blocksize0> > MatrixBlock00;
+    typedef Dune::BCRSMatrix<Dune::FieldMatrix<double, blocksize0, blocksize1> > MatrixBlock01;
+    typedef Dune::BCRSMatrix<Dune::FieldMatrix<double, blocksize1, blocksize0> > MatrixBlock10;
+    typedef Dune::BCRSMatrix<Dune::FieldMatrix<double, blocksize1, blocksize1> > MatrixBlock11;
+
+  public:
+    typedef Dune::MultiTypeBlockMatrix<Dune::MultiTypeBlockVector<MatrixBlock00,MatrixBlock01>,
+        Dune::MultiTypeBlockVector<MatrixBlock10,MatrixBlock11> > MatrixType;
+    const Basis basis_;
+
+    std::shared_ptr<LocalStiffness> localStiffness_;
+
+  public:
+
+    /** \brief Constructor for a given basis */
+    template <class LocalStiffnessT>
+    MixedGFEAssembler(const Basis& basis,
+                      LocalStiffnessT&& localStiffness)
+      : basis_(basis),
+      localStiffness_(Dune::Solvers::wrap_own_share<LocalStiffness>(std::forward<LocalStiffnessT>(localStiffness)))
+    {}
+
+    /** \brief Constructor for a given grid */
+    MixedGFEAssembler(const Basis& basis,
+                      LocalGeodesicFEStiffness<Basis, Dune::GFE::ProductManifold<TargetSpace0, TargetSpace1> >* localStiffness)
+      : basis_(basis),
+      localStiffness_(localStiffness)
+    {}
+
+    /** \brief Assemble the tangent stiffness matrix and the functional gradient together
+     *
+     * This is more efficient than computing them separately, because you need the gradient
+     * anyway to compute the Riemannian Hessian.
+     */
+    void assembleGradientAndHessian(const std::vector<TargetSpace0>& configuration0,
+                                    const std::vector<TargetSpace1>& configuration1,
+                                    Dune::BlockVector<Dune::FieldVector<double, blocksize0> >& gradient0,
+                                    Dune::BlockVector<Dune::FieldVector<double, blocksize1> >& gradient1,
+                                    MatrixType& hessian,
+                                    bool computeOccupationPattern=true) const;
+
+    /** \brief Compute the energy of a deformation state */
+    double computeEnergy(const std::vector<TargetSpace0>& configuration0,
+                         const std::vector<TargetSpace1>& configuration1) const;
+
+    //protected:
+    void getMatrixPattern(Dune::MatrixIndexSet& nb00,
+                          Dune::MatrixIndexSet& nb01,
+                          Dune::MatrixIndexSet& nb10,
+                          Dune::MatrixIndexSet& nb11) const;
+
+  }; // end class
+
+
+
+  template <class Basis, class TargetSpace>
+  void MixedGFEAssembler<Basis,TargetSpace>::
+  getMatrixPattern(Dune::MatrixIndexSet& nb00,
+                   Dune::MatrixIndexSet& nb01,
+                   Dune::MatrixIndexSet& nb10,
+                   Dune::MatrixIndexSet& nb11) const
   {
-    // Bind the local FE basis view to the current element
-    localView.bind(element);
-    // Add element stiffness matrix onto the global stiffness matrix
-    for (size_t i=0; i<localView.size(); i++)
-    {
-      // The global index of the i-th local degree of freedom of the element 'e'
-      auto row = localView.index(i);
+    nb00.resize(basis_.size({0}), basis_.size({0}));
+    nb01.resize(basis_.size({0}), basis_.size({1}));
+    nb10.resize(basis_.size({1}), basis_.size({0}));
+    nb11.resize(basis_.size({1}), basis_.size({1}));
+
+    // A view on the FE basis on a single element
+    auto localView = basis_.localView();
 
-      for (size_t j=0; j<localView.size(); j++ )
+    // Loop over grid elements
+    for (const auto& element : elements(basis_.gridView(), Dune::Partitions::interior))
+    {
+      // Bind the local FE basis view to the current element
+      localView.bind(element);
+      // Add element stiffness matrix onto the global stiffness matrix
+      for (size_t i=0; i<localView.size(); i++)
       {
-        // The global index of the j-th local degree of freedom of the element 'e'
-        auto col = localView.index(j);
-
-        if (row[0]==0 and col[0]==0)
-          nb00.add(row[1],col[1]);
-        if (row[0]==0 and col[0]==1)
-          nb01.add(row[1],col[1]);
-        if (row[0]==1 and col[0]==0)
-          nb10.add(row[1],col[1]);
-        if (row[0]==1 and col[0]==1)
-          nb11.add(row[1],col[1]);
+        // The global index of the i-th local degree of freedom of the element 'e'
+        auto row = localView.index(i);
+
+        for (size_t j=0; j<localView.size(); j++ )
+        {
+          // The global index of the j-th local degree of freedom of the element 'e'
+          auto col = localView.index(j);
+
+          if (row[0]==0 and col[0]==0)
+            nb00.add(row[1],col[1]);
+          if (row[0]==0 and col[0]==1)
+            nb01.add(row[1],col[1]);
+          if (row[0]==1 and col[0]==0)
+            nb10.add(row[1],col[1]);
+          if (row[0]==1 and col[0]==1)
+            nb11.add(row[1],col[1]);
+
+        }
 
       }
 
@@ -134,211 +139,211 @@ getMatrixPattern(Dune::MatrixIndexSet& nb00,
 
   }
 
-}
-
-template <class Basis, class TargetSpace>
-void MixedGFEAssembler<Basis,TargetSpace>::
-assembleGradientAndHessian(const std::vector<TargetSpace0>& configuration0,
-                           const std::vector<TargetSpace1>& configuration1,
-                           Dune::BlockVector<Dune::FieldVector<double, blocksize0> >& gradient0,
-                           Dune::BlockVector<Dune::FieldVector<double, blocksize1> >& gradient1,
-                           MatrixType& hessian,
-                           bool computeOccupationPattern) const
-{
-  if (computeOccupationPattern) {
+  template <class Basis, class TargetSpace>
+  void MixedGFEAssembler<Basis,TargetSpace>::
+  assembleGradientAndHessian(const std::vector<TargetSpace0>& configuration0,
+                             const std::vector<TargetSpace1>& configuration1,
+                             Dune::BlockVector<Dune::FieldVector<double, blocksize0> >& gradient0,
+                             Dune::BlockVector<Dune::FieldVector<double, blocksize1> >& gradient1,
+                             MatrixType& hessian,
+                             bool computeOccupationPattern) const
+  {
+    if (computeOccupationPattern) {
 
-    Dune::MatrixIndexSet pattern00;
-    Dune::MatrixIndexSet pattern01;
-    Dune::MatrixIndexSet pattern10;
-    Dune::MatrixIndexSet pattern11;
+      Dune::MatrixIndexSet pattern00;
+      Dune::MatrixIndexSet pattern01;
+      Dune::MatrixIndexSet pattern10;
+      Dune::MatrixIndexSet pattern11;
 
-    getMatrixPattern(pattern00, pattern01, pattern10, pattern11);
+      getMatrixPattern(pattern00, pattern01, pattern10, pattern11);
 
-    using namespace Dune::TypeTree::Indices;
-    pattern00.exportIdx(hessian[_0][_0]);
-    pattern01.exportIdx(hessian[_0][_1]);
-    pattern10.exportIdx(hessian[_1][_0]);
-    pattern11.exportIdx(hessian[_1][_1]);
+      using namespace Dune::Indices;
+      pattern00.exportIdx(hessian[_0][_0]);
+      pattern01.exportIdx(hessian[_0][_1]);
+      pattern10.exportIdx(hessian[_1][_0]);
+      pattern11.exportIdx(hessian[_1][_1]);
 
-  }
+    }
 
-  hessian = 0;
+    hessian = 0;
 
-  gradient0.resize(configuration0.size());
-  gradient0 = 0;
-  gradient1.resize(configuration1.size());
-  gradient1 = 0;
+    gradient0.resize(configuration0.size());
+    gradient0 = 0;
+    gradient1.resize(configuration1.size());
+    gradient1 = 0;
 
-  // A view on the FE basis on a single element
-  auto localView = basis_.localView();
-  for (const auto& element : elements(basis_.gridView(), Dune::Partitions::interior))
-  {
-    // Bind the local FE basis view to the current element
-    localView.bind(element);
-    using namespace Dune::TypeTree::Indices;
-
-    const int nDofs0 = localView.tree().child(_0,0).finiteElement().size();
-    const int nDofs1 = localView.tree().child(_1,0).finiteElement().size();
-    // This loop reads out the pattern for a local matrix; in each element, we have localView.size() degrees of freedom; from the composite and powerbasis layers
-    // nDofs0 are the degrees of freedom for *one* subspacebasis of the power basis of the displacement part;
-    // nDofs1 are the degrees of freedom for *one* subspacebasis of the power basis of the rotational part
-    // this is why the indices (_0,0) and (_1,0) are used: _0 takes the whole displacement part and _1 the whole rotational part; and 0 the first subspacebasis respectively
-    // Extract local solution
-    Dune::TupleVector<std::vector<TargetSpace0>, std::vector<TargetSpace1> > localConfiguration;
-    localConfiguration[_0].resize(nDofs0);
-    localConfiguration[_1].resize(nDofs1);
-
-    for (int i=0; i<nDofs0+nDofs1; i++)
+    // A view on the FE basis on a single element
+    auto localView = basis_.localView();
+    for (const auto& element : elements(basis_.gridView(), Dune::Partitions::interior))
     {
-      int localIndexI = 0;
-      if (i < nDofs0) {
-        auto& node = localView.tree().child(_0,0);
-        localIndexI = node.localIndex(i);
-      } else {
-        auto& node = localView.tree().child(_1,0);
-        localIndexI = node.localIndex(i-nDofs0);
+      // Bind the local FE basis view to the current element
+      localView.bind(element);
+      using namespace Dune::Indices;
+
+      const int nDofs0 = localView.tree().child(_0,0).finiteElement().size();
+      const int nDofs1 = localView.tree().child(_1,0).finiteElement().size();
+      // This loop reads out the pattern for a local matrix; in each element, we have localView.size() degrees of freedom; from the composite and powerbasis layers
+      // nDofs0 are the degrees of freedom for *one* subspacebasis of the power basis of the displacement part;
+      // nDofs1 are the degrees of freedom for *one* subspacebasis of the power basis of the rotational part
+      // this is why the indices (_0,0) and (_1,0) are used: _0 takes the whole displacement part and _1 the whole rotational part; and 0 the first subspacebasis respectively
+      // Extract local solution
+      Dune::TupleVector<std::vector<TargetSpace0>, std::vector<TargetSpace1> > localConfiguration;
+      localConfiguration[_0].resize(nDofs0);
+      localConfiguration[_1].resize(nDofs1);
+
+      for (int i=0; i<nDofs0+nDofs1; i++)
+      {
+        int localIndexI = 0;
+        if (i < nDofs0) {
+          auto& node = localView.tree().child(_0,0);
+          localIndexI = node.localIndex(i);
+        } else {
+          auto& node = localView.tree().child(_1,0);
+          localIndexI = node.localIndex(i-nDofs0);
+        }
+        auto multiIndex = localView.index(localIndexI);
+        //CompositeBasis number is contained in multiIndex[0], the Subspacebasis is contained in multiIndex[2]
+        //multiIndex[1] contains the actual index
+        if (multiIndex[0] == 0)
+          localConfiguration[_0][i] = configuration0[multiIndex[1]];
+        else if (multiIndex[0] == 1)
+          localConfiguration[_1][i-nDofs0] = configuration1[multiIndex[1]];
       }
-      auto multiIndex = localView.index(localIndexI);
-      //CompositeBasis number is contained in multiIndex[0], the Subspacebasis is contained in multiIndex[2]
-      //multiIndex[1] contains the actual index
-      if (multiIndex[0] == 0)
-        localConfiguration[_0][i] = configuration0[multiIndex[1]];
-      else if (multiIndex[0] == 1)
-        localConfiguration[_1][i-nDofs0] = configuration1[multiIndex[1]];
-    }
-
-    std::vector<double> localGradient(nDofs0*blocksize0 + nDofs1*blocksize1);
 
-    using HessianType = Dune::FieldMatrix<Dune::Matrix<double>,2,2>;
+      std::vector<double> localGradient(nDofs0*blocksize0 + nDofs1*blocksize1);
 
-    HessianType localHessian;
+      using HessianType = Dune::FieldMatrix<Dune::Matrix<double>,2,2>;
 
-    // setup local matrix and gradient
-    localStiffness_->assembleGradientAndHessian(localView,
-                                                localConfiguration,
-                                                localGradient,
-                                                localHessian);
-
-    // Add element matrix to global stiffness matrix
-    for (int i=0; i<nDofs0+nDofs1; i++)
-    {
-      int localIndexRow = 0;
-      if (i < nDofs0) {
-        auto& node = localView.tree().child(_0,0);
-        localIndexRow = node.localIndex(i);
-      } else {
-        auto& node = localView.tree().child(_1,0);
-        localIndexRow = node.localIndex(i-nDofs0);
-      }
+      HessianType localHessian;
 
-      auto row = localView.index(localIndexRow);
+      // setup local matrix and gradient
+      localStiffness_->assembleGradientAndHessian(localView,
+                                                  localConfiguration,
+                                                  localGradient,
+                                                  localHessian);
 
-      for (int j=0; j<nDofs0+nDofs1; j++ )
+      // Add element matrix to global stiffness matrix
+      for (int i=0; i<nDofs0+nDofs1; i++)
       {
-        int localIndexCol = 0;
-        if (j < nDofs0) {
+        int localIndexRow = 0;
+        if (i < nDofs0) {
           auto& node = localView.tree().child(_0,0);
-          localIndexCol = node.localIndex(j);
+          localIndexRow = node.localIndex(i);
         } else {
           auto& node = localView.tree().child(_1,0);
-          localIndexCol = node.localIndex(j-nDofs0);
+          localIndexRow = node.localIndex(i-nDofs0);
         }
 
-        auto col = localView.index(localIndexCol);
-
-        if (row[0]==0 and col[0]==0)
-          for (std::size_t ii=0; ii<blocksize0; ii++)
-            for (std::size_t jj=0; jj<blocksize0; jj++)
-              hessian[_0][_0][row[1]][col[1]][ii][jj] += localHessian[_0][_0][i*blocksize0+ii][j*blocksize0+jj];
-
-        if (row[0]==0 and col[0]==1)
-          for (std::size_t ii=0; ii<blocksize0; ii++)
-            for (std::size_t jj=0; jj<blocksize1; jj++)
-              hessian[_0][_1][row[1]][col[1]][ii][jj] += localHessian[_0][_1][i*blocksize0+ii][(j-nDofs0)*blocksize1+jj];
-
-        if (row[0]==1 and col[0]==0)
-          for (std::size_t ii=0; ii<blocksize1; ii++)
-            for (std::size_t jj=0; jj<blocksize0; jj++)
-              hessian[_1][_0][row[1]][col[1]][ii][jj] += localHessian[_1][_0][(i-nDofs0)*blocksize1+ii][j*blocksize0+jj];
+        auto row = localView.index(localIndexRow);
+
+        for (int j=0; j<nDofs0+nDofs1; j++ )
+        {
+          int localIndexCol = 0;
+          if (j < nDofs0) {
+            auto& node = localView.tree().child(_0,0);
+            localIndexCol = node.localIndex(j);
+          } else {
+            auto& node = localView.tree().child(_1,0);
+            localIndexCol = node.localIndex(j-nDofs0);
+          }
+
+          auto col = localView.index(localIndexCol);
+
+          if (row[0]==0 and col[0]==0)
+            for (std::size_t ii=0; ii<blocksize0; ii++)
+              for (std::size_t jj=0; jj<blocksize0; jj++)
+                hessian[_0][_0][row[1]][col[1]][ii][jj] += localHessian[_0][_0][i*blocksize0+ii][j*blocksize0+jj];
+
+          if (row[0]==0 and col[0]==1)
+            for (std::size_t ii=0; ii<blocksize0; ii++)
+              for (std::size_t jj=0; jj<blocksize1; jj++)
+                hessian[_0][_1][row[1]][col[1]][ii][jj] += localHessian[_0][_1][i*blocksize0+ii][(j-nDofs0)*blocksize1+jj];
+
+          if (row[0]==1 and col[0]==0)
+            for (std::size_t ii=0; ii<blocksize1; ii++)
+              for (std::size_t jj=0; jj<blocksize0; jj++)
+                hessian[_1][_0][row[1]][col[1]][ii][jj] += localHessian[_1][_0][(i-nDofs0)*blocksize1+ii][j*blocksize0+jj];
+
+          if (row[0]==1 and col[0]==1)
+            for (std::size_t ii=0; ii<blocksize1; ii++)
+              for (std::size_t jj=0; jj<blocksize1; jj++)
+                hessian[_1][_1][row[1]][col[1]][ii][jj] += localHessian[_1][_1][(i-nDofs0)*blocksize1+ii][(j-nDofs0)*
+                                                                                                          blocksize1+jj];
+        }
 
-        if (row[0]==1 and col[0]==1)
-          for (std::size_t ii=0; ii<blocksize1; ii++)
-            for (std::size_t jj=0; jj<blocksize1; jj++)
-              hessian[_1][_1][row[1]][col[1]][ii][jj] += localHessian[_1][_1][(i-nDofs0)*blocksize1+ii][(j-nDofs0)*
-                                                                                                        blocksize1+jj];
+        // Add local gradient to global gradient
+        if (row[0] == 0)
+          for (std::size_t j=0; j<blocksize0; j++)
+            gradient0[row[1]][j] += localGradient[i*blocksize0 + j];
+        else
+          for (std::size_t j=0; j<blocksize1; j++)
+            gradient1[row[1]][j] += localGradient[nDofs0*blocksize0 + (i-nDofs0)*blocksize1 + j];
       }
 
-      // Add local gradient to global gradient
-      if (row[0] == 0)
-        for (std::size_t j=0; j<blocksize0; j++)
-          gradient0[row[1]][j] += localGradient[i*blocksize0 + j];
-      else
-        for (std::size_t j=0; j<blocksize1; j++)
-          gradient1[row[1]][j] += localGradient[nDofs0*blocksize0 + (i-nDofs0)*blocksize1 + j];
     }
-
   }
-}
 
-template <class Basis, class TargetSpace>
-double MixedGFEAssembler<Basis, TargetSpace>::
-computeEnergy(const std::vector<TargetSpace0>& configuration0,
-              const std::vector<TargetSpace1>& configuration1) const
-{
-  double energy = 0;
+  template <class Basis, class TargetSpace>
+  double MixedGFEAssembler<Basis, TargetSpace>::
+  computeEnergy(const std::vector<TargetSpace0>& configuration0,
+                const std::vector<TargetSpace1>& configuration1) const
+  {
+    double energy = 0;
 
-  if (configuration0.size()!=basis_.size({0}))
-    DUNE_THROW(Dune::Exception, "Configuration vector 0 doesn't match the basis!");
+    if (configuration0.size()!=basis_.size({0}))
+      DUNE_THROW(Dune::Exception, "Configuration vector 0 doesn't match the basis!");
 
-  if (configuration1.size()!=basis_.size({1}))
-    DUNE_THROW(Dune::Exception, "Configuration vector 1 doesn't match the basis!");
+    if (configuration1.size()!=basis_.size({1}))
+      DUNE_THROW(Dune::Exception, "Configuration vector 1 doesn't match the basis!");
 
-  // A view on the FE basis on a single element
-  auto localView = basis_.localView();
+    // A view on the FE basis on a single element
+    auto localView = basis_.localView();
 
-  // Loop over all elements
-  for (const auto& element : elements(basis_.gridView(), Dune::Partitions::interior))
-  {
-    // Bind the local FE basis view to the current element
-    localView.bind(element);
+    // Loop over all elements
+    for (const auto& element : elements(basis_.gridView(), Dune::Partitions::interior))
+    {
+      // Bind the local FE basis view to the current element
+      localView.bind(element);
 
-    // Number of degrees of freedom on this element
-    using namespace Dune::Indices;
-    const int nDofs0 = localView.tree().child(_0,0).finiteElement().size();
-    const int nDofs1 = localView.tree().child(_1,0).finiteElement().size();
+      // Number of degrees of freedom on this element
+      using namespace Dune::Indices;
+      const int nDofs0 = localView.tree().child(_0,0).finiteElement().size();
+      const int nDofs1 = localView.tree().child(_1,0).finiteElement().size();
 
-    Dune::TupleVector<std::vector<TargetSpace0>, std::vector<TargetSpace1> > localConfiguration;
-    localConfiguration[_0].resize(nDofs0);
-    localConfiguration[_1].resize(nDofs1);
+      Dune::TupleVector<std::vector<TargetSpace0>, std::vector<TargetSpace1> > localConfiguration;
+      localConfiguration[_0].resize(nDofs0);
+      localConfiguration[_1].resize(nDofs1);
 
-    for (int i=0; i<nDofs0+nDofs1; i++)
-    {
-      int localIndexI = 0;
-      if (i < nDofs0) {
-        auto& node = localView.tree().child(_0,0);
-        localIndexI = node.localIndex(i);
-      } else {
-        auto& node = localView.tree().child(_1,0);
-        localIndexI = node.localIndex(i-nDofs0);
+      for (int i=0; i<nDofs0+nDofs1; i++)
+      {
+        int localIndexI = 0;
+        if (i < nDofs0) {
+          auto& node = localView.tree().child(_0,0);
+          localIndexI = node.localIndex(i);
+        } else {
+          auto& node = localView.tree().child(_1,0);
+          localIndexI = node.localIndex(i-nDofs0);
+        }
+
+        auto multiIndex = localView.index(localIndexI);
+
+        // The CompositeBasis number is contained in multiIndex[0]
+        // multiIndex[1] contains the actual index
+        if (multiIndex[0] == 0)
+          localConfiguration[_0][i] = configuration0[multiIndex[1]];
+        else if (multiIndex[0] == 1)
+          localConfiguration[_1][i-nDofs0] = configuration1[multiIndex[1]];
       }
 
-      auto multiIndex = localView.index(localIndexI);
+      energy += localStiffness_->energy(localView,
+                                        localConfiguration);
 
-      // The CompositeBasis number is contained in multiIndex[0]
-      // multiIndex[1] contains the actual index
-      if (multiIndex[0] == 0)
-        localConfiguration[_0][i] = configuration0[multiIndex[1]];
-      else if (multiIndex[0] == 1)
-        localConfiguration[_1][i-nDofs0] = configuration1[multiIndex[1]];
     }
 
-    energy += localStiffness_->energy(localView,
-                                      localConfiguration);
+    return energy;
 
   }
 
-  return energy;
-
-}
+}  // namespace Dune::GFE
 
 #endif

dune/gfe/assemblers/nonplanarcosseratshellenergy.hh

@@ -5,9 +5,9 @@
 #include <dune/common/parametertree.hh>
 #include <dune/geometry/quadraturerules.hh>
 
-#include <dune/matrix-vector/crossproduct.hh>
+#include <dune/functions/functionspacebases/subspacebasis.hh>
 
-#include <dune/fufem/boundarypatch.hh>
+#include <dune/matrix-vector/crossproduct.hh>
 
 #include <dune/functions/gridfunctions/discreteglobalbasisfunction.hh>
 
@@ -15,11 +15,8 @@
 #include <dune/gmsh4/gridcreators/lagrangegridcreator.hh>
 #endif
 
-#include <dune/gfe/assemblers/cosseratenergystiffness.hh>
 #include <dune/gfe/assemblers/localenergy.hh>
-#include <dune/gfe/localgeodesicfefunction.hh>
-#include <dune/gfe/tensor3.hh>
-#include <dune/gfe/localprojectedfefunction.hh>
+#include <dune/gfe/densities/cosseratshelldensity.hh>
 #include <dune/gfe/spaces/productmanifold.hh>
 #include <dune/gfe/spaces/realtuple.hh>
 #include <dune/gfe/spaces/rotation.hh>
@@ -29,650 +26,318 @@
 #include <dune/localfunctions/lagrange/lfecache.hh>
 #endif
 
-/** \brief Assembles the cosserat energy for a single element.
- *
- * \tparam Basis                       Type of the Basis used for assembling
- * \tparam dim                         Dimension of the Targetspace, 3
- * \tparam field_type                  The coordinate type of the TargetSpace
- * \tparam StressFreeStateGridFunction Type of the GridFunction representing the Cosserat shell in a stress free state
- */
-template<class Basis, int dim, class field_type, class StressFreeStateGridFunction>
-class NonplanarCosseratShellEnergy
-  : public Dune::GFE::LocalEnergy<Basis,Dune::GFE::ProductManifold<RealTuple<field_type,dim>,Rotation<field_type,dim> > >
+namespace Dune::GFE
 {
-  // grid types
-  typedef typename Basis::GridView GridView;
-  typedef typename GridView::ctype DT;
-  typedef Dune::GFE::ProductManifold<RealTuple<field_type,dim>,Rotation<field_type,dim> > TargetSpace;
-  typedef typename TargetSpace::ctype RT;
-  typedef typename GridView::template Codim<0>::Entity Entity;
-
-  // some other sizes
-  constexpr static int gridDim = GridView::dimension;
-  constexpr static int dimworld = GridView::dimensionworld;
-
-public:
-
-  /** \brief Constructor with a set of material parameters
-   * \param parameters                  The material parameters
-   * \param stressFreeStateGridFunction Pointer to a parametrization representing the Cosserat shell in a stress-free state
+  /** \brief Assembles the cosserat energy for a single element.
+   *
+   * \tparam Basis                       Type of the Basis used for assembling
+   * \tparam LocalGFEFunction The geometric finite element function that represents the shell configuration
+   * \tparam dim                         Dimension of the Targetspace, 3
+   * \tparam field_type                  The coordinate type of the TargetSpace
+   * \tparam StressFreeStateGridFunction Type of the GridFunction representing the Cosserat shell in a stress free state
    */
-  NonplanarCosseratShellEnergy(const Dune::ParameterTree& parameters,
-                               const StressFreeStateGridFunction* stressFreeStateGridFunction,
-                               const BoundaryPatch<GridView>* neumannBoundary,
-                               const std::function<Dune::FieldVector<double,3>(Dune::FieldVector<double,dimworld>)> neumannFunction,
-                               const std::function<Dune::FieldVector<double,3>(Dune::FieldVector<double,dimworld>)> volumeLoad)
-    : stressFreeStateGridFunction_(stressFreeStateGridFunction),
-    neumannBoundary_(neumannBoundary),
-    neumannFunction_(neumannFunction),
-    volumeLoad_(volumeLoad)
-  {
-    // The shell thickness
-    thickness_ = parameters.template get<double>("thickness");
-
-    // Lame constants
-    mu_ = parameters.template get<double>("mu");
-    lambda_ = parameters.template get<double>("lambda");
-
-    // Cosserat couple modulus
-    mu_c_ = parameters.template get<double>("mu_c");
-
-    // Length scale parameter
-    L_c_ = parameters.template get<double>("L_c");
-
-    // Curvature parameters
-    b1_ = parameters.template get<double>("b1");
-    b2_ = parameters.template get<double>("b2");
-    b3_ = parameters.template get<double>("b3");
-
-    // Indicator to use the alternative energy W_Coss from Birsan 2021:
-    useAlternativeEnergyWCoss_ = parameters.template get<bool>("useAlternativeEnergyWCoss", false);
-  }
-
-  /** \brief Assemble the energy for a single element */
-  RT energy (const typename Basis::LocalView& localView,
-             const std::vector<TargetSpace>& localSolution) const override;
-
-  RT energy (const typename Basis::LocalView& localView,
-             const typename Dune::GFE::Impl::LocalEnergyTypes<TargetSpace>::CompositeCoefficients& coefficients) const override;
-
-  /*  Sources:
-      Birsan 2019: Derivation of a refined six-parameter shell model, equation (111)
-      Birsan 2021: Alternative derivation of the higher-order constitudtive model for six-parameter elastic shells, equations (119) and (126)
-   */
-  RT W_Coss(const Dune::FieldMatrix<field_type,3,3>& S, const Dune::FieldMatrix<double,3,3>& a, const Dune::FieldVector<double,3>& n0) const
-  {
-    return W_Coss_mixt(S,S,a,n0);
-  }
-
-  RT W_Coss_mixt(const Dune::FieldMatrix<field_type,3,3>& S, const Dune::FieldMatrix<field_type,3,3>& T, const Dune::FieldMatrix<double,3,3>& a, const Dune::FieldVector<double,3>& n0) const
-  {
-    auto planarPart = W_mixt(a*S,a*T);
-    Dune::FieldVector<field_type, 3> n0S;
-    Dune::FieldVector<field_type, 3> n0T;
-    S.mtv(n0, n0S);
-    T.mtv(n0, n0T);
-    field_type normalPart = 2*mu_*mu_c_* n0S * n0T /(mu_ + mu_c_);
-    return planarPart + normalPart;
-  }
-
-  RT W_m(const Dune::FieldMatrix<field_type,3,3>& S) const
+  template<class Basis, class LocalGFEFunction, int dim, class field_type, class StressFreeStateGridFunction>
+  class NonplanarCosseratShellEnergy
+    : public Dune::GFE::LocalEnergy<Basis,Dune::GFE::ProductManifold<RealTuple<field_type,dim>,Rotation<field_type,dim> > >
   {
-    return W_mixt(S,S);
-  }
-
-  RT W_mixt(const Dune::FieldMatrix<field_type,3,3>& S, const Dune::FieldMatrix<field_type,3,3>& T) const
-  {
-    return mu_ * Dune::GFE::frobeniusProduct(Dune::GFE::sym(S), Dune::GFE::sym(T))
-           + mu_c_ * Dune::GFE::frobeniusProduct(Dune::GFE::skew(S), Dune::GFE::skew(T))
-           + lambda_ * mu_ / (lambda_ + 2*mu_) * Dune::GFE::trace(S) * Dune::GFE::trace(T);
-  }
-
-  RT W_mp(const Dune::FieldMatrix<field_type,3,3>& S) const
-  {
-    return mu_ * Dune::GFE::sym(S).frobenius_norm2() + mu_c_ * Dune::GFE::skew(S).frobenius_norm2() + lambda_ * 0.5 * Dune::GFE::traceSquared(S);
-  }
-
-  //For b1 = b2 = 1 and b3 = 1/3, this reduces to S.frobenius_norm2()
-  RT W_curv(const Dune::FieldMatrix<field_type,3,3>& S) const
-  {
-    return mu_ * L_c_ * L_c_ * (b1_ * Dune::GFE::dev(Dune::GFE::sym(S)).frobenius_norm2()
-                                + b2_ * Dune::GFE::skew(S).frobenius_norm2() + b3_ * Dune::GFE::traceSquared(S));
-  }
+    // grid types
+    typedef typename Basis::GridView GridView;
+    typedef typename GridView::ctype DT;
+    typedef Dune::GFE::ProductManifold<RealTuple<field_type,dim>,Rotation<field_type,dim> > TargetSpace;
+    typedef typename TargetSpace::ctype RT;
+    typedef typename GridView::template Codim<0>::Entity Element;
+
+    // some other sizes
+    constexpr static int gridDim = GridView::dimension;
+    constexpr static int dimworld = GridView::dimensionworld;
+
+  public:
+
+    /** \brief Constructor from a GFE function as a shared pointer
+     * \param stressFreeStateGridFunction Pointer to a parametrization representing the Cosserat shell in a stress-free state
+     */
+    NonplanarCosseratShellEnergy(std::shared_ptr<LocalGFEFunction> localGFEFunction,
+                                 const std::shared_ptr<Dune::GFE::CosseratShellDensity<Element,field_type> >& density,
+                                 const StressFreeStateGridFunction* stressFreeStateGridFunction)
+      : localGFEFunction_(localGFEFunction)
+      , stressFreeStateGridFunction_(stressFreeStateGridFunction)
+      , density_(density)
+    {}
+
+    /** \brief Constructor from a GFE function as an r-value reference
+     * \param stressFreeStateGridFunction Pointer to a parametrization representing the Cosserat shell in a stress-free state
+     */
+    NonplanarCosseratShellEnergy(LocalGFEFunction&& localGFEFunction,
+                                 const std::shared_ptr<Dune::GFE::CosseratShellDensity<Element,field_type> >& density,
+                                 const StressFreeStateGridFunction* stressFreeStateGridFunction)
+      : localGFEFunction_(std::make_shared<LocalGFEFunction>(std::move(localGFEFunction)))
+      , stressFreeStateGridFunction_(stressFreeStateGridFunction)
+      , density_(density)
+    {}
+
+    /** \brief Assemble the energy for a single element */
+    RT energy (const typename Basis::LocalView& localView,
+               const std::vector<TargetSpace>& localSolution) const override;
+
+    RT energy (const typename Basis::LocalView& localView,
+               const typename Dune::GFE::Impl::LocalEnergyTypes<TargetSpace>::CompositeCoefficients& coefficients) const override;
 
 #if HAVE_DUNE_GMSH4
-  static int getOrder(const Dune::Gmsh4::LagrangeGridCreator<typename GridView::Grid>* lagrangeGridCreator)
-  {
-    return lagrangeGridCreator->order();
-  }
+    static int getOrder(const Dune::Gmsh4::LagrangeGridCreator<typename GridView::Grid>* lagrangeGridCreator)
+    {
+      return lagrangeGridCreator->order();
+    }
 #endif
 
-  template<typename B, typename V, typename NTREM, typename R>
-  static int getOrder(const Dune::Functions::DiscreteGlobalBasisFunction<B,V, NTREM, R>* gridFunction)
-  {
-    return gridFunction->basis().preBasis().subPreBasis().order();
-  }
-
-  /** \brief The shell thickness */
-  double thickness_;
-
-  /** \brief Lame constants */
-  double mu_, lambda_;
-
-  /** \brief Cosserat couple modulus */
-  double mu_c_;
-
-  /** \brief Length scale parameter */
-  double L_c_;
-
-  /** \brief Curvature parameters */
-  double b1_, b2_, b3_;
-
-  /** \brief Indicator to use the alternative energy W_Coss from Birsan 2021:
-             Alternative derivation of the higher-order constitudtive model for six-parameter elastic shells, equations (119) and (126). */
-  bool useAlternativeEnergyWCoss_;
-
-  /** \brief The geometry of the reference deformation used for assembling */
-  const StressFreeStateGridFunction* stressFreeStateGridFunction_;
+    template<typename B, typename V, typename NTREM, typename R>
+    static int getOrder(const Dune::Functions::DiscreteGlobalBasisFunction<B,V, NTREM, R>* gridFunction)
+    {
+      return gridFunction->basis().preBasis().subPreBasis().order();
+    }
 
-  /** \brief The Neumann boundary */
-  const BoundaryPatch<GridView>* neumannBoundary_;
+    // The value and derivative of this function are evaluated at the quadrature points,
+    // and given to the density.
+    const std::shared_ptr<LocalGFEFunction> localGFEFunction_;
 
-  /** \brief The function implementing the Neumann data */
-  const std::function<Dune::FieldVector<double,3>(Dune::FieldVector<double,dimworld>)> neumannFunction_;
+    /** \brief The geometry of the reference deformation used for assembling */
+    const StressFreeStateGridFunction* stressFreeStateGridFunction_;
 
-  /** \brief The function implementing a volume load */
-  const std::function<Dune::FieldVector<double,3>(Dune::FieldVector<double,dimworld>)> volumeLoad_;
-};
+    /** \brief The energy density of a Cosserat shell with nonplanar reference configuration */
+    const std::shared_ptr<Dune::GFE::CosseratShellDensity<Element,field_type> > density_;
+  };
 
-template <class Basis, int dim, class field_type, class StressFreeStateGridFunction>
-[[deprecated("Use an std::vector<RealTuple<field_type,dim>> and an std::vector<Rotation<field_type,dim>> together with the MixedGFEAssembler and the GFEAssemblerWrapper instead of std::vector<Dune::GFE::ProductManifold<RealTuple<field_type,dim>,Rotation<field_type,dim> >>.")]]
-typename NonplanarCosseratShellEnergy<Basis, dim, field_type, StressFreeStateGridFunction>::RT
-NonplanarCosseratShellEnergy<Basis,dim,field_type, StressFreeStateGridFunction>::
-energy(const typename Basis::LocalView& localView,
-       const std::vector<Dune::GFE::ProductManifold<RealTuple<field_type,dim>,Rotation<field_type,dim> > >& localSolution) const
-{
-  using namespace Dune::Indices;
+  template <class Basis, class LocalGFEFunction, int dim, class field_type, class StressFreeStateGridFunction>
+  typename NonplanarCosseratShellEnergy<Basis, LocalGFEFunction, dim, field_type, StressFreeStateGridFunction>::RT
+  NonplanarCosseratShellEnergy<Basis,LocalGFEFunction,dim,field_type, StressFreeStateGridFunction>::
+  energy(const typename Basis::LocalView& localView,
+         const std::vector<Dune::GFE::ProductManifold<RealTuple<field_type,dim>,Rotation<field_type,dim> > >& localSolution) const
+  {
+    RT energy = 0;
 
-  // The element geometry
-  auto element = localView.element();
+    // The element geometry
+    auto element = localView.element();
 
-  // The set of shape functions on this element
-  using namespace Dune::TypeTree::Indices;
-  const auto& localFiniteElement = LocalFiniteElementFactory<Basis,0>::get(localView,_0);
+    if constexpr (Basis::LocalView::Tree::isLeaf || Basis::LocalView::Tree::isPower)
+    {
+      // The set of shape functions on this element
+      using namespace Dune::Indices;
 
 #if HAVE_DUNE_CURVEDGEOMETRY
-  // Construct a curved geometry of this element of the Cosserat shell in its stress-free state
-  // The variable local holds the local coordinates in the reference element
-  // and localGeometry.global maps them to the world coordinates
-  Dune::CurvedGeometry<DT, gridDim, dimworld, Dune::CurvedGeometryTraits<DT, Dune::LagrangeLFECache<DT,DT,gridDim> > > geometry(referenceElement(element),
-                                                                                                                                [this,element](const auto& local) {
-                                                                                                                                auto localGridFunction = localFunction(*stressFreeStateGridFunction_);
-                                                                                                                                localGridFunction.bind(element);
-                                                                                                                                return localGridFunction(local);
-    }, getOrder(stressFreeStateGridFunction_));
+      // Construct a curved geometry of this element of the Cosserat shell in its stress-free state
+      // The variable local holds the local coordinates in the reference element
+      // and localGeometry.global maps them to the world coordinates
+      Dune::CurvedGeometry<DT, gridDim, dimworld, Dune::CurvedGeometryTraits<DT, Dune::LagrangeLFECache<DT,DT,gridDim> > > geometry(referenceElement(element),
+                                                                                                                                    [this,element](const auto& local) {
+                                                                                                                                    auto localGridFunction = localFunction(*stressFreeStateGridFunction_);
+                                                                                                                                    localGridFunction.bind(element);
+                                                                                                                                    return localGridFunction(local);
+        }, getOrder(stressFreeStateGridFunction_));
 #else
-  // When using element.geometry(), the geometry of the element is flat
-  auto geometry = element.geometry();
+      // When using element.geometry(), the geometry of the element is flat
+      auto geometry = element.geometry();
 #endif
 
-  ////////////////////////////////////////////////////////////////////////////////////
-  //  Set up the local nonlinear finite element function
-  ////////////////////////////////////////////////////////////////////////////////////
-  typedef LocalGeodesicFEFunction<gridDim, DT, decltype(localFiniteElement), TargetSpace> LocalGFEFunctionType;
-  LocalGFEFunctionType localGeodesicFEFunction(localFiniteElement,localSolution);
-
-  RT energy = 0;
-
-  auto quadOrder = (element.type().isSimplex()) ? localFiniteElement.localBasis().order()
-                                                : localFiniteElement.localBasis().order() * gridDim;
+      localGFEFunction_->bind(element,localSolution);
 
-  const auto& quad = Dune::QuadratureRules<DT, gridDim>::rule(element.type(), quadOrder);
+      // Bind the density to the current element
+      density_->bind(element);
 
-  for (size_t pt=0; pt<quad.size(); pt++)
-  {
-    // Local position of the quadrature point
-    const Dune::FieldVector<DT,gridDim>& quadPos = quad[pt].position();
-
-    const DT integrationElement = geometry.integrationElement(quadPos);
+      const auto localGFEOrder = localGFEFunction_->localFiniteElement().localBasis().order();
+      const auto quadOrder = (element.type().isSimplex()) ? localGFEOrder
+                                                        : localGFEOrder * gridDim;
 
-    // The value of the local function
-    Dune::GFE::ProductManifold<RealTuple<field_type,dim>,Rotation<field_type,dim> > value = localGeodesicFEFunction.evaluate(quadPos);
+      const auto& quad = Dune::QuadratureRules<DT, gridDim>::rule(element.type(), quadOrder);
 
-    // The derivative of the local function w.r.t. the coordinate system of the tangent space
-    auto derivative = localGeodesicFEFunction.evaluateDerivative(quadPos,value);
+      for (size_t pt=0; pt<quad.size(); pt++)
+      {
+        // Local position of the quadrature point
+        const Dune::FieldVector<DT,gridDim>& quadPos = quad[pt].position();
 
-    //////////////////////////////////////////////////////////
-    //  The rotation and its derivative
-    //  Note: we need it in matrix coordinates
-    //////////////////////////////////////////////////////////
+        const DT integrationElement = geometry.integrationElement(quadPos);
 
-    Dune::FieldMatrix<field_type,dim,dim> R;
-    value[_1].matrix(R);
-    auto RT = Dune::GFE::transpose(R);
+        // The value of the local function
+        const auto value = localGFEFunction_->evaluate(quadPos);
 
-    // Extract the orientation derivative
-    Dune::FieldMatrix<field_type,4,gridDim> orientationDerivative;
-    for (size_t i=0; i<4; ++i)
-      for (size_t j=0; j<gridDim; ++j)
-        orientationDerivative[i][j] = derivative[3+i][j];
+        // The derivative of the local function w.r.t. the coordinate system of the tangent space
+        const auto derivative = localGFEFunction_->evaluateDerivative(quadPos,value);
 
-    //Derivative of the rotation w.r.t. the coordinate system of the tangent space
-    Tensor3<field_type,3,3,gridDim> DR = value[_1].quaternionTangentToMatrixTangent(orientationDerivative);
+        ////////////////////////////////
+        //  First fundamental form
+        ////////////////////////////////
 
-    //////////////////////////////////////////////////////////
-    //  Fundamental forms and curvature
-    //////////////////////////////////////////////////////////
+        Dune::FieldMatrix<double,3,3> aCovariant;
 
-    // First fundamental form
-    Dune::FieldMatrix<double,3,3> aCovariant;
+        // If dimworld==3, then the first two lines of aCovariant are simply the jacobianTransposed
+        // of the element.  If dimworld<3 (i.e., ==2), we have to explicitly enters 0.0 in the last column.
+        auto jacobianTransposed = geometry.jacobianTransposed(quadPos);
 
-    // If dimworld==3, then the first two lines of aCovariant are simply the jacobianTransposed
-    // of the element.  If dimworld<3 (i.e., ==2), we have to explicitly enters 0.0 in the last column.
-    auto jacobianTransposed = geometry.jacobianTransposed(quadPos);
+        for (int i=0; i<2; i++)
+        {
+          for (int j=0; j<dimworld; j++)
+            aCovariant[i][j] = jacobianTransposed[i][j];
+          for (int j=dimworld; j<3; j++)
+            aCovariant[i][j] = 0.0;
+        }
 
-    for (int i=0; i<2; i++)
-    {
-      for (int j=0; j<dimworld; j++)
-        aCovariant[i][j] = jacobianTransposed[i][j];
-      for (int j=dimworld; j<3; j++)
-        aCovariant[i][j] = 0.0;
-    }
+        aCovariant[2] = Dune::MatrixVector::crossProduct(aCovariant[0], aCovariant[1]);
+        aCovariant[2] /= aCovariant[2].two_norm();
 
-    aCovariant[2] = Dune::MatrixVector::crossProduct(aCovariant[0], aCovariant[1]);
-    aCovariant[2] /= aCovariant[2].two_norm();
-
-    auto aContravariant = aCovariant;
-    aContravariant.invert();
-    // The contravariant base vectors are the *columns* of the inverse of the covariant matrix
-    // To get an easier access to the columns, we use the transpose of the contravariant matrix
-    aContravariant = Dune::GFE::transpose(aContravariant);
-
-    Dune::FieldMatrix<double,3,3> a(0);
-    for (int alpha=0; alpha<gridDim; alpha++)
-      a += Dune::GFE::dyadicProduct(aCovariant[alpha], aContravariant[alpha]);
-
-    auto a00 = aCovariant[0] * aCovariant[0];
-    auto a01 = aCovariant[0] * aCovariant[1];
-    auto a10 = aCovariant[1] * aCovariant[0];
-    auto a11 = aCovariant[1] * aCovariant[1];
-    auto aScalar = std::sqrt(a00*a11 - a10*a01);
-
-    // Alternator tensor
-    Dune::FieldMatrix<int,2,2> eps = {{0,1},{-1,0}};
-    Dune::FieldMatrix<double,3,3> c(0);
-
-    for (int alpha=0; alpha<2; alpha++)
-      for (int beta=0; beta<2; beta++)
-        c += aScalar * eps[alpha][beta] * Dune::GFE::dyadicProduct(aContravariant[alpha], aContravariant[beta]);
-
-    Dune::FieldMatrix<double,3,3> b(0);
 #if HAVE_DUNE_CURVEDGEOMETRY
-    // In case dune-curvedgeometry is not installed, we assume the second fundamental form b
-    // ( (-1) * the derivative of the normal field, evaluated at each quadrature point) is zero!
-    auto grad_s_n = geometry.normalGradient(quad[pt].position());
-    grad_s_n *= (-1);
-    for (int i = 0; i < dimworld; i++)
-      for (int j = 0; j < dimworld; j++)
-        b[i][j] = grad_s_n[i][j];
+        const auto normalGradient = geometry.normalGradient(quad[pt].position());
+#else
+        // Assume that the geometry is flat if DUNE_CURVEDGEOMETRY is not present.
+        // TODO: This is not always true!
+        Dune::FieldMatrix<double,3,3> normalGradient(0);
 #endif
 
-    // Mean curvatue
-    auto H = 0.5 * Dune::GFE::trace(b);
+        //////////////////////////////////////////////////////////
+        // Add the local energy density
+        //////////////////////////////////////////////////////////
 
-    // Gauss curvature, calculated with the normalGradient in the euclidean coordinate system
-    // see e.g. formula (3.5) in "Refined dimensional reduction for isotropic elastic Cosserat shells with initial curvature"
-    auto bSquared = b*b;
-    auto K = 2*H*H - 0.5*Dune::GFE::trace(bSquared);
+        const auto energyDensity = (*density_)(quadPos,
+                                               aCovariant,
+                                               normalGradient,
+                                               value,
+                                               derivative);
 
-    //////////////////////////////////////////////////////////
-    //  Strain tensors
-    //////////////////////////////////////////////////////////
-
-    // Elastic shell strain
-    Dune::FieldMatrix<field_type,3,3> Ee(0);
-    Dune::FieldMatrix<field_type,3,3> grad_s_m(0);
-    for (int alpha=0; alpha<gridDim; alpha++)
-    {
-      Dune::FieldVector<field_type,3> vec;
-      for (int i=0; i<3; i++)
-        vec[i] = derivative[i][alpha];
-      grad_s_m += Dune::GFE::dyadicProduct(vec, aContravariant[alpha]);
+        // Add energy density
+        energy += quad[pt].weight() * integrationElement * energyDensity;
+      }
     }
-
-    Ee = RT * grad_s_m - a;
-
-    // Elastic shell bending-curvature strain, e.g. formula (4.56) in "Refined dimensional reduction for isotropic elastic Cosserat shells with initial curvature"
-    Dune::FieldMatrix<field_type,3,3> Ke(0);
-    for (int alpha=0; alpha<gridDim; alpha++)
+    else
     {
-      Dune::FieldMatrix<field_type,3,3> tmp;
-      for (int i=0; i<3; i++)
-        for (int j=0; j<3; j++)
-          tmp[i][j] = DR[i][j][alpha];
-      auto tmp2 = RT * tmp;
-      Ke += Dune::GFE::dyadicProduct(SkewMatrix<field_type,3>(tmp2).axial(), aContravariant[alpha]);
+      // You need a scalar basis or a power basis when calling this method.
+      std::abort();
     }
 
-    //////////////////////////////////////////////////////////
-    // Add the local energy density
-    //////////////////////////////////////////////////////////
-
-    // Add the membrane energy density
-    field_type energyDensity = 0;
-    if (useAlternativeEnergyWCoss_) {
-      energyDensity += (thickness_ - K*Dune::power(thickness_,3) / 12.0) * W_Coss(Ee, a, aContravariant[2])
-                       + (Dune::power(thickness_,3) / 12.0 - K * Dune::power(thickness_,5) / 80.0)*W_Coss(Ee*b + c*Ke, a, aContravariant[2])
-                       + Dune::power(thickness_,3) / 6.0 * W_Coss_mixt(Ee, c*Ke*b - 2*H*c*Ke, a, aContravariant[2])
-                       + Dune::power(thickness_,5) / 80.0 * W_Coss( (Ee*b + c*Ke)*b, a, aContravariant[2]);
-    } else {
-      energyDensity += (thickness_ - K*Dune::power(thickness_,3) / 12.0) * W_m(Ee)
-                       + (Dune::power(thickness_,3) / 12.0 - K * Dune::power(thickness_,5) / 80.0)*W_m(Ee*b + c*Ke)
-                       + Dune::power(thickness_,3) / 6.0 * W_mixt(Ee, c*Ke*b - 2*H*c*Ke)
-                       + Dune::power(thickness_,5) / 80.0 * W_mp( (Ee*b + c*Ke)*b);
-    }
-
-    // Add the bending energy density
-    energyDensity += (thickness_ - K*Dune::power(thickness_,3) / 12.0) * W_curv(Ke)
-                     + (Dune::power(thickness_,3) / 12.0 - K * Dune::power(thickness_,5) / 80.0)*W_curv(Ke*b)
-                     + Dune::power(thickness_,5) / 80.0 * W_curv(Ke*b*b);
-
-    // Add energy density
-    energy += quad[pt].weight() * integrationElement * energyDensity;
-
-    ///////////////////////////////////////////////////////////
-    // Volume load contribution
-    ///////////////////////////////////////////////////////////
-
-    if (not volumeLoad_)
-      continue;
-
-    // Value of the volume load density at the current position
-    Dune::FieldVector<double,3> volumeLoadDensity = volumeLoad_(geometry.global(quad[pt].position()));
-
-    // Only translational dofs are affected by the volume load
-    for (size_t i=0; i<volumeLoadDensity.size(); i++)
-      energy += (volumeLoadDensity[i] * value[_0].globalCoordinates()[i]) * quad[pt].weight() * integrationElement;
-  }
-
-
-  //////////////////////////////////////////////////////////////////////////////
-  //   Assemble boundary contributions
-  //////////////////////////////////////////////////////////////////////////////
-
-  if (not neumannFunction_)
     return energy;
+  }
 
-  for (auto&& it : intersections(neumannBoundary_->gridView(),element) )
+  template <class Basis, class LocalGFEFunction, int dim, class field_type, class StressFreeStateGridFunction>
+  typename NonplanarCosseratShellEnergy<Basis, LocalGFEFunction, dim, field_type, StressFreeStateGridFunction>::RT
+  NonplanarCosseratShellEnergy<Basis,LocalGFEFunction,dim,field_type, StressFreeStateGridFunction>::
+  energy(const typename Basis::LocalView& localView,
+         const typename Dune::GFE::Impl::LocalEnergyTypes<TargetSpace>::CompositeCoefficients& localConfiguration) const
   {
-    if (not neumannBoundary_ or not neumannBoundary_->contains(it))
-      continue;
+    // The element geometry
+    auto element = localView.element();
 
-    const auto& quad = Dune::QuadratureRules<DT, gridDim-1>::rule(it.type(), quadOrder);
+    RT energy = 0;
 
-    for (size_t pt=0; pt<quad.size(); pt++)
+    if constexpr (Impl::LocalEnergyTypes<TargetSpace>::isProductManifold
+                  && Basis::LocalView::Tree::isComposite)
     {
-      // Local position of the quadrature point
-      const Dune::FieldVector<DT,gridDim>& quadPos = it.geometryInInside().global(quad[pt].position());
-
-      const DT integrationElement = it.geometry().integrationElement(quad[pt].position());
-
-      // The value of the local function
-      Dune::GFE::ProductManifold<RealTuple<field_type,dim>,Rotation<field_type,dim> > value = localGeodesicFEFunction.evaluate(quadPos);
-
-      // Value of the Neumann data at the current position
-      Dune::FieldVector<double,3> neumannValue = neumannFunction_(it.geometry().global(quad[pt].position()));
 
-      // Only translational dofs are affected by the Neumann force
-      for (size_t i=0; i<neumannValue.size(); i++)
-        energy += (neumannValue[i] * value[_0].globalCoordinates()[i]) * quad[pt].weight() * integrationElement;
-    }
-  }
+      // The set of shape functions on this element
 
-  return energy;
-}
-
-template <class Basis, int dim, class field_type, class StressFreeStateGridFunction>
-typename NonplanarCosseratShellEnergy<Basis, dim, field_type, StressFreeStateGridFunction>::RT
-NonplanarCosseratShellEnergy<Basis,dim,field_type, StressFreeStateGridFunction>::
-energy(const typename Basis::LocalView& localView,
-       const typename Dune::GFE::Impl::LocalEnergyTypes<TargetSpace>::CompositeCoefficients& localConfiguration) const
-{
-  // The element geometry
-  auto element = localView.element();
-
-  // The set of shape functions on this element
-
-  using namespace Dune::Indices;
-  const auto& deformationLocalFiniteElement = LocalFiniteElementFactory<Basis,0>::get(localView,_0);
-  const auto& orientationLocalFiniteElement = LocalFiniteElementFactory<Basis,1>::get(localView,_1);
+      using namespace Dune::Indices;
 
 #if HAVE_DUNE_CURVEDGEOMETRY
-  // Construct a curved geometry of this element of the Cosserat shell in its stress-free state
-  // The variable local holds the local coordinates in the reference element
-  // and localGeometry.global maps them to the world coordinates
-  Dune::CurvedGeometry<DT, gridDim, dimworld, Dune::CurvedGeometryTraits<DT, Dune::LagrangeLFECache<DT,DT,gridDim> > > geometry(referenceElement(element),
-                                                                                                                                [this,element](const auto& local) {
-                                                                                                                                auto localGridFunction = localFunction(*stressFreeStateGridFunction_);
-                                                                                                                                localGridFunction.bind(element);
-                                                                                                                                return localGridFunction(local);
-    }, getOrder(stressFreeStateGridFunction_));
+      // Construct a curved geometry of this element of the Cosserat shell in its stress-free state
+      // The variable local holds the local coordinates in the reference element
+      // and localGeometry.global maps them to the world coordinates
+      Dune::CurvedGeometry<DT, gridDim, dimworld, Dune::CurvedGeometryTraits<DT, Dune::LagrangeLFECache<DT,DT,gridDim> > > geometry(referenceElement(element),
+                                                                                                                                    [this,element](const auto& local) {
+                                                                                                                                    auto localGridFunction = localFunction(*stressFreeStateGridFunction_);
+                                                                                                                                    localGridFunction.bind(element);
+                                                                                                                                    return localGridFunction(local);
+        }, getOrder(stressFreeStateGridFunction_));
 #else
-  // When using element.geometry(), the geometry of the element is flat
-  auto geometry = element.geometry();
+      // When using element.geometry(), the geometry of the element is flat
+      auto geometry = element.geometry();
 #endif
 
-  ////////////////////////////////////////////////////////////////////////////////////
-  //  Set up the local nonlinear finite element function
-  ////////////////////////////////////////////////////////////////////////////////////
-  typedef LocalGeodesicFEFunction<gridDim, DT, decltype(deformationLocalFiniteElement), RealTuple<field_type,dim> > LocalDeformationGFEFunctionType;
-  typedef LocalGeodesicFEFunction<gridDim, DT, decltype(orientationLocalFiniteElement), Rotation<field_type,dim> > LocalOrientationGFEFunctionType;
-  LocalDeformationGFEFunctionType localDeformationGFEFunction(deformationLocalFiniteElement,localConfiguration[_0]);
-  LocalOrientationGFEFunctionType localOrientationGFEFunction(orientationLocalFiniteElement,localConfiguration[_1]);
-
-  RT energy = 0;
-
-  auto quadOrder = (deformationLocalFiniteElement.type().isSimplex()) ? deformationLocalFiniteElement.localBasis().order()
-                                                : deformationLocalFiniteElement.localBasis().order() * gridDim;
-
-  const auto& quad = Dune::QuadratureRules<DT, gridDim>::rule(element.type(), quadOrder);
-
-  for (size_t pt=0; pt<quad.size(); pt++)
-  {
-    // Local position of the quadrature point
-    const Dune::FieldVector<DT,gridDim>& quadPos = quad[pt].position();
-
-    const DT integrationElement = geometry.integrationElement(quadPos);
+      ////////////////////////////////////////////////////////////////////////////////////
+      //  Set up the local nonlinear finite element function
+      ////////////////////////////////////////////////////////////////////////////////////
+      std::get<0>(*localGFEFunction_).bind(element,localConfiguration[_0]);
+      std::get<1>(*localGFEFunction_).bind(element,localConfiguration[_1]);
 
-    // The value of the local function
-    RealTuple<field_type,dim> deformationValue = localDeformationGFEFunction.evaluate(quadPos);
-    Rotation<field_type,dim>  orientationValue = localOrientationGFEFunction.evaluate(quadPos);
+      // Bind the density to the current element
+      density_->bind(element);
 
-    // The derivative of the local function w.r.t. the coordinate system of the tangent space
-    typename LocalDeformationGFEFunctionType::DerivativeType deformationDerivative = localDeformationGFEFunction.evaluateDerivative(quadPos,deformationValue);
-    typename LocalOrientationGFEFunctionType::DerivativeType orientationDerivative = localOrientationGFEFunction.evaluateDerivative(quadPos,orientationValue);
+      const auto deformationGFEOrder = std::get<0>(*localGFEFunction_).localFiniteElement().localBasis().order();
 
-    //////////////////////////////////////////////////////////
-    //  The rotation and its derivative
-    //  Note: we need it in matrix coordinates
-    //////////////////////////////////////////////////////////
+      auto quadOrder = (element.type().isSimplex()) ? deformationGFEOrder
+                                                : deformationGFEOrder * gridDim;
 
-    Dune::FieldMatrix<field_type,dim,dim> R;
-    orientationValue.matrix(R);
-    auto RT = Dune::GFE::transpose(R);
+      const auto& quad = Dune::QuadratureRules<DT, gridDim>::rule(element.type(), quadOrder);
 
-    //Derivative of the rotation w.r.t. the coordinate system of the tangent space
-    Tensor3<field_type,3,3,gridDim> DR = orientationValue.quaternionTangentToMatrixTangent(orientationDerivative);
+      for (size_t pt=0; pt<quad.size(); pt++)
+      {
+        // Local position of the quadrature point
+        const Dune::FieldVector<DT,gridDim>& quadPos = quad[pt].position();
 
-    //////////////////////////////////////////////////////////
-    //  Fundamental forms and curvature
-    //////////////////////////////////////////////////////////
+        const DT integrationElement = geometry.integrationElement(quadPos);
 
-    // First fundamental form
-    Dune::FieldMatrix<double,3,3> aCovariant;
+        // The value of the local function
+        TargetSpace value;
+        value[_0] = std::get<0>(*localGFEFunction_).evaluate(quadPos);
+        value[_1] = std::get<1>(*localGFEFunction_).evaluate(quadPos);
 
-    // If dimworld==3, then the first two lines of aCovariant are simply the jacobianTransposed
-    // of the element.  If dimworld<3 (i.e., ==2), we have to explicitly enters 0.0 in the last column.
-    auto jacobianTransposed = geometry.jacobianTransposed(quadPos);
-    for (int i=0; i<2; i++)
-    {
-      for (int j=0; j<dimworld; j++)
-        aCovariant[i][j] = jacobianTransposed[i][j];
-      for (int j=dimworld; j<3; j++)
-        aCovariant[i][j] = 0.0;
-    }
+        // The derivative of the local function w.r.t. the coordinate system of the tangent space
+        const auto deformationDerivative = std::get<0>(*localGFEFunction_).evaluateDerivative(quadPos,value[_0]);
+        const auto orientationDerivative = std::get<1>(*localGFEFunction_).evaluateDerivative(quadPos,value[_1]);
 
-    aCovariant[2] = Dune::MatrixVector::crossProduct(aCovariant[0], aCovariant[1]);
-    aCovariant[2] /= aCovariant[2].two_norm();
+        // Concatenate the two derivative matrices
+        Dune::FieldMatrix<RT, TargetSpace::embeddedDim, gridDim> derivative;
 
+        for (int i=0; i<deformationDerivative.rows; ++i)
+          derivative[i] = deformationDerivative[i];
 
-    auto aContravariant = aCovariant;
-    aContravariant.invert();
-    // The contravariant base vectors are the *columns* of the inverse of the covariant matrix
-    // To get an easier access to the columns, we use the transpose of the contravariant matrix
-    aContravariant = Dune::GFE::transpose(aContravariant);
+        for (int i=0; i<orientationDerivative.rows; ++i)
+          derivative[i+deformationDerivative.rows] = orientationDerivative[i];
 
-    Dune::FieldMatrix<double,3,3> a(0);
-    for (int alpha=0; alpha<gridDim; alpha++)
-      a += Dune::GFE::dyadicProduct(aCovariant[alpha], aContravariant[alpha]);
+        //////////////////////////////////////////////////////////
+        //  Fundamental forms and curvature
+        //////////////////////////////////////////////////////////
 
-    auto a00 = aCovariant[0] * aCovariant[0];
-    auto a01 = aCovariant[0] * aCovariant[1];
-    auto a10 = aCovariant[1] * aCovariant[0];
-    auto a11 = aCovariant[1] * aCovariant[1];
-    auto aScalar = std::sqrt(a00*a11 - a10*a01);
+        // First fundamental form
+        Dune::FieldMatrix<double,3,3> aCovariant;
 
-    // Alternator tensor
-    Dune::FieldMatrix<int,2,2> eps = {{0,1},{-1,0}};
-    Dune::FieldMatrix<double,3,3> c(0);
+        // If dimworld==3, then the first two lines of aCovariant are simply the jacobianTransposed
+        // of the element.  If dimworld<3 (i.e., ==2), we have to explicitly enters 0.0 in the last column.
+        auto jacobianTransposed = geometry.jacobianTransposed(quadPos);
+        for (int i=0; i<2; i++)
+        {
+          for (int j=0; j<dimworld; j++)
+            aCovariant[i][j] = jacobianTransposed[i][j];
+          for (int j=dimworld; j<3; j++)
+            aCovariant[i][j] = 0.0;
+        }
 
-    for (int alpha=0; alpha<2; alpha++)
-      for (int beta=0; beta<2; beta++)
-        c += aScalar * eps[alpha][beta] * Dune::GFE::dyadicProduct(aContravariant[alpha], aContravariant[beta]);
+        aCovariant[2] = Dune::MatrixVector::crossProduct(aCovariant[0], aCovariant[1]);
+        aCovariant[2] /= aCovariant[2].two_norm();
 
-    Dune::FieldMatrix<double,3,3> b(0);
 #if HAVE_DUNE_CURVEDGEOMETRY
-    // In case dune-curvedgeometry is not installed, we assume the second fundamental form b
-    // ( (-1) * the derivative of the normal field, evaluated at each quadrature point) is zero!
-    auto grad_s_n = geometry.normalGradient(quad[pt].position());
-    for (int i = 0; i < dimworld; i++)
-      for (int j = 0; j < dimworld; j++)
-        b[i][j] = grad_s_n[i][j];
+        const auto normalGradient = geometry.normalGradient(quad[pt].position());
+#else
+        // Assume that the geometry is flat if DUNE_CURVEDGEOMETRY is not present.
+        // TODO: This is not always true!
+        Dune::FieldMatrix<double,3,3> normalGradient(0);
 #endif
-    b *= (-1);
 
-    // Mean curvatue
-    auto H = 0.5 * Dune::GFE::trace(b);
+        //////////////////////////////////////////////////////////
+        // Add the local energy density
+        //////////////////////////////////////////////////////////
 
-    // Gauss curvature, calculated with the normalGradient in the euclidean coordinate system
-    // see e.g. formula (3.5) in "Refined dimensional reduction for isotropic elastic Cosserat shells with initial curvature"
-    auto bSquared = b*b;
-    auto K = 2*H*H - 0.5*Dune::GFE::trace(bSquared);
+        const auto energyDensity = (*density_)(quadPos,
+                                               aCovariant,
+                                               normalGradient,
+                                               value,
+                                               derivative);
 
-    //////////////////////////////////////////////////////////
-    //  Strain tensors
-    //////////////////////////////////////////////////////////
+        // Add energy density
+        energy += quad[pt].weight() * integrationElement * energyDensity;
+      }
 
-    // Elastic shell strain
-    Dune::FieldMatrix<field_type,3,3> Ee(0);
-    Dune::FieldMatrix<field_type,3,3> grad_s_m(0);
-    for (int alpha=0; alpha<gridDim; alpha++)
-    {
-      Dune::FieldVector<field_type,3> vec;
-      for (int i=0; i<3; i++)
-        vec[i] = deformationDerivative[i][alpha];
-      grad_s_m += Dune::GFE::dyadicProduct(vec, aContravariant[alpha]);
     }
+    else
+      DUNE_THROW(Dune::NotImplemented, "Non-product manifold or non-composite basis");
 
-    Ee = RT * grad_s_m - a;
-
-    // Elastic shell bending-curvature strain, e.g. formula (4.56) in "Refined dimensional reduction for isotropic elastic Cosserat shells with initial curvature"
-    Dune::FieldMatrix<field_type,3,3> Ke(0);
-    for (int alpha=0; alpha<gridDim; alpha++)
-    {
-      Dune::FieldMatrix<field_type,3,3> tmp;
-      for (int i=0; i<3; i++)
-        for (int j=0; j<3; j++)
-          tmp[i][j] = DR[i][j][alpha];
-      auto tmp2 = RT * tmp;
-      Ke += Dune::GFE::dyadicProduct(SkewMatrix<field_type,3>(tmp2).axial(), aContravariant[alpha]);
-    }
-
-    //////////////////////////////////////////////////////////
-    // Add the local energy density
-    //////////////////////////////////////////////////////////
-
-    // Add the membrane energy density
-    field_type energyDensity = 0;
-    if (useAlternativeEnergyWCoss_) {
-      energyDensity += (thickness_ - K*Dune::power(thickness_,3) / 12.0) * W_Coss(Ee, a, aContravariant[2])
-                       + (Dune::power(thickness_,3) / 12.0 - K * Dune::power(thickness_,5) / 80.0)*W_Coss(Ee*b + c*Ke, a, aContravariant[2])
-                       + Dune::power(thickness_,3) / 6.0 * W_Coss_mixt(Ee, c*Ke*b - 2*H*c*Ke, a, aContravariant[2])
-                       + Dune::power(thickness_,5) / 80.0 * W_Coss( (Ee*b + c*Ke)*b, a, aContravariant[2]);
-    } else {
-      energyDensity += (thickness_ - K*Dune::power(thickness_,3) / 12.0) * W_m(Ee)
-                       + (Dune::power(thickness_,3) / 12.0 - K * Dune::power(thickness_,5) / 80.0)*W_m(Ee*b + c*Ke)
-                       + Dune::power(thickness_,3) / 6.0 * W_mixt(Ee, c*Ke*b - 2*H*c*Ke)
-                       + Dune::power(thickness_,5) / 80.0 * W_mp( (Ee*b + c*Ke)*b);
-    }
-
-    energyDensity += (thickness_ - K*Dune::power(thickness_,3) / 12.0) * W_curv(Ke)
-                     + (Dune::power(thickness_,3) / 12.0 - K * Dune::power(thickness_,5) / 80.0)*W_curv(Ke*b)
-                     + Dune::power(thickness_,5) / 80.0 * W_curv(Ke*b*b);
-
-    // Add energy density
-    energy += quad[pt].weight() * integrationElement * energyDensity;
-
-    ///////////////////////////////////////////////////////////
-    // Volume load contribution
-    ///////////////////////////////////////////////////////////
-
-    if (not volumeLoad_)
-      continue;
-
-    // Value of the volume load density at the current position
-    Dune::FieldVector<double,3> volumeLoadDensity = volumeLoad_(geometry.global(quad[pt].position()));
-
-    // Only translational dofs are affected by the volume load
-    for (size_t i=0; i<volumeLoadDensity.size(); i++)
-      energy += (volumeLoadDensity[i] * deformationValue[i]) * quad[pt].weight() * integrationElement;
-  }
-
-
-  //////////////////////////////////////////////////////////////////////////////
-  //   Assemble boundary contributions
-  //////////////////////////////////////////////////////////////////////////////
-
-  if (not neumannFunction_)
     return energy;
-
-  for (auto&& it : intersections(neumannBoundary_->gridView(),element) )
-  {
-    if (not neumannBoundary_ or not neumannBoundary_->contains(it))
-      continue;
-
-    const auto& quad = Dune::QuadratureRules<DT, gridDim-1>::rule(it.type(), quadOrder);
-
-    for (size_t pt=0; pt<quad.size(); pt++)
-    {
-      // Local position of the quadrature point
-      const Dune::FieldVector<DT,gridDim>& quadPos = it.geometryInInside().global(quad[pt].position());
-
-      const DT integrationElement = it.geometry().integrationElement(quad[pt].position());
-
-      // The value of the local function
-      RealTuple<field_type,dim> deformationValue = localDeformationGFEFunction.evaluate(quadPos);
-
-      // Value of the Neumann data at the current position
-      Dune::FieldVector<double,3> neumannValue = neumannFunction_(it.geometry().global(quad[pt].position()));
-
-      // Only translational dofs are affected by the Neumann force
-      for (size_t i=0; i<neumannValue.size(); i++)
-        energy += (neumannValue[i] * deformationValue[i]) * quad[pt].weight() * integrationElement;
-    }
   }
 
-  return energy;
-}
+}  // namespace Dune::GFE
 
 #endif   //#ifndef DUNE_GFE_NONPLANARCOSSERATSHELLENERGY_HH

dune/gfe/assemblers/simofoxenergy.hh

@@ -5,39 +5,15 @@
 #include <dune/common/parametertree.hh>
 #include <dune/common/transpose.hh>
 #include <dune/common/tuplevector.hh>
-#include <dune/fufem/boundarypatch.hh>
+#include <dune/functions/functionspacebases/subspacebasis.hh>
 #include <dune/geometry/quadraturerules.hh>
 #include <dune/gfe/linearalgebra.hh>
 #include <dune/gfe/assemblers/localenergy.hh>
 #include <dune/gfe/spaces/realtuple.hh>
 #include <dune/gfe/spaces/unitvector.hh>
 
-namespace Dune::GFE {
-
-  /** \brief Get LocalFiniteElements from a localView, for different tree depths of the local view
-   *
-   * Copied from CosseratEnergyLocalStiffness.hh
-   */
-  template <class Basis, std::size_t i>
-  class LocalFiniteElementFactory {
-  public:
-    static auto get(const typename Basis::LocalView &localView, std::integral_constant<std::size_t, i> iType)
-    -> decltype(localView.tree().child(iType, 0).finiteElement()) {
-      return localView.tree().child(iType, 0).finiteElement();
-    }
-  };
-
-  /** \brief Specialize for scalar bases, here we cannot call tree().child() */
-  template <class GridView, int order, std::size_t i>
-  class LocalFiniteElementFactory<Dune::Functions::LagrangeBasis<GridView, order>, i> {
-  public:
-    static auto get(const typename Dune::Functions::LagrangeBasis<GridView, order>::LocalView &localView,
-                    std::integral_constant<std::size_t, i> iType) -> decltype(localView.tree().finiteElement()) {
-      return localView.tree().finiteElement();
-    }
-  };
-
-
+namespace Dune::GFE
+{
   /** \brief Implements the energy of a Simo-Fox shell
    *
    * Described in:
@@ -47,8 +23,8 @@ namespace Dune::GFE {
    * \tparam LocalFEFunction Decides which interpolation function is used for the interpolation of the midsurface
    * and director field.
    */
-  template <class Basis, template <int, typename, typename, typename> typename LocalFEFunction, typename field_type = double>
-  class SimoFoxEnergyLocalStiffness
+  template <class Basis, typename LocalFEFunction, typename ReferenceLocalGFEFunction, typename field_type = double>
+  class SimoFoxEnergy
     : public LocalEnergy<Basis, ProductManifold<RealTuple<field_type, 3>,
           UnitVector<field_type, 3> > >
   {
@@ -58,43 +34,29 @@ namespace Dune::GFE {
     typedef typename Basis::GridView GridView;
     typedef typename GridView::ctype DT;
     typedef field_type RT;
-    typedef typename GridView::template Codim<0>::Entity Entity;
 
     // some other sizes
     static constexpr int gridDim  = GridView::dimension;
     static constexpr int dimworld = GridView::dimensionworld;
 
-    // The local finite element type used for midsurface position and displacement interpolation
-    using MidSurfaceElement =
-      typename std::result_of<decltype (&LocalFiniteElementFactory<Basis, 0>::get)(typename Basis::LocalView, decltype(Dune::Indices::_0))>::type;
-    // The local finite element type used for the director interpolation
-    using DirectorElement =
-      typename std::result_of<decltype (&LocalFiniteElementFactory<Basis, 1>::get)(typename Basis::LocalView, decltype(Dune::Indices::_1))>::type;
-
-    // The local finite element function type to evaluate the midsurface position
-    using LocalMidSurfaceFunctionType = LocalFEFunction<gridDim, DT, MidSurfaceElement, RealTuple<field_type, 3> >;
-    // The local finite element function type to evaluate the director
-    using LocalDirectorFunctionType = LocalFEFunction<gridDim, DT, DirectorElement, UnitVector<field_type, 3> >;
-
-    // Extra function type for reference quantities since they are unconditionally doubles, i.e. no ADOL-C types
-    using LocalMidSurfaceReferenceFunctionType = LocalFEFunction<gridDim, DT, MidSurfaceElement, RealTuple<double, 3> >;
-    using LocalDirectorReferenceFunctionType   = LocalFEFunction<gridDim, DT, DirectorElement, UnitVector<double, 3> >;
-
   public:
     /** \brief Constructor with a set of material parameters
      * \param parameters The material parameters
-     * \param x0 reference configuration
+     * \param localGFEFunction The geometric FE function that describes the current configuration
+     * \param referenceLocalGFEFunction The GFE function that describes
+     *   the stress-free configuration
+     * \param x0 Coefficients of the stress-free configuration
      */
-    SimoFoxEnergyLocalStiffness(const Dune::ParameterTree &parameters, const BoundaryPatch<GridView> *neumannBoundary,
-                                const std::function<Dune::FieldVector<double, 3>(Dune::FieldVector<double, 2>)> neumannFunction,
-                                const std::function<Dune::FieldVector<double, 3>(Dune::FieldVector<double, 2>)> volumeLoad,
-                                const Dune::TupleVector<std::vector<RealTuple<double, 3> >, std::vector<UnitVector<double, 3> > > &x0)
-      : neumannBoundary_(neumannBoundary),
-      neumannFunction_(neumannFunction),
-      thickness_{parameters.template get<double>("thickness")},      // The sheeqll thickness
+    SimoFoxEnergy(const Dune::ParameterTree &parameters,
+                  LocalFEFunction &&localGFEFunction,
+                  ReferenceLocalGFEFunction &&referenceLocalGFEFunction,
+                  const Dune::TupleVector<std::vector<RealTuple<double, 3> >, std::vector<UnitVector<double, 3> > > &x0)
+      : thickness_{parameters.template get<double>("thickness")},      // The sheeqll thickness
       mu_{parameters.template get<double>("mu")},                    // Lame constant 1
       lambda_{parameters.template get<double>("lambda")},            // Lame constant 2
       kappa_{parameters.template get<double>("kappa")},              // Shear correction factor
+      localGFEFunction_(std::make_shared<LocalFEFunction>(std::move(localGFEFunction))),
+      referenceLocalGFEFunction_(std::make_shared<ReferenceLocalGFEFunction>(std::move(referenceLocalGFEFunction))),
       midSurfaceRefConfig{x0[Dune::Indices::_0]},
       directorRefConfig{x0[Dune::Indices::_1]}
     {
@@ -146,27 +108,9 @@ namespace Dune::GFE {
 
     auto getReferenceLocalConfigurations(const typename Basis::LocalView &localView) const;
 
-    /** \brief Calculates all kinematic quantities */
-    template <typename Element, typename LocalDirectorFunction, typename LocalMidSurfaceFunction, typename LocalDirectorReferenceFunction,
-        typename LocalMidSurfaceReferenceFunction, typename IntegrationPointPosition>
-    static auto kinematicVariablesFactory(const Element &element, const LocalDirectorFunction &directorFunction,
-                                          const LocalDirectorReferenceFunction &directorReferenceFunction,
-                                          const LocalMidSurfaceFunction &midSurfaceFunction,
-                                          const LocalMidSurfaceReferenceFunction &midSurfaceReferenceFunction,
-                                          const LocalMidSurfaceFunction &midSurfaceDisplacementFunction, const IntegrationPointPosition &quadPos);
-
     /** \brief Calculate the Green-Lagrange strain components */
     static Dune::FieldVector<RT, 8> calculateGreenLagrangianStrains(const KinematicVariables &kin);
 
-    /** \brief Save all tangent base matrices for all nodes in one place*/
-    Dune::BlockVector<Dune::FieldMatrix<field_type, 2, 3> > directorTangentSpaces;
-
-    /** \brief The Neumann boundary */
-    const BoundaryPatch<GridView> *neumannBoundary_;
-
-    /** \brief The function implementing the Neumann data */
-    const std::function<Dune::FieldVector<double, 3>(Dune::FieldVector<double, dimworld>)> neumannFunction_;
-
     /** \brief The shell thickness */
     double thickness_;
 
@@ -179,10 +123,10 @@ namespace Dune::GFE {
     /** \brief Material tangent matrix */
     Dune::FieldMatrix<double, 8, 8> CMat_;
 
-    /** \brief The function implementing a volume load */
-    const std::function<Dune::FieldVector<double, 3>(Dune::FieldVector<double, dimworld>)> volumeLoad_;
+    const std::shared_ptr<LocalFEFunction> localGFEFunction_;
 
     /** \brief Stores the reference configuration of the midsurface and the director field */
+    const std::shared_ptr<ReferenceLocalGFEFunction> referenceLocalGFEFunction_;
     const std::vector<RealTuple<double, 3> > &midSurfaceRefConfig;
     const std::vector<UnitVector<double, 3> > &directorRefConfig;
   };
@@ -194,8 +138,8 @@ namespace Dune::GFE {
    * b is similar to the second fundamental form of the surface, but the unit normal is replaced with the shell director
    * gamma_1 and gamma_2 are the transverse shear in the two parametric directions
    * Paper Equation 4.10 */
-  template <class Basis, template <int, typename, typename, typename> typename LocalFEFunction, typename field_type>
-  Dune::FieldVector<field_type, 8> SimoFoxEnergyLocalStiffness<Basis, LocalFEFunction, field_type>::calculateGreenLagrangianStrains(
+  template <class Basis, typename LocalFEFunction, typename ReferenceLocalGFEFunction, typename field_type>
+  Dune::FieldVector<field_type, 8> SimoFoxEnergy<Basis, LocalFEFunction, ReferenceLocalGFEFunction, field_type>::calculateGreenLagrangianStrains(
     const KinematicVariables &kin)
   {
     Dune::FieldVector<RT, 8> egl;
@@ -217,37 +161,10 @@ namespace Dune::GFE {
     return egl;
   }
 
-  /** \brief Calculates all kinematic quantities that are needed for strain calculation
-   *
-   * Paper Equations 6.4 - 6.8
-   */
-  template <class Basis, template <int, typename, typename, typename> typename LocalFEFunction, typename field_type>
-  template <typename Element, typename LocalDirectorFunction, typename LocalMidSurfaceFunction, typename LocalDirectorReferenceFunction,
-      typename LocalMidSurfaceReferenceFunction, typename IntegrationPointPosition>
-  auto SimoFoxEnergyLocalStiffness<Basis, LocalFEFunction, field_type>::kinematicVariablesFactory(
-    const Element &element, const LocalDirectorFunction &directorFunction, const LocalDirectorReferenceFunction &directorReferenceFunction,
-    const LocalMidSurfaceFunction &midSurfaceFunction, const LocalMidSurfaceReferenceFunction &midSurfaceReferenceFunction,
-    const LocalMidSurfaceFunction &midSurfaceDisplacementFunction, const IntegrationPointPosition &quadPos)
-  {
-    KinematicVariables kin{};
-
-    const auto jInvT = element.geometry().jacobianInverseTransposed(quadPos);
-
-    kin.t           = directorFunction.evaluate(quadPos).globalCoordinates();
-    kin.t0          = directorReferenceFunction.evaluate(quadPos).globalCoordinates();
-    kin.a1anda2     = transpose(midSurfaceFunction.evaluateDerivative(quadPos) * transpose(jInvT));
-    kin.A1andA2     = transpose(midSurfaceReferenceFunction.evaluateDerivative(quadPos) * transpose(jInvT));
-    kin.ud1andud2   = transpose(midSurfaceDisplacementFunction.evaluateDerivative(quadPos) * transpose(jInvT));
-    kin.t0d1Andt0d2 = transpose(directorReferenceFunction.evaluateDerivative(quadPos) * transpose(jInvT));
-    kin.td1Andtd2   = transpose(directorFunction.evaluateDerivative(quadPos) * transpose(jInvT));
-
-    return kin;
-  }
-
-  template <class Basis, template <int, typename, typename, typename> typename LocalFEFunction, typename field_type>
-  auto SimoFoxEnergyLocalStiffness<Basis, LocalFEFunction, field_type>::getReferenceLocalConfigurations(
+  template <class Basis, typename LocalFEFunction, typename ReferenceLocalGFEFunction, typename field_type>
+  auto SimoFoxEnergy<Basis, LocalFEFunction, ReferenceLocalGFEFunction, field_type>::getReferenceLocalConfigurations(
     const typename Basis::LocalView &localView) const {
-    using namespace Dune::TypeTree::Indices;
+    using namespace Dune::Indices;
     const int nDofs0 = localView.tree().child(_0, 0).finiteElement().size();
     const int nDofs1 = localView.tree().child(_1, 0).finiteElement().size();
 
@@ -282,9 +199,9 @@ namespace Dune::GFE {
    *  E are the components of the Green-Lagrangian strains [membrane strains, bending, transverse shear]
    *  see for details Paper Equation 4.11,4.10 and 10.1
    */
-  template <class Basis, template <int, typename, typename, typename> typename LocalFEFunction, typename field_type>
-  typename SimoFoxEnergyLocalStiffness<Basis, LocalFEFunction, field_type>::RT
-  SimoFoxEnergyLocalStiffness<Basis, LocalFEFunction, field_type>::energy(
+  template <class Basis, typename LocalFEFunction, typename ReferenceLocalGFEFunction, typename field_type>
+  typename SimoFoxEnergy<Basis, LocalFEFunction, ReferenceLocalGFEFunction, field_type>::RT
+  SimoFoxEnergy<Basis, LocalFEFunction, ReferenceLocalGFEFunction, field_type>::energy(
     const typename Basis::LocalView &localView,
     const typename Impl::LocalEnergyTypes<TargetSpace>::CompositeCoefficients &localConfiguration) const
   {
@@ -294,9 +211,6 @@ namespace Dune::GFE {
 
     auto element = localView.element();
 
-    const auto &midSurfaceElement = LocalFiniteElementFactory<Basis, 0>::get(localView, _0);
-    const auto &directorElement   = LocalFiniteElementFactory<Basis, 1>::get(localView, _1);
-
     const auto [localRefMidSurfaceConfiguration, localRefDirectorConfiguration] = getReferenceLocalConfigurations(localView);
 
     std::vector<RealTuple<field_type, 3> > displacements;
@@ -304,14 +218,26 @@ namespace Dune::GFE {
     for (size_t i = 0; i < localMidSurfaceConfiguration.size(); ++i)
       displacements.emplace_back(localMidSurfaceConfiguration[i].globalCoordinates() - localRefMidSurfaceConfiguration[i].globalCoordinates());
 
-    const LocalMidSurfaceFunctionType localMidSurfaceFunction(midSurfaceElement, localMidSurfaceConfiguration);
-    const LocalMidSurfaceFunctionType localMidSurfaceDisplacementFunction(midSurfaceElement, displacements);
-    const LocalMidSurfaceReferenceFunctionType localMidSurfaceReferenceFunction(midSurfaceElement, localRefMidSurfaceConfiguration);
-    const LocalDirectorFunctionType localDirectorFunction(directorElement, localDirectorConfiguration);
-    const LocalDirectorReferenceFunctionType localDirectorReferenceFunction(directorElement, localRefDirectorConfiguration);
+    auto& localMidSurfaceFunction = std::get<0>(*localGFEFunction_);
+    auto& localDirectorFunction = std::get<1>(*localGFEFunction_);
+    auto& localMidSurfaceReferenceFunction = std::get<0>(*referenceLocalGFEFunction_);
+    auto& localDirectorReferenceFunction = std::get<1>(*referenceLocalGFEFunction_);
+
+    // Copy of the midsurface deformation function to hold the displacement
+    using LocalMidSurfaceFunctionType = std::tuple_element_t<0,LocalFEFunction>;
+    LocalMidSurfaceFunctionType localMidSurfaceDisplacementFunction(localMidSurfaceFunction);
 
-    const int quadOrder = (element.type().isSimplex()) ? std::max(midSurfaceElement.localBasis().order(), directorElement.localBasis().order())
-                                                       : std::max(midSurfaceElement.localBasis().order(), directorElement.localBasis().order()) + 1;
+    localMidSurfaceFunction.bind(element, localMidSurfaceConfiguration);
+    localMidSurfaceDisplacementFunction.bind(element, displacements);
+    localMidSurfaceReferenceFunction.bind(element, localRefMidSurfaceConfiguration);
+    localDirectorFunction.bind(element, localDirectorConfiguration);
+    localDirectorReferenceFunction.bind(element, localRefDirectorConfiguration);
+
+    const auto midSurfaceOrder = localMidSurfaceFunction.localFiniteElement().localBasis().order();
+    const auto directorOrder = localDirectorFunction.localFiniteElement().localBasis().order();
+    const auto gfeOrder = std::max(midSurfaceOrder, directorOrder);
+
+    const int quadOrder = (element.type().isSimplex()) ? gfeOrder : gfeOrder + 1;
 
     const auto &quad = Dune::QuadratureRules<DT, gridDim>::rule(element.type(), quadOrder);
 
@@ -319,46 +245,29 @@ namespace Dune::GFE {
     for (const auto &curQuad : quad) {
       const Dune::FieldVector<DT, gridDim> &quadPos = curQuad.position();
 
-      const KinematicVariables kin
-        = kinematicVariablesFactory(element, localDirectorFunction, localDirectorReferenceFunction, localMidSurfaceFunction,
-                                    localMidSurfaceReferenceFunction, localMidSurfaceDisplacementFunction, quadPos);
+      KinematicVariables kin{};
+
+      const auto geometryJacobianInverse = element.geometry().jacobianInverse(quadPos);
+
+      // Compute all kinematic quantities needed for strain calculation
+      //
+      // Paper Equations (6.4)--(6.8)
+      kin.t           = localDirectorFunction.evaluate(quadPos).globalCoordinates();
+      kin.t0          = localDirectorReferenceFunction.evaluate(quadPos).globalCoordinates();
+      kin.a1anda2     = transpose(localMidSurfaceFunction.evaluateDerivative(quadPos) * geometryJacobianInverse);
+      kin.A1andA2     = transpose(localMidSurfaceReferenceFunction.evaluateDerivative(quadPos) * geometryJacobianInverse);
+      kin.ud1andud2   = transpose(localMidSurfaceDisplacementFunction.evaluateDerivative(quadPos) * geometryJacobianInverse);
+      kin.t0d1Andt0d2 = transpose(localDirectorReferenceFunction.evaluateDerivative(quadPos) * geometryJacobianInverse);
+      kin.td1Andtd2   = transpose(localDirectorFunction.evaluateDerivative(quadPos) * geometryJacobianInverse);
 
       const DT integrationElement = element.geometry().integrationElement(quadPos);
       const FieldVector<field_type, 8> Egl = calculateGreenLagrangianStrains(kin);
       energy +=  0.5 * Egl * (CMat_ * Egl) * curQuad.weight() * integrationElement;
     }
 
-    //////////////////////////////////////////////////////////////////////////////
-    //   Assemble boundary contributions
-    //////////////////////////////////////////////////////////////////////////////
-
-    if (not neumannFunction_) return energy;
-
-    for (auto &&it : intersections(neumannBoundary_->gridView(), element))
-    {
-      if (not neumannBoundary_ or not neumannBoundary_->contains(it)) continue;
-
-      const auto &quadLine = QuadratureRules<DT, gridDim - 1>::rule(it.type(), quadOrder);
-
-      for (const auto &curQuad : quadLine)
-      {
-        // Local position of the quadrature point
-        const FieldVector<DT, gridDim> &quadPos = it.geometryInInside().global(curQuad.position());
-
-        const DT integrationElement = it.geometry().integrationElement(curQuad.position());
-
-        // The value of the local function
-        RealTuple<field_type, 3> deformationValue = localMidSurfaceFunction.evaluate(quadPos);
-
-        // Value of the Neumann data at the current position
-        auto neumannValue = neumannFunction_(it.geometry().global(curQuad.position()));
-
-        // Only translational dofs are affected by the Neumann force
-        for (size_t i = 0; i < neumannValue.size(); i++)
-          energy -= (neumannValue[i] * deformationValue.globalCoordinates()[i]) * curQuad.weight() * integrationElement;
-      }
-    }
     return energy;
   }
+
 }  // namespace Dune::GFE
+
 #endif  // DUNE_GFE_SIMOFOX_ENERGY_HH

dune/gfe/assemblers/sumenergy.hh

@@ -5,8 +5,8 @@
 
 #include <dune/gfe/assemblers/localenergy.hh>
 
-namespace Dune::GFE {
-
+namespace Dune::GFE
+{
   /**
      \brief Assembles the a sum of energies for a single element by summing up the energies of each GFE::LocalEnergy.
 

dune/gfe/assemblers/surfacecosseratenergy.hh

@@ -6,10 +6,12 @@
 
 #include <dune/fufem/boundarypatch.hh>
 
+#include <dune/functions/functionspacebases/subspacebasis.hh>
+
 #include <dune/gfe/assemblers/localenergy.hh>
-#include <dune/gfe/localgeodesicfefunction.hh>
-#include <dune/gfe/localprojectedfefunction.hh>
 #include <dune/gfe/densities/cosseratshelldensity.hh>
+#include <dune/gfe/functions/localgeodesicfefunction.hh>
+#include <dune/gfe/functions/localprojectedfefunction.hh>
 #include <dune/gfe/spaces/productmanifold.hh>
 
 #if HAVE_DUNE_CURVEDGEOMETRY
@@ -37,24 +39,22 @@ namespace Dune::GFE
     constexpr static int gridDim = GridView::dimension;
     static constexpr int boundaryDim = gridDim - 1;
 
+    using Intersection = typename GridView::Intersection;
+
   public:
 
     /** \brief Constructor with a set of material parameters
-     * \param parameters The material parameters
+     * \param density The density that is being integrated over
      * \param shellBoundary The shellBoundary contains the faces where the cosserat energy is assembled
      * \param curvedGeometryGridFunction The curvedGeometryGridFunction gives the geometry of the shell in stress-free state.
               When assembling, we deform the intersections using the curvedGeometryGridFunction and then use the deformed geometries.
-     * \param thicknessF The shell thickness parameter, given as a function and evaluated at each quadrature point
-     * \param lameF The Lame parameters, given as a function and evaluated at each quadrature point
      */
-    SurfaceCosseratEnergy(const Dune::ParameterTree& parameters,
+    SurfaceCosseratEnergy(const std::shared_ptr<CosseratShellDensity<Intersection,RT> >& density,
                           const BoundaryPatch<GridView>* shellBoundary,
-                          const CurvedGeometryGridFunction& curvedGeometryGridFunction,
-                          const std::function<double(Dune::FieldVector<double,dimWorld>)> thicknessF,
-                          const std::function<Dune::FieldVector<double,2>(Dune::FieldVector<double,dimWorld>)> lameF)
+                          const CurvedGeometryGridFunction& curvedGeometryGridFunction)
       : shellBoundary_(shellBoundary),
       curvedGeometryGridFunction_(curvedGeometryGridFunction),
-      density_(parameters,thicknessF,lameF)
+      density_(density)
     {}
 
     /** \brief Assemble the energy for a single element */
@@ -71,6 +71,8 @@ namespace Dune::GFE
 
       using namespace Dune::Indices;
 
+      const auto element = localView.element();
+
       ////////////////////////////////////////////////////////////////////////////////////
       //  Set up the local nonlinear finite element function
       ////////////////////////////////////////////////////////////////////////////////////
@@ -78,27 +80,29 @@ namespace Dune::GFE
       using RBM0 = typename std::tuple_element<0,TargetSpace>::type;
       using RBM1 = typename std::tuple_element<1,TargetSpace>::type;
 
-      // The set of shape functions on this element
-      const auto& deformationLocalFiniteElement = localView.tree().child(_0,0).finiteElement();
-      const auto& orientationLocalFiniteElement = localView.tree().child(_1,0).finiteElement();
+      auto deformationScalarBasis = Functions::subspaceBasis(localView.globalBasis(), _0, 0);
+      auto rotationScalarBasis = Functions::subspaceBasis(localView.globalBasis(), _1, 0);
 
-      typedef LocalGeodesicFEFunction<gridDim, DT, decltype(deformationLocalFiniteElement), RBM0> LocalGFEFunctionType0;
-      typedef LocalGeodesicFEFunction<gridDim, DT, decltype(orientationLocalFiniteElement), RBM1> LocalGFEFunctionType1;
-      LocalGFEFunctionType0 localGeodesicFEFunction0(deformationLocalFiniteElement,localCoefficients[_0]);
-      LocalGFEFunctionType1 localGeodesicFEFunction1(orientationLocalFiniteElement,localCoefficients[_1]);
+      typedef LocalGeodesicFEFunction<decltype(deformationScalarBasis), RBM0> LocalGFEFunctionType0;
+      typedef LocalGeodesicFEFunction<decltype(rotationScalarBasis), RBM1> LocalGFEFunctionType1;
+      LocalGFEFunctionType0 localGeodesicFEFunction0(deformationScalarBasis);
+      LocalGFEFunctionType1 localGeodesicFEFunction1(rotationScalarBasis);
+      localGeodesicFEFunction0.bind(element,localCoefficients[_0]);
+      localGeodesicFEFunction1.bind(element,localCoefficients[_1]);
 
       RT energy = 0;
 
-      const auto element = localView.element();
-
       for (auto&& it : intersections(shellBoundary_->gridView(), element))
       {
         if (not shellBoundary_->contains(it))
           continue;
 
+        // Bind the density to the current intersection
+        density_->bind(it);
+
+        const auto deformationOrder = localView.tree().child(_0,0).finiteElement().localBasis().order();
 #if HAVE_DUNE_CURVEDGEOMETRY
         auto localGridFunction = localFunction(curvedGeometryGridFunction_);
-        auto curvedGeometryGridFunctionOrder = deformationLocalFiniteElement.localBasis().order();
         localGridFunction.bind(element);
         auto referenceElement = Dune::referenceElement<DT,boundaryDim>(it.type());
 
@@ -110,13 +114,12 @@ namespace Dune::GFE
         BoundaryGeometry boundaryGeometry(referenceElement,
                                           [localGridFunction, localGeometry=it.geometryInInside()](const auto& local) {
                                           return localGridFunction(localGeometry.global(local));
-          }, curvedGeometryGridFunctionOrder);
+          }, deformationOrder);
 #else
         const auto& boundaryGeometry = it.geometry();
 #endif
 
-        auto quadOrder = (it.type().isSimplex()) ? deformationLocalFiniteElement.localBasis().order()
-                                                  : deformationLocalFiniteElement.localBasis().order() * boundaryDim;
+        auto quadOrder = (it.type().isSimplex()) ? deformationOrder : deformationOrder * boundaryDim;
 
         const auto& quad = Dune::QuadratureRules<DT, boundaryDim>::rule(it.type(), quadOrder);
 
@@ -125,12 +128,8 @@ namespace Dune::GFE
           // Local position of the quadrature point
           const auto& quadPos = it.geometryInInside().global(qp.position());
 
-          // Global position of the quadrature point
-          auto quadPosGlobal = it.geometry().global(qp.position());
-
           const DT integrationElement = boundaryGeometry.integrationElement(qp.position());
 
-          FieldMatrix<RT, TargetSpace::embeddedDim, dimWorld> derivative;
           FieldMatrix<RT, TargetSpace::embeddedDim, boundaryDim> derivative2D;
 
           // The value of the local functions
@@ -145,12 +144,10 @@ namespace Dune::GFE
           // Put the value and the derivatives together from the separated values
           for (int j = 0; j < dimWorld; j++) {
             for (int i = 0; i < RBM0::embeddedDim; i++) {
-              derivative[i][j] = derivative0[i][j];
               if (j < boundaryDim)
                 derivative2D[i][j] = derivative0[i][j];
             }
             for (int i = 0; i < RBM1::embeddedDim; i++) {
-              derivative[RBM0::embeddedDim + i][j] = derivative1[i][j];
               if (j < boundaryDim)
                 derivative2D[RBM0::embeddedDim + i][j] = derivative1[i][j];
             }
@@ -186,11 +183,11 @@ namespace Dune::GFE
           // Add the local energy density
           //////////////////////////////////////////////////////////
 
-          const auto energyDensity = density_(quadPosGlobal,
-                                              aCovariant,
-                                              normalGradient,
-                                              value,
-                                              derivative2D);
+          const auto energyDensity = (*density_)(qp.position(),
+                                                 aCovariant,
+                                                 normalGradient,
+                                                 value,
+                                                 derivative2D);
           energy += qp.weight() * integrationElement * energyDensity;
         }
       }
@@ -206,7 +203,7 @@ namespace Dune::GFE
     const CurvedGeometryGridFunction curvedGeometryGridFunction_;
 
     /** \brief The density that is being integrated */
-    const CosseratShellDensity<FieldVector<DT,3>,RT> density_;
+    const std::shared_ptr<CosseratShellDensity<Intersection,RT> > density_;
 
   };
 }  // namespace Dune::GFE

dune/gfe/assemblers/surfacecosseratstressassembler.hh

@@ -4,20 +4,22 @@
 #include <dune/fufem/boundarypatch.hh>
 
 #include <dune/gfe/linearalgebra.hh>
-#include <dune/gfe/localgeodesicfefunction.hh>
 #include <dune/gfe/spaces/rotation.hh>
 
 #include <dune/matrix-vector/transpose.hh>
 
-namespace Dune::GFE {
+namespace Dune::GFE
+{
+
   /** \brief An assembler that can calculate the norms of specific stress tensors for each element for an output by film-on-substrate
 
      \tparam BasisOrderD Basis used for the displacement
      \tparam BasisOrderR Basis used for the rotation
+     \tparam LocalGFEFunctionR Geometric FE function for the rotations
      \tparam TargetSpaceD Target space for the Displacement
      \tparam TargetSpaceR Target space for the Rotation
    */
-  template <class BasisOrderD, class BasisOrderR, class TargetSpaceD, class TargetSpaceR>
+  template <class BasisOrderD, class BasisOrderR, class LocalGFEFunctionR, class TargetSpaceD, class TargetSpaceR>
   class SurfaceCosseratStressAssembler
   {
   public:
@@ -179,6 +181,7 @@ namespace Dune::GFE {
        \param stressShellBiotTensor Vector containing the Biot-Stress-Tensor for each element
      */
     void assembleShellStress(
+      LocalGFEFunctionR& localGeodesicFEFunction,
       const VectorR rot,
       const VectorD x,
       const VectorD xInitial,
@@ -233,8 +236,8 @@ namespace Dune::GFE {
           VectorR localConfigurationRot(lFEOrderR.size());
           for (std::size_t i=0; i<localConfigurationRot.size(); i++)
             localConfigurationRot[i] = rot[localViewOrderR.index(i)[0]];  //localViewOrderR.index(i) is a multiindex, its first entry is the actual index
-          typedef LocalGeodesicFEFunction<dim, double, decltype(lFEOrderR), TargetSpaceR> LocalGFEFunctionType;
-          LocalGFEFunctionType localGeodesicFEFunction(lFEOrderR,localConfigurationRot);
+
+          localGeodesicFEFunction.bind(element,localConfigurationRot);
 
           auto evaluateAtPoint = [&](FieldVector<double,3> pointGlobal, FieldVector<double,3> pointLocal3d) -> FieldMatrix<double,dim,dim> {
                                    Dune::FieldMatrix<double,dim,dim> nablaTheta;
@@ -301,5 +304,7 @@ namespace Dune::GFE {
       }
     }
   };
-}
+
+}  // namespace Dune::GFE
+
 #endif

dune/gfe/assemblers/weightedsumenergy.hh

@@ -7,52 +7,57 @@
 #include <dune/geometry/quadraturerules.hh>
 
 #include <dune/gfe/assemblers/localenergy.hh>
-#include <dune/gfe/localgeodesicfefunction.hh>
 
-template<class Basis, class TargetSpace>
-class WeightedSumEnergy
-  : public Dune::GFE::LocalEnergy<Basis,TargetSpace>
+
+namespace Dune::GFE
 {
-  // grid types
-  typedef typename Basis::GridView GridView;
-  typedef typename GridView::ctype DT;
-  typedef typename TargetSpace::ctype RT;
 
-  // some other sizes
-  constexpr static int gridDim = GridView::dimension;
+  template<class Basis, class TargetSpace>
+  class WeightedSumEnergy
+    : public Dune::GFE::LocalEnergy<Basis,TargetSpace>
+  {
+    // grid types
+    typedef typename Basis::GridView GridView;
+    typedef typename GridView::ctype DT;
+    typedef typename TargetSpace::ctype RT;
 
-public:
+    // some other sizes
+    constexpr static int gridDim = GridView::dimension;
 
-  std::vector<std::shared_ptr<Dune::GFE::LocalEnergy<Basis,TargetSpace> > > addends_;
+  public:
 
-  std::vector<double> weights_;
+    std::vector<std::shared_ptr<Dune::GFE::LocalEnergy<Basis,TargetSpace> > > addends_;
 
-  WeightedSumEnergy(std::vector<std::shared_ptr<Dune::GFE::LocalEnergy<Basis,TargetSpace> > > addends,
-                    std::vector<double> weights)
-    : addends_(addends),
-    weights_(weights)
-  {}
+    std::vector<double> weights_;
 
-  /** \brief Assemble the energy for a single element */
-  RT energy (const typename Basis::LocalView& localView,
-             const std::vector<TargetSpace>& localConfiguration) const override
+    WeightedSumEnergy(std::vector<std::shared_ptr<Dune::GFE::LocalEnergy<Basis,TargetSpace> > > addends,
+                      std::vector<double> weights)
+      : addends_(addends),
+      weights_(weights)
+    {}
 
-  {
-    RT energy = 0;
+    /** \brief Assemble the energy for a single element */
+    RT energy (const typename Basis::LocalView& localView,
+               const std::vector<TargetSpace>& localConfiguration) const override
 
-    assert(weights_.size() == addends_.size());
+    {
+      RT energy = 0;
 
-    for (size_t i=0; i<addends_.size(); i++)
-      energy += weights_[i] * addends_[i]->energy(localView, localConfiguration);
+      assert(weights_.size() == addends_.size());
 
-    return energy;
-  }
+      for (size_t i=0; i<addends_.size(); i++)
+        energy += weights_[i] * addends_[i]->energy(localView, localConfiguration);
 
-  RT energy (const typename Basis::LocalView& localView,
-             const typename Dune::GFE::Impl::LocalEnergyTypes<TargetSpace>::CompositeCoefficients& coefficients) const override
-  {
-    DUNE_THROW(Dune::NotImplemented, "!");
-  }
-};
+      return energy;
+    }
+
+    RT energy (const typename Basis::LocalView& localView,
+               const typename Dune::GFE::Impl::LocalEnergyTypes<TargetSpace>::CompositeCoefficients& coefficients) const override
+    {
+      DUNE_THROW(Dune::NotImplemented, "!");
+    }
+  };
+
+}  // namespace Dune::GFE
 
 #endif

dune/gfe/averagedistanceassembler.hh

@@ -5,97 +5,103 @@
 
 #include <dune/gfe/symmetricmatrix.hh>
 
-/** \tparam TargetSpace The manifold that we are mapping to */
-template <class TargetSpace, class WeightType=double>
-class AverageDistanceAssembler
+
+namespace Dune::GFE
 {
-  typedef typename TargetSpace::ctype ctype;
-  static const int size         = TargetSpace::TangentVector::dimension;
-  static const int embeddedSize = TargetSpace::EmbeddedTangentVector::dimension;
-
-public:
-
-  /** \brief Constructor with given coefficients \f$ v_i \f$ and weights \f$ w_i \f$
-   */
-  AverageDistanceAssembler(const std::vector<TargetSpace>& coefficients,
-                           const std::vector<WeightType>& weights)
-    : coefficients_(coefficients),
-    weights_(weights)
-  {}
-
-  /** \brief Constructor with given coefficients \f$ v_i \f$ and weights \f$ w_i \f$
-   *
-   * The weights are given as a vector of length-1 FieldVectors instead of as a
-   * vector of doubles.  The reason is that these weights are actually the values
-   * of Lagrange shape functions, and the dune-localfunction interface returns
-   * shape function values this way.
-   */
-  AverageDistanceAssembler(const std::vector<TargetSpace>& coefficients,
-                           const std::vector<Dune::FieldVector<WeightType,1> >& weights)
-    : coefficients_(coefficients),
-    weights_(weights.size())
+
+  /** \tparam TargetSpace The manifold that we are mapping to */
+  template <class TargetSpace, class WeightType=double>
+  class AverageDistanceAssembler
   {
-    for (size_t i=0; i<weights.size(); i++)
-      weights_[i] = weights[i][0];
-  }
+    typedef typename TargetSpace::ctype ctype;
+    static const int size         = TargetSpace::TangentVector::dimension;
+    static const int embeddedSize = TargetSpace::EmbeddedTangentVector::dimension;
+
+  public:
+
+    /** \brief Constructor with given coefficients \f$ v_i \f$ and weights \f$ w_i \f$
+     */
+    AverageDistanceAssembler(const std::vector<TargetSpace>& coefficients,
+                             const std::vector<WeightType>& weights)
+      : coefficients_(coefficients),
+      weights_(weights)
+    {}
+
+    /** \brief Constructor with given coefficients \f$ v_i \f$ and weights \f$ w_i \f$
+     *
+     * The weights are given as a vector of length-1 FieldVectors instead of as a
+     * vector of doubles.  The reason is that these weights are actually the values
+     * of Lagrange shape functions, and the dune-localfunction interface returns
+     * shape function values this way.
+     */
+    AverageDistanceAssembler(const std::vector<TargetSpace>& coefficients,
+                             const std::vector<Dune::FieldVector<WeightType,1> >& weights)
+      : coefficients_(coefficients),
+      weights_(weights.size())
+    {
+      for (size_t i=0; i<weights.size(); i++)
+        weights_[i] = weights[i][0];
+    }
+
+    ctype value(const TargetSpace& x) const {
 
-  ctype value(const TargetSpace& x) const {
+      ctype result = 0;
+      for (size_t i=0; i<coefficients_.size(); i++) {
+        ctype dist = TargetSpace::distance(coefficients_[i], x);
+        result += weights_[i]*dist*dist;
+      }
 
-    ctype result = 0;
-    for (size_t i=0; i<coefficients_.size(); i++) {
-      ctype dist = TargetSpace::distance(coefficients_[i], x);
-      result += weights_[i]*dist*dist;
+      return result;
     }
 
-    return result;
-  }
+    void assembleEmbeddedGradient(const TargetSpace& x,
+                                  typename TargetSpace::EmbeddedTangentVector& gradient) const
+    {
+      gradient = 0;
+      for (size_t i=0; i<coefficients_.size(); i++)
+        gradient.axpy(weights_[i],
+                      TargetSpace::derivativeOfDistanceSquaredWRTSecondArgument(coefficients_[i], x));
+    }
 
-  void assembleEmbeddedGradient(const TargetSpace& x,
-                                typename TargetSpace::EmbeddedTangentVector& gradient) const
-  {
-    gradient = 0;
-    for (size_t i=0; i<coefficients_.size(); i++)
-      gradient.axpy(weights_[i],
-                    TargetSpace::derivativeOfDistanceSquaredWRTSecondArgument(coefficients_[i], x));
-  }
-
-  void assembleGradient(const TargetSpace& x,
-                        typename TargetSpace::TangentVector& gradient) const
-  {
-    typename TargetSpace::EmbeddedTangentVector embeddedGradient;
-    assembleEmbeddedGradient(x,embeddedGradient);
+    void assembleGradient(const TargetSpace& x,
+                          typename TargetSpace::TangentVector& gradient) const
+    {
+      typename TargetSpace::EmbeddedTangentVector embeddedGradient;
+      assembleEmbeddedGradient(x,embeddedGradient);
 
-    Dune::FieldMatrix<ctype,size,embeddedSize> orthonormalFrame = x.orthonormalFrame();
-    orthonormalFrame.mv(embeddedGradient,gradient);
-  }
+      Dune::FieldMatrix<ctype,size,embeddedSize> orthonormalFrame = x.orthonormalFrame();
+      orthonormalFrame.mv(embeddedGradient,gradient);
+    }
 
-  void assembleEmbeddedHessian(const TargetSpace& x,
-                               Dune::SymmetricMatrix<ctype,embeddedSize>& matrix) const
-  {
-    matrix = 0;
-    for (size_t i=0; i<coefficients_.size(); i++)
-      matrix.axpy(weights_[i], TargetSpace::secondDerivativeOfDistanceSquaredWRTSecondArgument(coefficients_[i], x));
-  }
-
-  // TODO Use a Symmetric matrix for the result
-  void assembleHessian(const TargetSpace& x,
-                       Dune::FieldMatrix<ctype,size,size>& matrix) const
-  {
-    Dune::SymmetricMatrix<ctype,embeddedSize> embeddedHessian;
-    assembleEmbeddedHessian(x,embeddedHessian);
+    void assembleEmbeddedHessian(const TargetSpace& x,
+                                 SymmetricMatrix<ctype,embeddedSize>& matrix) const
+    {
+      matrix = 0;
+      for (size_t i=0; i<coefficients_.size(); i++)
+        matrix.axpy(weights_[i], TargetSpace::secondDerivativeOfDistanceSquaredWRTSecondArgument(coefficients_[i], x));
+    }
 
-    Dune::FieldMatrix<ctype,size,embeddedSize> frame = x.orthonormalFrame();
+    // TODO Use a Symmetric matrix for the result
+    void assembleHessian(const TargetSpace& x,
+                         Dune::FieldMatrix<ctype,size,size>& matrix) const
+    {
+      SymmetricMatrix<ctype,embeddedSize> embeddedHessian;
+      assembleEmbeddedHessian(x,embeddedHessian);
+
+      Dune::FieldMatrix<ctype,size,embeddedSize> frame = x.orthonormalFrame();
+
+      // this is frame * embeddedHessian * frame^T
+      for (int i=0; i<size; i++)
+        for (int j=0; j<size; j++)
+          matrix[i][j] = embeddedHessian.energyScalarProduct(frame[i], frame[j]);
+    }
 
-    // this is frame * embeddedHessian * frame^T
-    for (int i=0; i<size; i++)
-      for (int j=0; j<size; j++)
-        matrix[i][j] = embeddedHessian.energyScalarProduct(frame[i], frame[j]);
-  }
+    const std::vector<TargetSpace> coefficients_;
 
-  const std::vector<TargetSpace> coefficients_;
+    std::vector<WeightType> weights_;
 
-  std::vector<WeightType> weights_;
+  };
 
-};
+}  // namespace Dune::GFE
 
 #endif

dune/gfe/bendingisometryhelper.hh

+#ifndef DUNE_GFE_BENDINGISOMETRYHELPER_HH
+#define DUNE_GFE_BENDINGISOMETRYHELPER_HH
+
+#include <dune/common/parametertree.hh>
+#include <dune/common/parametertreeparser.hh>
+
+#include <dune/geometry/quadraturerules.hh>
+
+#include <dune/functions/functionspacebases/cubichermitebasis.hh>
+#include <dune/functions/common/differentiablefunction.hh>
+#include <dune/functions/common/differentiablefunctionfromcallables.hh>
+#include <dune/functions/common/functionconcepts.hh>
+
+
+#include <dune/gfe/spaces/productmanifold.hh>
+#include <dune/gfe/spaces/realtuple.hh>
+#include <dune/gfe/spaces/rotation.hh>
+#include <dune/gfe/functions/embeddedglobalgfefunction.hh>
+#include <dune/gfe/functions/localprojectedfefunction.hh>
+
+#include <dune/matrix-vector/crossproduct.hh>
+
+/** \file
+ * \brief Various helper functions related to bending isometries.
+ */
+
+
+namespace Dune::GFE::Impl
+{
+
+
+  /** \brief The identity function in R^d, and its derivative
+   *
+   * This is needed to compute the displacement from the deformation (for visualization).
+   * Since we want to represent it in a Hermite finite element space, we need
+   * the derivative as well.
+   *
+   * Is there no shorter way to implement this?
+   */
+  template<class K>
+  class IdentityGridEmbedding
+  {
+  public:
+    //! Evaluate identity
+    Dune::FieldVector<K,3> operator() (const Dune::FieldVector<K,2>& x) const
+    {
+      return {x[0], x[1], 0.0};
+    }
+
+    /**
+     * \brief Obtain derivative of identity function
+     */
+    friend auto derivative(const IdentityGridEmbedding& p)
+    {
+      return [](const Dune::FieldVector<K,2>& x) { return Dune::FieldMatrix<K,3,2>({{1,0}, {0,1}, {0,0}}); };
+    }
+  };
+
+
+  /**
+      Define some Generic lambda expressions (C++ 20)
+   */
+  auto cancelLastMatrixColumn = []<typename RT>(Dune::FieldMatrix<RT,3,3> matrix) -> Dune::FieldMatrix<RT,3,2>
+  {
+    return Dune::FieldMatrix<RT,3,2>{{matrix[0][0],matrix[0][1]},
+      {matrix[1][0],matrix[1][1]},
+      {matrix[2][0],matrix[2][1]}
+    };
+  };
+
+  auto rotationExtraction = []<typename RT>(Dune::FieldVector<RT,7> vec) -> Rotation<RT,3>
+  {
+    std::array<RT,4> quaternion {vec[3],vec[4],vec[5],vec[6]};
+    return Rotation<RT,3>(quaternion);
+  };
+
+  auto deformationExtraction = []<typename RT>(Dune::FieldVector<RT,7> vec) -> Dune::FieldVector<RT,3>
+  {
+    return Dune::FieldVector<RT,3>{vec[0],vec[1],vec[2]};
+  };
+
+
+
+
+
+
+
+  /**
+   * @brief Data structure conversion from 'VectorSpaceCoefficients' to 'IsometryCoefficients'.
+   */
+  template<class DiscreteKirchhoffBasis, class CoefficientBasis, class Coefficients, class IsometryCoefficients>
+  void vectorToIsometryCoefficientMap(const DiscreteKirchhoffBasis& discreteKirchhoffBasis,
+                                      const CoefficientBasis& coefficientBasis,
+                                      const Coefficients& vectorCoefficients,
+                                      IsometryCoefficients& isometryCoefficients)
+  {
+    isometryCoefficients.resize(coefficientBasis.size());
+    auto localView = discreteKirchhoffBasis.localView();
+    auto localViewCoefficients = coefficientBasis.localView();
+
+    using RT = typename IsometryCoefficients::value_type::ctype;
+
+    for (const auto &element :  elements(coefficientBasis.gridView()))
+    {
+      localView.bind(element);
+      localViewCoefficients.bind(element);
+
+      /**
+         Create data structures to store deformation values and the deformation gradient (Dofs)
+         on the current element. We store the values in an std::vector where each entry
+         contains the deformation vector (resp. x,y derivative vectors) at one node of the element.
+       */
+      std::vector<Dune::FieldVector<RT,3> > deformationDofs(3);
+      std::vector<Dune::FieldVector<RT,3> > xDerivativeDofs(3);
+      std::vector<Dune::FieldVector<RT,3> > yDerivativeDofs(3);
+
+      // Loop over components of power basis
+      for(std::size_t k=0; k<3 ; k++)
+      {
+        // This could be placed out of the loop if we assume a power basis.
+        const auto& hermiteLFE = localView.tree().child(k).finiteElement();
+
+        for(std::size_t i=0; i<hermiteLFE.size(); i++)
+        {
+          auto localIdx = localView.tree().child(k).localIndex(i);
+          auto globalIdx = localView.index(localIdx);
+          // Get current node index (Caution: This assumes that the hermite LocalFiniteElement only contains vertex dofs.)
+          auto nodeIdx = hermiteLFE.localCoefficients().localKey(i).subEntity();
+
+          /**
+             We use the functionalDescriptor of the (discreteKirchhoffBasis) hermite LocalFiniteELement
+             to determine whether a current DOF corresponds to a function evaluation or partial derivative w.r.t
+             the first or second variable.
+             (Alternatively to the - experimental - FunctionalDescriptor: use localKey(i).index(): 0~value, 1~xDerivative, 2~yDerivative.)
+           */
+          if(std::array<unsigned int,2> val {0,0} ; hermiteLFE.localInterpolation().functionalDescriptor(i).partialDerivativeOrder() == val  )
+          {
+            deformationDofs[nodeIdx][k] = vectorCoefficients[globalIdx[0]][globalIdx[1]];
+          }
+          if(std::array<unsigned int,2> val {1,0} ; hermiteLFE.localInterpolation().functionalDescriptor(i).partialDerivativeOrder() == val  )
+          {
+            xDerivativeDofs[nodeIdx][k] = vectorCoefficients[globalIdx[0]][globalIdx[1]];
+          }
+          if(std::array<unsigned int,2> val {0,1} ; hermiteLFE.localInterpolation().functionalDescriptor(i).partialDerivativeOrder() == val  )
+          {
+            yDerivativeDofs[nodeIdx][k] = vectorCoefficients[globalIdx[0]][globalIdx[1]];
+          }
+        }
+      }
+
+      // Loop over vertices of the current element
+      for(std::size_t c=0; c<3 ; c++)
+      {
+        /** Normalize derivative vectors (this is necessary due to the h-dependent scaling used in the hermite-Basis). */
+        xDerivativeDofs[c] /= xDerivativeDofs[c].two_norm();
+        yDerivativeDofs[c] /= yDerivativeDofs[c].two_norm();
+
+        /**
+            cross product to get last column of rotation matrix.
+         */
+        Dune::FieldVector<RT,3> cross = Dune::MatrixVector::crossProduct(xDerivativeDofs[c],yDerivativeDofs[c]);
+
+        Dune::FieldMatrix<RT,3,3> rotMatrix(0);
+
+        for(std::size_t k=0; k<3; k++)
+        {
+          rotMatrix[k][0] = xDerivativeDofs[c][k];
+          rotMatrix[k][1] = yDerivativeDofs[c][k];
+          rotMatrix[k][2] = cross[k];
+        }
+
+        /**
+           Check if derivative vectors are orthonormal.
+         */
+        Dune::FieldMatrix<double,3,3> I = Dune::ScaledIdentityMatrix<double,3>(1);
+
+        if(((rotMatrix.transposed()*rotMatrix) - I).frobenius_norm()>1e-8)
+          DUNE_THROW(Dune::Exception, "Error: input rotation is not orthonormal!"<< ((rotMatrix.transposed()*rotMatrix) - I).frobenius_norm());
+
+        size_t localIdxOut = localViewCoefficients.tree().localIndex(c);
+        size_t globalIdxOut = localViewCoefficients.index(localIdxOut);
+
+        using namespace Dune::Indices;
+        isometryCoefficients[globalIdxOut][_1].set(rotMatrix);
+        isometryCoefficients[globalIdxOut][_0] = (RealTuple<RT, 3>)deformationDofs[c];
+      }
+    }
+  }
+
+
+  /**
+   * @brief Data structure conversion from 'IsometryCoefficients' to 'VectorSpaceCoefficients'.
+   */
+  template<class DiscreteKirchhoffBasis, class CoefficientBasis, class Coefficients, class IsometryCoefficients>
+  void isometryToVectorCoefficientMap(const DiscreteKirchhoffBasis& discreteKirchhoffBasis,
+                                      const CoefficientBasis& coefficientBasis,
+                                      Coefficients& vectorCoefficients,
+                                      const IsometryCoefficients& isometryCoefficients)
+  {
+    auto localView = discreteKirchhoffBasis.localView();
+    auto localViewCoefficients = coefficientBasis.localView();
+
+    auto gridView = discreteKirchhoffBasis.gridView();
+
+    /** Create an EmbeddedGlobalGFEFunction that serves as a GridViewFunction later on.
+        We need a function to be passed into the interpolation method of the hermite basis
+        (although the values in between grid nodes do not matter).
+
+        The range type of EmbeddedGlobalFEFunction is a FieldVector<double,7>
+        where the first three entries correspond to the deformation and the
+        last four entries correspond to a (quaternion) rotation.
+     */
+    typedef GFE::LocalProjectedFEFunction<CoefficientBasis, GFE::ProductManifold<RealTuple<double,3>, Rotation<double,3> > > LocalInterpolationRule;
+    Dune::GFE::EmbeddedGlobalGFEFunction<CoefficientBasis, LocalInterpolationRule, Dune::GFE::ProductManifold<RealTuple<double,3>, Rotation<double,3> > > embeddedGlobalFunction(coefficientBasis, isometryCoefficients);
+
+    auto productSpaceGridViewFunction = Dune::Functions::makeAnalyticGridViewFunction(embeddedGlobalFunction, gridView);
+
+    /** Extract the deformation and rotation into separate functions. */
+    auto deformationGlobalFunction = Dune::Functions::makeComposedGridFunction(deformationExtraction,productSpaceGridViewFunction);
+
+
+
+    auto rotationGlobalFunction = Dune::Functions::makeComposedGridFunction(cancelLastMatrixColumn,
+                                                                            Dune::Functions::makeComposedGridFunction(Rotation<double,3>::quaternionToMatrix,
+                                                                                                                      Dune::Functions::makeComposedGridFunction(rotationExtraction,productSpaceGridViewFunction)));
+
+    using Domain = const Dune::FieldVector<double,2>&;
+    using Range = Dune::FieldVector<double,3>;
+    auto globalIsometryFunction = makeDifferentiableFunctionFromCallables(Dune::Functions::SignatureTag<Range(Domain)>(),deformationGlobalFunction,rotationGlobalFunction);
+
+    /** Interpolate into the hermite basis to get the coefficient vector. */
+    interpolate(discreteKirchhoffBasis, vectorCoefficients, globalIsometryFunction);
+  }
+
+
+  template<class LocalDiscreteKirchhoffFunction, class NormalBasis>
+  auto computeDiscreteSurfaceNormal(LocalDiscreteKirchhoffFunction& deformationFunction,
+                                    const NormalBasis& normalBasis)
+  {
+    std::vector<Dune::FieldVector<double,3> > discreteNormalVectorCoefficients(normalBasis.size());
+
+    auto localView = normalBasis.localView();
+    for (auto&& element : elements(normalBasis.gridView()))
+    {
+      auto geometry = element.geometry();
+
+      localView.bind(element);
+      const auto nSf = localView.tree().child(0).finiteElement().localBasis().size();
+
+      deformationFunction.bind(element);
+
+      for (int i=0; i<geometry.corners(); i++)
+      {
+        auto numDir = deformationFunction.evaluateDerivative(geometry.local(geometry.corner(i)));
+
+        /**
+            Get the normal vector via cross-product.
+         */
+        Dune::FieldVector<double,3> normal = {numDir[1][0]*numDir[2][1] - numDir[1][1]*numDir[2][0],
+                                              numDir[2][0]*numDir[0][1] - numDir[0][0]*numDir[2][1],
+                                              numDir[0][0]*numDir[1][1] - numDir[1][0]*numDir[0][1]};
+        // printvector(std::cout, normal, "discrete normal:" , "--");
+        for (size_t k=0; k < 3; k++)
+          for (size_t i=0; i < nSf; i++)
+          {
+            size_t localIdx = localView.tree().child(k).localIndex(i); // hier i:leafIdx
+            size_t globalIdx = localView.index(localIdx);
+            discreteNormalVectorCoefficients[globalIdx] = normal[k];
+          }
+      }
+    }
+
+    return discreteNormalVectorCoefficients;
+  }
+
+
+  /**
+   * @brief  Compute Coefficients of the local discrete gradient operator.
+   *
+   * The discrete Gradient is a special linear combination represented in a [P2]^2  space (locally by a representation local finite element)
+   * The coefficients of this linear combination correspond to certain linear combinations of the Gradients of localfunction_ .
+   * The coefficients are stored in the form [Basisfunctions x components x gridDim]
+   * in a BlockVector<FieldMatrix<RT, 3, gridDim>> .
+   */
+  template<class Coefficient,class LocalP2Element, class DeformationFunction >
+  auto discreteGradientCoefficients(Coefficient& discreteGradientCoefficients,
+                                    LocalP2Element& representationLFE,
+                                    DeformationFunction& deformationFunction,
+                                    auto& element)
+  {
+    using RT = typename DeformationFunction::RT;
+
+    auto geometry = element.geometry();
+    constexpr static int gridDim = DeformationFunction::gridDim;
+    auto gridView = deformationFunction.gridView();
+    const auto &indexSet = gridView.indexSet();
+    discreteGradientCoefficients.resize(representationLFE.size());
+    /**
+     * @brief On the current element we need to assign the coefficients of the discrete Gradient to the right P2-Lagrange-Basisfunctions.
+     *        this is done by iterating over all P2-Lagrange-Basisfunctions and determine the corresponding coefficients.
+     *        We distinguish two situations:
+     *        - [1] If the Lagrange-basisfunction is assigned to a node: the corresponding coefficient is just the evaluation of the deformation gradient at that specific node.
+     *        - [2] If the Lagrange-basisfunction is assigned to an edge center: the corresponding coefficient is decomposed into a normal-& tangential part,
+     *              where the normal-part: is given by the affine combination of the deformation gradient values at the edge-vertices
+     *              and the tangential-part: is given by the tangential part of the deformation gradient at the edge-center.
+     *
+     *        * We do this for all (k=3) components of the deformation-function simulatneously.
+     */
+    for(std::size_t i=0; i<representationLFE.size(); i++)
+    {
+      // [1] Determine coefficients for basis functions assigned to a vertex
+      if (representationLFE.localCoefficients().localKey(i).codim() == 2)
+      {
+        size_t nodeIndex = representationLFE.localCoefficients().localKey(i).subEntity();     // This gives index on reference Element ?
+        // Alternative: use subEntity to get (local) node position
+
+        typename DeformationFunction::DerivativeType whJacobianValue = deformationFunction.evaluateDerivative(geometry.local(geometry.corner(nodeIndex)));
+        discreteGradientCoefficients[i] = whJacobianValue;
+      }
+
+      // [2] Determine coefficients for basis functions assigned to an edge
+      if (representationLFE.localCoefficients().localKey(i).codim() == 1)
+      {
+        /**
+         * @brief On each edge the (componentwise) value of the discete jacobian on the edge midpoint
+         * is decomposed into a normal and tangent component.
+         * The normal component at the edge center is given as the affine combination of the
+         * normal component of the jacobian of localfunctions_ at the two vertices on each edge.
+         * The tangential component is the tangent component of the jacobian of localfunctions_ on
+         * the edge center.
+         */
+        size_t edgeIndex =  representationLFE.localCoefficients().localKey(i).subEntity();
+        auto edgeGeometry = element.template subEntity<1>(edgeIndex).geometry();
+
+        /**
+         * @brief Get the right orientation for tangent and normal vectors on current edge.
+         */
+        const auto &refElement = referenceElement(element);
+        // Local vertex indices within the element
+        auto localV0 = refElement.subEntity(edgeIndex, gridDim-1, 0, gridDim);
+        auto localV1 = refElement.subEntity(edgeIndex, gridDim-1, 1, gridDim);
+
+        // Global vertex indices within the grid
+        auto globalV0 = indexSet.subIndex(element, localV0, gridDim);
+        auto globalV1 = indexSet.subIndex(element, localV1, gridDim);
+
+        double edgeOrientation = 1.0;
+
+        if ((localV0<localV1 && globalV0>globalV1) || (localV0>localV1 && globalV0<globalV1))
+          edgeOrientation = -1.0;
+
+        // Get tangent vector on current edge
+        Dune::FieldVector<RT,2> tmp = (geometry.corner(localV1) - geometry.corner(localV0)) / abs(edgeGeometry.volume());
+
+        //convert to FieldMatrix
+        Dune::FieldMatrix<RT,2,1> edgeTangent = {tmp[0], tmp[1]};
+        edgeTangent *= edgeOrientation;
+
+        /**
+         * @brief Get edge normal by rotating the edge-tangent by by pi/2 clockwise
+         */
+        Dune::FieldMatrix<RT,2,1> edgeNormal = {edgeTangent[1], -1.0*edgeTangent[0]};
+
+        /**
+         * @brief Evaluate Jacobian of localfunctions_ at edge-vertices and edge-midpoints.
+         */
+        typename DeformationFunction::DerivativeType whJacobianValueCenter = deformationFunction.evaluateDerivative(geometry.local(edgeGeometry.center()));  // Q: Are these transformations correct?
+        typename DeformationFunction::DerivativeType whJacobianValueFirstVertex = deformationFunction.evaluateDerivative(geometry.local(geometry.corner(localV0)));
+        typename DeformationFunction::DerivativeType whJacobianValueSecondVertex = deformationFunction.evaluateDerivative(geometry.local(geometry.corner(localV1)));
+
+        // @OS: Is there a elementwise/hadamard matrix multiplication in Dune?
+        Dune::FieldMatrix<RT,3,gridDim> normalComponent(0);
+        Dune::FieldMatrix<RT,3,gridDim> tangentialComponent(0);
+
+        auto normalComponentFactor = 0.5 * (whJacobianValueFirstVertex + whJacobianValueSecondVertex) * edgeNormal;
+        auto tangentialComponentFactor = whJacobianValueCenter * edgeTangent;
+
+        /**
+         * @brief Compute normal and tangential component on current edge.
+         */
+        for (int k=0; k<3; k++)
+          for (int l=0; l<gridDim; l++)
+          {
+            normalComponent[k][l]   = normalComponentFactor[k][0]*edgeNormal[l];
+            tangentialComponent[k][l] = tangentialComponentFactor[k][0]*edgeTangent[l];
+          }
+        discreteGradientCoefficients[i] = normalComponent + tangentialComponent;
+
+      }
+    } //end of P2-indices
+  }
+
+
+
+} // end namespace Dune::GFE::Impl
+
+#endif

dune/gfe/cosseratstrain.hh

@@ -3,106 +3,109 @@
 
 #include <dune/common/fmatrix.hh>
 
-namespace Dune {
-
-  namespace GFE {
-
-    /** \brief Strain tensor of a Cosserat material
-     *
-     * Mathematically, this object is a dimworld x dimworld matrix.
-     * However, the last dimworld-dim columns consists of the corresponding
-     * columns of the identity matrix.  Knowing this allows more efficient
-     * implementations of several frequently occurring operations.
+namespace Dune::GFE
+{
+
+  /** \brief Strain tensor of a Cosserat material
+   *
+   * Mathematically, this object is a dimworld x dimworld matrix.
+   * However, the last dimworld-dim columns consists of the corresponding
+   * columns of the identity matrix.  Knowing this allows more efficient
+   * implementations of several frequently occurring operations.
+   *
+   * \tparam T Type used for real numbers
+   * \tparam dimworld Dimension of the world space
+   * \tparam dim Dimension of the domain
+   */
+  template <class T, int dimworld, int dim>
+  class CosseratStrain
+  {
+  public:
+    /** \brief Compute strain from the deformation gradient and the microrotation
      *
-     * \tparam T Type used for real numbers
-     * \tparam dimworld Dimension of the world space
-     * \tparam dim Dimension of the domain
+     * \param R The microrotation
      */
-    template <class T, int dimworld, int dim>
-    class CosseratStrain
+    CosseratStrain(const FieldMatrix<T,dimworld+4,dim>& deformationGradient,
+                   const FieldMatrix<T,dimworld,dimworld>& R)
     {
-    public:
-      /** \brief Compute strain from the deformation gradient and the microrotation
-       *
-       * \param R The microrotation
+      /* Compute F, the deformation gradient.
+         In the continuum case (domain dimension == world dimension) this is
+         \f$ \hat{F} = \nabla m  \f$
+         In the case of a shell it is
+         \f$ \hat{F} = (\nabla m | \overline{R}_3) \f$
        */
-      CosseratStrain(const FieldMatrix<T,dimworld+4,dim>& deformationGradient,
-                     const FieldMatrix<T,dimworld,dimworld>& R)
-      {
-        /* Compute F, the deformation gradient.
-           In the continuum case (domain dimension == world dimension) this is
-           \f$ \hat{F} = \nabla m  \f$
-           In the case of a shell it is
-           \f$ \hat{F} = (\nabla m | \overline{R}_3) \f$
-         */
-        FieldMatrix<T,dimworld,dimworld> F;
-        for (int i=0; i<dimworld; i++)
-          for (int j=0; j<dim; j++)
-            F[i][j] = deformationGradient[i][j];
-
-        for (int i=0; i<dimworld; i++)
-          for (int j=dim; j<dimworld; j++)
-            F[i][j] = R[i][j];
-
-        // U = R^T F
-        for (int i=0; i<dimworld; i++)
-          for (int j=0; j<dimworld; j++) {
-            data_[i][j] = 0;
-            for (int k=0; k<dimworld; k++)
-              data_[i][j] += R[k][i] * F[k][j];
-          }
-
-      }
-
-      /** \brief Compute strain from the deformation gradient and the microrotation
-       *
-       * \param R The microrotation
+      FieldMatrix<T,dimworld,dimworld> F;
+      for (int i=0; i<dimworld; i++)
+        for (int j=0; j<dim; j++)
+          F[i][j] = deformationGradient[i][j];
+
+      for (int i=0; i<dimworld; i++)
+        for (int j=dim; j<dimworld; j++)
+          F[i][j] = R[i][j];
+
+      // U = R^T F
+      for (int i=0; i<dimworld; i++)
+        for (int j=0; j<dimworld; j++) {
+          data_[i][j] = 0;
+          for (int k=0; k<dimworld; k++)
+            data_[i][j] += R[k][i] * F[k][j];
+        }
+
+    }
+
+    /** \brief Compute strain from the deformation gradient and the microrotation
+     *
+     * \param R The microrotation
+     */
+    CosseratStrain(const FieldMatrix<T,dimworld,dim>& deformationGradient,
+                   const FieldMatrix<T,dimworld,dimworld>& R)
+    {
+      /* Compute F, the deformation gradient.
+         In the continuum case (domain dimension == world dimension) this is
+         \f$ \hat{F} = \nabla m  \f$
+         In the case of a shell it is
+         \f$ \hat{F} = (\nabla m | \overline{R}_3) \f$
        */
-      CosseratStrain(const FieldMatrix<T,dimworld,dim>& deformationGradient,
-                     const FieldMatrix<T,dimworld,dimworld>& R)
-      {
-        /* Compute F, the deformation gradient.
-           In the continuum case (domain dimension == world dimension) this is
-           \f$ \hat{F} = \nabla m  \f$
-           In the case of a shell it is
-           \f$ \hat{F} = (\nabla m | \overline{R}_3) \f$
-         */
-        FieldMatrix<T,dimworld,dimworld> F;
-        for (int i=0; i<dimworld; i++)
-          for (int j=0; j<dim; j++)
-            F[i][j] = deformationGradient[i][j];
-
-        for (int i=0; i<dimworld; i++)
-          for (int j=dim; j<dimworld; j++)
-            F[i][j] = R[i][j];
-
-        // U = R^T F
-        for (int i=0; i<dimworld; i++)
-          for (int j=0; j<dimworld; j++) {
-            data_[i][j] = 0;
-            for (int k=0; k<dimworld; k++)
-              data_[i][j] += R[k][i] * F[k][j];
-          }
-
-      }
-
-      T determinant() const
-      {
-        return data_.determinant();
-      }
-
-      /** \brief Temporary: get the raw matrix data */
-      const FieldMatrix<T,dimworld,dimworld>& matrix() const
-      {
-        return data_;
-      }
-
-    private:
-      FieldMatrix<T,dimworld,dimworld> data_;
-    };
-
-  }   // namespace GFE
-
-}  // namespace Dune
+      FieldMatrix<T,dimworld,dimworld> F;
+      for (int i=0; i<dimworld; i++)
+        for (int j=0; j<dim; j++)
+          F[i][j] = deformationGradient[i][j];
+
+      for (int i=0; i<dimworld; i++)
+        for (int j=dim; j<dimworld; j++)
+          F[i][j] = R[i][j];
+
+      // U = R^T F
+      for (int i=0; i<dimworld; i++)
+        for (int j=0; j<dimworld; j++) {
+          data_[i][j] = 0;
+          for (int k=0; k<dimworld; k++)
+            data_[i][j] += R[k][i] * F[k][j];
+        }
+
+    }
+
+    T determinant() const
+    {
+      // If dim==2 we know that the last column is the identity matrix.
+      // Hence the determinant is equal to the determinant of the
+      // upper left 2x2 sub-block.
+      if (dim==2)
+        return data_[0][0]*data_[1][1] - data_[1][0]*data_[0][1];
+
+      return data_.determinant();
+    }
+
+    /** \brief Temporary: get the raw matrix data */
+    const FieldMatrix<T,dimworld,dimworld>& matrix() const
+    {
+      return data_;
+    }
+
+  private:
+    FieldMatrix<T,dimworld,dimworld> data_;
+  };
+
+}  // namespace Dune::GFE
 
 #endif   // DUNE_GFE_COSSERATSTRAIN_HH

dune/gfe/cosseratvtkreader.hh

@@ -5,42 +5,38 @@
 #include <dune/gfe/spaces/realtuple.hh>
 #include <dune/gfe/spaces/rotation.hh>
 
-namespace Dune
+namespace Dune::GFE
 {
-  namespace GFE
+
+  /** \brief Read configurations of Cosserat models from VTK files into memory */
+  class CosseratVTKReader
   {
+  public:
 
-    /** \brief Read configurations of Cosserat models from VTK files into memory */
-    class CosseratVTKReader
+    static void read(std::vector<ProductManifold<RealTuple<double,3>,Rotation<double,3> > >& configuration,
+                     const std::string& filename)
     {
-    public:
-
-      static void read(std::vector<ProductManifold<RealTuple<double,3>,Rotation<double,3> > >& configuration,
-                       const std::string& filename)
-      {
-        VTKFile vtkFile;
-        vtkFile.read(filename);
+      VTKFile vtkFile;
+      vtkFile.read(filename);
 
-        configuration.resize(vtkFile.points_.size());
+      configuration.resize(vtkFile.points_.size());
 
-        for (size_t i=0; i<configuration.size(); i++)
-        {
-          configuration[i][Indices::_0].globalCoordinates() = vtkFile.points_[i];
-
-          FieldMatrix<double,3,3> R;
-          for (int j=0; j<3; j++)
-            for (int k=0; k<3; k++)
-              R[j][k] = vtkFile.directors_[k][i][j];
+      for (size_t i=0; i<configuration.size(); i++)
+      {
+        configuration[i][Indices::_0].globalCoordinates() = vtkFile.points_[i];
 
-          configuration[i][Indices::_1].set(R);
-        }
+        FieldMatrix<double,3,3> R;
+        for (int j=0; j<3; j++)
+          for (int k=0; k<3; k++)
+            R[j][k] = vtkFile.directors_[k][i][j];
 
+        configuration[i][Indices::_1].set(R);
       }
 
-    };
+    }
 
-  }
+  };
 
-}
+}  // namespace Dune::GFE
 
 #endif

dune/gfe/cosseratvtkwriter.hh

@@ -6,9 +6,15 @@
 #include <dune/grid/io/file/vtk/vtkwriter.hh>
 #include <dune/grid/io/file/vtk/pvtuwriter.hh>
 
+#include <dune/istl/bvector.hh>
+
 #include <dune/functions/functionspacebases/lagrangebasis.hh>
 #include <dune/functions/functionspacebases/interpolate.hh>
 #include <dune/functions/gridfunctions/discreteglobalbasisfunction.hh>
+#include <dune/functions/gridfunctions/composedgridfunction.hh>
+
+#include <dune/vtk/vtkwriter.hh>
+#include <dune/vtk/datacollectors/lagrangedatacollector.hh>
 
 #include <dune/gfe/vtkfile.hh>
 #include <dune/gfe/spaces/productmanifold.hh>
@@ -16,553 +22,435 @@
 #include <dune/gfe/spaces/rotation.hh>
 
 
-/** \brief Write the configuration of a Cosserat material in VTK format */
-template <class GridType>
-class CosseratVTKWriter
+namespace Dune::GFE
 {
 
-  static const int dim = GridType::dimension;
-
-  template <typename Basis1, typename Basis2>
-  static void downsample(const Basis1& basis1, const std::vector<Dune::GFE::ProductManifold<RealTuple<double,3>,Rotation<double,3> > >& v1,
-                         const Basis2& basis2,       std::vector<Dune::GFE::ProductManifold<RealTuple<double,3>,Rotation<double,3> > >& v2)
-  {
-    // Embed v1 into R^7
-    std::vector<Dune::FieldVector<double,7> > v1Embedded(v1.size());
-    for (size_t i=0; i<v1.size(); i++)
-      v1Embedded[i] = v1[i].globalCoordinates();
-
-    // Interpolate
-    auto function = Dune::Functions::makeDiscreteGlobalBasisFunction<Dune::FieldVector<double,7> >(basis1, v1Embedded);
-    std::vector<Dune::FieldVector<double,7> > v2Embedded;
-    Dune::Functions::interpolate(basis2, v2Embedded, function);
-
-    // Copy back from R^7 into ProductManifold
-    v2.resize(v2Embedded.size());
-    for (size_t i=0; i<v2.size(); i++)
-      v2[i] = Dune::GFE::ProductManifold<RealTuple<double,3>,Rotation<double,3> >(v2Embedded[i]);
-  }
-
-  template <typename Basis1, typename Basis2>
-  static void downsample(const Basis1& basis1, const std::vector<RealTuple<double,3> >& v1,
-                         const Basis2& basis2,       std::vector<RealTuple<double,3> >& v2)
-  {
-    // Copy from RealTuple to FieldVector
-    std::vector<Dune::FieldVector<double,3> > v1Embedded(v1.size());
-    for (size_t i=0; i<v1.size(); i++)
-      v1Embedded[i] = v1[i].globalCoordinates();
-
-    // Interpolate
-    auto function = Dune::Functions::makeDiscreteGlobalBasisFunction<Dune::FieldVector<double,3> >(basis1, v1Embedded);
-    std::vector<Dune::FieldVector<double,3> > v2Embedded;
-    Dune::Functions::interpolate(basis2, v2Embedded, function);
-
-    // Copy back from FieldVector to RealTuple
-    v2.resize(v2Embedded.size());
-    for (size_t i=0; i<v2.size(); i++)
-      v2[i] = RealTuple<double,3>(v2Embedded[i]);
-  }
-
-  /** \brief Extend filename to contain communicator rank and size
-   *
-   * Copied from dune-grid vtkwriter.hh
-   */
-  static std::string getParallelPieceName(const std::string& name,
-                                          const std::string& path,
-                                          int commRank, int commSize)
-  {
-    std::ostringstream s;
-    if(path.size() > 0) {
-      s << path;
-      if(path[path.size()-1] != '/')
-        s << '/';
-    }
-    s << 's' << std::setw(4) << std::setfill('0') << commSize << '-';
-    s << 'p' << std::setw(4) << std::setfill('0') << commRank << '-';
-    s << name;
-    if(GridType::dimension > 1)
-      s << ".vtu";
-    else
-      s << ".vtp";
-    return s.str();
-  }
-
-  /** \brief Extend filename to contain communicator rank and size
-   *
-   * Copied from dune-grid vtkwriter.hh
-   */
-  static std::string getParallelName(const std::string& name,
-                                     const std::string& path,
-                                     int commSize)
-  {
-    std::ostringstream s;
-    if(path.size() > 0) {
-      s << path;
-      if(path[path.size()-1] != '/')
-        s << '/';
-    }
-    s << 's' << std::setw(4) << std::setfill('0') << commSize << '-';
-    s << name;
-    if(GridType::dimension > 1)
-      s << ".pvtu";
-    else
-      s << ".pvtp";
-    return s.str();
-  }
-
-public:
-  /** \brief Write a configuration given with respect to a scalar function space basis
-   */
-  template <typename Basis>
-  static void write(const Basis& basis,
-                    const Dune::TupleVector<std::vector<RealTuple<double,3> >,
-                        std::vector<Rotation<double,3> > >& configuration,
-                    const std::string& filename)
-  {
-    using namespace Dune::TypeTree::Indices;
-    std::vector<Dune::GFE::ProductManifold<RealTuple<double,3>,Rotation<double,3> > > xRBM(basis.size());
-    for (std::size_t i = 0; i < basis.size(); i++) {
-      for (int j = 0; j < 3; j ++)   // Displacement part
-        xRBM[i][_0].globalCoordinates()[j] = configuration[_0][i][j];
-      xRBM[i][_1] = configuration[_1][i];      // Rotation part
-    }
-    write(basis,xRBM,filename);
-  }
-  /** \brief Write a configuration given with respect to a scalar function space basis
-   */
-  template <typename Basis, typename VectorType>
-  static void write(const Basis& basis,
-                    const VectorType& configuration,
-                    const std::string& filename)
-  {
-    using namespace Dune::TypeTree::Indices;
-    std::vector<Dune::GFE::ProductManifold<RealTuple<double,3>,Rotation<double,3> > > xRBM(basis.size());
-    for (std::size_t i = 0; i < basis.size(); i++) {
-      for (int j = 0; j < 3; j ++)   // Displacement part
-        xRBM[i][_0].globalCoordinates()[j] = configuration[i][j];
-    }
-    write(basis,xRBM,filename);
-  }
-
-  /** \brief Write a configuration given with respect to a scalar function space basis
-   */
-  template <typename Basis>
-  static void write(const Basis& basis,
-                    const std::vector<Dune::GFE::ProductManifold<RealTuple<double,3>,Rotation<double,3> > >& configuration,
-                    const std::string& filename)
+  /** \brief Write the configuration of a Cosserat material in VTK format */
+  template <class GridView>
+  class CosseratVTKWriter
   {
-    using namespace Dune::Indices;
 
-    assert(basis.size() == configuration.size());
-    auto gridView = basis.gridView();
+    static const int dim = GridView::dimension;
 
-    // Determine order of the basis
-    // We check for the order of the first element, and assume it is the same for all others
-    auto localView = basis.localView();
-    localView.bind(*gridView.template begin<0>());
-    const int order = localView.tree().finiteElement().localBasis().order();
-    // order of the approximation of the VTK file -- can only be two or one
-    const auto vtkOrder = std::min(2,order);
-
-    //  Downsample 3rd-order functions onto a P2-space.  That's all VTK can visualize today.
-    if (order>=3)
+    template <typename Basis1, typename Basis2>
+    static void downsample(const Basis1& basis1, const std::vector<Dune::GFE::ProductManifold<RealTuple<double,3>,Rotation<double,3> > >& v1,
+                           const Basis2& basis2,       std::vector<Dune::GFE::ProductManifold<RealTuple<double,3>,Rotation<double,3> > >& v2)
     {
-      using namespace Dune::Functions::BasisFactory;
-      auto p2Basis = makeBasis(gridView, lagrange<2>());
-
-      auto blockedP2Basis = makeBasis(
-        gridView,
-        power<3>(
-          lagrange<2>(),
-          blockedInterleaved()
-          ));
-
-      std::vector<Dune::GFE::ProductManifold<RealTuple<double,3>,Rotation<double,3> > > downsampledConfig;
-
-      downsample(basis, configuration, blockedP2Basis, downsampledConfig);
-
-      write(p2Basis, downsampledConfig, filename);
-      return;
+      // Embed v1 into R^7
+      std::vector<Dune::FieldVector<double,7> > v1Embedded(v1.size());
+      for (size_t i=0; i<v1.size(); i++)
+        v1Embedded[i] = v1[i].globalCoordinates();
+
+      // Interpolate
+      auto function = Dune::Functions::makeDiscreteGlobalBasisFunction<Dune::FieldVector<double,7> >(basis1, v1Embedded);
+      std::vector<Dune::FieldVector<double,7> > v2Embedded;
+      Dune::Functions::interpolate(basis2, v2Embedded, function);
+
+      // Copy back from R^7 into ProductManifold
+      v2.resize(v2Embedded.size());
+      for (size_t i=0; i<v2.size(); i++)
+        v2[i] = Dune::GFE::ProductManifold<RealTuple<double,3>,Rotation<double,3> >(v2Embedded[i]);
     }
 
-    Dune::GFE::VTKFile vtkFile;
-
-    // Count the number of elements of the different types
-    std::map<Dune::GeometryType,std::size_t> numElements;
-    for (const auto t : gridView.indexSet().types(0))
-      numElements[t] = 0;
-
-    for (auto&& t : elements(gridView, Dune::Partitions::interior))
-      numElements[t.type()]++;
-
-    std::size_t totalNumElements = 0;
-    for (const auto nE : numElements)
-      totalNumElements += nE.second;
-
-    // Enter vertex coordinates
-    std::vector<Dune::FieldVector<double, 3> > points(configuration.size());
-    for (size_t i=0; i<configuration.size(); i++)
-      points[i] = configuration[i][_0].globalCoordinates();
-
-    vtkFile.points_ = points;
+    template <typename Basis1, typename Basis2>
+    static void downsample(const Basis1& basis1, const std::vector<RealTuple<double,3> >& v1,
+                           const Basis2& basis2,       std::vector<RealTuple<double,3> >& v2)
+    {
+      // Copy from RealTuple to FieldVector
+      std::vector<Dune::FieldVector<double,3> > v1Embedded(v1.size());
+      for (size_t i=0; i<v1.size(); i++)
+        v1Embedded[i] = v1[i].globalCoordinates();
+
+      // Interpolate
+      auto function = Dune::Functions::makeDiscreteGlobalBasisFunction<Dune::FieldVector<double,3> >(basis1, v1Embedded);
+      std::vector<Dune::FieldVector<double,3> > v2Embedded;
+      Dune::Functions::interpolate(basis2, v2Embedded, function);
+
+      // Copy back from FieldVector to RealTuple
+      v2.resize(v2Embedded.size());
+      for (size_t i=0; i<v2.size(); i++)
+        v2[i] = RealTuple<double,3>(v2Embedded[i]);
+    }
 
-    // Enter elements
-    std::size_t connectivitySize = 0;
-    for (const auto nE : numElements)
+    /** \brief Extend filename to contain communicator rank and size
+     *
+     * Copied from dune-grid vtkwriter.hh
+     */
+    static std::string getParallelPieceName(const std::string& name,
+                                            const std::string& path,
+                                            int commRank, int commSize)
     {
-      if (nE.first.isQuadrilateral())
-        connectivitySize += ((vtkOrder==2) ? 8 : 4) * nE.second;
-      else if (nE.first.isTriangle())
-        connectivitySize += ((vtkOrder==2) ? 6 : 3) * nE.second;
-      else if (nE.first.isHexahedron())
-        connectivitySize += ((vtkOrder==2) ? 20 : 8) * nE.second;
-      else if (nE.first.isLine())
-        connectivitySize += ((vtkOrder==2) ? 3 : 2) * nE.second;
+      std::ostringstream s;
+      if(path.size() > 0) {
+        s << path;
+        if(path[path.size()-1] != '/')
+          s << '/';
+      }
+      s << 's' << std::setw(4) << std::setfill('0') << commSize << '-';
+      s << 'p' << std::setw(4) << std::setfill('0') << commRank << '-';
+      s << name;
+      if(GridView::dimension > 1)
+        s << ".vtu";
       else
-        DUNE_THROW(Dune::IOError, "Unsupported element type '" << nE.first << "' found!");
+        s << ".vtp";
+      return s.str();
     }
-    std::vector<unsigned int> connectivity(connectivitySize);
 
-    size_t i=0;
-    for (const auto& element : elements(gridView, Dune::Partitions::interior))
+    /** \brief Extend filename to contain communicator rank and size
+     *
+     * Copied from dune-grid vtkwriter.hh
+     */
+    static std::string getParallelName(const std::string& name,
+                                       const std::string& path,
+                                       int commSize)
     {
-      localView.bind(element);
-
-      if (element.type().isQuadrilateral())
-      {
-        if (vtkOrder==2)
-        {
-          connectivity[i++] = localView.index(0);
-          connectivity[i++] = localView.index(2);
-          connectivity[i++] = localView.index(8);
-          connectivity[i++] = localView.index(6);
-
-          connectivity[i++] = localView.index(1);
-          connectivity[i++] = localView.index(5);
-          connectivity[i++] = localView.index(7);
-          connectivity[i++] = localView.index(3);
-        }
-        else      // first order
-        {
-          connectivity[i++] = localView.index(0);
-          connectivity[i++] = localView.index(1);
-          connectivity[i++] = localView.index(3);
-          connectivity[i++] = localView.index(2);
-        }
-      }
-      if (element.type().isTriangle())
-      {
-        if (vtkOrder==2)
-        {
-          connectivity[i++] = localView.index(0);
-          connectivity[i++] = localView.index(2);
-          connectivity[i++] = localView.index(5);
-          connectivity[i++] = localView.index(1);
-          connectivity[i++] = localView.index(4);
-          connectivity[i++] = localView.index(3);
-        }
-        else      // first order
-        {
-          connectivity[i++] = localView.index(0);
-          connectivity[i++] = localView.index(1);
-          connectivity[i++] = localView.index(2);
-        }
-      }
-      if (element.type().isHexahedron())
-      {
-        if (vtkOrder==2)
-        {
-          // Corner dofs
-          connectivity[i++] = localView.index(0);
-          connectivity[i++] = localView.index(2);
-          connectivity[i++] = localView.index(8);
-          connectivity[i++] = localView.index(6);
-
-          connectivity[i++] = localView.index(18);
-          connectivity[i++] = localView.index(20);
-          connectivity[i++] = localView.index(26);
-          connectivity[i++] = localView.index(24);
-
-          // Edge dofs
-          connectivity[i++] = localView.index(1);
-          connectivity[i++] = localView.index(5);
-          connectivity[i++] = localView.index(7);
-          connectivity[i++] = localView.index(3);
-
-          connectivity[i++] = localView.index(19);
-          connectivity[i++] = localView.index(23);
-          connectivity[i++] = localView.index(25);
-          connectivity[i++] = localView.index(21);
-
-          connectivity[i++] = localView.index(9);
-          connectivity[i++] = localView.index(11);
-          connectivity[i++] = localView.index(17);
-          connectivity[i++] = localView.index(15);
-        }
-        else      // first order
-        {
-          connectivity[i++] = localView.index(0);
-          connectivity[i++] = localView.index(1);
-          connectivity[i++] = localView.index(3);
-          connectivity[i++] = localView.index(2);
-          connectivity[i++] = localView.index(4);
-          connectivity[i++] = localView.index(5);
-          connectivity[i++] = localView.index(7);
-          connectivity[i++] = localView.index(6);
-        }
+      std::ostringstream s;
+      if(path.size() > 0) {
+        s << path;
+        if(path[path.size()-1] != '/')
+          s << '/';
       }
+      s << 's' << std::setw(4) << std::setfill('0') << commSize << '-';
+      s << name;
+      if(GridView::dimension > 1)
+        s << ".pvtu";
+      else
+        s << ".pvtp";
+      return s.str();
+    }
 
-      if (element.type().isLine())
-      {
-        if (vtkOrder==2)
-        {
-          connectivity[i++] = localView.index(0);
-          connectivity[i++] = localView.index(2);
-          connectivity[i++] = localView.index(1);
-        }
-        else      // first order
-        {
-          connectivity[i++] = localView.index(0);
-          connectivity[i++] = localView.index(1);
-        }
+  public:
+    /** \brief Write a configuration given with respect to a scalar function space basis
+     */
+    template <typename Basis>
+    static void write(const Basis& basis,
+                      const Dune::TupleVector<std::vector<RealTuple<double,3> >,
+                          std::vector<Rotation<double,3> > >& configuration,
+                      const std::string& filename)
+    {
+      using namespace Dune::Indices;
+      std::vector<Dune::GFE::ProductManifold<RealTuple<double,3>,Rotation<double,3> > > xRBM(basis.size());
+      for (std::size_t i = 0; i < basis.size(); i++) {
+        for (int j = 0; j < 3; j ++) // Displacement part
+          xRBM[i][_0].globalCoordinates()[j] = configuration[_0][i][j];
+        xRBM[i][_1] = configuration[_1][i];    // Rotation part
       }
+      write(basis,xRBM,filename);
     }
-
-    vtkFile.cellConnectivity_ = connectivity;
-
-    std::vector<int> offsets(totalNumElements);
-    i = 0;
-    int offsetCounter = 0;
-    for (const auto& element : elements(gridView, Dune::Partitions::interior))
+    /** \brief Write a configuration given with respect to a scalar function space basis
+     */
+    template <typename Basis, typename VectorType>
+    static void write(const Basis& basis,
+                      const VectorType& configuration,
+                      const std::string& filename)
     {
-      if (element.type().isQuadrilateral())
-        offsetCounter += (vtkOrder==2) ? 8 : 4;
-
-      if (element.type().isTriangle())
-        offsetCounter += (vtkOrder==2) ? 6 : 3;
-
-      if (element.type().isHexahedron())
-        offsetCounter += (vtkOrder==2) ? 20 : 8;
-
-      offsets[i++] += offsetCounter;
+      using namespace Dune::Indices;
+      std::vector<Dune::GFE::ProductManifold<RealTuple<double,3>,Rotation<double,3> > > xRBM(basis.size());
+      for (std::size_t i = 0; i < basis.size(); i++) {
+        for (int j = 0; j < 3; j ++) // Displacement part
+          xRBM[i][_0].globalCoordinates()[j] = configuration[i][j];
+      }
+      write(basis,xRBM,filename);
     }
 
-    vtkFile.cellOffsets_ = offsets;
-
-    std::vector<int> cellTypes(totalNumElements);
-    i = 0;
-    for (const auto& element : elements(gridView, Dune::Partitions::interior))
+    /** \brief Write a configuration given with respect to a scalar function space basis
+     */
+    template <typename Basis>
+    static void write(const Basis& basis,
+                      const std::vector<Dune::GFE::ProductManifold<RealTuple<double,3>,Rotation<double,3> > >& configuration,
+                      const std::string& filename)
     {
-      if (element.type().isQuadrilateral())
-        cellTypes[i++] = (vtkOrder==2) ? 23 : 9;
-
-      if (element.type().isTriangle())
-        cellTypes[i++] = (vtkOrder==2) ? 22 : 5;
-
-      if (element.type().isHexahedron())
-        cellTypes[i++] = (vtkOrder==2) ? 25 : 12;
-    }
-    vtkFile.cellTypes_ = cellTypes;
+      using namespace Dune::Indices;
 
-    // Z coordinate for better visualization of wrinkles
-    std::vector<double> zCoord(points.size());
-    for (size_t i=0; i<configuration.size(); i++)
-      zCoord[i] = configuration[i][_0].globalCoordinates()[2];
+      assert(basis.size() == configuration.size());
+      auto gridView = basis.gridView();
 
-    vtkFile.zCoord_ = zCoord;
+      // Determine order of the basis
+      // We check for the order of the first element, and assume it is the same for all others
+      auto localView = basis.localView();
+      localView.bind(*gridView.template begin<0>());
+      const int order = localView.tree().finiteElement().localBasis().order();
+      // order of the approximation of the VTK file -- can only be two or one
+      const auto vtkOrder = std::min(2,order);
 
-    // The three director fields
-    for (size_t i=0; i<3; i++)
-    {
-      vtkFile.directors_[i].resize(configuration.size());
-      for (size_t j=0; j<configuration.size(); j++)
-        vtkFile.directors_[i][j] = configuration[j][_1].director(i);
-    }
+      //  Downsample 3rd-order functions onto a P2-space.  That's all VTK can visualize today.
+      if (order>=3)
+      {
+        using namespace Dune::Functions::BasisFactory;
+        auto p2Basis = makeBasis(gridView, lagrange<2>());
 
-    // Actually write the VTK file to disk
-    vtkFile.write(filename);
-  }
-
-  /** \brief Write a configuration given with respect to a scalar function space basis
-   */
-  template <typename DisplacementBasis, typename OrientationBasis>
-  static void writeMixed(const DisplacementBasis& displacementBasis,
-                         const std::vector<RealTuple<double,3> >& deformationConfiguration,
-                         const OrientationBasis& orientationBasis,
-                         const std::vector<Rotation<double,3> >& orientationConfiguration,
-                         const std::string& filename)
-  {
-    assert(displacementBasis.size() == deformationConfiguration.size());
-    assert(orientationBasis.size()  == orientationConfiguration.size());
-    auto gridView = displacementBasis.gridView();
+        auto blockedP2Basis = makeBasis(
+          gridView,
+          power<3>(
+            lagrange<2>(),
+            blockedInterleaved()
+            ));
 
-    // Determine order of the displacement basis
-    auto localView = displacementBasis.localView();
-    localView.bind(*gridView.template begin<0>());
-    int order = localView.tree().finiteElement().localBasis().order();
+        std::vector<Dune::GFE::ProductManifold<RealTuple<double,3>,Rotation<double,3> > > downsampledConfig;
 
-    // We only do P2 spaces at the moment
-    if (order != 2 and order != 3)
-    {
-      std::cout << "Warning: CosseratVTKWriter only supports P2 spaces -- skipping" << std::endl;
-      return;
-    }
+        downsample(basis, configuration, blockedP2Basis, downsampledConfig);
 
-    std::vector<RealTuple<double,3> > displacementConfiguration = deformationConfiguration;
+        write(p2Basis, downsampledConfig, filename);
+        return;
+      }
 
-    if (order == 3)
-    {
-      using namespace Dune::Functions::BasisFactory;
+      Dune::GFE::VTKFile vtkFile;
 
-      auto p2DeformationBasis = makeBasis(
-        gridView,
-        power<3>(
-          lagrange<2>(),
-          blockedInterleaved()
-          ));
+      // Count the number of elements of the different types
+      std::map<Dune::GeometryType,std::size_t> numElements;
+      for (const auto t : gridView.indexSet().types(0))
+        numElements[t] = 0;
 
-      // resample to 2nd order -- vtk can't do anything higher
-      std::vector<RealTuple<double,3> > p2DeformationConfiguration;
+      for (auto&& t : elements(gridView, Dune::Partitions::interior))
+        numElements[t.type()]++;
 
-      downsample(displacementBasis, displacementConfiguration,
-                 p2DeformationBasis, p2DeformationConfiguration);
+      std::size_t totalNumElements = 0;
+      for (const auto nE : numElements)
+        totalNumElements += nE.second;
 
-      displacementConfiguration = p2DeformationConfiguration;
-    }
+      // Enter vertex coordinates
+      std::vector<Dune::FieldVector<double, 3> > points(configuration.size());
+      for (size_t i=0; i<configuration.size(); i++)
+        points[i] = configuration[i][_0].globalCoordinates();
 
-    std::string fullfilename = filename + ".vtu";
+      vtkFile.points_ = points;
 
-    // Prepend rank and communicator size to the filename, if there are more than one process
-    if (gridView.comm().size() > 1)
-      fullfilename = getParallelPieceName(filename, "", gridView.comm().rank(), gridView.comm().size());
+      // Enter elements
+      std::size_t connectivitySize = 0;
+      for (const auto nE : numElements)
+      {
+        if (nE.first.isQuadrilateral())
+          connectivitySize += ((vtkOrder==2) ? 8 : 4) * nE.second;
+        else if (nE.first.isTriangle())
+          connectivitySize += ((vtkOrder==2) ? 6 : 3) * nE.second;
+        else if (nE.first.isHexahedron())
+          connectivitySize += ((vtkOrder==2) ? 20 : 8) * nE.second;
+        else if (nE.first.isLine())
+          connectivitySize += ((vtkOrder==2) ? 3 : 2) * nE.second;
+        else
+          DUNE_THROW(Dune::IOError, "Unsupported element type '" << nE.first << "' found!");
+      }
+      std::vector<unsigned int> connectivity(connectivitySize);
 
-    // Write the pvtu file that ties together the different parts
-    if (gridView.comm().size() > 1 && gridView.comm().rank()==0)
-    {
-      std::ofstream pvtuOutFile(getParallelName(filename, "", gridView.comm().size()));
-      Dune::VTK::PVTUWriter writer(pvtuOutFile, Dune::VTK::unstructuredGrid);
+      size_t i=0;
+      for (const auto& element : elements(gridView, Dune::Partitions::interior))
+      {
+        localView.bind(element);
 
-      writer.beginMain();
+        if (element.type().isQuadrilateral())
+        {
+          if (vtkOrder==2)
+          {
+            connectivity[i++] = localView.index(0);
+            connectivity[i++] = localView.index(2);
+            connectivity[i++] = localView.index(8);
+            connectivity[i++] = localView.index(6);
+
+            connectivity[i++] = localView.index(1);
+            connectivity[i++] = localView.index(5);
+            connectivity[i++] = localView.index(7);
+            connectivity[i++] = localView.index(3);
+          }
+          else    // first order
+          {
+            connectivity[i++] = localView.index(0);
+            connectivity[i++] = localView.index(1);
+            connectivity[i++] = localView.index(3);
+            connectivity[i++] = localView.index(2);
+          }
+        }
+        if (element.type().isTriangle())
+        {
+          if (vtkOrder==2)
+          {
+            connectivity[i++] = localView.index(0);
+            connectivity[i++] = localView.index(2);
+            connectivity[i++] = localView.index(5);
+            connectivity[i++] = localView.index(1);
+            connectivity[i++] = localView.index(4);
+            connectivity[i++] = localView.index(3);
+          }
+          else    // first order
+          {
+            connectivity[i++] = localView.index(0);
+            connectivity[i++] = localView.index(1);
+            connectivity[i++] = localView.index(2);
+          }
+        }
+        if (element.type().isHexahedron())
+        {
+          if (vtkOrder==2)
+          {
+            // Corner dofs
+            connectivity[i++] = localView.index(0);
+            connectivity[i++] = localView.index(2);
+            connectivity[i++] = localView.index(8);
+            connectivity[i++] = localView.index(6);
+
+            connectivity[i++] = localView.index(18);
+            connectivity[i++] = localView.index(20);
+            connectivity[i++] = localView.index(26);
+            connectivity[i++] = localView.index(24);
+
+            // Edge dofs
+            connectivity[i++] = localView.index(1);
+            connectivity[i++] = localView.index(5);
+            connectivity[i++] = localView.index(7);
+            connectivity[i++] = localView.index(3);
+
+            connectivity[i++] = localView.index(19);
+            connectivity[i++] = localView.index(23);
+            connectivity[i++] = localView.index(25);
+            connectivity[i++] = localView.index(21);
+
+            connectivity[i++] = localView.index(9);
+            connectivity[i++] = localView.index(11);
+            connectivity[i++] = localView.index(17);
+            connectivity[i++] = localView.index(15);
+          }
+          else    // first order
+          {
+            connectivity[i++] = localView.index(0);
+            connectivity[i++] = localView.index(1);
+            connectivity[i++] = localView.index(3);
+            connectivity[i++] = localView.index(2);
+            connectivity[i++] = localView.index(4);
+            connectivity[i++] = localView.index(5);
+            connectivity[i++] = localView.index(7);
+            connectivity[i++] = localView.index(6);
+          }
+        }
 
-      //writer.beginPointData();
-      //writer.addArray<float>("director0", 3);
-      //writer.addArray<float>("director1", 3);
-      //writer.addArray<float>("director2", 3);
-      //writer.addArray<float>("zCoord", 1);
-      //writer.endPointData();
+        if (element.type().isLine())
+        {
+          if (vtkOrder==2)
+          {
+            connectivity[i++] = localView.index(0);
+            connectivity[i++] = localView.index(2);
+            connectivity[i++] = localView.index(1);
+          }
+          else    // first order
+          {
+            connectivity[i++] = localView.index(0);
+            connectivity[i++] = localView.index(1);
+          }
+        }
+      }
 
-      // dump point coordinates
-      writer.beginPoints();
-      writer.addArray("Coordinates", 3, Dune::VTK::Precision::float32);
-      writer.endPoints();
+      vtkFile.cellConnectivity_ = connectivity;
 
-      for (int i=0; i<gridView.comm().size(); i++)
-        writer.addPiece(getParallelPieceName(filename, "", i, gridView.comm().size()));
+      std::vector<int> offsets(totalNumElements);
+      i = 0;
+      int offsetCounter = 0;
+      for (const auto& element : elements(gridView, Dune::Partitions::interior))
+      {
+        if (element.type().isQuadrilateral())
+          offsetCounter += (vtkOrder==2) ? 8 : 4;
 
-      // finish main section
-      writer.endMain();
-    }
+        if (element.type().isTriangle())
+          offsetCounter += (vtkOrder==2) ? 6 : 3;
 
-    /////////////////////////////////////////////////////////////////////////////////
-    //  Write the actual vtu file
-    /////////////////////////////////////////////////////////////////////////////////
+        if (element.type().isHexahedron())
+          offsetCounter += (vtkOrder==2) ? 20 : 8;
 
-    // Stupid: I can't directly get the number of Interior_Partition elements
-    size_t numElements = std::distance(gridView.template begin<0, Dune::Interior_Partition>(),
-                                       gridView.template end<0, Dune::Interior_Partition>());
+        offsets[i++] += offsetCounter;
+      }
 
-    std::ofstream outFile(fullfilename);
+      vtkFile.cellOffsets_ = offsets;
 
-    // Write header
-    outFile << "<?xml version=\"1.0\"?>" << std::endl;
-    outFile << "<VTKFile type=\"UnstructuredGrid\" version=\"0.1\" byte_order=\"LittleEndian\">" << std::endl;
-    outFile << "  <UnstructuredGrid>" << std::endl;
-    outFile << "    <Piece NumberOfCells=\"" << numElements << "\" NumberOfPoints=\"" << displacementConfiguration.size() << "\">" << std::endl;
+      std::vector<int> cellTypes(totalNumElements);
+      i = 0;
+      for (const auto& element : elements(gridView, Dune::Partitions::interior))
+      {
+        if (element.type().isQuadrilateral())
+          cellTypes[i++] = (vtkOrder==2) ? 23 : 9;
 
-    // Write vertex coordinates
-    outFile << "      <Points>" << std::endl;
-    outFile << "        <DataArray type=\"Float32\" Name=\"Coordinates\" NumberOfComponents=\"3\" format=\"ascii\">" << std::endl;
-    for (size_t i=0; i<displacementConfiguration.size(); i++)
-      outFile << "          " << displacementConfiguration[i] << std::endl;
-    outFile << "        </DataArray>" << std::endl;
-    outFile << "      </Points>" << std::endl;
+        if (element.type().isTriangle())
+          cellTypes[i++] = (vtkOrder==2) ? 22 : 5;
 
-    // Write elements
-    outFile << "      <Cells>" << std::endl;
+        if (element.type().isHexahedron())
+          cellTypes[i++] = (vtkOrder==2) ? 25 : 12;
+      }
+      vtkFile.cellTypes_ = cellTypes;
 
-    outFile << "         <DataArray type=\"Int32\" Name=\"connectivity\" NumberOfComponents=\"1\" format=\"ascii\">" << std::endl;
+      // Z coordinate for better visualization of wrinkles
+      std::vector<double> zCoord(points.size());
+      for (size_t i=0; i<configuration.size(); i++)
+        zCoord[i] = configuration[i][_0].globalCoordinates()[2];
 
-    for (const auto& element : elements(gridView, Dune::Partitions::interior))
-    {
-      localView.bind(element);
+      vtkFile.zCoord_ = zCoord;
 
-      outFile << "          ";
-      if (element.type().isQuadrilateral())
+      // The three director fields
+      for (size_t i=0; i<3; i++)
       {
-        outFile << localView.index(0) << " ";
-        outFile << localView.index(2) << " ";
-        outFile << localView.index(8) << " ";
-        outFile << localView.index(6) << " ";
-
-        outFile << localView.index(1) << " ";
-        outFile << localView.index(5) << " ";
-        outFile << localView.index(7) << " ";
-        outFile << localView.index(3) << " ";
-
-        outFile << std::endl;
+        vtkFile.directors_[i].resize(configuration.size());
+        for (size_t j=0; j<configuration.size(); j++)
+          vtkFile.directors_[i][j] = configuration[j][_1].director(i);
       }
-    }
-    outFile << "         </DataArray>" << std::endl;
 
-    outFile << "         <DataArray type=\"Int32\" Name=\"offsets\" NumberOfComponents=\"1\" format=\"ascii\">" << std::endl;
-    size_t offsetCounter = 0;
-    for (const auto& element : elements(gridView))
-    {
-      outFile << "          ";
-      if (element.type().isQuadrilateral())
-        offsetCounter += 8;
-      outFile << offsetCounter << std::endl;
+      // Actually write the VTK file to disk
+      vtkFile.write(filename);
     }
-    outFile << "         </DataArray>" << std::endl;
 
-    outFile << "         <DataArray type=\"UInt8\" Name=\"types\" NumberOfComponents=\"1\" format=\"ascii\">" << std::endl;
-    for (const auto& element : elements(gridView))
+    /** \brief Write a Cosserat configuration as VTK file
+     *
+     * The microrotation field will be represented as the three director fields
+     *
+     * \param directorBasis The basis that will be used to represent director fields
+     * \param order The polynomial order that the data will have in the VTK file
+     */
+    template <typename DisplacementFunction, typename OrientationFunction>
+    static void write(const GridView& gridView,
+                      const DisplacementFunction& displacement,
+                      const OrientationFunction& orientation,
+                      int order,
+                      const std::string& filename)
     {
-      outFile << "          ";
-      if (element.type().isQuadrilateral())
-        outFile << "23" << std::endl;
-    }
-    outFile << "         </DataArray>" << std::endl;
-
-    outFile << "      </Cells>" << std::endl;
-#if 0
-    // Point data
-    outFile << "      <PointData Scalars=\"zCoord\" Vectors=\"director0\">" << std::endl;
-
-    // Z coordinate for better visualization of wrinkles
-    outFile << "        <DataArray type=\"Float32\" Name=\"zCoord\" NumberOfComponents=\"1\" format=\"ascii\">" << std::endl;
-    for (size_t i=0; i<configuration.size(); i++)
-      outFile << "          " << configuration[i].r[2] << std::endl;
-    outFile << "        </DataArray>" << std::endl;
+      using namespace Dune;
+
+      // Create a writer object
+#if DUNE_VERSION_GTE(DUNE_VTK, 2, 10)
+      auto vtkWriter = Vtk::UnstructuredGridWriter(Vtk::LagrangeDataCollector(gridView,order));
+#else
+      Vtk::LagrangeDataCollector<GridView> dataCollector(gridView, order);
+      auto vtkWriter = VtkUnstructuredGridWriter(dataCollector);
+#endif
 
-    // The three director fields
-    for (size_t i=0; i<3; i++)
-    {
-      outFile << "        <DataArray type=\"Float32\" Name=\"director" << i <<"\" NumberOfComponents=\"3\" format=\"ascii\">" << std::endl;
-      for (size_t j=0; j<configuration.size(); j++)
-        outFile << "          " << configuration[j].q.director(i) << std::endl;
-      outFile << "        </DataArray>" << std::endl;
+      // Attach the displacement field
+      vtkWriter.addPointData(displacement, Dune::VTK::FieldInfo("displacement", Dune::VTK::FieldInfo::Type::vector, 3));
+
+      // Attach the director fields
+      // This lambda takes a unit quaternion and extracts one column
+      // of the corresponding rotation matrix
+      auto directorExtractor = [](FieldVector<double,4> q,int columnNumber) -> FieldVector<double,3>
+                               {
+                                 FieldMatrix<double,3,3> matrix;
+                                 Rotation<double,3>(q).matrix(matrix);
+                                 FieldVector<double,3> column;
+                                 for (size_t i=0; i<3; ++i)
+                                   column[i] = matrix[i][columnNumber];
+                                 return column;
+                               };
+
+      auto director0Function = Functions::makeComposedGridFunction(std::bind(directorExtractor,std::placeholders::_1,0),
+                                                                   orientation);
+
+      auto director1Function = Functions::makeComposedGridFunction(std::bind(directorExtractor,std::placeholders::_1,1),
+                                                                   orientation);
+
+      auto director2Function = Functions::makeComposedGridFunction(std::bind(directorExtractor,std::placeholders::_1,2),
+                                                                   orientation);
+
+      vtkWriter.addPointData(director0Function, Dune::VTK::FieldInfo("director0", Dune::VTK::FieldInfo::Type::vector, 3));
+      vtkWriter.addPointData(director1Function, Dune::VTK::FieldInfo("director1", Dune::VTK::FieldInfo::Type::vector, 3));
+      vtkWriter.addPointData(director2Function, Dune::VTK::FieldInfo("director2", Dune::VTK::FieldInfo::Type::vector, 3));
+
+      // Write the file
+      vtkWriter.write(filename);
     }
 
-    outFile << "      </PointData>" << std::endl;
-#endif
-    // Write footer
-    outFile << "    </Piece>" << std::endl;
-    outFile << "  </UnstructuredGrid>" << std::endl;
-    outFile << "</VTKFile>" << std::endl;
-
-  }
+  };
 
-};
+}  // namespace Dune::GFE
 
 #endif

dune/gfe/densities/CMakeLists.txt

@@ -3,7 +3,9 @@ install(FILES
         bulkcosseratdensity.hh
         chiralskyrmiondensity.hh
         cosseratshelldensity.hh
+        cosseratvolumeloaddensity.hh
         duneelasticitydensity.hh
         harmonicdensity.hh
         localdensity.hh
+        planarcosseratshelldensity.hh
         DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/dune/gfe/densities)

dune/gfe/densities/bulkcosseratdensity.hh

 #ifndef DUNE_GFE_DENSITIES_BULKCOSSERATDENSITY_HH
 #define DUNE_GFE_DENSITIES_BULKCOSSERATDENSITY_HH
 
-#define DONT_USE_CURL
-#define CURVATURE_WITH_WRYNESS
-
 #include <dune/common/fmatrix.hh>
 
 #include <dune/gfe/cosseratstrain.hh>
@@ -13,16 +10,19 @@
 #include <dune/gfe/spaces/realtuple.hh>
 #include <dune/gfe/spaces/rotation.hh>
 
-namespace Dune::GFE {
+namespace Dune::GFE
+{
 
-  template<class Position, class field_type>
+  template<class Element, class field_type>
   class BulkCosseratDensity final
-    : public GFE::LocalDensity<Position, GFE::ProductManifold<RealTuple<field_type,3>,Rotation<field_type,3> > >
+    : public GFE::LocalDensity<Element, GFE::ProductManifold<RealTuple<field_type,3>,Rotation<field_type,3> > >
   {
   private:
+    using LocalCoordinate = typename Element::Geometry::LocalCoordinate;
+
     // BulkCosseratDensity only works for 3d->3d problems
-    static_assert(Position::size()==3);
-    static const int gridDim = Position::size();
+    static_assert(LocalCoordinate::size()==3);
+    static const int gridDim = LocalCoordinate::size();
     static const int embeddedDim = Rotation<field_type,gridDim>::embeddedDim;
 
     // The target space with 'adouble' as the number type
@@ -41,8 +41,48 @@ namespace Dune::GFE {
       return result;
     }
 
+    /** \brief Compute the wryness tensor from the microrotation and its derivative
+     *
+     * The wryness tensor implemented here is defined in Eq. (19) in:
+     *
+     *   M. Bîrsan, I.D. Ghiba, R.J. Martin, P. Neff: "Refined dimensional reduction
+     *   for isotropic elastic Cosserat shells with initial curvature",
+     *   Mathematics and Mechanics of Solids, 24(12), 2019
+     *   DOI: https://doi.org/10.1177/1081286519856061
+     *
+     */
+    static FieldMatrix<field_type,3,3> wryness(const FieldMatrix<field_type,gridDim,gridDim>& R,
+                                               const Tensor3<field_type,3,3,3>& DR)
+    {
+      FieldMatrix<field_type,3,3> dRx1(0);  //Derivative of R with respect to x1
+      FieldMatrix<field_type,3,3> dRx2(0);  //Derivative of R with respect to x2
+      FieldMatrix<field_type,3,3> dRx3(0);  //Derivative of R with respect to x3
+      for (int i=0; i<3; i++)
+        for (int j=0; j<3; j++) {
+          dRx1[i][j] = DR[i][j][0];
+          dRx2[i][j] = DR[i][j][1];
+          dRx3[i][j] = DR[i][j][2];
+        }
+
+      FieldMatrix<field_type,3,3> wrynessTensor(0);
+
+      auto axialVectorx1 = SkewMatrix<field_type,3>(transpose(R)*dRx1).axial();
+      auto axialVectorx2 = SkewMatrix<field_type,3>(transpose(R)*dRx2).axial();
+      auto axialVectorx3 = SkewMatrix<field_type,3>(transpose(R)*dRx3).axial();
+      for (int i=0; i<3; i++) {
+        wrynessTensor[i][0] = axialVectorx1[i];
+        wrynessTensor[i][1] = axialVectorx2[i];
+        wrynessTensor[i][2] = axialVectorx3[i];
+      }
+
+      return wrynessTensor;
+    }
+
   public:
 
+    /** \brief The different types of curvature tensors */
+    enum CurvatureType {Norm, Curl, Wryness};
+
     /** \brief Constructor with a set of material parameters
      * \param parameters The material parameters
      */
@@ -57,6 +97,17 @@ namespace Dune::GFE {
       // Length scale parameter
       L_c_ = parameters.template get<double>("L_c");
 
+      // Curvature type
+      const auto& curvatureString = parameters.template get<std::string>("curvatureType");
+      if (curvatureString=="norm")
+        curvatureType_ = CurvatureType::Norm;
+      else if (curvatureString=="curl")
+        curvatureType_ = CurvatureType::Curl;
+      else if (curvatureString=="wryness")
+        curvatureType_ = CurvatureType::Wryness;
+      else
+        DUNE_THROW(NotImplemented, "Curvature type '" << curvatureString << "' is not implemented!");
+
       // Curvature exponent
       q_ = parameters.template get<double>("q");
 
@@ -68,8 +119,8 @@ namespace Dune::GFE {
 
     /** \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])
+    BulkCosseratDensity(double mu, double lambda, double mu_c, double L_c, CurvatureType curvatureType, double q, const std::array<double,3>& b)
+      : mu_(mu), lambda_(lambda), mu_c_(mu_c), L_c_(L_c), curvatureType_(curvatureType), q_(q), b1_(b[0]), b2_(b[1]), b3_(b[2])
     {}
 
     /** \brief The energy \f$ W_{mp}(\overline{U}) \f$, as written in
@@ -82,52 +133,26 @@ namespace Dune::GFE {
       for (int i=0; i<gridDim; i++)
         UMinus1[i][i] -= 1;
 
-      return mu_ * GFE::sym(UMinus1).frobenius_norm2()
-             + mu_c_ * GFE::skew(UMinus1).frobenius_norm2()
 #ifdef QUADRATIC_2006
-             + (mu_*lambda_)/(2*mu_ + lambda_) * GFE::traceSquared(UMinus1);  // GFE::traceSquared(UMinus1) = GFE::traceSquared(GFE::sym(UMinus1))
+      field_type materialFactor = (mu_*lambda_)/(2*mu_ + lambda_);
 #else
-             + lambda_/2 * GFE::traceSquared(UMinus1);  // GFE::traceSquared(UMinus1) = GFE::traceSquared(GFE::sym(UMinus1))
+      field_type materialFactor = lambda_/2;
 #endif
-    }
-
-    field_type curvatureWithWryness(const FieldMatrix<field_type,gridDim,gridDim>& R, const Tensor3<field_type,3,3,3>& DR) const
-    {
-      // construct Wryness tensor \Gamma as in "Refined dimensional reduction for isotropic elastic Cosserat shells with initial curvature"
-      FieldMatrix<field_type,3,3> dRx1(0); //Derivative of R with respect to x1
-      FieldMatrix<field_type,3,3> dRx2(0); //Derivative of R with respect to x2
-      FieldMatrix<field_type,3,3> dRx3(0); //Derivative of R with respect to x3
-      for (int i=0; i<3; i++)
-        for (int j=0; j<3; j++) {
-          dRx1[i][j] = DR[i][j][0];
-          dRx2[i][j] = DR[i][j][1];
-          dRx3[i][j] = DR[i][j][2];
-        }
-
-      FieldMatrix<field_type,3,3> wryness(0);
 
-      auto axialVectorx1 = SkewMatrix<field_type,3>(transpose(R)*dRx1).axial();
-      auto axialVectorx2 = SkewMatrix<field_type,3>(transpose(R)*dRx2).axial();
-      auto axialVectorx3 = SkewMatrix<field_type,3>(transpose(R)*dRx3).axial();
-      for (int i=0; i<3; i++) {
-        wryness[i][0] = axialVectorx1[i];
-        wryness[i][1] = axialVectorx2[i];
-        wryness[i][2] = axialVectorx3[i];
-      }
-
-      // The choice for the Frobenius norm here is b1=b2=1 and b3 = 1/3
-      return mu_ * L_c_ * L_c_ * (b1_ * GFE::dev(GFE::sym(wryness)).frobenius_norm2()
-                                  + b2_ * GFE::skew(wryness).frobenius_norm2() + b3_ * GFE::traceSquared(wryness));
+      // In various papers the last term is traceSquared(sym(UMinus1)),
+      // but UMinus1 is symmetric and hence
+      //   traceSquared(UMinus1) = traceSquared(sym(UMinus1)).
+      return mu_ * GFE::sym(UMinus1).frobenius_norm2()
+             + mu_c_ * GFE::skew(UMinus1).frobenius_norm2()
+             + materialFactor * GFE::traceSquared(UMinus1);
     }
 
     /** \brief Evaluate the density
      */
-    virtual field_type operator() (const Position& x,
+    virtual field_type operator() (const LocalCoordinate& x,
                                    const typename GFE::ProductManifold<RealTuple<field_type,3>,Rotation<field_type,3> >::CoordinateType& value,
                                    const FieldMatrix<field_type,7,gridDim>& derivative) const override
     {
-      using namespace Dune::Indices;
-
       field_type strainEnergyDensity = 0;
 
       /////////////////////////////////////////////////////////
@@ -136,48 +161,87 @@ namespace Dune::GFE {
 
       Rotation<field_type,3> rotationValue(FieldVector<field_type,4>{value[3], value[4], value[5], value[6]});
 
-      FieldMatrix<field_type,3,3> deformationDerivative = {derivative[0], derivative[1], derivative[2]};
-      FieldMatrix<field_type,4,3> orientationDerivative = {derivative[3], derivative[4], derivative[5], derivative[6]};
+      FieldMatrix<field_type,3,gridDim> deformationDerivative = {derivative[0], derivative[1], derivative[2]};
+      FieldMatrix<field_type,4,gridDim> orientationDerivative = {derivative[3], derivative[4], derivative[5], derivative[6]};
 
       /////////////////////////////////////////////////////////
       // compute U, the Cosserat strain
       /////////////////////////////////////////////////////////
 
-      FieldMatrix<field_type,gridDim,gridDim> R;
+      FieldMatrix<field_type,3,3> R;
       rotationValue.matrix(R);
 
-      GFE::CosseratStrain<field_type,gridDim,gridDim> U(deformationDerivative,R);
+      GFE::CosseratStrain<field_type,3,gridDim> U(deformationDerivative,R);
 
       /////////////////////////////////////////////////////////////////////
       //  Transfer the derivative of the rotation into matrix coordinates
       ////////////////////////////////////////////////////////////////////
 
       Tensor3<field_type,3,3,gridDim> DR = rotationValue.quaternionTangentToMatrixTangent(orientationDerivative);
-      strainEnergyDensity = quadraticEnergy(U);
-
-    #ifdef CURVATURE_WITH_WRYNESS
-      strainEnergyDensity += curvatureWithWryness(R,DR);
-    #else
 
-    #ifdef DONT_USE_CURL
-      auto argument = DR;
-    #else
-      auto argument = curl(DR);
-    #endif
+      // Compute the bending energy density
+      strainEnergyDensity = quadraticEnergy(U);
 
-      auto norm = (b1_ * GFE::dev(GFE::sym(argument)).frobenius_norm2()
-                   + b2_ * GFE::skew(argument).frobenius_norm2() + b3_ * GFE::traceSquared(argument));
+      // Compute the curvature energy density
+      field_type norm = 0;
+
+      switch (curvatureType_)
+      {
+      case CurvatureType::Norm :
+      {
+        /* This is a simplified version of the curvature energy that appears in
+         *
+         *   P. Neff: "A geometrically exact Cosserat shell-model including size effects,
+         *   avoiding degeneracy in the thin shell limit. Part I: Formal dimensional reduction
+         *   for elastic plates and existence of minimizers for positive Cosserat couple modulus"
+         *   Continuum Mech. Thermodyn. (2004) 16: 577–628
+         *   DOI: 10.1007/s00161-004-0182-4
+         *
+         * The curvature energy given there is much more complicated, but it
+         * does use the full norm of DR to define the curvature.
+         */
+        norm = DR.frobenius_norm2();
+        break;
+      }
+      case CurvatureType::Curl :
+      {
+        auto curvatureTensor = curl(DR);
+
+        // The choice for the Frobenius norm here is b1=b2=1 and b3 = 1/3
+        norm = (b1_ * GFE::dev(GFE::sym(curvatureTensor)).frobenius_norm2()
+                + b2_ * GFE::skew(curvatureTensor).frobenius_norm2()
+                + b3_ * GFE::traceSquared(curvatureTensor));
+        break;
+      }
+      case CurvatureType::Wryness :
+      {
+        auto curvatureTensor = wryness(R,DR);
+
+        // The choice for the Frobenius norm here is b1=b2=1 and b3 = 1/3
+        norm = (b1_ * GFE::dev(GFE::sym(curvatureTensor)).frobenius_norm2()
+                + b2_ * GFE::skew(curvatureTensor).frobenius_norm2()
+                + b3_ * GFE::traceSquared(curvatureTensor));
+        break;
+      }
+      }
 
       strainEnergyDensity += mu_ * std::pow(L_c_ * L_c_ * norm,q_/2.0);
-    #endif
 
       return strainEnergyDensity;
     }
 
     // Construct a copy of this density but using 'adouble' as the number type
-    virtual std::unique_ptr<LocalDensity<Position,ATargetSpace> > makeActiveDensity() const
+    virtual std::unique_ptr<LocalDensity<Element,ATargetSpace> > makeActiveDensity() const override
     {
-      return std::make_unique<BulkCosseratDensity<Position,adouble> >(mu_, lambda_, mu_c_, L_c_, q_, std::array<double,3>{b1_, b2_, b3_});
+      // curvatureType_ is a local enum type, and therefore its type changes
+      // together with the type of the surrounding class.
+      auto activeCurvatureType = static_cast<typename BulkCosseratDensity<Element,adouble>::CurvatureType>(curvatureType_);
+
+      auto result = std::make_unique<BulkCosseratDensity<Element,adouble> >(mu_, lambda_, mu_c_, L_c_, activeCurvatureType, q_, std::array<double,3>{b1_, b2_, b3_});
+
+      if (this->elementOrIntersection_)
+        result->bind(*this->elementOrIntersection_);
+      return result;
     }
 
     /** \brief The density depends on the microrotation value, but not on the deformation
@@ -194,15 +258,18 @@ namespace Dune::GFE {
       return true;
     }
 
-    /** \brief Lame constants */
+    /** \brief Lamé constants */
     double mu_, lambda_;
 
-    /** \brief Cosserat couple modulus, preferably 0 */
+    /** \brief Cosserat couple modulus */
     double mu_c_;
 
     /** \brief Length scale parameter */
     double L_c_;
 
+    /** \brief The precise way to compute the curvature tensor */
+    enum CurvatureType curvatureType_;
+
     /** \brief Curvature exponent */
     double q_;
 

dune/gfe/densities/chiralskyrmiondensity.hh

@@ -7,6 +7,7 @@
 #include <dune/common/parametertree.hh>
 
 #include <dune/gfe/densities/localdensity.hh>
+#include <dune/gfe/spaces/unitvector.hh>
 
 namespace Dune::GFE
 {
@@ -16,13 +17,14 @@ namespace Dune::GFE
    * The energy is discussed in:
    * - Christof Melcher, "Chiral skyrmions in the plane", Proc. of the Royal Society, online DOI DOI: 10.1098/rspa.2014.0394
    */
-  template<class Position, class field_type>
+  template<class ElementOrIntersection, class field_type>
   class ChiralSkyrmionDensity
-    : public GFE::LocalDensity<Position,UnitVector<field_type,3> >
+    : public GFE::LocalDensity<ElementOrIntersection,UnitVector<field_type,3> >
   {
     // various useful types
+    using LocalCoordinate = typename ElementOrIntersection::Geometry::LocalCoordinate;
     using TargetSpace = UnitVector<field_type,3>;
-    using Derivative = FieldMatrix<field_type, TargetSpace::embeddedDim, Position::size()>;
+    using Derivative = FieldMatrix<field_type, TargetSpace::embeddedDim, LocalCoordinate::size ()>;
 
     using ATargetSpace = typename TargetSpace::template rebind<adouble>::other;
 
@@ -47,7 +49,7 @@ namespace Dune::GFE
      * \param value The deformation at the current position
      * \param derivative The derivative of the deformation at the current position
      */
-    virtual field_type operator() (const Position& x,
+    virtual field_type operator() (const LocalCoordinate& x,
                                    const typename TargetSpace::CoordinateType& value,
                                    const Derivative& derivative) const override
     {
@@ -77,9 +79,12 @@ namespace Dune::GFE
     }
 
     // Construct a copy of this density but using 'adouble' as the number type
-    virtual std::unique_ptr<LocalDensity<Position,ATargetSpace> > makeActiveDensity() const
+    virtual std::unique_ptr<LocalDensity<ElementOrIntersection,ATargetSpace> > makeActiveDensity() const
     {
-      return std::make_unique<ChiralSkyrmionDensity<Position,adouble> >(h_, kappa_);
+      auto result = std::make_unique<ChiralSkyrmionDensity<ElementOrIntersection,adouble> >(h_, kappa_);
+      if (this->elementOrIntersection_)
+        result->bind(*this->elementOrIntersection_);
+      return result;
     }
 
     /** \brief The density depends on the value */

dune/gfe/densities/cosseratshelldensity.hh

 #ifndef DUNE_GFE_DENSITIES_COSSERATSHELLDENSITY_HH
 #define DUNE_GFE_DENSITIES_COSSERATSHELLDENSITY_HH
 
+#include <adolc/adouble.h>
+
 #include <dune/common/fmatrix.hh>
 #include <dune/common/parametertree.hh>
 
+#if HAVE_DUNE_CURVEDGRID
+#include <dune/curvedgrid/curvedgrid.hh>
+#endif
+
+#include <dune/matrix-vector/crossproduct.hh>
+
+#include <dune/gfe/linearalgebra.hh>
 #include <dune/gfe/tensor3.hh>
+#include <dune/gfe/densities/localdensity.hh>
 #include <dune/gfe/spaces/productmanifold.hh>
 #include <dune/gfe/spaces/realtuple.hh>
 #include <dune/gfe/spaces/rotation.hh>
@@ -24,21 +34,62 @@ namespace Dune::GFE
    *
    * (and earlier papers by Neff and Bîrsan).
    *
-   * \tparam Position Type for the points where the density will be evaluated
+   * \tparam ElementOrIntersection The domain of the density.
+   *   Can be either a grid element (i.e., a codimension-0 Entity)
+   *   or an Intersection.
    * \tparam field_type Type used for numbers
    */
-  template<class Position, class field_type>
-  class CosseratShellDensity
+  template<class ElementOrIntersection, class field_type>
+  class CosseratShellDensity final
+    : public LocalDensity<ElementOrIntersection,
+          ProductManifold<RealTuple<field_type,3>, Rotation<field_type,3> > >
   {
     constexpr static int dimWorld = 3;
     constexpr static int domainDim = 2;
 
-    // TODO: This is likely not a good long-term solution
-    static_assert(Position::size()==dimWorld, "Position must be in world-coordinates.");
+    using Geometry = typename ElementOrIntersection::Geometry;
+    using ctype = typename Geometry::ctype;
+    using LocalCoordinate = typename Geometry::LocalCoordinate;
 
     using RigidBodyMotion = ProductManifold<RealTuple<field_type,dimWorld>, Rotation<field_type,dimWorld> >;
+    // The target space with 'adouble' as the number type
+    using ARigidBodyMotion = ProductManifold<RealTuple<adouble,3>,Rotation<adouble,3> >;
+
+    // The type used for derivatives of geometric FE functions
     using Derivative = FieldMatrix<field_type,RigidBodyMotion::embeddedDim,domainDim>;
 
+#if HAVE_DUNE_CURVEDGEOMETRY
+    /** \brief Get the derivative of the element normal vector
+     *
+     * Currently only CurvedGrid provides this information, and therefore
+     * we have this specialization.
+     *
+     * The first argument is a geometry implementation rather than
+     * the geometry itself because it makes the template specialization simpler.
+     */
+    template <class LF, class LG, bool useInterpolation>
+    static auto normalGradient(const Dune::Curved::Geometry<ctype,domainDim,dimWorld,LF,LG,useInterpolation>& geometryImpl,
+                               const typename Geometry::LocalCoordinate& pos)
+    {
+      return geometryImpl.normalGradient(pos);
+    }
+#endif
+
+    /** \brief Derivative of the element normal if the geometry is not Dune::Curved::Geometry
+     *
+     * In this case, the element is assumed to be affine, and the derivative
+     * of the normal vector is simply the zero matrix.
+     *
+     * \todo This may lead to wrong results, for example when using
+     * bilinear quadrilateral elements in a 3d world.
+     */
+    template <typename GeometryImpl>
+    static auto normalGradient(const GeometryImpl& geometryImpl,
+                               const typename Geometry::LocalCoordinate& pos)
+    {
+      return Dune::FieldMatrix<ctype,3,3>(0);
+    }
+
     /** \brief Compute the derivative of the rotation (with respect to x), but wrt matrix coordinates
         \param value Value of the gfe function at a certain point
         \param derivative First derivative of the gfe function wrt x at that point, in quaternion coordinates
@@ -87,23 +138,82 @@ namespace Dune::GFE
       return mu * GFE::sym(S).frobenius_norm2() + mu_c_ * GFE::skew(S).frobenius_norm2() + lambda * 0.5 * GFE::traceSquared(S);
     }
 
+    // For b1 = b2 = 1 and b3 = 1/3, this reduces to S.frobenius_norm2()
     field_type W_curv(const Dune::FieldMatrix<field_type,3,3>& S, double mu) const
     {
       return mu * L_c_ * L_c_ * (b1_ * GFE::dev(Dune::GFE::sym(S)).frobenius_norm2()
                                  + b2_ * GFE::skew(S).frobenius_norm2() + b3_ * GFE::traceSquared(S));
     }
 
+    /*  Sources:
+        Birsan 2019: Derivation of a refined six-parameter shell model, equation (111)
+        Birsan 2021: Alternative derivation of the higher-order constitudtive model for six-parameter elastic shells, equations (119) and (126)
+     */
+    field_type W_Coss(const FieldMatrix<field_type,3,3>& S,
+                      const FieldMatrix<double,3,3>& a,
+                      const Dune::FieldVector<double,3>& n0,
+                      double mu, double lambda) const
+    {
+      return W_Coss_mixt(S,S,a,n0,mu,lambda);
+    }
+
+    field_type W_Coss_mixt(const FieldMatrix<field_type,3,3>& S,
+                           const FieldMatrix<field_type,3,3>& T,
+                           const FieldMatrix<double,3,3>& a,
+                           const FieldVector<double,3>& n0,
+                           double mu, double lambda) const
+    {
+      auto planarPart = W_mixt(a*S,a*T,mu,lambda);
+      FieldVector<field_type, 3> n0S;
+      FieldVector<field_type, 3> n0T;
+      S.mtv(n0, n0S);
+      T.mtv(n0, n0T);
+      field_type normalPart = 2*mu*mu_c_* n0S * n0T /(mu + mu_c_);
+      return planarPart + normalPart;
+    }
 
 
   public:
 
     /** \brief Constructor with a set of material parameters
-     * \param parameters The material parameters
-     * \param shellBoundary The shellBoundary contains the faces where the cosserat energy is assembled
-     * \param curvedGeometryGridFunction The curvedGeometryGridFunction gives the geometry of the shell in stress-free state.
-              When assembling, we deform the intersections using the curvedGeometryGridFunction and then use the deformed geometries.
-     * \param thicknessF The shell thickness parameter, given as a function and evaluated at each quadrature point
-     * \param lameF The Lame parameters, given as a function and evaluated at each quadrature point
+     * \param parameters The material parameters, including the shell thickness
+     */
+    CosseratShellDensity(const Dune::ParameterTree& parameters)
+    {
+      // Make a constant function of the given scalar value for the thickness
+      auto thicknessValue = parameters.template get<double>("thickness");
+
+      thicknessF_ = [thicknessValue](const FieldVector<double,dimWorld>& x){return thicknessValue;};
+
+      // Same for the Lamé parameters
+      double mu = parameters.template get<double>("mu");
+      double lambda = parameters.template get<double>("lambda");
+
+      lameF_ = [mu,lambda](const FieldVector<double,dimWorld>& x)
+               -> FieldVector<double,2>
+               {
+                 return {mu,lambda};
+               };
+
+      // Cosserat couple modulus
+      mu_c_ = parameters.template get<double>("mu_c");
+
+      // Length scale parameter
+      L_c_ = parameters.template get<double>("L_c");
+
+      // Curvature parameters
+      b1_ = parameters.template get<double>("b1");
+      b2_ = parameters.template get<double>("b2");
+      b3_ = parameters.template get<double>("b3");
+
+      // Indicator to use the alternative energy W_Coss from Birsan 2021:
+      useAlternativeEnergyWCoss_ = parameters.template get<bool>("useAlternativeEnergyWCoss", false);
+    }
+
+    /** \brief Constructor with space-dependent thickness and Lamé parameters
+     * \param parameters The constant-in-space material parameters
+     * \param thickness The shell thickness parameter, as a function of 3d space
+     * \param lame The Lamé parameters, as a function of 3d space
      */
     CosseratShellDensity(const Dune::ParameterTree& parameters,
                          const std::function<double(Dune::FieldVector<double,dimWorld>)> thickness,
@@ -121,6 +231,75 @@ namespace Dune::GFE
       b1_ = parameters.template get<double>("b1");
       b2_ = parameters.template get<double>("b2");
       b3_ = parameters.template get<double>("b3");
+
+      // Indicator to use the alternative energy W_Coss from Birsan 2021:
+      useAlternativeEnergyWCoss_ = parameters.template get<bool>("useAlternativeEnergyWCoss", false);
+    }
+
+    /** \brief Constructor with space-dependent thickness and Lamé parameters
+     * \param thickness The shell thickness parameter, as a function of 3d space
+     * \param lame The Lamé parameters, as a function of 3d space
+     */
+    CosseratShellDensity(double mu_c,
+                         double L_c,
+                         double b1, double b2, double b3,
+                         bool useAlternativeEnergyWCoss,
+                         const std::function<double(Dune::FieldVector<double,dimWorld>)> thickness,
+                         const std::function<Dune::FieldVector<double,2>(Dune::FieldVector<double,dimWorld>)> lame)
+      : thicknessF_(thickness),
+      lameF_(lame),
+      mu_c_(mu_c), L_c_(L_c), b1_(b1), b2_(b2), b3_(b3),
+      useAlternativeEnergyWCoss_(useAlternativeEnergyWCoss)
+    {}
+
+    /** \brief Evaluate the density for a given value and first derivative
+     *
+     * \param x The current position
+     * \param value The value of the integrand at x
+     * \param derivative The derivative of the integrand at x
+     */
+    field_type operator() (const LocalCoordinate& x,
+                           const typename RigidBodyMotion::CoordinateType& value,
+                           const Derivative& derivative) const override
+    {
+      // TODO: Using this method via LocalIntegralEnergy will not produce
+      // correct results: LocalIntegralEnergy transforms the 'derivative'
+      // argument from the reference element to the grid element, while
+      // this density expects the derivatives from the reference element.
+      // At least that's how I think it is.  Until I have found the time to check,
+      // let's abort here.
+      DUNE_THROW(NotImplemented, "Check whether this method produces the correct results!");
+
+      ////////////////////////////////
+      //  First fundamental form
+      ////////////////////////////////
+
+      FieldMatrix<ctype,3,3> aCovariant;
+
+      // If dimworld==3, then the first two lines of aCovariant are simply the jacobianTransposed
+      // of the element.  If dimworld<3 (i.e., ==2), we have to explicitly enters 0.0 in the last column.
+      const auto jacobianTransposed = this->elementOrIntersection_->geometry().jacobianTransposed(x);
+
+      for (int i=0; i<2; i++)
+      {
+        for (int j=0; j<dimWorld; j++)
+          aCovariant[i][j] = jacobianTransposed[i][j];
+        for (int j=dimWorld; j<3; j++)
+          aCovariant[i][j] = 0.0;
+      }
+
+      aCovariant[2] = MatrixVector::crossProduct(aCovariant[0], aCovariant[1]);
+      aCovariant[2] /= aCovariant[2].two_norm();
+
+      //////////////////////////////////////////////////////////
+      // Add the local energy density
+      //////////////////////////////////////////////////////////
+
+      return operator()(x,
+                        aCovariant,
+                        normalGradient(this->elementOrIntersection_->geometry().impl(), x),
+                        RigidBodyMotion(value),
+                        derivative);
     }
 
     /** \brief Evaluation with the current position, the deformation function, the deformation gradient, the rotation and the rotation gradient
@@ -130,15 +309,16 @@ namespace Dune::GFE
      * \param value The deformation at the current position
      * \param derivative The derivative of the deformation at the current position
      */
-    field_type operator() (const Position& x,
+    field_type operator() (const LocalCoordinate& x,
                            const FieldMatrix<double,dimWorld,dimWorld>& aCovariant,
                            // TODO: Fix the following type
                            const FieldMatrix<double,dimWorld,dimWorld>& normalGradient,
                            const RigidBodyMotion& value,
                            const Derivative& derivative) const
     {
-      double thickness = thicknessF_(x);
-      auto lameConstants = lameF_(x);
+      auto xGlobal = this->elementOrIntersection_->geometry().global(x);
+      double thickness = thicknessF_(xGlobal);
+      auto lameConstants = lameF_(xGlobal);
       auto mu = lameConstants[0];
       auto lambda = lameConstants[1];
       // TODO: Use structured binding here
@@ -234,10 +414,21 @@ namespace Dune::GFE
       //////////////////////////////////////////////////////////
 
       // The membrane density
-      auto density = (thickness - K*Dune::power(thickness,3) / 12.0) * W_m(Ee, mu, lambda);
-      density += (Dune::power(thickness,3) / 12.0 - K * Dune::power(thickness,5) / 80.0)*W_m(Ee*b + c*Ke, mu, lambda);
-      density += Dune::power(thickness,3) / 6.0 * W_mixt(Ee, c*Ke*b - 2*H*c*Ke, mu, lambda);
-      density += Dune::power(thickness,5) / 80.0 * W_mp( (Ee*b + c*Ke)*b, mu, lambda);
+      field_type density = 0;
+      if (useAlternativeEnergyWCoss_)
+      {
+        density += (thickness - K*Dune::power(thickness,3) / 12.0) * W_Coss(Ee, a, aContravariant[2], mu, lambda)
+                   + (Dune::power(thickness,3) / 12.0 - K * Dune::power(thickness,5) / 80.0)*W_Coss(Ee*b + c*Ke, a, aContravariant[2], mu, lambda)
+                   + Dune::power(thickness,3) / 6.0 * W_Coss_mixt(Ee, c*Ke*b - 2*H*c*Ke, a, aContravariant[2], mu, lambda)
+                   + Dune::power(thickness,5) / 80.0 * W_Coss( (Ee*b + c*Ke)*b, a, aContravariant[2], mu, lambda);
+      }
+      else
+      {
+        density += (thickness - K*Dune::power(thickness,3) / 12.0) * W_m(Ee, mu, lambda);
+        density += (Dune::power(thickness,3) / 12.0 - K * Dune::power(thickness,5) / 80.0)*W_m(Ee*b + c*Ke, mu, lambda);
+        density += Dune::power(thickness,3) / 6.0 * W_mixt(Ee, c*Ke*b - 2*H*c*Ke, mu, lambda);
+        density += Dune::power(thickness,5) / 80.0 * W_mp( (Ee*b + c*Ke)*b, mu, lambda);
+      }
 
       // The bending density
       density += (thickness - K*Dune::power(thickness,3) / 12.0) * W_curv(Ke, mu)
@@ -247,6 +438,41 @@ namespace Dune::GFE
       return density;
     }
 
+    // Construct a copy of this density but using 'adouble' as the number type
+    virtual std::unique_ptr<LocalDensity<ElementOrIntersection,ARigidBodyMotion> > makeActiveDensity() const override
+    {
+      auto result = std::make_unique<CosseratShellDensity<ElementOrIntersection,adouble> >(mu_c_,
+                                                                                           L_c_,
+                                                                                           b1_, b2_, b3_,
+                                                                                           useAlternativeEnergyWCoss_,
+                                                                                           thicknessF_,
+                                                                                           lameF_);
+
+      if (this->elementOrIntersection_)
+        result->bind(*this->elementOrIntersection_);
+      return result;
+    }
+
+    /** \brief Whether the density depends on the 'value' argument
+     *
+     * \return The Cosserat shell density depends on the microrotation value,
+     * but not on the deformation (only on the deformation gradient).
+     */
+    virtual bool dependsOnValue(int factor=-1) const override
+    {
+      return factor != 0;
+    }
+
+    /** \brief Whether the density depends on the 'derivative' argument
+     *
+     * \return Always `true` because the Cosserat shell density depends
+     * on the derivatives of both factors: The deformation and the microrotation.
+     */
+    virtual bool dependsOnDerivative([[maybe_unused]] int factor=-1) const override
+    {
+      return true;
+    }
+
   private:
     /** \brief The shell thickness as a function*/
     std::function<double(Dune::FieldVector<double,dimWorld>)> thicknessF_;
@@ -254,10 +480,7 @@ namespace Dune::GFE
     /** \brief The Lamé-parameters as a function*/
     std::function<Dune::FieldVector<double,2>(Dune::FieldVector<double,dimWorld>)> lameF_;
 
-    /** \brief Lame constants */
-    double mu_, lambda_;
-
-    /** \brief Cosserat couple modulus, preferably 0 */
+    /** \brief Cosserat couple modulus */
     double mu_c_;
 
     /** \brief Length scale parameter */
@@ -266,6 +489,10 @@ namespace Dune::GFE
     /** \brief Curvature parameters */
     double b1_, b2_, b3_;
 
+
+    /** \brief Whether to use the alternative energy W_Coss from Birsan 2021:
+               "Alternative derivation of the higher-order constitutive model for six-parameter elastic shells", equations (119) and (126). */
+    bool useAlternativeEnergyWCoss_;
   };
 
 }  // namespace Dune::GFE

dune/gfe/densities/cosseratvolumeloaddensity.hh

+#ifndef DUNE_GFE_DENSITIES_COSSERATVOLUMELOADDENSITY_HH
+#define DUNE_GFE_DENSITIES_COSSERATVOLUMELOADDENSITY_HH
+
+#include <dune/common/fmatrix.hh>
+
+#include <dune/gfe/densities/localdensity.hh>
+#include <dune/gfe/spaces/productmanifold.hh>
+#include <dune/gfe/spaces/realtuple.hh>
+#include <dune/gfe/spaces/rotation.hh>
+
+namespace Dune::GFE
+{
+  template<class ElementOrIntersection, class field_type>
+  class CosseratVolumeLoadDensity final
+    : public GFE::LocalDensity<ElementOrIntersection, GFE::ProductManifold<RealTuple<field_type,3>,Rotation<field_type,3> > >
+  {
+    using LocalCoordinate = typename ElementOrIntersection::Geometry::LocalCoordinate;
+
+    static constexpr int gridDim = LocalCoordinate::size();
+    static constexpr auto dimworld = ElementOrIntersection::Geometry::coorddimension;
+
+    // The target space with 'adouble' as the number type
+    using ATargetSpace = GFE::ProductManifold<RealTuple<adouble,3>,Rotation<adouble,3> >;
+
+    /** \brief The function implementing a volume load */
+    const std::function<FieldVector<double,3>(FieldVector<double,dimworld>)> volumeLoad_;
+
+  public:
+
+    /** \brief Constructor with a given load density
+     */
+    CosseratVolumeLoadDensity(const std::function<FieldVector<double,3>(FieldVector<double,dimworld>)>& volumeLoad)
+      : volumeLoad_(volumeLoad)
+    {}
+
+    /** \brief Evaluate the density
+     */
+    virtual field_type operator() (const LocalCoordinate& x,
+                                   const typename GFE::ProductManifold<RealTuple<field_type,3>,Rotation<field_type,3> >::CoordinateType& value,
+                                   const FieldMatrix<field_type,7,gridDim>& derivative) const override
+    {
+      field_type density = 0;
+
+      // Value of the volume load density at the current position
+      auto loadVector = volumeLoad_(this->elementOrIntersection_->geometry().global(x));
+
+      // In this implementation, only translational dofs are affected by the volume load
+      for (size_t i=0; i<loadVector.size(); i++)
+        density += loadVector[i] * value[i];
+
+      return density;
+    }
+
+    // Construct a copy of this density but using 'adouble' as the number type
+    virtual std::unique_ptr<LocalDensity<ElementOrIntersection,ATargetSpace> > makeActiveDensity() const override
+    {
+      auto result = std::make_unique<CosseratVolumeLoadDensity<ElementOrIntersection,adouble> >(volumeLoad_);
+
+      if (this->elementOrIntersection_)
+        result->bind(*this->elementOrIntersection_);
+      return result;
+    }
+
+    /** \brief The density depends on the deformation value, but not on the microrotation value
+     */
+    virtual bool dependsOnValue(int factor=-1) const override
+    {
+      return factor==-1 || factor==0;
+    }
+
+    /** \brief The density neither depends on the deformation gradient nor on the microrotation gradient
+     */
+    virtual bool dependsOnDerivative([[maybe_unused]] int factor=-1) const override
+    {
+      return false;
+    }
+  };
+
+} // namespace GFE
+
+#endif   //#ifndef DUNE_GFE_DENSITIES_COSSERATVOLUMELOADDENSITY_HH

dune/gfe/densities/duneelasticitydensity.hh

@@ -5,25 +5,33 @@
 
 #include <dune/gfe/densities/localdensity.hh>
 
+#if DUNE_VERSION_GTE(DUNE_ELASTICITY, 2, 11)
+#include <dune/elasticity/densities/localdensity.hh>
+#else
 #include <dune/elasticity/materials/localdensity.hh>
+#endif
 
 namespace Dune::GFE
 {
 
   /** \brief Adapter for densities from the dune-elasticity module
    *
+   * \tparam ElementOrIntersection The domain of the density.
+   *   Can be either a grid element (i.e., a codimension-0 Entity)
+   *   or an Intersection.
    * \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>
+  template<class ElementOrIntersection, class TargetSpace, std::size_t index=0>
   class DuneElasticityDensity final
-    : public LocalDensity<Position,TargetSpace>
+    : public LocalDensity<ElementOrIntersection,TargetSpace>
   {
-    using ctype = typename Position::value_type;
+    using LocalCoordinate = typename ElementOrIntersection::Geometry::LocalCoordinate;
+    using ctype = typename ElementOrIntersection::Geometry::ctype;
     using field_type = typename TargetSpace::field_type;
 
-    constexpr static auto dim = Position::size();
+    constexpr static auto dim = LocalCoordinate::size();
     constexpr static auto embeddedBlocksize = TargetSpace::EmbeddedTangentVector::dimension;
 
     using ATargetSpace = typename TargetSpace::template rebind<adouble>::other;
@@ -48,7 +56,7 @@ namespace Dune::GFE
      * \param value The deformation at the current position
      * \param derivative The derivative at the current position
      */
-    virtual field_type operator() (const Position& x,
+    virtual field_type operator() (const LocalCoordinate& x,
                                    const typename TargetSpace::CoordinateType& value,
                                    const FieldMatrix<field_type,embeddedBlocksize,dim>& derivative) const override
     {
@@ -65,13 +73,17 @@ namespace Dune::GFE
     }
 
     // Construct a copy of this density but using 'adouble' as the number type
-    virtual std::unique_ptr<LocalDensity<Position,ATargetSpace> > makeActiveDensity() const
+    virtual std::unique_ptr<LocalDensity<ElementOrIntersection,ATargetSpace> > makeActiveDensity() const override
     {
       // 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));
+      auto result = std::make_unique<DuneElasticityDensity<ElementOrIntersection,ATargetSpace,index> >(std::move(activeDensity));
+
+      if (this->elementOrIntersection_)
+        result->bind(*this->elementOrIntersection_);
+      return result;
     }
 
     /** \brief The density does not depend on the value */

dune/gfe/densities/harmonicdensity.hh

@@ -7,13 +7,14 @@
 
 namespace Dune::GFE
 {
-  template<class Position, class TargetSpace>
+  template<class ElementOrIntersection, class TargetSpace>
   class HarmonicDensity final
-    : public LocalDensity<Position,TargetSpace>
+    : public LocalDensity<ElementOrIntersection,TargetSpace>
   {
+    using LocalCoordinate = typename ElementOrIntersection::Geometry::LocalCoordinate;
     using field_type = typename TargetSpace::field_type;
 
-    constexpr static auto dim = Position::size();
+    constexpr static auto dim = LocalCoordinate::size();
     constexpr static auto embeddedBlocksize = TargetSpace::EmbeddedTangentVector::dimension;
 
     using ATargetSpace = typename TargetSpace::template rebind<adouble>::other;
@@ -26,7 +27,7 @@ namespace Dune::GFE
      * \param value The deformation at the current position
      * \param derivative The derivative of the deformation at the current position
      */
-    virtual field_type operator() (const Position& x,
+    virtual field_type operator() (const LocalCoordinate& x,
                                    const typename TargetSpace::CoordinateType& value,
                                    const FieldMatrix<field_type,embeddedBlocksize,dim>& derivative) const override
     {
@@ -35,7 +36,7 @@ namespace Dune::GFE
 
     /** \brief Compute value, first and second derivatives of the density
      */
-    virtual void derivatives(const Position& x,
+    virtual void derivatives(const LocalCoordinate& x,
                              const TargetSpace& value,
                              const FieldMatrix<field_type,embeddedBlocksize,dim>& derivative,
                              field_type& densityValue,
@@ -62,9 +63,12 @@ namespace Dune::GFE
     }
 
     // Construct a copy of this density but using 'adouble' as the number type
-    virtual std::unique_ptr<LocalDensity<Position,ATargetSpace> > makeActiveDensity() const
+    virtual std::unique_ptr<LocalDensity<ElementOrIntersection,ATargetSpace> > makeActiveDensity() const override
     {
-      return std::make_unique<HarmonicDensity<Position,ATargetSpace> >();
+      auto result = std::make_unique<HarmonicDensity<ElementOrIntersection,ATargetSpace> >();
+      if (this->elementOrIntersection_)
+        result->bind(*this->elementOrIntersection_);
+      return result;
     }
 
     /** \brief The density does not depend on the value */