Material class for Neo-Hookean materials

[[Imported from SVN: r11295]]

dune/elasticity/materials/neohookeanmaterial.hh

+#ifndef NEO_HOOKIAN_MATERIAL
+#define NEO_HOOKIAN_MATERIAL
+
+#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/elasticity/assemblers/neohookefunctionalassembler.hh>
+#include <dune/elasticity/assemblers/neohookeoperatorassembler.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 :-( )
+ */
+template <class Basis>
+class NeoHookeMaterial
+{
+public:
+    typedef typename Basis::GridView::Grid GridType;
+    typedef typename Basis::LocalFiniteElement Lfe;    
+
+    typedef NeoHookeFunctionalAssembler<GridType> LocalLinearization;
+    typedef NeoHookeOperatorAssembler<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:
+    NeoHookeMaterial() :
+        lambda_(1.0),
+        mu_(0.3)
+    {}
+
+    NeoHookeMaterial(const Basis& basis, ctype E, ctype nu) :
+        basis_(basis)
+    {      
+        lambda_ = E*nu / ((1+nu)*(1-2*nu));
+        mu_     = E / (2*(1+nu));
+
+        localLinearization_ = Dune::shared_ptr<LocalLinearization>(new LocalLinearization(lambda_,mu_));
+        localHessian_ = Dune::shared_ptr<LocalHessian>(new LocalHessian(lambda_,mu_));
+    }
+
+
+    void setup(ctype E, ctype nu, const Basis& basis)
+    {
+        lambda_ = E*nu / ((1+nu)*(1-2*nu));
+        mu_     = E / (2*(1+nu));
+
+        if (localLinearization_)
+            localLinearization_.reset();       
+        if (localHessian_)
+            localHessian_.reset();       
+        
+        localLinearization_ = Dune::shared_ptr<LocalLinearization>(new LocalLinearization(lambda_,mu_));
+        localHessian_ = Dune::shared_ptr<LocalHessian>(new LocalHessian(lambda_,mu_));
+
+        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) {
+
+            // TODO Get proper quadrature rule
+            // get quadrature rule
+            const int order = (eIt->type().isSimplex()) ? 4 : 4*dim;
+
+            const Dune::template QuadratureRule<ctype, dim>& quad = QuadratureRuleCache<ctype, dim>::rule(order);
+
+            // 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);
+
+                SymmetricTensor<dim> strain;
+                computeStrain(localConfGrad,strain);
+            
+                // the trace
+                ctype trE(0);
+                for (int i=0; i<dim; i++)
+                    trE += strain(i,i);                
+
+                // make deformation gradient out of the discplacement
+                for (int i=0;i<dim;i++)
+                    localConfGrad[i][i] += 1;        
+
+                // evaluate the derminante of the deformation gradient
+                const ctype J = localConfGrad.determinant();
+
+                ctype z = quad[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 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_;
+
+    //! Lame constant
+    ctype lambda_;
+    //! Lame constant
+   ctype mu_;
+    
+
+    //! 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];
+            }
+        }
+    }
+
+};
+
+#endif