Skip to content
Snippets Groups Projects
Select Git revision
  • 4a335a72c2527832093e40430c61c7c6fa5bb313
  • master default protected
  • dune-tkr-article
  • patrizio-convexity-test
  • releases/2.6-1
  • releases/2.5-1
  • releases/2.4-1
  • releases/2.3-1
  • releases/2.2-1
  • releases/2.1-1
  • releases/2.0-1
  • dune-tkr-article-base
  • dune-tkr-article-patched
  • subversion->git
14 results

neohookeanmaterial.hh

Blame
  • Forked from agnumpde / dune-elasticity
    Source project has a limited visibility.
    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