Merge branch 'feature/cosserat-continuum-nonplanar' into 'master'

Cosserat-Continuum-Nonplanar

See merge request !72

dune.module

@@ -8,4 +8,4 @@ Version: svn
 Maintainer: oliver.sander@tu-dresden.de
 #depending on
 Depends: dune-common(>=2.7) dune-grid(>=2.7) dune-uggrid dune-istl dune-localfunctions dune-geometry (>=2.7) dune-functions (>=2.7) dune-solvers dune-fufem dune-elasticity (>=2.7)
-Suggests: dune-foamgrid dune-parmg dune-vtk dune-curvedgeometry
+Suggests: dune-foamgrid dune-parmg dune-vtk dune-curvedgrid

dune/gfe/cosseratvtkwriter.hh

@@ -109,6 +109,23 @@ class CosseratVTKWriter
     }
 
 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<RigidBodyMotion<double,3>> xRBM(basis.size());
+      for (int i = 0; i < basis.size(); i++) {
+        for (int j = 0; j < 3; j ++) // Displacement part
+          xRBM[i].r[j] = configuration[_0][i][j];
+        xRBM[i].q = configuration[_1][i];    // Rotation part
+      }
+      write(basis,xRBM,filename);
+    }
 
     /** \brief Write a configuration given with respect to a scalar function space basis
      */

dune/gfe/geometries/moebiusstrip.hh

+#ifndef DUNE_GFE_MOEBIUSSTRIP_GRIDFUNCTION_HH
+#define DUNE_GFE_MOEBIUSSTRIP_GRIDFUNCTION_HH
+
+#include <cmath>
+
+#include <dune/common/fmatrix.hh>
+#include <dune/common/fvector.hh>
+#include <dune/curvedgrid/gridfunctions/analyticgridfunction.hh>
+
+namespace Dune {
+
+/// \brief Functor representing a Möbius strip in 3D, where the base circle is the circle in the x-y-plane with radius_ around 0,0,0
+template <class T = double>
+class MoebiusStripProjection
+{
+  static const int dim = 3;
+  using Domain = FieldVector<T,dim>;
+  using Jacobian = FieldMatrix<T,dim,dim>;
+
+  T radius_;
+
+public:
+  MoebiusStripProjection (T radius)
+    : radius_(radius)
+  {}
+
+  /// \brief Project the coordinate to the MS using closest point projection
+  Domain operator() (const Domain& x) const
+  {
+    double nrm = std::sqrt(x[0]*x[0] + x[1]*x[1]);
+    //center for the point x - it lies on the circle around (0,0,0) with radius radius_ in the x-y-plane
+    Domain center{x[0] * radius_/nrm, x[1] * radius_/nrm, 0}; 
+    double cosu = x[0]/nrm;
+    double sinu = x[1]/nrm;
+    double cosuhalf = std::sqrt((1+cosu)/2);
+    cosuhalf *= (sinu < 0) ? (-1) : 1 ; // u goes from 0 to 2*pi, so cosuhalf is negative if sinu < 0, we multiply that here
+    double sinuhalf = std::sqrt((1-cosu)/2); //u goes from 0 to 2*pi, so sinuhalf is always >= 0 
+
+    // now calculate line from center to the new point
+    // the direction is (cosuhalf*cosu,cosuhalf*sinu,sinuhalf), as can be seen from the formula to construct the MS, we normalize this vector
+    Domain centerToPointOnMS{cosuhalf*cosu,cosuhalf*sinu,sinuhalf};
+    centerToPointOnMS /= centerToPointOnMS.two_norm();
+
+    Domain centerToX = center - x;
+    // We need the length, let theta be the angle between the vector centerToX and center to centerToPointOnMS
+    // then cos(theta) = centerToX*centerToPointOnMS/(len(centerToX)*len(centerToPointOnMS))
+    // We want to project the point to MS, s.t. the angle between the MS and the projection is 90deg
+    // Then, length = cos(theta) * len(centerToX) = centerToX*centerToPointOnMS/len(centerToPointOnMS) = centerToX*centerToPointOnMS
+    double length = - centerToX * centerToPointOnMS;
+    centerToPointOnMS *= length;
+    return center + centerToPointOnMS;
+  }
+
+
+  /// \brief derivative of the projection to the MS
+  friend auto derivative (const MoebiusStripProjection& moebiusStrip)
+  {
+    DUNE_THROW(NotImplemented,"The derivative of the projection to the Möbius strip is not implemented yet!");
+    return [radius = moebiusStrip.radius_](const Domain& x)
+    {
+      Jacobian out;
+      return out;
+    };
+  }
+
+  /// \brief Normal Vector of the MS
+  Domain normal (const Domain& x) const
+  {
+    using std::sqrt;
+    Domain nVec = {x[0],x[1],0};
+    double nrm = std::sqrt(x[0]*x[0] + x[1]*x[1]);
+    double cosu = x[0]/nrm;
+    double sinu = x[1]/nrm;
+    double cosuhalf = std::sqrt((1+cosu)/2);
+    cosuhalf *= (sinu < 0) ? (-1) : 1 ; // u goes from 0 to 2*pi, so cosuhalf is negative if sinu < 0, we multiply that here
+    double sinuhalf = std::sqrt((1-cosu)/2);
+    nVec[2] = (-1)*cosuhalf*(cosu*x[0]+sinu*x[1])/sinuhalf;
+    nVec /= nVec.two_norm();
+    return nVec;
+  }
+
+  /// \brief The mean curvature of the MS
+  T mean_curvature (const Domain& /*x*/) const
+  {
+    DUNE_THROW(NotImplemented,"The mean curvature of the Möbius strip is not implemented yet!");
+    return 1;
+  }
+
+  /// \brief The area of the MS
+  T area () const
+  {
+    DUNE_THROW(NotImplemented,"The area of the Möbius strip is not implemented yet!");
+    return 1;
+  }
+};
+
+/// \brief construct a grid function representing the parametrization of a MS
+template <class Grid, class T>
+auto moebiusStripGridFunction (T radius)
+{
+  return analyticGridFunction<Grid>(MoebiusStripProjection<T>{radius});
+}
+
+} // end namespace Dune
+
+#endif // DUNE_GFE_MOEBIUSSTRIP_GRIDFUNCTION_HH

