Experimental Monney-Rivlin Material

I don't remember anymore if it was working perfectly when I implemented it

dune/elasticity/assemblers/mooneyrivlinfunctionalassembler.hh

+#ifndef MOONEY_RIVLIN_FUNCTIONAL_ASSEMBLER_HH
+#define MOONEY_RIVLIN_FUNCTIONAL_ASSEMBLER_HH
+
+#include <dune/common/fvector.hh>
+#include <dune/istl/bvector.hh>
+
+#include <dune/fufem/quadraturerules/quadraturerulecache.hh>
+#include "dune/fufem/assemblers/localfunctionalassembler.hh"
+#include <dune/fufem/symmetrictensor.hh>
+#include <dune/fufem/functions/virtualgridfunction.hh>
+
+#include <dune/elasticity/common/elasticityhelpers.hh>
+
+/** \brief Local assembler for the linearization of a nonlinear Mooney-Rivlin energy at a displacement u
+ *
+ *    \f[
+ *      W(u)= a tr(E(u)) + b tr(E(u))^2 + c ||E(u)||^2 + d J^2 + e J^-k,  k>=2  
+ *    \f]
+ *
+ *   where
+ *      - \f$E\f$: the nonlinear strain tensor
+ *      - \f$tr \f$: the trace operator
+ *      - \f$J\f$ the determinant of the deformation gradient
+ *      - \f$\a\f$,...,\f$\e\f$ material parameters
+ *      - \f$k\f$ controlls orientation preservation (k \approx \floor(1/(1-2nu) -1))
+*/  
+template <class GridType, class TrialLocalFE>
+class MooneyRivlinFunctionalAssembler : 
+    public LocalFunctionalAssembler<GridType,TrialLocalFE, Dune::FieldVector<typename GridType::ctype ,GridType::dimension> >
+
+{
+    static const int dim = GridType::dimension;
+    typedef typename GridType::ctype ctype;
+    typedef typename Dune::FieldVector<ctype,GridType::dimension> T;
+
+public:
+    typedef typename LocalFunctionalAssembler<GridType,TrialLocalFE,T>::Element Element;
+    typedef typename LocalFunctionalAssembler<GridType,TrialLocalFE,T>::LocalVector LocalVector;
+    typedef VirtualGridFunction<GridType, T > GridFunction;
+
+    //! Create assembler from material parameters and discplacement
+    MooneyRivlinFunctionalAssembler(const Dune::shared_ptr<GridFunction> displacement,
+                                    ctype a, ctype b, ctype c, ctype d, ctype e, int k) :
+        a_(a),b_(b),c_(c),d_(d),e_(e),k_(k),
+        displacement_(displacement)
+    {}
+
+    //! Create assembler from material parameters
+    MooneyRivlinFunctionalAssembler(ctype a, ctype b, ctype c, ctype d, ctype e, int k) :
+        a_(a),b_(b),c_(c),d_(d),e_(e),k_(k)
+    {}
+
+    void assemble(const Element& element, LocalVector& localVector, const TrialLocalFE& tFE) const
+    {
+
+        typedef typename Dune::FieldMatrix<ctype,dim,dim> FMdimdim;
+        typedef typename TrialLocalFE::Traits::LocalBasisType::Traits::JacobianType JacobianType;    
+        typedef typename GridType::template Codim<0>::Geometry::LocalCoordinate LocalCoordinate; 
+
+        int rows = tFE.localBasis().size();
+
+        localVector.resize(rows);
+        localVector = 0.0;
+
+        // get quadrature rule
+        const int order = (element.type().isSimplex()) 
+            ? 6*(tFE.localBasis().order())
+            : 6*(tFE.localBasis().order()*dim);
+
+        // get quadrature rule
+        const Dune::template QuadratureRule<ctype, dim>& quad = QuadratureRuleCache<ctype, dim>::rule(element.type(), order, IsRefinedLocalFiniteElement<TrialLocalFE>::value(tFE) );
+
+        // store gradients of shape functions and base functions
+        std::vector<JacobianType> referenceGradients(tFE.localBasis().size());
+        std::vector<T> gradients(tFE.localBasis().size());
+
+        // store geometry mapping of the entity
+        const typename ElementPointer::Entity::Geometry geometry = element.geometry();
+
+        // loop over quadrature points
+        for (size_t pt=0; pt < quad.size(); ++pt) {
+
+            // get quadrature point
+            const LocalCoordinate& quadPos = quad[pt].position();
+
+            // get transposed inverse of Jacobian of transformation
+            const FMdimdim& invJacobian = geometry.jacobianInverseTransposed(quadPos);
+
+            // get integration factor
+            const ctype integrationElement = geometry.integrationElement(quadPos);
+
+            // get gradients of shape functions
+            tFE.localBasis().evaluateJacobian(quadPos, referenceGradients);
+
+            // compute gradients of base functions
+            for (size_t i=0; i<gradients.size(); ++i) {
+
+                // transform gradients
+                gradients[i] = 0.0;
+                invJacobian.umv(referenceGradients[i][0], gradients[i]);
+
+            }
+
+            // evaluate the displacement gradient at the quadrature point
+            typename GridFunction::DerivativeType localConfGrad;
+            if (displacement_->isDefinedOn(element))
+                displacement_->evaluateDerivativeLocal(element, quadPos, localConfGrad);
+            else
+                displacement_->evaluateDerivative(geometry.global(quadPos),localConfGrad);
+
+            // compute linearization of the determinante of the deformation gradient
+            FMdimdim linDefDet;
+            Dune::Elasticity::linearisedDefDet(localConfGrad,linDefDet);
+
+            // the trace of the strain tensor
+            SymmetricTensor<dim> strain;
+            computeStrain(localConfGrad,strain);
+            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();
+
+            //compute linearized strain for all shapefunctions
+            std::vector<Dune::array<SymmetricTensor<dim>,dim> > linStrain(rows);
+
+            for (int i=0; i<rows; i++)
+                for (int k=0; k<dim; k++) {
+
+                    // The deformation gradient of the shape function
+                    FMdimdim deformationGradient(0);
+                    deformationGradient[k] = gradients[i];
+
+                    // The linearized strain is the symmetric product of defGrad and (Id+localConf)
+                    linStrain[i][k]=symmetricProduct(deformationGradient,localConfGrad);
+                }
+
+
+            // collect terms 
+            FMdimdim fu(0);
+
+            // the linearization of the trace is just deformation gradient
+            fu.axpy(a_+2*b_*trE,localConfGrad);
+
+            // linearization of the compressibility function
+            fu.axpy(2*d_*J-e_*k_*std::pow(J,-k_-1),linDefDet);
+
+            ctype z = quad[pt].weight()*integrationElement;
+
+            // add vector entries
+            for(int row=0; row<rows; row++) {
+                fu.usmv(z,gradients[row],localVector[row]);
+
+                // linearization of the tensor contraction
+                for (int rcomp=0;rcomp<dim;rcomp++)
+                    localVector[row][rcomp] +=(strain*linStrain[row][rcomp])*z*2*c_;
+            }
+        }
+    }
+
+    /** \brief Set new configuration. In Newton iterations this needs to be assembled more than one time.*/
+    void setConfiguration(Dune::shared_ptr<GridFunction> newDisplacement) {
+        displacement_ = newDisplacement;
+    }
+
+private:
+    //! Material parameters
+    ctype a_; ctype b_; ctype c_; ctype d_; ctype e_;
+    //! Compressibility
+    int k_;
+
+    /** \brief The configuration at which the functional is evaluated.*/
+    Dune::shared_ptr<GridFunction> displacement_;
+
+    //! 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];
+            }
+        }
+    }
+
+    /** \brief Compute the symmetric product of two matrices, i.e.
+     *      0.5*(A^T * B + B^T * A )
+     */
+    SymmetricTensor<dim> symmetricProduct(const Dune::FieldMatrix<ctype, dim, dim>& a, const Dune::FieldMatrix<ctype, dim, dim>& b) const {
+
+        SymmetricTensor<dim> result(0.0);
+        for (int i=0;i<dim; i++)
+            for (int j=i;j<dim;j++)
+                for (int k=0;k<dim;k++)
+                    result(i,j)+= a[k][i]*b[k][j]+b[k][i]*a[k][j];   
+        result *= 0.5;                
+        
+        return result;
+    }
+};
+
+#endif
+
+

