Add Simo-Fox shell model

This commit adds the energy implementat for the nonlinear Reissner-Mindlin
shell model known as the Simo and Fox shell model.

Details may be found in

  "A consistent finite element formulation of the geometrically
   non-linear Reissner-Mindlin shell model [Müller, Bischoff]"

dune/gfe/simofoxenergy.hh

+#ifndef DUNE_GFE_SIMOFOX_ENERGY_HH
+#define DUNE_GFE_SIMOFOX_ENERGY_HH
+
+#include <dune/common/fmatrix.hh>
+#include <dune/common/parametertree.hh>
+#if DUNE_VERSION_GTE(DUNE_COMMON, 2, 8)
+#include <dune/common/transpose.hh>
+#endif
+#include <dune/common/tuplevector.hh>
+#include <dune/fufem/boundarypatch.hh>
+#include <dune/geometry/quadraturerules.hh>
+#include <dune/gfe/linearalgebra.hh>
+#include <dune/gfe/localenergy.hh>
+#include <dune/gfe/mixedlocalgeodesicfestiffness.hh>
+#include <dune/gfe/realtuple.hh>
+#include <dune/gfe/unitvector.hh>
+
+namespace Dune::GFE {
+
+  /** \brief Get LocalFiniteElements from a localView, for different tree depths of the local view
+   *
+   * Copied from CosseratEnergyLocalStiffness.hh
+   */
+  template <class Basis, std::size_t i>
+  class LocalFiniteElementFactory {
+  public:
+    static auto get(const typename Basis::LocalView &localView, std::integral_constant<std::size_t, i> iType)
+        -> decltype(localView.tree().child(iType, 0).finiteElement()) {
+      return localView.tree().child(iType, 0).finiteElement();
+    }
+  };
+
+  /** \brief Specialize for scalar bases, here we cannot call tree().child() */
+  template <class GridView, int order, std::size_t i>
+  class LocalFiniteElementFactory<Dune::Functions::LagrangeBasis<GridView, order>, i> {
+  public:
+    static auto get(const typename Dune::Functions::LagrangeBasis<GridView, order>::LocalView &localView,
+                    std::integral_constant<std::size_t, i> iType) -> decltype(localView.tree().finiteElement()) {
+      return localView.tree().finiteElement();
+    }
+  };
+
+
+  /** \brief Implements the energy of a Simo-Fox shell
+   *
+   * Described in:
+   *  "A consistent finite element formulation of the geometrically non-linear Reissner-Mindlin shell model [Müller,
+   * Bischoff]"
+   *
+   * \tparam LocalFEFunction Decides which interpolation function is used for the interpolation of the midsurface
+   * and director field.
+   */
+  template <class Basis, template <int, typename, typename, typename> typename LocalFEFunction, typename field_type = double>
+  class SimoFoxEnergyLocalStiffness
+      : public Dune::GFE::LocalEnergy<Basis, RealTuple<field_type, 3>,
+                                      UnitVector<field_type, 3>>,  // inheritance to allow usage with LocalGeodesicFEADOLCStiffness
+        public MixedLocalGeodesicFEStiffness<Basis, RealTuple<field_type, 3>,
+                                             UnitVector<field_type, 3>>  // inheritance to allow usage with MixedGFEAssembler
+  {
+    // grid types
+    typedef typename Basis::GridView GridView;
+    typedef typename GridView::ctype DT;
+    typedef field_type RT;
+    typedef typename GridView::template Codim<0>::Entity Entity;
+
+    // some other sizes
+    static constexpr int gridDim  = GridView::dimension;
+    static constexpr int dimworld = GridView::dimensionworld;
+
+    // The local finite element type used for midsurface position and displacement interpolation
+    using MidSurfaceElement =
+        typename std::result_of<decltype (&LocalFiniteElementFactory<Basis, 0>::get)(typename Basis::LocalView, decltype(Dune::Indices::_0))>::type;
+    // The local finite element type used for the director interpolation
+    using DirectorElement =
+        typename std::result_of<decltype (&LocalFiniteElementFactory<Basis, 1>::get)(typename Basis::LocalView, decltype(Dune::Indices::_1))>::type;
+
+    // The local finite element function type to evaluate the midsurface position
+    using LocalMidSurfaceFunctionType = LocalFEFunction<gridDim, DT, MidSurfaceElement, RealTuple<field_type, 3>>;
+    // The local finite element function type to evaluate the director
+    using LocalDirectorFunctionType = LocalFEFunction<gridDim, DT, DirectorElement, UnitVector<field_type, 3>>;
+
+    // Extra function type for reference quantities since they are unconditionally doubles, i.e. no ADOL-C types
+    using LocalMidSurfaceReferenceFunctionType = LocalFEFunction<gridDim, DT, MidSurfaceElement, RealTuple<double, 3>>;
+    using LocalDirectorReferenceFunctionType   = LocalFEFunction<gridDim, DT, DirectorElement, UnitVector<double, 3>>;
+
+  public:
+    /** \brief Constructor with a set of material parameters
+     * \param parameters The material parameters
+     * \param x0 reference configuration
+     */
+    SimoFoxEnergyLocalStiffness(const Dune::ParameterTree &parameters, const BoundaryPatch<GridView> *neumannBoundary,
+                                const std::function<Dune::FieldVector<double, 3>(Dune::FieldVector<double, 2>)> neumannFunction,
+                                const std::function<Dune::FieldVector<double, 3>(Dune::FieldVector<double, 2>)> volumeLoad,
+                                const Dune::TupleVector<std::vector<RealTuple<double, 3>>, std::vector<UnitVector<double, 3>>> &x0)
+        : neumannBoundary_(neumannBoundary),
+          neumannFunction_(neumannFunction),
+          thickness_{parameters.template get<double>("thickness")},  // The sheeqll thickness
+          mu_{parameters.template get<double>("mu")},                // Lame constant 1
+          lambda_{parameters.template get<double>("lambda")},        // Lame constant 2
+          kappa_{parameters.template get<double>("kappa")},          // Shear correction factor
+          midSurfaceRefConfig{x0[Dune::Indices::_0]},
+          directorRefConfig{x0[Dune::Indices::_1]}
+    {
+      // Calculate the over the thickness preintegrated St.Venant Kirchhoff material matrix, Paper equation 10.1 */
+      const double Emodul = mu_ * (3 * lambda_ + 2 * mu_) / (lambda_ + mu_);  // Young's modulus
+      const double nu     = lambda_ / (2 * (lambda_ + mu_));                  // Poisson ratio
+
+      // membrane
+      const double fac1 = thickness_ * Emodul / (1 - nu * nu);
+      CMat_[0][0] = CMat_[1][1] = fac1;
+      CMat_[2][2]               = fac1 * (1 - nu) * 0.5;
+      CMat_[1][0] = CMat_[0][1] = fac1 * nu;
+
+      // bending
+      const double fac2 = thickness_ * thickness_ * thickness_ / 12 * Emodul / (1 - nu * nu);
+      CMat_[3][3] = CMat_[4][4] = fac2;
+      CMat_[5][5]               = fac2 * (1 - nu) * 0.5;
+      CMat_[3][4] = CMat_[4][3] = fac2 * nu;
+
+      // transverse shear
+      const double fac3 = kappa_ * thickness_ * Emodul * 0.5 / (1 + nu);
+      CMat_[6][6] = CMat_[7][7] = fac3;
+    }
+
+    /** \brief Assemble the energy for a single element */
+    RT energy(const typename Basis::LocalView &localView, const std::vector<RealTuple<field_type, 3>> &localMidSurfaceConfiguration,
+              const std::vector<UnitVector<field_type, 3>> &localDirectorConfiguration) const override;
+
+  private:
+    /** \brief A structure that contains all quantities to calculate the Lagrangian strains */
+    struct KinematicVariables {
+      // current configuration
+      FieldMatrix<field_type, 2, 3> a1anda2;    // first and second tangent base vector on the current geometry
+      FieldMatrix<field_type, 2, 3> td1Andtd2;  // partial derivative of the director on the current geometry in first and second direction
+      FieldMatrix<field_type, 2, 3> ud1andud2;  // partial derivative of the displacement function in first and second direction
+      FieldVector<field_type, 3> t;             // unit director on the current geometry
+
+      // reference configuration
+      FieldMatrix<double, 2, 3> A1andA2;      // first and second tangent base vector on the reference geometry
+      FieldMatrix<double, 2, 3> t0d1Andt0d2;  // partial derivative of the director on the reference geometry in first and second direction
+      FieldVector<double, 3> t0;              // unit director on the reference geometry
+    };
+
+    auto getReferenceLocalConfigurations(const typename Basis::LocalView &localView) const;
+
+    /** \brief Calculates all kinematic quantities */
+    template <typename Element, typename LocalDirectorFunction, typename LocalMidSurfaceFunction, typename LocalDirectorReferenceFunction,
+              typename LocalMidSurfaceReferenceFunction, typename IntegrationPointPosition>
+    static auto kinematicVariablesFactory(const Element &element, const LocalDirectorFunction &directorFunction,
+                                          const LocalDirectorReferenceFunction &directorReferenceFunction,
+                                          const LocalMidSurfaceFunction &midSurfaceFunction,
+                                          const LocalMidSurfaceReferenceFunction &midSurfaceReferenceFunction,
+                                          const LocalMidSurfaceFunction &midSurfaceDisplacementFunction, const IntegrationPointPosition &quadPos);
+
+    /** \brief Calculate the Green-Lagrange strain components */
+    static Dune::FieldVector<RT, 8> calculateGreenLagrangianStrains(const KinematicVariables &kin);
+
+    /** \brief Save all tangent base matrices for all nodes in one place*/
+    Dune::BlockVector<Dune::FieldMatrix<field_type, 2, 3>> directorTangentSpaces;
+
+    /** \brief The Neumann boundary */
+    const BoundaryPatch<GridView> *neumannBoundary_;
+
+    /** \brief The function implementing the Neumann data */
+    const std::function<Dune::FieldVector<double, 3>(Dune::FieldVector<double, dimworld>)> neumannFunction_;
+
+    /** \brief The shell thickness */
+    double thickness_;
+
+    /** \brief Lame constants */
+    double mu_, lambda_;
+
+    /** \brief Shear correction factor */
+    double kappa_;
+
+    /** \brief Material tangent matrix */
+    Dune::FieldMatrix<double, 8, 8> CMat_;
+
+    /** \brief The function implementing a volume load */
+    const std::function<Dune::FieldVector<double, 3>(Dune::FieldVector<double, dimworld>)> volumeLoad_;
+
+    /** \brief Stores the reference configuration of the midsurface and the director field */
+    const std::vector<RealTuple<double, 3>> &midSurfaceRefConfig;
+    const std::vector<UnitVector<double, 3>> &directorRefConfig;
+  };
+
+  /** \brief Calculate the Green-Lagrange strain components for the thickness-integrated setting
+   *
+   * The components are returned as [a_11, a_22, a_12, b_11, b_22, b_12, gamma_1, gamma_2]
+   * where a is the midsurface metric,
+   * b is similar to the second fundamental form of the surface, but the unit normal is replaced with the shell director
+   * gamma_1 and gamma_2 are the transverse shear in the two parametric directions
+   * Paper Equation 4.10 */
+  template <class Basis, template <int, typename, typename, typename> typename LocalFEFunction, typename field_type>
+  Dune::FieldVector<field_type, 8> SimoFoxEnergyLocalStiffness<Basis, LocalFEFunction, field_type>::calculateGreenLagrangianStrains(
+      const KinematicVariables &kin)
+  {
+    Dune::FieldVector<RT, 8> egl;
+    // membrane
+    egl[0] = kin.A1andA2[0] * kin.ud1andud2[0] + 0.5 * kin.ud1andud2[0].two_norm2();
+    egl[1] = kin.A1andA2[1] * kin.ud1andud2[1] + 0.5 * kin.ud1andud2[1].two_norm2();
+    egl[2] = kin.a1anda2[0] * kin.a1anda2[1];
+
+    // bending
+    egl[3] = kin.a1anda2[0] * kin.td1Andtd2[0] - kin.A1andA2[0] * kin.t0d1Andt0d2[0];
+    egl[4] = kin.a1anda2[1] * kin.td1Andtd2[1] - kin.A1andA2[1] * kin.t0d1Andt0d2[1];
+    egl[5] = kin.a1anda2[0] * kin.td1Andtd2[1] + kin.a1anda2[1] * kin.td1Andtd2[0] - kin.A1andA2[0] * kin.t0d1Andt0d2[1]
+             - kin.A1andA2[1] * kin.t0d1Andt0d2[0];
+
+    // transverse shear
+    egl[6] = kin.a1anda2[0] * kin.t - kin.A1andA2[0] * kin.t0;
+    egl[7] = kin.a1anda2[1] * kin.t - kin.A1andA2[1] * kin.t0;
+
+    return egl;
+  }
+
+  /** \brief Calculates all kinematic quantities that are needed for strain calculation
+   *
+   * Paper Equations 6.4 - 6.8
+   */
+  template <class Basis, template <int, typename, typename, typename> typename LocalFEFunction, typename field_type>
+  template <typename Element, typename LocalDirectorFunction, typename LocalMidSurfaceFunction, typename LocalDirectorReferenceFunction,
+            typename LocalMidSurfaceReferenceFunction, typename IntegrationPointPosition>
+  auto SimoFoxEnergyLocalStiffness<Basis, LocalFEFunction, field_type>::kinematicVariablesFactory(
+      const Element &element, const LocalDirectorFunction &directorFunction, const LocalDirectorReferenceFunction &directorReferenceFunction,
+      const LocalMidSurfaceFunction &midSurfaceFunction, const LocalMidSurfaceReferenceFunction &midSurfaceReferenceFunction,
+      const LocalMidSurfaceFunction &midSurfaceDisplacementFunction, const IntegrationPointPosition &quadPos)
+  {
+    KinematicVariables kin{};
+
+    const auto jInvT = element.geometry().jacobianInverseTransposed(quadPos);
+
+    kin.t           = directorFunction.evaluate(quadPos).globalCoordinates();
+    kin.t0          = directorReferenceFunction.evaluate(quadPos).globalCoordinates();
+    kin.a1anda2     = transpose(midSurfaceFunction.evaluateDerivative(quadPos) * transpose(jInvT));
+    kin.A1andA2     = transpose(midSurfaceReferenceFunction.evaluateDerivative(quadPos) * transpose(jInvT));
+    kin.ud1andud2   = transpose(midSurfaceDisplacementFunction.evaluateDerivative(quadPos) * transpose(jInvT));
+    kin.t0d1Andt0d2 = transpose(directorReferenceFunction.evaluateDerivative(quadPos) * transpose(jInvT));
+    kin.td1Andtd2   = transpose(directorFunction.evaluateDerivative(quadPos) * transpose(jInvT));
+
+    return kin;
+  }
+
+  template <class Basis, template <int, typename, typename, typename> typename LocalFEFunction, typename field_type>
+  auto SimoFoxEnergyLocalStiffness<Basis, LocalFEFunction, field_type>::getReferenceLocalConfigurations(
+      const typename Basis::LocalView &localView) const {
+    using namespace Dune::TypeTree::Indices;
+    const int nDofs0 = localView.tree().child(_0, 0).finiteElement().size();
+    const int nDofs1 = localView.tree().child(_1, 0).finiteElement().size();
+
+    std::vector<RealTuple<double, 3>> localConfiguration0(nDofs0);
+    std::vector<UnitVector<double, 3>> localConfiguration1(nDofs1);
+
+    for (int i = 0; i < nDofs0 + nDofs1; i++) {
+      int localIndexI = 0;
+      if (i < nDofs0) {
+        auto &node  = localView.tree().child(_0, 0);
+        localIndexI = node.localIndex(i);
+      } else {
+        auto &node  = localView.tree().child(_1, 0);
+        localIndexI = node.localIndex(i - nDofs0);
+      }
+      auto multiIndex = localView.index(localIndexI);
+
+      // The CompositeBasis number is contained in multiIndex[0]
+      // multiIndex[1] contains the actual index
+      if (multiIndex[0] == 0)
+        localConfiguration0[i] = midSurfaceRefConfig[multiIndex[1]];
+      else if (multiIndex[0] == 1)
+        localConfiguration1[i - nDofs0] = directorRefConfig[multiIndex[1]];
+    }
+    return std::make_tuple(localConfiguration0, localConfiguration1);
+  }
+
+  /** \brief Computes the Simo-Fox shell energy
+   *
+   *  The energy is 0.5 * S * E = 0.5 * transpose(E) * Cmat * E, where
+   *  S are the stress resultants [membrane forces, bending moments, transverse shear forces]
+   *  E are the components of the Green-Lagrangian strains [membrane strains, bending , tranverse shear]
+   *  see for details Paper Equation 4.11,4.10 and 10.1
+   */
+  template <class Basis, template <int, typename, typename, typename> typename LocalFEFunction, typename field_type>
+  typename SimoFoxEnergyLocalStiffness<Basis, LocalFEFunction, field_type>::RT
+  SimoFoxEnergyLocalStiffness<Basis, LocalFEFunction, field_type>::energy(
+      const typename Basis::LocalView &localView, const std::vector<RealTuple<field_type, 3>> &localMidSurfaceConfiguration,
+      const std::vector<UnitVector<field_type, 3>> &localDirectorConfiguration) const
+  {
+    auto element = localView.element();
+
+    using namespace Dune::Indices;
+    const auto &midSurfaceElement = LocalFiniteElementFactory<Basis, 0>::get(localView, _0);
+    const auto &directorElement   = LocalFiniteElementFactory<Basis, 1>::get(localView, _1);
+
+    const auto [localRefMidSurfaceConfiguration, localRefDirectorConfiguration] = getReferenceLocalConfigurations(localView);
+
+    std::vector<RealTuple<field_type, 3>> displacements;
+    displacements.reserve(localMidSurfaceConfiguration.size());
+    for (size_t i = 0; i < localMidSurfaceConfiguration.size(); ++i)
+      displacements.emplace_back(localMidSurfaceConfiguration[i].globalCoordinates() - localRefMidSurfaceConfiguration[i].globalCoordinates());
+
+    const LocalMidSurfaceFunctionType localMidSurfaceFunction(midSurfaceElement, localMidSurfaceConfiguration);
+    const LocalMidSurfaceFunctionType localMidSurfaceDisplacementFunction(midSurfaceElement, displacements);
+    const LocalMidSurfaceReferenceFunctionType localMidSurfaceReferenceFunction(midSurfaceElement, localRefMidSurfaceConfiguration);
+    const LocalDirectorFunctionType localDirectorFunction(directorElement, localDirectorConfiguration);
+    const LocalDirectorReferenceFunctionType localDirectorReferenceFunction(directorElement, localRefDirectorConfiguration);
+
+    const int quadOrder = (element.type().isSimplex()) ? std::max(midSurfaceElement.localBasis().order(), directorElement.localBasis().order())
+                                                       : std::max(midSurfaceElement.localBasis().order(), directorElement.localBasis().order()) + 1;
+
+    const auto &quad = Dune::QuadratureRules<DT, gridDim>::rule(element.type(), quadOrder);
+
+    RT energy = 0;
+    for (const auto &curQuad : quad) {
+      const Dune::FieldVector<DT, gridDim> &quadPos = curQuad.position();
+
+      const KinematicVariables kin
+          = kinematicVariablesFactory(element, localDirectorFunction, localDirectorReferenceFunction, localMidSurfaceFunction,
+                                      localMidSurfaceReferenceFunction, localMidSurfaceDisplacementFunction, quadPos);
+
+      const DT integrationElement = element.geometry().integrationElement(quadPos);
+      const FieldVector<field_type, 8> Egl = calculateGreenLagrangianStrains(kin);
+      energy +=  0.5 * Egl * (CMat_ * Egl) * curQuad.weight() * integrationElement;
+    }
+
+    //////////////////////////////////////////////////////////////////////////////
+    //   Assemble boundary contributions
+    //////////////////////////////////////////////////////////////////////////////
+
+    if (not neumannFunction_) return energy;
+
+    for (auto &&it : intersections(neumannBoundary_->gridView(), element))
+    {
+      if (not neumannBoundary_ or not neumannBoundary_->contains(it)) continue;
+
+      const auto &quadLine = QuadratureRules<DT, gridDim - 1>::rule(it.type(), quadOrder);
+
+      for (const auto &curQuad : quadLine)
+      {
+        // Local position of the quadrature point
+        const FieldVector<DT, gridDim> &quadPos = it.geometryInInside().global(curQuad.position());
+
+        const DT integrationElement = it.geometry().integrationElement(curQuad.position());
+
+        // The value of the local function
+        RealTuple<field_type, 3> deformationValue = localMidSurfaceFunction.evaluate(quadPos);
+
+        // Value of the Neumann data at the current position
+        auto neumannValue = neumannFunction_(it.geometry().global(curQuad.position()));
+
+        // Only translational dofs are affected by the Neumann force
+        for (size_t i = 0; i < neumannValue.size(); i++)
+          energy -= (neumannValue[i] * deformationValue.globalCoordinates()[i]) * curQuad.weight() * integrationElement;
+      }
+    }
+    return energy;
+  }
+}  // namespace Dune::GFE
+#endif  // DUNE_GFE_SIMOFOX_ENERGY_HH

