dune/gfe/sumenergy.hh

+#ifndef DUNE_GFE_SUMENERGY_HH
+#define DUNE_GFE_SUMENERGY_HH
+
+#include <vector>
+
+#include <dune/gfe/localenergy.hh>
+#include <dune/gfe/mixedlocalgeodesicfestiffness.hh>
+
+namespace Dune::GFE {
+
+  /**
+    \brief Assembles the a sum of energies for a single element by summing up the energies of each GFE::LocalEnergy.
+
+           This class works similarly to the class Dune::Elasticity::SumEnergy, where Dune::Elasticity::SumEnergy extends
+           Dune::Elasticity::LocalEnergy and Dune::GFE::SumEnergy extends Dune::GFE::LocalEnergy.
+  */
+
+template<class Basis, class... TargetSpaces>
+class SumEnergy
+  : public Dune::GFE::LocalEnergy<Basis, TargetSpaces...>,
+    public MixedLocalGeodesicFEStiffness<Basis, TargetSpaces...>
+    //Inheriting from MixedLocalGeodesicFEStiffness is hack, and will be replaced eventually; once MixedLocalGFEADOLCStiffness
+    //will be removed and its functionality will be included in LocalGeodesicFEADOLCStiffness this is not needed anymore!
+
+{
+
+  using LocalView = typename Basis::LocalView;
+  using GridView = typename LocalView::GridView;
+  using DT = typename GridView::Grid::ctype;
+  using RT = typename Dune::GFE::LocalEnergy<Basis,TargetSpaces...>::RT;
+
+public:
+  
+  /** \brief Default Constructor*/
+  SumEnergy()
+  {}
+
+  void addLocalEnergy(std::shared_ptr<GFE::LocalEnergy<Basis,TargetSpaces...>> localEnergy)
+  {
+    localEnergies_.push_back(localEnergy);
+  }
+
+  RT energy(const typename Basis::LocalView& localView,
+       const std::vector<TargetSpaces>&... localSolution) const
+  {
+    RT sum = 0.;
+    for ( const auto& localEnergy : localEnergies_ )
+      sum += localEnergy->energy( localView, localSolution... );
+
+    return sum;
+  }
+
+private:
+  // vector containing pointers to the local energies
+  std::vector<std::shared_ptr<GFE::LocalEnergy<Basis, TargetSpaces...>>> localEnergies_;
+};
+
+}  // namespace Dune::GFE
+
+#endif   //#ifndef DUNE_GFE_SUMENERGY_HH

dune/gfe/surfacecosseratenergy.hh

 #ifndef DUNE_GFE_SURFACECOSSERATENERGY_HH
 #define DUNE_GFE_SURFACECOSSERATENERGY_HH
 
+#include <dune/common/indices.hh>
 #include <dune/geometry/quadraturerules.hh>
 
 #include <dune/fufem/functions/virtualgridfunction.hh>
@@ -16,172 +17,32 @@
 #include <dune/gfe/tensor3.hh>
 #include <dune/gfe/vertexnormals.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).finiteElement())