dune/elasticity/assemblers/mooneyrivlinoperatorassembler.hh

+#ifndef MOONEY_RIVLIN_OPERATOR_ASSEMBLER_HH
+#define MOONEY_RIVLIN_OPERATOR_ASSEMBLER_HH
+
+
+#include <dune/common/fvector.hh>
+#include <dune/common/fmatrix.hh>
+
+#include <dune/istl/matrix.hh>
+
+#include <dune/fufem/quadraturerules/quadraturerulecache.hh>
+#include <dune/fufem/assemblers/localoperatorassembler.hh>
+#include <dune/fufem/staticmatrixtools.hh>
+#include <dune/fufem/symmetrictensor.hh>
+#include <dune/fufem/functions/virtualgridfunction.hh>
+
+#include <dune/elasticity/common/elasticityhelpers.hh>
+
+/** \brief Assemble the second derivative of a Mooney-Rivlin material.
+ *
+ *    \f[
+ *      W(u)= a tr(E(u)) + b tr(E(u))^2 + c ||E(u)||^2 + d J^2 + e J^-k,  k>=2  
+ *    \f]
+ *
+ *   where
+ *      - \f$E\f$: the nonlinear strain tensor
+ *      - \f$tr \f$: the trace operator
+ *      - \f$J\f$ the determinant of the deformation gradient
+ *      - \f$\a\f$,...,\f$\e\f$ material parameters
+ *      - \f$k\f$ controlls orientation preservation (k \approx \floor(1/(1-2nu) -1))
+*/  
+template <class GridType, class TrialLocalFE, class AnsatzLocalFE>
+class MooneyRivlinOperatorAssembler 
+    : public LocalOperatorAssembler < GridType, TrialLocalFE, AnsatzLocalFE, Dune::FieldMatrix<typename GridType::ctype ,GridType::dimension,GridType::dimension> >
+{
+
+
+    static const int dim = GridType::dimension;
+    typedef typename GridType::ctype ctype;
+    
+    typedef VirtualGridFunction<GridType, Dune::FieldVector<ctype,GridType::dimension> > GridFunction;
+
+public:
+    typedef typename Dune::FieldMatrix<ctype,GridType::dimension,GridType::dimension> T;
+
+    typedef typename LocalOperatorAssembler < GridType, TrialLocalFE, AnsatzLocalFE, T >::Element Element;
+    typedef typename LocalOperatorAssembler < GridType, TrialLocalFE, AnsatzLocalFE,T >::BoolMatrix BoolMatrix;
+    typedef typename LocalOperatorAssembler < GridType, TrialLocalFE, AnsatzLocalFE,T >::LocalMatrix LocalMatrix;
+
+
+    MooneyRivlinOperatorAssembler(const Dune::shared_ptr<GridFunction> displacement,
+                                    ctype a, ctype b, ctype c, ctype d, ctype e, int k) :
+        a_(a),b_(b),c_(c),d_(d),e_(e),k_(k), displacement_(displacement)
+    {}
+
+    MooneyRivlinOperatorAssembler(ctype a, ctype b, ctype c, ctype d, ctype e, int k) :
+        a_(a),b_(b),c_(c),d_(d),e_(e),k_(k)
+    {}
+
+    void indices(const Element& element, BoolMatrix& isNonZero, const TrialLocalFE& tFE, const AnsatzLocalFE& aFE) const
+    {
+        isNonZero = true;
+    }
+
+    void assemble(const Element& element, LocalMatrix& localMatrix, const TrialLocalFE& tFE, const AnsatzLocalFE& aFE) const
+    {
+        typedef Dune::FieldVector<ctype,dim> FVdim;
+        typedef Dune::FieldMatrix<ctype,dim,dim> FMdimdim;
+        typedef typename TrialLocalFE::Traits::LocalBasisType::Traits::JacobianType JacobianType;
+
+        int rows = tFE.localBasis().size();
+        int cols = aFE.localBasis().size();
+
+        localMatrix.setSize(rows,cols);
+        localMatrix = 0.0;
+
+        // get quadrature rule
+        const int order = (element.type().isSimplex()) 
+            ? 6*(tFE.localBasis().order())
+            : 6*(tFE.localBasis().order()*dim);
+        const Dune::template QuadratureRule<ctype, dim>& quad = QuadratureRuleCache<ctype, dim>::rule(element.type(), order, IsRefinedLocalFiniteElement<TrialLocalFE>::value(tFE) );
+
+        // store gradients of shape functions and base functions
+        std::vector<JacobianType> referenceGradients(tFE.localBasis().size());
+        std::vector<FVdim> gradients(tFE.localBasis().size());
+
+        // store entity geometry mapping
+        const typename Element::Geometry geometry = element.geometry();
+
+        // loop over quadrature points
+        for (size_t pt=0; pt < quad.size(); ++pt) {
+
+            // get quadrature point
+            const FVdim& quadPos = quad[pt].position();
+
+            // get transposed inverse of Jacobian of transformation
+            const FMdimdim& invJacobian = geometry.jacobianInverseTransposed(quadPos);
+
+            // get integration factor
+            const ctype integrationElement = geometry.integrationElement(quadPos);
+
+            // get gradients of shape functions
+            tFE.localBasis().evaluateJacobian(quadPos, referenceGradients);
+
+            // compute gradients of base functions
+            for (size_t i=0; i<gradients.size(); ++i) {
+
+                // transform gradients
+                gradients[i] = 0.0;
+                invJacobian.umv(referenceGradients[i][0], gradients[i]);
+
+            }
+
+            // evaluate displacement gradients the configuration at the quadrature point
+            typename GridFunction::DerivativeType localConfGrad;
+            if (displacement_->isDefinedOn(element))
+                displacement_->evaluateDerivativeLocal(element, quadPos, localConfGrad);
+            else
+                displacement_->evaluateDerivative(geometry.global(quadPos),localConfGrad);
+
+            // compute linearization of the determinante of the deformation gradient
+            FMdimdim linDefDet;
+            Dune::Elasticity::linearisedDefDet(localConfGrad,linDefDet);
+
+            // compute deformation gradient of the configuration
+            for (int i=0;i<dim;i++)
+                localConfGrad[i][i] += 1;                
+
+            // evaluate the derminante of the deformation gradient
+            const ctype J = localConfGrad.determinant();
+
+            //compute linearized strain for all shapefunctions
+            std::vector<Dune::array<FMdimdim, dim > > defGrad(rows);
+            std::vector<Dune::array<SymmetricTensor<dim>,dim> > linStrain(rows);
+
+            for (int i=0; i<rows; i++)
+                for (int k=0; k<dim; k++) {
+
+                    // The deformation gradient of the shape function
+                    defGrad[i][k] = 0;
+                    defGrad[i][k][k] = gradients[i];
+
+                    // The linearized strain is the symmetric product of defGrad and (Id+localConf)
+                    linStrain[i][k]=symmetricProduct(defGrad[i][k],localConfGrad);
+                }
+
+            // the trace of the strain tensor
+            SymmetricTensor<dim> strain;
+            computeStrain(localConfGrad,strain);
+            ctype trE(0);
+            for (int i=0; i<dim; i++)
+                trE += strain(i,i);                
+
+            // /////////////////////////////////////////////////
+            //   Assemble matrix
+            // /////////////////////////////////////////////////
+            ctype z = quad[pt].weight() * integrationElement;
+
+            ctype coeff= (2*d_*J-k_*e_*std::pow(J,-k_-1))*z;
+            ctype coeff2 = (2*d_+k_*(k_+1)*e_*std::pow(J,-k_-2))*z;
+
+            for (int row=0; row<rows; row++)
+                    for (int col=0; col<cols; col++) {
+                        
+                        // second derivative of the trace terms 
+                        StaticMatrix::addToDiagonal(localMatrix[row][col],z*(a_ + 2*b_*trE)*((gradients[row]*gradients[col])));
+                        // second derivative of the determinat term 
+                        localMatrix[row][col].axpy(coeff,hessDefDet(localConfGrad,gradients[row],gradients[col]));
+                      
+                        FVdim rowTerm(0),colTerm(0);
+                        linDefDet.mv(gradients[row],rowTerm); 
+                        linDefDet.mv(gradients[col],colTerm); 
+
+                        FVdim traceRowTerm(0),traceColTerm(0);
+                        localConfGrad.mv(gradients[row],traceRowTerm); 
+                        localConfGrad.mv(gradients[col],traceColTerm); 
+
+                        // second derivative of the scalar product term
+                        for (int rcomp=0; rcomp<dim; rcomp++) {
+                            // derivative of the coefficient term in the compressibility function
+                            localMatrix[row][col][rcomp].axpy(coeff2*rowTerm[rcomp],colTerm);
+
+                            // second derivative of the trace terms part 2
+                            localMatrix[row][col][rcomp].axpy(2*z*b_*traceRowTerm[rcomp],traceColTerm);
+
+                            // derivative of the scalar product terms
+                            for (int ccomp=0; ccomp<dim; ccomp++) {
+                                localMatrix[row][col][rcomp][ccomp] += (linStrain[row][rcomp]*linStrain[col][ccomp])*2*c_*z;
+                                localMatrix[row][col][rcomp][ccomp] += (strain*symmetricProduct(defGrad[row][rcomp],defGrad[col][ccomp]))*z*2*c_;
+                            } 
+                        }
+
+            }
+        }
+
+    }
+
+    /** \brief Set new configuration to be assembled at. Needed e.g. during Newton iteration.*/
+    void setConfiguration(Dune::shared_ptr<GridFunction> newConfig) {
+        displacement_ = newConfig;
+    }
+
+private:
+
+    //! Material parameters
+    ctype a_; ctype b_; ctype c_; ctype d_; ctype e_;
+    int k_;
+
+    /** \brief The configuration at which the linearized operator is evaluated.*/
+    Dune::shared_ptr<GridFunction> displacement_;
+
+
+    /** \brief Compute the matrix entry that arises from the second derivative of the determinant. 
+     *  Note that it is assumed the defGrad is the deformation gradient and not the displacement gradient   
+     */
+    Dune::FieldMatrix<ctype,dim,dim> hessDefDet(const Dune::FieldMatrix<ctype,dim,dim>& defGrad,
+                            const Dune::FieldVector<ctype,dim>& testGrad, const Dune::FieldVector<ctype,dim>& ansatzGrad) const
+    {
+        Dune::FieldMatrix<ctype, dim, dim> res(0);
+        
+        if (dim==2) {
+            res[0][1] = testGrad[0]*ansatzGrad[1]-ansatzGrad[0]*testGrad[1];
+            res[1][0] = -res[0][1];
+            return res;
+        }
+
+        //compute the cross product
+        Dune::FieldVector<ctype,dim> cross(0);    
+        for (int k=0; k<dim; k++)
+            cross[k]=testGrad[(k+1)%dim]*ansatzGrad[(k+2)%dim] -testGrad[(k+2)%dim]*ansatzGrad[(k+1)%dim];
+
+        // the resulting matrix is skew symmetric with entries <cross,degGrad[i]>
+        for (int i=0; i<dim; i++)
+            for (int j=i+1; j<dim; j++) {
+                int k= (dim-(i+j))%dim;
+            res[i][j] = (cross*defGrad[k])*std::pow(-1,k);
+            res[j][i] = -res[i][j];
+        }
+        
+        return res;
+    }
+
+    //! 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];
+            }
+        }
+    }
+
+    /** \brief Compute the symmetric product of two matrices, i.e.
+     *      0.5*(A^T * B + B^T * A )
+     */
+    SymmetricTensor<dim> symmetricProduct(const Dune::FieldMatrix<ctype, dim, dim>& a, const Dune::FieldMatrix<ctype, dim, dim>& b) const {
+
+        SymmetricTensor<dim> result(0.0);
+        for (int i=0;i<dim; i++)
+            for (int j=i;j<dim;j++)
+                for (int k=0;k<dim;k++)
+                    result(i,j)+= a[k][i]*b[k][j]+b[k][i]*a[k][j];   
+        result *= 0.5;                
+        
+        return result;
+    }
+
+};
+
+
+#endif
+

dune/elasticity/materials/mooneyrivlinmaterial.hh

+#ifndef MOONEY_RIVLIN_MATERIAL
+#define MOONEY_RIVLIN_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/mooneyrivlinfunctionalassembler.hh>
+#include <dune/elasticity/assemblers/mooneyrivlinoperatorassembler.hh>
+
+#include <dune/elasticity/materials/material.hh>
+#warning Mooney-Rivlin Might not work properly yet
+
+/** \brief Class representing a hyperelastic Mooney-Rivlin 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 :-( )
+ *
+ *    \f[
+ *      W(u)= a tr(E(u)) + b tr(E(u))^2 + c ||E(u)||^2 + gamma(J)
+ *    \f]
+ *
+ *  where
+ *      - \f$E\f$: the nonlinear strain tensor
+ *      - \f$tr \f$: the trace operator
+ *      - \f$J\f$ the determinant of the deformation gradient
+ *      - \f$a\f$,\f$b\f$,\f$c$\f material parameters
+ *      - \f$\gamma$\f A 'penalty' function which enforces orientation preservation
+ *
+ *  with
+ *      \f$\gamma(J)$:= dJ^2 + e J^-k,  k>=2\f
+ *      \f$k\f$ controlls orientation preservation (k \approx \floor(1/(1-2nu) -1))
+ */
+template <class Basis>
+class MooneyRivlinMaterial : public Material<Basis>
+{
+public:
+    typedef Material<Basis> Base;
+    typedef typename Basis::GridView::Grid GridType;
+    using typename Base::GlobalBasis;
+    using typename Base::Lfe;
+    using typename Base::LocalLinearization;
+    using typename Base::LocalHessian;
+    using typename Base::VectorType;
+
+    typedef MooneyRivlinFunctionalAssembler<GridType,Lfe> MonRivLinearization;
+    typedef MooneyRivlinOperatorAssembler<GridType, Lfe, Lfe> MonRivHessian;
+
+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:
+    MooneyRivlinMaterial(int k=3) :
+        E_(100), nu_(0.3), k_(k)
+    {
+        setupCoefficients();
+    }
+
+    template <BasisT>
+    MooneyRivlinMaterial(BasisT&& basis, ctype E, ctype nu, int k=3) :
+        Base(std::forward<BasisT>(basis)), E_(E), nu_(nu), k_(k)
+    {
+
+        setupCoefficients();
+
+        localLinearization_ = Dune::shared_ptr<MonRivLinearization>(new MonRivLinearization(a_, b_,c_ ,d_ ,e_,k_));
+        localHessian_ = Dune::shared_ptr<MonRivHessian>(new MonRivHessian(a_ , b_, c_, d_, e_,k_));
+    }
+
+    template <BasisT>
+    void setup(BasisT&& basis, ctype E, ctype nu)
+    {
+        E_ = E; nu_ = nu;
+
+        setupCoefficients();
+
+        if (localLinearization_)
+            localLinearization_.reset();
+        if (localHessian_)
+            localHessian_.reset();
+
+        localLinearization_ = std::shared_ptr<MonRivLinearization>(new MonRivLinearization(a_,b_,c_,d_,e_,k_));
+        localHessian_ = std::shared_ptr<MonRivHessian>(new MonRivHessian(a_,b_,c_,d_,e_,k_));
+
+        this->setBasis(std::forward<BasisT>(basis));
+    }
+
+    void getMaterialParameters(ctype& E, ctype& nu) {
+        E = E_; nu = nu_;
+    }
+
+    //! Evaluate the strain energy
+    ctype energy(const VectorType& coeff)
+    {
+        // make grid function
+        BasisGridFunction<Basis,VectorType> configuration(*this->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()) ? 5 : 5*dim;
+
+            const Dune::template QuadratureRule<ctype, dim>& quad = QuadratureRuleCache<ctype, dim>::rule(eIt->type(),order,0);
+
+            // 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,VectorType>::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;
+
+                // W(u)= a tr(E(u)) + b tr(E(u))^2 + c ||E(u)||^2 + gamma(J)
+                energy += z*(a_*trE + b_*trE*trE + c_*(strain*strain) + d_*J*J + e_*std::pow(J,-k_));
+            }
+        }
+
+        return 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_;}
+
+private:
+
+    //! Compute coefficients s.t. polyconvexity holds
+    void setupCoefficients() {
+
+        ctype lambda = E_*nu_ / ((1+nu_)*(1-2*nu_));
+        ctype mu     = E_ / (2*(1+nu_));
+
+        // Choose coefficients in such a way that we have polyconvexity
+        // First/Second derivative of the compressible function part of gamma (=-kJ^(-k-1)) at 1.0
+        ctype ld1 = -k_;
+        ctype ld2 = k_*(k_+1);
+        if (ld1 >=0 || ld2 <= 0)
+            std::cout<<"Coerciveness failed\n";
+
+        // Check if lame constants admit coerciveness
+        ctype rho= -ld1/(ld2-ld1);
+
+        if( ( (rho < 0.5 && lambda < (1.0/rho-2.0)*mu) || lambda <= 0 || mu <= 0 ))
+            std::cout<<"Coerciveness failed\n";
+
+        //const ctype somePositiveConstant = - mu_ + rho*(lambda_+2.*mu_);
+        //if(somePositiveConstant <= 0)
+
+        e_ = (lambda+2.0*mu)/(ld2-ld1);
+        if(e_ <= 0 )
+            std::cout<<"Coerciveness failed\n";
+
+        // parameters a, and b are used if expressed in terms of the Green-StVenant strain tensor
+        // the choice of b and d is admissible but heuristic -> @TODO justified choice
+        d_ = (0.5*rho-0.25)*mu+0.25*rho*lambda;
+        if(d_ > 0.25*mu)
+            d_ = (rho-0.75)*mu+0.5*rho*lambda;
+
+        ctype b_ = -mu + rho*(lambda+2.*mu)-2.*d_;
+        ctype c_ = -b_;
+        ctype a_ = b_ + mu;
+
+        // last check if I didn't miss a condition
+        ctype alpha = 0.5*(mu-b_);
+        if(alpha <= 0 || b_ <= 0 || d_ <= 0)
+            std::cout<<"Coerciveness failed\n";
+    }
+
+    //! First derivative of the strain energy
+    std::shared_ptr<MonRivLinearization> localLinearization_;
+
+    //! Second derivative of the strain energy
+    std::shared_ptr<MonRivHessian> localHessian_;
+
+    //! Elasticity modulus
+    ctype E_;
+    //! Shear modulus
+    ctype nu_;
+    //! Exponent of the compressibility function (>= 2)
+    int k_;
+    //! Material parameters
+    ctype a_; ctype b_; ctype c_; ctype d_; ctype e_;
+
+    //! 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