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

Minor cleanup

parent d3d666b1
No related branches found
No related tags found
No related merge requests found
......@@ -90,15 +90,15 @@ public:
//! Evaluate the strain energy
field_type energy(std::shared_ptr<GridFunction> displace) const
{
ctype energy=0;
const auto& leafView = this->basis().getGridView().grid().leafGridView();
field_type energy=0;
const auto& leafView = this->basis().getGridView();
for (const auto& e : elements(leafView)) {
// TODO Get proper quadrature rule
// get quadrature rule (should depend on k?)
const int order = (e.type().isSimplex()) ? 5 : 5*dim;
const auto& quad = QuadratureRuleCache<ctype, dim>::rule(e.type(), order, 0);
const auto& quad = QuadratureRuleCache<field_type, dim>::rule(e.type(), order, 0);
const auto& geometry = e.geometry();
......@@ -109,7 +109,7 @@ public:
auto integrationElement = geometry.integrationElement(quadPos);
// evaluate displacement gradient at the quadrature point
typename BasisGridFunction<Basis,VectorType>::DerivativeType localConfGrad;
typename BasisGridFunction<Basis, typename Base::VectorType>::DerivativeType localConfGrad;
if (displace->isDefinedOn(e))
displace->evaluateDerivativeLocal(e, quadPos, localConfGrad);
......
......@@ -107,11 +107,8 @@ public:
// loop over quadrature points
for (const auto& pt : quad) {
// get quadrature point
const auto& quadPos = pt.position();
// get integration factor
const ctype integrationElement = geometry.integrationElement(quadPos);
const auto integrationElement = geometry.integrationElement(quadPos);
// evaluate displacement gradient at the quadrature point
typename BasisGridFunction<Basis, typename Base::VectorType>::DerivativeType localDispGrad;
......@@ -119,21 +116,17 @@ public:
if (displace->isDefinedOn(e))
displace->evaluateDerivativeLocal(e, quadPos, localDispGrad);
else
displace->evaluateDerivative(geometry.global(quadPos), localDispGrad);
SymmetricTensor<dim,ReturnType> strain;
Dune::Elasticity::strain(localDispGrad,strain);
displace->evaluateDerivative(geometry.global(quadPos),localDispGrad);
// the trace
auto strain = Dune::Elasticity::strain(localDispGrad);
auto trE = strain.trace();
// turn displacement gradient into deformation gradient
Dune::MatrixVector::addToDiagonal(localDispGrad, 1.0);
// evaluate the derminante of the deformation gradient
const ReturnType J = localDispGrad.determinant();
ctype z = pt.weight()*integrationElement;
auto J = localDispGrad.determinant();
auto z = pt.weight()*integrationElement;
#ifdef LAURSEN
energy += z*(0.25*lambda_*(J*J-1)-(lambda_*0.5+mu_)*std::log(J)+mu_*trE);
......@@ -142,7 +135,6 @@ public:
#endif
}
}
return energy;
}
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment