Commit b80f2e49 authored by jonathan.froehlich's avatar jonathan.froehlich
Browse files

Merge branch '143-use-optimized-elements' into 'feature'

Resolve "Use optimized elements"

Closes #143

See merge request !164
parents 90b32d4f 4b870dd0
Pipeline #221153 failed
......@@ -6,6 +6,7 @@ build/
# == Build Directories from CLion ==
.idea/
cmake-build-debug/
cmake-build-debug-coverage/
.build-clion/
# == Python directories ==
tools/venv/
......
#include <problem/PressureProblem.hpp>
#include <LagrangeDiscretization.hpp>
#include "VectorFieldElement.hpp"
#include "VectorFieldFaceElement.hpp"
#include "LagrangeElasticity.hpp"
#include "LagrangeTransfer.hpp"
......@@ -123,10 +122,10 @@ void LagrangeElasticity::Residual(const cell &c, const Vector &U, Vector &r) con
for (int i = 0; i < E.size(); ++i) {
for (int k = 0; k < c.dim(); ++k) {
VectorField u_ik = E.VectorValue(q, i, k);
auto u_ik = E.VectorComponentValue(q, i, k);
// Should this be rotated?
Tensor F_ik = E.VectorGradient(q, i, k);
Tensor FA_ik = F_ik * inverseFA;
auto F_ik = E.VectorRowGradient(q, i, k);
auto FA_ik = F_ik * inverseFA;
// Acceleration
r_c(i, k) += w * rho * (a * u_ik);
......@@ -178,7 +177,7 @@ void LagrangeElasticity::PressureBoundary(const cell &c, int face, int bc, const
for (int i = 0; i < E.size(); ++i) {
for (int k = 0; k < c.dim(); ++k) {
VectorField u_ik = E.VectorValue(q, i, k);
auto u_ik = E.VectorComponentValue(q, i, k);
r_c(i, k) -= w * (localPressure) * N * u_ik;
}
}
......@@ -203,7 +202,7 @@ void LagrangeElasticity::TractionBoundary(const cell &c, int face, int bc, const
for (int i = 0; i < E.size(); ++i) {
for (int k = 0; k < c.dim(); ++k) {
VectorField u_ik = E.VectorValue(q, i, k);
auto u_ik = E.VectorComponentValue(q, i, k);
r_c(i, k) += w * t * u_ik;
}
......@@ -255,14 +254,14 @@ void LagrangeElasticity::Jacobi(const cell &c, const Vector &U, Matrix &A) const
for (int i = 0; i < E.size(); ++i) {
for (int k = 0; k < c.dim(); ++k) {
VectorField u_ik = E.VectorValue(q, i, k);
Tensor F_ik = E.VectorGradient(q, i, k);
Tensor FA_ik = F_ik * inverseFA;
auto u_ik = E.VectorComponentValue(q, i, k);
auto F_ik = E.VectorRowGradient(q, i, k);
auto FA_ik = F_ik * inverseFA;
for (int j = 0; j < E.size(); ++j) {
for (int l = 0; l < c.dim(); ++l) {
VectorField u_jl = E.VectorValue(q, j, l);
Tensor F_jl = E.VectorGradient(q, j, l);
Tensor FA_jl = F_jl * inverseFA;
auto u_jl = E.VectorComponentValue(q, j, l);
auto F_jl = E.VectorRowGradient(q, j, l);
auto FA_jl = F_jl * inverseFA;
// Acceleration
A_c(i, j, k, l) += w * rho * timeScheme->Jabobi() * (u_ik * u_jl);
......@@ -518,7 +517,7 @@ void LagrangeElasticity::InitializePrestress(const Vector &U) {
for (int i = 0; i < E.size(); ++i) {
for (int k = 0; k < c.dim(); ++k) {
Tensor F_ik = E.VectorGradient(q, i, k);
auto F_ik = E.VectorRowGradient(q, i, k);
// Passive material response
r_c(i, k) += w * cellMat.IsoDerivative(Fiso, F_ik);
......
......@@ -40,12 +40,6 @@ GuccioneMaterial::FirstPiolaKirchhoff(const Tensor &F, const Tensor &Fiso, const
return F * firstPK(E);
}
std::function<double(const Tensor &)>
GuccioneMaterial::ConstitutiveMatrix(const Tensor &F, const Tensor &Q, const Tensor &H) const {
using namespace std::placeholders;
return std::bind(&GuccioneMaterial::IsotropicSecondDerivative, this, F, H, _1);
}
double GuccioneMaterial::IsotropicEnergy(const Tensor &F) const {
return 0.0;
}
......@@ -103,6 +97,17 @@ GuccioneMaterial::AnisotropicSecondDerivative(const Tensor &F, const Tensor &Q,
secondDerivative(E, sym(FT * HQ), sym(FT * GQ));
}
double GuccioneMaterial::IsotropicSecondDerivative(const Tensor &F, const TensorRow &H,
const TensorRow &G) const {
return 0.0;
}
double
GuccioneMaterial::AnisotropicSecondDerivative(const Tensor &F, const Tensor &Q, const TensorRow &H,
const TensorRow &G) const {
return AnisotropicSecondDerivative(F, Q, Tensor(H), Tensor(G));
}
/*
quad GuccioneBonetMaterial::expQ(const Tensor &E, const VectorField &f) const {
......@@ -168,12 +173,6 @@ Tensor GuccioneBonetMaterial::FirstPiolaKirchhoff(const Tensor &F, const Tensor
return F * (firstPK(E, f));
}
std::function<double(const Tensor &)>
GuccioneBonetMaterial::ConstitutiveMatrix(const Tensor &F, const Tensor &Q, const Tensor &H) const {
using namespace std::placeholders;
return std::bind(&GuccioneBonetMaterial::IsotropicSecondDerivative, this, F, H, _1);
}
double GuccioneBonetMaterial::IsotropicEnergy(const Tensor &F) const {
return 0.0;
}
......@@ -189,6 +188,18 @@ GuccioneBonetMaterial::IsotropicSecondDerivative(const Tensor &F, const Tensor &
return 0.0;
}
double
GuccioneBonetMaterial::IsotropicDerivative(const Tensor &F, const TensorRow &H) const {
return 0.0;
}
double
GuccioneBonetMaterial::IsotropicSecondDerivative(const Tensor &F, const TensorRow &H,
const TensorRow &G) const {
return 0.0;
}
string GuccioneBonetMaterial::Name() const {
return "GuccioneBonet";
}
......@@ -218,3 +229,25 @@ double GuccioneBonetMaterial::AnisotropicSecondDerivative(const Tensor &F, const
return Frobenius(firstPK(E, Q[0]), sym(transpose(G) * H)) +
Frobenius(constitutiveMatrix(E, Q[0], sym(FT * H)), sym(FT * G));
}
double GuccioneBonetMaterial::AnisotropicDerivative(const Tensor &F, const Tensor &Q,
const TensorRow &H) const {
Tensor FT = transpose(F);
Tensor E = 0.5 * (FT * F - One);
return Frobenius(F * firstPK(E, Q[0]), H);
}
double GuccioneBonetMaterial::AnisotropicSecondDerivative(const Tensor &F, const Tensor &Q,
const TensorRow &H, const TensorRow &G) const {
// Using numerical derivative
/*using namespace ::std::placeholders;
auto derivativeFunc = std::bind(&GuccioneBonetMaterial::AnisotropicDerivative, this, _1, Q, G);
return numericalDerivative(derivativeFunc, F, H);*/
Tensor FT = transpose(F);
Tensor E = 0.5 * (FT * F - One);
return Frobenius(firstPK(E, Q[0]), sym(ColRowProduct(G, H))) +
Frobenius(constitutiveMatrix(E, Q[0], sym(FT * H)), sym(FT * G));
}
......@@ -22,8 +22,8 @@ class GuccioneMaterial : public HyperelasticMaterial {
double bfs() const { return GetParameter(3); }
std::pair<size_t, size_t> GetMaxParameterRange() const {
// start, length
return {0, 1};
// start, length
return {0, 1};
}
Tensor fiberWeights(const Tensor &E) const {
......@@ -54,15 +54,14 @@ public:
Tensor FirstPiolaKirchhoff(const Tensor &F, const Tensor &Fiso, const Tensor &Q) const override;
std::function<double(const Tensor &)>
ConstitutiveMatrix(const Tensor &F, const Tensor &Q, const Tensor &H) const override;
double IsotropicEnergy(const Tensor &F) const override;
double IsotropicDerivative(const Tensor &F, const Tensor &H) const override;
double
IsotropicSecondDerivative(const Tensor &F, const Tensor &H, const Tensor &G) const override;
double
IsotropicSecondDerivative(const Tensor &F, const TensorRow &H, const TensorRow &G) const override;
double AnisotropicEnergy(const Tensor &F, const Tensor &Q) const override;
......@@ -71,6 +70,9 @@ public:
double AnisotropicSecondDerivative(const Tensor &F, const Tensor &Q, const Tensor &H,
const Tensor &G) const override;
double AnisotropicSecondDerivative(const Tensor &F, const Tensor &Q, const TensorRow &H,
const TensorRow &G) const override;
int ParameterCount() const override { return 4; }
string Name() const override;
......@@ -105,8 +107,6 @@ class GuccioneBonetMaterial : public HyperelasticMaterial {
double bfs() const { return GetParameter(3); }
double q(const Tensor &E, const VectorField &f) const;
Tensor dq(const Tensor &E, const VectorField &f) const;
......@@ -141,9 +141,6 @@ public:
Tensor FirstPiolaKirchhoff(const Tensor &E, const Tensor &Fiso, const Tensor &Q) const override;
std::function<double(const Tensor &)>
ConstitutiveMatrix(const Tensor &E, const Tensor &Q, const Tensor &H) const override;
double IsotropicEnergy(const Tensor &F) const override;
double IsotropicDerivative(const Tensor &F, const Tensor &H) const override;
......@@ -151,6 +148,11 @@ public:
double
IsotropicSecondDerivative(const Tensor &F, const Tensor &H, const Tensor &G) const override;
double IsotropicDerivative(const Tensor &F, const TensorRow &H) const override;
double
IsotropicSecondDerivative(const Tensor &F, const TensorRow &H, const TensorRow &G) const override;
double AnisotropicEnergy(const Tensor &F, const Tensor &Q) const override;
double AnisotropicDerivative(const Tensor &F, const Tensor &Q, const Tensor &H) const override;
......@@ -158,6 +160,13 @@ public:
double AnisotropicSecondDerivative(const Tensor &F, const Tensor &Q, const Tensor &H,
const Tensor &G) const override;
double AnisotropicDerivative(const Tensor &F, const Tensor &Q, const TensorRow &H) const override;
double AnisotropicSecondDerivative(const Tensor &F, const Tensor &Q, const TensorRow &H,
const TensorRow &G) const override;
int ParameterCount() const override { return 4; }
string Name() const override;
......
......@@ -70,12 +70,6 @@ HolzapfelMaterial::FirstPiolaKirchhoff(const Tensor &F, const Tensor &Fiso, cons
}
std::function<double(const Tensor &)>
HolzapfelMaterial::ConstitutiveMatrix(const Tensor &F, const Tensor &Q, const Tensor &H) const {
using namespace std::placeholders;
return std::bind(&HolzapfelMaterial::IsotropicSecondDerivative, this, F, H, _1);
}
double HolzapfelMaterial::IsotropicEnergy(const Tensor &F) const {
Tensor C = transpose(F) * F;
return energyIso(C);
......@@ -123,3 +117,33 @@ string HolzapfelMaterial::Name() const {
return "Ogden-Holzapfel";
}
double HolzapfelMaterial::IsotropicDerivative(const Tensor &F, const TensorRow &H) const {
Tensor C = transpose(F) * F;
return 2 * Frobenius(F * firstPKIso(C), H);
}
double HolzapfelMaterial::IsotropicSecondDerivative(const Tensor &F, const TensorRow &H,
const TensorRow &G) const {
Tensor FT = transpose(F);
Tensor C = FT * F;
return 2 * Frobenius(firstPKIso(C), sym(ColRowProduct(G, H))) +
4 * Frobenius(constitutiveMatrixIso(C, sym(FT*H)), sym(FT * G));
}
double HolzapfelMaterial::AnisotropicDerivative(const Tensor &F, const Tensor &Q,
const TensorRow &H) const {
Tensor C = transpose(F) * F;
return 2 * Frobenius(F * firstPKAniso(C, Q), H);
}
double
HolzapfelMaterial::AnisotropicSecondDerivative(const Tensor &F, const Tensor &Q, const TensorRow &H,
const TensorRow &G) const {
Tensor FT = transpose(F);
Tensor C = FT * F;
return 2 * Frobenius(firstPKAniso(C, Q), sym(ColRowProduct(G, H))) +
4 * Frobenius(constitutiveMatrixAniso(C, Q, sym(FT*H)), sym(FT * G));
}
......@@ -67,9 +67,6 @@ public:
Tensor FirstPiolaKirchhoff(const Tensor &F, const Tensor &Fiso, const Tensor &Q) const override;
std::function<double(const Tensor &)>
ConstitutiveMatrix(const Tensor &F, const Tensor &Q, const Tensor &H) const override;
double IsotropicEnergy(const Tensor &F) const override;
double IsotropicDerivative(const Tensor &F, const Tensor &H) const override;
......@@ -77,6 +74,7 @@ public:
double
IsotropicSecondDerivative(const Tensor &F, const Tensor &H, const Tensor &G) const override;
double AnisotropicEnergy(const Tensor &F, const Tensor &Q) const override;
double AnisotropicDerivative(const Tensor &F, const Tensor &Q, const Tensor &H) const override;
......@@ -85,6 +83,17 @@ public:
AnisotropicSecondDerivative(const Tensor &F, const Tensor &Q, const Tensor &H,
const Tensor &G) const override;
double IsotropicDerivative(const Tensor &F, const TensorRow &H) const override;
double IsotropicSecondDerivative(const Tensor &F, const TensorRow &H, const TensorRow &G) const override;
double AnisotropicDerivative(const Tensor &F, const Tensor &Q, const TensorRow &H) const override;
double AnisotropicSecondDerivative(const Tensor &F, const Tensor &Q, const TensorRow &H, const TensorRow &G) const override;
int ParameterCount() const override { return 8; }
string Name() const override;
......
......@@ -15,18 +15,17 @@ LaplaceMaterial::constitutiveMatrix(const Tensor &Du, const Tensor &H) const {
return permeability * H;
}
Tensor
LaplaceMaterial::constitutiveMatrix(const Tensor &Du, const TensorRow &H) const {
return permeability * H;
}
Tensor
LaplaceMaterial::FirstPiolaKirchhoff(const Tensor &F, const Tensor &Fiso, const Tensor &Q) const {
return firstPK(Fiso - One);
}
std::function<double(const Tensor &)>
LaplaceMaterial::ConstitutiveMatrix(const Tensor &F, const Tensor &Q, const Tensor &H) const {
using namespace std::placeholders;
return std::bind(&LaplaceMaterial::IsotropicSecondDerivative, this, F, H, _1);
}
double LaplaceMaterial::IsotropicEnergy(const Tensor &F) const {
return energy(F - One);
}
......@@ -45,4 +44,13 @@ string LaplaceMaterial::Name() const {
return "Laplace";
}
double LaplaceMaterial::IsotropicDerivative(const Tensor &F, const TensorRow &H) const {
return Frobenius(firstPK(F - One), H);
}
double LaplaceMaterial::IsotropicSecondDerivative(const Tensor &F, const TensorRow &H,
const TensorRow &G) const {
return Frobenius(constitutiveMatrix(F - One, H), G);
}
......@@ -11,6 +11,7 @@ class LaplaceMaterial : public HyperelasticMaterial {
Tensor firstPK(const Tensor &Du) const;
Tensor constitutiveMatrix(const Tensor &Du, const Tensor &H) const;
Tensor constitutiveMatrix(const Tensor &Du, const TensorRow &H) const;
public:
LaplaceMaterial() noexcept: HyperelasticMaterial(1) {
......@@ -25,9 +26,6 @@ public:
Tensor FirstPiolaKirchhoff(const Tensor &F, const Tensor &Fiso, const Tensor &Q) const override;
std::function<double(const Tensor &)>
ConstitutiveMatrix(const Tensor &F, const Tensor &Q, const Tensor &H) const override;
double IsotropicEnergy(const Tensor &F) const override;
double IsotropicDerivative(const Tensor &F, const Tensor &H) const override;
......@@ -35,9 +33,15 @@ public:
double
IsotropicSecondDerivative(const Tensor &F, const Tensor &H, const Tensor &G) const override;
double IsotropicDerivative(const Tensor &F, const TensorRow &H) const override;
double
IsotropicSecondDerivative(const Tensor &F, const TensorRow &H, const TensorRow &G) const override;
int ParameterCount() const override { return 1; }
string Name() const override;
};
#endif //LAPLACEMATERIAL_HPP
......@@ -21,12 +21,6 @@ LinearMaterial::FirstPiolaKirchhoff(const Tensor &F, const Tensor &Fiso, const T
return firstPK(sym(Fiso - One));
}
std::function<double(const Tensor &)>
LinearMaterial::ConstitutiveMatrix(const Tensor &F, const Tensor &Q, const Tensor &H) const {
using namespace std::placeholders;
return std::bind(&LinearMaterial::IsotropicSecondDerivative, this, F, H, _1);
}
double LinearMaterial::IsotropicEnergy(const Tensor &F) const {
return energy(sym(F - One));
}
......@@ -43,3 +37,12 @@ LinearMaterial::IsotropicSecondDerivative(const Tensor &F, const Tensor &H, cons
string LinearMaterial::Name() const {
return "Linearized";
}
double LinearMaterial::IsotropicDerivative(const Tensor &F, const TensorRow &H) const {
return Frobenius(firstPK(F - One), sym(H));
}
double LinearMaterial::IsotropicSecondDerivative(const Tensor &F, const TensorRow &H,
const TensorRow &G) const {
return Frobenius(constitutiveMatrix(sym(F - One), sym(H)), sym(G));
}
......@@ -35,9 +35,6 @@ public:
Tensor FirstPiolaKirchhoff(const Tensor &F, const Tensor &Fiso, const Tensor &Q) const override;
std::function<double(const Tensor &)>
ConstitutiveMatrix(const Tensor &F, const Tensor &Q, const Tensor &H) const override;
double IsotropicEnergy(const Tensor &F) const override;
double IsotropicDerivative(const Tensor &F, const Tensor &H) const override;
......@@ -45,6 +42,11 @@ public:
double
IsotropicSecondDerivative(const Tensor &F, const Tensor &H, const Tensor &G) const override;
double IsotropicDerivative(const Tensor &F, const TensorRow &H) const override;
double
IsotropicSecondDerivative(const Tensor &F, const TensorRow &H, const TensorRow &G) const override;
int ParameterCount() const override { return 2; }
string Name() const override;
......
......@@ -33,6 +33,9 @@ public:
* @param F = (Du + One) for the derivative Du of a given displacement u.
*/
virtual double IsotropicDerivative(const Tensor &F, const Tensor &H) const = 0;
virtual double IsotropicDerivative(const Tensor &F, const TensorRow &H) const {
return IsotropicDerivative(F, Tensor(H));
}
/**
* Returns the second strain derivative of the material given the deformation Gradient F
......@@ -41,6 +44,10 @@ public:
virtual double
IsotropicSecondDerivative(const Tensor &F, const Tensor &H, const Tensor &G) const = 0;
virtual double IsotropicSecondDerivative(const Tensor &F, const TensorRow &H, const TensorRow &G) const{
return IsotropicSecondDerivative(F, Tensor(H), Tensor(G));
}
/**
* Returns the strain energy of the material given the deformation Gradient F
......@@ -53,6 +60,7 @@ public:
* @param F = (Du + One) for the derivative Du of a given displacement u.
*/
virtual double AnisotropicDerivative(const Tensor &F, const Tensor &Q, const Tensor &H) const {return 0.0; };
virtual double AnisotropicDerivative(const Tensor &F, const Tensor &Q, const TensorRow &H) const{return AnisotropicDerivative(F, Q,Tensor(H));}
/**
* Returns the second strain derivative of the material given the deformation Gradient F
......@@ -61,6 +69,7 @@ public:
virtual double
AnisotropicSecondDerivative(const Tensor &F, const Tensor &Q, const Tensor &H,
const Tensor &G) const {return 0.0;};
virtual double AnisotropicSecondDerivative(const Tensor &F, const Tensor &Q, const TensorRow &H, const TensorRow &G) const{return AnisotropicSecondDerivative(F, Q, Tensor(H), Tensor(G));}
/**
* Returns the derivative of the strain energy given the appropriate tensor.
......@@ -70,16 +79,6 @@ public:
virtual Tensor
FirstPiolaKirchhoff(const Tensor &F, const Tensor &Fiso, const Tensor &Q) const = 0;
/**
* Returns the second derivative of the strain energy given the appropriate tensor
* The derivative of the first Piola Kirchhoff tensor along H is computed.
* @param F = (Du + One) for the derivative Du of a given displacement u.
* @param H : Directional derivative, usually a test function in the given FE space.
* @param Q : Rotation matrix containing fibre directions, i.e. Q[0]=f etc.
*/
virtual std::function<double(const Tensor &)>
ConstitutiveMatrix(const Tensor &F, const Tensor &Q, const Tensor &H) const = 0;
virtual int ParameterCount() const = 0;
double GetParameter(int i) const { return parameters[i]; }
......@@ -148,12 +147,14 @@ public:
return hyperelMat->IsotropicEnergy(F);
}
double IsoDerivative(const Tensor &F, const Tensor &H) const {
template<typename T>
double IsoDerivative(const Tensor &F, const T &H) const {
return hyperelMat->IsotropicDerivative(F, H);
}
template<typename T>
double
IsoSecondDerivative(const Tensor &F, const Tensor &H, const Tensor &G) const {
IsoSecondDerivative(const Tensor &F, const T &H, const T &G) const {
return hyperelMat->IsotropicSecondDerivative(F, H, G);
}
......@@ -162,12 +163,14 @@ public:
return hyperelMat->AnisotropicEnergy(F, Q);
}
double AnisoDerivative(const Tensor &F, const Tensor &Q, const Tensor &H) const {
template<typename T>
double AnisoDerivative(const Tensor &F, const Tensor &Q, const T &H) const {
return hyperelMat->AnisotropicDerivative(F, Q, H);
}
template<typename T>
double
AnisoSecondDerivative(const Tensor &F, const Tensor &Q, const Tensor &H, const Tensor &G) const {
AnisoSecondDerivative(const Tensor &F, const Tensor &Q, const T &H, const T &G) const {
return hyperelMat->AnisotropicSecondDerivative(F, Q, H, G);
}
......@@ -175,11 +178,6 @@ public:
return hyperelMat->FirstPiolaKirchhoff(F, Fiso, Q);
}
std::function<double(const Tensor &)>
IsoConstitutive(const Tensor &F, const Tensor &Q, const Tensor &H) const {
return hyperelMat->ConstitutiveMatrix(F, Q, H);
}
double Penalty() const {
return volPenalty->Penalty();
}
......@@ -190,8 +188,11 @@ public:
double VolDerivative(const Tensor &F, const Tensor &H) const;
double
VolSecondDerivative(const Tensor &F, const Tensor &H, const Tensor &G) const;
double VolSecondDerivative(const Tensor &F, const Tensor &H, const Tensor &G) const;
double VolDerivative(const Tensor &F, const TensorRow &H) const {return VolDerivative(F, Tensor(H));}
double VolSecondDerivative(const Tensor &F, const TensorRow &H, const TensorRow &G) const{return VolSecondDerivative(F, Tensor(H), Tensor(G));};
Tensor VolFPK(double J) const {
if (volPenalty->Penalty() > 0) {
......@@ -214,11 +215,13 @@ public:
return Zero;
}
double TotalDerivative(const Tensor &Fiso, const Tensor &F, const Tensor &Q, const Tensor &H) const{
template<typename T>
double TotalDerivative(const Tensor &Fiso, const Tensor &F, const Tensor &Q, const T &H) const{
return IsoDerivative(Fiso, H) + AnisoDerivative(F, Q, H) + Penalty() * VolDerivative(F, H);
}
double TotalSecondDerivative(const Tensor &Fiso, const Tensor &F, const Tensor &Q, const Tensor &H, const Tensor &G) const{
template<typename T>
double TotalSecondDerivative(const Tensor &Fiso, const Tensor &F, const Tensor &Q, const T &H, const T &G) const{
return
IsoSecondDerivative(Fiso, H, G)
+ AnisoSecondDerivat