Merge branch '99-extract-volumetric-function-from-neo-hooke' into 'debug'

Moved penalty from Neo Hooke

Closes #99

See merge request mpp/cardmech!114

cardmech/src/elasticity/materials/Material.cpp

@@ -12,8 +12,12 @@ Material::Material(const string &materialName) {
   // === Create HyperElasticMaterial ==== //
   if (materialName == "Laplace") hyperelMat = std::make_unique<LaplaceMaterial>();
   else if (materialName == "Linear") hyperelMat = std::make_unique<LinearMaterial>();
-  else if (materialName == "NeoHooke") hyperelMat = std::make_unique<NeoHookeMaterial>();
-  else if (materialName == "MooneyRivlin") hyperelMat = std::make_unique<MooneyRivlinMaterial>();
+  else if (materialName == "NeoHooke") {
+    hyperelMat = std::make_unique<NeoHookeMaterial>();
+    volPenalty = std::make_unique<CiarletPenalty>(hyperelMat->GetParameter(0),
+                                                  hyperelMat->GetParameter(1));
+    return;
+  } else if (materialName == "MooneyRivlin") hyperelMat = std::make_unique<MooneyRivlinMaterial>();
   else if (materialName == "Guccione") hyperelMat = std::make_unique<GuccioneMaterial>();
   else if (materialName == "Bonet") hyperelMat = std::make_unique<GuccioneBonetMaterial>();
   else if (materialName == "Holzapfel") hyperelMat = std::make_unique<HolzapfelMaterial>();
@@ -26,9 +30,11 @@ Material::Material(const string &materialName) {
 
 Material::Material(const string &materialName, const vector<double> &materialParameters) {
   if (materialName == "Linear") hyperelMat = std::make_unique<LinearMaterial>(materialParameters);
-  else if (materialName == "NeoHooke")
+  else if (materialName == "NeoHooke") {
     hyperelMat = std::make_unique<NeoHookeMaterial>(materialParameters);
-  else if (materialName == "MooneyRivlin")
+    volPenalty = std::make_unique<CiarletPenalty>(materialParameters[0], materialParameters[1]);
+    return;
+  } else if (materialName == "MooneyRivlin")
     hyperelMat = std::make_unique<MooneyRivlinMaterial>(materialParameters);
   else if (materialName == "Guccione")
     hyperelMat = std::make_unique<GuccioneMaterial>(materialParameters);
@@ -45,10 +51,7 @@ Material::Material(const string &materialName, const vector<double> &materialPar
   initVolumetricPenalty();
 }
 
-void Material::initVolumetricPenalty() {
-  string penaltyType = "None";
-  config.get("QuasiCompressiblePenalty", penaltyType);
-
+void Material::initVolumetricPenalty(const std::string &penaltyType) {
   if (penaltyType == "None") volPenalty = std::make_unique<NoPenalty>();
   else if (penaltyType == "Ciarlet")
     volPenalty = std::make_unique<CiarletPenalty>();
@@ -56,3 +59,10 @@ void Material::initVolumetricPenalty() {
     volPenalty = std::make_unique<LogarithmicCompressiblePenalty>();
   else if (penaltyType == "Mixed") volPenalty = std::make_unique<MixedCompressiblePenalty>();
 }
+
+void Material::initVolumetricPenalty() {
+  string penaltyType = "None";
+  config.get("QuasiCompressiblePenalty", penaltyType);
+  initVolumetricPenalty(penaltyType);
+
+}

cardmech/src/elasticity/materials/Material.hpp

@@ -96,10 +96,7 @@ class Material {
 
   void initVolumetricPenalty();
 
-  std::function<double(const Tensor &, const Tensor &, const Tensor &)> derivativeFunc =
-      [this](const Tensor &a, const Tensor &b, const Tensor &c) {
-        return this->hyperelMat->StrainDerivative(a, b, c);
-      };
+  void initVolumetricPenalty(const string &penaltyType);
 
 public:
   Material() = default;
@@ -182,6 +179,7 @@ public:
   const HyperelasticMaterial &Isochoric() const { return *hyperelMat; }
 
   const VolumetricPenalty &Volumetric() const { return *volPenalty; }
+
 };
 
 namespace mat {

cardmech/src/elasticity/materials/NeoHookeMaterial.cpp

@@ -2,54 +2,31 @@
 
 
 Tensor NeoHookeMaterial::FirstPiolaKirchhoff(const Tensor &F, const Tensor &Q) const {
-  double J = det(F);
-  return mu() * sym(F) + gammaDerivative(J) * mat::cofactor(F);
+  return mu() * sym(F);
 }
 
 Tensor
 NeoHookeMaterial::constitutiveMatrix(const Tensor &F, const Tensor &Q, const Tensor &H) const {
-  double J = det(F);
-  const Tensor FI = Invert(F);
-  const Tensor FIT = transpose(FI);
-  return mu() * sym(H)
-         + gammaSecondDerivative(J) * mat::d_det(F, H) * mat::cofactor(F)
-         + gammaDerivative(J) * mat::d_cofactor(F, H);
+  return FirstPiolaKirchhoff(H, Q);
 }
 
 double NeoHookeMaterial::StrainEnergy(const Tensor &F, const Tensor &Q) const {
-  return mu() * 0.5 * (Frobenius(F, F) - 3) + gammaFunction(det(F));
+  return mu() * 0.5 * (Frobenius(F, F) - 3);
 }
 
 double NeoHookeMaterial::StrainDerivative(const Tensor &F, const Tensor &Q, const Tensor &H) const {
-  double J = det(F);
-  Tensor CofF = J * transposeInvert(F);
-  return mu() * Frobenius(sym(F), H) + gammaDerivative(J) * mat::d_det(F, H);
+  return mu() * Frobenius(sym(F), H);
 }
 
 double NeoHookeMaterial::StrainSecondDerivative(const Tensor &F, const Tensor &Q, const Tensor &H,
                                                 const Tensor &G) const {
-  double J = det(F);
-  Tensor CofF = J * transposeInvert(F);
-  return mu() * Frobenius(sym(H), G)
-         + gammaSecondDerivative(J) * mat::d_det(F, H) * mat::d_det(F, G)
-         + gammaDerivative(J) * mat::dd_det(F, G, H);
+  return mu() * Frobenius(sym(H), G);
 }
 
 string NeoHookeMaterial::Name() const {
   return "Neo Hooke";
 }
 
-double NeoHookeMaterial::gammaFunction(double J) const {
-  return lambda() * 0.5 * (J - 1) * (J - 1) - mu() * log(J);
-}
-
-double NeoHookeMaterial::gammaDerivative(double J) const {
-  return lambda() * (J - 1) - mu() / J;
-}
-
-double NeoHookeMaterial::gammaSecondDerivative(double J) const {
-  return lambda() + mu() / (J * J);
-}
 
 std::function<double(const Tensor &)>
 NeoHookeMaterial::ConstitutiveMatrix(const Tensor &F, const Tensor &Q, const Tensor &H) const {

cardmech/src/elasticity/materials/NeoHookeMaterial.hpp

@@ -8,13 +8,6 @@ class NeoHookeMaterial : public HyperelasticMaterial {
 
   double mu() const { return GetParameter(1); }
 
-
-  double gammaFunction(double J) const;;
-
-  double gammaDerivative(double J) const;;
-
-  double gammaSecondDerivative(double J) const;
-
   Tensor constitutiveMatrix(const Tensor &F, const Tensor &Q, const Tensor &H) const;
 
 public:

cardmech/src/elasticity/materials/VolumetricPenalty.hpp

@@ -5,7 +5,9 @@ class VolumetricPenalty {
 protected:
   double penalty = 0.0;
 public:
-  VolumetricPenalty() = default;
+  VolumetricPenalty() : penalty(20000.0) {
+    config.get("VolumetricPenalty", penalty);
+  };
 
   explicit VolumetricPenalty(double lambda) : penalty(lambda) {};
 
@@ -32,23 +34,32 @@ public:
 };
 
 class CiarletPenalty : public VolumetricPenalty {
+  double mu{1.0};
 public:
-  CiarletPenalty() : VolumetricPenalty(20000.0) {
-    config.get("VolumetricPenalty", penalty);
+  CiarletPenalty(bool mixed = false) : VolumetricPenalty() {
+    if (mixed && penalty > 0.0) {
+      config.get("PermeabilityPenalty", mu);
+      mu /= penalty;
+    }
   }
 
-  double Functional(double J) const override { return 0.5 * (J - 1.0) * (J - 1.0) - log(J); };
+  CiarletPenalty(double p) : VolumetricPenalty(p) {}
 
-  double Derivative(double J) const override { return (J - 1.0) - 1 / J; };
+  CiarletPenalty(double p, double m) : VolumetricPenalty(p), mu(p > 0.0 ? m / p : 1.0) {
+  }
 
-  double SecondDerivative(double J) const override { return 1.0 + 1 / (J * J); };
+  double Functional(double J) const override { return 0.5 * (J - 1.0) * (J - 1.0) - mu * log(J); };
+
+  double Derivative(double J) const override { return (J - 1.0) - mu / J; };
+
+  double SecondDerivative(double J) const override { return 1.0 + mu / (J * J); };
 };
 
 class LogarithmicCompressiblePenalty : public VolumetricPenalty {
 public:
-  LogarithmicCompressiblePenalty() : VolumetricPenalty(20000.0) {
-    config.get("VolumetricPenalty", penalty);
-  }
+  LogarithmicCompressiblePenalty() : VolumetricPenalty() {}
+
+  LogarithmicCompressiblePenalty(double p) : VolumetricPenalty(p) {}
 
   double Functional(double J) const override { return 0.5 * (J - 1.0) * log(J); };
 
@@ -58,18 +69,27 @@ public:
 };
 
 class MixedCompressiblePenalty : public VolumetricPenalty {
+  double mu{0.0};
 public:
-  MixedCompressiblePenalty() : VolumetricPenalty(20000.0) {
-    config.get("VolumetricPenalty", penalty);
+  MixedCompressiblePenalty() : VolumetricPenalty() {
+    if (penalty > 0.0) {
+      config.get("PermeabilityPenalty", mu);
+      mu /= penalty;
+    }
+
+  }
+
+  MixedCompressiblePenalty(double l, double m) : VolumetricPenalty(l), mu(l > 0.0 ? m / l : 0.0) {
   }
 
   double Functional(double J) const override {
-    return 0.25 * ((J - 1.0) * (J - 1.0) + log(J) * log(J));
+    return 0.25 * ((J - 1.0) * (J - 1.0) + mu * log(J) * log(J));
   };
 
-  double Derivative(double J) const override { return 0.5 * ((J - 1.0) + log(J) / J); };
+  double Derivative(double J) const override { return 0.5 * ((J - 1.0) + mu * log(J) / J); };
+
+  double SecondDerivative(double J) const override { return 0.5 * (1.0 - mu * log(J)) / J / J; };
 
-  double SecondDerivative(double J) const override { return 0.5 * (1.0 - log(J)) / J / J; };
 };
 
 #endif //VOLUMETRICPENALTY_H

cardmech/test/elasticity/materials/TestNeoHookeMaterial.hpp

-#ifndef CARDMECH_TESTNEOHOOKEMATERIAL_HPP
-#define CARDMECH_TESTNEOHOOKEMATERIAL_HPP
+#ifndef TESTNEOHOOKEMATERIAL_HPP
+#define TESTNEOHOOKEMATERIAL_HPP
 
 #include "materials/NeoHookeMaterial.hpp"
 #include "TestDefaultMaterial.hpp"
 
 INSTANTIATE_TYPED_TEST_SUITE_P(NeoHooke, MaterialTest, NeoHookeMaterial);
+
+
+std::vector<MaterialTestParameters> neohookeTestParameters{
+    {
+        std::vector{1.0, 1.0},
+        Tensor{3, 0, 0, 0, 1, 0, 0, 0, 1},
+        exp(1.0) - 1.0,
+        exp(1) * 0.75,
+        5.2666710 //exp(1) * (0.25 + 81.0 / 16.0)
+    }
+};
+
+class NeoHookeMaterialTest : public TestWithParam<MaterialTestParameters> {
+protected:
+  Material mat;
+public:
+  NeoHookeMaterialTest() : mat("NeoHooke", GetParam().matParams) {
+
+  }
+};
+
+TEST_P(NeoHookeMaterialTest, PenaltyInclusion) {
+  double lambda = GetParam().matParams[0];
+  double mu = GetParam().matParams[1];
+
+  Tensor F = 2 * One;
+
+  EXPECT_DOUBLE_EQ(lambda * mat.VolEnergy(det(F)), lambda * 0.5 * 49.0 - mu * log(8.0));
+}
+
+INSTANTIATE_TEST_SUITE_P(GuccioneMaterial, NeoHookeMaterialTest, ValuesIn(neohookeTestParameters)
+);
+
 #endif //CARDMECH_TESTNEOHOOKEMATERIAL_HPP