Commit 1a3eef03 authored by jonathan.froehlich's avatar jonathan.froehlich
Browse files

Added and tested currently used material methods

parent 1987b989
Pipeline #220278 canceled with stages
......@@ -122,9 +122,9 @@ 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) {
auto u_ik = E.VectorValue(q, i, k);
auto u_ik = E.VectorComponentValue(q, i, k);
// Should this be rotated?
auto F_ik = E.VectorGradient(q, i, k);
auto F_ik = E.VectorRowGradient(q, i, k);
auto FA_ik = F_ik * inverseFA;
// Acceleration
......@@ -177,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;
}
}
......@@ -202,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;
}
......@@ -254,13 +254,13 @@ 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) {
auto u_ik = E.VectorValue(q, i, k);
auto F_ik = E.VectorGradient(q, i, k);
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) {
auto u_jl = E.VectorValue(q, j, l);
auto F_jl = E.VectorGradient(q, j, l);
auto u_jl = E.VectorComponentValue(q, j, l);
auto F_jl = E.VectorRowGradient(q, j, l);
auto FA_jl = F_jl * inverseFA;
// Acceleration
......@@ -484,7 +484,7 @@ void LagrangeElasticity::InitializePrestress(const Vector &U) {
for (int i = 0; i < E.size(); ++i) {
for (int k = 0; k < c.dim(); ++k) {
auto 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;
TensorR 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,7 +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));}
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,7 +43,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));}
virtual double IsotropicSecondDerivative(const Tensor &F, const TensorRow &H, const TensorRow &G) const{
return IsotropicSecondDerivative(F, Tensor(H), Tensor(G));
}
/**
......@@ -74,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]; }
......@@ -183,12 +178,6 @@ public:
return hyperelMat->FirstPiolaKirchhoff(F, Fiso, Q);
}
template<typename T>
std::function<double(const Tensor &)>
IsoConstitutive(const Tensor &F, const Tensor &Q, const T &H) const {
return hyperelMat->ConstitutiveMatrix(F, Q, H);
}
double Penalty() const {
return volPenalty->Penalty();
}
......
......@@ -23,12 +23,6 @@ Tensor MooneyRivlinMaterial::FirstPiolaKirchhoff(const Tensor &F, const Tensor &
return Zero;
}
std::function<double(const Tensor &)>
MooneyRivlinMaterial::ConstitutiveMatrix(const Tensor &F, const Tensor &Q, const Tensor &H) const {
using namespace std::placeholders;
return std::bind(&MooneyRivlinMaterial::IsotropicSecondDerivative, this, F, H, _1);
}
double MooneyRivlinMaterial::IsotropicEnergy(const Tensor &F) const {
Tensor H = mat::cofactor(F);
return alpha() * (Frobenius(F, F) - 3) + beta() * Frobenius(H, H) +
......
......@@ -33,9 +33,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;
......
......@@ -29,8 +29,3 @@ string NeoHookeMaterial::Name() const {
}
std::function<double(const Tensor &)>
NeoHookeMaterial::ConstitutiveMatrix(const Tensor &F, const Tensor &Q, const Tensor &H) const {
using namespace std::placeholders;
return std::bind(&NeoHookeMaterial::IsotropicSecondDerivative, this, F, H, _1);
}
......@@ -25,9 +25,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;
......
......@@ -20,13 +20,6 @@ Tensor StVenantKirchhoffMaterial::FirstPiolaKirchhoff(const Tensor &F, const Ten
return Fiso * firstPK(E);
}
std::function<double(const Tensor &)>
StVenantKirchhoffMaterial::ConstitutiveMatrix(const Tensor &F, const Tensor &Q,
const Tensor &H) const {
using namespace std::placeholders;
return std::bind(&StVenantKirchhoffMaterial::IsotropicSecondDerivative, this, F, H, _1);
}
double StVenantKirchhoffMaterial::IsotropicEnergy(const Tensor &F) const {
Tensor E = 0.5 * (transpose(F) * F - One);
return energy(E);
......
......@@ -30,9 +30,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,18 +35,19 @@ struct MaterialTestParameters {
static double numMatDerivative(const HyperelasticMaterial &mat, const Tensor &F, const Tensor &Q,
const Tensor &H) {
const Tensor &H) {
double h = 1e-8;
double dwi = mat.IsotropicEnergy(F + h * H) + mat.AnisotropicEnergy(F + h * H, Q)
- mat.IsotropicEnergy(F - h * H) - mat.AnisotropicEnergy(F - h * H, Q);
return 0.5 * dwi / h;
}
static double numMatSecondDerivative(const HyperelasticMaterial &mat, const Tensor &F, const Tensor &Q,
const Tensor &H, const Tensor& G) {
static double
numMatSecondDerivative(const HyperelasticMaterial &mat, const Tensor &F, const Tensor &Q,
const Tensor &H, const Tensor &G) {