diff --git a/cardmech/src/elasticity/assemble/LagrangeElasticity.cpp b/cardmech/src/elasticity/assemble/LagrangeElasticity.cpp
index c7e2c19ce00928ef11645e9ca542402cb2c34028..8b302a2ede4f7fa617625082e47d6a22d0e03544 100644
--- a/cardmech/src/elasticity/assemble/LagrangeElasticity.cpp
+++ b/cardmech/src/elasticity/assemble/LagrangeElasticity.cpp
@@ -108,9 +108,7 @@ void LagrangeElasticity::Residual(const cell &c, const Vector &U, Vector &r) con
VectorField u = E.VectorValue(q, U);
Tensor F = DeformationGradient(E, q, U);
- auto flory = eProblem.FlorySplit(F);
- auto Fiso = flory.first;
- auto J = flory.second;
+ auto Fiso = eProblem.IsochoricPart(F);
Tensor inverseFA = mat::active_deformation(Q, E.VectorValue(q, *gamma));
F = F * inverseFA;
@@ -132,7 +130,7 @@ void LagrangeElasticity::Residual(const cell &c, const Vector &U, Vector &r) con
r_c(i, k) += w * cellMat.IsoDerivative(Fiso, F_ik);
r_c(i, k) += w * cellMat.AnisoDerivative(F, Q, F_ik);
r_c(i, k) +=
- w * mat::penalty(cellMat) * cellMat.VolDerivative(J, F_ik /* *transpose(inverseFA) */);
+ w * cellMat.Penalty() * cellMat.VolDerivative(F, F_ik /* *transpose(inverseFA) */);
// Body load
r_c(i, k) -= w * (f * u_ik);
@@ -142,8 +140,8 @@ void LagrangeElasticity::Residual(const cell &c, const Vector &U, Vector &r) con
// Prestress
r_c(i, k) += w * cellMat.IsoDerivative(F_0, F_ik);
r_c(i, k) +=
- w * mat::penalty(cellMat) *
- cellMat.VolDerivative(det(F_0), F_ik /* *transpose(inverseFA) */);
+ w * cellMat.Penalty() *
+ cellMat.VolDerivative(F_0, F_ik /* *transpose(inverseFA) */);
}
}
}
@@ -248,9 +246,7 @@ void LagrangeElasticity::Jacobi(const cell &c, const Vector &U, Matrix &A) const
double w = E.QWeight(q);
Tensor F = DeformationGradient(E, q, U);
- auto flory = eProblem.FlorySplit(F);
- auto Fiso = flory.first;
- auto J = flory.second;
+ auto Fiso = eProblem.IsochoricPart(F);
Tensor inverseFA = mat::active_deformation(Q, E.VectorValue(q, *gamma));
F = F * inverseFA;
@@ -273,7 +269,7 @@ void LagrangeElasticity::Jacobi(const cell &c, const Vector &U, Matrix &A) const
A_c(i, j, k, l) += w * cellMat.IsoSecondDerivative(Fiso, F_ik, F_jl);
A_c(i, j, k, l) += w * cellMat.AnisoSecondDerivative(F, Q, F_ik, F_jl);
A_c(i, j, k, l) +=
- w * mat::penalty(cellMat) * cellMat.VolSecondDerivative(J, F_ik, F_jl);
+ w * cellMat.Penalty() * cellMat.VolSecondDerivative(F, F_ik, F_jl);
// Active stress
//A_c(i,k) += w * JacobiStress(T, QT[0], DeformationGradient(Du), E.VectorGradient(q, i, k), E.VectorGradient(q, j, l));
diff --git a/cardmech/src/elasticity/materials/Material.cpp b/cardmech/src/elasticity/materials/Material.cpp
index aa62dcb26675d9c5b82abec28446ef7abacd56b5..e022c56f693606381db731ae96db2fd189ad00fa 100644
--- a/cardmech/src/elasticity/materials/Material.cpp
+++ b/cardmech/src/elasticity/materials/Material.cpp
@@ -66,3 +66,20 @@ void Material::initVolumetricPenalty() {
initVolumetricPenalty(penaltyType);
}
+
+double Material::VolDerivative(const Tensor &F, const Tensor &H) const {
+ if (volPenalty->Penalty() == 0.0) return 0.0;
+
+ double dwj = volPenalty->Derivative(det(F));
+ return dwj != 0.0 ? dwj * Frobenius(mat::cofactor(F), H) : 0.0;
+}
+
+double Material::VolSecondDerivative(const Tensor &F, const Tensor &H, const Tensor &G) const {
+ if (volPenalty->Penalty() == 0.0) return 0.0;
+
+ double J = det(F);
+ Tensor CofF = mat::cofactor(F);
+ double ddwj = volPenalty->SecondDerivative(J) * Frobenius(CofF, H) * Frobenius(CofF, G);
+ ddwj += volPenalty->Derivative(J) * Frobenius(mat::d_cofactor(F, H), G);
+ return ddwj;
+}
diff --git a/cardmech/src/elasticity/materials/Material.hpp b/cardmech/src/elasticity/materials/Material.hpp
index aa56b4128fd4532d7bce5b45ea28f8a62f0639e2..3cbdb93933b6f7b9079f5b333bdd4595bfaeaea2 100644
--- a/cardmech/src/elasticity/materials/Material.hpp
+++ b/cardmech/src/elasticity/materials/Material.hpp
@@ -3,6 +3,7 @@
#include "utility/Config.hpp"
#include "VolumetricPenalty.hpp"
+#include "materials.hpp"
#include <Tensor.hpp>
#include <functional>
@@ -161,32 +162,18 @@ public:
return hyperelMat->ConstitutiveMatrix(F, Q, H);
}
-
- double VolEnergy(double J) const {
- return volPenalty->Penalty() > 0 ? volPenalty->Functional(J) : 0.0;
+ double Penalty() const {
+ return volPenalty->Penalty();
}
- double VolDerivative(double J, const Tensor &H) const {
- if (volPenalty->Penalty() > 0) {
- double dwj = volPenalty->Functional(J) * volPenalty->Derivative(J) * J;
- return dwj != 0 ? dwj * trace(H) : 0.0;
- }
- return 0.0;
+ double VolEnergy(const Tensor &F) const {
+ return volPenalty->Penalty() > 0 ? volPenalty->Functional(det(F)) : 0.0;
}
+ double VolDerivative(const Tensor &F, const Tensor &H) const;
+
double
- VolSecondDerivative(double J, const Tensor &H, const Tensor &G) const {
- if (volPenalty->Penalty() > 0) {
- double volPart = 0.0;
- double dwj = volPenalty->Derivative(J) * J;
- double ddwj = volPenalty->SecondDerivative(J) * J;
- if (dwj != 0.0 || ddwj != 0.0)
- volPart = (ddwj + dwj) * trace(H) * trace(G)
- - dwj * Frobenius(G, H); // = trace((FEinv * H * FEinv * G)^T)
- return volPart;
- }
- return 0.0;
- }
+ VolSecondDerivative(const Tensor &F, const Tensor &H, const Tensor &G) const;
Tensor VolFPK(double J) const {
if (volPenalty->Penalty() > 0) {
@@ -216,104 +203,6 @@ public:
};
-namespace mat {
-
-
- static Tensor cofactor(const Tensor &F) {
- return 0.5 * Cross(F, F);
- }
-
- static Tensor d_cofactor(const Tensor &F, const Tensor &Du) {
- return Cross(F, Du);
- }
-
- static Tensor dd_cofactor(const Tensor &F, const Tensor &Du, const Tensor &Dv) {
- return Cross(Du, Dv);
- }
-
- static double d_det(const Tensor &F, const Tensor &Du) {
- return Frobenius(cofactor(F), Du);
- }
-
- static double dd_det(const Tensor &F, const Tensor &Du, const Tensor &Dv) {
- return Frobenius(F, Cross(Du, Dv));
- }
-
- static double macaulay(double value) {
- return (value > 0) * value;
- }
-
- static double
- invariant(int k, const Tensor &F, const VectorField &f = zero, const VectorField &l = zero) {
- switch (k) {
- case 1:
- return Frobenius(F, F);
- case 2:
- return Frobenius(cofactor(F), cofactor(F));
- case 3:
- return det(F);
- case 4:
- return (F * f) * (F * f);
- case 8:
- return abs((F * f) * (F * l));
- default:
- return 1.0;
- }
- }
-
- static Tensor
- dinvariant(int k, const Tensor &F, const VectorField &f = zero, const VectorField &l = zero) {
- switch (k) {
- case 1:
- return 2.0 * F;
- case 2:
- return invariant(1, F) * dinvariant(1, F) - 4.0 * F * transpose(F) * F;
- case 3:
- return cofactor(F);
- case 4:
- return 2.0 * F * Product(f, f);
- case 8:
- return F * (Product(f, l) + Product(l, f));
- default:
- return 1.0;
- }
- }
-
- static Tensor
- ddinvariant(int k, const Tensor &F, const Tensor &H, const VectorField &f = zero,
- const VectorField &l = zero) {
- switch (k) {
- case 1:
- return 2.0 * H;
- case 2:
- return dinvariant(1, F) * dinvariant(1, F) + invariant(1, F) * ddinvariant(1, F, H)
- - 4.0 * (F * transpose(F) * H + F * transpose(H) * F + H * transpose(F) * F);
- case 3:
- // https://math.stackexchange.com/questions/2136747/derivate-of-the-cofactor-and-the-determinant
- return transposeInvert(F) * (Frobenius(cofactor(F), H) * One - transpose(H) * cofactor(F));
- case 4:
- return 2.0 * H * Product(f, f);
- case 8:
- return H * (Product(f, l) + Product(l, f));
- default:
- return 1.0;
- }
- }
-
- static Tensor active_deformation(const Tensor &Q, const VectorField &gamma) {
- if (gamma[0] == 0.0)
- return One;
-
- Tensor AD = One + gamma[0] * Product(Q[0], Q[0])
- + gamma[1] * Product(Q[1], Q[1])
- + gamma[2] * Product(Q[2], Q[2]);
- return AD;
- }
-
- static double penalty(const Material &material) {
- return material.Volumetric().Penalty();
- }
-}
#endif
diff --git a/cardmech/src/elasticity/materials/VolumetricPenalty.hpp b/cardmech/src/elasticity/materials/VolumetricPenalty.hpp
index f2d26d089358f81cd1553b27a9097fb5573042c6..8daf183e6965c78179d832d23eecf2888fc63282 100644
--- a/cardmech/src/elasticity/materials/VolumetricPenalty.hpp
+++ b/cardmech/src/elasticity/materials/VolumetricPenalty.hpp
@@ -34,7 +34,7 @@ public:
};
class CiarletPenalty : public VolumetricPenalty {
- double mu{1.0};
+ double mu{0.0};
public:
CiarletPenalty(bool mixed = false) : VolumetricPenalty() {
if (mixed && penalty > 0.0) {
diff --git a/cardmech/src/elasticity/materials/materials.hpp b/cardmech/src/elasticity/materials/materials.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..e820095c56604b08d24098ff92eb77c58209be63
--- /dev/null
+++ b/cardmech/src/elasticity/materials/materials.hpp
@@ -0,0 +1,100 @@
+#ifndef MATERIALS_HPP
+#define MATERIALS_HPP
+
+
+#include "Tensor.hpp"
+
+namespace mat {
+ static Tensor cofactor(const Tensor &F) {
+ return 0.5 * Cross(F, F);
+ }
+
+ static Tensor d_cofactor(const Tensor &F, const Tensor &Du) {
+ return Cross(F, Du);
+ }
+
+ static Tensor dd_cofactor(const Tensor &F, const Tensor &Du, const Tensor &Dv) {
+ return Cross(Du, Dv);
+ }
+
+ static double d_det(const Tensor &F, const Tensor &Du) {
+ return Frobenius(cofactor(F), Du);
+ }
+
+ static double dd_det(const Tensor &F, const Tensor &Du, const Tensor &Dv) {
+ return Frobenius(F, Cross(Du, Dv));
+ }
+
+ static double macaulay(double value) {
+ return (value > 0) * value;
+ }
+
+ static double
+ invariant(int k, const Tensor &F, const VectorField &f = zero, const VectorField &l = zero) {
+ switch (k) {
+ case 1:
+ return Frobenius(F, F);
+ case 2:
+ return Frobenius(cofactor(F), cofactor(F));
+ case 3:
+ return det(F);
+ case 4:
+ return (F * f) * (F * f);
+ case 8:
+ return abs((F * f) * (F * l));
+ default:
+ return 1.0;
+ }
+ }
+
+ static Tensor
+ dinvariant(int k, const Tensor &F, const VectorField &f = zero, const VectorField &l = zero) {
+ switch (k) {
+ case 1:
+ return 2.0 * F;
+ case 2:
+ return invariant(1, F) * dinvariant(1, F) - 4.0 * F * transpose(F) * F;
+ case 3:
+ return cofactor(F);
+ case 4:
+ return 2.0 * F * Product(f, f);
+ case 8:
+ return F * (Product(f, l) + Product(l, f));
+ default:
+ return 1.0;
+ }
+ }
+
+ static Tensor
+ ddinvariant(int k, const Tensor &F, const Tensor &H, const VectorField &f = zero,
+ const VectorField &l = zero) {
+ switch (k) {
+ case 1:
+ return 2.0 * H;
+ case 2:
+ return dinvariant(1, F) * dinvariant(1, F) + invariant(1, F) * ddinvariant(1, F, H)
+ - 4.0 * (F * transpose(F) * H + F * transpose(H) * F + H * transpose(F) * F);
+ case 3:
+ // https://math.stackexchange.com/questions/2136747/derivate-of-the-cofactor-and-the-determinant
+ return transposeInvert(F) * (Frobenius(cofactor(F), H) * One - transpose(H) * cofactor(F));
+ case 4:
+ return 2.0 * H * Product(f, f);
+ case 8:
+ return H * (Product(f, l) + Product(l, f));
+ default:
+ return 1.0;
+ }
+ }
+
+ static Tensor active_deformation(const Tensor &Q, const VectorField &gamma) {
+ if (gamma[0] == 0.0)
+ return One;
+
+ Tensor AD = One + gamma[0] * Product(Q[0], Q[0])
+ + gamma[1] * Product(Q[1], Q[1])
+ + gamma[2] * Product(Q[2], Q[2]);
+ return AD;
+ }
+}
+
+#endif //MATERIALS_HPP
diff --git a/cardmech/src/elasticity/problem/ElasticityProblem.hpp b/cardmech/src/elasticity/problem/ElasticityProblem.hpp
index 6c53019f970471b23395c29ff56b582ff3229044..dc45332b12717c47adf37197eaffbbbde75bfe4f 100644
--- a/cardmech/src/elasticity/problem/ElasticityProblem.hpp
+++ b/cardmech/src/elasticity/problem/ElasticityProblem.hpp
@@ -68,9 +68,8 @@ public:
const VectorField &n) const { return 0.0; };
- std::pair<Tensor, double> FlorySplit(const Tensor &F) const {
- double J = det(F);
- return std::pair<Tensor, double>{useFlorySplit ? pow(J, -1.0 / 3.0) * F : F, J};
+ Tensor IsochoricPart(const Tensor &F) const {
+ return useFlorySplit ? pow(det(F), -1.0 / 3.0) * F : F;
}
/**