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

Add the old implementation again, it is much much faster than Adol-C

parent b698c106
No related branches found
No related tags found
No related merge requests found
Show whitespace changes
Inline Side-by-side
dune/elasticity/materials/neohookeanmaterial.hh
+ 101
− 84
View file @ 4a910cf6
......@@ -3,15 +3,15 @@
#ifndef DUNE_ELASTICITY_MATERIALS_NEO_HOOKEAN_MATERIAL_HH
#define DUNE_ELASTICITY_MATERIALS_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/elasticity/materials/material.hh>
#include <dune/elasticity/assemblers/neohookefunctionalassembler.hh>
#include <dune/elasticity/assemblers/neohookeoperatorassembler.hh>
/** \brief Class representing a hyperelastic Neo-hookian material.
*
......@@ -34,127 +34,144 @@
* - \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>
template <class Basis>
class NeoHookeanMaterial : public Material<Basis>
{
public:
typedef Material<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;
enum {dim = GridType::dimension};
private:
using Base::dim;
typedef typename GridType::ctype ctype;
typedef typename GridType::template Codim<0>::Entity Element;
typedef NeoHookeFunctionalAssembler<GridType,Lfe> NeoLinearization;
typedef NeoHookeOperatorAssembler<GridType, Lfe, Lfe> NeoHessian;
typedef typename GridType::template Codim<0>::Geometry::LocalCoordinate LocalCoordinate;
typedef typename GridType::template Codim<0>::LeafIterator ElementIterator;
typedef typename LocalFiniteElement::Traits::LocalBasisType::Traits::RangeFieldType field_type;
typedef typename LocalFiniteElement::Traits::LocalBasisType::Traits::JacobianType JacobianType;
public:
NeoHookeanMaterial() :
lambda_(1.0),
mu_(0.3)
{}
typedef Adolc::LocalEnergy<GridType, LocalFiniteElement, dim> Base;
typedef typename Base::CoefficientVectorType CoefficientVectorType;
typedef typename Base::ReturnType ReturnType;
NeoHookeanMaterial(const Basis& basis, ReturnType E, ReturnType nu) :
Base(basis)
{
lambda_ = E*nu / ((1+nu)*(1-2*nu));
mu_ = E / (2*(1+nu));
LocalNeoHookeanEnergy(field_type E, field_type nu)
localLinearization_ = std::shared_ptr<NeoLinearization>(new NeoLinearization(lambda_,mu_));
localHessian_ = std::shared_ptr<NeoHessian>(new NeoHessian(lambda_,mu_));
}
void setup(ReturnType E, ReturnType nu, const Basis& basis)
{
lambda_ = E*nu / ((1+nu)*(1-2*nu));
mu_ = E / (2*(1+nu));
localLinearization_ = std::shared_ptr<NeoLinearization>(new NeoLinearization(lambda_,mu_));
localHessian_ = std::shared_ptr<NeoHessian>(new NeoHessian(lambda_,mu_));
this->basis_ = &basis;
}
ReturnType energy(const Element& element, const LocalFiniteElement& lfe,
const CoefficientVectorType& localCoeff) const
//! Evaluate the strain energy
ReturnType energy(std::shared_ptr<GridFunction> displace)
{
ReturnType energy=0;
const GridType& grid = this->basis().getGridView().grid();
// TODO Get proper quadrature rule
// get quadrature rule
const int order = (element.type().isSimplex()) ? 4 : 4*dim;
ElementIterator eIt = grid.template leafbegin<0>();
ElementIterator eItEnd = grid.template leafend<0>();
const auto& quad = QuadratureRuleCache<ctype, dim>::rule(element.type(),order,
IsRefinedLocalFiniteElement<LocalFiniteElement>::value(lfe));
for (;eIt != eItEnd; ++eIt) {
// the element geometry mapping
const typename Element::Geometry geometry = element.geometry();
// TODO Get proper quadrature rule
// get quadrature rule
const int order = (eIt->type().isSimplex()) ? 4 : 4*dim;
// store gradients of shape functions and base functions
std::vector<JacobianType> referenceGradients(lfe.localBasis().size());
const Dune::template QuadratureRule<ctype, dim>& quad = QuadratureRuleCache<ctype, dim>::rule(eIt->type(),order,0);
// 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);
const LocalCoordinate& quadPos = quad[pt].position();
// 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];
}
const ctype integrationElement = eIt->geometry().integrationElement(quadPos);
SymmetricTensor<dim, ReturnType> strain(0.0);
Dune::Elasticity::computeNonlinearStrain(localGradient,strain);
// evaluate displacement gradient at the quadrature point
typename BasisGridFunction<Basis,VectorType>::DerivativeType localDispGrad;
if (displace->isDefinedOn(*eIt))
displace->evaluateDerivativeLocal(*eIt, quadPos, localDispGrad);
else
displace->evaluateDerivative(eIt->geometry().global(quadPos),localDispGrad);
SymmetricTensor<dim> strain;
Dune::Elasticity::computeNonlinearStrain(localDispGrad,strain);
// the trace
ReturnType trE(0);
for (int i=0; i<dim; i++)
trE += strain(i,i);
// turn displacement gradient into deformation gradient
for (int i=0;i<dim;i++)
localGradient[i][i] += 1;
localDispGrad[i][i] += 1;
// evaluate the derminante of the deformation gradient
const ReturnType J = localGradient.determinant();
const ReturnType J = localDispGrad.determinant();
ReturnType z = quad[pt].weight()*integrationElement;
ctype 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());
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_*strain.trace());
energy += z*(0.5*lambda_*(J-1)*(J-1)-2*mu_*std::log(J)+mu_*trE);
#endif
}
}
return energy;
}
private:
//! First Lame constant
field_type lambda_;
//! Second Lame constant
field_type mu_;
};
//! Return the local assembler of the first derivative of the strain energy
LocalLinearization& firstDerivative(std::shared_ptr<GridFunction> displace) {
template <class Basis>
class NeoHookeMaterial : 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;
localLinearization_->setConfiguration(displace);
return *localLinearization_;
}
typedef LocalNeoHookeanEnergy<GridType,Lfe> LocalEnergy;
using Base::dim;
//! Return the local assembler of the second derivative of the strain energy
LocalHessian& secondDerivative(std::shared_ptr<GridFunction> displace) {
NeoHookeMaterial(const Basis& basis, ReturnType E, ReturnType nu, bool vectorMode=true) :
localEnergy_(E,nu)
{
this->setup(basis, localEnergy_,vectorMode);
localHessian_->setConfiguration(displace);
return *localHessian_;
}
using Base::energy;
private:
LocalEnergy localEnergy_;
//! First derivative of the strain energy
std::shared_ptr<NeoLinearization> localLinearization_;
//! Second derivative of the strain energy
std::shared_ptr<NeoHessian> localHessian_;
//! First Lame constant
ReturnType lambda_;
//! Second Lame constant
ReturnType mu_;
};
#endif
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment