Inherit from new abstract base class + cleanup and renaming.

[[Imported from SVN: r12666]]

dune/elasticity/materials/geomexactstvenantkirchhoffmaterial.hh

-#ifndef GEOMETRICALLY_EXACT_ST_VENANT_KIRCHHOFF_MATERIAL
-#define GEOMETRICALLY_EXACT_ST_VENANT_KIRCHHOFF_MATERIAL
+// -*- tab-width: 8; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set ts=8 sw=4 et sts=4:
+#ifndef DUNE_ELASTICITY_MATERIALS_GEOMETRICALLY_EXACT_ST_VENANT_KIRCHHOFF_MATERIAL_HH
+#define DUNE_ELASTICITY_MATERIALS_GEOMETRICALLY_EXACT_ST_VENANT_KIRCHHOFF_MATERIAL_HH
 
 #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/fufem/assemblers/localassemblers/geononlinstvenantfunctionalassembler.hh>
 #include <dune/fufem/assemblers/localassemblers/geononlinlinearizedstvenantassembler.hh>
 
+#include <dune/elasticity/materials/material.hh>
+#include <dune/elasticity/common/elasticityhelpers.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
+class GeomExactStVenantMaterial : public Material<Basis>
 {
+    typedef Material<Basis> Base;
 public:
     typedef typename Basis::GridView::Grid GridType;
-    typedef typename Basis::LocalFiniteElement Lfe;    
+    typedef typename Base::GlobalBasis GlobalBasis;
+    typedef typename Base::Lfe Lfe;
+    typedef typename Base::LocalLinearization LocalLinearization;
+    typedef typename Base::LocalHessian LocalHessian;
+    typedef typename Base::VectorType VectorType;
+    typedef typename Base::GridFunction GridFunction;
+    typedef typename Base::ReturnType ReturnType;
 
-    typedef GeomNonlinElasticityFunctionalAssembler<GridType,Lfe> LocalLinearization;
-    typedef GeomNonlinLinearizedStVenantAssembler<GridType, Lfe, Lfe> LocalHessian;
-    typedef Basis GlobalBasis;
 
 private:
+    using Base::dim;
     typedef typename GridType::ctype ctype;
-
-    static const int dim = GridType::dimension;
-    static const int dimworld = GridType::dimensionworld;
-
+    typedef GeomNonlinElasticityFunctionalAssembler<GridType,Lfe> GeomLinearization;
+    typedef GeomNonlinLinearizedStVenantAssembler<GridType, Lfe, Lfe> GeomHessian;
     typedef typename GridType::template Codim<0>::Geometry::LocalCoordinate LocalCoordinate; 
     typedef typename GridType::template Codim<0>::LeafIterator ElementIterator; 
 
@@ -41,51 +46,31 @@ public:
         nu_(0.3)
     {}
 
-    GeomExactStVenantMaterial(const Basis& basis, ctype E, ctype nu) :
-        localLinearization_(E,nu),
-        localHessian_(E,nu),
-        basis_(basis),
+    GeomExactStVenantMaterial(const Basis& basis, ReturnType E, ReturnType nu) :
+        Base(basis),
         E_(E),
         nu_(nu)
-    {}
-
-    void getMaterialParameters(ctype& E, ctype& nu)
-    {
-        E = E_;
-        nu = nu_;
-    }
-
-
-    void setMaterialParameters(ctype E, ctype nu)
     {
-        E_ = E;
-        nu_ = nu;
+        localLinearization_ = Dune::shared_ptr<GeomLinearization>(new GeomLinearization(E,nu));
+        localHessian_ = Dune::shared_ptr<GeomHessian>(new GeomHessian(E,nu));
     }
 
-    void setup(ctype E, ctype nu, const Basis& basis)
+    void setup(ReturnType E, ReturnType 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));
+        localLinearization_ = Dune::shared_ptr<GeomLinearization>(new GeomLinearization(E,nu));
+        localHessian_ = Dune::shared_ptr<GeomHessian>(new GeomHessian(E,nu));
 
-        basis_ = &basis;
+        this->basis_ = &basis;
     }
 
     //! Evaluate the strain energy
