geomexactstvenantkirchhoffmaterial.hh

#ifndef GEOMETRICALLY_EXACT_ST_VENANT_KIRCHHOFF_MATERIAL
#define GEOMETRICALLY_EXACT_ST_VENANT_KIRCHHOFF_MATERIAL

#include <dune/fufem/functiontools/quadraturerulecache.hh>

#include <dune/fufem/symmetrictensor.hh>
#include <dune/fufem/functions/virtualgridfunction.hh>
#include <dune/fufem/functions/basisgridfunction.hh>

#include <dune/fufem/assemblers/localassemblers/geononlinstvenantfunctionalassembler.hh>
#include <dune/fufem/assemblers/localassemblers/geononlinlinearizedstvenantassembler.hh>

/* \brief Class representing a hyperelastic geometrically exact St.Venant Kirchhoff 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 :-( )
*/
template <class Basis>
class GeomExactStVenantMaterial
{
public:
    typedef typename Basis::GridView::Grid GridType;
    typedef typename Basis::LocalFiniteElement Lfe;    

    typedef GeomNonlinElasticityFunctionalAssembler<GridType> LocalLinearization;
    typedef GeomNonlinLinearizedStVenantAssembler<GridType, Lfe, Lfe> LocalHessian;
    typedef Basis GlobalBasis;

private:
    typedef typename GridType::ctype ctype;

    static const int dim = GridType::dimension;
    static const int dimworld = GridType::dimensionworld;

    typedef typename GridType::template Codim<0>::Geometry::LocalCoordinate LocalCoordinate; 
    typedef typename GridType::template Codim<0>::LeafIterator ElementIterator; 

public:
    GeomExactStVenantMaterial() :
        E_(1.0),
        nu_(0.3)
    {}

    GeomExactStVenantMaterial(const Basis& basis, ctype E, ctype nu) :
        localLinearization_(E,nu),
        localHessian_(E,nu),
        basis_(basis),
        E_(E),
        nu_(nu)
    {}

    void setMaterialParameter(ctype E, ctype nu)
    {
        E_ = E;
        nu_ = nu;
    }

    void setup(ctype E, ctype nu, const Basis& basis)
    {
        E_ = E; 
        nu_ = nu;
        if (localLinearization_)
            localLinearization_.reset();       
        if (localHessian_)
            localHessian_.reset();       
        
        localLinearization_ = Dune::shared_ptr<LocalLinearization>(new LocalLinearization(E,nu));
        localHessian_ = Dune::shared_ptr<LocalHessian>(new LocalHessian(E,nu));

        basis_ = &basis;
    }

    //! Evaluate the strain energy
    template <class CoeffType>
    ctype energy(const CoeffType& coeff)
    {
        // make grid function
        BasisGridFunction<Basis,CoeffType> configuration(*basis_,coeff);

        ctype energy=0;
        const GridType& grid = configuration.grid();

        ElementIterator eIt = grid.template leafbegin<0>();
        ElementIterator eItEnd = grid.template leafend<0>();

        for (;eIt != eItEnd; ++eIt) {

            // get quadrature rule
            QuadratureRuleKey quad1(basis_->getLocalFiniteElement(*eIt));
            QuadratureRuleKey quadKey = quad1.derivative().square().square();

            const Dune::template QuadratureRule<ctype, dim>& quad = QuadratureRuleCache<ctype, dim>::rule(quadKey);

            // loop over quadrature points
            for (size_t pt=0; pt < quad.size(); ++pt) {

                // get quadrature point
                const LocalCoordinate& quadPos = quad[pt].position();

                // get integration factor
                const ctype integrationElement = eIt->geometry().integrationElement(quadPos);

                // evaluate displacement gradient at the quadrature point
                typename BasisGridFunction<Basis,CoeffType>::DerivativeType localConfGrad;

                if (configuration.isDefinedOn(*eIt))
                    configuration.evaluateDerivativeLocal(*eIt, quadPos, localConfGrad);
                else
                    configuration.evaluateDerivative(eIt->geometry().global(quadPos),localConfGrad);

                // Compute the nonlinear strain tensor from the deformation gradient
                SymmetricTensor<dim> strain;
                computeStrain(localConfGrad,strain);

                // and the stress
                SymmetricTensor<dim> stress = hookeTimesStrain(strain);

                ctype z = quad[pt].weight()*integrationElement;
                energy += (stress*strain)*z;

            }
        }
            
        return 0.5*energy;
    }
   
    //! Return the local assembler of the first derivative of the strain energy 
    LocalLinearization& firstDerivative() {return *localLinearization_;}

    //! Return the local assembler of the second derivative of the strain energy 
    LocalHessian& secondDerivative() {return *localHessian_;}

    //! Return the global basis
    const Basis& basis() {return *basis_;}

private:
    //! First derivative of the strain energy
    Dune::shared_ptr<LocalLinearization> localLinearization_;

    //! Second derivative of the strain energy
    Dune::shared_ptr<LocalHessian> localHessian_;

    //! Global basis used for the spatial discretization
    const Basis* basis_;

    //! Elasticity modulus
    ctype E_;
    //! Poisson ratio
    ctype nu_;
    

    //! Compute nonlinear strain tensor from the deformation gradient.
    void computeStrain(const Dune::FieldMatrix<ctype, dim, dim>& grad, SymmetricTensor<dim>& strain) const {
        strain = 0;
        for (int i=0; i<dim ; ++i) {
            strain(i,i) +=grad[i][i];

            for (int k=0;k<dim;++k)
                strain(i,i) += 0.5*grad[k][i]*grad[k][i];

            for (int j=i+1; j<dim; ++j) {
                strain(i,j) += 0.5*(grad[i][j] + grad[j][i]);
                for (int k=0;k<dim;++k)
                    strain(i,j) += 0.5*grad[k][i]*grad[k][j];
            }
        }
    }

    //! Compute linear elastic stress from the strain
    SymmetricTensor<dim> hookeTimesStrain(const SymmetricTensor<dim>& strain) const {

        SymmetricTensor<dim> stress = strain;
        stress *= E_/(1+nu_);

        // Compute the trace of the strain
        ctype trace = 0;
        for (int i=0; i<dim; i++)
            trace += strain(i,i);

        ctype f = E_*nu_ / ((1+nu_)*(1-2*nu_)) * trace;

        for (int i=0; i<dim; i++)
            stress(i,i) += f;

        return stress;
    }


};

#endif