Merge branch 'fix/mooneyrivlin' into 'master'

Fix in Mooney-Rivlin density: Use proper 4th term for Ciarlet

See merge request !69

dune/elasticity/materials/mooneyrivlindensity.hh

@@ -64,13 +64,14 @@ public:
     /* The Mooney-Rivlin-Density is given as a function of the eigenvalues of the right Cauchy-Green-Deformation tensor C = F^TF
        or the right Cauchy-Green-Deformation tensor B = FF^T.
        C = F^TF and B = FF^T have the same eigenvalues - we can use either one of them.
-       
+
        There are three Mooney-Rivlin-Variants:
        ciarlet: W(F) = mooneyrivlin_a*(normF)^2 + mooneyrivlin_b*(normFinv)^2*detF^2 + mooneyrivlin_c*(detF)^2 -
-                       ((dim-1)*mooneyrivlin_a + mooneyrivlin_b + 2*mooneyrivlin_c)*ln(detF)
+                       (2*mooneyrivlin_a + 4*mooneyrivlin_b + 2*mooneyrivlin_c)*ln(detF),
+                       where the last term is chosen s.t. W( t*I ) is minimal for t=1
        log:     W(F) = \sum_{i,j=0}^{i+j<=n} mooneyrivlin_ij * (I1 - 3)^i * (I2 - 3)^j + mooneyrivlin_k * ln(det(F))^2
        square:  W(F) = \sum_{i,j=0}^{i+j<=n} mooneyrivlin_ij * (I1 - 3)^i * (I2 - 3)^j + mooneyrivlin_k * 0.5 * (det(F) - 1)^2
-       
+
        For log and square: I1 and I2 are the first two invariants of C (or B), multiplied with a factor depending on det(F):
        I1 = (det(F)^(-2/dim)) * [ first invariant of C ]
           = (det(F)^(-2/dim)) * (sum of all eigenvalues of C)
@@ -99,7 +100,8 @@ public:
       gradientInverse.invert();
       field_type frobeinusNormFInverseSquared = gradientInverse.frobenius_norm2();
       using std::log;
-      return mooneyrivlin_a*frobeniusNormFsquared + mooneyrivlin_b*frobeinusNormFInverseSquared*detF*detF + mooneyrivlin_c*detF*detF - ((dim-1)*mooneyrivlin_a + mooneyrivlin_b + 2*mooneyrivlin_c)*log(detF);
+      return mooneyrivlin_a*frobeniusNormFsquared + mooneyrivlin_b*frobeinusNormFInverseSquared*detF*detF + mooneyrivlin_c*detF*detF
+             - (2.0*mooneyrivlin_a + 4.0*mooneyrivlin_b + 2.0*mooneyrivlin_c)*log(detF);
     }
     else
     {

test/mooneyrivlintest.cc

@@ -52,7 +52,7 @@ int assembleAndCompare (const Basis basis, const ParameterTree parameters, const
   //////////////////////////////////////////////////////////////
   //  Compute the energy, assemble and compare!
   //////////////////////////////////////////////////////////////
-  
+
   VectorType gradient;
   MatrixType hessianMatrix;
   double energy = assembler.computeEnergy(x);
@@ -90,7 +90,7 @@ int assembleAndCompare (const Basis basis, const ParameterTree parameters, const
 int main (int argc, char *argv[])
 {
 
-  Dune::MPIHelper& mpiHelper = MPIHelper::instance(argc, argv);
+  [[maybe_unused]] Dune::MPIHelper& mpiHelper = MPIHelper::instance(argc, argv);
 
   // ///////////////////////////////////////
   //    Create the grid
@@ -126,8 +126,6 @@ int main (int argc, char *argv[])
       blockedInterleaved()
   ));
 
-  using PowerBasis = decltype(basis);
-
   //////////////////////////////////////////////////////////////
   //  Prepare the iterate x where we want to assemble
   //////////////////////////////////////////////////////////////
@@ -167,10 +165,10 @@ int main (int argc, char *argv[])
   parameters["mooneyrivlin_a"] = "2.42e+6";
   parameters["mooneyrivlin_b"] = "6.52e+6";
   parameters["mooneyrivlin_c"] = "-7.34e+6";
-  expectedEnergy = 68753239.6;
-  expectedGradientTwoNorm = 31244643.1;
-  expectedGradientInfinityNorm = 9673572.39;
-  expectedMatrixFrobeniusNorm = 1.67660965e+09;
+  expectedEnergy = 76302830.4;
+  expectedGradientTwoNorm = 40670527.3;
+  expectedGradientInfinityNorm = 11116511.8;
+  expectedMatrixFrobeniusNorm = 2.1978108e+09;
   int testCiarlet = assembleAndCompare(basis, parameters, x, expectedEnergy, expectedGradientTwoNorm, expectedGradientInfinityNorm, expectedMatrixFrobeniusNorm);
 
   return testCiarlet + testLog + testSquare;