-  {
-    return localView.tree().child(iType).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 {
 
-template<class Basis, class TargetSpace, class field_type=double, class GradientRT=double>
+template<class Basis, class... TargetSpaces>
 class SurfaceCosseratEnergy
-: public Dune::GFE::LocalEnergy<Basis,TargetSpace>
+: public Dune::GFE::LocalEnergy<Basis, TargetSpaces...>
 {
- // grid types
-  typedef typename Basis::GridView GridView;
-  typedef typename GridView::ctype ctype;
-  typedef typename Basis::LocalView::Tree::FiniteElement LocalFiniteElement;
-  typedef typename GridView::ctype DT;
-  typedef typename TargetSpace::ctype RT;
-  typedef typename GridView::template Codim<0>::Entity Entity;
-  typedef typename TargetSpace::template rebind<adouble>::other ATargetSpace;
-
-  enum {dim=GridView::dimension};
-  enum {dimworld=GridView::dimensionworld};
+  using GridView = typename Basis::GridView;
+  using DT = typename GridView::ctype ;
+  using RT = typename Dune::GFE::LocalEnergy<Basis, TargetSpaces...>::RT ;
+  using Entity = typename GridView::template Codim<0>::Entity ;
+  using RBM0 = RealTuple<RT,GridView::dimensionworld> ;
+  using RBM1 = Rotation<RT,GridView::dimensionworld> ;
+  using RBM = RigidBodyMotion<RT,GridView::dimensionworld> ;
+
+  enum {dimWorld=GridView::dimensionworld};
   enum {gridDim=GridView::dimension};
   static constexpr int boundaryDim = gridDim - 1;
 
-
-  /** \brief Compute the symmetric part of a matrix A, i.e. \f$ \frac 12 (A + A^T) \f$ */
-  static Dune::FieldMatrix<field_type,dim,dim> sym(const Dune::FieldMatrix<field_type,dim,dim>& A)
-  {
-    Dune::FieldMatrix<field_type,dim,dim> result;
-    for (int i=0; i<dim; i++)
-      for (int j=0; j<dim; j++)
-        result[i][j] = 0.5 * (A[i][j] + A[j][i]);
-    return result;
-  }
-
-  /** \brief Compute the antisymmetric part of a matrix A, i.e. \f$ \frac 12 (A - A^T) \f$ */
-  static Dune::FieldMatrix<field_type,dim,dim> skew(const Dune::FieldMatrix<field_type,dim,dim>& A)
-  {
-    Dune::FieldMatrix<field_type,dim,dim> result;
-    for (int i=0; i<dim; i++)
-      for (int j=0; j<dim; j++)
-        result[i][j] = 0.5 * (A[i][j] - A[j][i]);
-    return result;
-  }
-
-  /** \brief Compute the deviator of a matrix A */
-  static Dune::FieldMatrix<field_type,dim,dim> dev(const Dune::FieldMatrix<field_type,dim,dim>& A)
-  {
-    Dune::FieldMatrix<field_type,dim,dim> result = A;
-    auto t = trace(A);
-    for (int i=0; i<dim; i++)
-      result[i][i] -= t / dim;
-    return result;
-  }
-
-  /** \brief Return the trace of a matrix */
-  template <class T, int N>
-  static T trace(const Dune::FieldMatrix<T,N,N>& A)
-  {
-    T trace = 0;
-    for (int i=0; i<N; i++)
-      trace += A[i][i];
-    return trace;
-  }
-
-  /** \brief Return the square of the trace of a matrix */
-  template <int N>
-  static field_type traceSquared(const Dune::FieldMatrix<field_type,N,N>& A)
-  {
-    field_type trace = 0;
-    for (int i=0; i<N; i++)
-      trace += A[i][i];
-    return trace*trace;
-  }
-
-  template <class T, int N>
-  static T frobeniusProduct(const Dune::FieldMatrix<T,N,N>& A, const Dune::FieldMatrix<T,N,N>& B)
-  {
-    T result(0.0);
-
-    for (int i=0; i<N; i++)
-      for (int j=0; j<N; j++)
-        result += A[i][j] * B[i][j];
-
-    return result;
-  }
-
-  template <int N>
-  static auto dyadicProduct(const Dune::FieldVector<adouble,N>& A, const Dune::FieldVector<adouble,N>& B)
-    -> Dune::FieldMatrix<adouble,N,N>
-  {
-    Dune::FieldMatrix<adouble,N,N> result;
-
-    for (int i=0; i<N; i++)
-      for (int j=0; j<N; j++)
-        result[i][j] = A[i]*B[j];
-
-    return result;
-  }
-
-  template <int N>
-  static auto dyadicProduct(const Dune::FieldVector<double,N>& A, const Dune::FieldVector<double,N>& B)
-    -> Dune::FieldMatrix<double,N,N>
-  {
-    Dune::FieldMatrix<double,N,N> result;
-
-    for (int i=0; i<N; i++)
-      for (int j=0; j<N; j++)
-        result[i][j] = A[i]*B[j];
-
-    return result;
-  }
-
-  template <int N>
-  static auto dyadicProduct(const Dune::FieldVector<adouble,N>& A, const Dune::FieldVector<double,N>& B)
-    -> Dune::FieldMatrix<adouble,N,N>
-  {
-    Dune::FieldMatrix<adouble,N,N> result;
-
-    for (int i=0; i<N; i++)
-      for (int j=0; j<N; j++)
-        result[i][j] = A[i]*B[j];
-
-    return result;
-  }
-
-  template <class T, int N>
-  static Dune::FieldMatrix<T,N,N> transpose(const Dune::FieldMatrix<T,N,N>& A)
-  {
-    Dune::FieldMatrix<T,N,N> result;
-
-    for (int i=0; i<N; i++)
-      for (int j=0; j<N; j++)
-        result[i][j] = A[j][i];
-
-    return result;
-  }
-
-
   /** \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
       \param DR First derivative of the gfe function wrt x at that point, in matrix coordinates
    */
-  static void computeDR(const RigidBodyMotion<field_type,dimworld>& value,
-                        const Dune::FieldMatrix<field_type,7,boundaryDim>& derivative,
-                        Tensor3<field_type,3,3,boundaryDim>& DR)
+  static void computeDR(const RBM& value,
+                        const Dune::FieldMatrix<RT,RBM::embeddedDim,boundaryDim>& derivative,
+                        Tensor3<RT,dimWorld,dimWorld,boundaryDim>& derivative_rotation)
   {
     // The LocalGFEFunction class gives us the derivatives of the orientation variable,
     // but as a map into quaternion space.  To obtain matrix coordinates we use the
@@ -192,15 +53,15 @@ class SurfaceCosseratEnergy
     // corresponding orthogonal matrix, we need to invert the i and j indices
     //
     // So, DR[i][j][k] contains \partial R_ij / \partial k
-    Tensor3<field_type,3 , 3, 4> dd_dq;
+    Tensor3<RT, dimWorld, dimWorld, 4> dd_dq;
     value.q.getFirstDerivativesOfDirectors(dd_dq);
 
-    DR = field_type(0);
-    for (int i=0; i<3; i++)
-      for (int j=0; j<3; j++)
+    derivative_rotation = RT(0);
+    for (int i=0; i<dimWorld; i++)
+      for (int j=0; j<dimWorld; j++)
         for (int k=0; k<boundaryDim; k++)
           for (int l=0; l<4; l++)
-            DR[i][j][k] += dd_dq[j][i][l] * derivative[l+3][k];
+            derivative_rotation[i][j][k] += dd_dq[j][i][l] * derivative[l+3][k];
   }
 
 
@@ -211,11 +72,11 @@ public:
    * \param parameters The material parameters
    */
   SurfaceCosseratEnergy(const Dune::ParameterTree& parameters,
-    const std::vector<UnitVector<double,3> >& vertexNormals,
+    const std::vector<UnitVector<double,dimWorld> >& vertexNormals,
     const BoundaryPatch<GridView>* shellBoundary,
-    const std::unordered_map<typename GridView::Grid::GlobalIdSet::IdType,Dune::MultiLinearGeometry<double, dim-1, dim>>& geometriesOnShellBoundary,
-    const std::function<double(Dune::FieldVector<double,dim>)> thicknessF,
-    const std::function<Dune::FieldVector<double,2>(Dune::FieldVector<double,dim>)> lameF)
+    const std::unordered_map<typename GridView::Grid::GlobalIdSet::IdType,Dune::MultiLinearGeometry<double, dimWorld-1, dimWorld>>& geometriesOnShellBoundary,
+    const std::function<double(Dune::FieldVector<double,dimWorld>)> thicknessF,
+    const std::function<Dune::FieldVector<double,2>(Dune::FieldVector<double,dimWorld>)> lameF)
   : shellBoundary_(shellBoundary),
     vertexNormals_(vertexNormals),
     geometriesOnShellBoundary_(geometriesOnShellBoundary),
@@ -235,36 +96,66 @@ public:
   }
 
 RT energy(const typename Basis::LocalView& localView,
-       const std::vector<RigidBodyMotion<field_type,dim> >& localSolution) const
-{
+          const std::vector<TargetSpaces>&... localSolutions) const
+{ 
+  static_assert(sizeof...(TargetSpaces) == 2, "SurfaceCosseratEnergy needs exactly two TargetSpace!");
+
+  using namespace Dune::Indices;
+  using TargetSpace0 = typename std::tuple_element<0, std::tuple<TargetSpaces...> >::type;
+  const std::vector<TargetSpace0>& localSolution0 = std::get<0>(std::forward_as_tuple(localSolutions...));
+  using TargetSpace1 = typename std::tuple_element<1, std::tuple<TargetSpaces...> >::type;
+  const std::vector<TargetSpace1>& localSolution1 = std::get<1>(std::forward_as_tuple(localSolutions...));
+
   // The element geometry
   auto element = localView.element();
-
-  // The set of shape functions on this element
-  const auto& localFiniteElement = localView.tree().finiteElement();
+  auto gridView = localView.globalBasis().gridView();
 
   ////////////////////////////////////////////////////////////////////////////////////
   //  Construct a linear (i.e., non-constant!) normal field on each element
   ////////////////////////////////////////////////////////////////////////////////////
-  auto gridView = localView.globalBasis().gridView();
-
-  assert(vertexNormals_.size() == gridView.indexSet().size(gridDim));
-  std::vector<UnitVector<double,3> > cornerNormals(element.subEntities(gridDim));
-  for (size_t i=0; i<cornerNormals.size(); i++)
-    cornerNormals[i] = vertexNormals_[gridView.indexSet().subIndex(element,i,gridDim)];
-
   typedef typename Dune::PQkLocalFiniteElementCache<DT, double, gridDim, 1> P1FiniteElementCache;
   typedef typename P1FiniteElementCache::FiniteElementType P1LocalFiniteElement;
   //it.type()?
   P1FiniteElementCache p1FiniteElementCache;
   const auto& p1LocalFiniteElement = p1FiniteElementCache.get(element.type());
+
+  assert(vertexNormals_.size() == gridView.indexSet().size(gridDim));
+  std::vector<UnitVector<double,3> > cornerNormals(element.subEntities(gridDim));
+  for (size_t i=0; i<cornerNormals.size(); i++)
+    cornerNormals[i] = vertexNormals_[gridView.indexSet().subIndex(element,i,gridDim)];
   Dune::GFE::LocalProjectedFEFunction<gridDim, DT, P1LocalFiniteElement, UnitVector<double,3> > unitNormals(p1LocalFiniteElement, cornerNormals);
 
   ////////////////////////////////////////////////////////////////////////////////////
   //  Set up the local nonlinear finite element function
   ////////////////////////////////////////////////////////////////////////////////////
-  typedef LocalGeodesicFEFunction<gridDim, DT, LocalFiniteElement, TargetSpace> LocalGFEFunctionType;
-  LocalGFEFunctionType localGeodesicFEFunction(localFiniteElement,localSolution);
+  using namespace Dune::TypeTree::Indices;
+
+  // 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();
+    
+  // to construt a local GFE function, in case they are the shape functions are the same, we can use use one GFE-Function
+#if MIXED_SPACE
+    std::vector<RBM0> localSolutionRBM0(localSolution0.size());
+    std::vector<RBM1> localSolutionRBM1(localSolution1.size());
+    for (int i = 0; i < localSolution0.size(); i++)
+      localSolutionRBM0[i] = localSolution0[i];
+    for (int i = 0; i< localSolution1.size(); i++)
+      localSolutionRBM1[i] = localSolution1[i];
+    typedef LocalGeodesicFEFunction<gridDim, DT, decltype(deformationLocalFiniteElement), RBM0> LocalGFEFunctionType0;
+    typedef LocalGeodesicFEFunction<gridDim, DT, decltype(orientationLocalFiniteElement), RBM1> LocalGFEFunctionType1;
+    LocalGFEFunctionType0 localGeodesicFEFunction0(deformationLocalFiniteElement,localSolutionRBM0);
+    LocalGFEFunctionType1 localGeodesicFEFunction1(orientationLocalFiniteElement,localSolutionRBM1);
+#else
+    std::vector<RBM> localSolutionRBM(localSolution0.size());
+    for (int i = 0; i < localSolution0.size(); i++) {
+      for (int j = 0; j < dimWorld; j++)
+        localSolutionRBM[i].r[j] = localSolution0[i][j];
+      localSolutionRBM[i].q = localSolution1[i];
+    }
+    typedef LocalGeodesicFEFunction<gridDim, DT, decltype(deformationLocalFiniteElement), RBM> LocalGFEFunctionType;
+    LocalGFEFunctionType localGeodesicFEFunction(deformationLocalFiniteElement,localSolutionRBM);
+#endif
 
   RT energy = 0;
 
@@ -276,8 +167,8 @@ RT energy(const typename Basis::LocalView& localView,
     
     auto id = idSet.subId(it.inside(), it.indexInInside(), 1);
     auto boundaryGeometry = geometriesOnShellBoundary_.at(id);
-    auto quadOrder = (it.type().isSimplex()) ? localFiniteElement.localBasis().order()
-                                                  : localFiniteElement.localBasis().order() * boundaryDim;
+    auto quadOrder = (it.type().isSimplex()) ? deformationLocalFiniteElement.localBasis().order()
+                                                  : deformationLocalFiniteElement.localBasis().order() * boundaryDim;
 
     const auto& quad = Dune::QuadratureRules<DT, boundaryDim>::rule(it.type(), quadOrder);
     for (size_t pt=0; pt<quad.size(); pt++) {
@@ -294,47 +185,71 @@ RT energy(const typename Basis::LocalView& localView,
 
       const DT integrationElement = boundaryGeometry.integrationElement(quad[pt].position());
 
-      // The value of the local function
-      RigidBodyMotion<field_type,dim> value = localGeodesicFEFunction.evaluate(quadPos);
-
-      // The derivative of the local function
-      auto derivative = localGeodesicFEFunction.evaluateDerivative(quadPos,value);
+      RBM value;
+      Dune::FieldMatrix<RT, RBM::embeddedDim, dimWorld> derivative;
+      Dune::FieldMatrix<RT, RBM::embeddedDim, boundaryDim> derivative2D;
 
-      Dune::FieldMatrix<RT,7,2> derivative2D;
-      for (int i = 0; i < 7; i++) {
-        for (int j = 0; j < 2; j++) {
-          derivative2D[i][j] = derivative[i][j];
+#if MIXED_SPACE
+      // The value of the local function
+      RBM0 value0 = localGeodesicFEFunction0.evaluate(quadPos);
+      RBM1 value1 = localGeodesicFEFunction1.evaluate(quadPos);
+      // The derivatives of the local functions
+      auto derivative0 = localGeodesicFEFunction0.evaluateDerivative(quadPos,value0);
+      auto derivative1 = localGeodesicFEFunction1.evaluateDerivative(quadPos,value1);
+
+      // Put the value and the derivatives together from the separated values
+      for (int i = 0; i < RBM0::dim; i++)
+          value.r[i] = value0[i];
+      value.q = value1;
+      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];
         }
       }
+#else
+      value = localGeodesicFEFunction.evaluate(quadPos);
+      derivative = localGeodesicFEFunction.evaluateDerivative(quadPos,value);
+      for (int i = 0; i < RBM::embeddedDim; i++)
+        for (int j = 0; j < boundaryDim; j++)
+          derivative2D[i][j] = derivative[i][j];
+#endif
+
       //////////////////////////////////////////////////////////
       //  The rotation and its derivative
       //  Note: we need it in matrix coordinates
       //////////////////////////////////////////////////////////
 
-      Dune::FieldMatrix<field_type,dim,dim> R;
+      Dune::FieldMatrix<RT,dimWorld,dimWorld> R;
       value.q.matrix(R);
-      auto RT = transpose(R);
+      auto rt = Dune::GFE::transpose(R);
 
-      Tensor3<field_type,3,3,boundaryDim> DR;
-      computeDR(value, derivative2D, DR);
+      Tensor3<RT,dimWorld,dimWorld,boundaryDim> derivative_rotation;
+      computeDR(value, derivative2D, derivative_rotation);
 
       //////////////////////////////////////////////////////////
       //  Fundamental forms and curvature
       //////////////////////////////////////////////////////////
 
       // First fundamental form
-      Dune::FieldMatrix<double,3,3> aCovariant;
+      Dune::FieldMatrix<double,dimWorld,dimWorld> 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.
+      // 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 = boundaryGeometry.jacobianTransposed(quad[pt].position());
       // auto jacobianTransposed = geometry.jacobianTransposed(quadPos);
 
       for (int i=0; i<2; i++)
       {
-        for (int j=0; j<dimworld; j++)
+        for (int j=0; j<dimWorld; j++)
           aCovariant[i][j] = jacobianTransposed[i][j];
-        for (int j=dimworld; j<3; j++)
+        for (int j=dimWorld; j<3; j++)
           aCovariant[i][j] = 0.0;
       }
 
@@ -346,7 +261,7 @@ RT energy(const typename Basis::LocalView& localView,
 
       Dune::FieldMatrix<double,3,3> a(0);
       for (int alpha=0; alpha<boundaryDim; alpha++)
-        a += dyadicProduct(aCovariant[alpha], aContravariant[alpha]);
+        a += Dune::GFE::dyadicProduct(aCovariant[alpha], aContravariant[alpha]);
 
       auto a00 = aCovariant[0] * aCovariant[0];
       auto a01 = aCovariant[0] * aCovariant[1];
@@ -360,7 +275,7 @@ RT energy(const typename Basis::LocalView& localView,
 
       for (int alpha=0; alpha<2; alpha++)
         for (int beta=0; beta<2; beta++)
-          c += sqrt(aScalar) * eps[alpha][beta] * dyadicProduct(aContravariant[alpha], aContravariant[beta]);
+          c += aScalar * eps[alpha][beta] * Dune::GFE::dyadicProduct(aContravariant[alpha], aContravariant[beta]);
 
       // Second fundamental form
       // The derivative of the normal field
@@ -372,42 +287,42 @@ RT energy(const typename Basis::LocalView& localView,
         Dune::FieldVector<double,3> vec;
         for (int i=0; i<3; i++)
           vec[i] = normalDerivative[i][alpha];
-        b -= dyadicProduct(vec, aContravariant[alpha]);
+        b -= Dune::GFE::dyadicProduct(vec, aContravariant[alpha]);
       }
 
       // Gauss curvature
       auto K = b.determinant();
 
       // Mean curvatue
-      auto H = 0.5 * trace(b);
+      auto H = 0.5 * Dune::GFE::trace(b);
 
       //////////////////////////////////////////////////////////
       //  Strain tensors
       //////////////////////////////////////////////////////////
 
       // Elastic shell strain
-      Dune::FieldMatrix<field_type,3,3> Ee(0);
-      Dune::FieldMatrix<field_type,3,3> grad_s_m(0);
+      Dune::FieldMatrix<RT,3,3> Ee(0);
+      Dune::FieldMatrix<RT,3,3> grad_s_m(0);
       for (int alpha=0; alpha<boundaryDim; alpha++)
       {
-        Dune::FieldVector<field_type,3> vec;
+        Dune::FieldVector<RT,3> vec;
         for (int i=0; i<3; i++)
           vec[i] = derivative[i][alpha];
-        grad_s_m += dyadicProduct(vec, aContravariant[alpha]);
+        grad_s_m += Dune::GFE::dyadicProduct(vec, aContravariant[alpha]);
       }
 
-      Ee = RT * grad_s_m - a;
+      Ee = rt * grad_s_m - a;
 
       // Elastic shell bending-curvature strain
-      Dune::FieldMatrix<field_type,3,3> Ke(0);
+      Dune::FieldMatrix<RT,3,3> Ke(0);
       for (int alpha=0; alpha<boundaryDim; alpha++)
       {
-        Dune::FieldMatrix<field_type,3,3> tmp;
+        Dune::FieldMatrix<RT,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 += dyadicProduct(SkewMatrix<field_type,3>(tmp2).axial(), aContravariant[alpha]);
+            tmp[i][j] = derivative_rotation[i][j][alpha];
+        auto tmp2 = rt * tmp;
+        Ke += Dune::GFE::dyadicProduct(SkewMatrix<RT,3>(tmp2).axial(), aContravariant[alpha]);
       }
 
       //////////////////////////////////////////////////////////
@@ -433,26 +348,27 @@ RT energy(const typename Basis::LocalView& localView,
   return energy;
 }
 
-  RT W_m(const Dune::FieldMatrix<field_type,3,3>& S, double mu, double lambda) const
+  RT W_m(const Dune::FieldMatrix<RT,3,3>& S, double mu, double lambda) const
   {
     return W_mixt(S,S, mu, lambda);
   }
 
-  RT W_mixt(const Dune::FieldMatrix<field_type,3,3>& S, const Dune::FieldMatrix<field_type,3,3>& T, double mu, double lambda) const
+  RT W_mixt(const Dune::FieldMatrix<RT,3,3>& S, const Dune::FieldMatrix<RT,3,3>& T, double mu, double lambda) const
   {
-    return mu * frobeniusProduct(sym(S), sym(T))
-         + mu_c_ * frobeniusProduct(skew(S), skew(T))
-         + lambda * mu / (lambda + 2*mu) * trace(S) * trace(T);
+    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, double mu, double lambda) const
+  RT W_mp(const Dune::FieldMatrix<RT,3,3>& S, double mu, double lambda) const
   {
-    return mu * sym(S).frobenius_norm2() + mu_c_ * skew(S).frobenius_norm2() + lambda * 0.5 * traceSquared(S);
+    return mu * Dune::GFE::sym(S).frobenius_norm2() + mu_c_ * Dune::GFE::skew(S).frobenius_norm2() + lambda * 0.5 * Dune::GFE::traceSquared(S);
   }
 
-  RT W_curv(const Dune::FieldMatrix<field_type,3,3>& S, double mu) const
+  RT W_curv(const Dune::FieldMatrix<RT,3,3>& S, double mu) const
   {
-    return mu * L_c_ * L_c_ * (b1_ * dev(sym(S)).frobenius_norm2() + b2_ * skew(S).frobenius_norm2() + b3_ * traceSquared(S));
+    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));
   }
 
 private:
@@ -460,16 +376,16 @@ private:
   const BoundaryPatch<GridView>* shellBoundary_;
 
   /** \brief Stress-free geometries of the shell elements*/
-  const std::unordered_map<typename GridView::Grid::GlobalIdSet::IdType, Dune::MultiLinearGeometry<double, dim-1, dim>> geometriesOnShellBoundary_;
+  const std::unordered_map<typename GridView::Grid::GlobalIdSet::IdType, Dune::MultiLinearGeometry<double, dimWorld-1, dimWorld>> geometriesOnShellBoundary_;
 
   /** \brief The normal vectors at the grid vertices. They are used to compute the reference surface curvature. */
   std::vector<UnitVector<double,3> > vertexNormals_;
 
   /** \brief The shell thickness as a function*/
-  std::function<double(Dune::FieldVector<double,dim>)> thicknessF_;
+  std::function<double(Dune::FieldVector<double,dimWorld>)> thicknessF_;
 
   /** \brief The Lamé-parameters as a function*/
-  std::function<Dune::FieldVector<double,2>(Dune::FieldVector<double,dim>)> lameF_;
+  std::function<Dune::FieldVector<double,2>(Dune::FieldVector<double,dimWorld>)> lameF_;
 
   /** \brief Lame constants */
   double mu_, lambda_;
@@ -484,5 +400,6 @@ private:
   double b1_, b2_, b3_;
 
 };
+}  // namespace Dune::GFE
 
 #endif   //#ifndef DUNE_GFE_SURFACECOSSERATENERGY_HH

dune/gfe/vtkreader.hh

-// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
-// vi: set et ts=4 sw=2 sts=2:
-#ifndef DUNE_GFE_VTKREADER_HH
-#define DUNE_GFE_VTKREADER_HH
-
-#include <memory>
-
-#include <dune/gfe/vtkfile.hh>
-
-/** \brief Read a grid from a VTK file
- */
-template <class GridType>
-class VTKReader
-{
-public:
-
-  /** \brief Read a grid from a VTK file */
-  static std::unique_ptr<GridType> read(std::string filename)
-  {
-    constexpr auto dimworld = GridType::dimensionworld;
-
-    Dune::GFE::VTKFile vtkFile;
-
-    // Read test file from disk
-    vtkFile.read(filename);
-
-    Dune::GridFactory<GridType> factory;
-
-    for (const auto& v : vtkFile.points_)
-    {
-      // use the first dimworld components as vertex coordinate; discard the rest
-      Dune::FieldVector<typename GridType::ctype, dimworld> pos;
-      for (int i=0; i<dimworld; i++)
-        pos[i] = v[i];
-      factory.insertVertex(pos);
-    }
-
-    for (size_t i=0; i<vtkFile.cellConnectivity_.size(); i+=3)
-    {
-      factory.insertElement(Dune::GeometryTypes::triangle, {vtkFile.cellConnectivity_[i],
-                                       vtkFile.cellConnectivity_[i+1],
-                                       vtkFile.cellConnectivity_[i+2]});
-
-    }
-
-    return std::unique_ptr<GridType>(factory.createGrid());
-  }
-
-};
-
-#endif

src/film-on-substrate.parset → problems/film-on-substrate.parset

@@ -9,13 +9,13 @@ structuredGrid = true
 lower = 0 0 0
 upper = 200 100 200
 
-elements = 10 5 5
+elements = 4 2 2
 
 # Number of grid levels, all elements containing surfaceshell grid vertices will get adaptively refined
-numLevels = 2
+numLevels = 1
 
 # When starting from a file, the stress-free configuration of the surfaceShell is read from a file, this file needs to match the *finest* grid level!
-startFromFile = true
+startFromFile = false
 pathToGridDeformationFile = ./
 gridDeformationFile = deformation
 
@@ -26,7 +26,7 @@ gridDeformation="[1.3*x[0], x[1], x[2]]"
 #  Boundary values
 #############################################
 
-problem = cantilever
+dirichletValues = identity-dirichlet-values
 
 ###  Python predicate specifying all Dirichlet grid vertices
 # x is the vertex coordinate
@@ -34,7 +34,7 @@ dirichletVerticesPredicate = "x[0] < 0.01"
 
 ###  Python predicate specifying all surfaceshell grid vertices, elements conataining these vertices will get adaptively refined
 # x is the vertex coordinate
-surfaceShellVerticesPredicate = "x[2] > 199.99"
+surfaceShellVerticesPredicate = "x[2] > 199.99 and x[0] > 49.99 and x[0] < 150.01"
 
 ###  Python predicate specifying all Neumann grid vertices
 # x is the vertex coordinate
@@ -50,14 +50,31 @@ initialDeformation = "x"
 #  Solver parameters
 #############################################
 
+# Inner solver, cholmod or multigrid
+solvertype = multigrid
+
 # Number of homotopy steps for the Dirichlet boundary conditions
 numHomotopySteps = 1
 
-# Tolerance of the trust region solver
+# Tolerance of the solver
 tolerance = 1e-3
 
-# Max number of steps of the trust region solver
-maxTrustRegionSteps = 20
+# Max number of solver steps
+maxSolverSteps = 1
+
+# Measure convergence
+instrumented = 0
+
+#############################################
+#  Solver parameters specific for proximal newton solver using cholmod
+#############################################
+
+# initial regularization parameter
+initialRegularization = 1000000
+
+#############################################
+#  Solver parameters specific for trust-region solver using multigrid solver 
+#############################################
 
 trustRegionScaling = 1 1 1 0.01 0.01 0.01
 
@@ -85,14 +102,11 @@ mgTolerance = 1e-5
 # Tolerance of the base grid solver
 baseTolerance = 1e-8
 
-# Measure convergence
-instrumented = 0
-
 ############################
 #   Material parameters
 ############################
 
-energy = stvenantkirchhoff
+energy = mooneyrivlin
 
 ## For the Wriggers L-shape example
 [materialParameters]
@@ -108,13 +122,7 @@ surfaceShellParameters = surface-shell-parameters-1-3
 mu_c = 0
 
 # Length scale parameter
-L_c = 0.2 
-
-# Curvature exponent
-q = 2
-
-# Shear correction factor
-kappa = 1
+L_c = 0.2
 
 b1 = 1
 b2 = 1
@@ -140,4 +148,4 @@ mooneyrivlin_energy = log # log, square or ciarlet; different ways to compute th
 # log: Generalized Rivlin model or polynomial hyperelastic model, using  0.5*mooneyrivlin_k*log(det(∇φ)) as the volume-preserving penalty term
 # square: Generalized Rivlin model or polynomial hyperelastic model, using mooneyrivlin_k*(det(∇φ)-1)² as the volume-preserving penalty term
 
-[]
+[]
\ No newline at end of file

src/cantilever-dirichlet-values.py → problems/identity-dirichlet-values.py

@@ -5,13 +5,8 @@ class DirichletValues:
         self.homotopyParameter = homotopyParameter
 
     def deformation(self, x):
-        # Dirichlet b.c. simply clamp the shell in the reference configuration
-        out = x
-
-        return out
-
+        return x
 
     def orientation(self, x):
-        rotation = [[1,0,0], [0, 1, 0], [0, 0, 1]]
+        rotation = [[1,0,0], [0,1,0], [0,0,1]]
         return rotation
-

src/CMakeLists.txt

@@ -5,7 +5,7 @@ set(programs compute-disc-error
              rod3d)
 #            rodobstacle)
 
-if (dune-parmg_FOUND)
+if (dune-parmg_FOUND AND ${DUNE_ELASTICITY_VERSION} VERSION_GREATER_EQUAL 2.8)
 	set(programs film-on-substrate ${programs})
 endif()
 foreach(_program ${programs})

src/cosserat-continuum.cc

@@ -46,7 +46,6 @@
 #include <dune/gfe/nonplanarcosseratshellenergy.hh>
 #include <dune/gfe/cosseratvtkwriter.hh>
 #include <dune/gfe/cosseratvtkreader.hh>
-#include <dune/gfe/vtkreader.hh>
 #include <dune/gfe/geodesicfeassembler.hh>
 #include <dune/gfe/riemanniantrsolver.hh>
 #include <dune/gfe/vertexnormals.hh>
@@ -54,6 +53,10 @@
 #include <dune/gfe/mixedgfeassembler.hh>
 #include <dune/gfe/mixedriemanniantrsolver.hh>
 
+#if HAVE_DUNE_VTK
+#include <dune/vtk/vtkreader.hh>
+#endif
+
 // grid dimension
 const int dim = 2;
 const int dimworld = 2;
@@ -78,45 +81,6 @@ const int blocksize = TargetSpace::TangentVector::dimension;
 
 using namespace Dune;
 
-/** \brief A constant vector-valued function, for simple Neumann boundary values */
-struct NeumannFunction
-    : public Dune::VirtualFunction<FieldVector<double,dimworld>, FieldVector<double,3> >
-{
-    NeumannFunction(const FieldVector<double,3> values,
-                    double homotopyParameter)
-    : values_(values),
-      homotopyParameter_(homotopyParameter)
-    {}
-
-    void evaluate(const FieldVector<double, dimworld>& x, FieldVector<double,3>& out) const {
-        out = 0;
-        out.axpy(homotopyParameter_, values_);
-    }
-
-    FieldVector<double,3> values_;
-    double homotopyParameter_;
-};
-
-/** \brief A constant vector-valued function, for simple volume loads */
-struct VolumeLoad
-    : public Dune::VirtualFunction<FieldVector<double,dimworld>, FieldVector<double,3> >
-{
-    VolumeLoad(const FieldVector<double,3> values,
-               double homotopyParameter)
-    : values_(values),
-      homotopyParameter_(homotopyParameter)
-    {}
-
-    void evaluate(const FieldVector<double, dimworld>& x, FieldVector<double,3>& out) const {
-        out = 0;
-        out.axpy(homotopyParameter_, values_);
-    }
-
-    FieldVector<double,3> values_;
-    double homotopyParameter_;
-};
-
-
 int main (int argc, char *argv[]) try
 {
     // initialize MPI, finalize is done automatically on exit
@@ -203,7 +167,11 @@ int main (int argc, char *argv[]) try
         if (suffix == ".msh")
             grid = std::shared_ptr<GridType>(GmshReader<GridType>::read(path + "/" + gridFile));
         else if (suffix == ".vtu" or suffix == ".vtp")
-            grid = VTKReader<GridType>::read(path + "/" + gridFile);
+#if HAVE_DUNE_VTK
+            grid = VtkReader<GridType>::createGridFromFile(path + "/" + gridFile);
+#else
+            DUNE_THROW(NotImplemented, "Please install dune-vtk for VTK reading support!");
+#endif
     }
 
     grid->globalRefine(numLevels-1);
@@ -218,13 +186,17 @@ int main (int argc, char *argv[]) try
 
     using namespace Dune::Functions::BasisFactory;
 #ifdef MIXED_SPACE
+    const int dimRotation = Rotation<double,dim>::embeddedDim;
     auto compositeBasis = makeBasis(
       gridView,
       composite(
-          lagrange<displacementOrder>(),
-          lagrange<rotationOrder>()
-      )
-    );
+        power<dim>(
+            lagrange<displacementOrder>()
+        ),
+        power<dimRotation>(
+            lagrange<rotationOrder>()
+        )
+    ));
 
     typedef Dune::Functions::LagrangeBasis<GridView,displacementOrder> DeformationFEBasis;
     typedef Dune::Functions::LagrangeBasis<GridView,rotationOrder> OrientationFEBasis;