dune/gfe/geometries/twistedstrip.hh

+#ifndef DUNE_GFE_TWISTEDSTRIP_GRIDFUNCTION_HH
+#define DUNE_GFE_TWISTEDSTRIP_GRIDFUNCTION_HH
+
+#include <cmath>
+
+#include <dune/common/fmatrix.hh>
+#include <dune/common/fvector.hh>
+#include <dune/curvedgrid/gridfunctions/analyticgridfunction.hh>
+
+namespace Dune {
+
+/// \brief Functor representing a twisted strip in 3D, with length length_ and nTwists twists
+template <class T = double>
+class TwistedStripProjection
+{
+  static const int dim = 3;
+  using Domain = FieldVector<T,dim>;
+  using Jacobian = FieldMatrix<T,dim,dim>;
+
+  T length_;
+  double nTwists_;
+
+public:
+  TwistedStripProjection (T length, double nTwists)
+    : length_(length),  
+      nTwists_(nTwists)
+  {}
+
+  /// \brief Project the coordinate to the twisted strip using closest point projection
+  Domain operator() (const Domain& x) const
+  {
+    // Angle for the x - position of the point
+    double alpha = std::acos(-1)*nTwists_*x[0]/length_;
+    double cosalpha = std::cos(alpha);
+    double sinalpha = std::sin(alpha);
+
+    // Add the line from middle to the new point - the direction is (0,cosalpha,sinalpha)
+    // We need the distance, calculated by (y^2 + z^2)^0.5
+    double dist = std::sqrt(x[1]*x[1] + x[2]*x[2]);
+    double y = std::abs(cosalpha*dist);
+    if (x[1] < 0 and std::abs(x[1]) > std::abs(x[2])) { // only trust the sign of x[1] if it is larger than x[2]
+      y *= (-1);
+    } else if (std::abs(x[1]) <= std::abs(x[2])) {
+      if (cosalpha * sinalpha >= 0 and x[2] < 0) // if cosalpha*sinalpha > 0, x[1] and x[2] must have the same sign
+        y *= (-1);
+      if (cosalpha * sinalpha <= 0 and x[2] > 0) // if cosalpha*sinalpha < 0, x[1] and x[2] must have opposite signs
+        y *= (-1);
+    }
+    double z = std::abs(sinalpha*dist);
+    if (x[2] < 0 and std::abs(x[2]) > std::abs(x[1])) {
+      z *= (-1);
+    } else if (std::abs(x[2]) <= std::abs(x[1])) {
+      if (cosalpha * sinalpha >= 0 and x[1] < 0)
+        z *= (-1);
+      if (cosalpha * sinalpha <= 0 and x[1] > 0)
+        z *= (-1);
+    }
+    return {x[0], y , z};
+  }
+
+  /// \brief derivative of the projection to the twistedstrip
+  friend auto derivative (const TwistedStripProjection& twistedStrip)
+  {
+    DUNE_THROW(NotImplemented,"The derivative of the projection to the twisted strip is not implemented yet!");
+    return [length = twistedStrip.length_, nTwists = twistedStrip.nTwists_](const Domain& x)
+    {
+      Jacobian out;
+      return out;
+    };
+  }
+
+  /// \brief Normal Vector of the twisted strip
+  Domain normal (const Domain& x) const
+  {
+    // Angle for the x - position of the point
+    double alpha = std::acos(-1)*nTwists_*x[0]/length_;
+    std::cout << "normal at " << x[0] << "is " << -std::sin(alpha) << ", " << std::cos(alpha) << std::endl; 
+    return {0,-std::sin(alpha), std::cos(alpha)};
+  }
+
+  /// \brief The mean curvature of the twisted strip
+  T mean_curvature (const Domain& /*x*/) const
+  {
+    DUNE_THROW(NotImplemented,"The mean curvature of the twisted strip is not implemented yet!");
+    return 1;
+  }
+
+  /// \brief The area of the twisted strip
+  T area (T width) const
+  {
+    return length_*width;
+  }
+};
+
+/// \brief construct a grid function representing the parametrization of a twisted strip
+template <class Grid, class T>
+auto twistedStripGridFunction (T length, double nTwists)
+{
+  return analyticGridFunction<Grid>(TwistedStripProjection<T>{length, nTwists});
+}
+
+} // end namespace Dune
+
+#endif // DUNE_GFE_TWISTEDSTRIP_GRIDFUNCTION_HH

dune/gfe/nonplanarcosseratshellenergy.hh

@@ -9,6 +9,8 @@
 
 #include <dune/fufem/boundarypatch.hh>
 
+#include <dune/functions/gridfunctions/discreteglobalbasisfunction.hh>
+
 #include <dune/gfe/localenergy.hh>
 #include <dune/gfe/localgeodesicfefunction.hh>
 #include <dune/gfe/rigidbodymotion.hh>
@@ -21,7 +23,15 @@
 #include <dune/localfunctions/lagrange/lfecache.hh>
 #endif
 
-template<class Basis, int dim, class field_type=double>
+/** \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=double, class StressFreeStateGridFunction = 
+  Dune::Functions::DiscreteGlobalBasisFunction<Basis,std::vector<Dune::FieldVector<double, Basis::GridView::dimensionworld>> > >
 class NonplanarCosseratShellEnergy
   : public Dune::GFE::LocalEnergy<Basis,RigidBodyMotion<field_type,dim> >
 {
@@ -40,13 +50,16 @@ class NonplanarCosseratShellEnergy
 public:
 
   /** \brief Constructor with a set of material parameters
-   * \param parameters The material parameters
+   * \param parameters                  The material parameters
+   * \param stressFreeStateGridFunction Pointer to a parametrization 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)
-  : neumannBoundary_(neumannBoundary),
+  : stressFreeStateGridFunction_(stressFreeStateGridFunction),
+    neumannBoundary_(neumannBoundary),
     neumannFunction_(neumannFunction),
     volumeLoad_(volumeLoad)
   {
@@ -111,6 +124,9 @@ public:
   /** \brief Curvature parameters */
   double b1_, b2_, b3_;
 
