GetMaxParameter for DG

cardmech/conf/elasticity_benchmarks/beam.conf

@@ -57,7 +57,7 @@ StretchFactorNormal = 4.0
 ######################
 # Material Parameters#
 ######################
-#ActiveMaterial = Linear
+ActiveMaterial = Linear
 #ActiveMaterial = NeoHooke
 #ActiveMaterial = Bonet
 ActiveMaterial = Holzapfel
@@ -71,7 +71,7 @@ Density = 0.00
 # Volumetric Penalty
 #------------------------
 Incompressible = false
-QuasiCompressiblePenalty = Ciarlet
+QuasiCompressiblePenalty = None
 VolumetricPenalty = 2
 PermeabilityPenalty = 2
 
@@ -148,7 +148,7 @@ presmoothing = 3;
 postsmoothing = 3;
 
 DGSign=-1.0
-DGPenalty=100.0
+DGPenalty= 20.0
 Overlap = dG1
 #Overlap = NoOverlap
 Overlap_Distribution = 1

cardmech/conf/m++.conf

@@ -16,7 +16,7 @@ GeoPath = ../geo/
 # === Cardiac Mechanics ================
 #loadconf = coupled.conf
 #loadconf = klotz.conf
-#loadconf = elasticity_benchmarks/beam.conf
+loadconf = elasticity_benchmarks/beam.conf
 loadconf = elasticity_benchmarks/ellipsoid.conf
 loadconf = elasticity.conf
 # ======================================

cardmech/src/elasticity/assemble/DGElasticity.cpp

@@ -157,7 +157,8 @@ void DGElasticity::Residual(const cell &c, const Vector &U, Vector &r) const {
       }
   }
   for (int f = 0; f < c.Faces(); ++f) {
-    Scalar s = cellMat.Isochoric().GetParameter(0) *( pow(degree,2) * penalty) / U.GetMesh().MaxMeshWidth();
+    Scalar s = cellMat.Isochoric().GetMaxParameter() *( pow(degree,2) * penalty) / U.GetMesh().MaxMeshWidth();
+
     DGVectorFieldFaceElement FE(U, c, f);
     if (E.Bnd(f) >= 230 && E.Bnd(f) <= 233) { // Neumann-Randbedingung
       for (int q = 0; q < FE.nQ(); ++q) {
@@ -234,8 +235,8 @@ void DGElasticity::Residual(const cell &c, const Vector &U, Vector &r) const {
       DGVectorFieldElement E_nghbr(U, cf);
       DGVectorFieldFaceElement FE_nghbr(r, cf, f1);
 
-      s = cellMat.Isochoric().GetParameter(0) * ( pow(degree,2) * penalty) / U.GetMesh().MaxMeshWidth();
-      for (int q = 0; q < FE.nQ(); ++q) {
+      s = cellMat.Isochoric().GetMaxParameter() * ( pow(degree,2) * penalty) / U.GetMesh().MaxMeshWidth();
+        for (int q = 0; q < FE.nQ(); ++q) {
         double w = FE.QWeight(q);
         Point z = FE.QPoint(q);
         VectorField N = FE.QNormal(q);
@@ -313,7 +314,7 @@ void DGElasticity::DirichletBoundary(const cell &c, int face, int bc, const Vect
   DGVectorFieldElement E(U, c);
   DGVectorFieldFaceElement FE(U, c, face);
   const auto &cellMat = eProblem.GetMaterial(c);
-  Scalar s = cellMat.Isochoric().GetParameter(0) * ( pow(degree,2) * penalty) / U.GetMesh().MaxMeshWidth();
+  Scalar s = cellMat.Isochoric().GetMaxParameter() * ( pow(degree,2) * penalty) / U.GetMesh().MaxMeshWidth();
 
   for (int q = 0; q < FE.nQ(); ++q) {
     double w = FE.QWeight(q);
@@ -339,7 +340,7 @@ void DGElasticity::FixationBoundary(const cell &c, int face, int bc, const Vecto
   DGVectorFieldElement E(U, c);
   DGVectorFieldFaceElement FE(U, c, face);
   const auto &cellMat = eProblem.GetMaterial(c);
-  Scalar s = cellMat.Isochoric().GetParameter(0) * ( pow(degree,2) * penalty) / U.GetMesh().MaxMeshWidth();
+  Scalar s = cellMat.Isochoric().GetMaxParameter() * ( pow(degree,2) * penalty) / U.GetMesh().MaxMeshWidth();
 
   for (int q = 0; q < FE.nQ(); ++q) {
     double w = FE.QWeight(q);
@@ -478,7 +479,7 @@ void DGElasticity::Jacobi(const cell &c, const Vector &U, Matrix &A) const {
   }
   //BFParts bnd(U.GetMesh(), c);
   for (int face = 0; face < c.Faces(); ++face) {
-    Scalar s = cellMat.Isochoric().GetParameter(0) * ( pow(degree,2) * penalty) / U.GetMesh().MaxMeshWidth();
+    Scalar s = cellMat.Isochoric().GetMaxParameter() * ( pow(degree,2) * penalty) / U.GetMesh().MaxMeshWidth();
     DGVectorFieldFaceElement FE(U, c, face);
     if (E.Bnd(face) == DIRICHLET_BC || E.Bnd(face) == DIRICHLET_FIX_BC) {
       for (int q = 0; q < FE.nQ(); ++q) {

cardmech/src/elasticity/materials/GuccioneMaterial.hpp

@@ -21,6 +21,11 @@ class GuccioneMaterial : public HyperelasticMaterial {
 
   double bfs() const { return GetParameter(3); }
 
+  std::pair<size_t, size_t> GetMaxParameterRange() const {
+      // start, length
+      return {0, 4};
+  }
+
   Tensor fiberWeights(const Tensor &E) const {
     return {bf() * E[0][0], bfs() * E[0][1], bfs() * E[0][2],
             bfs() * E[1][0], bs() * E[1][1], bs() * E[1][2],

cardmech/src/elasticity/materials/HolzapfelMaterial.hpp

@@ -26,6 +26,11 @@ class HolzapfelMaterial : public HyperelasticMaterial {
 
   double bfs() const { return GetParameter(7); }
 
+  std::pair<size_t, size_t> GetMaxParameterRange() const {
+    // start, length
+    return {0, 4};
+  }
+
   /// Returns d if d is positive and 0 instead.
   double macaulay(double d) const { return (d > 0 ? d : 0.0); };
 

cardmech/src/elasticity/materials/LinearMaterial.hpp

@@ -9,6 +9,11 @@ class LinearMaterial : public HyperelasticMaterial {
 
   double mu() const { return GetParameter(1); }
 
+  std::pair<size_t, size_t> GetMaxParameterRange() const {
+      // start, length
+      return {0, 2};
+  }
+
 
   double energy(const Tensor &e_u) const;
 

cardmech/src/elasticity/materials/Material.hpp

@@ -84,6 +84,23 @@ public:
 
   double GetParameter(int i) const { return parameters[i]; }
 
+  virtual std::pair<size_t, size_t> GetMaxParameterRange() const {
+      // start, length
+      return {0, 1};
+  }
+
+  double GetMaxParameter() const{
+      auto range = GetMaxParameterRange();
+      double maxParameter = parameters[range.first];
+      for (size_t i = range.first; i <= range.second; ++i){
+          if(parameters[i] > maxParameter){
+              maxParameter = parameters[i];
+          }
+      }
+      mout << "maxParameter " << maxParameter << endl;
+      return maxParameter;
+  }
+
   void UpdateParameters(const std::vector<double> &params) {
     std::copy(params.begin(), params.end(), parameters.begin());
   }