Select Git revision
WhiteHand.mat
Code owners
Assign users and groups as approvers for specific file changes. Learn more.
neohookeanmaterial.hh 9.06 KiB
// -*- tab-width: 8; indent-tabs-mode: nil; c-basic-offset: 4 -*-
// vi: set ts=8 sw=4 et sts=4:
#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>
#include <dune/fufem/functions/virtualgridfunction.hh>
#include <dune/fufem/functions/basisgridfunction.hh>
#include <dune/fufem/quadraturerules/quadraturerulecache.hh>
#include <dune/elasticity/materials/material.hh>
#include <dune/elasticity/assemblers/neohookefunctionalassembler.hh>
#include <dune/elasticity/assemblers/neohookeoperatorassembler.hh>
#include <dune/matrix-vector/addtodiagonal.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 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;
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;
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;
typedef typename GridType::ctype ctype;
typedef NeoHookeFunctionalAssembler<GridType,Lfe> NeoLinearization;
typedef NeoHookeOperatorAssembler<GridType, Lfe, Lfe> NeoHessian;
typedef typename GridType::template Codim<0>::Geometry::LocalCoordinate LocalCoordinate;
typedef typename GridType::template Codim<0>::LeafIterator ElementIterator;
public:
NeoHookeanMaterial() :
lambda_(1.0),
mu_(0.3)
{}
template <class BasisT>
NeoHookeanMaterial(BasisT&& basis, ReturnType E, ReturnType nu) :
Base(std::forward<BasisT>(basis))
{
lambda_ = E*nu / ((1+nu)*(1-2*nu));
mu_ = E / (2*(1+nu));
localLinearization_ = std::make_shared<NeoLinearization>(lambda_,mu_);
localHessian_ = std::make_shared<NeoHessian>(lambda_,mu_);
}
template <class BasisT>
void setup(BasisT&& basis, ReturnType E, ReturnType nu)
{
lambda_ = E*nu / ((1+nu)*(1-2*nu));
mu_ = E / (2*(1+nu));
localLinearization_ = std::make_shared<NeoLinearization>(lambda_,mu_);
localHessian_ = std::make_shared<NeoHessian>(lambda_,mu_);
this->setBasis(std::forward<BasisT>(basis));
}
//! Evaluate the strain energy
ReturnType energy(std::shared_ptr<GridFunction> displace) const
{
ReturnType energy=0;
const auto& leafView = this->basis().getGridView().grid().leafGridView();
for (const auto& e : elements(leafView)) {
// get quadrature rule
QuadratureRuleKey feKey(this->basis_->getLocalFiniteElement(e));
// the determinant involves derivative^dim
// further use 2nd order to approximate the logarithm
//feKey.setOrder(std::pow(feKey.derivative().order(),2*dim));
const int order = (e.type().isSimplex()) ? 6 : 4*dim;
const auto& quad = QuadratureRuleCache<ctype, dim>::rule(e.type(), order, 0);
//const auto& quad = QuadratureRuleCache<ctype, dim>::rule(feKey);
const auto& geometry = e.geometry();
// loop over quadrature points
for (const auto& pt : quad) {
// get quadrature point
const auto& quadPos = pt.position();
// get integration factor
const ctype integrationElement = geometry.integrationElement(quadPos);
// evaluate displacement gradient at the quadrature point
typename BasisGridFunction<Basis,VectorType>::DerivativeType localDispGrad;
if (displace->isDefinedOn(e))
displace->evaluateDerivativeLocal(e, quadPos, localDispGrad);
else
displace->evaluateDerivative(geometry.global(quadPos), localDispGrad);
SymmetricTensor<dim,ReturnType> strain;
Dune::Elasticity::strain(localDispGrad,strain);
// the trace
auto trE = strain.trace();
// turn displacement gradient into deformation gradient
Dune::MatrixVector::addToDiagonal(localDispGrad, 1.0);
// evaluate the derminante of the deformation gradient
const ReturnType J = localDispGrad.determinant();
ctype z = pt.weight()*integrationElement;
#ifdef LAURSEN
energy += z*(0.25*lambda_*(J*J-1)-(lambda_*0.5+mu_)*std::log(J)+mu_*trE);
#else
energy += z*(0.5*lambda_*(J-1)*(J-1)-2*mu_*std::log(J)+mu_*trE);
#endif
}
}
return energy;
}
AdolcEnergy energy(const Element& element, const Lfe& lfe,
const AdolcCoefficients& localCoeff) const
{
AdolcEnergy energy=0;
// TODO Get proper quadrature rule
// get quadrature rule
const int order = (element.type().isSimplex()) ? 5 : 4*dim;
const auto& quad = QuadratureRuleCache<typename GridType::ctype, dim>::rule(element.type(), order,
IsRefinedLocalFiniteElement<Lfe>::value(lfe));
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();
const auto& invJacobian = geometry.jacobianInverseTransposed(quadPos);
const auto integrationElement = geometry.integrationElement(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);
// 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_ * strain.trace());
#else
energy += z*(0.5 * lambda_ * (J-1)*(J-1) - 2 * mu_ * std::log(J) + mu_ * strain.trace());
#endif
}
return energy;
}
//! Return the local assembler of the first derivative of the strain energy
LocalLinearization& firstDerivative(std::shared_ptr<GridFunction> displace) {
localLinearization_->setConfiguration(displace);
return *localLinearization_;
}
//! Return the local assembler of the second derivative of the strain energy
LocalHessian& secondDerivative(std::shared_ptr<GridFunction> displace) {
localHessian_->setConfiguration(displace);
return *localHessian_;
}
private:
//! First derivative of the strain energy
std::shared_ptr<NeoLinearization> localLinearization_;
//! Second derivative of the strain energy
std::shared_ptr<NeoHessian> localHessian_;
//! First Lame constant
ReturnType lambda_;
//! Second Lame constant
ReturnType mu_;
};
#endif