+  /** \brief The geometry used for assembling */
+  const StressFreeStateGridFunction* stressFreeStateGridFunction_;
+
   /** \brief The Neumann boundary */
   const BoundaryPatch<GridView>* neumannBoundary_;
 
@@ -121,15 +137,34 @@ public:
   const std::function<Dune::FieldVector<double,3>(Dune::FieldVector<double,dimworld>)> volumeLoad_;
 };
 
-template <class Basis, int dim, class field_type>
-typename NonplanarCosseratShellEnergy<Basis,dim,field_type>::RT
-NonplanarCosseratShellEnergy<Basis,dim,field_type>::
+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 std::vector<RigidBodyMotion<field_type,dim> >& localSolution) const
 {
   // The element geometry
   auto element = localView.element();
+
+#if HAVE_DUNE_CURVEDGEOMETRY
+  // Construct a curved geometry of this element of the Cosserat shell in stress-free state
+  // When using element.geometry(), then the curvatures on the element are zero, when using a curved geometry, they are not
+  // If a parametrization representing the Cosserat shell in a stress-free state is given,
+  // this is used for the curved geometry approximation.
+  // 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) {
+      if (not stressFreeStateGridFunction_) {
+        return element.geometry().global(local);
+      }
+      auto localGridFunction = localFunction(*stressFreeStateGridFunction_);
+      localGridFunction.bind(element);
+      return localGridFunction(local);
+    }, 2); /*order*/
+#else
   auto geometry = element.geometry();
+#endif
 
   // The set of shape functions on this element
   const auto& localFiniteElement = localView.tree().finiteElement();
