Commit 085ebe8b authored by akbib's avatar akbib Committed by akbib@FU-BERLIN.DE
Browse files

Add new folder for hyperelastic materials.

A material has a function "energy" to compute the strain energy and methods first/secondDerivative which return local assemblers of the linearization and the hessian of the strain energy. 

[[Imported from SVN: r11174]]
parent e4acbdd2
materialsdir = $(includedir)/dune/elasticity/materials
materials_HEADERS = geomexactstvenantkirchhoffmaterial.hh
include $(top_srcdir)/am/global-rules
#include <dune/fufem/functiontools/quadraturerulecache.hh>
#include <dune/fufem/symmetrictensor.hh>
#include <dune/fufem/functions/virtualgridfunction.hh>
#include <dune/fufem/assemblers/localassemblers/geononlinstvenantfunctionalassembler.hh>
#include <dune/fufem/assemblers/localassemblers/geononlinlinearizedstvenantassembler.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
typedef typename Basis::GridView::Grid GridType;
typedef typename GridType::ctype ctype;
static const int dim = GridType::dimension;
static const int dimworld = GridType::dimensionworld;
typedef typename Basis::LocalFiniteElement Lfe;
typedef VirtualGridFunction<GridType, Dune::FieldVector<double,GridType::dimension> > GridFunction;
typedef typename GridType::template Codim<0>::Geometry::LocalCoordinate LocalCoordinate;
typedef typename GridType::template Codim<0>::LeafIterator ElementIterator;
typedef GeomNonlinElasticityFunctionalAssembler<GridType> LocalLinearization;
typedef GeomNonlinLinearizedStVenantAssembler<GridType, Lfe, Lfe> LocalHessian;
GeomExactStVenantMaterial(const Basis& basis, ctype E, ctype nu) :
//! 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) {
// get quadrature rule
QuadratureRuleKey quad1(basis_.getLocalFiniteElement(*eIt));
QuadratureRuleKey quadKey = quad1.derivative().square().square();
const Dune::template QuadratureRule<ctype, dim>& quad = QuadratureRuleCache<ctype, dim>::rule(quadKey);
// 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 GridFunction::DerivativeType localConfGrad;
if (configuration.isDefinedOn(*eIt))
configuration.evaluateDerivativeLocal(*eIt, quadPos, localConfGrad);
// Compute the nonlinear strain tensor from the deformation gradient
SymmetricTensor<dim> strain;
// and the stress
SymmetricTensor<dim> stress = hookeTimesStrain(strain);
ctype z = quad[pt].weight()*integrationElement;
energy += (stress*strain)*z;
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_;}
//! First derivative of the strain energy
LocalLinearization localLinearization_;
//! Second derivative of the strain energy
LocalHessian localHessian_;
//! Global basis used for the spatial discretization
const Basis& basis_;
//! Elasticity modulus
const ctype E_;
//! Poisson ratio
const 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];
//! Compute linear elastic stress from the strain
SymmetricTensor<dim> hookeTimesStrain(const SymmetricTensor<dim>& strain) const {
SymmetricTensor<dim> 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;
for (int i=0; i<dim; i++)
stress(i,i) += f;
return stress;
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment