stvenantkirchhoffassembler.hh

#ifndef ST_VENANT_KIRCHHOFF_ASSEMBLER_HH
#define ST_VENANT_KIRCHHOFF_ASSEMBLER_HH


#include <dune/common/fvector.hh>
#include <dune/common/fmatrix.hh>

#include <dune/istl/matrix.hh>

#include "dune/fufem/quadraturerules/quadraturerulecache.hh"

#include "dune/fufem/assemblers/localoperatorassembler.hh"
#include "dune/fufem/staticmatrixtools.hh"

#include <dune/fufem/symmetrictensor.hh>


/** \brief Assembler for a St.-Venant-Kirchhoff material (i.e., linear elasticity)
 */
template <class GridType, class TrialLocalFE, class AnsatzLocalFE>
class StVenantKirchhoffAssembler 
    : public LocalOperatorAssembler < GridType, TrialLocalFE, AnsatzLocalFE, Dune::FieldMatrix<typename GridType::ctype,GridType::dimension,GridType::dimension> >
{
    private:
        static const int dim = GridType::dimension;
        typedef typename GridType::ctype ctype;

    /** \brief Young's modulus */
    double E_;

    /** \brief The Poisson ratio */
    double nu_;

    public:
        typedef typename Dune::FieldMatrix<ctype,GridType::dimension,GridType::dimension> T;

        typedef typename LocalOperatorAssembler < GridType, TrialLocalFE, AnsatzLocalFE, T >::Element Element;
        typedef typename LocalOperatorAssembler < GridType, TrialLocalFE, AnsatzLocalFE,T >::BoolMatrix BoolMatrix;
        typedef typename LocalOperatorAssembler < GridType, TrialLocalFE, AnsatzLocalFE,T >::LocalMatrix LocalMatrix;

    StVenantKirchhoffAssembler(double E, double nu):
        E_(E), nu_(nu)
    {}
    
    void indices(const Element& element, BoolMatrix& isNonZero, const TrialLocalFE& tFE, const AnsatzLocalFE& aFE) const
    {
        isNonZero = true;
    }

    void assemble(const Element& element, LocalMatrix& localMatrix, const TrialLocalFE& tFE, const AnsatzLocalFE& aFE) const
    {
        typedef typename Dune::template FieldVector<ctype,dim> FVdim;
        typedef typename Dune::template FieldMatrix<ctype,dim,dim> FMdimdim;
        typedef typename TrialLocalFE::Traits::LocalBasisType::Traits::JacobianType JacobianType;
        
        int rows = tFE.localBasis().size();
        int cols = aFE.localBasis().size();


        localMatrix.setSize(rows,cols);
        localMatrix = 0.0;
        
        // get quadrature rule
        const int order = (element.type().isSimplex()) 
            ? 2*(tFE.localBasis().order()-1)
            : 2*(tFE.localBasis().order()*dim-1);
        //            const Dune::template QuadratureRule<double, dim>& quad = Dune::template QuadratureRules<double, dim>::rule(element.type(), order);
        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());

        // the element geometry mapping 
        const typename Element::Geometry geometry = element.geometry();

        // loop over quadrature points
        for (size_t pt=0; pt < quad.size(); ++pt) {
            
            // get quadrature point
            const FVdim& 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]);

            }

            // /////////////////////////////////////////////
            //   Compute strain for all shape functions
            // /////////////////////////////////////////////
            std::vector<Dune::array<SymmetricTensor<dim>,dim> > strain(rows);

            for (int i=0; i<rows; i++) {

                for (int k=0; k<dim; k++) {

                    // The deformation gradient of the shape function
                    FMdimdim deformationGradient(0);
                    deformationGradient[k] = gradients[i];
                    
                    /* Computes the linear strain tensor from the deformation gradient*/
                    computeStrain(deformationGradient,strain[i][k]);
                    
                }
                
            }

            // /////////////////////////////////////////////////
            //   Assemble matrix
            // /////////////////////////////////////////////////
            ctype z = quad[pt].weight() * integrationElement;
            for (int row=0; row<rows; ++row) {

                for (int rcomp=0; rcomp<dim; rcomp++) {

                    // Compute stress
                    SymmetricTensor<dim> stress = hookeTimesStrain(strain[row][rcomp]);

                    for (int col=0; col<=row; col++) {
                        for (int ccomp=0; ccomp<dim; ccomp++) {

                            ctype zij = stress*strain[col][ccomp] * z;
                            localMatrix[row][col][rcomp][ccomp] += zij;
                            if (col!=row)
                                localMatrix[col][row][ccomp][rcomp] += zij;
                        }
                    }
                }
            }
            
        }

    }

    void computeStrain(const Dune::FieldMatrix<double, dim, dim>& grad, SymmetricTensor<dim>& strain) const
    {
        for (int i=0; i<dim ; ++i)
            {
                strain(i,i) = grad[i][i];
                for (int j=i+1; j<dim; ++j)
                    strain(i,j) = 0.5*(grad[i][j] + grad[j][i]);
            }
        
    }

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

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

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

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

        return stress;
      }
};


#endif