local assembler for the linearization of neo-hookean energy. Basically a...

local assembler for the linearization of neo-hookean energy. Basically a prettier version of yet existing code

[[Imported from SVN: r11289]]

dune/elasticity/assemblers/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
+
+