diff --git a/dune/elasticity/materials/geomexactstvenantkirchhoffmaterial.hh b/dune/elasticity/materials/geomexactstvenantkirchhoffmaterial.hh
index 476162d93ccccd12344bb159a107772489376c6a..c4714622bc745e11bf09f61332c51e5c2eac2d05 100644
--- a/dune/elasticity/materials/geomexactstvenantkirchhoffmaterial.hh
+++ b/dune/elasticity/materials/geomexactstvenantkirchhoffmaterial.hh
@@ -1,37 +1,42 @@
-#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_;
+ localLinearization_ = Dune::shared_ptr<GeomLinearization>(new GeomLinearization(E,nu));
+ localHessian_ = Dune::shared_ptr<GeomHessian>(new GeomHessian(E,nu));
}
-
- void setMaterialParameters(ctype E, ctype nu)
- {
- E_ = E;
- nu_ = 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));
-
- basis_ = &basis;
+
+ localLinearization_ = Dune::shared_ptr<GeomLinearization>(new GeomLinearization(E,nu));
+ localHessian_ = Dune::shared_ptr<GeomHessian>(new GeomHessian(E,nu));
+
+ 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_;
+ ReturnType 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,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