problems/simofox-cantilever-dirichlet-values.py

+import math
+
+class DirichletValues:
+    def __init__(self, homotopyParameter):
+        self.homotopyParameter = homotopyParameter
+
+    def deformation(self, x):
+        # Dirichlet b.c. simply clamp the shell in the reference configuration
+        out = [x[0], x[1], 0]
+
+        return out
+
+    def director(self, x):
+        dir = [0, 0, 1]
+        return dir
+

problems/simofox-cantilever.parset

+#############################################
+#  Grid parameters
+#############################################
+
+structuredGrid = true
+
+# bounding box
+lower = 0 0
+upper = 10 5
+
+elements = 10 1
+
+# Number of grid levels
+numLevels = 1
+
+#############################################
+#  Solver parameters
+#############################################
+
+# Number of homotopy steps for the Dirichlet boundary conditions
+numHomotopySteps = 10
+
+# Tolerance of the trust region solver
+tolerance = 1e-6
+
+# Max number of steps of the trust region solver
+maxTrustRegionSteps = 200
+
+trustRegionScaling = 1 1 1 0.01 0.01
+
+# Initial trust-region radius
+initialTrustRegionRadius = 10
+
+# Number of multigrid iterations per trust-region step
+numIt = 400
+
+# Number of presmoothing steps
+nu1 = 3
+
+# Number of postsmoothing steps
+nu2 = 3
+
+# Number of coarse grid corrections
+mu = 1
+
+# Number of base solver iterations
+baseIt = 1
+
+# Tolerance of the multigrid solver
+mgTolerance = 1e-6
+
+# Tolerance of the base grid solver
+baseTolerance = 1e-6
+
+# Measure convergence
+instrumented = 0
+
+############################
+#   Material parameters
+############################
+
+
+## For the Wriggers L-shape example
+[materialParameters]
+
+# shell thickness
+thickness = 1
+
+# Lame parameters
+mu = 100000
+lambda = 0
+
+# Shear correction factor
+kappa = 1
+
+[]
+
+#############################################
+#  Boundary values
+#############################################
+
+problem = simofox-cantilever
+
+###  Python predicate specifying all Dirichlet grid vertices
+# x is the vertex coordinate
+dirichletVerticesPredicate = "x[0] < 0.01"
+
+###  Python predicate specifying all Neumann grid vertices
+# x is the vertex coordinate
+neumannVerticesPredicate = "x[0] > 9.99"
+
+###  Neumann values
+neumannValues = 0 0 320
+
+# Initial deformation
+initialDeformation = "[x[0], x[1], 0]"
+
+#startFromFile = yes
+#initialIterateFilename = initial_iterate.vtu

