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

Major cleanup

parent 709c603d
Branches
Tags
No related merge requests found
...@@ -34,9 +34,9 @@ public: ...@@ -34,9 +34,9 @@ public:
private: private:
using Base::dim; using Base::dim;
typedef typename GridType::ctype ctype; typedef typename GridType::ctype ctype;
typedef GeomExactStVenantKirchhoffFunctionalAssembler<GridType,Lfe> GeomLinearization; typedef GeomExactStVenantKirchhoffFunctionalAssembler<GridType,Lfe> LinearisationAssembler;
typedef GeomExactStVenantKirchhoffOperatorAssembler<GridType, Lfe, Lfe> GeomHessian; typedef GeomExactStVenantKirchhoffOperatorAssembler<GridType, Lfe, Lfe> HessianAssembler;
typedef typename GridType::template Codim<0>::Geometry::LocalCoordinate LocalCoordinate;
using Derivative = typename GridFunction::DerivativeType; using Derivative = typename GridFunction::DerivativeType;
using StressFunction = std::function<SymmetricTensor<dim, ReturnType>(Derivative)>; using StressFunction = std::function<SymmetricTensor<dim, ReturnType>(Derivative)>;
...@@ -51,8 +51,8 @@ public: ...@@ -51,8 +51,8 @@ public:
lambda_(E*nu/((1+nu)*(1-2*nu))), lambda_(E*nu/((1+nu)*(1-2*nu))),
mu_(E/(2+2*nu)) mu_(E/(2+2*nu))
{ {
localLinearization_ = std::make_shared<GeomLinearization>(E,nu); localLinearisation_.setStiffnessParameter(E, nu);
localHessian_ = std::make_shared<GeomHessian>(E,nu); localHessian_.setStiffnessParameter(E, nu);
} }
void setup(ReturnType E, ReturnType nu, const Basis& basis) void setup(ReturnType E, ReturnType nu, const Basis& basis)
...@@ -61,8 +61,8 @@ public: ...@@ -61,8 +61,8 @@ public:
lambda_ = E*nu/((1+nu)*(1-2*nu)); lambda_ = E*nu/((1+nu)*(1-2*nu));
mu_ = E/(2+2*nu); mu_ = E/(2+2*nu);
localLinearization_ = std::make_shared<GeomLinearization>(E,nu); localLinearisation_.setStiffnessParameter(E, nu);
localHessian_ = std::make_shared<GeomHessian>(E,nu); localHessian_.setStiffnessParameter(E, nu);
this->basis_ = &basis; this->basis_ = &basis;
} }
...@@ -70,9 +70,9 @@ public: ...@@ -70,9 +70,9 @@ public:
//! Evaluate the strain energy //! Evaluate the strain energy
ReturnType energy(std::shared_ptr<GridFunction> displace) const ReturnType energy(std::shared_ptr<GridFunction> displace) const
{ {
ReturnType energy=0; ReturnType energy{0};
const auto gridView = this->basis_->getGridView(); const auto& gridView = this->basis_->getGridView();
for (const auto& e : elements(gridView)) { for (const auto& e : elements(gridView)) {
...@@ -82,55 +82,45 @@ public: ...@@ -82,55 +82,45 @@ public:
const auto& quad = QuadratureRuleCache<ctype, dim>::rule(quadKey); const auto& quad = QuadratureRuleCache<ctype, dim>::rule(quadKey);
// get the element geometry
const auto& geometry = e.geometry(); const auto& geometry = e.geometry();
// loop over quadrature points for (const auto& pt : quad) {
for (size_t pt=0; pt < quad.size(); ++pt) {
// get quadrature point
const LocalCoordinate& quadPos = quad[pt].position();
// get integration factor
const ctype integrationElement = geometry.integrationElement(quadPos);
// evaluate displacement gradient at the quadrature point const auto& quadPos = pt.position();
typename GridFunction::DerivativeType localDispGrad;
Derivative localDispGrad;
if (displace->isDefinedOn(e)) if (displace->isDefinedOn(e))
displace->evaluateDerivativeLocal(e, quadPos, localDispGrad); displace->evaluateDerivativeLocal(e, quadPos, localDispGrad);
else else
displace->evaluateDerivative(geometry.global(quadPos), localDispGrad); displace->evaluateDerivative(geometry.global(quadPos), localDispGrad);
// Compute the nonlinear strain tensor from the deformation gradient
SymmetricTensor<dim, ReturnType> strain; SymmetricTensor<dim, ReturnType> strain;
Dune::Elasticity::strain(localDispGrad, strain); Dune::Elasticity::strain(localDispGrad, strain);
// and the stress
SymmetricTensor<dim, ReturnType> stress = strain; SymmetricTensor<dim, ReturnType> stress = strain;
stress *= 2*mu_; stress *= 2*mu_;
stress.addToDiag(lambda_ * strain.trace()); stress.addToDiag(lambda_ * strain.trace());
ctype z = quad[pt].weight()*integrationElement; ReturnType z = pt.weight() * geometry.integrationElement(quadPos);
energy += (stress*strain)*z; energy += (stress*strain)*z;
} }
} }
return 0.5*energy; return 0.5*energy;
} }
//! Return the local assembler of the first derivative of the strain energy //! Return the local assembler of the first derivative of the strain energy
LocalLinearization& firstDerivative(std::shared_ptr<GridFunction> displace) LinearisationAssembler& firstDerivative(std::shared_ptr<GridFunction> displace)
{ {
localLinearization_->setConfiguration(displace); localLinearisation_.setConfiguration(displace);
return *localLinearization_; return localLinearisation_;
} }
//! Return the local assembler of the second derivative of the strain energy //! Return the local assembler of the second derivative of the strain energy
LocalHessian& secondDerivative(std::shared_ptr<GridFunction> displace) HessianAssembler& secondDerivative(std::shared_ptr<GridFunction> displace)
{ {
localHessian_->setConfiguration(displace); localHessian_.setConfiguration(displace);
return *localHessian_; return localHessian_;
}
static StressFunction stressFunction(ReturnType E, ReturnType nu) { static StressFunction stressFunction(ReturnType E, ReturnType nu) {
ReturnType lambda{E * nu / (1.0+nu) / (1.0-2.0*nu)}; ReturnType lambda{E * nu / (1.0+nu) / (1.0-2.0*nu)};
...@@ -151,10 +141,10 @@ public: ...@@ -151,10 +141,10 @@ public:
private: private:
//! First derivative of the strain energy //! First derivative of the strain energy
std::shared_ptr<GeomLinearization> localLinearization_; LinearisationAssembler localLinearisation_;
//! Second derivative of the strain energy //! Second derivative of the strain energy
std::shared_ptr<GeomHessian> localHessian_; HessianAssembler localHessian_;
//! 1. Lam\'e constant //! 1. Lam\'e constant
ReturnType lambda_; ReturnType lambda_;
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment