diff --git a/dune/elasticity/materials/adolcneohookeanmaterial.hh b/dune/elasticity/materials/adolcneohookeanmaterial.hh deleted file mode 100644 index 2a6f6d4f6e10110110b7ecb3ff7a25fd99a2b8f6..0000000000000000000000000000000000000000 --- a/dune/elasticity/materials/adolcneohookeanmaterial.hh +++ /dev/null @@ -1,176 +0,0 @@ -// -*- 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> -#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 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() = default; - - LocalNeoHookeanEnergy(field_type E, field_type nu) - { - lambda_ = E*nu / ((1+nu)*(1-2*nu)); - mu_ = E / (2*(1+nu)); - } - - LocalNeoHookeanEnergy& operator=(LocalNeoHookeanEnergy&& other) { - lambda_ = other.lambda_; - mu_ = other.mu_; - } - - 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::strain(localGradient,strain); - - // turn displacement gradient into deformation gradient - Dune::MatrixVector::addToDiagonal(localGradient, 1.0); - - // 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() = default; - - template <class BasisT> - AdolcNeoHookeanMaterial(BasisT&& basis, ReturnType E, ReturnType nu, bool vectorMode=true) : - localEnergy_(E,nu) - { - Base::setup(std::forward<BasisT>(basis), localEnergy_,vectorMode); - } - - - template <class BasisT> - void setup(BasisT&& basis, ReturnType E, ReturnType nu, bool vectorMode=true) { - localEnergy_ = LocalEnergy(E,nu); - Base::setup(std::forward<BasisT>(basis), localEnergy_,vectorMode); - } - - using Base::energy; - private: - LocalEnergy localEnergy_; -}; - -#endif