diff --git a/dune/elasticity/materials/adolcneohookeanmaterial.hh b/dune/elasticity/materials/adolcneohookeanmaterial.hh new file mode 100644 index 0000000000000000000000000000000000000000..29b5b295b83521396c0b2da292f367cf9277e04c --- /dev/null +++ b/dune/elasticity/materials/adolcneohookeanmaterial.hh @@ -0,0 +1,160 @@ +// -*- 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