@@ -215,17 +250,8 @@ energy(const typename Basis::LocalView& localView,
         c += aScalar * eps[alpha][beta] * Dune::GFE::dyadicProduct(aContravariant[alpha], aContravariant[beta]);
 
 #if HAVE_DUNE_CURVEDGEOMETRY
-    // Construct a curved geometry to evaluate the derivative of the normal field on each quadrature point
-    // The variable local holds the local coordinates in the reference element
-    // and localGeometry.global maps them to the world coordinates
-    // we want to take the derivative of the normal field on the element in world coordinates
-    Dune::CurvedGeometry<DT, gridDim, dimworld, Dune::CurvedGeometryTraits<DT, Dune::LagrangeLFECache<DT,DT,gridDim>>> curvedGeometry(referenceElement(element),
-      [localGeometry=element.geometry()](const auto& local) {
-        return localGeometry.global(local);
-      }, 1); //order = 1
-
-    // Second fundamental form: The derivative of the normal field
-    auto normalDerivative = curvedGeometry.normalGradient(quad[pt].position());
+    // Second fundamental form: The derivative of the normal field, on each quadrature point
+    auto normalDerivative = geometry.normalGradient(quad[pt].position());
 #else
     //In case dune-curvedgeometry is not installed, the normal derivative is set to zero.
     Dune::FieldMatrix<double,3,3> normalDerivative(0);

dune/gfe/riemannianpnsolver.cc

@@ -67,7 +67,7 @@ setup(const GridType& grid,
     // //////////////////////////////////////////////////////////////////////////////////////
 
     typedef DuneFunctionsBasis<Basis> FufemBasis;
-    FufemBasis basis(assembler_->basis_);
+    FufemBasis basis(assembler_->getBasis());
     OperatorAssembler<FufemBasis,FufemBasis> operatorAssembler(basis, basis);
 
     LaplaceAssembler<GridType, typename FufemBasis::LocalFiniteElement, typename FufemBasis::LocalFiniteElement> laplaceStiffness;

problems/cosserat-continuum-cantilever.parset

@@ -24,7 +24,7 @@ numHomotopySteps = 1
 tolerance = 1e-3
 
 # Max number of steps of the trust region solver
-maxTrustRegionSteps = 200
+maxSolverSteps = 200
 
 trustRegionScaling = 1 1 1 0.01 0.01 0.01
 

problems/cosserat-continuum-wriggers-l-shape.parset

@@ -20,7 +20,7 @@ numHomotopySteps = 1
 tolerance = 1e-8
 
 # Max number of steps of the trust region solver
-maxTrustRegionSteps = 1000
+maxSolverSteps = 1000
 
 trustRegionScaling = 1 1 1 0.01 0.01 0.01
 

problems/film-on-substrate.parset

@@ -59,7 +59,7 @@ initialDeformation = "x"
 #############################################
 
 # Inner solver, cholmod or multigrid
-solvertype = multigrid
+solvertype = trustRegion
 
 # Number of homotopy steps for the Dirichlet boundary conditions
 numHomotopySteps = 1

src/cosserat-continuum.cc

+#define MIXED_SPACE 0
+
 #include <config.h>
 
 #include <fenv.h>
@@ -37,7 +39,6 @@
 #include <dune/solvers/solvers/iterativesolver.hh>
 #include <dune/solvers/norms/energynorm.hh>
 
-#include <dune/gfe/rigidbodymotion.hh>
 #include <dune/gfe/localgeodesicfeadolcstiffness.hh>
 #include <dune/gfe/mixedlocalgfeadolcstiffness.hh>
 #include <dune/gfe/cosseratenergystiffness.hh>
@@ -45,10 +46,17 @@
 #include <dune/gfe/cosseratvtkwriter.hh>
 #include <dune/gfe/cosseratvtkreader.hh>
 #include <dune/gfe/geodesicfeassembler.hh>
-#include <dune/gfe/riemanniantrsolver.hh>
 #include <dune/gfe/embeddedglobalgfefunction.hh>
 #include <dune/gfe/mixedgfeassembler.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
 
 #if HAVE_DUNE_VTK
 #include <dune/vtk/vtkreader.hh>
@@ -59,22 +67,17 @@
 const int dim = 2;
 const int dimworld = 2;
 
-//#define MIXED_SPACE
-
 // Order of the approximation space for the displacement
 const int displacementOrder = 2;
 
 // Order of the approximation space for the microrotations
 const int rotationOrder = 2;
 
-#ifndef MIXED_SPACE
+#if !MIXED_SPACE
 static_assert(displacementOrder==rotationOrder, "displacement and rotation order do not match!");
 
 // Image space of the geodesic fe functions
-typedef RigidBodyMotion<double,3> TargetSpace;
-
-// Tangent vector of the image space
-const int blocksize = TargetSpace::TangentVector::dimension;
+using TargetSpace = RigidBodyMotion<double, 3>;
 #endif
 
 using namespace Dune;
@@ -84,6 +87,15 @@ 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 << "DISPLACEMENTORDER = " << displacementOrder << ", ROTATIONORDER = " << rotationOrder << std::endl;
+        std::cout << "dim = " << dim << ", dimworld = " << dimworld << 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();
     Python::Reference main = Python::import("__main__");
@@ -97,12 +109,8 @@ int main (int argc, char *argv[]) try
         << std::endl;
 
     using namespace TypeTree::Indices;
-#ifdef MIXED_SPACE
     using SolutionType = TupleVector<std::vector<RealTuple<double,3> >,
                                      std::vector<Rotation<double,3> > >;
-#else
-    typedef std::vector<TargetSpace> SolutionType;
-#endif
 
     // parse data file
     ParameterTree parameterSet;
@@ -117,7 +125,8 @@ int main (int argc, char *argv[]) try
     const 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 initialRegularization    = parameterSet.get<double>("initialRegularization", 1000);
     const double initialTrustRegionRadius = parameterSet.get<double>("initialTrustRegionRadius");
     const int multigridIterations         = parameterSet.get<int>("numIt");
     const int nu1                         = parameterSet.get<int>("nu1");
@@ -145,6 +154,8 @@ int main (int argc, char *argv[]) try
 
     std::string structuredGridType = parameterSet["structuredGrid"];
     if (structuredGridType == "cube") {
+        if (dim!=dimworld)
+            DUNE_THROW(GridError, "Please use FoamGrid and read in a grid for problems with dim != dimworld.");
 
         lower = parameterSet.get<FieldVector<double,dimworld> >("lower");
         upper = parameterSet.get<FieldVector<double,dimworld> >("upper");
@@ -183,7 +194,7 @@ int main (int argc, char *argv[]) try
     GridView gridView = grid->leafGridView();
 
     using namespace Dune::Functions::BasisFactory;
-#ifdef MIXED_SPACE
+
     const int dimRotation = Rotation<double,dim>::embeddedDim;
     auto compositeBasis = makeBasis(
       gridView,
@@ -195,6 +206,13 @@ int main (int argc, char *argv[]) try
             lagrange<rotationOrder>()
         )
     ));
+    using CompositeBasis = decltype(compositeBasis);
+
+    auto deformationPowerBasis = makeBasis(
+        gridView,
+        power<3>(
+            lagrange<displacementOrder>()
+    ));
 
     typedef Dune::Functions::LagrangeBasis<GridView,displacementOrder> DeformationFEBasis;
     typedef Dune::Functions::LagrangeBasis<GridView,rotationOrder> OrientationFEBasis;
@@ -202,10 +220,6 @@ int main (int argc, char *argv[]) try
     DeformationFEBasis deformationFEBasis(gridView);
     OrientationFEBasis orientationFEBasis(gridView);
 
-#else
-    typedef Dune::Functions::LagrangeBasis<typename GridType::LeafGridView, displacementOrder> FEBasis;
-    FEBasis feBasis(gridView);
-#endif
 
     // /////////////////////////////////////////
     //   Read Dirichlet values
@@ -237,10 +251,10 @@ int main (int argc, char *argv[]) try
     BoundaryPatch<GridView> dirichletBoundary(gridView, dirichletVertices);
     BoundaryPatch<GridView> neumannBoundary(gridView, neumannVertices);
 
-    if (mpiHelper.rank()==0)
-      std::cout << "Neumann boundary has " << neumannBoundary.numFaces() << " faces\n";
-
-#ifdef MIXED_SPACE
+    std::cout << "On rank " << mpiHelper.rank() << ": Dirichlet boundary has " << dirichletBoundary.numFaces()
+              << " faces and " << dirichletVertices.count() << " nodes.\n";
+    std::cout << "On rank " << mpiHelper.rank() << ": Neumann boundary has " << neumannBoundary.numFaces() << " faces.\n";
+  
     BitSetVector<1> deformationDirichletNodes(deformationFEBasis.size(), false);
     constructBoundaryDofs(dirichletBoundary,deformationFEBasis,deformationDirichletNodes);
 
@@ -261,43 +275,26 @@ int main (int argc, char *argv[]) try
       if (orientationDirichletNodes[i][0])
         for (int j=0; j<3; j++)
           orientationDirichletDofs[i][j] = true;
-#else
-    BitSetVector<1> dirichletNodes(feBasis.size(), false);
-    constructBoundaryDofs(dirichletBoundary,feBasis,dirichletNodes);
-
-    BitSetVector<1> neumannNodes(feBasis.size(), false);
-    constructBoundaryDofs(neumannBoundary,feBasis,neumannNodes);
-
-
-    BitSetVector<blocksize> dirichletDofs(feBasis.size(), false);
-    for (size_t i=0; i<feBasis.size(); i++)
-      if (dirichletNodes[i][0])
-        for (int j=0; j<5; j++)
-          dirichletDofs[i][j] = true;
-#endif
 
     // //////////////////////////
     //   Initial iterate
     // //////////////////////////
 
-#ifdef MIXED_SPACE
     SolutionType x;
 
-    x[_0].resize(deformationFEBasis.size());
+    x[_0].resize(compositeBasis.size({0}));
 
     lambda = std::string("lambda x: (") + parameterSet.get<std::string>("initialDeformation") + std::string(")");
     auto pythonInitialDeformation = Python::make_function<FieldVector<double,3> >(Python::evaluate(lambda));
 
     std::vector<FieldVector<double,3> > v;
-    Functions::interpolate(deformationFEBasis, v, pythonInitialDeformation);
+    Functions::interpolate(deformationPowerBasis, v, pythonInitialDeformation);
 
     for (size_t i=0; i<x[_0].size(); i++)
       x[_0][i] = v[i];
 
-    x[_1].resize(orientationFEBasis.size());
-#else
-    SolutionType x(feBasis.size());
-
+    x[_1].resize(compositeBasis.size({1}));
+#if !MIXED_SPACE
     if (parameterSet.hasKey("startFromFile"))
     {
       std::shared_ptr<GridType> initialIterateGrid;
@@ -314,35 +311,32 @@ int main (int argc, char *argv[]) try
         initialIterateGrid = std::shared_ptr<GridType>(GmshReader<GridType>::read(path + "/" + initialIterateGridFilename));
       }
 
-      SolutionType initialIterate;
+      std::vector<TargetSpace> initialIterate;
       GFE::CosseratVTKReader::read(initialIterate, parameterSet.get<std::string>("initialIterateFilename"));
 
       typedef Dune::Functions::LagrangeBasis<typename GridType::LeafGridView, 2> InitialBasis;
       InitialBasis initialBasis(initialIterateGrid->leafGridView());
 
 #ifdef PROJECTED_INTERPOLATION
-      using LocalInterpolationRule  = LocalProjectedFEFunction<dim, double, FEBasis::LocalView::Tree::FiniteElement, TargetSpace>;
+      using LocalInterpolationRule  = LocalProjectedFEFunction<dim, double, DeformationFEBasis::LocalView::Tree::FiniteElement, TargetSpace>;
 #else
-      using LocalInterpolationRule  = LocalGeodesicFEFunction<dim, double, FEBasis::LocalView::Tree::FiniteElement, TargetSpace>;
+      using LocalInterpolationRule  = LocalGeodesicFEFunction<dim, double, DeformationFEBasis::LocalView::Tree::FiniteElement, TargetSpace>;
 #endif
       GFE::EmbeddedGlobalGFEFunction<InitialBasis,LocalInterpolationRule,TargetSpace> initialFunction(initialBasis,initialIterate);
-
+      auto powerBasis = makeBasis(
+          gridView,
+          power<7>(
+              lagrange<displacementOrder>()
+      ));
       std::vector<FieldVector<double,7> > v;
-      Dune::Functions::interpolate(feBasis,v,initialFunction);
-      DUNE_THROW(NotImplemented, "Replace scalar basis by power basis!");
-
-      for (size_t i=0; i<x.size(); i++)
-        x[i] = TargetSpace(v[i]);
-
-    } else {
-    lambda = std::string("lambda x: (") + parameterSet.get<std::string>("initialDeformation") + std::string(")");
-      auto pythonInitialDeformation = Python::make_function<FieldVector<double,3> >(Python::evaluate(lambda));
+      //TODO: Interpolate does not work with an GFE:EmbeddedGlobalGFEFunction
+      //Dune::Functions::interpolate(powerBasis,v,initialFunction);
 
-    std::vector<FieldVector<double,3> > v;
-      Functions::interpolate(feBasis, v, pythonInitialDeformation);
-
-    for (size_t i=0; i<x.size(); i++)
-      x[i].r = v[i];
+      for (size_t i=0; i<x.size(); i++) {
+        auto vTargetSpace = TargetSpace(v[i]);
+        x[_0][i] = vTargetSpace.r;
+        x[_1][i] = vTargetSpace.q;
+      }
     }
 #endif
 
