diff --git a/dune/elasticity/assemblers/mooneyrivlinfunctionalassembler.hh b/dune/elasticity/assemblers/mooneyrivlinfunctionalassembler.hh
new file mode 100644
index 0000000000000000000000000000000000000000..d68ce37392ea8bc653fe71549783ea9a8fcd7e06
--- /dev/null
+++ b/dune/elasticity/assemblers/mooneyrivlinfunctionalassembler.hh
@@ -0,0 +1,215 @@
+#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
+
+
diff --git a/dune/elasticity/assemblers/mooneyrivlinoperatorassembler.hh b/dune/elasticity/assemblers/mooneyrivlinoperatorassembler.hh
new file mode 100644
index 0000000000000000000000000000000000000000..353ebc3078c8c874fb8982c9888bac748984a582
--- /dev/null
+++ b/dune/elasticity/assemblers/mooneyrivlinoperatorassembler.hh
@@ -0,0 +1,278 @@
+#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
+
diff --git a/dune/elasticity/materials/mooneyrivlinmaterial.hh b/dune/elasticity/materials/mooneyrivlinmaterial.hh
new file mode 100644
index 0000000000000000000000000000000000000000..57d2a465fc9dfd2b2a5f5d6b71f20ccad1f7633e
--- /dev/null
+++ b/dune/elasticity/materials/mooneyrivlinmaterial.hh
@@ -0,0 +1,246 @@
+#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