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

Cleanup and reduce typedefs

parent 2922693c
Branches
Tags
No related merge requests found
......@@ -48,25 +48,19 @@ public:
using typename Base::GridType;
using typename Base::GlobalBasis;
using typename Base::Lfe;
using typename Base::LocalLinearization;
using typename Base::LocalHessian;
using typename Base::VectorType;
using typename Base::ReturnType;
using typename Base::GridFunction;
using field_type = typename Base::ReturnType;
using Base::dim;
using AdolCBase = Adolc::LocalEnergy<GridType, Lfe, dim>;
using Element = typename GridType::template Codim<0>::Entity;
using AdolcCoefficients = typename AdolCBase::CoefficientVectorType;
using AdolcEnergy = typename AdolCBase::ReturnType;
using AdolcFieldType = typename AdolCBase::ReturnType;
using MonRivLinearization = MooneyRivlinFunctionalAssembler<GridType, Lfe>;
using MonRivHessian = MooneyRivlinOperatorAssembler<GridType, Lfe, Lfe>;
private:
using ctype = typename GridType::ctype;
static constexpr int dim = GridType::dimension;
static constexpr int dimworld = GridType::dimensionworld;
using MonRivLinearisation = MooneyRivlinFunctionalAssembler<GridType, Lfe>;
using MonRivHessian = MooneyRivlinOperatorAssembler<GridType, Lfe, Lfe>;
public:
MooneyRivlinMaterial(int k=3) :
......@@ -77,7 +71,7 @@ public:
}
template <class BasisT>
MooneyRivlinMaterial(BasisT&& basis, ctype E, ctype nu, int k=3) :
MooneyRivlinMaterial(BasisT&& basis, field_type E, field_type nu, int k=3) :
Base(std::forward<BasisT>(basis)), E_(E), nu_(nu), k_(k)
{
setupCoefficients();
......@@ -85,7 +79,7 @@ public:
}
template <class BasisT>
void setup(BasisT&& basis, ctype E, ctype nu, int k=3)
void setup(BasisT&& basis, field_type E, field_type nu, int k=3)
{
this->setBasis(std::forward<BasisT>(basis));
......@@ -94,12 +88,12 @@ public:
constructAssemblers();
}
void getMaterialParameters(ctype& E, ctype& nu) {
void getMaterialParameters(field_type& E, field_type& nu) {
E = E_; nu = nu_;
}
//! Evaluate the strain energy
ReturnType energy(std::shared_ptr<GridFunction> displace) const
field_type energy(std::shared_ptr<GridFunction> displace) const
{
ctype energy=0;
const auto& leafView = this->basis().getGridView().grid().leafGridView();
......@@ -200,13 +194,13 @@ public:
}
//! Return the local assembler of the first derivative of the strain energy
LocalLinearization& firstDerivative(std::shared_ptr<GridFunction> displace) {
typename Base::LocalLinearization& firstDerivative(std::shared_ptr<GridFunction> displace) {
localLinearization_->setConfiguration(displace);
return *localLinearization_;
}
//! Return the local assembler of the second derivative of the strain energy
LocalHessian& secondDerivative(std::shared_ptr<GridFunction> displace) {
typename Base::LocalHessian& secondDerivative(std::shared_ptr<GridFunction> displace) {
localHessian_->setConfiguration(displace);
return *localHessian_;
}
......@@ -215,23 +209,23 @@ private:
//! Compute coefficients s.t. polyconvexity holds
void setupCoefficients() {
ctype lambda = E_*nu_ / ((1 +nu_)*(1 - 2*nu_));
ctype mu = E_ / (2*(1 + nu_));
field_type lambda = E_*nu_ / ((1 +nu_)*(1 - 2*nu_));
field_type mu = E_ / (2*(1 + nu_));
// Choose coefficients in such a way that we have polyconvexity
// First/Second derivative of the compressible function part of gamma (=-kJ^(-k-1)) at 1.0
ctype ld1 = -k_;
ctype ld2 = k_*(k_+1);
field_type ld1 = -k_;
field_type ld2 = k_*(k_+1);
if (ld1 >=0 || ld2 <= 0)
std::cout<<"Coerciveness failed\n";
// Check if lame constants admit coerciveness
ctype rho= -ld1/(ld2-ld1);
field_type rho= -ld1/(ld2-ld1);
if( ( (rho < 0.5 && lambda < (1.0/rho-2.0)*mu) || lambda <= 0 || mu <= 0 ))
std::cout<<"Coerciveness failed\n";
//const ctype somePositiveConstant = - mu_ + rho*(lambda_+2.*mu_);
//const field_type somePositiveConstant = - mu_ + rho*(lambda_+2.*mu_);
//if(somePositiveConstant <= 0)
e_ = (lambda+2.0*mu)/(ld2-ld1);
......@@ -249,7 +243,7 @@ private:
a_ = b_ + mu;
// last check if I didn't miss a condition
ctype alpha = 0.5*(mu - b_);
field_type alpha = 0.5*(mu - b_);
if(alpha <= 0 || b_ <= 0 || d_ <= 0)
std::cout<<"Coerciveness failed\n";
}
......@@ -267,13 +261,13 @@ private:
std::shared_ptr<MonRivHessian> localHessian_;
//! Elasticity modulus
ctype E_;
field_type E_;
//! Shear modulus
ctype nu_;
field_type nu_;
//! Exponent of the compressibility function (>= 2)
int k_;
//! Material parameters
ctype a_; ctype b_; ctype c_; ctype d_; ctype e_;
field_type a_; field_type b_; field_type c_; field_type d_; field_type e_;
};
#endif
......@@ -45,15 +45,13 @@ class NeoHookeanMaterial : public Material<Basis>,
Material<Basis>::GridType::dimension>
{
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;
using Base = Material<Basis>;
using typename Base::GridType;
using typename Base::GlobalBasis;
using typename Base::Lfe;
using typename Base::GridFunction;
using field_type = typename Base::ReturnType;
using Base::dim;
using AdolCBase = Adolc::LocalEnergy<GridType, Lfe, dim>;
using Element = typename GridType::template Codim<0>::Entity;
......@@ -62,12 +60,8 @@ public:
using AdolcFieldType = typename AdolCBase::ReturnType;
private:
using Base::dim;
typedef typename GridType::ctype ctype;
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;
using NeoLinearisation = NeoHookeFunctionalAssembler<GridType,Lfe>;
using NeoHessian = NeoHookeOperatorAssembler<GridType, Lfe, Lfe>;
public:
NeoHookeanMaterial() :
......@@ -87,7 +81,7 @@ public:
}
template <class BasisT>
void setup(BasisT&& basis, ReturnType E, ReturnType nu)
void setup(BasisT&& basis, field_type E, field_type nu)
{
lambda_ = E*nu / ((1+nu)*(1-2*nu));
......@@ -100,10 +94,10 @@ public:
}
//! Evaluate the strain energy
ReturnType energy(std::shared_ptr<GridFunction> displace) const
field_type energy(std::shared_ptr<GridFunction> displace) const
{
ReturnType 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)) {
......@@ -129,7 +123,7 @@ public:
const ctype integrationElement = geometry.integrationElement(quadPos);
// evaluate displacement gradient at the quadrature point
typename BasisGridFunction<Basis,VectorType>::DerivativeType localDispGrad;
typename BasisGridFunction<Basis, typename Base::VectorType>::DerivativeType localDispGrad;
if (displace->isDefinedOn(e))
displace->evaluateDerivativeLocal(e, quadPos, localDispGrad);
......@@ -219,14 +213,14 @@ public:
}
//! Return the local assembler of the first derivative of the strain energy
LocalLinearization& firstDerivative(std::shared_ptr<GridFunction> displace) {
typename Base::LocalLinearization& firstDerivative(std::shared_ptr<GridFunction> displace) {
localLinearization_->setConfiguration(displace);
return *localLinearization_;
}
//! Return the local assembler of the second derivative of the strain energy
LocalHessian& secondDerivative(std::shared_ptr<GridFunction> displace) {
typename Base::LocalHessian& secondDerivative(std::shared_ptr<GridFunction> displace) {
localHessian_->setConfiguration(displace);
return *localHessian_;
......@@ -240,9 +234,9 @@ private:
std::shared_ptr<NeoHessian> localHessian_;
//! First Lame constant
ReturnType lambda_;
field_type lambda_;
//! Second Lame constant
ReturnType mu_;
field_type mu_;
};
#endif
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment