diff --git a/dune/elasticity/materials/mooneyrivlinmaterial.hh b/dune/elasticity/materials/mooneyrivlinmaterial.hh
index 7450b8cd767b2d6a00bdcb9c0913f5d533a0d9c6..675e9150bdbd861f58f03aacb6d15c73c4d2d0c7 100644
--- a/dune/elasticity/materials/mooneyrivlinmaterial.hh
+++ b/dune/elasticity/materials/mooneyrivlinmaterial.hh
@@ -1,9 +1,9 @@
#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);
diff --git a/dune/elasticity/materials/neohookeanmaterial.hh b/dune/elasticity/materials/neohookeanmaterial.hh
index a7c0204eeca00b52342f8bd8b554b1a5b7d6d6be..8b090b1cff587d963393d6d93f183c17f611e84a 100644
--- a/dune/elasticity/materials/neohookeanmaterial.hh
+++ b/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;
}