Skip to content
Snippets Groups Projects
Commit 86c702d5 authored by Jonathan Youett's avatar Jonathan Youett
Browse files

Remove redundant class

The normal material classes can now fulfill this purpuse
parent 4a335a72
Branches
Tags
No related merge requests found
// -*- tab-width: 8; indent-tabs-mode: nil; c-basic-offset: 4 -*-
// vi: set ts=8 sw=4 et sts=4:
#ifndef DUNE_ELASTICITY_MATERIALS_ADOLC_NEO_HOOKEAN_MATERIAL_HH
#define DUNE_ELASTICITY_MATERIALS_ADOLC_NEO_HOOKEAN_MATERIAL_HH
#include <dune/fufem/assemblers/localassemblers/adolclocalenergy.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/elasticity/materials/adolcmaterial.hh>
#include <dune/elasticity/common/elasticityhelpers.hh>
#include <dune/matrix-vector/addtodiagonal.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 :-( )
*
* * 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]
*
* 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 (first lame and shear modulus)
*/
template <class GridType, class LocalFiniteElement>
class LocalNeoHookeanEnergy : public Adolc::LocalEnergy<GridType, LocalFiniteElement,GridType::dimension>
{
public:
enum {dim = GridType::dimension};
typedef typename GridType::ctype ctype;
typedef typename GridType::template Codim<0>::Entity Element;
typedef typename LocalFiniteElement::Traits::LocalBasisType::Traits::RangeFieldType field_type;
typedef typename LocalFiniteElement::Traits::LocalBasisType::Traits::JacobianType JacobianType;
typedef Adolc::LocalEnergy<GridType, LocalFiniteElement, dim> Base;
typedef typename Base::CoefficientVectorType CoefficientVectorType;
typedef typename Base::ReturnType ReturnType;
LocalNeoHookeanEnergy() = default;
LocalNeoHookeanEnergy(field_type E, field_type nu)
{
lambda_ = E*nu / ((1+nu)*(1-2*nu));
mu_ = E / (2*(1+nu));
}
LocalNeoHookeanEnergy& operator=(LocalNeoHookeanEnergy&& other) {
lambda_ = other.lambda_;
mu_ = other.mu_;
}
ReturnType energy(const Element& element, const LocalFiniteElement& lfe,
const CoefficientVectorType& localCoeff) const
{
ReturnType energy=0;
// TODO Get proper quadrature rule
// get quadrature rule
const int order = (element.type().isSimplex()) ? 4 : 4*dim;
const auto& quad = QuadratureRuleCache<ctype, dim>::rule(element.type(),order,
IsRefinedLocalFiniteElement<LocalFiniteElement>::value(lfe));
// the element geometry mapping
const typename Element::Geometry geometry = element.geometry();
// store gradients of shape functions and base functions
std::vector<JacobianType> referenceGradients(lfe.localBasis().size());
// loop over quadrature points
for (size_t pt=0; pt < quad.size(); ++pt) {
// get quadrature point
const auto& quadPos = quad[pt].position();
// get transposed inverse of Jacobian of transformation
const auto& invJacobian = geometry.jacobianInverseTransposed(quadPos);
// get integration factor
const ctype integrationElement = geometry.integrationElement(quadPos);
// get gradients of shape functions
lfe.localBasis().evaluateJacobian(quadPos, referenceGradients);
// compute gradient of the configuration
Dune::FieldMatrix<ReturnType,dim,dim> localGradient(0);
for (size_t k=0; k<referenceGradients.size(); ++k) {
Dune::FieldVector<ReturnType,dim> gradient(0);
invJacobian.umv(referenceGradients[k][0], gradient);
for (int i=0; i<dim; ++i)
for (int j=0; j<dim; ++j)
localGradient[i][j] += localCoeff[k][i] * gradient[j];
}
SymmetricTensor<dim, ReturnType> strain(0.0);
Dune::Elasticity::strain(localGradient,strain);
// turn displacement gradient into deformation gradient
Dune::MatrixVector::addToDiagonal(localGradient, 1.0);
// evaluate the derminante of the deformation gradient
const ReturnType J = localGradient.determinant();
ReturnType z = quad[pt].weight()*integrationElement;
#ifdef LAURSEN
energy += z*(0.25*lambda_*(J*J-1)-(lambda_*0.5+mu_)*std::log(J)+mu_*strain.trace());
#else
energy += z*(0.5*lambda_*(J-1)*(J-1)-2*mu_*std::log(J)+mu_*strain.trace());
#endif
}
return energy;
}
private:
//! First Lame constant
field_type lambda_;
//! Second Lame constant
field_type mu_;
};
template <class Basis>
class AdolcNeoHookeanMaterial : public AdolcMaterial<Basis>
{
public:
typedef AdolcMaterial<Basis> Base;
typedef typename Base::GridType GridType;
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 LocalNeoHookeanEnergy<GridType,Lfe> LocalEnergy;
using Base::dim;
AdolcNeoHookeanMaterial() = default;
template <class BasisT>
AdolcNeoHookeanMaterial(BasisT&& basis, ReturnType E, ReturnType nu, bool vectorMode=true) :
localEnergy_(E,nu)
{
Base::setup(std::forward<BasisT>(basis), localEnergy_,vectorMode);
}
template <class BasisT>
void setup(BasisT&& basis, ReturnType E, ReturnType nu, bool vectorMode=true) {
localEnergy_ = LocalEnergy(E,nu);
Base::setup(std::forward<BasisT>(basis), localEnergy_,vectorMode);
}
using Base::energy;
private:
LocalEnergy localEnergy_;
};
#endif
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment