Make the material also act like a local energy so we can use Adol-C

This way we dont need a seperate class for that

dune/elasticity/materials/mooneyrivlinmaterial.hh

 #ifndef MOONEY_RIVLIN_MATERIAL
 #define MOONEY_RIVLIN_MATERIAL
 
+#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>
 
@@ -36,7 +36,10 @@
  *      \f$k\f$ controlls orientation preservation (k \approx \floor(1/(1-2nu) -1))
  */
 template <class Basis>
-class MooneyRivlinMaterial : public Material<Basis>
+class MooneyRivlinMaterial : public Material<Basis>,
+            public Adolc::LocalEnergy<typename Material<Basis>::GridType,
+                                      typename Material<Basis>::Lfe,
+                                      Material<Basis>::GridType::dimension>
 {
 public:
     typedef Material<Basis> Base;
@@ -49,6 +52,11 @@ public:
     using typename Base::ReturnType;
     using typename Base::GridFunction;
 
+    using AdolCBase = Adolc::LocalEnergy<GridType, Lfe, dim>;
+    using Element = typename GridType::template Codim<0>::Entity;
+    using AdolcCoefficients = typename AdolCBase::CoefficientVectorType;
+    using AdolcEnergy = typename AdolCBase::ReturnType;
+    using AdolcFieldType = typename AdolCBase::ReturnType;
     using MonRivLinearization = MooneyRivlinFunctionalAssembler<GridType, Lfe>;
     using MonRivHessian = MooneyRivlinOperatorAssembler<GridType, Lfe, Lfe>;
 
@@ -134,6 +142,61 @@ public:
         return energy;
     }
 
+    //! Compute local energy, required for the Adol-C interface
+    AdolcEnergy energy(const Element& element, const Lfe& lfe,
+            const AdolcCoefficients& localCoeff) const
+    {
+        AdolcEnergy energy=0;
+
+        // TODO Get proper quadrature rule
+        // get quadrature rule (should depend on k?)
+        const int order = (element.type().isSimplex()) ? 5 : 5*dim;
+        const auto& quad = QuadratureRuleCache<typename GridType::ctype, dim>::rule(element.type(), order, 0);
+
+        const auto& geometry = element.geometry();
+
+        using LfeTraits = typename Lfe::Traits::LocalBasisType::Traits;
+        using Jacobian = typename LfeTraits::JacobianType;
+        std::vector<Jacobian> referenceGradients(lfe.localBasis().size());
+
+        for (const auto& pt : quad) {
+
+            const auto& quadPos = pt.position();
+            auto integrationElement = geometry.integrationElement(quadPos);
+            const auto& invJacobian = geometry.jacobianInverseTransposed(quadPos);
+
+
+            // get gradients of shape functions
+            lfe.localBasis().evaluateJacobian(quadPos, referenceGradients);
+
+            // compute gradient of the configuration
+            Dune::FieldMatrix<AdolcFieldType, dim, dim> localGradient(0);
+            for (size_t k=0; k < referenceGradients.size(); ++k) {
+
+              Dune::FieldVector<AdolcFieldType, 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];
+            }
+
+            auto strain = Dune::Elasticity::strain(localGradient);
+            auto trE = strain.trace();
+
+            // make deformation gradient out of the discplacement
+            Dune::MatrixVector::addToDiagonal(localGradient, 1.0);
+
+            auto J = localGradient.determinant();
+
+            auto z = 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(std::shared_ptr<GridFunction> displace) {
         localLinearization_->setConfiguration(displace);

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>
@@ -38,7 +39,10 @@
  *      - \f$\lambda\f$,\f$\mu\f$ material parameters (first lame and shear modulus)
  */
 template <class Basis>
-class NeoHookeanMaterial : public Material<Basis>
+class NeoHookeanMaterial : public Material<Basis>,
+                           public Adolc::LocalEnergy<typename Material<Basis>::GridType,
+                                      typename Material<Basis>::Lfe,
+                                      Material<Basis>::GridType::dimension>
 {
 public:
     typedef Material<Basis> Base;
@@ -51,6 +55,11 @@ public:
     typedef typename Base::GridFunction GridFunction;
     typedef typename Base::ReturnType ReturnType;
 
+    using AdolCBase = Adolc::LocalEnergy<GridType, Lfe, dim>;
+    using Element = typename GridType::template Codim<0>::Entity;
+    using AdolcCoefficients = typename AdolCBase::CoefficientVectorType;
+    using AdolcEnergy = typename AdolCBase::ReturnType;
+    using AdolcFieldType = typename AdolCBase::ReturnType;
 
 private:
     using Base::dim;
@@ -152,56 +161,60 @@ public:
         return energy;
     }
 
-    //! Evaluate the strain energy
-    template <class Element>
-    ReturnType localEnergy(const GridFunction& displace, const Element& e) const
+    AdolcEnergy energy(const Element& element, const Lfe& lfe,
+            const AdolcCoefficients& localCoeff) const
     {
-        ReturnType energy=0;
+        AdolcEnergy energy=0;
 
-            // TODO Get proper quadrature rule
-            // get quadrature rule
-            const int order = (e.type().isSimplex()) ? 5 : 4*dim;
+        // TODO Get proper quadrature rule
+        // get quadrature rule
+        const int order = (element.type().isSimplex()) ? 5 : 4*dim;
 
-            const auto& quad = QuadratureRuleCache<ctype, dim>::rule(e.type(), order, 0);
+        const auto& quad = QuadratureRuleCache<typename GridType::ctype, dim>::rule(element.type(), order,
+                 IsRefinedLocalFiniteElement<Lfe>::value(lfe));
 
-            const auto& geometry = e.geometry();
+        const auto& geometry = element.geometry();
 
-            // loop over quadrature points
-            for (const auto& pt : quad) {
+        using LfeTraits = typename Lfe::Traits::LocalBasisType::Traits;
+        using Jacobian = typename LfeTraits::JacobianType;
+        std::vector<Jacobian> referenceGradients(lfe.localBasis().size());
 
-                // get quadrature point
-                const auto& quadPos = pt.position();
+        for (const auto& pt : quad) {
 
-                // get integration factor
-                const auto integrationElement = geometry.integrationElement(quadPos);
+            const auto& quadPos = pt.position();
 
-                // evaluate displacement gradient at the quadrature point
-                typename BasisGridFunction<Basis,VectorType>::DerivativeType localDispGrad;
+            const auto& invJacobian = geometry.jacobianInverseTransposed(quadPos);
+            const auto integrationElement = geometry.integrationElement(quadPos);
 
-                if (displace.isDefinedOn(e))
-                    displace.evaluateDerivativeLocal(e, quadPos, localDispGrad);
-                else
-                    displace.evaluateDerivative(geometry.global(quadPos),localDispGrad);
+            // get gradients of shape functions
+            lfe.localBasis().evaluateJacobian(quadPos, referenceGradients);
 
-                SymmetricTensor<dim, ReturnType> strain;
-                Dune::Elasticity::strain(localDispGrad, strain);
+            // compute gradient of the configuration
+            Dune::FieldMatrix<AdolcFieldType, dim, dim> localGradient(0);
+            for (size_t k=0; k < referenceGradients.size(); ++k) {
 
-                // the trace
-                auto trE = strain.trace();
+                Dune::FieldVector<AdolcFieldType, dim> gradient(0);
+                invJacobian.umv(referenceGradients[k][0], gradient);
 
-                // turn displacement gradient into deformation gradient
-                Dune::MatrixVector::addToDiagonal(localDispGrad, 1.0);
+                for (int i=0; i < dim; ++i)
+                    for (int j=0; j < dim; ++j)
+                        localGradient[i][j] += localCoeff[k][i] * gradient[j];
+            }
 
-                // evaluate the derminante of the deformation gradient
-                ReturnType J = localDispGrad.determinant();
-                ctype z = pt.weight()*integrationElement;
+            auto strain = Dune::Elasticity::strain(localGradient);
+
+            // turn displacement gradient into deformation gradient
+            Dune::MatrixVector::addToDiagonal(localGradient, 1.0);
+
+            AdolcFieldType J = localGradient.determinant();
+            AdolcFieldType z = pt.weight()*integrationElement;
 
 #ifdef LAURSEN
-                energy += z*(0.25*lambda_*(J*J-1)-(lambda_*0.5+mu_)*std::log(J)+mu_*trE);
+            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_*trE);
+            energy += z*(0.5 * lambda_ * (J-1)*(J-1) - 2 * mu_ * std::log(J) + mu_ * strain.trace());
 #endif
-            }
+        }
         return energy;
     }