diff --git a/dune/elasticity/assemblers/neohookefunctionalassembler.hh b/dune/elasticity/assemblers/neohookefunctionalassembler.hh new file mode 100644 index 0000000000000000000000000000000000000000..2292b8deb083c772f08089aa155681c72b976619 --- /dev/null +++ b/dune/elasticity/assemblers/neohookefunctionalassembler.hh @@ -0,0 +1,182 @@ +#ifndef NEO_HOOKE_FUNCTIONAL_ASSEMBLER_HH +#define NEO_HOOKE_FUNCTIONAL_ASSEMBLER_HH + +#include <dune/common/fvector.hh> +#include <dune/istl/bvector.hh> + +#include <dune/fufem/quadraturerules/quadraturerulecache.hh> +#include "dune/fufem/assemblers/localfunctionalassembler.hh" +#include <dune/fufem/symmetrictensor.hh> +#include <dune/fufem/functions/virtualgridfunction.hh> + +#include <dune/elasticity/common/elasticityhelper.hh> + +/** \brief Local assembler for the linearization of a nonlinear Neo-Hookean energy at a displacement u + * + * 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] + * + * The linearization of this energy at \f$ u\f$ as a functional applied to a test function \f$ v\f$] reads + * + * + * Laursen: + * \f[ + * l_u(v) := (\frac{\lambda}2 J - \frac{\frac{\lambda}2 + \mu}{J})\frac{\partial J}{\partial u}v + + \mu\frac{\partial tr(E(u))}{\partial u} v + * \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 +*/ +template <class GridType> +class NeoHookeFunctionalAssembler : + public LocalFunctionalAssembler<GridType, Dune::FieldVector<typename GridType::ctype ,GridType::dimension> > + +{ + static const int dim = GridType::dimension; + typedef typename GridType::ctype ctype; + typedef typename Dune::FieldVector<ctype,GridType::dimension> T; + +public: + typedef typename LocalFunctionalAssembler<GridType,T>::ElementPointer ElementPointer; + typedef typename LocalFunctionalAssembler<GridType,T>::LocalVector LocalVector; + typedef VirtualGridFunction<GridType, Dune::FieldVector<ctype,GridType::dimension> > GridFunction; + + //! Create assembler from material parameters and a grid + NeoHookeFunctionalAssembler(ctype lambda, ctype mu, const Dune::shared_ptr<GridFunction> displacement) : + lambda_(lambda), + mu_(mu), + configuration_(configuration) + {} + + //! Create assembler from material parameters + NeoHookeFunctionalAssembler(ctype lambda, ctype mu) : + lambda_(lambda), + mu_(mu) + {} + + template <class TrialLocalFE> + void assemble(const ElementPointer& element, LocalVector& localVector, const TrialLocalFE& tFE) { + + typedef typename Dune::FieldVector<ctype,dim> FVdim; + typedef typename Dune::FieldMatrix<ctype,dim,dim> FMdimdim; + typedef typename TrialLocalFE::Traits::LocalBasisType::Traits::JacobianType JacobianType; + typedef typename GridType::template Codim<0>::Geometry::LocalCoordinate LocalCoordinate; + + int rows = tFE.localBasis().size(); + + localVector.resize(rows); + localVector = 0.0; + + // get quadrature rule + const int order = (element->type().isSimplex()) + ? 4*(tFE.localBasis().order()) + : 4*(tFE.localBasis().order()*dim); + + // get quadrature rule + const Dune::template QuadratureRule<ctype, dim>& quad = QuadratureRuleCache<ctype, dim>::rule(element->type(), order, IsRefinedLocalFiniteElement<TrialLocalFE>::value(tFE) ); + + // store gradients of shape functions and base functions + std::vector<JacobianType> referenceGradients(tFE.localBasis().size()); + std::vector<FVdim> gradients(tFE.localBasis().size()); + + // store geometry mapping of the entity + const typename ElementPointer::Entity::Geometry geometry = element->geometry(); + + // loop over quadrature points + for (size_t pt=0; pt < quad.size(); ++pt) { + + // get quadrature point + const LocalCoordinate& quadPos = quad[pt].position(); + + // get transposed inverse of Jacobian of transformation + const FMdimdim& invJacobian = geometry.jacobianInverseTransposed(quadPos); + + // get integration factor + const ctype integrationElement = geometry.integrationElement(quadPos); + + // get gradients of shape functions + tFE.localBasis().evaluateJacobian(quadPos, referenceGradients); + + // compute gradients of base functions + for (size_t i=0; i<gradients.size(); ++i) { + + // transform gradients + gradients[i] = 0.0; + invJacobian.umv(referenceGradients[i][0], gradients[i]); + + } + + // evaluate the displacement gradients of the configuration at the quadrature point + typename GridFunction::DerivativeType localConfGrad; + if (configuration_->isDefinedOn(*element)) + configuration_->evaluateDerivativeLocal(*element, quadPos, localConfGrad); + else + configuration_->evaluateDerivative(geometry.global(quadPos),localConfGrad); + + // evaluate the derminante of the deformation gradient + const ctype J = localConfGrad.determinant(); + + // compute linearization of the determinante of the deformation gradient + FMdimdim linDefDet; + Dune::Elasticity::linearisedDefDet(localConfGrad,linDefDet); + + + /* Compute the nonlinear strain tensor from the displacement gradient*/ + SymmetricTensor<dim> strain; + computeStrain(localConfGrad,strain); + + // make deformation gradient out of the discplacement + for (int i=0;i<dim;i++) + localConfGrad[i][i] += 1; + + // collect terms + FMdimdim fu(0); + +#ifdef LAURSEN + // Laursen + fu.axpy(J*lambda_/2-(lambda_/2+mu_)/J,linDefDet); +#else + // Fischer + fu.axpy((J-1)*lambda_-2*mu_/J,linDefDet); +#endif + // the linearization of the trace is just deformation gradient + fu.axpy(mu_,localConfGrad); + + ctype z = quad[pt].weight()*integrationElement; + + // add vector entries + for(int row=0; row<rows; row++) + fu.usmv(z,gradients[row],localVector[row]); + } + } + + /** \brief Set new configuration. In Newton iterations this needs to be assembled more than one time.*/ + void setConfiguration(Dune::shared_ptr<GridFunction> newConfig) { + configuration_ = newConfig; + } + +private: + /** \brief Lame constant */ + ctype lambda_; + + /** \brief Lame constant */ + ctype mu_; + + /** \brief The configuration at which the functional is evaluated.*/ + Dune::shared_ptr<GridFunction> configuration_; +}; + +#endif + +