The NeoHookean material using Adol-C. This can be used as an example

dune/elasticity/materials/adolcneohookeanmaterial.hh

+// -*- tab-width: 8; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set ts=8 sw=4 et sts=4:
+#ifndef DUNE_ELASTICITY_MATERIALS_ADOLC_NEO_HOOKEAN_MATERIAL_HH
+#define DUNE_ELASTICITY_MATERIALS_ADOLC_NEO_HOOKEAN_MATERIAL_HH
+
+#include <dune/fufem/assemblers/localassemblers/adolclocalenergy.hh>
+#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/materials/adolcmaterial.hh>
+#include <dune/elasticity/common/elasticityhelpers.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 :-( )
+ *
+ * *  Laursen:
+ *    \f[
+ *      W(u)= \frac{\lambda}4(J^2-1)-(\frac{\lambda}2+\mu)log(J)+\mu tr(E(u))
+ *    \f]
+ *
+ *  Fischer/Wriggers:
+ *    \f[
+ *      W(u)= \frac{\lambda}2(J-1)^2 - 2 \mu*log(J)+\mu tr(E(u))
+ *    \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$\lambda\f$,\f$\mu\f$ material parameters (first lame and shear modulus)
+ */
+template <class GridType, class LocalFiniteElement>
+class LocalNeoHookeanEnergy : public Adolc::LocalEnergy<GridType, LocalFiniteElement,GridType::dimension>
+{
+  public:
+
+    enum {dim = GridType::dimension};
+    typedef typename GridType::ctype ctype;
+    typedef typename GridType::template Codim<0>::Entity Element;
+
+    typedef typename LocalFiniteElement::Traits::LocalBasisType::Traits::RangeFieldType field_type;
+    typedef typename LocalFiniteElement::Traits::LocalBasisType::Traits::JacobianType JacobianType;
+
+    typedef Adolc::LocalEnergy<GridType, LocalFiniteElement, dim> Base;
+    typedef typename Base::CoefficientVectorType CoefficientVectorType;
+    typedef typename Base::ReturnType ReturnType;
+
+    LocalNeoHookeanEnergy(field_type E, field_type nu)
+    {
+        lambda_ = E*nu / ((1+nu)*(1-2*nu));
+        mu_     = E / (2*(1+nu));
+    }
+
+    ReturnType energy(const Element& element, const LocalFiniteElement& lfe,
+            const CoefficientVectorType& localCoeff) const
+    {
+        ReturnType energy=0;
+
+        // TODO Get proper quadrature rule
+        // get quadrature rule
+        const int order = (element.type().isSimplex()) ? 4 : 4*dim;
+
+        const auto& quad = QuadratureRuleCache<ctype, dim>::rule(element.type(),order,
+                 IsRefinedLocalFiniteElement<LocalFiniteElement>::value(lfe));
+
+        // the element geometry mapping
+        const typename Element::Geometry geometry = element.geometry();
+
+        // store gradients of shape functions and base functions
+        std::vector<JacobianType> referenceGradients(lfe.localBasis().size());
+
+        // loop over quadrature points
+        for (size_t pt=0; pt < quad.size(); ++pt) {
+
+            // get quadrature point
+            const auto& quadPos = quad[pt].position();
+
+            // get transposed inverse of Jacobian of transformation
+            const auto& invJacobian = geometry.jacobianInverseTransposed(quadPos);
+
+            // get integration factor
+            const ctype integrationElement = geometry.integrationElement(quadPos);
+
+            // get gradients of shape functions
+            lfe.localBasis().evaluateJacobian(quadPos, referenceGradients);
+
+            // compute gradient of the configuration
+            Dune::FieldMatrix<ReturnType,dim,dim> localGradient(0);
+            for (size_t k=0; k<referenceGradients.size(); ++k) {
+                Dune::FieldVector<ReturnType,dim> gradient(0);
+                invJacobian.umv(referenceGradients[k][0], gradient);
+                for (int i=0; i<dim; ++i)
+                    for (int j=0; j<dim; ++j)
+                        localGradient[i][j] += localCoeff[k][i] * gradient[j];
+            }
+
+            SymmetricTensor<dim, ReturnType> strain(0.0);
+            Dune::Elasticity::computeNonlinearStrain(localGradient,strain);
+
+            // turn displacement gradient into deformation gradient
+            for (int i=0;i<dim;i++)
+                localGradient[i][i] += 1;
+
+            // evaluate the derminante of the deformation gradient
+            const ReturnType J = localGradient.determinant();
+
+            ReturnType z = quad[pt].weight()*integrationElement;
+
+#ifdef LAURSEN
+            energy += z*(0.25*lambda_*(J*J-1)-(lambda_*0.5+mu_)*std::log(J)+mu_*strain.trace());
+#else
+            energy += z*(0.5*lambda_*(J-1)*(J-1)-2*mu_*std::log(J)+mu_*strain.trace());
+#endif
+        }
+
+        return energy;
+    }
+
+  private:
+    //! First Lame constant
+    field_type lambda_;
+    //! Second Lame constant
+    field_type mu_;
+
+};
+
+template <class Basis>
+class AdolcNeoHookeanMaterial : public AdolcMaterial<Basis>
+{
+  public:
+    typedef AdolcMaterial<Basis> Base;
+    typedef typename Base::GridType GridType;
+    typedef typename Base::GlobalBasis GlobalBasis;
+    typedef typename Base::Lfe Lfe;
+    typedef typename Base::LocalLinearization LocalLinearization;
+    typedef typename Base::LocalHessian LocalHessian;
+    typedef typename Base::VectorType VectorType;
+    typedef typename Base::GridFunction GridFunction;
+    typedef typename Base::ReturnType ReturnType;
+
+    typedef LocalNeoHookeanEnergy<GridType,Lfe> LocalEnergy;
+    using Base::dim;
+
+    AdolcNeoHookeanMaterial(const Basis& basis, ReturnType E, ReturnType nu, bool vectorMode=true) :
+        localEnergy_(E,nu)
+    {
+        this->setup(basis, localEnergy_,vectorMode);
+    }
+
+    using Base::energy;
+  private:
+    LocalEnergy localEnergy_;
+};
+
+#endif