src/CMakeLists.txt

 set(programs compute-disc-error
              cosserat-continuum
              gradient-flow
-             harmonicmaps)
+             harmonicmaps
+             simofoxshell)
 #            rodobstacle)
 if (${DUNE_COMMON_VERSION} VERSION_GREATER_EQUAL 2.8)
   set (programs rod3d ${programs})

src/simofoxshell.cc

+#include <config.h>
+
+#include <array>
+#include <dune/fufem/utilities/adolcnamespaceinjections.hh>
+
+#include <dune/common/typetraits.hh>
+#include <dune/common/bitsetvector.hh>
+#include <dune/common/parametertree.hh>
+#include <dune/common/parametertreeparser.hh>
+#include <dune/common/tuplevector.hh>
+
+#include <dune/grid/utility/structuredgridfactory.hh>
+#include <dune/grid/io/file/gmshreader.hh>
+#include <dune/grid/uggrid.hh>
+
+#include <dune/functions/functionspacebases/lagrangebasis.hh>
+#include <dune/functions/functionspacebases/compositebasis.hh>
+#include <dune/functions/functionspacebases/powerbasis.hh>
+
+#include <dune/fufem/functiontools/boundarydofs.hh>
+#include <dune/fufem/dunepython.hh>
+
+#include <dune/gfe/mixedlocalgfeadolcstiffness.hh>
+#include <dune/gfe/simofoxenergy.hh>
+#include <dune/gfe/cosseratvtkwriter.hh>
+#include <dune/gfe/embeddedglobalgfefunction.hh>
+#include <dune/gfe/mixedgfeassembler.hh>
+#include <dune/gfe/mixedriemanniantrsolver.hh>
+#include <dune/gfe/unitvector.hh>
+#include <dune/gfe/localgeodesicfefunction.hh>
+#include <dune/gfe/localprojectedfefunction.hh>
+
+#if HAVE_DUNE_VTK
+#include <dune/vtk/vtkreader.hh>
+#endif
+
+template <int dim, class ctype, class LocalFiniteElement, class TS>
+using LocalFEFunction = Dune::GFE::LocalProjectedFEFunction<dim,ctype,LocalFiniteElement,TS>;
+//using LocalFEFunction = LocalGeodesicFEFunction<dim,ctype,LocalFiniteElement,TS>;
+
+// Order of the approximation space for the midsurface position
+const int midsurfaceOrder = 1;
+
+// Order of the approximation space for the director
+const int directorOrder = 1;
+
+using namespace Dune;
+
+int main(int argc, char *argv[]) try
+{
+  // Initialize MPI, finalize is done automatically on exit
+  Dune::MPIHelper &mpiHelper = MPIHelper::instance(argc, argv);
+
+  // Start Python interpreter
+  Python::start();
+  Python::Reference main = Python::import("__main__");
+  Python::run("import math");
+
+  Python::runStream()
+    << std::endl << "import sys"
+    << std::endl << "import os"
+    << std::endl << "sys.path.append(os.getcwd() + '/../../problems/')"
+    << std::endl;
+
+  using namespace TypeTree::Indices;
+  using SolutionType = TupleVector<std::vector<RealTuple<double,3> >, std::vector<UnitVector<double,3> > >;
+
+  // parse data file
+  ParameterTree parameterSet;
+  if (argc < 2)
+    DUNE_THROW(Exception, "Usage: ./simofoxshell inputfile.dat");
+
+  ParameterTreeParser::readINITree(argv[1], parameterSet);
+
+  ParameterTreeParser::readOptions(argc, argv, parameterSet);
+
+  // read solver settings
+  const auto numLevels = parameterSet.get<int>("numLevels");
+  const auto totalLoadSteps = parameterSet.get<int>("numHomotopySteps");
+  const auto tolerance = parameterSet.get<double>("tolerance");
+  const auto maxTrustRegionSteps = parameterSet.get<int>("maxTrustRegionSteps");
+  const auto initialTrustRegionRadius = parameterSet.get<double>("initialTrustRegionRadius");
+  const auto multigridIterations = parameterSet.get<int>("numIt");
+  const auto nu1 = parameterSet.get<int>("nu1");
+  const auto nu2 = parameterSet.get<int>("nu2");
+  const auto mu = parameterSet.get<int>("mu");
+  const auto baseIterations = parameterSet.get<int>("baseIt");
+  const auto mgTolerance = parameterSet.get<double>("mgTolerance");
+  const auto baseTolerance = parameterSet.get<double>("baseTolerance");
+  const auto instrumented = parameterSet.get<bool>("instrumented");
+  std::string resultPath = parameterSet.get("resultPath", "");
+
+  /////////////////////////////////////////
+  //    Create the grid
+  /////////////////////////////////////////
+  using Grid = UGGrid<2>;
+
+  std::shared_ptr<Grid> grid;
+
+  FieldVector<double, 2> lower(0), upper(1);
+
+  std::string structuredGridType = parameterSet["structuredGrid"];
+  if (parameterSet.get<bool>("structuredGrid"))
+  {
+    lower = parameterSet.get<FieldVector<double, 2> >("lower");
+    upper = parameterSet.get<FieldVector<double, 2> >("upper");
+
+    auto elements = parameterSet.get<std::array<unsigned int, 2> >("elements");
+    grid = StructuredGridFactory<Grid>::createCubeGrid(lower, upper, elements);
+
+  } else {
+    auto path = parameterSet.get<std::string>("path");
+    auto gridFile = parameterSet.get<std::string>("gridFile");
+
+    // Guess the grid file format by looking at the file name suffix
+    auto dotPos = gridFile.rfind('.');
+    if (dotPos == std::string::npos)
+      DUNE_THROW(IOError, "Could not determine grid input file format");
+    std::string suffix = gridFile.substr(dotPos, gridFile.length() - dotPos);
+
+    if (suffix == ".msh")
+      grid = std::shared_ptr<Grid>(GmshReader<Grid>::read(path + "/" + gridFile));
+    else if (suffix == ".vtu" or suffix == ".vtp")
+#if HAVE_DUNE_VTK
+      grid = VtkReader<Grid>::createGridFromFile(path + "/" + gridFile);
+#else
+      DUNE_THROW(NotImplemented, "Please install dune-vtk for VTK reading support!");
+#endif
+  }
+
+  grid->globalRefine(numLevels - 1);
+  grid->loadBalance();
+
+  if (mpiHelper.rank() == 0)
+    std::cout << "There are " << grid->leafGridView().comm().size() << " processes" << std::endl;
+
+  using GridView = Grid::LeafGridView;
+  GridView gridView = grid->leafGridView();
+
+  using namespace Dune::Functions::BasisFactory;
+
+  auto compositeBasis = makeBasis(gridView,
+                                  composite(
+                                    power<2>(lagrange<midsurfaceOrder>()),
+                                    power<2>(lagrange<directorOrder>())));
+
+  typedef Dune::Functions::LagrangeBasis<GridView, midsurfaceOrder> MidsurfaceFEBasis;
+  typedef Dune::Functions::LagrangeBasis<GridView, directorOrder> DirectorFEBasis;
+
+  MidsurfaceFEBasis midsurfaceFEBasis(gridView);
+  DirectorFEBasis directorFEBasis(gridView);
+
+  ///////////////////////////////////////////
+  //   Read Dirichlet values
+  ///////////////////////////////////////////
+
+  BitSetVector<1> dirichletVertices(gridView.size(2), false);
+  BitSetVector<1> neumannVertices(gridView.size(2), false);
+
+  const GridView::IndexSet &indexSet = gridView.indexSet();
+
+  // Make Python function that computes which vertices are on the Dirichlet boundary,
+  // based on the vertex positions.
+  std::string lambda = std::string("lambda x: (") + parameterSet.get<std::string>("dirichletVerticesPredicate") + std::string(")");
+  auto pythonDirichletVertices = Python::make_function<bool>(Python::evaluate(lambda));
+
+  // Same for the Neumann boundary
+  lambda = std::string("lambda x: (") + parameterSet.get<std::string>("neumannVerticesPredicate", "0") + std::string(")");
+  auto pythonNeumannVertices = Python::make_function<bool>(Python::evaluate(lambda));
+
+  for (auto &&vertex: vertices(gridView))
+  {
+    bool isDirichlet = pythonDirichletVertices(vertex.geometry().corner(0));
+    dirichletVertices[indexSet.index(vertex)] = isDirichlet;
+
+    bool isNeumann = pythonNeumannVertices(vertex.geometry().corner(0));
+    neumannVertices[indexSet.index(vertex)] = isNeumann;
+  }
+
+  BoundaryPatch<GridView> dirichletBoundary(gridView, dirichletVertices);
+  BoundaryPatch<GridView> neumannBoundary(gridView, neumannVertices);
+
+  if (mpiHelper.rank() == 0)
+    std::cout << "Neumann boundary has " << neumannBoundary.numFaces() << " faces\n";
+
+  BitSetVector<1> deformationDirichletNodes(midsurfaceFEBasis.size(), false);
+  constructBoundaryDofs(dirichletBoundary, midsurfaceFEBasis, deformationDirichletNodes);
+
+  BitSetVector<1> neumannNodes(midsurfaceFEBasis.size(), false);
+  constructBoundaryDofs(neumannBoundary, directorFEBasis, neumannNodes);
+
+  BitSetVector<3> deformationDirichletDofs(midsurfaceFEBasis.size(), false);
+  for (size_t i = 0; i < midsurfaceFEBasis.size(); i++)
+    if (deformationDirichletNodes[i][0])
+      for (int j = 0; j < 3; j++)
+        deformationDirichletDofs[i][j] = true;
+
+  BitSetVector<1> orientationDirichletNodes(directorFEBasis.size(), false);
+  constructBoundaryDofs(dirichletBoundary, directorFEBasis, orientationDirichletNodes);
+
+  BitSetVector<2> orientationDirichletDofs(directorFEBasis.size(), false);
+  for (size_t i = 0; i < directorFEBasis.size(); i++)
+    if (orientationDirichletNodes[i][0])
+      for (int j = 0; j < 2; j++)
+        orientationDirichletDofs[i][j] = true;
+
+  ////////////////////////////
+  //   Initial iterate
+  ////////////////////////////
+
+  SolutionType x;
+
+  x[_0].resize(midsurfaceFEBasis.size());
+
+  lambda = std::string("lambda x: (") + parameterSet.get<std::string>("initialDeformation") + std::string(")");
+  auto pythonInitialDeformation = Python::make_function<FieldVector<double, 3> >(Python::evaluate(lambda));
+
+  auto deformationPowerBasis = makeBasis(gridView,power<3>(lagrange<midsurfaceOrder>()));
+
+  std::vector<FieldVector<double, 3> > v;
+  Functions::interpolate(deformationPowerBasis, v, pythonInitialDeformation);
+  std::copy(v.begin(), v.end(), x[_0].begin());
+
+  x[_1].resize(directorFEBasis.size());
+
+  // The code currently assumes that the grid resides in the x-y plane
+  // Therefore, all references directors point into z direction
+  // If the reference is not planar one has to calculate the reference nodal directors.
+  // E.g. for bilinear grid elements this can be just the average of the 4 element normals
+  const FieldVector<double, 3> referenceDirectorVector({0, 0, 1});   //only plain reference!
+  std::fill(x[_1].begin(), x[_1].end(), referenceDirectorVector);
+
+  const SolutionType x0 = x;
+
+  ////////////////////////////////////////////////////////
+  //   Main homotopy loop
+  ////////////////////////////////////////////////////////
+
+  for (int loadStep = 0; loadStep < totalLoadSteps; loadStep++)
+  {
+    const double loadFactor = (loadStep + 1) * (1.0 / totalLoadSteps);
+    if (mpiHelper.rank() == 0)
+      std::cout << "Homotopy step: " << loadStep << ",    parameter: " << loadFactor << std::endl;
+
+    //////////////////////////////////////////////////////////////
+    //   Create an assembler for the energy functional
+    //////////////////////////////////////////////////////////////
+
+    const ParameterTree &materialParameters = parameterSet.sub("materialParameters");
+    FieldVector<double, 3> neumannValues{0, 0, 0};
+    if (parameterSet.hasKey("neumannValues"))
+      neumannValues = parameterSet.get<FieldVector<double, 3> >("neumannValues");
+
+    auto neumannFunction = [&](FieldVector<double, 2>) {
+                             auto nV = neumannValues;
+                             nV *= loadFactor;
+                             return nV;
+                           };
+
+    // Output initial iterate (of homotopy loop)
+    auto directorPowerBasis = makeBasis(gridView, power<3>(lagrange<directorOrder>()));
+
+    BlockVector<FieldVector<double, 3> > displacementInitial(compositeBasis.size({0}));
+    //calculates the global coordinates of midsurface displacements into displacementInitial
+    std::transform(x[_0].begin(),x[_0].end(),x0[_0].begin(),displacementInitial.begin(),[](const auto& x_0Entry,const auto& x0_0Entry){
+      return x_0Entry.globalCoordinates()-x0_0Entry.globalCoordinates();
+    });
+
+    //copies the global coordinates of the directors to directorInitial for a VTKWriter usable format
+    BlockVector<FieldVector<double, 3> > directorInitial(compositeBasis.size({1}));
+    std::transform(x[_1].begin(),x[_1].end(),directorInitial.begin(),[](const auto& x_1Entry){
+      return x_1Entry.globalCoordinates();
+    });
+
+    auto displacementFunctionInitial = Dune::Functions::makeDiscreteGlobalBasisFunction<FieldVector<double, 3> >(deformationPowerBasis,
+                                                                                                                displacementInitial);
+    auto directorFunctionInitial = Dune::Functions::makeDiscreteGlobalBasisFunction<FieldVector<double, 3> >(directorPowerBasis,
+                                                                                                            directorInitial);
+    //  We need to subsample, because VTK cannot natively display real second-order functions
+    SubsamplingVTKWriter<GridView> vtkWriter(gridView, Dune::refinementLevels(midsurfaceOrder - 1));
+    vtkWriter.addVertexData(displacementFunctionInitial, VTK::FieldInfo("displacement", VTK::FieldInfo::Type::scalar, 3));
+    vtkWriter.addVertexData(directorFunctionInitial, VTK::FieldInfo("director", VTK::FieldInfo::Type::scalar, 3));
+    vtkWriter.write(resultPath + "simo-fox_homotopy_" + "_" + std::to_string(neumannValues[2]) + "_" + std::to_string(0));
+
+    if (mpiHelper.rank() == 0) {
+      std::cout << "Material parameters:" << std::endl;
+      materialParameters.report();
+    }
+
+    // Assembler using ADOL-C
+    Dune::GFE::SimoFoxEnergyLocalStiffness<decltype(compositeBasis), LocalFEFunction,adouble> simoFoxEnergyADOLCLocalStiffness(materialParameters,
+                                                                                                    &neumannBoundary,
+                                                                                                    neumannFunction,
+                                                                                                    nullptr, x0);
+
+    MixedLocalGFEADOLCStiffness<decltype(compositeBasis),
+        RealTuple<double,3>,
+        UnitVector<double,3> > localGFEADOLCStiffness(&simoFoxEnergyADOLCLocalStiffness);
+
+    MixedGFEAssembler<decltype(compositeBasis),
+        RealTuple<double,3>, UnitVector<double,3> > assembler(compositeBasis, &localGFEADOLCStiffness);
+
+    // /////////////////////////////////////////////////
+    //   Create a Riemannian trust-region solver
+    // /////////////////////////////////////////////////
+
+    MixedRiemannianTrustRegionSolver<Grid,
+        decltype(compositeBasis),
+        MidsurfaceFEBasis, RealTuple<double,3>,
+        DirectorFEBasis, UnitVector<double,3> > solver;
+
+    solver.setup(*grid,
+                 &assembler,
+                 midsurfaceFEBasis,
+                 directorFEBasis,
+                 x,
+                 deformationDirichletDofs,
+                 orientationDirichletDofs,
+                 tolerance,
+                 maxTrustRegionSteps,
+                 initialTrustRegionRadius,
+                 multigridIterations,
+                 mgTolerance,
+                 mu, nu1, nu2,
+                 baseIterations,
+                 baseTolerance,
+                 instrumented);
+
+    solver.setScaling(parameterSet.get<FieldVector<double, 5> >("trustRegionScaling"));
+
+    ////////////////////////////////////////////////////////
+    //   Set Dirichlet values
+    ////////////////////////////////////////////////////////
+
+    Python::Reference dirichletValuesClass = Python::import(parameterSet.get<std::string>("problem") + "-dirichlet-values");
+
+    Python::Callable C = dirichletValuesClass.get("DirichletValues");
+
+    // Call a constructor.
+    Python::Reference dirichletValuesPythonObject = C(loadFactor);
+
+    // Extract object member functions as Dune functions
+    auto deformationDirichletValues = Python::make_function<FieldVector<double, 3> >(dirichletValuesPythonObject.get("deformation"));
+
+    std::vector<FieldVector<double, 3> > ddV;
+    Functions::interpolate(deformationPowerBasis, ddV, deformationDirichletValues, deformationDirichletDofs);
+
+    for (size_t j = 0; j < x[_0].size(); j++)
+      if (deformationDirichletNodes[j][0])
+        x[_0][j] = ddV[j];
+
+    // /////////////////////////////////////////////////////
+    //   Solve!
+    // /////////////////////////////////////////////////////
+
+    solver.setInitialIterate(x);
+    solver.solve();
+
+    x = solver.getSol();
+
+    // Output result of each load step
+    std::stringstream iAsAscii;
+    iAsAscii << loadStep + 1;
+
+    //calculates the global coordinates of midsurface displacements into displacementInitial
+    BlockVector<FieldVector<double, 3> > displacement(compositeBasis.size({0}));
+    std::transform(x[_0].begin(),x[_0].end(),x0[_0].begin(),displacement.begin(),[](const auto& x_0Entry,const auto& x0_0Entry){
+      return x_0Entry.globalCoordinates()-x0_0Entry.globalCoordinates();
+    });
+
+    //copies the global coordinates of the directors to directorInitial for a VTKWriter usable format
+    BlockVector<FieldVector<double, 3> > director(compositeBasis.size({1}));
+    std::transform(x[_1].begin(),x[_1].end(),director.begin(),[](const auto& x_1Entry){
+      return x_1Entry.globalCoordinates();
+    });
+
+    auto displacementFunction = Dune::Functions::makeDiscreteGlobalBasisFunction<FieldVector<double, 3> >(deformationPowerBasis,
+                                                                                                          displacement);
+    auto directorFunction = Dune::Functions::makeDiscreteGlobalBasisFunction<FieldVector<double, 3> >(directorPowerBasis,
+                                                                                                      director);
+    //  We need to subsample, because VTK cannot natively display real second-order functions
+    //  SubsamplingVTKWriter<GridView> vtkWriter(gridView, Dune::refinementLevels(midsurfaceOrder - 1));
+    vtkWriter.addVertexData(displacementFunction, VTK::FieldInfo("displacement", VTK::FieldInfo::Type::scalar, 3));
+    vtkWriter.addVertexData(directorFunction, VTK::FieldInfo("director", VTK::FieldInfo::Type::scalar, 3));
+    vtkWriter.write(resultPath + "simo-fox_homotopy_" + "_" + std::to_string(neumannValues[2]) + "_" + std::to_string(loadStep + 1));
+  }
+}
+catch (Exception &e) {
+  std::cout << e.what() << std::endl;
+}