Add a material class that use Adol-C to compute

the linearization and hessian.

dune/elasticity/materials/neohookeanmaterial.hh

@@ -3,6 +3,7 @@
 #ifndef DUNE_ELASTICITY_MATERIALS_NEO_HOOKEAN_MATERIAL_HH
 #define DUNE_ELASTICITY_MATERIALS_NEO_HOOKEAN_MATERIAL_HH
 
+#include <dune/fufem/assemblers/localassemblers/adolclocalenergy.hh>
 #include <dune/fufem/quadraturerules/quadraturerulecache.hh>
 
 #include <dune/fufem/symmetrictensor.hh>
@@ -174,4 +175,129 @@ private:
     ReturnType mu_;
 };
 
+
+//! Local energy for a geometric exact St. Venant--Kirchhoff material
+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 NeoHookeMaterial : 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;
+
+    NeoHookeMaterial(const Basis& basis, ReturnType E, ReturnType nu) :
+        localEnergy_(E,nu)
+    {
+        this->setup(basis, localEnergy_);
+    }
+
+    using Base::energy;
+  private:
+    LocalEnergy localEnergy_;
+};
+
 #endif