neohookefunctionalassembler.hh

#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