@@ -427,15 +399,25 @@ int main (int argc, char *argv[]) try
     // ////////////////////////////////////////////////////////////
 
     const ParameterTree& materialParameters = parameterSet.sub("materialParameters");
-    std::shared_ptr<NeumannFunction> neumannFunction;
+    FieldVector<double,3> neumannValues {0,0,0};
     if (parameterSet.hasKey("neumannValues"))
-        neumannFunction = std::make_shared<NeumannFunction>(parameterSet.get<FieldVector<double,3> >("neumannValues"),
-                                                       homotopyParameter);
+        neumannValues = parameterSet.get<FieldVector<double,3> >("neumannValues");
+
+    auto neumannFunction = [&]( FieldVector<double,dim> ) {
+      auto nV = neumannValues;
+      nV *= homotopyParameter;
+      return nV;
+    };
 
-    std::shared_ptr<VolumeLoad> volumeLoad;
+    FieldVector<double,3> volumeLoadValues {0,0,0};
     if (parameterSet.hasKey("volumeLoad"))
-        volumeLoad = std::make_shared<VolumeLoad>(parameterSet.get<FieldVector<double,3> >("volumeLoad"),
-                                                                                          homotopyParameter);
+        volumeLoadValues = parameterSet.get<FieldVector<double,3> >("volumeLoad");
+
+    auto volumeLoad = [&]( FieldVector<double,dim>) {
+      auto vL = volumeLoadValues;
+      vL *= homotopyParameter;
+      return vL;
+    };
 
     if (mpiHelper.rank() == 0) {
         std::cout << "Material parameters:" << std::endl;

src/film-on-substrate.cc

+#define MIXED_SPACE 0
 #include <iostream>
 #include <fstream>
 
@@ -12,71 +13,70 @@
 #include <dune/fufem/utilities/adolcnamespaceinjections.hh>
 
 #include <dune/common/bitsetvector.hh>
+#include <dune/common/indices.hh>
 #include <dune/common/parametertree.hh>
 #include <dune/common/parametertreeparser.hh>
 #include <dune/common/timer.hh>
+#include <dune/common/tuplevector.hh>
 #include <dune/common/version.hh>
 
-#include <dune/grid/uggrid.hh>
-#include <dune/grid/utility/structuredgridfactory.hh>
-
-#include <dune/grid/io/file/gmshreader.hh>
-#include <dune/grid/io/file/vtk.hh>
-
-#if DUNE_VERSION_GTE(DUNE_ELASTICITY, 2, 8)
 #include <dune/elasticity/materials/exphenckydensity.hh>
 #include <dune/elasticity/materials/henckydensity.hh>
 #include <dune/elasticity/materials/mooneyrivlindensity.hh>
 #include <dune/elasticity/materials/neohookedensity.hh>
 #include <dune/elasticity/materials/stvenantkirchhoffdensity.hh>
-#include <dune/elasticity/materials/sumenergy.hh>
-#include <dune/elasticity/materials/localintegralenergy.hh>
-#include <dune/elasticity/materials/neumannenergy.hh>
-#else
-#include <dune/elasticity/materials/exphenckyenergy.hh>
-#include <dune/elasticity/materials/henckyenergy.hh>
-#if !DUNE_VERSION_LT(DUNE_ELASTICITY, 2, 7)
-#include <dune/elasticity/materials/mooneyrivlinenergy.hh>
-#endif
-#include <dune/elasticity/materials/neohookeenergy.hh>
-#include <dune/elasticity/materials/neumannenergy.hh>
-#include <dune/elasticity/materials/sumenergy.hh>
-#include <dune/elasticity/materials/stvenantkirchhoffenergy.hh>
-#endif
-
 
 #include <dune/functions/functionspacebases/interpolate.hh>
 #include <dune/functions/functionspacebases/lagrangebasis.hh>
 #include <dune/functions/gridfunctions/discreteglobalbasisfunction.hh>
+#include <dune/functions/functionspacebases/compositebasis.hh>
 #include <dune/functions/functionspacebases/powerbasis.hh>
 
-#include <dune/fufem/boundarypatch.hh>
 #include <dune/fufem/functiontools/boundarydofs.hh>
-#include <dune/fufem/functiontools/basisinterpolator.hh>
-#include <dune/fufem/functionspacebases/dunefunctionsbasis.hh>
 #include <dune/fufem/dunepython.hh>
 
-#include <dune/gfe/rigidbodymotion.hh>
-#include <dune/gfe/localgeodesicfeadolcstiffness.hh>
+#include <dune/grid/uggrid.hh>
+#include <dune/grid/utility/structuredgridfactory.hh>
+#include <dune/grid/io/file/gmshreader.hh>
+#include <dune/grid/io/file/vtk.hh>
+
 #include <dune/gfe/cosseratvtkwriter.hh>
-#include <dune/gfe/geodesicfeassembler.hh>
-#include <dune/gfe/riemanniantrsolver.hh>
-#include <dune/gfe/vertexnormals.hh>
+#include <dune/gfe/localintegralenergy.hh>
+#include <dune/gfe/mixedgfeassembler.hh>
+#include <dune/gfe/mixedlocalgfeadolcstiffness.hh>
+#include <dune/gfe/neumannenergy.hh>
 #include <dune/gfe/surfacecosseratenergy.hh>
-#include <dune/gfe/sumcosseratenergy.hh>
+#include <dune/gfe/sumenergy.hh>
+#include <dune/gfe/vertexnormals.hh>
+
+#if MIXED_SPACE
+#include <dune/gfe/mixedriemanniantrsolver.hh>
+#else
+#include <dune/gfe/geodesicfeassemblerwrapper.hh>
+#include <dune/gfe/riemannianpnsolver.hh>
+#include <dune/gfe/riemanniantrsolver.hh>
+#include <dune/gfe/rigidbodymotion.hh>
+#endif
+
+#include <dune/istl/multitypeblockvector.hh>
 
 #include <dune/solvers/solvers/iterativesolver.hh>
 #include <dune/solvers/norms/energynorm.hh>
 
-#include <iostream>
-#include <fstream>
-
 // grid dimension
 #ifndef WORLD_DIM
 #  define WORLD_DIM 3
 #endif
 const int dim = WORLD_DIM;
-const int order = 2;
+
+const int targetDim = WORLD_DIM;
+
+const int displacementOrder = 2;
+const int rotationOrder = 2;
+
+#if !MIXED_SPACE
+static_assert(displacementOrder==rotationOrder, "displacement and rotation order do not match!");
+#endif
 
 #if DUNE_VERSION_LT(DUNE_COMMON, 2, 7)
 template<>
@@ -95,27 +95,6 @@ struct Dune::MathematicalConstants<adouble>
 typedef adouble ValueType;
 
 using namespace Dune;
-#if DUNE_VERSION_LT(DUNE_ELASTICITY, 2, 8)
-/** \brief A constant vector-valued function, for simple Neumann boundary values */
-struct NeumannFunction
-    : public Dune::VirtualFunction<FieldVector<double,dim>, FieldVector<double,dim> >
-{
-  NeumannFunction(const FieldVector<double,dim> values,
-                  double homotopyParameter)
-  : values_(values),
-    homotopyParameter_(homotopyParameter)
-  {}
-
-  void evaluate(const FieldVector<double, dim>& x, FieldVector<double,dim>& out) const
-  {
-    out = 0;
-    out.axpy(-homotopyParameter_, values_);
-  }
-
-  FieldVector<double,dim> values_;
-  double homotopyParameter_;
-};
-#endif
 
 int main (int argc, char *argv[]) try
 {
@@ -123,8 +102,14 @@ int main (int argc, char *argv[]) try
   // initialize MPI, finalize is done automatically on exit
   Dune::MPIHelper& mpiHelper = MPIHelper::instance(argc, argv);
 
-  if (mpiHelper.rank()==0)
-    std::cout << "ORDER = " << order << std::endl;
+  if (mpiHelper.rank()==0) {
+    std::cout << "DISPLACEMENTORDER = " << displacementOrder << ", ROTATIONORDER = " << rotationOrder << std::endl;
+#if MIXED_SPACE
+    std::cout << "MIXED_SPACE = 1" << std::endl;
+#else
+    std::cout << "MIXED_SPACE = 0" << std::endl;
+#endif
+  }
 
   // Start Python interpreter
   Python::start();
@@ -135,15 +120,9 @@ int main (int argc, char *argv[]) try
   Python::runStream()
         << std::endl << "import sys"
         << std::endl << "import os"
-        << std::endl << "sys.path.append(os.getcwd() + '/../../src/')"
+        << std::endl << "sys.path.append(os.getcwd() + '/../../problems/')"
         << std::endl;
 
-  typedef BlockVector<FieldVector<double,dim> > SolutionType;
-  typedef RigidBodyMotion<double,dim> TargetSpace;
-  typedef std::vector<TargetSpace> SolutionTypeCosserat;
-
-  const int blocksize = TargetSpace::TangentVector::dimension;
-
   // parse data file
   ParameterTree parameterSet;
 
@@ -155,8 +134,9 @@ int main (int argc, char *argv[]) try
   int numLevels                         = parameterSet.get<int>("numLevels");
   int numHomotopySteps                  = parameterSet.get<int>("numHomotopySteps");
   const double tolerance                = parameterSet.get<double>("tolerance");
-  const int maxTrustRegionSteps         = parameterSet.get<int>("maxTrustRegionSteps");
+  const int maxSolverSteps              = parameterSet.get<int>("maxSolverSteps");
   const double initialTrustRegionRadius = parameterSet.get<double>("initialTrustRegionRadius");
+  const double initialRegularization    = parameterSet.get<double>("initialRegularization");
   const int multigridIterations         = parameterSet.get<int>("numIt");
   const int nu1                         = parameterSet.get<int>("nu1");
   const int nu2                         = parameterSet.get<int>("nu2");
@@ -168,15 +148,16 @@ int main (int argc, char *argv[]) try
   const bool startFromFile              = parameterSet.get<bool>("startFromFile");
   const std::string resultPath          = parameterSet.get("resultPath", "");
 
-  // ///////////////////////////////////////
-  //    Create the grid
-  // ///////////////////////////////////////
+  /////////////////////////////////////////////////////////////
+  //                      CREATE THE GRID
+  /////////////////////////////////////////////////////////////
   typedef UGGrid<dim> GridType;
 
   std::shared_ptr<GridType> grid;
 
   FieldVector<double,dim> lower(0), upper(1);
 
+
   if (parameterSet.get<bool>("structuredGrid")) {
 
     lower = parameterSet.get<FieldVector<double,dim> >("lower");
@@ -217,28 +198,60 @@ int main (int argc, char *argv[]) try
 
     grid->adapt();
 
-    grid->loadBalance();
-
     numLevels--;
   }
 
+  grid->loadBalance();
+
   if (mpiHelper.rank()==0)
     std::cout << "There are " << grid->leafGridView().comm().size() << " processes" << std::endl;
 
   typedef GridType::LeafGridView GridView;
   GridView gridView = grid->leafGridView();
 
-  // FE basis spanning the FE space that we are working in
-  typedef Dune::Functions::LagrangeBasis<GridView, order> FEBasis;
-  FEBasis feBasis(gridView);
+  /////////////////////////////////////////////////////////////
+  //                      DATA TYPES 
+  /////////////////////////////////////////////////////////////
 
-  // dune-fufem-style FE basis for the transition from dune-fufem to dune-functions
-  typedef DuneFunctionsBasis<FEBasis> FufemFEBasis;
-  FufemFEBasis fufemFEBasis(feBasis);
+  using namespace Dune::Indices;
+
+  typedef std::vector<RealTuple<double,dim> > DisplacementVector;
+  typedef std::vector<Rotation<double,dim>> RotationVector;
+  const int dimRotation = Rotation<double,dim>::embeddedDim;
+  typedef TupleVector<DisplacementVector, RotationVector> SolutionType;
+
+  /////////////////////////////////////////////////////////////
+  //                      FUNCTION SPACE
+  /////////////////////////////////////////////////////////////
+
+  using namespace Dune::Functions::BasisFactory;
+  auto compositeBasis = makeBasis(
+    gridView,
+    composite(
+      power<dim>(
+        lagrange<displacementOrder>()
+      ),
+      power<dimRotation>(
+        lagrange<rotationOrder>()
+      )
+  ));
+
+  auto deformationPowerBasis = makeBasis(
+    gridView,
+    power<dim>(
+        lagrange<displacementOrder>()
+  ));
+  typedef Dune::Functions::LagrangeBasis<GridView,displacementOrder> DeformationFEBasis;
+  typedef Dune::Functions::LagrangeBasis<GridView,rotationOrder> OrientationFEBasis;
+  DeformationFEBasis deformationFEBasis(gridView);
+  OrientationFEBasis orientationFEBasis(gridView);
 
-  // /////////////////////////////////////////
-  //   Read Dirichlet values
-  // /////////////////////////////////////////
+  using CompositeBasis = decltype(compositeBasis);
+  using LocalView = typename CompositeBasis::LocalView;
+
+  /////////////////////////////////////////////////////////////
+  //                      BOUNDARY DATA
+  /////////////////////////////////////////////////////////////
 
   BitSetVector<1> dirichletVertices(gridView.size(dim), false);
   BitSetVector<1> neumannVertices(gridView.size(dim), false);
@@ -262,99 +275,64 @@ int main (int argc, char *argv[]) try
   auto neumannBoundary = std::make_shared<BoundaryPatch<GridView>>(gridView, neumannVertices);
   BoundaryPatch<GridView> surfaceShellBoundary(gridView, surfaceShellVertices);
 
-  if (mpiHelper.rank()==0) {
-    std::cout << "Neumann boundary has " << neumannBoundary->numFaces() << " faces\n";
-    std::cout << "Shell boundary has " << surfaceShellBoundary.numFaces() << " faces\n";
-  }
-
-
-  BitSetVector<1> dirichletNodes(feBasis.size(), false);
-#if DUNE_VERSION_LT(DUNE_FUFEM, 2, 7)
-  using FufemFEBasis = DuneFunctionsBasis<FEBasis>;
-  FufemFEBasis fufemFeBasis(feBasis);
-  constructBoundaryDofs(dirichletBoundary,fufemFeBasis,dirichletNodes);
-#else
-  constructBoundaryDofs(dirichletBoundary,feBasis,dirichletNodes);
-#endif
-
-  BitSetVector<1> neumannNodes(feBasis.size(), false);
-#if DUNE_VERSION_LT(DUNE_FUFEM, 2, 7)
-  constructBoundaryDofs(*neumannBoundary,fufemFeBasis,neumannNodes);
-#else
-  constructBoundaryDofs(*neumannBoundary,feBasis,neumannNodes);
-#endif
-
-  BitSetVector<1> surfaceShellNodes(feBasis.size(), false);
-#if DUNE_VERSION_LT(DUNE_FUFEM, 2, 7)
-  constructBoundaryDofs(surfaceShellBoundary,fufemFeBasis,surfaceShellNodes);
-#else
-  constructBoundaryDofs(surfaceShellBoundary,feBasis,surfaceShellNodes);
-#endif
-
-  BitSetVector<blocksize> dirichletDofs(feBasis.size(), false);
-  for (size_t i=0; i<feBasis.size(); i++)
-    for (int j=0; j<blocksize; j++) {
-      if (dirichletNodes[i][0])
-        dirichletDofs[i][j] = true;
+  std::cout << "On rank " << mpiHelper.rank() << ": Dirichlet boundary has " << dirichletBoundary.numFaces() << " faces\n";
+  std::cout << "On rank " << mpiHelper.rank() << ": Neumann boundary has " << neumannBoundary->numFaces() << " faces\n";
+  std::cout << "On rank " << mpiHelper.rank() << ": Shell boundary has " << surfaceShellBoundary.numFaces() << " faces\n";
+
+  BitSetVector<1> dirichletNodes(compositeBasis.size({0}),false);
+  BitSetVector<1> surfaceShellNodes(compositeBasis.size({1}),false);
+  constructBoundaryDofs(dirichletBoundary,deformationFEBasis,dirichletNodes);
+  constructBoundaryDofs(surfaceShellBoundary,orientationFEBasis,surfaceShellNodes);
+
+  //Create BitVector matching the tangential space
+  const int dimRotationTangent = Rotation<double,dim>::TangentVector::dimension;
+  typedef MultiTypeBlockVector<std::vector<FieldVector<double,dim> >, std::vector<FieldVector<double,dimRotationTangent>> > VectorForBit;
+  typedef Solvers::DefaultBitVector_t<VectorForBit> BitVector;
+
+  BitVector dirichletDofs;
+  dirichletDofs[_0].resize(compositeBasis.size({0}));
+  dirichletDofs[_1].resize(compositeBasis.size({1}));
+  for (size_t i = 0; i < compositeBasis.size({0}); i++){
+    for (int j = 0; j < dim; j++){
+      dirichletDofs[_0][i][j] = dirichletNodes[i][0];
     }
-  for (size_t i=0; i<feBasis.size(); i++)
-    for (int j=3; j<blocksize; j++) {
-      if (not surfaceShellNodes[i][0])
-        dirichletDofs[i][j] = true;
+  }
+  for (size_t i = 0; i < compositeBasis.size({1}); i++) {
+    for (int j = 0; j < dimRotationTangent; j++){
+      dirichletDofs[_1][i][j] = (not surfaceShellNodes[i][0] or dirichletNodes[i][0]);
     }
+  }
 
-  // //////////////////////////
-  //   Initial iterate
-  // //////////////////////////
+  /////////////////////////////////////////////////////////////
+  //                      INITIAL DATA
+  /////////////////////////////////////////////////////////////
   
-  SolutionTypeCosserat x(feBasis.size());
+  SolutionType x;
+  x[_0].resize(compositeBasis.size({0}));
+  x[_1].resize(compositeBasis.size({1}));
 
-  BlockVector<FieldVector<double,3> > v;
-  
   //Initial deformation of the underlying substrate
+  BlockVector<FieldVector<double,dim> > displacement(compositeBasis.size({0}));
   lambda = std::string("lambda x: (") + parameterSet.get<std::string>("initialDeformation") + std::string(")");
   PythonFunction<FieldVector<double,dim>, FieldVector<double,dim> > pythonInitialDeformation(Python::evaluate(lambda));
-  ::Functions::interpolate(fufemFEBasis, v, pythonInitialDeformation);
-
-  //Copy over the initial deformation
-  for (size_t i=0; i<x.size(); i++) {
-      x[i].r = v[i];
-  }
+  Dune::Functions::interpolate(deformationPowerBasis, displacement, pythonInitialDeformation);
 
-  ////////////////////////////////////////////////////////
-  //   Main homotopy loop
-  ////////////////////////////////////////////////////////
+  BlockVector<FieldVector<double,dim> > identity(compositeBasis.size({0}));
+  Dune::Functions::interpolate(deformationPowerBasis, identity, [](FieldVector<double,dim> x){ return x; });
 
-  using namespace Dune::Functions::BasisFactory;
-  auto powerBasis = makeBasis(
-    gridView,
-    power<dim>(
-      lagrange<order>(),
-      blockedInterleaved()
-  ));
-
-  BlockVector<FieldVector<double,dim>> identity;
-  Dune::Functions::interpolate(powerBasis, identity, [](FieldVector<double,dim> x){ return x; });
-
-  BlockVector<FieldVector<double,dim> > displacement(v.size());
-  for (int i = 0; i < v.size(); i++) {
-    displacement[i] = v[i];
+  for (int i = 0; i < displacement.size(); i++) {
+    x[_0][i] = displacement[i]; //Copy over the initial deformation to the deformation part of x
+    displacement[i] -= identity[i]; //Subtract identity to get the initial displacement as a function
   }
-  displacement -= identity;
-
-  auto displacementFunction = Dune::Functions::makeDiscreteGlobalBasisFunction<FieldVector<double,dim>>(feBasis, displacement);
-  auto localDisplacementFunction = localFunction(displacementFunction);
-
-  FieldVector<double,dim> neumannValues {0,0,0};
-
-  if (parameterSet.hasKey("neumannValues"))
-    neumannValues = parameterSet.get<FieldVector<double,dim> >("neumannValues");
-  std::cout << "Neumann values: " << neumannValues << std::endl;
-
+  auto displacementFunction = Dune::Functions::makeDiscreteGlobalBasisFunction<FieldVector<double,dim>>(deformationPowerBasis, displacement);
   //  We need to subsample, because VTK cannot natively display real second-order functions
-  SubsamplingVTKWriter<GridView> vtkWriter(gridView, Dune::refinementLevels(order-1));
-  vtkWriter.addVertexData(localDisplacementFunction, VTK::FieldInfo("displacement", VTK::FieldInfo::Type::scalar, dim));
+  SubsamplingVTKWriter<GridView> vtkWriter(gridView, Dune::refinementLevels(displacementOrder-1));
+  vtkWriter.addVertexData(displacementFunction, VTK::FieldInfo("displacement", VTK::FieldInfo::Type::scalar, dim));
   vtkWriter.write(resultPath + "finite-strain_homotopy_" + parameterSet.get<std::string>("energy") + "_0");
+  
+  /////////////////////////////////////////////////////////////
+  //               INITIAL SURFACE SHELL DATA
+  /////////////////////////////////////////////////////////////
 
   typedef MultiLinearGeometry<double, dim-1, dim> ML;
   std::unordered_map<GridType::GlobalIdSet::IdType, ML> geometriesOnShellBoundary;
@@ -363,34 +341,52 @@ int main (int argc, char *argv[]) try
 
   // Read in the grid deformation
   if (startFromFile) {
+    // Create a basis of order 1 in order to deform the geometries on the surface shell boundary
     const std::string pathToGridDeformationFile = parameterSet.get("pathToGridDeformationFile", "");
-
     // for this, we create a basis of order 1 in order to deform the geometries on the surface shell boundary
-    typedef Dune::Functions::LagrangeBasis<GridView, 1> FEBasisOrder1;
-    FEBasisOrder1 feBasisOrder1(gridView);
-  
+    auto feBasisOrder1 = makeBasis(
+      gridView,
+      power<dim>(
+        lagrange<1>(),
+        blockedInterleaved()
+    ));
+    GlobalIndexSet<GridView> globalVertexIndexSet(gridView,dim);
+
+    std::unordered_map<std::string, FieldVector<double,3>> deformationMap;
+    std::string line, displacement, entry;
+    if (mpiHelper.rank() == 0) 
+      std::cout << "Reading in deformation file: " << pathToGridDeformationFile + parameterSet.get<std::string>("gridDeformationFile") << std::endl;
     // Read grid deformation information from the file specified in the parameter set via gridDeformationFile
-    BlockVector<FieldVector<double,3> > gridDeformationFromFile(feBasisOrder1.size());
-    std::string line;
-    std::ifstream file(pathToGridDeformationFile + parameterSet.get<std::string>("gridDeformationFile"));
+    std::ifstream file(pathToGridDeformationFile + parameterSet.get<std::string>("gridDeformationFile"), std::ios::in);
     if (file.is_open()) {
-      size_t i = 0;
       while (std::getline(file, line)) {
         size_t j = 0;
-        std::stringstream entries(line);
-        std::string entry;
-        FieldVector<double,3> coord(0);
-        while(entries >> entry) {
-          coord[j++] = std::stod(entry);
-        }
-        gridDeformationFromFile[i++] = coord;
+        size_t pos = line.find(":");
+        displacement = line.substr(pos + 1);
+        line.erase(pos);
+        std::stringstream entries(displacement);
+        FieldVector<double,3> displacementVector(0);
+        while(entries >> entry)
+          displacementVector[j++] = std::stod(entry);
+        deformationMap.insert({line,displacementVector});
       }
-      if (i != feBasisOrder1.size())
+      if (mpiHelper.rank() == 0)
+        std::cout << "... done: The grid has " << globalVertexIndexSet.size(dim) << " vertices and the defomation file has " << deformationMap.size() << " entries." << std::endl;
+      if (deformationMap.size() != globalVertexIndexSet.size(dim))
         DUNE_THROW(Exception, "Error: Grid and deformation vector do not match!");
       file.close();
     } else {
       DUNE_THROW(Exception, "Error: Could not open the file containing the deformation vector!");
     }
+  
+    BlockVector<FieldVector<double,dim>> gridDeformationFromFile;
+    Dune::Functions::interpolate(feBasisOrder1, gridDeformationFromFile, [](FieldVector<double,dim> x){ return x; });
+
+    for (auto& entry : gridDeformationFromFile) {
+      std::stringstream stream;
+      stream << entry;
+      entry = deformationMap.at(stream.str()); //Look up the deformation for this vertex in the deformationMap
+    }
 
     auto gridDeformationFromFileFunction = Dune::Functions::makeDiscreteGlobalBasisFunction<FieldVector<double,dim>>(feBasisOrder1, gridDeformationFromFile);
     auto localGridDeformationFromFileFunction = localFunction(gridDeformationFromFileFunction);
@@ -434,15 +430,30 @@ int main (int argc, char *argv[]) try
     }
   }
 
+  /////////////////////////////////////////////////////////////
+  //                      MAIN HOMOTOPY LOOP
+  /////////////////////////////////////////////////////////////
+  FieldVector<double,dim> neumannValues {0,0,0};
+  if (parameterSet.hasKey("neumannValues"))
+    neumannValues = parameterSet.get<FieldVector<double,dim> >("neumannValues");
+  std::cout << "Neumann values: " << neumannValues << std::endl;
+
+  // Vertex Normals for the 3D-Part
+  std::vector<UnitVector<double,dim> > vertexNormals(gridView.size(dim));
+  Dune::FieldVector<double,dim> vertexNormalRaw = {0,0,1};
+  for (int i = 0; i< vertexNormals.size(); i++) {
+    UnitVector vertexNormal(vertexNormalRaw);
+    vertexNormals[i] = vertexNormal;
+  }
+  //Function for the Cosserat material parameters
   const ParameterTree& materialParameters = parameterSet.sub("materialParameters");
   Python::Reference surfaceShellClass = Python::import(materialParameters.get<std::string>("surfaceShellParameters"));
   Python::Callable surfaceShellCallable = surfaceShellClass.get("SurfaceShellParameters");
-
   Python::Reference pythonObject = surfaceShellCallable();
   PythonFunction<Dune::FieldVector<double, dim>, double> fThickness(pythonObject.get("thickness"));
   PythonFunction<Dune::FieldVector<double, dim>, Dune::FieldVector<double, 2>> fLame(pythonObject.get("lame"));
 
-  for (int i=0; i<numHomotopySteps; i++)
+  for (int i = 0; i < numHomotopySteps; i++)
   {
     double homotopyParameter = (i+1)*(1.0/numHomotopySteps);
 
@@ -451,19 +462,13 @@ int main (int argc, char *argv[]) try
     // ////////////////////////////////////////////////////////////
 
 
-#if DUNE_VERSION_LT(DUNE_ELASTICITY, 2, 8)
-    std::shared_ptr<NeumannFunction> neumannFunction;
-    neumannFunction = std::make_shared<NeumannFunction>(neumannValues, homotopyParameter);
-#else
       // A constant vector-valued function, for simple Neumann boundary values
     std::shared_ptr<std::function<Dune::FieldVector<double,dim>(Dune::FieldVector<double,dim>)>> neumannFunctionPtr;
     neumannFunctionPtr = std::make_shared<std::function<Dune::FieldVector<double,dim>(Dune::FieldVector<double,dim>)>>([&](FieldVector<double,dim> ) {
-      auto nv = neumannValues;
-      nv *= (-homotopyParameter);
-      return nv;
+      return neumannValues * (-homotopyParameter);
     });
-#endif
 
+    const ParameterTree& materialParameters = parameterSet.sub("materialParameters");
     if (mpiHelper.rank() == 0)
     {
       std::cout << "Homotopy step: " << i << ",    parameter: " << homotopyParameter << std::endl;
@@ -471,9 +476,12 @@ int main (int argc, char *argv[]) try
       materialParameters.report();
     }
 
-    // Assembler using ADOL-C
     std::cout << "Selected energy is: " << parameterSet.get<std::string>("energy") << std::endl;
-#if DUNE_VERSION_GTE(DUNE_ELASTICITY, 2,8)
+
+    // /////////////////////////////////////////////////
+    //   Create the energy functional
+    // /////////////////////////////////////////////////
+
     std::shared_ptr<Elasticity::LocalDensity<dim,ValueType>> elasticDensity;
     if (parameterSet.get<std::string>("energy") == "stvenantkirchhoff")
       elasticDensity = std::make_shared<Elasticity::StVenantKirchhoffDensity<dim,ValueType>>(materialParameters);
@@ -486,164 +494,178 @@ int main (int argc, char *argv[]) try
     if (parameterSet.get<std::string>("energy") == "mooneyrivlin")
       elasticDensity = std::make_shared<Elasticity::MooneyRivlinDensity<dim,ValueType>>(materialParameters);
 
-    if(!elasticDensity)
-      DUNE_THROW(Exception, "Error: Selected energy not available!");
-
-    auto elasticEnergy = std::make_shared<LocalIntegralEnergy<GridView, FEBasis::LocalView::Tree::FiniteElement, ValueType>>(elasticDensity);
-    auto neumannEnergy = std::make_shared<NeumannEnergy<GridView,
-                                                        FEBasis::LocalView::Tree::FiniteElement,
-                                                        ValueType> >(neumannBoundary,neumannFunctionPtr);
-
-    auto elasticAndNeumann = std::make_shared<Dune::SumEnergy<GridView,
-          FEBasis::LocalView::Tree::FiniteElement,
-          ValueType>>(elasticEnergy, neumannEnergy);
-#else
-#if DUNE_VERSION_LT(DUNE_ELASTICITY, 2, 7)
-    std::shared_ptr<LocalFEStiffness<GridView,
-#else
-    std::shared_ptr<Elasticity::LocalEnergy<GridView,
-#endif
-                                     FEBasis::LocalView::Tree::FiniteElement,
-                                     std::vector<Dune::FieldVector<ValueType, dim>> > > elasticEnergy;
-
-    if (parameterSet.get<std::string>("energy") == "stvenantkirchhoff")
-      elasticEnergy = std::make_shared<StVenantKirchhoffEnergy<GridView,
-                                                               FEBasis::LocalView::Tree::FiniteElement,
-                                                               ValueType> >(materialParameters);
-#if !DUNE_VERSION_LT(DUNE_ELASTICITY, 2, 7)
-    if (parameterSet.get<std::string>("energy") == "mooneyrivlin")
-      elasticEnergy = std::make_shared<MooneyRivlinEnergy<GridView,
-                                                               FEBasis::LocalView::Tree::FiniteElement,
-                                                               ValueType> >(materialParameters);
-#endif
-
-    if (parameterSet.get<std::string>("energy") == "neohooke")
-      elasticEnergy = std::make_shared<NeoHookeEnergy<GridView,
-                                                      FEBasis::LocalView::Tree::FiniteElement,
-                                                      ValueType> >(materialParameters);
-
-    if (parameterSet.get<std::string>("energy") == "hencky")
-      elasticEnergy = std::make_shared<HenckyEnergy<GridView,
-                                                    FEBasis::LocalView::Tree::FiniteElement,
-                                                    ValueType> >(materialParameters);
-
-    if (parameterSet.get<std::string>("energy") == "exphencky")
-      elasticEnergy = std::make_shared<ExpHenckyEnergy<GridView,
-                                                       FEBasis::LocalView::Tree::FiniteElement,
-                                                       ValueType> >(materialParameters);
 
-    if(!elasticEnergy)
+    if(!elasticDensity)
       DUNE_THROW(Exception, "Error: Selected energy not available!");
 
-    auto neumannEnergy = std::make_shared<NeumannEnergy<GridView,
-                                                        FEBasis::LocalView::Tree::FiniteElement,
-                                                        ValueType> >(neumannBoundary.get(),neumannFunction.get());
-
-    auto elasticAndNeumann = std::make_shared<SumEnergy<GridView,
-          FEBasis::LocalView::Tree::FiniteElement,
-          ValueType>>(elasticEnergy, neumannEnergy);
-#endif
-
-    using LocalEnergyBase = GFE::LocalEnergy<FEBasis,RigidBodyMotion<adouble, dim> >;
-
-    std::shared_ptr<LocalEnergyBase> surfaceCosseratEnergy;
-    std::vector<UnitVector<double,3> > vertexNormals(gridView.size(3));
-    Dune::FieldVector<double,3> vertexNormalRaw = {0,0,1};
-    for (int i = 0; i< vertexNormals.size(); i++) {
-      UnitVector vertexNormal(vertexNormalRaw);
-      vertexNormals[i] = vertexNormal;
-    }
-
-    surfaceCosseratEnergy = std::make_shared<SurfaceCosseratEnergy<FEBasis,RigidBodyMotion<adouble, dim>, adouble, adouble>>(materialParameters, std::move(vertexNormals), &surfaceShellBoundary, std::move(geometriesOnShellBoundary), fThickness, fLame);
-
-    std::shared_ptr<LocalEnergyBase> totalEnergy;
-    totalEnergy = std::make_shared<GFE::SumCosseratEnergy<FEBasis,RigidBodyMotion<adouble, dim>, adouble>> (elasticAndNeumann, surfaceCosseratEnergy);
-
+    auto elasticEnergy = std::make_shared<GFE::LocalIntegralEnergy<CompositeBasis, RealTuple<ValueType,targetDim>, Rotation<ValueType,dim>>>(elasticDensity);
+    auto neumannEnergy = std::make_shared<GFE::NeumannEnergy<CompositeBasis, RealTuple<ValueType,targetDim>, Rotation<ValueType,dim>>>(neumannBoundary,*neumannFunctionPtr);
+    auto surfaceCosseratEnergy = std::make_shared<GFE::SurfaceCosseratEnergy<CompositeBasis, RealTuple<ValueType,dim>, Rotation<ValueType,dim> >>(materialParameters, std::move(vertexNormals), &surfaceShellBoundary, std::move(geometriesOnShellBoundary), fThickness, fLame);
 
-    LocalGeodesicFEADOLCStiffness<FEBasis,
-                                  TargetSpace> localGFEADOLCStiffness(totalEnergy.get());
+    GFE::SumEnergy<CompositeBasis, RealTuple<ValueType,targetDim>, Rotation<ValueType,targetDim>> sumEnergy;
+    sumEnergy.addLocalEnergy(neumannEnergy);
+    sumEnergy.addLocalEnergy(elasticEnergy);
+    sumEnergy.addLocalEnergy(surfaceCosseratEnergy);
 
-    GeodesicFEAssembler<FEBasis,TargetSpace> assembler(gridView, &localGFEADOLCStiffness);
-
-    // /////////////////////////////////////////////////
-    //   Create a Riemannian trust-region solver
-    // /////////////////////////////////////////////////
-
-    RiemannianTrustRegionSolver<FEBasis,TargetSpace> solver;
-    solver.setup(*grid,
-                 &assembler,
-                 x,
-                 dirichletDofs,
-                 tolerance,
-                 maxTrustRegionSteps,
-                 initialTrustRegionRadius,
-                 multigridIterations,
-                 mgTolerance,
-                 mu, nu1, nu2,
-                 baseIterations,
-                 baseTolerance,
-                 instrumented);
-
-    solver.setScaling(parameterSet.get<FieldVector<double,6> >("trustRegionScaling"));
+    MixedLocalGFEADOLCStiffness<CompositeBasis,
+                                RealTuple<double,dim>,
+                                Rotation<double,dim> > localGFEADOLCStiffness(&sumEnergy);
+    MixedGFEAssembler<CompositeBasis,
+                      RealTuple<double,dim>,
+                      Rotation<double,dim> > mixedAssembler(compositeBasis, &localGFEADOLCStiffness);
 
     ////////////////////////////////////////////////////////
     //   Set Dirichlet values
     ////////////////////////////////////////////////////////
 
-    Python::Reference dirichletValuesClass = Python::import(parameterSet.get<std::string>("problem") + "-dirichlet-values");
+    BitSetVector<RealTuple<double,dim>::TangentVector::dimension> deformationDirichletDofs(dirichletDofs[_0].size(), false);
+    for(int i = 0; i < dirichletDofs[_0].size(); i++)
+      for (int j = 0; j < RealTuple<double,dim>::TangentVector::dimension; j++)
+        deformationDirichletDofs[i][j] = dirichletDofs[_0][i][j];
+    BitSetVector<Rotation<double,dim>::TangentVector::dimension> orientationDirichletDofs(dirichletDofs[_1].size(), false);
+    for(int i = 0; i < dirichletDofs[_1].size(); i++)
+      for (int j = 0; j < Rotation<double,dim>::TangentVector::dimension; j++)
+        orientationDirichletDofs[i][j] = dirichletDofs[_1][i][j];
 
+    Python::Reference dirichletValuesClass = Python::import(parameterSet.get<std::string>("dirichletValues"));
     Python::Callable C = dirichletValuesClass.get("DirichletValues");
 
     // Call a constructor.
     Python::Reference dirichletValuesPythonObject = C(homotopyParameter);
 
     // Extract object member functions as Dune functions
-    PythonFunction<FieldVector<double,dim>, FieldVector<double,3> >   deformationDirichletValues(dirichletValuesPythonObject.get("deformation"));
-    PythonFunction<FieldVector<double,dim>, FieldMatrix<double,3,3> > orientationDirichletValues(dirichletValuesPythonObject.get("orientation"));
-
-    std::vector<FieldVector<double,3> > ddV;
-    ::Functions::interpolate(fufemFEBasis, ddV, deformationDirichletValues, dirichletDofs);
-
-    std::vector<FieldMatrix<double,3,3> > dOV;
-    ::Functions::interpolate(fufemFEBasis, dOV, orientationDirichletValues, dirichletDofs);
-
-    for (size_t j=0; j<x.size(); j++)
-      if (dirichletNodes[j][0])
-      {
-        x[j].r = ddV[j];
-        x[j].q.set(dOV[j]);
+    PythonFunction<FieldVector<double,dim>, FieldVector<double,targetDim> >   deformationDirichletValues(dirichletValuesPythonObject.get("deformation"));
+    PythonFunction<FieldVector<double,dim>, FieldMatrix<double,targetDim,targetDim> > rotationalDirichletValues(dirichletValuesPythonObject.get("orientation"));
+
+    BlockVector<FieldVector<double,targetDim> > ddV;
+    Dune::Functions::interpolate(deformationPowerBasis, ddV, deformationDirichletValues, deformationDirichletDofs);
+
+    BlockVector<FieldMatrix<double,targetDim,targetDim> > dOV;
+    Dune::Functions::interpolate(orientationFEBasis, dOV, rotationalDirichletValues, orientationDirichletDofs);
+
+    for (int i = 0; i < compositeBasis.size({0}); i++)
+      if (dirichletDofs[_0][i][0])
+        x[_0][i] = ddV[i];
+    for (int i = 0; i < compositeBasis.size({1}); i++)
+      if (dirichletDofs[_1][i][0])
+        x[_1][i].set(dOV[i]);
+
+#if !MIXED_SPACE
+    //The MixedRiemannianTrustRegionSolver can treat the Displacement and Orientation Space as separate ones
+    //The RiemannianTrustRegionSolver can only treat the Displacement and Rotation together in a RigidBodyMotion
+    //Therefore, x and the dirichletDofs are converted to a RigidBodyMotion structure, as well as the Hessian and Gradient that are returned by the assembler
+    typedef RigidBodyMotion<double, dim> RBM;
+    std::vector<RBM> xRBM(compositeBasis.size({0}));
+    BitSetVector<RBM::TangentVector::dimension> dirichletDofsRBM(compositeBasis.size({0}), false);
+    for (int i = 0; i < compositeBasis.size({0}); i++) {
+      for (int j = 0; j < dim; j ++) { // Displacement part
+        xRBM[i].r[j] = x[_0][i][j];
+        dirichletDofsRBM[i][j] = dirichletDofs[_0][i][j];
       }
+      xRBM[i].q = x[_1][i]; // Rotation part
+      for (int j = dim; j < RBM::TangentVector::dimension; j ++)
+        dirichletDofsRBM[i][j] = dirichletDofs[_1][i][j-dim];
+    }
+    typedef Dune::GFE::GeodesicFEAssemblerWrapper<CompositeBasis, DeformationFEBasis, RBM, RealTuple<double, dim>, Rotation<double,dim>> GFEAssemblerWrapper;
+    GFEAssemblerWrapper assembler(&mixedAssembler, deformationFEBasis);
+#endif
 
-    // /////////////////////////////////////////////////////
-    //   Solve!
-    // /////////////////////////////////////////////////////
-
-    solver.setInitialIterate(x);
-    solver.solve();
+    ///////////////////////////////////////////////////
+    //   Create the chosen solver and solve!
+    ///////////////////////////////////////////////////
+
+    if (parameterSet.get<std::string>("solvertype") == "multigrid") {
+#if MIXED_SPACE
+      MixedRiemannianTrustRegionSolver<GridType, CompositeBasis, DeformationFEBasis, RealTuple<double,dim>, OrientationFEBasis, Rotation<double,dim>> solver;
+      solver.setup(*grid,
+                   &mixedAssembler,
+                   deformationFEBasis,
+                   orientationFEBasis,
+                   x,
+                   deformationDirichletDofs,
+                   orientationDirichletDofs,
+                   tolerance,
+                   maxSolverSteps,
+                   initialTrustRegionRadius,
+                   multigridIterations,
+                   mgTolerance,
+                   mu, nu1, nu2,
+                   baseIterations,
+                   baseTolerance,
+                   instrumented);
+
+      solver.setScaling(parameterSet.get<FieldVector<double,6> >("trustRegionScaling"));
+      solver.setInitialIterate(x);
+      solver.solve();
+      x = solver.getSol();
+#else
+      RiemannianTrustRegionSolver<DeformationFEBasis, RBM, GFEAssemblerWrapper> solver;
+      solver.setup(*grid,
+                   &assembler,
+                   xRBM,
+                   dirichletDofsRBM,
+                   tolerance,
+                   maxSolverSteps,
+                   initialTrustRegionRadius,
+                   multigridIterations,
+                   mgTolerance,
+                   mu, nu1, nu2,
+                   baseIterations,
+                   baseTolerance,
+                   instrumented);
+
+      solver.setScaling(parameterSet.get<FieldVector<double,6> >("trustRegionScaling"));
+      solver.setInitialIterate(xRBM);
+      solver.solve();
+      xRBM = solver.getSol();
+      for (int i = 0; i < xRBM.size(); i++) {
+        x[_0][i] = xRBM[i].r;
+        x[_1][i] = xRBM[i].q;
+      }
+#endif
+    } else { //parameterSet.get<std::string>("solvertype") == "cholmod"
 
-    x = solver.getSol();
+#if MIXED_SPACE
+    DUNE_THROW(Exception, "Error: There is no MixedRiemannianProximalNewtonSolver!");
+#else
+      RiemannianProximalNewtonSolver<DeformationFEBasis, RBM, GFEAssemblerWrapper> solver;
+      solver.setup(*grid,
+                   &assembler,
+                   xRBM,
+                   dirichletDofsRBM,
+                   tolerance,
+                   maxSolverSteps,
+                   initialRegularization,
+                   instrumented);
+      solver.setInitialIterate(xRBM);
+      solver.solve();
+      xRBM = solver.getSol();
+      for (int i = 0; i < xRBM.size(); i++) {
+        x[_0][i] = xRBM[i].r;
+        x[_1][i] = xRBM[i].q;
+      }
+#endif
+    }
 
     std::cout << "Overall calculation took " << overallTimer.elapsed() << " sec." << std::endl;
 
     /////////////////////////////////
     //   Output result
     /////////////////////////////////
-    for (size_t i=0; i<x.size(); i++)
-      v[i] = x[i].r;
 
     // Compute the displacement
-    BlockVector<FieldVector<double,dim> > displacement(v.size());
-    for (int i = 0; i < v.size(); i++) {
-      displacement[i] = v[i];
+    BlockVector<FieldVector<double,dim> > displacement(compositeBasis.size({0}));
+    for (int i = 0; i < compositeBasis.size({0}); i++) {
+       for (int j = 0; j  < dim; j++) {
+        displacement[i][j] = x[_0][i][j];
+      }
+      displacement[i] -= identity[i];
     }
-    displacement -= identity;
-
-    auto displacementFunction = Dune::Functions::makeDiscreteGlobalBasisFunction<FieldVector<double,dim>>(feBasis, displacement);
-    auto localDisplacementFunction = localFunction(displacementFunction);
+    auto displacementFunction = Dune::Functions::makeDiscreteGlobalBasisFunction<FieldVector<double,dim>>(deformationPowerBasis, displacement);
 
     //  We need to subsample, because VTK cannot natively display real second-order functions
-    SubsamplingVTKWriter<GridView> vtkWriter(gridView, Dune::refinementLevels(order-1));
-    vtkWriter.addVertexData(localDisplacementFunction, VTK::FieldInfo("displacement", VTK::FieldInfo::Type::scalar, dim));
+    SubsamplingVTKWriter<GridView> vtkWriter(gridView, Dune::refinementLevels(displacementOrder-1));
+    vtkWriter.addVertexData(displacementFunction, VTK::FieldInfo("displacement", VTK::FieldInfo::Type::scalar, dim));
     vtkWriter.write(resultPath + "finite-strain_homotopy_" + parameterSet.get<std::string>("energy") + "_" + std::to_string(neumannValues[0]) + "_" + std::to_string(i+1));
   }
 } catch (Exception& e) {

test/CMakeLists.txt

@@ -26,5 +26,12 @@ dune_add_test(SOURCES harmonicmaptest.cc
               TIMEOUT 600
               CMAKE_GUARD MPI_FOUND)
 
+if (${DUNE_ELASTICITY_VERSION} VERSION_GREATER_EQUAL 2.7)
+dune_add_test(SOURCES geodesicfeassemblerwrappertest.cc
+              MPI_RANKS 1 4
+              TIMEOUT 600
+              CMAKE_GUARD MPI_FOUND)
+endif()
+
 # Copy the example grid used for testing into the build dir
 file(COPY grids/irregular-square.msh DESTINATION ${CMAKE_CURRENT_BINARY_DIR}/grids)

test/geodesicfeassemblerwrappertest.cc

+#include <config.h>
+
+#include <array>
+
+// Includes for the ADOL-C automatic differentiation library
+// Need to come before (almost) all others.
+#include <adolc/adouble.h>
+#include <dune/fufem/utilities/adolcnamespaceinjections.hh>
+
+#include <dune/common/typetraits.hh>
+#include <dune/common/bitsetvector.hh>
+
+#include <dune/functions/functionspacebases/compositebasis.hh>
+#include <dune/functions/functionspacebases/interpolate.hh>
+#include <dune/functions/functionspacebases/lagrangebasis.hh>
+#include <dune/functions/functionspacebases/powerbasis.hh>
+#include <dune/functions/gridfunctions/discreteglobalbasisfunction.hh>
+
+#include <dune/fufem/boundarypatch.hh>
+
+#include <dune/grid/utility/structuredgridfactory.hh>
+#include <dune/grid/uggrid.hh>
+
+#include <dune/gfe/cosseratenergystiffness.hh>
+#include <dune/gfe/geodesicfeassembler.hh>
+#include <dune/gfe/harmonicenergy.hh>
+#include <dune/gfe/localgeodesicfeadolcstiffness.hh>
+#include <dune/gfe/localprojectedfefunction.hh>
+#include <dune/gfe/mixedgfeassembler.hh>
+#include <dune/gfe/mixedlocalgfeadolcstiffness.hh>
+#include <dune/gfe/riemanniantrsolver.hh>
+
+#include <dune/gfe/geodesicfeassemblerwrapper.hh>
+
+#include <dune/istl/multitypeblockmatrix.hh>
+#include <dune/istl/multitypeblockvector.hh>
+
+// grid dimension
+const int gridDim = 2;
+
+// target dimension
+const int dim = 3;
+
+//order of the finite element space
+const int displacementOrder = 2;
+const int rotationOrder = 2;
+
+using namespace Dune;
+using namespace Indices;
+
+//differentiation method: ADOL-C
+using ValueType = adouble;
+
+//Types for the mixed space
+using DisplacementVector = std::vector<RealTuple<double,dim>>;
+using RotationVector =  std::vector<Rotation<double,dim>>;
+using Vector = MultiTypeBlockVector<DisplacementVector, RotationVector>;
+const int dimCR = Rotation<double,dim>::TangentVector::dimension; //dimCorrectionRotation = Dimension of the correction for rotations
+using CorrectionType = MultiTypeBlockVector<BlockVector<FieldVector<double,dim> >, BlockVector<FieldVector<double,dimCR>>>;
+
+using MatrixRow0 = MultiTypeBlockVector<BCRSMatrix<FieldMatrix<double,dim,dim>>,  BCRSMatrix<FieldMatrix<double,dim,dimCR>>>;
+using MatrixRow1 = MultiTypeBlockVector<BCRSMatrix<FieldMatrix<double,dimCR,dim>>, BCRSMatrix<FieldMatrix<double,dimCR,dimCR>>>;
+using MatrixType = MultiTypeBlockMatrix<MatrixRow0,MatrixRow1>;
+
+//Types for the Non-mixed space
+using RBM = RigidBodyMotion<double, dim>;
+const static int blocksize = RBM::TangentVector::dimension;
+using CorrectionTypeWrapped = BlockVector<FieldVector<double, blocksize> >;
+using MatrixTypeWrapped = BCRSMatrix<FieldMatrix<double, blocksize, blocksize> >;
+
+int main (int argc, char *argv[])
+{
+  MPIHelper::instance(argc, argv);
+
+  /////////////////////////////////////////////////////////////////////////
+  //    Create the grid
+  /////////////////////////////////////////////////////////////////////////
+  using GridType = UGGrid<gridDim>;
+  auto grid = StructuredGridFactory<GridType>::createCubeGrid({0,0}, {1,1}, {2,2});
+  grid->globalRefine(2);
+  grid->loadBalance();
+
+  using GridView = GridType::LeafGridView;
+  GridView gridView = grid->leafGridView();
+
+  std::function<bool(FieldVector<double,gridDim>)> isNeumann = [](FieldVector<double,gridDim> coordinate) {
+    return coordinate[0] > 0.99;
+  };
+
+  BitSetVector<1> neumannVertices(gridView.size(gridDim), false);
+  const GridView::IndexSet& indexSet = gridView.indexSet();
+  
+  for (auto&& vertex : vertices(gridView))
+      neumannVertices[indexSet.index(vertex)] = isNeumann(vertex.geometry().corner(0));
+
+  BoundaryPatch<GridView> neumannBoundary(gridView, neumannVertices);
+
+
+  /////////////////////////////////////////////////////////////////////////
+  //  Create a composite basis
+  /////////////////////////////////////////////////////////////////////////
+
+  using namespace Functions::BasisFactory;
+
+  auto compositeBasis = makeBasis(
+    gridView,
+    composite(
+      power<dim>(
+        lagrange<displacementOrder>()
+      ),
+      power<dim>(
+        lagrange<rotationOrder>()
+      )
+  ));
+
+  using CompositeBasis = decltype(compositeBasis);
+  using LocalView = typename CompositeBasis::LocalView;
+
+  /////////////////////////////////////////////////////////////////////////
+  //  Create the energy functions with their parameters
+  /////////////////////////////////////////////////////////////////////////
+
+  //Surface-Cosserat-Energy-Parameters
+  ParameterTree parameters;
+  parameters["thickness"] = "1";
+  parameters["mu"] = "2.7191e+4";
+  parameters["lambda"] = "4.4364e+4";
+  parameters["mu_c"] = "0";
+  parameters["L_c"] = "0.01";
+  parameters["q"] = "2";
+  parameters["kappa"] = "1";
+
+
+  FieldVector<double,dim> values_ = {3e4,2e4,1e4};
+  auto neumannFunction = [&](FieldVector<double, gridDim>){
+    return values_;
+  };
+
+  CosseratEnergyLocalStiffness<decltype(compositeBasis), dim,adouble> cosseratEnergy(parameters,
+                                                                     &neumannBoundary,
+                                                                     neumannFunction,
+                                                                     nullptr);
+  MixedLocalGFEADOLCStiffness<CompositeBasis,
+                              RealTuple<double,dim>,
+                              Rotation<double,dim> > mixedLocalGFEADOLCStiffness(&cosseratEnergy);
+  MixedGFEAssembler<CompositeBasis,
+                    RealTuple<double,dim>,
+                    Rotation<double,dim> > mixedAssembler(compositeBasis, &mixedLocalGFEADOLCStiffness);
+
+  using RBM = RigidBodyMotion<double, dim>;
+  using DeformationFEBasis = Functions::LagrangeBasis<GridView,displacementOrder>;
+  DeformationFEBasis deformationFEBasis(gridView);
+  using GFEAssemblerWrapper = GFE::GeodesicFEAssemblerWrapper<CompositeBasis, DeformationFEBasis, RBM, RealTuple<double, dim>, Rotation<double,dim>>;
+  GFEAssemblerWrapper assembler(&mixedAssembler, deformationFEBasis);
+  
+  /////////////////////////////////////////////////////////////////////////
+  //  Prepare the iterate x where we want to assemble - identity in 2D with z = 0
+  /////////////////////////////////////////////////////////////////////////  
+  auto deformationPowerBasis = makeBasis(
+    gridView,
+    power<gridDim>(
+        lagrange<displacementOrder>()
+  ));
+  BlockVector<FieldVector<double,gridDim> > identity(compositeBasis.size({0}));
+  Functions::interpolate(deformationPowerBasis, identity, [](FieldVector<double,gridDim> x){ return x; });
+  BlockVector<FieldVector<double,dim> > initialDeformation(compositeBasis.size({0}));
+  initialDeformation = 0;
+
+  Vector x;
+  x[_0].resize(compositeBasis.size({0}));
+  x[_1].resize(compositeBasis.size({1}));
+  std::vector<RBM> xRBM(compositeBasis.size({0}));
+  for (int i = 0; i < compositeBasis.size({0}); i++) {
+    for (int j = 0; j < gridDim; j++)
+      initialDeformation[i][j] = identity[i][j];
+    x[_0][i] = initialDeformation[i];
+    for (int j = 0; j < dim; j ++) { // Displacement part
+      xRBM[i].r[j] = x[_0][i][j];
+    }
+    xRBM[i].q = x[_1][i]; // Rotation part
+  }
+
+  //////////////////////////////////////////////////////////////////////////////
+  //  Compute the energy, assemble the Gradient and Hessian using 
+  //  the GeodesicFEAssemblerWrapper and the MixedGFEAssembler and compare!
+  //////////////////////////////////////////////////////////////////////////////
+  CorrectionTypeWrapped gradient;
+  MatrixTypeWrapped hessianMatrix;
+  double energy = assembler.computeEnergy(xRBM);
+  assembler.assembleGradientAndHessian(xRBM, gradient, hessianMatrix, true);
+  double gradientTwoNorm = gradient.two_norm();
+  double gradientInfinityNorm = gradient.infinity_norm();
+  double matrixFrobeniusNorm = hessianMatrix.frobenius_norm();
+
+  CorrectionType gradientMixed;
+  gradientMixed[_0].resize(x[_0].size());
+  gradientMixed[_1].resize(x[_1].size());
+  MatrixType hessianMatrixMixed;
+  double energyMixed = mixedAssembler.computeEnergy(x[_0], x[_1]);
+  mixedAssembler.assembleGradientAndHessian(x[_0], x[_1], gradientMixed[_0], gradientMixed[_1], hessianMatrixMixed, true);
+  double gradientMixedTwoNorm = gradientMixed.two_norm();
+  double gradientMixedInfinityNorm = gradientMixed.infinity_norm();
+  double matrixMixedFrobeniusNorm = hessianMatrixMixed.frobenius_norm();
+
+  if (std::abs(energy - energyMixed)/energyMixed > 1e-8)
+  {
+      std::cerr << std::setprecision(9);
+      std::cerr << "The energy calculated by the GeodesicFEAssemblerWrapper is " << energy << " but "
+                << energyMixed << " (calculated by the MixedGFEAssembler) was expected!" << std::endl;
+      return 1;
+  }
+  if ( std::abs(gradientTwoNorm - gradientMixedTwoNorm)/gradientMixedTwoNorm > 1e-8 ||
+       std::abs(gradientInfinityNorm - gradientMixedInfinityNorm)/gradientMixedInfinityNorm > 1e-8)
+  {
+      std::cerr << std::setprecision(9);
+      std::cerr << "The gradient infinity norm calculated by the GeodesicFEAssemblerWrapper is " << gradientInfinityNorm << " but "
+                << gradientMixedInfinityNorm << " (calculated by the MixedGFEAssembler) was expected!" << std::endl;
+      std::cerr << "The gradient norm calculated by the GeodesicFEAssemblerWrapper is " << gradientTwoNorm << " but "
+                << gradientMixedTwoNorm << " (calculated by the MixedGFEAssembler) was expected!" << std::endl;
+      return 1;
+  }
+
+  if (std::abs(matrixFrobeniusNorm - matrixMixedFrobeniusNorm)/matrixMixedFrobeniusNorm > 1e-8)
+  {
+      std::cerr << std::setprecision(9);
+      std::cerr << "The matrix norm calculated by the GeodesicFEAssemblerWrapper is " << matrixFrobeniusNorm << " but "
+                << matrixMixedFrobeniusNorm << " (calculated by the MixedGFEAssembler) was expected!" << std::endl;
+      return 1;
+  }
+}