diff --git a/CMakeLists.txt b/CMakeLists.txt
index 055570791442f4e9f74198d8d350c3aaf2ffb323..417cba3f8ebe40a760e67cf22e0d0169fd19c86d 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -14,6 +14,7 @@ dune_project()
 add_subdirectory("test")
 add_subdirectory("dune")
 add_subdirectory("m4")
+add_subdirectory("src")
 
 set(programs linelast viscoelast nonlinelast)
 foreach(_program ${programs})
diff --git a/dune/elasticity/assemblers/CMakeLists.txt b/dune/elasticity/assemblers/CMakeLists.txt
index 60b0eefc751ca3b848f82315360ae5e2129b14b3..c510e61e630ed74c1eedf022285e38e95f6c5ef4 100644
--- a/dune/elasticity/assemblers/CMakeLists.txt
+++ b/dune/elasticity/assemblers/CMakeLists.txt
@@ -1,4 +1,7 @@
 install(FILES
+    feassembler.hh
+    localadolcstiffness.hh
+    localfestiffness.hh
     neohookeassembler.cc
     neohookeassembler.hh
     neohookefunctionalassembler.hh
diff --git a/dune/elasticity/assemblers/feassembler.hh b/dune/elasticity/assemblers/feassembler.hh
new file mode 100644
index 0000000000000000000000000000000000000000..e10bde18314eca3ea1f421ee161fb36533c7463d
--- /dev/null
+++ b/dune/elasticity/assemblers/feassembler.hh
@@ -0,0 +1,190 @@
+#ifndef DUNE_ELASTICITY_ASSEMBLERS_FEASSEMBLER_HH
+#define DUNE_ELASTICITY_ASSEMBLERS_FEASSEMBLER_HH
+
+#include <dune/istl/bcrsmatrix.hh>
+#include <dune/common/fmatrix.hh>
+#include <dune/istl/matrixindexset.hh>
+#include <dune/istl/matrix.hh>
+
+#include "localfestiffness.hh"
+
+
+/** \brief A global FE assembler for problems involving functions that map into non-Euclidean spaces
+ */
+template <class Basis, class VectorType>
+class FEAssembler {
+
+    typedef typename Basis::GridView GridView;
+    typedef typename GridView::template Codim<0>::template Partition<Dune::Interior_Partition>::Iterator ElementIterator;
+
+    //! Dimension of the grid.
+    enum { gridDim = GridView::dimension };
+
+    //! Dimension of a tangent space
+    enum { blocksize = VectorType::value_type::dimension };
+
+    //!
+    typedef Dune::FieldMatrix<double, blocksize, blocksize> MatrixBlock;
+
+public:
+    const Basis basis_;
+
+protected:
+
+    LocalFEStiffness<GridView,
+                             typename Basis::LocalFiniteElement,
+                             VectorType>* localStiffness_;
+
+public:
+
+    /** \brief Constructor for a given grid */
+    FEAssembler(const Basis& basis,
+                LocalFEStiffness<GridView,typename Basis::LocalFiniteElement, VectorType>* localStiffness)
+        : basis_(basis),
+          localStiffness_(localStiffness)
+    {}
+
+    /** \brief Assemble the tangent stiffness matrix and the functional gradient together
+     *
+     * This is more efficient than computing them separately, because you need the gradient
+     * anyway to compute the Riemannian Hessian.
+     */
+    virtual void assembleGradientAndHessian(const VectorType& sol,
+                                            Dune::BlockVector<Dune::FieldVector<double, blocksize> >& gradient,
+                                            Dune::BCRSMatrix<MatrixBlock>& hessian,
+                                            bool computeOccupationPattern=true) const;
+
+    /** \brief Compute the energy of a deformation state */
+    virtual double computeEnergy(const VectorType& sol) const;
+
+    //protected:
+    void getNeighborsPerVertex(Dune::MatrixIndexSet& nb) const;
+
+}; // end class
+
+
+
+template <class Basis, class TargetSpace>
+void FEAssembler<Basis,TargetSpace>::
+getNeighborsPerVertex(Dune::MatrixIndexSet& nb) const
+{
+    int n = basis_.size();
+
+    nb.resize(n, n);
+
+    ElementIterator it    = basis_.getGridView().template begin<0,Dune::Interior_Partition>();
+    ElementIterator endit = basis_.getGridView().template end<0,Dune::Interior_Partition>  ();
+
+    for (; it!=endit; ++it) {
+
+        const typename Basis::LocalFiniteElement& lfe = basis_.getLocalFiniteElement(*it);
+
+        for (size_t i=0; i<lfe.localBasis().size(); i++) {
+
+            for (size_t j=0; j<lfe.localBasis().size(); j++) {
+
+                int iIdx = basis_.index(*it,i);
+                int jIdx = basis_.index(*it,j);
+
+                nb.add(iIdx, jIdx);
+
+            }
+
+        }
+
+    }
+
+}
+
+template <class Basis, class VectorType>
+void FEAssembler<Basis,VectorType>::
+assembleGradientAndHessian(const VectorType& sol,
+                           Dune::BlockVector<Dune::FieldVector<double, blocksize> >& gradient,
+                           Dune::BCRSMatrix<MatrixBlock>& hessian,
+                           bool computeOccupationPattern) const
+{
+    if (computeOccupationPattern) {
+
+        Dune::MatrixIndexSet neighborsPerVertex;
+        getNeighborsPerVertex(neighborsPerVertex);
+        neighborsPerVertex.exportIdx(hessian);
+
+    }
+
+    hessian = 0;
+
+    gradient.resize(sol.size());
+    gradient = 0;
+
+    ElementIterator it    = basis_.getGridView().template begin<0,Dune::Interior_Partition>();
+    ElementIterator endit = basis_.getGridView().template end<0,Dune::Interior_Partition>  ();
+
+    for( ; it != endit; ++it ) {
+
+        const int numOfBaseFct = basis_.getLocalFiniteElement(*it).localBasis().size();
+
+        // Extract local solution
+        VectorType localSolution(numOfBaseFct);
+        VectorType localPointLoads(numOfBaseFct);
+
+        for (int i=0; i<numOfBaseFct; i++)
+            localSolution[i]   = sol[basis_.index(*it,i)];
+
+        VectorType localGradient(numOfBaseFct);
+
+        // setup local matrix and gradient
+        localStiffness_->assembleGradientAndHessian(*it, basis_.getLocalFiniteElement(*it), localSolution, localGradient);
+
+        // Add element matrix to global stiffness matrix
+        for(int i=0; i<numOfBaseFct; i++) {
+
+            int row = basis_.index(*it,i);
+
+            for (int j=0; j<numOfBaseFct; j++ ) {
+
+                int col = basis_.index(*it,j);
+                hessian[row][col] += localStiffness_->A_[i][j];
+
+            }
+        }
+
+        // Add local gradient to global gradient
+        for (int i=0; i<numOfBaseFct; i++)
+            gradient[basis_.index(*it,i)] += localGradient[i];
+
+    }
+
+}
+
+
+template <class Basis, class VectorType>
+double FEAssembler<Basis, VectorType>::
+computeEnergy(const VectorType& sol) const
+{
+    double energy = 0;
+
+    if (sol.size()!=basis_.size())
+        DUNE_THROW(Dune::Exception, "Solution vector doesn't match the basis!");
+
+    ElementIterator it    = basis_.getGridView().template begin<0,Dune::Interior_Partition>();
+    ElementIterator endIt = basis_.getGridView().template end<0,Dune::Interior_Partition>();
+
+    // Loop over all elements
+    for (; it!=endIt; ++it) {
+
+        // Number of degrees of freedom on this element
+        size_t nDofs = basis_.getLocalFiniteElement(*it).localBasis().size();
+
+        VectorType localSolution(nDofs);
+
+        for (size_t i=0; i<nDofs; i++)
+            localSolution[i]   = sol[basis_.index(*it,i)];
+
+        energy += localStiffness_->energy(*it, basis_.getLocalFiniteElement(*it), localSolution);
+
+    }
+
+    return energy;
+
+}
+#endif
diff --git a/dune/elasticity/assemblers/localadolcstiffness.hh b/dune/elasticity/assemblers/localadolcstiffness.hh
new file mode 100644
index 0000000000000000000000000000000000000000..fcc0f2dfc5071c7cc9bf758b7df5d90a5dea8ab6
--- /dev/null
+++ b/dune/elasticity/assemblers/localadolcstiffness.hh
@@ -0,0 +1,198 @@
+#ifndef DUNE_ELASTICITY_ASSEMBLERS_LOCALADOLCSTIFFNESS_HH
+#define DUNE_ELASTICITY_ASSEMBLERS_LOCALADOLCSTIFFNESS_HH
+
+#include <adolc/adouble.h>            // use of active doubles
+#include <adolc/drivers/drivers.h>    // use of "Easy to Use" drivers
+// gradient(.) and hessian(.)
+#include <adolc/interfaces.h>    // use of "Easy to Use" drivers
+#include <adolc/taping.h>             // use of taping
+
+#include <dune/fufem/utilities/adolcnamespaceinjections.hh>
+
+#include <dune/common/fmatrix.hh>
+#include <dune/istl/matrix.hh>
+
+#include <dune/elasticity/assemblers/localfestiffness.hh>
+
+//#define ADOLC_VECTOR_MODE
+
+/** \brief Assembles energy gradient and Hessian with ADOL-C (automatic differentiation)
+ */
+template<class GridView, class LocalFiniteElement, class VectorType>
+class LocalADOLCStiffness
+    : public LocalFEStiffness<GridView,LocalFiniteElement,VectorType>
+{
+    // grid types
+    typedef typename GridView::Grid::ctype DT;
+    typedef typename VectorType::value_type::field_type RT;
+    typedef typename GridView::template Codim<0>::Entity Entity;
+
+    // some other sizes
+    enum {gridDim=GridView::dimension};
+
+    // Hack
+    typedef std::vector<Dune::FieldVector<adouble,gridDim> > AVectorType;
+
+public:
+
+    //! Dimension of a tangent space
+    enum { blocksize = VectorType::value_type::dimension };
+
+    LocalADOLCStiffness(const LocalFEStiffness<GridView, LocalFiniteElement, AVectorType>* energy)
+    : localEnergy_(energy)
+    {}
+
+    /** \brief Compute the energy at the current configuration */
+    virtual RT energy (const Entity& e,
+               const LocalFiniteElement& localFiniteElement,
+               const VectorType& localConfiguration) const;
+
+    /** \brief Assemble the local stiffness matrix at the current position
+
+    This uses the automatic differentiation toolbox ADOL_C.
+    */
+    virtual void assembleGradientAndHessian(const Entity& e,
+                         const LocalFiniteElement& localFiniteElement,
+                         const VectorType& localConfiguration,
+                         VectorType& localGradient);
+
+    const LocalFEStiffness<GridView, LocalFiniteElement, AVectorType>* localEnergy_;
+
+};
+
+
+template <class GridView, class LocalFiniteElement, class VectorType>
+typename LocalADOLCStiffness<GridView, LocalFiniteElement, VectorType>::RT
+LocalADOLCStiffness<GridView, LocalFiniteElement, VectorType>::
+energy(const Entity& element,
+       const LocalFiniteElement& localFiniteElement,
+       const VectorType& localSolution) const
+{
+    double pureEnergy;
+
+    std::vector<Dune::FieldVector<adouble,blocksize> > localASolution(localSolution.size());
+
+    trace_on(1);
+
+    adouble energy = 0;
+
+    for (size_t i=0; i<localSolution.size(); i++)
+      for (size_t j=0; j<localSolution[i].size(); j++)
+        localASolution[i][j] <<= localSolution[i][j];
+
+    energy = localEnergy_->energy(element,localFiniteElement,localASolution);
+
+    energy >>= pureEnergy;
+
+    trace_off();
+
+    return pureEnergy;
+}
+
+
+
+// ///////////////////////////////////////////////////////////
+//   Compute gradient and Hessian together
+//   To compute the Hessian we need to compute the gradient anyway, so we may
+//   as well return it.  This saves assembly time.
+// ///////////////////////////////////////////////////////////
+template <class GridType, class LocalFiniteElement, class VectorType>
+void LocalADOLCStiffness<GridType, LocalFiniteElement, VectorType>::
+assembleGradientAndHessian(const Entity& element,
+                const LocalFiniteElement& localFiniteElement,
+                const VectorType& localSolution,
+                VectorType& localGradient)
+{
+    // Tape energy computation.  We may not have to do this every time, but it's comparatively cheap.
+    energy(element, localFiniteElement, localSolution);
+
+    /////////////////////////////////////////////////////////////////
+    // Compute the energy gradient
+    /////////////////////////////////////////////////////////////////
+
+    // Compute the actual gradient
+    size_t nDofs = localSolution.size();
+    size_t nDoubles = nDofs*blocksize;
+    std::vector<double> xp(nDoubles);
+    int idx=0;
+    for (size_t i=0; i<nDofs; i++)
+        for (size_t j=0; j<blocksize; j++)
+            xp[idx++] = localSolution[i][j];
+
+  // Compute gradient
+    std::vector<double> g(nDoubles);
+    gradient(1,nDoubles,xp.data(),g.data());                  // gradient evaluation
+
+    // Copy into Dune type
+    localGradient.resize(localSolution.size());
+
+    idx=0;
+    for (size_t i=0; i<nDofs; i++)
+        for (size_t j=0; j<blocksize; j++)
+            localGradient[i][j] = g[idx++];
+
+    /////////////////////////////////////////////////////////////////
+    // Compute Hessian
+    /////////////////////////////////////////////////////////////////
+
+    this->A_.setSize(nDofs,nDofs);
+
+#ifndef ADOLC_VECTOR_MODE
+    std::vector<double> v(nDoubles);
+    std::vector<double> w(nDoubles);
+
+    std::fill(v.begin(), v.end(), 0.0);
+
+    for (size_t i=0; i<nDofs; i++)
+      for (size_t ii=0; ii<blocksize; ii++)
+      {
+        // Evaluate Hessian in the direction of each vector of the orthonormal frame
+        for (size_t k=0; k<blocksize; k++)
+          v[i*blocksize + k] = (k==ii);
+
+        int rc= 3;
+        MINDEC(rc, hess_vec(1, nDoubles, xp.data(), v.data(), w.data()));
+        if (rc < 0)
+          DUNE_THROW(Dune::Exception, "ADOL-C has returned with error code " << rc << "!");
+
+        for (size_t j=0; j<nDoubles; j++)
+          this->A_[i][j/blocksize][ii][j%blocksize] = w[j];
+
+        // Make v the null vector again
+        std::fill(&v[i*blocksize], &v[(i+1)*blocksize], 0.0);
+      }
+#else
+    int n = nDoubles;
+    int nDirections = nDofs * blocksize;
+    double* tangent[nDoubles];
+    for(size_t i=0; i<nDoubles; i++)
+        tangent[i] = (double*)malloc(nDirections*sizeof(double));
+
+    double* rawHessian[nDoubles];
+    for(size_t i=0; i<nDoubles; i++)
+        rawHessian[i] = (double*)malloc(nDirections*sizeof(double));
+
+    for (int j=0; j<nDirections; j++)
+    {
+      for (int i=0; i<n; i++)
+        tangent[i][j] = 0.0;
+
+      for (int i=0; i<embeddedBlocksize; i++)
+        tangent[(j/blocksize)*embeddedBlocksize+i][j] = orthonormalFrame[j/blocksize][j%blocksize][i];
+    }
+
+    hess_mat(1,nDoubles,nDirections,xp.data(),tangent,rawHessian);
+
+    // Copy Hessian into Dune data type
+    for(size_t i=0; i<nDoubles; i++)
+      for (size_t j=0; j<nDirections; j++)
+        this->A_[j/blocksize][i/embeddedBlocksize][j%blocksize][i%embeddedBlocksize] = rawHessian[i][j];
+
+    for(size_t i=0; i<nDoubles; i++) {
+        free(rawHessian[i]);
+        free(tangent[i]);
+    }
+#endif
+}
+
+#endif
diff --git a/dune/elasticity/assemblers/localfestiffness.hh b/dune/elasticity/assemblers/localfestiffness.hh
new file mode 100644
index 0000000000000000000000000000000000000000..413d0eb21abb97c375d88ee2e5a65ff1a9629a5e
--- /dev/null
+++ b/dune/elasticity/assemblers/localfestiffness.hh
@@ -0,0 +1,58 @@
+#ifndef LOCAL_FE_STIFFNESS_HH
+#define LOCAL_FE_STIFFNESS_HH
+
+#include <dune/common/fmatrix.hh>
+#include <dune/istl/matrix.hh>
+
+
+template<class GridView, class LocalFiniteElement, class VectorType>
+class LocalFEStiffness
+{
+    // grid types
+    typedef typename GridView::Grid::ctype DT;
+    typedef typename VectorType::value_type::field_type RT;
+    typedef typename GridView::template Codim<0>::Entity Entity;
+
+    // some other sizes
+    enum {gridDim=GridView::dimension};
+
+public:
+
+    //! Dimension of a tangent space
+    enum { blocksize = VectorType::value_type::dimension };
+
+    /** \brief Assemble the local stiffness matrix at the current position
+
+    This default implementation used finite-difference approximations to compute the second derivatives
+    */
+    virtual void assembleGradientAndHessian(const Entity& e,
+                                 const LocalFiniteElement& localFiniteElement,
+                                 const VectorType& localConfiguration,
+                                 VectorType& localGradient);
+
+    /** \brief Compute the energy at the current configuration */
+    virtual RT energy (const Entity& e,
+                       const LocalFiniteElement& localFiniteElement,
+                       const VectorType& localConfiguration) const = 0;
+
+    // assembled data
+    Dune::Matrix<Dune::FieldMatrix<RT,blocksize,blocksize> > A_;
+
+};
+
+
+// ///////////////////////////////////////////////////////////
+//   Compute gradient by finite-difference approximation
+// ///////////////////////////////////////////////////////////
+template <class GridType, class LocalFiniteElement, class VectorType>
+void LocalFEStiffness<GridType, LocalFiniteElement, VectorType>::
+assembleGradientAndHessian(const Entity& element,
+                const LocalFiniteElement& localFiniteElement,
+                const VectorType& localConfiguration,
+                VectorType& localGradient)
+{
+  DUNE_THROW(Dune::NotImplemented, "!");
+}
+
+#endif
+
diff --git a/dune/elasticity/common/CMakeLists.txt b/dune/elasticity/common/CMakeLists.txt
index 09932dd0daf3bf0d89c4d0ace0b736e2c7f1ffd6..712cd0849103a43a8238624fd6fcc14459f8d80d 100644
--- a/dune/elasticity/common/CMakeLists.txt
+++ b/dune/elasticity/common/CMakeLists.txt
@@ -1,5 +1,8 @@
 install(FILES
     elasticityhelpers.hh
+    maxnormtrustregion.hh
     nonlinearelasticityproblem.hh
     nonlinearelasticityproblem.cc
+    trustregionsolver.hh
+    trustregionsolver.cc
     DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/dune/elasticity/common)
diff --git a/dune/elasticity/common/maxnormtrustregion.hh b/dune/elasticity/common/maxnormtrustregion.hh
new file mode 100644
index 0000000000000000000000000000000000000000..5145cb507137655723bfaf0967721e2d5a913034
--- /dev/null
+++ b/dune/elasticity/common/maxnormtrustregion.hh
@@ -0,0 +1,99 @@
+#ifndef DUNE_ELASTICITY_COMMON_MAXNORMTRUSTREGION_HH
+#define DUNE_ELASTICITY_COMMON_MAXNORMTRUSTREGION_HH
+
+#include <vector>
+
+#include <dune/solvers/common/boxconstraint.hh>
+
+template <int blocksize, class field_type=double>
+class MaxNormTrustRegion
+{
+public:
+
+    MaxNormTrustRegion(size_t size, field_type initialRadius)
+        : obstacles_(size)
+    {
+        set(initialRadius);
+    }
+
+    void set(field_type radius) {
+
+        radius_ = radius;
+
+        for (size_t i=0; i<obstacles_.size(); i++) {
+
+            for (int k=0; k<blocksize; k++) {
+
+                obstacles_[i].lower(k) = -radius;
+                obstacles_[i].upper(k) =  radius;
+
+            }
+
+        }
+
+    }
+
+    /** \brief Set trust region radius with a separate scaling for each vector block component
+     */
+    void set(field_type radius, const Dune::FieldVector<field_type, blocksize>& scaling) {
+
+      radius_ = radius;
+
+        for (size_t i=0; i<obstacles_.size(); i++) {
+
+            for (int k=0; k<blocksize; k++) {
+
+                obstacles_[i].lower(k) = -radius*scaling[k];
+                obstacles_[i].upper(k) =  radius*scaling[k];
+
+            }
+        }
+
+    }
+
+    /** \brief Return true if the given vector is not contained in the trust region */
+    bool violates(const Dune::BlockVector<Dune::FieldVector<double,blocksize> >& v) const
+    {
+      assert(v.size() == obstacles_.size());
+      for (size_t i=0; i<v.size(); i++)
+        for (size_t j=0; j<blocksize; j++)
+          if (v[i][j] < obstacles_[i].lower(j) or v[i][j] > obstacles_[i].upper(j))
+            return true;
+
+      return false;
+    }
+
+    field_type radius() const {
+        return radius_;
+    }
+
+    void scale(field_type factor) {
+
+        radius_ *= factor;
+
+        for (size_t i=0; i<obstacles_.size(); i++) {
+
+            for (int k=0; k<blocksize; k++) {
+
+                obstacles_[i].lower(k) *= factor;
+                obstacles_[i].upper(k) *= factor;
+
+            }
+
+        }
+
+    }
+
+    const std::vector<BoxConstraint<field_type,blocksize> >& obstacles() const {
+        return obstacles_;
+    }
+
+private:
+
+    std::vector<BoxConstraint<field_type,blocksize> > obstacles_;
+
+    field_type radius_;
+
+};
+
+#endif
diff --git a/dune/elasticity/common/trustregionsolver.cc b/dune/elasticity/common/trustregionsolver.cc
new file mode 100644
index 0000000000000000000000000000000000000000..f91551faff56e66713a6428ac84880873985f309
--- /dev/null
+++ b/dune/elasticity/common/trustregionsolver.cc
@@ -0,0 +1,397 @@
+#include <dune/common/bitsetvector.hh>
+#include <dune/common/timer.hh>
+
+#include <dune/istl/io.hh>
+
+#include <dune/functions/functionspacebases/pq1nodalbasis.hh>
+#include <dune/functions/functionspacebases/pqknodalbasis.hh>
+
+#include <dune/fufem/functionspacebases/p1nodalbasis.hh>
+#include <dune/fufem/assemblers/operatorassembler.hh>
+#include <dune/fufem/assemblers/localassemblers/laplaceassembler.hh>
+#include <dune/fufem/assemblers/localassemblers/massassembler.hh>
+#include <dune/fufem/assemblers/basisinterpolationmatrixassembler.hh>
+
+// Using a monotone multigrid as the inner solver
+#include <dune/solvers/iterationsteps/trustregiongsstep.hh>
+#include <dune/solvers/iterationsteps/mmgstep.hh>
+#include <dune/solvers/transferoperators/truncatedcompressedmgtransfer.hh>
+#if defined THIRD_ORDER || defined SECOND_ORDER
+#include <dune/gfe/pktop1mgtransfer.hh>
+#endif
+#include <dune/solvers/transferoperators/mandelobsrestrictor.hh>
+#include <dune/solvers/solvers/iterativesolver.hh>
+#include <dune/solvers/norms/twonorm.hh>
+#include <dune/solvers/norms/h1seminorm.hh>
+
+#include <dune/elasticity/common/maxnormtrustregion.hh>
+
+
+template <class BasisType, class VectorType>
+void TrustRegionSolver<BasisType,VectorType>::
+setup(const typename BasisType::GridView::Grid& grid,
+      const FEAssembler<BasisType, VectorType>* assembler,
+         const SolutionType& x,
+         const Dune::BitSetVector<blocksize>& dirichletNodes,
+         double tolerance,
+         int maxTrustRegionSteps,
+         double initialTrustRegionRadius,
+         int multigridIterations,
+         double mgTolerance,
+         int mu,
+         int nu1,
+         int nu2,
+         int baseIterations,
+         double baseTolerance)
+{
+    grid_                     = &grid;
+    assembler_                = assembler;
+    x_                        = x;
+    this->tolerance_          = tolerance;
+    maxTrustRegionSteps_      = maxTrustRegionSteps;
+    initialTrustRegionRadius_ = initialTrustRegionRadius;
+    innerIterations_          = multigridIterations;
+    innerTolerance_           = mgTolerance;
+    ignoreNodes_              = &dirichletNodes;
+
+    int numLevels = grid_->maxLevel()+1;
+
+    // ////////////////////////////////
+    //   Create a multigrid solver
+    // ////////////////////////////////
+
+#ifdef HAVE_IPOPT
+    // First create an IPOpt base solver
+    QuadraticIPOptSolver<MatrixType, CorrectionType>* baseSolver = new QuadraticIPOptSolver<MatrixType,CorrectionType>;
+    baseSolver->verbosity_ = NumProc::QUIET;
+    baseSolver->tolerance_ = baseTolerance;
+    baseSolver->linearSolverType_ = "mumps";
+#else
+    // First create a Gauss-seidel base solver
+    TrustRegionGSStep<MatrixType, CorrectionType>* baseSolverStep = new TrustRegionGSStep<MatrixType, CorrectionType>;
+
+    // Hack: the two-norm may not scale all that well, but it is fast!
+    TwoNorm<CorrectionType>* baseNorm = new TwoNorm<CorrectionType>;
+
+    ::LoopSolver<CorrectionType>* baseSolver = new ::LoopSolver<CorrectionType>(baseSolverStep,
+                                                                            baseIterations,
+                                                                            baseTolerance,
+                                                                            baseNorm,
+                                                                            Solver::QUIET);
+#endif
+
+    // Make pre and postsmoothers
+    TrustRegionGSStep<MatrixType, CorrectionType>* presmoother  = new TrustRegionGSStep<MatrixType, CorrectionType>;
+    TrustRegionGSStep<MatrixType, CorrectionType>* postsmoother = new TrustRegionGSStep<MatrixType, CorrectionType>;
+
+    MonotoneMGStep<MatrixType, CorrectionType>* mmgStep = new MonotoneMGStep<MatrixType, CorrectionType>;
+
+    mmgStep->setMGType(mu, nu1, nu2);
+    mmgStep->ignoreNodes_ = &dirichletNodes;
+    mmgStep->basesolver_        = baseSolver;
+    mmgStep->setSmoother(presmoother, postsmoother);
+    mmgStep->obstacleRestrictor_= new MandelObstacleRestrictor<CorrectionType>();
+    mmgStep->verbosity_         = Solver::QUIET;
+
+    // //////////////////////////////////////////////////////////////////////////////////////
+    //   Assemble a Laplace matrix to create a norm that's equivalent to the H1-norm
+    // //////////////////////////////////////////////////////////////////////////////////////
+
+    BasisType basis(grid.leafGridView());
+    OperatorAssembler<BasisType,BasisType> operatorAssembler(basis, basis);
+
+    LaplaceAssembler<GridType, typename BasisType::LocalFiniteElement, typename BasisType::LocalFiniteElement> laplaceStiffness;
+    typedef Dune::BCRSMatrix<Dune::FieldMatrix<double,1,1> > ScalarMatrixType;
+    ScalarMatrixType localA;
+
+    operatorAssembler.assemble(laplaceStiffness, localA);
+
+    if (h1SemiNorm_)
+        delete h1SemiNorm_;
+
+    ScalarMatrixType* A = new ScalarMatrixType(localA);
+
+    h1SemiNorm_ = new H1SemiNorm<CorrectionType>(*A);
+
+    innerSolver_ = std::shared_ptr<LoopSolver<CorrectionType> >(new ::LoopSolver<CorrectionType>(mmgStep,
+                                                                                                   innerIterations_,
+                                                                                                   innerTolerance_,
+                                                                                                   h1SemiNorm_,
+                                                                                                 Solver::REDUCED));
+
+    // //////////////////////////////////////////////////////////////////////////////////////
+    //   Assemble a mass matrix to create a norm that's equivalent to the L2-norm
+    //   This will be used to monitor the gradient
+    // //////////////////////////////////////////////////////////////////////////////////////
+
+    MassAssembler<GridType, typename BasisType::LocalFiniteElement, typename BasisType::LocalFiniteElement> massStiffness;
+    ScalarMatrixType localMassMatrix;
+
+    operatorAssembler.assemble(massStiffness, localMassMatrix);
+
+    ScalarMatrixType* massMatrix = new ScalarMatrixType(localMassMatrix);
+    l2Norm_ = std::make_shared<H1SemiNorm<CorrectionType> >(*massMatrix);
+
+    // ////////////////////////////////////////////////////////////
+    //    Create Hessian matrix and its occupation structure
+    // ////////////////////////////////////////////////////////////
+
+    hessianMatrix_ = std::auto_ptr<MatrixType>(new MatrixType);
+    Dune::MatrixIndexSet indices(grid_->size(1), grid_->size(1));
+    assembler_->getNeighborsPerVertex(indices);
+    indices.exportIdx(*hessianMatrix_);
+
+    // ////////////////////////////////////
+    //   Create the transfer operators
+    // ////////////////////////////////////
+
+    for (size_t k=0; k<mmgStep->mgTransfer_.size(); k++)
+        delete(mmgStep->mgTransfer_[k]);
+
+    ////////////////////////////////////////////////////////////////////////
+    //  The P1 space (actually P1/Q1, depending on the grid) is special:
+    //  If we work in such a space, then the multigrid hierarchy of spaces
+    //  is constructed in the usual way.  For all other space, there is
+    //  an additional restriction operator on the top of the hierarchy, which
+    //  restricts the FE space to the P1/Q1 space on the same grid.
+    //  On the lower grid levels a hierarchy of P1/Q1 spaces is used again.
+    ////////////////////////////////////////////////////////////////////////
+
+    typedef BasisType Basis;
+    bool isP1Basis = std::is_same<Basis,DuneFunctionsBasis<Dune::Functions::PQ1NodalBasis<typename Basis::GridView> > >::value
+                  || std::is_same<Basis,DuneFunctionsBasis<Dune::Functions::PQKNodalBasis<typename Basis::GridView, 1> > >::value;
+
+    if (isP1Basis)
+      mmgStep->mgTransfer_.resize(numLevels-1);
+    else
+      mmgStep->mgTransfer_.resize(numLevels);
+
+    // Here we set up the restriction of the leaf grid space into the leaf grid P1/Q1 space
+    if (not isP1Basis)
+    {
+        typedef typename TruncatedCompressedMGTransfer<CorrectionType>::TransferOperatorType TransferOperatorType;
+        P1NodalBasis<typename GridType::LeafGridView,double> p1Basis(grid_->leafGridView());
+
+        TransferOperatorType pkToP1TransferMatrix;
+        assembleBasisInterpolationMatrix<TransferOperatorType,
+                                         P1NodalBasis<typename GridType::LeafGridView,double>,
+                                         BasisType>(pkToP1TransferMatrix,p1Basis,basis);
+        mmgStep->mgTransfer_.back() = new TruncatedCompressedMGTransfer<CorrectionType>;
+        Dune::shared_ptr<TransferOperatorType> topTransferOperator = Dune::make_shared<TransferOperatorType>(pkToP1TransferMatrix);
+        dynamic_cast<TruncatedCompressedMGTransfer<CorrectionType>*>(mmgStep->mgTransfer_.back())->setMatrix(topTransferOperator);
+    }
+
+    // Now the P1/Q1 restriction operators for the remaining levels
+    for (int i=0; i<numLevels-1; i++) {
+
+        // Construct the local multigrid transfer matrix
+        TruncatedCompressedMGTransfer<CorrectionType>* newTransferOp = new TruncatedCompressedMGTransfer<CorrectionType>;
+        newTransferOp->setup(*grid_,i,i+1);
+
+        typedef typename TruncatedCompressedMGTransfer<CorrectionType>::TransferOperatorType TransferOperatorType;
+        mmgStep->mgTransfer_[i] = new TruncatedCompressedMGTransfer<CorrectionType>;
+        std::shared_ptr<TransferOperatorType> transferOperatorMatrix = Dune::make_shared<TransferOperatorType>(newTransferOp->getMatrix());
+        dynamic_cast<TruncatedCompressedMGTransfer<CorrectionType>*>(mmgStep->mgTransfer_[i])->setMatrix(transferOperatorMatrix);
+    }
+
+    // //////////////////////////////////////////////////////////
+    //   Create obstacles
+    // //////////////////////////////////////////////////////////
+
+    hasObstacle_.resize(basis.size(), true);
+    mmgStep->hasObstacle_ = &hasObstacle_;
+
+}
+
+
+template <class BasisType, class VectorType>
+void TrustRegionSolver<BasisType,VectorType>::solve()
+{
+    MonotoneMGStep<MatrixType,CorrectionType>* mgStep = NULL;
+
+    // if the inner solver is a monotone multigrid set up a max-norm trust-region
+    if (dynamic_cast<LoopSolver<CorrectionType>*>(innerSolver_.get())) {
+        mgStep = dynamic_cast<MonotoneMGStep<MatrixType,CorrectionType>*>(dynamic_cast<LoopSolver<CorrectionType>*>(innerSolver_.get())->iterationStep_);
+
+    }
+
+    BasisType basis(grid_->leafGridView());
+    MaxNormTrustRegion<blocksize> trustRegion(basis.size(), initialTrustRegionRadius_);
+
+    std::vector<BoxConstraint<field_type,blocksize> > trustRegionObstacles;
+
+    // /////////////////////////////////////////////////////
+    //   Trust-Region Solver
+    // /////////////////////////////////////////////////////
+
+    double oldEnergy = assembler_->computeEnergy(x_);
+
+    bool recomputeGradientHessian = true;
+    CorrectionType rhs;
+    MatrixType stiffnessMatrix;
+
+    for (int i=0; i<maxTrustRegionSteps_; i++) {
+
+        Dune::Timer totalTimer;
+        if (this->verbosity_ == Solver::FULL) {
+            std::cout << "----------------------------------------------------" << std::endl;
+            std::cout << "      Trust-Region Step Number: " << i
+                      << ",     radius: " << trustRegion.radius()
+                      << ",     energy: " << oldEnergy << std::endl;
+            std::cout << "----------------------------------------------------" << std::endl;
+        }
+
+        Dune::Timer gradientTimer;
+
+        if (recomputeGradientHessian) {
+
+            assembler_->assembleGradientAndHessian(x_,
+                                                   rhs,
+                                                   *hessianMatrix_,
+                                                   i==0    // assemble occupation pattern only for the first call
+                                                   );
+
+            rhs *= -1;        // The right hand side is the _negative_ gradient
+
+            // Compute gradient norm to monitor convergence
+            CorrectionType gradient = rhs;
+            for (size_t j=0; j<gradient.size(); j++)
+              for (size_t k=0; k<gradient[j].size(); k++)
+                if ((*ignoreNodes_)[j][k])
+                  gradient[j][k] = 0;
+
+            if (this->verbosity_ == Solver::FULL)
+              std::cout << "Gradient norm: " << l2Norm_->operator()(gradient) << std::endl;
+
+            if (this->verbosity_ == Solver::FULL)
+              std::cout << "Assembly took " << gradientTimer.elapsed() << " sec." << std::endl;
+
+            // Transfer matrix data
+            stiffnessMatrix = *hessianMatrix_;
+
+            recomputeGradientHessian = false;
+
+        }
+
+        CorrectionType corr(rhs.size());
+        corr = 0;
+
+        mgStep->setProblem(stiffnessMatrix, corr, rhs);
+
+        trustRegionObstacles = trustRegion.obstacles();
+        mgStep->obstacles_ = &trustRegionObstacles;
+
+        innerSolver_->preprocess();
+
+        ///////////////////////////////
+        //    Solve !
+        ///////////////////////////////
+
+        std::cout << "Solve quadratic problem..." << std::endl;
+
+        Dune::Timer solutionTimer;
+        innerSolver_->solve();
+        std::cout << "Solving the quadratic problem took " << solutionTimer.elapsed() << " seconds." << std::endl;
+
+        if (mgStep)
+            corr = mgStep->getSol();
+
+        //std::cout << "Correction: " << std::endl << corr << std::endl;
+
+        if (this->verbosity_ == NumProc::FULL)
+            std::cout << "Infinity norm of the correction: " << corr.infinity_norm() << std::endl;
+
+        if (corr.infinity_norm() < this->tolerance_) {
+            if (this->verbosity_ == NumProc::FULL)
+                std::cout << "CORRECTION IS SMALL ENOUGH" << std::endl;
+
+            if (this->verbosity_ != NumProc::QUIET)
+                std::cout << i+1 << " trust-region steps were taken." << std::endl;
+            break;
+        }
+
+        // ////////////////////////////////////////////////////
+        //   Check whether trust-region step can be accepted
+        // ////////////////////////////////////////////////////
+
+        SolutionType newIterate = x_;
+        for (size_t j=0; j<newIterate.size(); j++)
+            newIterate[j] += corr[j];
+
+        double energy    = assembler_->computeEnergy(newIterate);
+
+        // compute the model decrease
+        // It is $ m(x) - m(x+s) = -<g,s> - 0.5 <s, Hs>
+        // Note that rhs = -g
+        CorrectionType tmp(corr.size());
+        tmp = 0;
+        hessianMatrix_->umv(corr, tmp);
+        double modelDecrease = (rhs*corr) - 0.5 * (corr*tmp);
+
+        double relativeModelDecrease = modelDecrease / std::fabs(energy);
+
+        if (this->verbosity_ == NumProc::FULL) {
+            std::cout << "Absolute model decrease: " << modelDecrease
+                      << ",  functional decrease: " << oldEnergy - energy << std::endl;
+            std::cout << "Relative model decrease: " << relativeModelDecrease
+                      << ",  functional decrease: " << (oldEnergy - energy)/energy << std::endl;
+        }
+
+        assert(modelDecrease >= 0);
+
+        if (energy >= oldEnergy) {
+            if (this->verbosity_ == NumProc::FULL)
+                printf("Richtung ist keine Abstiegsrichtung!\n");
+        }
+
+        if (energy >= oldEnergy &&
+            (std::abs((oldEnergy-energy)/energy) < 1e-9 || relativeModelDecrease < 1e-9)) {
+            if (this->verbosity_ == NumProc::FULL)
+                std::cout << "Suspecting rounding problems" << std::endl;
+
+            if (this->verbosity_ != NumProc::QUIET)
+                std::cout << i+1 << " trust-region steps were taken." << std::endl;
+
+            x_ = newIterate;
+            break;
+        }
+
+        // //////////////////////////////////////////////
+        //   Check for acceptance of the step
+        // //////////////////////////////////////////////
+        if ( (oldEnergy-energy) / modelDecrease > 0.9) {
+            // very successful iteration
+
+            x_ = newIterate;
+            trustRegion.scale(2);
+
+            // current energy becomes 'oldEnergy' for the next iteration
+            oldEnergy = energy;
+
+            recomputeGradientHessian = true;
+
+        } else if ( (oldEnergy-energy) / modelDecrease > 0.01
+                    || std::abs(oldEnergy-energy) < 1e-12) {
+            // successful iteration
+            x_ = newIterate;
+
+            // current energy becomes 'oldEnergy' for the next iteration
+            oldEnergy = energy;
+
+            recomputeGradientHessian = true;
+
+        } else {
+
+            // unsuccessful iteration
+
+            // Decrease the trust-region radius
+            trustRegion.scale(0.5);
+
+            if (this->verbosity_ == NumProc::FULL)
+                std::cout << "Unsuccessful iteration!" << std::endl;
+        }
+
+        std::cout << "iteration took " << totalTimer.elapsed() << " sec." << std::endl;
+    }
+
+}
diff --git a/dune/elasticity/common/trustregionsolver.hh b/dune/elasticity/common/trustregionsolver.hh
new file mode 100644
index 0000000000000000000000000000000000000000..88dac98445e5c82a3d101db11beb4e418f470a4f
--- /dev/null
+++ b/dune/elasticity/common/trustregionsolver.hh
@@ -0,0 +1,130 @@
+#ifndef DUNE_ELASTICITY_COMMON_TRUSTREGIONSOLVER_HH
+#define DUNE_ELASTICITY_COMMON_TRUSTREGIONSOLVER_HH
+
+#include <vector>
+
+#include <dune/common/bitsetvector.hh>
+
+#include <dune/istl/bcrsmatrix.hh>
+#include <dune/istl/bvector.hh>
+
+#include <dune/solvers/common/boxconstraint.hh>
+#include <dune/solvers/norms/h1seminorm.hh>
+#include <dune/solvers/solvers/iterativesolver.hh>
+#include <dune/solvers/solvers/loopsolver.hh>
+
+#include <dune/fufem/functionspacebases/p1nodalbasis.hh>
+#include <dune/fufem/functionspacebases/p2nodalbasis.hh>
+#include <dune/fufem/functionspacebases/p3nodalbasis.hh>
+
+#include <dune/elasticity/assemblers/feassembler.hh>
+
+/** \brief Trust-region solver  */
+template <class BasisType, class VectorType>
+class TrustRegionSolver
+    : public IterativeSolver<VectorType,
+                             Dune::BitSetVector<VectorType::value_type::dimension> >
+{
+    typedef typename BasisType::GridView::Grid GridType;
+
+    const static int blocksize = VectorType::value_type::dimension;
+
+    const static int gridDim = GridType::dimension;
+
+    // Centralize the field type here
+    typedef double field_type;
+
+    // Some types that I need
+    typedef Dune::BCRSMatrix<Dune::FieldMatrix<field_type, blocksize, blocksize> > MatrixType;
+    typedef Dune::BlockVector<Dune::FieldVector<field_type, blocksize> >           CorrectionType;
+    typedef VectorType                                                             SolutionType;
+
+public:
+
+    TrustRegionSolver()
+        : IterativeSolver<VectorType, Dune::BitSetVector<blocksize> >(0,100,NumProc::FULL),
+          hessianMatrix_(std::auto_ptr<MatrixType>(NULL)), h1SemiNorm_(NULL)
+    {}
+
+    /** \brief Set up the solver using a monotone multigrid method as the inner solver */
+    void setup(const GridType& grid,
+               const FEAssembler<BasisType,VectorType>* assembler,
+               const SolutionType& x,
+               const Dune::BitSetVector<blocksize>& dirichletNodes,
+               double tolerance,
+               int maxTrustRegionSteps,
+               double initialTrustRegionRadius,
+               int multigridIterations,
+               double mgTolerance,
+               int mu,
+               int nu1,
+               int nu2,
+               int baseIterations,
+               double baseTolerance);
+
+    void setIgnoreNodes(const Dune::BitSetVector<blocksize>& ignoreNodes)
+    {
+        ignoreNodes_ = &ignoreNodes;
+        Dune::shared_ptr<LoopSolver<CorrectionType> > loopSolver = std::dynamic_pointer_cast<LoopSolver<CorrectionType> >(innerSolver_);
+        assert(loopSolver);
+        loopSolver->iterationStep_->ignoreNodes_ = ignoreNodes_;
+    }
+
+    void solve();
+
+    void setInitialSolution(const SolutionType& x) DUNE_DEPRECATED {
+        x_ = x;
+    }
+
+    void setInitialIterate(const SolutionType& x) {
+        x_ = x;
+    }
+
+    SolutionType getSol() const {return x_;}
+
+protected:
+
+    /** \brief The grid */
+    const GridType* grid_;
+
+    /** \brief The solution vector */
+    SolutionType x_;
+
+    /** \brief The initial trust-region radius in the maximum-norm */
+    double initialTrustRegionRadius_;
+
+    /** \brief Maximum number of trust-region steps */
+    int maxTrustRegionSteps_;
+
+    /** \brief Maximum number of multigrid iterations */
+    int innerIterations_;
+
+    /** \brief Error tolerance of the multigrid QP solver */
+    double innerTolerance_;
+
+    /** \brief Hessian matrix */
+    std::auto_ptr<MatrixType> hessianMatrix_;
+
+    /** \brief The assembler for the material law */
+    const FEAssembler<BasisType, VectorType>* assembler_;
+
+    /** \brief The solver for the quadratic inner problems */
+    std::shared_ptr<Solver> innerSolver_;
+
+    /** \brief Contains 'true' everywhere -- the trust-region is bounded */
+    Dune::BitSetVector<blocksize> hasObstacle_;
+
+    /** \brief The Dirichlet nodes */
+    const Dune::BitSetVector<blocksize>* ignoreNodes_;
+
+    /** \brief The norm used to measure multigrid convergence */
+    H1SemiNorm<CorrectionType>* h1SemiNorm_;
+
+    /** \brief An L2-norm, really.  The H1SemiNorm class is badly named */
+    std::shared_ptr<H1SemiNorm<CorrectionType> > l2Norm_;
+
+};
+
+#include "trustregionsolver.cc"
+
+#endif
diff --git a/dune/elasticity/materials/CMakeLists.txt b/dune/elasticity/materials/CMakeLists.txt
index 1748a21886665c48ead2b193b4b67fcbce9f19ec..0eda46c467e377ed475a4d89e67b4ec29701a494 100644
--- a/dune/elasticity/materials/CMakeLists.txt
+++ b/dune/elasticity/materials/CMakeLists.txt
@@ -1,6 +1,12 @@
 install(FILES
     adolcmaterial.hh
     adolcneohookeanmaterial.hh
+    exphenckyenergy.hh
     geomexactstvenantkirchhoffmaterial.hh
+    henckyenergy.hh
     neohookeanmaterial.hh
+    neohookeenergy.hh
+    neumannenergy.hh
+    stvenantkirchhoffenergy.hh
+    sumenergy.hh
     DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/dune/elasticity/materials)
diff --git a/dune/elasticity/materials/exphenckyenergy.hh b/dune/elasticity/materials/exphenckyenergy.hh
new file mode 100644
index 0000000000000000000000000000000000000000..d684227803067be93313c03fb6bc20aabb17aaca
--- /dev/null
+++ b/dune/elasticity/materials/exphenckyenergy.hh
@@ -0,0 +1,150 @@
+#ifndef DUNE_ELASTICITY_MATERIALS_EXPHENCKYENERGY_HH
+#define DUNE_ELASTICITY_MATERIALS_EXPHENCKYENERGY_HH
+
+#include <dune/common/fmatrix.hh>
+#include <dune/common/fmatrixev.hh>
+
+#include <dune/geometry/quadraturerules.hh>
+
+#include <dune/elasticity/assemblers/localfestiffness.hh>
+
+namespace Dune {
+
+template<class GridView, class LocalFiniteElement, class field_type=double>
+class ExpHenckyEnergy
+  : public LocalFEStiffness<GridView,
+                            LocalFiniteElement,
+                            std::vector<Dune::FieldVector<field_type,3> > >
+{
+  // grid types
+  typedef typename GridView::Grid::ctype DT;
+  typedef typename GridView::template Codim<0>::Entity Entity;
+
+  // some other sizes
+  enum {gridDim=GridView::dimension};
+  enum {dim=GridView::dimension};
+
+public:
+
+  /** \brief Constructor with a set of material parameters
+   * \param parameters The material parameters
+   */
+  ExpHenckyEnergy(const Dune::ParameterTree& parameters)
+  {
+    // Lame constants
+    mu_     = parameters.get<double>("mu");
+    lambda_ = parameters.get<double>("lambda");
+    kappa_  = (2.0 * mu_ + 3.0 * lambda_) / 3.0;
+
+    // Exponential Hencky constants with naming convention from Neff, Ghiba and Lankeit,
+    // "The exponentiated Hencky-logarithmic strain energy. Part I: Constitutive issues
+    //  and rank one convexity"
+    k_      = parameters.get<double>("k");
+    khat_   = parameters.get<double>("khat");
+
+    // weights for distortion and dilatation parts of the energy
+    alpha_  = mu_ / k_;
+    beta_   = kappa_ / (2.0 * khat_);
+  }
+
+  /** \brief Assemble the energy for a single element */
+  field_type energy (const Entity& e,
+                     const LocalFiniteElement& localFiniteElement,
+                     const std::vector<Dune::FieldVector<field_type, gridDim> >& localConfiguration) const;
+
+  /** \brief Lame constants */
+  field_type mu_, lambda_, kappa_, k_, khat_, alpha_, beta_;
+};
+
+template <class GridView, class LocalFiniteElement, class field_type>
+field_type
+ExpHenckyEnergy<GridView, LocalFiniteElement, field_type>::
+energy(const Entity& element,
+       const LocalFiniteElement& localFiniteElement,
+       const std::vector<Dune::FieldVector<field_type, gridDim> >& localConfiguration) const
+{
+  assert(element.type() == localFiniteElement.type());
+  typedef typename GridView::template Codim<0>::Entity::Geometry Geometry;
+
+  // store gradients of shape functions and base functions
+  std::vector<Dune::FieldMatrix<DT,1,gridDim> > referenceGradients(localFiniteElement.size());
+  std::vector<Dune::FieldVector<DT,gridDim> > gradients(localFiniteElement.size());
+
+  int quadOrder = (element.type().isSimplex()) ? localFiniteElement.localBasis().order()
+                                               : localFiniteElement.localBasis().order() * gridDim;
+
+  const Dune::QuadratureRule<DT, gridDim>& quad
+      = Dune::QuadratureRules<DT, gridDim>::rule(element.type(), quadOrder);
+
+  field_type distortionEnergy = 0, dilatationEnergy = 0;
+  for (size_t pt=0; pt<quad.size(); pt++)
+  {
+    // Local position of the quadrature point
+    const Dune::FieldVector<DT,gridDim>& quadPos = quad[pt].position();
+
+    const DT integrationElement = element.geometry().integrationElement(quadPos);
+
+    const typename Geometry::JacobianInverseTransposed& jacobianInverseTransposed = element.geometry().jacobianInverseTransposed(quadPos);
+
+    DT weight = quad[pt].weight() * integrationElement;
+
+    // get gradients of shape functions
+    localFiniteElement.localBasis().evaluateJacobian(quadPos, referenceGradients);
+
+    // compute gradients of base functions
+    for (size_t i=0; i<gradients.size(); ++i)
+      jacobianInverseTransposed.mv(referenceGradients[i][0], gradients[i]);
+
+    Dune::FieldMatrix<field_type,gridDim,gridDim> derivative(0);
+    for (size_t i=0; i<gradients.size(); i++)
+      for (int j=0; j<gridDim; j++)
+        derivative[j].axpy(localConfiguration[i][j], gradients[i]);
+
+    /////////////////////////////////////////////////////////
+    // compute C = F^T F
+    /////////////////////////////////////////////////////////
+
+    Dune::FieldMatrix<field_type,gridDim,gridDim> C(0);
+    for (int i=0; i<gridDim; i++)
+      for (int j=0; j<gridDim; j++)
+        for (int k=0; k<gridDim; k++)
+          C[i][j] += derivative[k][i] * derivative[k][j];
+
+    //////////////////////////////////////////////////////////
+    //  Eigenvalues of FTF
+    //////////////////////////////////////////////////////////
+
+    // compute eigenvalues of C, i.e., squared singular values \sigma_i of F
+    Dune::FieldVector<field_type, dim> sigmaSquared;
+    FMatrixHelp::eigenValues(C, sigmaSquared);
+
+    // logarithms of the singular values of F, i.e., eigenvalues of U
+    std::array<field_type, dim> lnSigma;
+    for (int i = 0; i < dim; i++)
+      lnSigma[i] = 0.5 * log(sigmaSquared[i]);
+
+    field_type trLogU = 0;
+    for (int i = 0; i < dim; i++)
+      trLogU += lnSigma[i];
+
+    field_type normDevLogUSquared = 0;
+    for (int i = 0; i < dim; i++)
+    {
+      // the deviator shifts the
+      field_type lnDevSigma_i = lnSigma[i] - (1.0 / double(dim)) * trLogU;
+      normDevLogUSquared += lnDevSigma_i * lnDevSigma_i;
+    }
+
+    // Add the local energy density
+    distortionEnergy += weight * alpha_ * exp(khat_ * normDevLogUSquared);
+    dilatationEnergy += weight * beta_  * exp(k_ * (trLogU * trLogU));
+  }
+
+  return distortionEnergy + dilatationEnergy;
+}
+
+}  // namespace Dune
+
+#endif   //#ifndef DUNE_ELASTICITY_MATERIALS_EXPHENCKYENERGY_HH
+
+
diff --git a/dune/elasticity/materials/henckyenergy.hh b/dune/elasticity/materials/henckyenergy.hh
new file mode 100644
index 0000000000000000000000000000000000000000..c603a3451957e107a5360b894082f1e1a192ab27
--- /dev/null
+++ b/dune/elasticity/materials/henckyenergy.hh
@@ -0,0 +1,139 @@
+#ifndef DUNE_ELASTICITY_MATERIALS_HENCKYENERGY_HH
+#define DUNE_ELASTICITY_MATERIALS_HENCKYENERGY_HH
+
+#include <dune/common/fmatrix.hh>
+#include <dune/common/fmatrixev.hh>
+
+#include <dune/geometry/quadraturerules.hh>
+
+#include <dune/elasticity/assemblers/localfestiffness.hh>
+
+namespace Dune {
+
+template<class GridView, class LocalFiniteElement, class field_type=double>
+class HenckyEnergy
+    : public LocalFEStiffness<GridView,LocalFiniteElement,std::vector<Dune::FieldVector<field_type,3> > >
+{
+    // grid types
+    typedef typename GridView::Grid::ctype DT;
+    typedef typename GridView::template Codim<0>::Entity Entity;
+
+    // some other sizes
+    enum {gridDim=GridView::dimension};
+    enum {dim=GridView::dimension};
+
+public:  // for testing
+
+    /** \brief Constructor with a set of material parameters
+     * \param parameters The material parameters
+     */
+    HenckyEnergy(const Dune::ParameterTree& parameters)
+    {
+      // Lame constants
+      mu_ = parameters.template get<double>("mu");
+      lambda_ = parameters.template get<double>("lambda");
+      kappa_  = (2.0 * mu_ + 3.0 * lambda_) / 3.0;
+    }
+
+    /** \brief Assemble the energy for a single element */
+    field_type energy (const Entity& e,
+               const LocalFiniteElement& localFiniteElement,
+               const std::vector<Dune::FieldVector<field_type,gridDim> >& localConfiguration) const;
+
+    /** \brief Lame constants */
+    double mu_, lambda_, kappa_;
+};
+
+template <class GridView, class LocalFiniteElement, class field_type>
+field_type
+HenckyEnergy<GridView,LocalFiniteElement,field_type>::
+energy(const Entity& element,
+       const LocalFiniteElement& localFiniteElement,
+       const std::vector<Dune::FieldVector<field_type,gridDim> >& localConfiguration) const
+{
+    assert(element.type() == localFiniteElement.type());
+    typedef typename GridView::template Codim<0>::Entity::Geometry Geometry;
+
+    field_type energy = 0;
+
+    // store gradients of shape functions and base functions
+    std::vector<Dune::FieldMatrix<DT,1,gridDim> > referenceGradients(localFiniteElement.size());
+    std::vector<Dune::FieldVector<DT,gridDim> > gradients(localFiniteElement.size());
+
+    int quadOrder = (element.type().isSimplex()) ? localFiniteElement.localBasis().order()
+                                                 : localFiniteElement.localBasis().order() * gridDim;
+
+    const Dune::QuadratureRule<DT, gridDim>& quad
+        = Dune::QuadratureRules<DT, gridDim>::rule(element.type(), quadOrder);
+
+    field_type distortionEnergy = 0, dilatationEnergy = 0;
+    for (size_t pt=0; pt<quad.size(); pt++) {
+
+        // Local position of the quadrature point
+        const Dune::FieldVector<DT,gridDim>& quadPos = quad[pt].position();
+
+        const DT integrationElement = element.geometry().integrationElement(quadPos);
+
+        const typename Geometry::JacobianInverseTransposed& jacobianInverseTransposed = element.geometry().jacobianInverseTransposed(quadPos);
+
+        DT weight = quad[pt].weight() * integrationElement;
+
+        // get gradients of shape functions
+        localFiniteElement.localBasis().evaluateJacobian(quadPos, referenceGradients);
+
+        // compute gradients of base functions
+        for (size_t i=0; i<gradients.size(); ++i)
+          jacobianInverseTransposed.mv(referenceGradients[i][0], gradients[i]);
+
+        Dune::FieldMatrix<field_type,gridDim,gridDim> derivative(0);
+        for (size_t i=0; i<gradients.size(); i++)
+          for (int j=0; j<gridDim; j++)
+            derivative[j].axpy(localConfiguration[i][j], gradients[i]);
+
+        /////////////////////////////////////////////////////////
+        // compute C = F^TF
+        /////////////////////////////////////////////////////////
+
+        Dune::FieldMatrix<field_type,gridDim,gridDim> C(0);
+        for (int i=0; i<gridDim; i++)
+          for (int j=0; j<gridDim; j++)
+            for (int k=0; k<gridDim; k++)
+              C[i][j] += derivative[k][i] * derivative[k][j];
+
+        //////////////////////////////////////////////////////////
+        //  Eigenvalues of C = F^TF
+        //////////////////////////////////////////////////////////
+
+        Dune::FieldVector<field_type,dim> sigmaSquared;
+        FMatrixHelp::eigenValues(C, sigmaSquared);
+
+        // logarithms of the singular values of F, i.e., eigenvalues of U
+        std::array<field_type, dim> lnSigma;
+        for (int i = 0; i < dim; i++)
+          lnSigma[i] = 0.5 * log(sigmaSquared[i]);
+
+        field_type trLogU = 0;
+        for (int i = 0; i < dim; i++)
+          trLogU += lnSigma[i];
+
+        field_type normDevLogUSquared = 0;
+        for (int i = 0; i < dim; i++)
+        {
+          // the deviator shifts the spectrum by the trace
+          field_type lnDevSigma_i = lnSigma[i] - (1.0 / double(dim)) * trLogU;
+          normDevLogUSquared += lnDevSigma_i * lnDevSigma_i;
+        }
+
+        // Add the local energy density
+        distortionEnergy += weight * mu_ * normDevLogUSquared;
+        dilatationEnergy += weight * 0.5 * kappa_ * (trLogU * trLogU);
+    }
+
+    return distortionEnergy + dilatationEnergy;
+}
+
+}  // namespace Dune
+
+#endif   //#ifndef DUNE_ELASTICITY_MATERIALS_HENCKYENERGY_HH
+
+
diff --git a/dune/elasticity/materials/neohookeenergy.hh b/dune/elasticity/materials/neohookeenergy.hh
new file mode 100644
index 0000000000000000000000000000000000000000..09118d0756843d5ad9853ff31b72d876e18854b8
--- /dev/null
+++ b/dune/elasticity/materials/neohookeenergy.hh
@@ -0,0 +1,138 @@
+#ifndef DUNE_ELASTICITY_MATERIALS_NEOHOOKEENERGY_HH
+#define DUNE_ELASTICITY_MATERIALS_NEOHOOKEENERGY_HH
+
+#include <dune/common/fmatrix.hh>
+#include <dune/common/fmatrixev.hh>
+
+#include <dune/geometry/quadraturerules.hh>
+
+#include <dune/elasticity/assemblers/localfestiffness.hh>
+
+namespace Dune {
+
+template<class GridView, class LocalFiniteElement, class field_type=double>
+class NeoHookeEnergy
+  : public LocalFEStiffness<GridView,
+                            LocalFiniteElement,
+                            std::vector<Dune::FieldVector<field_type,3> > >
+{
+  // grid types
+  typedef typename GridView::Grid::ctype DT;
+  typedef typename GridView::template Codim<0>::Entity Entity;
+
+  // some other sizes
+  enum {gridDim=GridView::dimension};
+  enum {dim=GridView::dimension};
+
+public:
+
+  /** \brief Constructor with a set of material parameters
+   * \param parameters The material parameters
+   */
+  NeoHookeEnergy(const Dune::ParameterTree& parameters)
+  {
+    // Lame constants
+    mu_     = parameters.get<double>("mu");
+    lambda_ = parameters.get<double>("lambda");
+    kappa_  = (2.0 * mu_ + 3.0 * lambda_) / 3.0;
+  }
+
+  /** \brief Assemble the energy for a single element */
+  field_type energy (const Entity& e,
+                     const LocalFiniteElement& localFiniteElement,
+                     const std::vector<Dune::FieldVector<field_type, gridDim> >& localConfiguration) const;
+
+  /** \brief Lame constants */
+  field_type mu_, lambda_, kappa_;
+};
+
+template <class GridView, class LocalFiniteElement, class field_type>
+field_type
+NeoHookeEnergy<GridView, LocalFiniteElement, field_type>::
+energy(const Entity& element,
+       const LocalFiniteElement& localFiniteElement,
+       const std::vector<Dune::FieldVector<field_type, gridDim> >& localConfiguration) const
+{
+  assert(element.type() == localFiniteElement.type());
+  typedef typename GridView::template Codim<0>::Entity::Geometry Geometry;
+
+  field_type distortionEnergy = 0, dilatationEnergy = 0;
+
+  // store gradients of shape functions and base functions
+  std::vector<Dune::FieldMatrix<DT,1,gridDim> > referenceGradients(localFiniteElement.size());
+  std::vector<Dune::FieldVector<DT,gridDim> > gradients(localFiniteElement.size());
+
+  int quadOrder = (element.type().isSimplex()) ? localFiniteElement.localBasis().order()
+                                               : localFiniteElement.localBasis().order() * gridDim;
+
+  const Dune::QuadratureRule<DT, gridDim>& quad
+      = Dune::QuadratureRules<DT, gridDim>::rule(element.type(), quadOrder);
+
+  for (size_t pt=0; pt<quad.size(); pt++)
+  {
+    // Local position of the quadrature point
+    const Dune::FieldVector<DT,gridDim>& quadPos = quad[pt].position();
+
+    const DT integrationElement = element.geometry().integrationElement(quadPos);
+
+    const typename Geometry::JacobianInverseTransposed& jacobianInverseTransposed = element.geometry().jacobianInverseTransposed(quadPos);
+
+    DT weight = quad[pt].weight() * integrationElement;
+
+    // get gradients of shape functions
+    localFiniteElement.localBasis().evaluateJacobian(quadPos, referenceGradients);
+
+    // compute gradients of base functions
+    for (size_t i=0; i<gradients.size(); ++i)
+      jacobianInverseTransposed.mv(referenceGradients[i][0], gradients[i]);
+
+    Dune::FieldMatrix<field_type,gridDim,gridDim> derivative(0);
+    for (size_t i=0; i<gradients.size(); i++)
+      for (int j=0; j<gridDim; j++)
+        derivative[j].axpy(localConfiguration[i][j], gradients[i]);
+
+    /////////////////////////////////////////////////////////
+    // compute C = F^T F
+    /////////////////////////////////////////////////////////
+
+    Dune::FieldMatrix<field_type,gridDim,gridDim> C(0);
+    for (int i=0; i<gridDim; i++)
+      for (int j=0; j<gridDim; j++)
+        for (int k=0; k<gridDim; k++)
+          C[i][j] += derivative[k][i] * derivative[k][j];
+
+    //////////////////////////////////////////////////////////
+    //  Eigenvalues of FTF
+    //////////////////////////////////////////////////////////
+
+    // eigenvalues of C, i.e., squared singular values \sigma_i of F
+    Dune::FieldVector<field_type, dim> sigmaSquared;
+    FMatrixHelp::eigenValues(C, sigmaSquared);
+
+    // singular values of F, i.e., eigenvalues of U
+    std::array<field_type, dim> sigma;
+    for (int i = 0; i < dim; i++)
+      sigma[i] = std::sqrt(sigmaSquared[i]);
+
+    field_type detC = 1.0;
+    for (int i = 0; i < dim; i++)
+      detC *= sigmaSquared[i];
+    field_type detF = std::sqrt(detC);
+
+    // \tilde{C} = \tilde{F}^T\tilde{F} = \frac{1}{\det{F}^{2/3}}C
+    field_type trCTilde = 0;
+    for (int i = 0; i < dim; i++)
+      trCTilde += sigmaSquared[i];
+    trCTilde /= pow(detC, 1.0/3.0);
+
+    // Add the local energy density
+    distortionEnergy += weight * (0.5 * mu_)    * (trCTilde - 3);
+    dilatationEnergy += weight * (0.5 * kappa_) * ((detF - 1) * (detF - 1));
+  }
+
+  return distortionEnergy + dilatationEnergy;
+}
+
+}  // namespace Dune
+
+#endif   //#ifndef DUNE_ELASTICITY_MATERIALS_NEOHOOKEENERGY_HH
diff --git a/dune/elasticity/materials/neumannenergy.hh b/dune/elasticity/materials/neumannenergy.hh
new file mode 100644
index 0000000000000000000000000000000000000000..539418d3fb6edaa0bbb0d719da292680594d04ce
--- /dev/null
+++ b/dune/elasticity/materials/neumannenergy.hh
@@ -0,0 +1,101 @@
+#ifndef DUNE_ELASTICITY_MATERIALS_NEUMANNENERGY_HH
+#define DUNE_ELASTICITY_MATERIALS_NEUMANNENERGY_HH
+
+#include <dune/common/fmatrix.hh>
+#include <dune/common/fmatrixev.hh>
+
+#include <dune/geometry/quadraturerules.hh>
+
+#include <dune/fufem/functions/virtualgridfunction.hh>
+#include <dune/fufem/boundarypatch.hh>
+
+#include <dune/elasticity/assemblers/localfestiffness.hh>
+
+namespace Dune {
+
+template<class GridView, class LocalFiniteElement, class field_type=double>
+class NeumannEnergy
+: public LocalFEStiffness<GridView,LocalFiniteElement,std::vector<Dune::FieldVector<field_type,3> > >
+{
+ // grid types
+  typedef typename GridView::ctype ctype;
+  typedef typename GridView::template Codim<0>::Entity Entity;
+
+  enum {dim=GridView::dimension};
+
+public:
+
+  /** \brief Constructor with a set of material parameters
+   * \param parameters The material parameters
+   */
+  NeumannEnergy(const BoundaryPatch<GridView>* neumannBoundary,
+                const Dune::VirtualFunction<Dune::FieldVector<ctype,dim>, Dune::FieldVector<double,3> >* neumannFunction)
+  : neumannBoundary_(neumannBoundary),
+    neumannFunction_(neumannFunction)
+  {}
+
+  /** \brief Assemble the energy for a single element */
+  field_type energy (const Entity& element,
+                     const LocalFiniteElement& localFiniteElement,
+                     const std::vector<Dune::FieldVector<field_type,dim> >& localConfiguration) const
+  {
+    assert(element.type() == localFiniteElement.type());
+
+    field_type energy = 0;
+
+    for (auto&& it : intersections(neumannBoundary_->gridView(), element)) {
+
+      if (not neumannBoundary_ or not neumannBoundary_->contains(it))
+        continue;
+
+      int quadOrder = localFiniteElement.localBasis().order();
+
+      const Dune::QuadratureRule<ctype, dim-1>& quad
+          = Dune::QuadratureRules<ctype, dim-1>::rule(it.type(), quadOrder);
+
+      for (size_t pt=0; pt<quad.size(); pt++) {
+
+        // Local position of the quadrature point
+        const Dune::FieldVector<ctype,dim>& quadPos = it.geometryInInside().global(quad[pt].position());
+
+        const auto integrationElement = it.geometry().integrationElement(quad[pt].position());
+
+        // The value of the local function
+        std::vector<Dune::FieldVector<ctype,1> > shapeFunctionValues;
+        localFiniteElement.localBasis().evaluateFunction(quadPos, shapeFunctionValues);
+
+        Dune::FieldVector<field_type,dim> value(0);
+        for (size_t i=0; i<localFiniteElement.size(); i++)
+          for (int j=0; j<dim; j++)
+            value[j] += shapeFunctionValues[i] * localConfiguration[i][j];
+
+        // Value of the Neumann data at the current position
+        Dune::FieldVector<double,3> neumannValue;
+
+        if (dynamic_cast<const VirtualGridViewFunction<GridView,Dune::FieldVector<double,3> >*>(neumannFunction_))
+            dynamic_cast<const VirtualGridViewFunction<GridView,Dune::FieldVector<double,3> >*>(neumannFunction_)->evaluateLocal(element, quadPos, neumannValue);
+        else
+            neumannFunction_->evaluate(it.geometry().global(quad[pt].position()), neumannValue);
+
+        // Only translational dofs are affected by the Neumann force
+        for (size_t i=0; i<neumannValue.size(); i++)
+          energy += (neumannValue[i] * value[i]) * quad[pt].weight() * integrationElement;
+
+      }
+
+    }
+
+    return energy;
+  }
+
+private:
+  /** \brief The Neumann boundary */
+  const BoundaryPatch<GridView>* neumannBoundary_;
+
+  /** \brief The function implementing the Neumann data */
+  const Dune::VirtualFunction<Dune::FieldVector<double,dim>, Dune::FieldVector<double,3> >* neumannFunction_;
+};
+
+}
+
+#endif   //#ifndef DUNE_ELASTICITY_MATERIALS_NEUMANNENERGY_HH
diff --git a/dune/elasticity/materials/stvenantkirchhoffenergy.hh b/dune/elasticity/materials/stvenantkirchhoffenergy.hh
new file mode 100644
index 0000000000000000000000000000000000000000..3edae89aa2f2a716452cb49e05f3a22b95722684
--- /dev/null
+++ b/dune/elasticity/materials/stvenantkirchhoffenergy.hh
@@ -0,0 +1,128 @@
+#ifndef DUNE_ELASTICITY_MATERIALS_STVENANTKIRCHHOFFENERGY_HH
+#define DUNE_ELASTICITY_MATERIALS_STVENANTKIRCHHOFFENERGY_HH
+
+#include <dune/common/fmatrix.hh>
+#include <dune/geometry/quadraturerules.hh>
+
+#include <dune/elasticity/assemblers/localfestiffness.hh>
+
+namespace Dune {
+
+template<class GridView, class LocalFiniteElement, class field_type=double>
+class StVenantKirchhoffEnergy
+  : public LocalFEStiffness<GridView,LocalFiniteElement,std::vector<Dune::FieldVector<field_type,GridView::dimension> > >
+{
+    // grid types
+    typedef typename GridView::Grid::ctype DT;
+    typedef typename GridView::template Codim<0>::Entity Entity;
+
+    // some other sizes
+    enum {gridDim=GridView::dimension};
+    enum {dim=GridView::dimension};
+
+public:
+
+    /** \brief Constructor with a set of material parameters
+     * \param parameters The material parameters
+     */
+    StVenantKirchhoffEnergy(const Dune::ParameterTree& parameters)
+    {
+      // Lame constants
+      mu_ = parameters.template get<double>("mu");
+      lambda_ = parameters.template get<double>("lambda");
+    }
+
+    /** \brief Assemble the energy for a single element */
+    field_type energy (const Entity& e,
+               const LocalFiniteElement& localFiniteElement,
+               const std::vector<Dune::FieldVector<field_type,gridDim> >& localConfiguration) const;
+
+    /** \brief Lame constants */
+    double mu_, lambda_;
+};
+
+template <class GridView, class LocalFiniteElement, class field_type>
+field_type
+StVenantKirchhoffEnergy<GridView,LocalFiniteElement,field_type>::
+energy(const Entity& element,
+       const LocalFiniteElement& localFiniteElement,
+       const std::vector<Dune::FieldVector<field_type,gridDim> >& localConfiguration) const
+{
+    assert(element.type() == localFiniteElement.type());
+    typedef typename GridView::template Codim<0>::Entity::Geometry Geometry;
+
+    field_type energy = 0;
+
+    // store gradients of shape functions and base functions
+    std::vector<Dune::FieldMatrix<DT,1,gridDim> > referenceGradients(localFiniteElement.localBasis().size());
+    std::vector<Dune::FieldVector<DT,gridDim> > gradients(localFiniteElement.localBasis().size());
+
+    int quadOrder = (element.type().isSimplex()) ? localFiniteElement.localBasis().order()
+                                                 : localFiniteElement.localBasis().order() * gridDim;
+
+    const Dune::QuadratureRule<DT, gridDim>& quad
+        = Dune::QuadratureRules<DT, gridDim>::rule(element.type(), quadOrder);
+
+    for (size_t pt=0; pt<quad.size(); pt++) {
+
+        // Local position of the quadrature point
+        const Dune::FieldVector<DT,gridDim>& quadPos = quad[pt].position();
+
+        const DT integrationElement = element.geometry().integrationElement(quadPos);
+
+        const typename Geometry::JacobianInverseTransposed& jacobianInverseTransposed = element.geometry().jacobianInverseTransposed(quadPos);
+
+        DT weight = quad[pt].weight() * integrationElement;
+
+        // get gradients of shape functions
+        localFiniteElement.localBasis().evaluateJacobian(quadPos, referenceGradients);
+
+        // compute gradients of base functions
+        for (size_t i=0; i<gradients.size(); ++i) {
+
+          // transform gradients
+          jacobianInverseTransposed.mv(referenceGradients[i][0], gradients[i]);
+
+        }
+
+        Dune::FieldMatrix<field_type,gridDim,gridDim> derivative(0);
+        for (size_t i=0; i<gradients.size(); i++)
+          for (int j=0; j<gridDim; j++)
+            derivative[j].axpy(localConfiguration[i][j], gradients[i]);
+
+        /////////////////////////////////////////////////////////
+        // compute strain E = 1/2 *( F^T F - I)
+        /////////////////////////////////////////////////////////
+
+        Dune::FieldMatrix<field_type,gridDim,gridDim> FTF(0);
+        for (int i=0; i<gridDim; i++)
+          for (int j=0; j<gridDim; j++)
+            for (int k=0; k<gridDim; k++)
+              FTF[i][j] += derivative[k][i] * derivative[k][j];
+
+        Dune::FieldMatrix<field_type,dim,dim> E = FTF;
+        for (int i=0; i<dim; i++)
+          E[i][i] -= 1.0;
+        E *= 0.5;
+
+        /////////////////////////////////////////////////////////
+        //  Compute energy
+        /////////////////////////////////////////////////////////
+
+        field_type trE = field_type(0);
+        for (int i=0; i<dim; i++)
+          trE += E[i][i];
+
+        // The trace of E^2, happens to be its squared Frobenius norm
+        field_type trESquared = E.frobenius_norm2();
+
+        energy += weight * mu_ * trESquared + weight * 0.5 * lambda_ * trE * trE;
+
+    }
+
+    return energy;
+}
+
+}  // namespace Dune
+
+#endif   //#ifndef DUNE_ELASTICITY_MATERIALS_STVENANTKIRCHHOFFENERGY_HH
diff --git a/dune/elasticity/materials/sumenergy.hh b/dune/elasticity/materials/sumenergy.hh
new file mode 100644
index 0000000000000000000000000000000000000000..c316a163bea1918606b85e2662e3500f10a1b7ce
--- /dev/null
+++ b/dune/elasticity/materials/sumenergy.hh
@@ -0,0 +1,55 @@
+#ifndef DUNE_ELASTICITY_MATERIALS_SUMENERGY_HH
+#define DUNE_ELASTICITY_MATERIALS_SUMENERGY_HH
+
+#include <dune/common/fmatrix.hh>
+#include <dune/common/fmatrixev.hh>
+
+#include <dune/geometry/quadraturerules.hh>
+
+#include <dune/fufem/functions/virtualgridfunction.hh>
+#include <dune/fufem/boundarypatch.hh>
+
+#include <dune/elasticity/assemblers/localfestiffness.hh>
+
+namespace Dune {
+
+template<class GridView, class LocalFiniteElement, class field_type=double>
+class SumEnergy
+: public LocalFEStiffness<GridView,LocalFiniteElement,std::vector<Dune::FieldVector<field_type,3> > >
+{
+ // grid types
+  typedef typename GridView::ctype ctype;
+  typedef typename GridView::template Codim<0>::Entity Entity;
+
+  enum {dim=GridView::dimension};
+
+public:
+
+  /** \brief Constructor with a set of material parameters
+   * \param parameters The material parameters
+   */
+  SumEnergy(std::shared_ptr<LocalFEStiffness<GridView,LocalFiniteElement,std::vector<Dune::FieldVector<field_type,3> > > > a,
+            std::shared_ptr<LocalFEStiffness<GridView,LocalFiniteElement,std::vector<Dune::FieldVector<field_type,3> > > > b)
+  : a_(a),
+    b_(b)
+  {}
+
+  /** \brief Assemble the energy for a single element */
+  field_type energy (const Entity& element,
+                     const LocalFiniteElement& localFiniteElement,
+                     const std::vector<Dune::FieldVector<field_type,dim> >& localConfiguration) const
+  {
+    return a_->energy(element, localFiniteElement, localConfiguration)
+         + b_->energy(element, localFiniteElement, localConfiguration);
+  }
+
+private:
+
+  std::shared_ptr<LocalFEStiffness<GridView,LocalFiniteElement,std::vector<Dune::FieldVector<field_type,3> > > > a_;
+
+  std::shared_ptr<LocalFEStiffness<GridView,LocalFiniteElement,std::vector<Dune::FieldVector<field_type,3> > > > b_;
+};
+
+}
+
+#endif   //#ifndef DUNE_ELASTICITY_MATERIALS_SUMENERGY_HH
diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt
new file mode 100644
index 0000000000000000000000000000000000000000..c2f16aac3d19f049b7639290c0df860b51008297
--- /dev/null
+++ b/src/CMakeLists.txt
@@ -0,0 +1,11 @@
+set(programs finite-strain-elasticity)
+
+foreach(_program ${programs})
+  add_executable(${_program} ${_program}.cc)
+  add_dune_adolc_flags(${_program})
+  add_dune_ipopt_flags(${_program})
+  add_dune_pythonlibs_flags(${_program})
+  add_dune_ug_flags(${_program})
+  add_dune_mpi_flags(${_program})
+#  target_compile_options(${_program} PRIVATE "-fpermissive")
+endforeach()
diff --git a/src/finite-strain-elasticity-bending.parset b/src/finite-strain-elasticity-bending.parset
new file mode 100644
index 0000000000000000000000000000000000000000..c05f7828ba35399306471e2b6a575116f0dcb599
--- /dev/null
+++ b/src/finite-strain-elasticity-bending.parset
@@ -0,0 +1,85 @@
+#############################################
+#  Grid parameters
+#############################################
+
+structuredGrid = true
+lower = 0 -0.5 -0.5
+upper = 10 0.5 0.5
+elements = 10 1 1
+
+# Number of grid levels
+numLevels = 2
+
+#############################################
+#  Solver parameters
+#############################################
+
+# Number of homotopy steps for the Dirichlet boundary conditions
+numHomotopySteps = 1
+
+# Tolerance of the trust region solver
+tolerance = 1e-6
+
+# Max number of steps of the trust region solver
+maxTrustRegionSteps = 500
+
+# Initial trust-region radius
+initialTrustRegionRadius = 0.1
+
+# Number of multigrid iterations per trust-region step
+numIt = 1000
+
+# 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 = 100
+
+# Tolerance of the multigrid solver
+mgTolerance = 1e-7
+
+# Tolerance of the base grid solver
+baseTolerance = 1e-8
+
+############################
+#   Material parameters
+############################
+
+energy = stvenantkirchhoff
+
+## For the Wriggers L-shape example
+[materialParameters]
+
+# Lame parameters
+# corresponds to E = 71240 N/mm^2, nu=0.31
+# However, we use units N/m^2
+mu = 2.7191e+6
+lambda = 4.4364e+6
+
+[]
+
+#############################################
+#  Boundary values
+#############################################
+
+problem = wriggers-l-shape
+
+###  Python predicate specifying all Dirichlet grid vertices
+# x is the vertex coordinate
+dirichletVerticesPredicate = "x[0] < 0.001"
+
+###  Python predicate specifying all Dirichlet grid vertices
+# x is the vertex coordinate
+neumannVerticesPredicate = "x[1] < -0.239"
+
+###  Neumann values, if needed
+neumannValues = 0 0 1e4
+
+# Initial deformation
+initialDeformation = "[x[0], x[1], x[2]]"
diff --git a/src/finite-strain-elasticity.cc b/src/finite-strain-elasticity.cc
new file mode 100644
index 0000000000000000000000000000000000000000..586872b2e6954898f6957fd991ee334352f1fe83
--- /dev/null
+++ b/src/finite-strain-elasticity.cc
@@ -0,0 +1,362 @@
+#include <config.h>
+
+// Includes for the ADOL-C automatic differentiation library
+// Need to come before (almost) all others.
+#include <adolc/adouble.h>
+#include <adolc/drivers/drivers.h>    // use of "Easy to Use" drivers
+#include <adolc/taping.h>
+
+#include <dune/fufem/utilities/adolcnamespaceinjections.hh>
+
+#include <dune/common/bitsetvector.hh>
+#include <dune/common/parametertree.hh>
+#include <dune/common/parametertreeparser.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/functions/functionspacebases/pq2nodalbasis.hh>
+#include <dune/functions/functionspacebases/interpolate.hh>
+#include <dune/functions/gridfunctions/discretescalarglobalbasisfunction.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/solvers/solvers/iterativesolver.hh>
+#include <dune/solvers/norms/energynorm.hh>
+
+#include <dune/elasticity/common/trustregionsolver.hh>
+#include <dune/elasticity/assemblers/localadolcstiffness.hh>
+#include <dune/elasticity/assemblers/feassembler.hh>
+#include <dune/elasticity/materials/stvenantkirchhoffenergy.hh>
+#include <dune/elasticity/materials/henckyenergy.hh>
+#include <dune/elasticity/materials/exphenckyenergy.hh>
+#include <dune/elasticity/materials/neohookeenergy.hh>
+#include <dune/elasticity/materials/neumannenergy.hh>
+#include <dune/elasticity/materials/sumenergy.hh>
+
+// grid dimension
+const int dim = 3;
+
+using namespace Dune;
+
+/** \brief A constant vector-valued function, for simple Neumann boundary values */
+struct NeumannFunction
+    : public Dune::VirtualFunction<FieldVector<double,dim>, FieldVector<double,3> >
+{
+  NeumannFunction(const FieldVector<double,3> values,
+                  double homotopyParameter)
+  : values_(values),
+    homotopyParameter_(homotopyParameter)
+  {}
+
+  void evaluate(const FieldVector<double, dim>& 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
+  Dune::MPIHelper& mpiHelper = MPIHelper::instance(argc, argv);
+
+  // Start Python interpreter
+  Python::start();
+  Python::Reference main = Python::import("__main__");
+  Python::run("import math");
+
+  //feenableexcept(FE_INVALID);
+  Python::runStream()
+        << std::endl << "import sys"
+        << std::endl << "sys.path.append('/home/sander/dune/dune-gfe/')"
+        << std::endl;
+
+  typedef BlockVector<FieldVector<double,dim> > SolutionType;
+
+  // parse data file
+  ParameterTree parameterSet;
+
+  ParameterTreeParser::readINITree(argv[1], parameterSet);
+
+  ParameterTreeParser::readOptions(argc, argv, parameterSet);
+
+  // read solver settings
+  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 double initialTrustRegionRadius = parameterSet.get<double>("initialTrustRegionRadius");
+  const int multigridIterations         = parameterSet.get<int>("numIt");
+  const int nu1                         = parameterSet.get<int>("nu1");
+  const int nu2                         = parameterSet.get<int>("nu2");
+  const int mu                          = parameterSet.get<int>("mu");
+  const int baseIterations              = parameterSet.get<int>("baseIt");
+  const double mgTolerance              = parameterSet.get<double>("mgTolerance");
+  const double baseTolerance            = parameterSet.get<double>("baseTolerance");
+  std::string resultPath                = parameterSet.get("resultPath", "");
+
+  // ///////////////////////////////////////
+  //    Create the grid
+  // ///////////////////////////////////////
+  typedef UGGrid<dim> GridType;
+
+  shared_ptr<GridType> grid;
+
+  FieldVector<double,dim> lower(0), upper(1);
+
+  if (parameterSet.get<bool>("structuredGrid")) {
+
+    lower = parameterSet.get<FieldVector<double,dim> >("lower");
+    upper = parameterSet.get<FieldVector<double,dim> >("upper");
+
+    array<unsigned int,dim> elements = parameterSet.get<array<unsigned int,dim> >("elements");
+    grid = StructuredGridFactory<GridType>::createCubeGrid(lower, upper, elements);
+
+  } else {
+    std::string path                = parameterSet.get<std::string>("path");
+    std::string gridFile            = parameterSet.get<std::string>("gridFile");
+    grid = shared_ptr<GridType>(GmshReader<GridType>::read(path + "/" + gridFile));
+  }
+
+  grid->globalRefine(numLevels-1);
+
+  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::PQ2NodalBasis<GridView> FEBasis;
+  FEBasis feBasis(gridView);
+
+  // dune-fufem-style FE basis for the transition from dune-fufem to dune-functions
+  typedef DuneFunctionsBasis<FEBasis> FufemFEBasis;
+  FufemFEBasis fufemFEBasis(feBasis);
+
+  // /////////////////////////////////////////
+  //   Read Dirichlet values
+  // /////////////////////////////////////////
+
+  BitSetVector<1> dirichletVertices(gridView.size(dim), false);
+  BitSetVector<1> neumannVertices(gridView.size(dim), 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(")");
+  PythonFunction<FieldVector<double,dim>, bool> pythonDirichletVertices(Python::evaluate(lambda));
+
+  // Same for the Neumann boundary
+  lambda = std::string("lambda x: (") + parameterSet.get<std::string>("neumannVerticesPredicate", "0") + std::string(")");
+  PythonFunction<FieldVector<double,dim>, bool> pythonNeumannVertices(Python::evaluate(lambda));
+
+  for (auto&& v : vertices(gridView))
+  {
+    bool isDirichlet;
+    pythonDirichletVertices.evaluate(v.geometry().corner(0), isDirichlet);
+    dirichletVertices[indexSet.index(v)] = isDirichlet;
+
+    bool isNeumann;
+    pythonNeumannVertices.evaluate(v.geometry().corner(0), isNeumann);
+    neumannVertices[indexSet.index(v)] = 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> dirichletNodes(feBasis.indexSet().size(), false);
+  constructBoundaryDofs(dirichletBoundary,fufemFEBasis,dirichletNodes);
+
+  BitSetVector<1> neumannNodes(feBasis.indexSet().size(), false);
+  constructBoundaryDofs(neumannBoundary,fufemFEBasis,neumannNodes);
+
+  BitSetVector<dim> dirichletDofs(feBasis.indexSet().size(), false);
+  for (size_t i=0; i<feBasis.indexSet().size(); i++)
+    if (dirichletNodes[i][0])
+      for (int j=0; j<dim; j++)
+        dirichletDofs[i][j] = true;
+
+  // //////////////////////////
+  //   Initial iterate
+  // //////////////////////////
+
+  SolutionType x(feBasis.indexSet().size());
+
+  lambda = std::string("lambda x: (") + parameterSet.get<std::string>("initialDeformation") + std::string(")");
+  PythonFunction<FieldVector<double,dim>, FieldVector<double,3> > pythonInitialDeformation(Python::evaluate(lambda));
+
+  ::Functions::interpolate(fufemFEBasis, x, pythonInitialDeformation);
+
+  ////////////////////////////////////////////////////////
+  //   Main homotopy loop
+  ////////////////////////////////////////////////////////
+
+  // Output initial iterate (of homotopy loop)
+  SolutionType identity;
+  Dune::Functions::interpolate(feBasis, identity, [](FieldVector<double,dim> x){ return x; });
+
+  SolutionType displacement = x;
+  displacement -= identity;
+
+  Dune::Functions::DiscreteScalarGlobalBasisFunction<FEBasis,SolutionType> displacementFunction(feBasis,displacement);
+  auto localDisplacementFunction = localFunction(displacementFunction);
+
+  //  We need to subsample, because VTK cannot natively display real second-order functions
+  SubsamplingVTKWriter<GridView> vtkWriter(gridView,2);
+  vtkWriter.addVertexData(localDisplacementFunction, VTK::FieldInfo("displacement", VTK::FieldInfo::Type::scalar, 3));
+  vtkWriter.write(resultPath + "finite-strain_homotopy_0");
+
+  for (int i=0; i<numHomotopySteps; i++)
+  {
+    double homotopyParameter = (i+1)*(1.0/numHomotopySteps);
+    if (mpiHelper.rank()==0)
+      std::cout << "Homotopy step: " << i << ",    parameter: " << homotopyParameter << std::endl;
+
+
+    // ////////////////////////////////////////////////////////////
+    //   Create an assembler for the energy functional
+    // ////////////////////////////////////////////////////////////
+
+    const ParameterTree& materialParameters = parameterSet.sub("materialParameters");
+    shared_ptr<NeumannFunction> neumannFunction;
+    if (parameterSet.hasKey("neumannValues"))
+      neumannFunction = make_shared<NeumannFunction>(parameterSet.get<FieldVector<double,3> >("neumannValues"),
+                                                     homotopyParameter);
+
+    std::cout << "Neumann values: " << parameterSet.get<FieldVector<double,3> >("neumannValues") << std::endl;
+
+    if (mpiHelper.rank() == 0)
+    {
+      std::cout << "Material parameters:" << std::endl;
+      materialParameters.report();
+    }
+
+    // Assembler using ADOL-C
+    std::cout << "Selected energy is: " << parameterSet.get<std::string>("energy") << std::endl;
+    std::shared_ptr<LocalFEStiffness<GridView,
+                                     FEBasis::LocalView::Tree::FiniteElement,
+                                     std::vector<Dune::FieldVector<adouble, 3> > > > elasticEnergy;
+
+    if (parameterSet.get<std::string>("energy") == "stvenantkirchhoff")
+      elasticEnergy = std::make_shared<StVenantKirchhoffEnergy<GridView,
+                                                               FEBasis::LocalView::Tree::FiniteElement,
+                                                               adouble> >(materialParameters);
+
+    if (parameterSet.get<std::string>("energy") == "neohooke")
+      elasticEnergy = std::make_shared<NeoHookeEnergy<GridView,
+                                                      FEBasis::LocalView::Tree::FiniteElement,
+                                                      adouble> >(materialParameters);
+
+    if (parameterSet.get<std::string>("energy") == "hencky")
+      elasticEnergy = std::make_shared<HenckyEnergy<GridView,
+                                                    FEBasis::LocalView::Tree::FiniteElement,
+                                                    adouble> >(materialParameters);
+
+    if (parameterSet.get<std::string>("energy") == "exphencky")
+      elasticEnergy = std::make_shared<ExpHenckyEnergy<GridView,
+                                                       FEBasis::LocalView::Tree::FiniteElement,
+                                                       adouble> >(materialParameters);
+
+    if(!elasticEnergy)
+      DUNE_THROW(Exception, "Error: Selected energy not available!");
+
+    auto neumannEnergy = std::make_shared<NeumannEnergy<GridView,
+                                                        FEBasis::LocalView::Tree::FiniteElement,
+                                                        adouble> >(&neumannBoundary,neumannFunction.get());
+
+    SumEnergy<GridView,
+              FEBasis::LocalView::Tree::FiniteElement,
+              adouble> totalEnergy(elasticEnergy, neumannEnergy);
+
+    LocalADOLCStiffness<GridView,
+                        FEBasis::LocalView::Tree::FiniteElement,
+                        SolutionType> localADOLCStiffness(&totalEnergy);
+
+    FEAssembler<FufemFEBasis,SolutionType> assembler(fufemFEBasis, &localADOLCStiffness);
+
+    // /////////////////////////////////////////////////
+    //   Create a Riemannian trust-region solver
+    // /////////////////////////////////////////////////
+
+    TrustRegionSolver<FufemFEBasis,SolutionType> solver;
+    solver.setup(*grid,
+                 &assembler,
+                 x,
+                 dirichletDofs,
+                 tolerance,
+                 maxTrustRegionSteps,
+                 initialTrustRegionRadius,
+                 multigridIterations,
+                 mgTolerance,
+                 mu, nu1, nu2,
+                 baseIterations,
+                 baseTolerance
+                );
+
+    ////////////////////////////////////////////////////////
+    //   Set Dirichlet values
+    ////////////////////////////////////////////////////////
+
+    // The python class that implements the Dirichlet values
+    Python::Reference dirichletValuesClass = Python::import(parameterSet.get<std::string>("problem") + "-dirichlet-values");
+
+    // The member method 'DirichletValues' of that class
+    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> >   dirichletValues(dirichletValuesPythonObject.get("dirichletValues"));
+
+    ::Functions::interpolate(fufemFEBasis, x, dirichletValues, dirichletDofs);
+
+    // /////////////////////////////////////////////////////
+    //   Solve!
+    // /////////////////////////////////////////////////////
+
+    solver.setInitialIterate(x);
+    solver.solve();
+
+    x = solver.getSol();
+
+    /////////////////////////////////
+    //   Output result
+    /////////////////////////////////
+
+    // Compute the displacement
+    auto displacement = x;
+    displacement -= identity;
+
+    Dune::Functions::DiscreteScalarGlobalBasisFunction<FEBasis,SolutionType> displacementFunction(feBasis,displacement);
+    auto localDisplacementFunction = localFunction(displacementFunction);
+
+    //  We need to subsample, because VTK cannot natively display real second-order functions
+    SubsamplingVTKWriter<GridView> vtkWriter(gridView,2);
+    vtkWriter.addVertexData(localDisplacementFunction, VTK::FieldInfo("displacement", VTK::FieldInfo::Type::scalar, 3));
+    vtkWriter.write(resultPath + "finite-strain_homotopy_" + std::to_string(i+1));
+  }
+
+} catch (Exception e) {
+    std::cout << e << std::endl;
+}