-    template <class CoeffType>
-    ctype energy(const CoeffType& coeff)
+    ReturnType energy(const GridFunction& displace)
     {
-        // make grid function
-        BasisGridFunction<Basis,CoeffType> configuration(*basis_,coeff);
-
         ctype energy=0;
-        const GridType& grid = configuration.grid();
+        const GridType& grid = this->basis_->getGridView().grid();
 
         ElementIterator eIt = grid.template leafbegin<0>();
         ElementIterator eItEnd = grid.template leafend<0>();
@@ -93,7 +78,7 @@ public:
         for (;eIt != eItEnd; ++eIt) {
 
             // get quadrature rule
-            QuadratureRuleKey quad1(basis_->getLocalFiniteElement(*eIt));
+            QuadratureRuleKey quad1(this->basis_->getLocalFiniteElement(*eIt));
             QuadratureRuleKey quadKey = quad1.derivative().square().square();
 
             const Dune::template QuadratureRule<ctype, dim>& quad = QuadratureRuleCache<ctype, dim>::rule(quadKey);
@@ -108,23 +93,22 @@ public:
                 const ctype integrationElement = eIt->geometry().integrationElement(quadPos);
 
                 // evaluate displacement gradient at the quadrature point
-                typename BasisGridFunction<Basis,CoeffType>::DerivativeType localConfGrad;
+                typename GridFunction::DerivativeType localDispGrad;
 
-                if (configuration.isDefinedOn(*eIt))
-                    configuration.evaluateDerivativeLocal(*eIt, quadPos, localConfGrad);
+                if (displace->isDefinedOn(*eIt))
+                    displace->evaluateDerivativeLocal(*eIt, quadPos, localDispGrad);
                 else
-                    configuration.evaluateDerivative(eIt->geometry().global(quadPos),localConfGrad);
+                    displace->evaluateDerivative(eIt->geometry().global(quadPos),localDispGrad);
 
                 // Compute the nonlinear strain tensor from the deformation gradient
-                SymmetricTensor<dim> strain;
-                computeStrain(localConfGrad,strain);
+                SymmetricTensor<dim,ReturnType> strain;
+                Dune::Elasticity::computeNonlinearStrain(localDispGrad,strain);
 
                 // and the stress
-                SymmetricTensor<dim> stress = hookeTimesStrain(strain);
+                SymmetricTensor<dim,ReturnType> stress = hookeTimesStrain(strain);
 
                 ctype z = quad[pt].weight()*integrationElement;
                 energy += (stress*strain)*z;
-
             }
         }
             
@@ -132,67 +116,44 @@ public:
     }
    
     //! Return the local assembler of the first derivative of the strain energy 
-    LocalLinearization& firstDerivative() {return *localLinearization_;}
+    LocalLinearization& firstDerivative(const GridFunction& displace)
+    {
+        localLinearization_->setConfiguration(displace);
+        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_;}
+    LocalHessian& secondDerivative(const GridFunction& displace)
+    {
+        localHessian_->setConfiguration(displace);
+        return *localHessian_;
+    }
 
 private:
     //! First derivative of the strain energy
-    Dune::shared_ptr<LocalLinearization> localLinearization_;
+    Dune::shared_ptr<GeomLinearization> localLinearization_;
 
     //! Second derivative of the strain energy
-    Dune::shared_ptr<LocalHessian> localHessian_;
-
-    //! Global basis used for the spatial discretization
-    const Basis* basis_;
+    Dune::shared_ptr<GeomHessian> localHessian_;
 
     //! Elasticity modulus
-    ctype E_;
+    ReturnType 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];
-            }
-        }
-    }
+    ReturnType nu_;
     
     //! Compute linear elastic stress from the strain
-    SymmetricTensor<dim> hookeTimesStrain(const SymmetricTensor<dim>& strain) const {
+    SymmetricTensor<dim,ReturnType> hookeTimesStrain(const SymmetricTensor<dim,ReturnType>& strain) const {
 
-        SymmetricTensor<dim> stress = strain;
+        SymmetricTensor<dim,ReturnType> 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;
+        ReturnType f = E_*nu_ / ((1+nu_)*(1-2*nu_)) * strain.trace();
 
         for (int i=0; i<dim; i++)
             stress(i,i) += f;
 
         return stress;
     }
-
-
 };
 
 #endif