@@ -351,14 +345,30 @@ int main (int argc, char *argv[]) try
     ////////////////////////////////////////////////////////
 
     // Output initial iterate (of homotopy loop)
-#ifdef MIXED_SPACE
+    BlockVector<FieldVector<double,3> > identity(compositeBasis.size({0}));
+    Dune::Functions::interpolate(deformationPowerBasis, identity, [&](FieldVector<double,dimworld> x){ return x;});
+    BlockVector<FieldVector<double,3> > displacement(compositeBasis.size({0}));
+
+    if (dim == dimworld) {
     CosseratVTKWriter<GridType>::writeMixed<DeformationFEBasis,OrientationFEBasis>(deformationFEBasis,x[_0],
                                                                                    orientationFEBasis,x[_1],
                                                                                    resultPath + "mixed-cosserat_homotopy_0");
+    } else {
+#if MIXED_SPACE
+        for (int i = 0; i < displacement.size(); i++) {
+            for (int j = 0; j < 3; j++)
+                displacement[i][j] = x[_0][i][j];
+            displacement[i] -= identity[i];
+        }
+        auto displacementFunction = Dune::Functions::makeDiscreteGlobalBasisFunction<FieldVector<double,3>>(deformationPowerBasis, displacement);
+    
+        SubsamplingVTKWriter<GridView> vtkWriter(gridView, Dune::refinementLevels(0));
+        vtkWriter.addVertexData(displacementFunction, VTK::FieldInfo("displacement", VTK::FieldInfo::Type::scalar, dimworld));
+        vtkWriter.write(resultPath + "mixed-cosserat_homotopy_0");
 #else
-    CosseratVTKWriter<GridType>::write<FEBasis>(feBasis,x, resultPath + "cosserat_homotopy_0");
-#endif
-
+    CosseratVTKWriter<GridType>::write<DeformationFEBasis>(deformationFEBasis, x, resultPath + "cosserat_homotopy_0");
+#endif   
+    }
     for (int i=0; i<numHomotopySteps; i++) {
 
         double homotopyParameter = (i+1)*(1.0/numHomotopySteps);
@@ -375,7 +385,7 @@ int main (int argc, char *argv[]) try
     if (parameterSet.hasKey("neumannValues"))
         neumannValues = parameterSet.get<FieldVector<double,3> >("neumannValues");
 
-    auto neumannFunction = [&]( FieldVector<double,dim> ) {
+    auto neumannFunction = [&]( FieldVector<double,dimworld> ) {
       auto nV = neumannValues;
       nV *= homotopyParameter;
       return nV;
@@ -385,7 +395,7 @@ int main (int argc, char *argv[]) try
     if (parameterSet.hasKey("volumeLoad"))
         volumeLoadValues = parameterSet.get<FieldVector<double,3> >("volumeLoad");
 
-    auto volumeLoad = [&]( FieldVector<double,dim>) {
+    auto volumeLoad = [&]( FieldVector<double,dimworld>) {
       auto vL = volumeLoadValues;
       vL *= homotopyParameter;
       return vL;
@@ -396,82 +406,6 @@ int main (int argc, char *argv[]) try
         materialParameters.report();
     }
 
-    // Assembler using ADOL-C
-#ifdef MIXED_SPACE
-    CosseratEnergyLocalStiffness<decltype(compositeBasis),
-                        3,adouble> cosseratEnergyADOLCLocalStiffness(materialParameters,
-                                                                     &neumannBoundary,
-                                                                     neumannFunction,
-                                                                     volumeLoad);
-
-    MixedLocalGFEADOLCStiffness<decltype(compositeBasis),
-                                RealTuple<double,3>,
-                                Rotation<double,3> > localGFEADOLCStiffness(&cosseratEnergyADOLCLocalStiffness);
-
-    MixedGFEAssembler<decltype(compositeBasis),
-                      RealTuple<double,3>,
-                      Rotation<double,3> > assembler(compositeBasis, &localGFEADOLCStiffness);
-#else
-    std::shared_ptr<GFE::LocalEnergy<FEBasis,RigidBodyMotion<adouble,3> > > localCosseratEnergy;
-
-    if (dim==dimworld)
-    {
-      localCosseratEnergy = std::make_shared<CosseratEnergyLocalStiffness<FEBasis,3,adouble> >(materialParameters,
-                                                                                               &neumannBoundary,
-                                                                                               neumannFunction,
-                                                                                               volumeLoad);
-    }
-    else
-    {
-      localCosseratEnergy = std::make_shared<NonplanarCosseratShellEnergy<FEBasis,3,adouble> >(materialParameters,
-                                                                                               &neumannBoundary,
-                                                                                               neumannFunction,
-                                                                                               volumeLoad);
-    }
-
-    LocalGeodesicFEADOLCStiffness<FEBasis,
-                                  TargetSpace> localGFEADOLCStiffness(localCosseratEnergy.get());
-
-    GeodesicFEAssembler<FEBasis,TargetSpace> assembler(gridView, localGFEADOLCStiffness);
-#endif
-
-    // /////////////////////////////////////////////////
-    //   Create a Riemannian trust-region solver
-    // /////////////////////////////////////////////////
-
-#ifdef MIXED_SPACE
-    MixedRiemannianTrustRegionSolver<GridType,
-                                     decltype(compositeBasis),
-                                     DeformationFEBasis, RealTuple<double,3>,
-                                     OrientationFEBasis, Rotation<double,3> > solver;
-#else
-    RiemannianTrustRegionSolver<FEBasis,TargetSpace> solver;
-#endif
-    solver.setup(*grid,
-                 &assembler,
-#ifdef MIXED_SPACE
-                 deformationFEBasis,
-                 orientationFEBasis,
-#endif
-                 x,
-#ifdef MIXED_SPACE
-                 deformationDirichletDofs,
-                 orientationDirichletDofs,
-#else
-                 dirichletDofs,
-#endif
-                 tolerance,
-                 maxTrustRegionSteps,
-                 initialTrustRegionRadius,
-                 multigridIterations,
-                 mgTolerance,
-                 mu, nu1, nu2,
-                 baseIterations,
-                 baseTolerance,
-                 instrumented);
-
-        solver.setScaling(parameterSet.get<FieldVector<double,6> >("trustRegionScaling"));
-
         ////////////////////////////////////////////////////////
         //   Set Dirichlet values
         ////////////////////////////////////////////////////////
@@ -484,85 +418,227 @@ int main (int argc, char *argv[]) try
         Python::Reference dirichletValuesPythonObject = C(homotopyParameter);
 
         // Extract object member functions as Dune functions
-        auto deformationDirichletValues = Python::make_function<FieldVector<double,3> >(dirichletValuesPythonObject.get("deformation"));
-        auto orientationDirichletValues = Python::make_function<FieldMatrix<double,3,3> >(dirichletValuesPythonObject.get("orientation"));
-
-        std::vector<FieldVector<double,3> > ddV;
-        std::vector<FieldMatrix<double,3,3> > dOV;
-
-#ifdef MIXED_SPACE
-        Functions::interpolate(deformationFEBasis, ddV, deformationDirichletValues, deformationDirichletDofs);
-        Functions::interpolate(orientationFEBasis, dOV, orientationDirichletValues, orientationDirichletDofs);
-
-        for (size_t j=0; j<x[_0].size(); j++)
-          if (deformationDirichletNodes[j][0])
-            x[_0][j] = ddV[j];
-
-        for (size_t j=0; j<x[_1].size(); j++)
-          if (orientationDirichletNodes[j][0])
-            x[_1][j].set(dOV[j]);
+        auto deformationDirichletValues = Python::make_function<FieldVector<double,3> >   (dirichletValuesPythonObject.get("deformation"));
+        auto orientationDirichletValues = Python::make_function<FieldMatrix<double,3,3> > (dirichletValuesPythonObject.get("orientation"));
+    
+        BlockVector<FieldVector<double,3> > ddV;
+        Dune::Functions::interpolate(deformationPowerBasis, ddV, deformationDirichletValues, deformationDirichletDofs);
+    
+        BlockVector<FieldMatrix<double,3,3> > dOV;
+        Dune::Functions::interpolate(orientationFEBasis, dOV, orientationDirichletValues, orientationDirichletDofs);
+    
+        for (int i = 0; i < compositeBasis.size({0}); i++) {
+          if (deformationDirichletDofs[i][0])
+            x[_0][i] = ddV[i];
+        }
+        for (int i = 0; i < compositeBasis.size({1}); i++)
+          if (orientationDirichletDofs[i][0])
+            x[_1][i].set(dOV[i]);
+
+        if (dim==dimworld) {
+            CosseratEnergyLocalStiffness<CompositeBasis, 3,adouble> localCosseratEnergy(materialParameters,
+                                                                                            &neumannBoundary,
+                                                                                            neumannFunction,
+                                                                                            volumeLoad);
+            MixedLocalGFEADOLCStiffness<CompositeBasis,
+                                RealTuple<double,3>,
+                                Rotation<double,3> > localGFEADOLCStiffness(&localCosseratEnergy);
+            MixedGFEAssembler<CompositeBasis,
+                      RealTuple<double,3>,
+                      Rotation<double,3> > mixedAssembler(compositeBasis, &localGFEADOLCStiffness);
+#if MIXED_SPACE
+            MixedRiemannianTrustRegionSolver<GridType,
+                                         CompositeBasis,
+                                         DeformationFEBasis, RealTuple<double,3>,
+                                         OrientationFEBasis, Rotation<double,3> > 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
-        Functions::interpolate(feBasis, ddV, deformationDirichletValues, dirichletDofs);
-        Functions::interpolate(feBasis, 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]);
-          }
+            //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
+            std::vector<TargetSpace> xTargetSpace(compositeBasis.size({0}));
+            BitSetVector<TargetSpace::TangentVector::dimension> dirichletDofsTargetSpace(compositeBasis.size({0}), false);
+            for (int i = 0; i < compositeBasis.size({0}); i++) {
+              for (int j = 0; j < 3; j ++) { // Displacement part
+                xTargetSpace[i].r[j] = x[_0][i][j];
+                dirichletDofsTargetSpace[i][j] = deformationDirichletDofs[i][j];
+              }
+              xTargetSpace[i].q = x[_1][i]; // Rotation part
+              for (int j = 3; j < TargetSpace::TangentVector::dimension; j ++)
+                dirichletDofsTargetSpace[i][j] = orientationDirichletDofs[i][j-3];
+            }
+
+            using GFEAssemblerWrapper = Dune::GFE::GeodesicFEAssemblerWrapper<CompositeBasis, DeformationFEBasis, TargetSpace, RealTuple<double, 3>, Rotation<double,3>>;
+            GFEAssemblerWrapper assembler(&mixedAssembler, deformationFEBasis);
+            RiemannianTrustRegionSolver<DeformationFEBasis, TargetSpace, GFEAssemblerWrapper> solver;
+            solver.setup(*grid,
+                     &assembler,
+                     xTargetSpace,
+                     dirichletDofsTargetSpace,
+                     tolerance,
+                     maxSolverSteps,
+                     initialTrustRegionRadius,
+                     multigridIterations,
+                     mgTolerance,
+                     mu, nu1, nu2,
+                     baseIterations,
+                     baseTolerance,
+                     instrumented);
+
+            solver.setScaling(parameterSet.get<FieldVector<double,6> >("trustRegionScaling"));
+            solver.setInitialIterate(xTargetSpace);
+            solver.solve();
+            xTargetSpace = solver.getSol();
+            for (int i = 0; i < xTargetSpace.size(); i++) {
+              x[_0][i] = xTargetSpace[i].r;
+              x[_1][i] = xTargetSpace[i].q;
+            }
 #endif
+        } else {
+#if MIXED_SPACE
+            DUNE_THROW(Exception, "Problems with dim != dimworld do not work for MIXED_SPACE = 1!");
+#else
+            std::shared_ptr<LocalGeodesicFEADOLCStiffness<DeformationFEBasis, TargetSpace>> localGFEADOLCStiffness;
+
+            NonplanarCosseratShellEnergy<DeformationFEBasis, 3, adouble> localCosseratEnergyPlanar(materialParameters,
+                                                                                                    nullptr,
+                                                                                                    &neumannBoundary,
+                                                                                                    neumannFunction,
+                                                                                                    volumeLoad);
+              localGFEADOLCStiffness = std::make_shared<LocalGeodesicFEADOLCStiffness<DeformationFEBasis, TargetSpace>>(&localCosseratEnergyPlanar);
+
+            GeodesicFEAssembler<DeformationFEBasis,TargetSpace> assembler(gridView, *localGFEADOLCStiffness);
+            //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
+            std::vector<TargetSpace> xTargetSpace(compositeBasis.size({0}));
+            BitSetVector<TargetSpace::TangentVector::dimension> dirichletDofsTargetSpace(compositeBasis.size({0}), false);
+            for (int i = 0; i < compositeBasis.size({0}); i++) {
+              for (int j = 0; j < 3; j ++) { // Displacement part
+                xTargetSpace[i].r[j] = x[_0][i][j];
+                dirichletDofsTargetSpace[i][j] = deformationDirichletDofs[i][j];
+              }
+              xTargetSpace[i].q = x[_1][i]; // Rotation part
+              for (int j = 3; j < TargetSpace::TangentVector::dimension; j ++)
+                dirichletDofsTargetSpace[i][j] = orientationDirichletDofs[i][j-3];
+            }
+            if (parameterSet.get<std::string>("solvertype", "trustRegion") == "trustRegion") {
+
+                RiemannianTrustRegionSolver<DeformationFEBasis, TargetSpace> solver;
+                solver.setup(*grid,
+                         &assembler,
+                         xTargetSpace,
+                         dirichletDofsTargetSpace,
+                         tolerance,
+                         maxSolverSteps,
+                         initialTrustRegionRadius,
+                         multigridIterations,
+                         mgTolerance,
+                         mu, nu1, nu2,
+                         baseIterations,
+                         baseTolerance,
+                         instrumented);
+
+                solver.setScaling(parameterSet.get<FieldVector<double,6> >("trustRegionScaling"));
+                solver.setInitialIterate(xTargetSpace);
+                solver.solve();
+                xTargetSpace = solver.getSol();
+            } else { //parameterSet.get<std::string>("solvertype") == "proximalNewton"
+#if DUNE_VERSION_LT(DUNE_COMMON, 2, 8)
+                DUNE_THROW(Exception, "Please install dune-solvers >= 2.8 to use the Proximal Newton Solver with Cholmod!");
+#else
+                RiemannianProximalNewtonSolver<DeformationFEBasis, TargetSpace> solver;
+                solver.setup(*grid,
+                             &assembler,
+                             xTargetSpace,
+                             dirichletDofsTargetSpace,
+                             tolerance,
+                             maxSolverSteps,
+                             initialRegularization,
+                             instrumented);
+                solver.setInitialIterate(xTargetSpace);
+                solver.solve();
+                xTargetSpace = solver.getSol();
+#endif
+            }
 
-        // /////////////////////////////////////////////////////
-        //   Solve!
-        // /////////////////////////////////////////////////////
-
-        solver.setInitialIterate(x);
-        solver.solve();
-
-        x = solver.getSol();
-
+            for (int i = 0; i < xTargetSpace.size(); i++) {
+              x[_0][i] = xTargetSpace[i].r;
+              x[_1][i] = xTargetSpace[i].q;
+            }
+#endif
+        }
         // Output result of each homotopy step
         std::stringstream iAsAscii;
         iAsAscii << i+1;
-#ifdef MIXED_SPACE
+        if (dim == dimworld) {
         CosseratVTKWriter<GridType>::writeMixed<DeformationFEBasis,OrientationFEBasis>(deformationFEBasis,x[_0],
                                                                                        orientationFEBasis,x[_1],
                                                                                        resultPath + "mixed-cosserat_homotopy_" + iAsAscii.str());
-#else
-        CosseratVTKWriter<GridType>::write<FEBasis>(feBasis,x, resultPath + "cosserat_homotopy_" + iAsAscii.str());
+        } else {
+#if MIXED_SPACE
+            for (int i = 0; i < displacement.size(); i++) {
+               for (int j = 0; j  < 3; j++) {
+                displacement[i][j] = x[_0][i][j];
+              }
+              displacement[i] -= identity[i];
+            }
+            auto displacementFunction = Dune::Functions::makeDiscreteGlobalBasisFunction<FieldVector<double,3>>(deformationPowerBasis, displacement);
+    
+            //  We need to subsample, because VTK cannot natively display real second-order functions
+            SubsamplingVTKWriter<GridView> vtkWriter(gridView, Dune::refinementLevels(displacementOrder-1));
+            vtkWriter.addVertexData(displacementFunction, VTK::FieldInfo("displacement", VTK::FieldInfo::Type::scalar, 3));
+            vtkWriter.write(resultPath + "cosserat_homotopy_" + std::to_string(i+1));
+#else 
+        CosseratVTKWriter<GridType>::write<DeformationFEBasis>(deformationFEBasis, x, resultPath + "cosserat_homotopy_" + std::to_string(i+1));
 #endif
-
+        }
     }
 
     // //////////////////////////////
     //   Output result
     // //////////////////////////////
 
-#ifndef MIXED_SPACE
+#if !MIXED_SPACE
     // Write the corresponding coefficient vector: verbatim in binary, to be completely lossless
     // This data may be used by other applications measuring the discretization error
-    BlockVector<TargetSpace::CoordinateType> xEmbedded(x.size());
-    for (size_t i=0; i<x.size(); i++)
-      xEmbedded[i] = x[i].globalCoordinates();
+    BlockVector<TargetSpace::CoordinateType> xEmbedded(compositeBasis.size({0}));
+    for (size_t i=0; i<compositeBasis.size({0}); i++)
+        for (size_t j=0; j<TargetSpace::TangentVector::dimension; j++)
+            xEmbedded[i][j] = j < 3 ? x[_0][i][j] : x[_1][i].globalCoordinates()[j-3];
 
     std::ofstream outFile("cosserat-continuum-result-" + std::to_string(numLevels) + ".data", std::ios_base::binary);
     MatrixVector::Generic::writeBinary(outFile, xEmbedded);
     outFile.close();
 #endif
 
+    if (parameterSet.get<bool>("computeAverageNeumannDeflection", false) and parameterSet.hasKey("neumannValues")) {
     // finally: compute the average deformation of the Neumann boundary
     // That is what we need for the locking tests
     FieldVector<double,3> averageDef(0);
-#ifdef MIXED_SPACE
     for (size_t i=0; i<x[_0].size(); i++)
         if (neumannNodes[i][0])
             averageDef += x[_0][i].globalCoordinates();
-#else
-    for (size_t i=0; i<x.size(); i++)
-        if (neumannNodes[i][0])
-            averageDef += x[i].r;
-#endif
     averageDef /= neumannNodes.count();
 
     if (mpiHelper.rank()==0 and parameterSet.hasKey("neumannValues"))
@@ -570,6 +646,7 @@ int main (int argc, char *argv[]) try
       std::cout << "Neumann values = " << parameterSet.get<FieldVector<double, 3> >("neumannValues") << "  "
                 << ",  average deflection: " << averageDef << std::endl;
     }
+    }
 
 } catch (Exception& e)
 {

src/film-on-substrate.cc

@@ -554,7 +554,7 @@ int main (int argc, char *argv[]) try
     //   Create the chosen solver and solve!
     ///////////////////////////////////////////////////
 
-    if (parameterSet.get<std::string>("solvertype") == "multigrid") {
+    if (parameterSet.get<std::string>("solvertype", "trustRegion") == "trustRegion") {
 #if MIXED_SPACE
       MixedRiemannianTrustRegionSolver<GridType, CompositeBasis, DeformationFEBasis, RealTuple<double,dim>, OrientationFEBasis, Rotation<double,dim>> solver;
       solver.setup(*grid,
@@ -603,7 +603,7 @@ int main (int argc, char *argv[]) try
         x[_1][i] = xRBM[i].q;
       }
 #endif
-    } else { //parameterSet.get<std::string>("solvertype") == "cholmod"
+    } else { //parameterSet.get<std::string>("solvertype") == "proximalNewton"
 
 #if MIXED_SPACE
     DUNE_THROW(Exception, "Error: There is no MixedRiemannianProximalNewtonSolver!");