diff --git a/dune/elasticity/materials/neohookeanmaterial.hh b/dune/elasticity/materials/neohookeanmaterial.hh
new file mode 100644
index 0000000000000000000000000000000000000000..c9f82935050a8189d03c6821ee85f0bf192b2b54
--- /dev/null
+++ b/dune/elasticity/materials/neohookeanmaterial.hh
@@ -0,0 +1,181 @@
+#ifndef NEO_HOOKIAN_MATERIAL
+#define NEO_HOOKIAN_MATERIAL
+
+#include <dune/fufem/quadraturerules/quadraturerulecache.hh>
+
+#include <dune/fufem/symmetrictensor.hh>
+#include <dune/fufem/functions/virtualgridfunction.hh>
+#include <dune/fufem/functions/basisgridfunction.hh>
+
+#include <dune/elasticity/assemblers/neohookefunctionalassembler.hh>
+#include <dune/elasticity/assemblers/neohookeoperatorassembler.hh>
+
+/** \brief Class representing a hyperelastic Neo-hookian material.
+ *
+ * \tparam Basis Global basis that is used for the spatial discretization.
+ * (This is not nice but I need a LocalFiniteElement for the local Hessian assembler :-( )
+ */
+template <class Basis>
+class NeoHookeMaterial
+{
+public:
+ typedef typename Basis::GridView::Grid GridType;
+ typedef typename Basis::LocalFiniteElement Lfe;
+
+ typedef NeoHookeFunctionalAssembler<GridType> LocalLinearization;
+ typedef NeoHookeOperatorAssembler<GridType, Lfe, Lfe> LocalHessian;
+ typedef Basis GlobalBasis;
+
+private:
+ typedef typename GridType::ctype ctype;
+
+ static const int dim = GridType::dimension;
+ static const int dimworld = GridType::dimensionworld;
+
+ typedef typename GridType::template Codim<0>::Geometry::LocalCoordinate LocalCoordinate;
+ typedef typename GridType::template Codim<0>::LeafIterator ElementIterator;
+
+public:
+ NeoHookeMaterial() :
+ lambda_(1.0),
+ mu_(0.3)
+ {}
+
+ NeoHookeMaterial(const Basis& basis, ctype E, ctype nu) :
+ basis_(basis)
+ {
+ lambda_ = E*nu / ((1+nu)*(1-2*nu));
+ mu_ = E / (2*(1+nu));
+
+ localLinearization_ = Dune::shared_ptr<LocalLinearization>(new LocalLinearization(lambda_,mu_));
+ localHessian_ = Dune::shared_ptr<LocalHessian>(new LocalHessian(lambda_,mu_));
+ }
+
+
+ void setup(ctype E, ctype nu, const Basis& basis)
+ {
+ lambda_ = E*nu / ((1+nu)*(1-2*nu));
+ mu_ = E / (2*(1+nu));
+
+ if (localLinearization_)
+ localLinearization_.reset();
+ if (localHessian_)
+ localHessian_.reset();
+
+ localLinearization_ = Dune::shared_ptr<LocalLinearization>(new LocalLinearization(lambda_,mu_));
+ localHessian_ = Dune::shared_ptr<LocalHessian>(new LocalHessian(lambda_,mu_));
+
+ basis_ = &basis;
+ }
+
+ //! Evaluate the strain energy
+ template <class CoeffType>
+ ctype energy(const CoeffType& coeff)
+ {
+ // make grid function
+ BasisGridFunction<Basis,CoeffType> configuration(*basis_,coeff);
+
+ ctype energy=0;
+ const GridType& grid = configuration.grid();
+
+ ElementIterator eIt = grid.template leafbegin<0>();
+ ElementIterator eItEnd = grid.template leafend<0>();
+
+ for (;eIt != eItEnd; ++eIt) {
+
+ // TODO Get proper quadrature rule
+ // get quadrature rule
+ const int order = (eIt->type().isSimplex()) ? 4 : 4*dim;
+
+ const Dune::template QuadratureRule<ctype, dim>& quad = QuadratureRuleCache<ctype, dim>::rule(order);
+
+ // loop over quadrature points
+ for (size_t pt=0; pt < quad.size(); ++pt) {
+
+ // get quadrature point
+ const LocalCoordinate& quadPos = quad[pt].position();
+
+ // get integration factor
+ const ctype integrationElement = eIt->geometry().integrationElement(quadPos);
+
+ // evaluate displacement gradient at the quadrature point
+ typename BasisGridFunction<Basis,CoeffType>::DerivativeType localConfGrad;
+
+ if (configuration.isDefinedOn(*eIt))
+ configuration.evaluateDerivativeLocal(*eIt, quadPos, localConfGrad);
+ else
+ configuration.evaluateDerivative(eIt->geometry().global(quadPos),localConfGrad);
+
+ SymmetricTensor<dim> strain;
+ computeStrain(localConfGrad,strain);
+
+ // the trace
+ ctype trE(0);
+ for (int i=0; i<dim; i++)
+ trE += strain(i,i);
+
+ // make deformation gradient out of the discplacement
+ for (int i=0;i<dim;i++)
+ localConfGrad[i][i] += 1;
+
+ // evaluate the derminante of the deformation gradient
+ const ctype J = localConfGrad.determinant();
+
+ ctype z = quad[pt].weight()*integrationElement;
+
+#ifdef LAURSEN
+ energy += z*(0.25*lambda_*(J*J-1)-(lambda_*0.5+mu_)*std::log(J)+mu_*trE);
+#else
+ energy += z*(0.5*lambda_*(J-1)*(J-1)-2*mu_*std::log(J)+mu_*trE);
+#endif
+ }
+ }
+
+ return 0.5*energy;
+ }
+
+ //! Return the local assembler of the first derivative of the strain energy
+ LocalLinearization& firstDerivative() {return *localLinearization_;}
+
+ //! Return the local assembler of the second derivative of the strain energy
+ LocalHessian& secondDerivative() {return *localHessian_;}
+
+ //! Return the global basis
+ const Basis& basis() {return *basis_;}
+
+private:
+ //! First derivative of the strain energy
+ Dune::shared_ptr<LocalLinearization> localLinearization_;
+
+ //! Second derivative of the strain energy
+ Dune::shared_ptr<LocalHessian> localHessian_;
+
+ //! Global basis used for the spatial discretization
+ const Basis* basis_;
+
+ //! Lame constant
+ ctype lambda_;
+ //! Lame constant
+ ctype mu_;
+
+
+ //! Compute nonlinear strain tensor from the deformation gradient.
+ void computeStrain(const Dune::FieldMatrix<ctype, dim, dim>& grad, SymmetricTensor<dim>& strain) const {
+ strain = 0;
+ for (int i=0; i<dim ; ++i) {
+ strain(i,i) +=grad[i][i];
+
+ for (int k=0;k<dim;++k)
+ strain(i,i) += 0.5*grad[k][i]*grad[k][i];
+
+ for (int j=i+1; j<dim; ++j) {
+ strain(i,j) += 0.5*(grad[i][j] + grad[j][i]);
+ for (int k=0;k<dim;++k)
+ strain(i,j) += 0.5*grad[k][i]*grad[k][j];
+ }
+ }
+ }
+
+};
+
+#endif