Simplify the calculation in mooneyrivlindensity.hh: Do not calulate the...

Simplify the calculation in mooneyrivlindensity.hh: Do not calulate the eigenvalues of F^TF and replace it by a function of trace(F^TF)

dune/elasticity/materials/mooneyrivlindensity.hh

@@ -58,53 +58,62 @@ public:
         for (int k=0; k<dim; k++)
           C[i][j] += gradient[k][i] * gradient[k][j];
 
-    //////////////////////////////////////////////////////////
-    //  Eigenvalues of C = F^TF
-    //////////////////////////////////////////////////////////
-
-
-    // eigenvalues of C, i.e., squared singular values \sigma_i of F
-    Dune::FieldVector<field_type, dim> sigmaSquared;
-    FMatrixHelp::eigenValues(C, sigmaSquared);
+    /* 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 + mooneyrivlin_c*(detF)^2 -
+                       ((dim-1)*mooneyrivlin_a + mooneyrivlin_b + 2*mooneyrivlin_c)*ln(detF)
+       log:     W(F) = \sum_{i,j=0}^{i+j<=n} mooneyrivlin_ij * (I1 - 3)^i * (I2 - 3)^j + mooneyrivlin_k * ln(det(F))
+       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)
+       
+       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)
+          = (det(F)^(-2/dim)) * trace(C)
+          = (det(F)^(-2/dim)) * trace(F^TF)
+          = (det(F)^(-2/dim)) * (frobenius_norm(F))^2
+       I2 = (det(F)^(-4/dim)) * [ second invariant of C ]
+          = (det(F)^(-4/dim)) * 0.5 * [ trace(C)^2 - tr(C^2) ]
+          = (det(F)^(-4/dim)) * [C_11C_22 + C_11C_33 + C_22C_33 - C_12C_21 - C_13C_31 - C_23C_32]
+          = (det(F)^(-4/dim)) * [C_11C_22 + C_11C_33 + C_22C_33 - C_12^2 - C_13^2 - C_23^2]
+    */
 
-    field_type normFSquared = gradient.frobenius_norm2();
+    field_type frobeniusNormFsquared = gradient.frobenius_norm2();
     field_type detF = gradient.determinant();
 
     if (detF < 0)
       DUNE_THROW( Dune::Exception, "det(F) < 0, so it is not possible to calculate the MooneyRivlinEnergy.");
 
-    field_type normFinvSquared = 0;
-
-    field_type c2Tilde = 0;
-    for (int i = 0; i < dim; i++) {
-      normFinvSquared += 1/sigmaSquared[i];
-      // compute D, which is the sum of the squared eigenvalues
-      for (int j = i+1; j < dim; j++)
-        c2Tilde += sigmaSquared[j]*sigmaSquared[i];
-      }
-
-    field_type trCTildeMinus3 = normFSquared/pow(detF, 2.0/dim) - 3;
-    // \tilde{D} = \frac{1}{\det{F}^{4/3}}D -> divide by det(F)^(4/3)= det(C)^(2/3)
-    c2Tilde /= pow(detF, 4.0/dim);
-    field_type c2TildeMinus3 = c2Tilde - 3;
-
     // Add the local energy density
     field_type strainEnergy = 0;
 
     if (mooneyrivlin_energy == "ciarlet")
     {
+      //To calculate the squared frobeniusnorm of F^(-1), we actually invert F^(-1)
+      auto gradientInverse = gradient;
+      gradientInverse.invert();
+      field_type frobeinusNormFInverseSquared = gradientInverse.frobenius_norm2();
       using std::log;
-      return mooneyrivlin_a*normFSquared + mooneyrivlin_b*normFinvSquared*detF + mooneyrivlin_c*detF*detF - ((dim-1)*mooneyrivlin_a + mooneyrivlin_b + 2*mooneyrivlin_c)*log(detF);
-    } else {
-      strainEnergy = mooneyrivlin_10 * trCTildeMinus3 +
-                        mooneyrivlin_01 * c2TildeMinus3  +
-                        mooneyrivlin_20 * trCTildeMinus3 * trCTildeMinus3 +
-                        mooneyrivlin_02 * c2TildeMinus3  * c2TildeMinus3  +
-                        mooneyrivlin_11 * trCTildeMinus3 * c2TildeMinus3  +
-                        mooneyrivlin_30 * trCTildeMinus3 * trCTildeMinus3 * trCTildeMinus3 +
-                        mooneyrivlin_21 * trCTildeMinus3 * trCTildeMinus3 * c2TildeMinus3  +
-                        mooneyrivlin_12 * trCTildeMinus3 * c2TildeMinus3  * c2TildeMinus3  +
-                        mooneyrivlin_03 * c2TildeMinus3  * c2TildeMinus3  * c2TildeMinus3;
+      return mooneyrivlin_a*frobeniusNormFsquared + mooneyrivlin_b*frobeinusNormFInverseSquared*detF + mooneyrivlin_c*detF*detF - ((dim-1)*mooneyrivlin_a + mooneyrivlin_b + 2*mooneyrivlin_c)*log(detF);
+    } else { // mooneyrivlin_energy is "log" or "square"
+      field_type a = pow(detF, 2.0/dim);
+      field_type invariant1Minus3 = frobeniusNormFsquared/a - 3;
+      field_type secondInvariantOfC = 0;
+      for (int i=0; i<dim-1; i++)
+        for (int j=i+1; j<dim; j++)
+          secondInvariantOfC += C[i][i]*C[j][j] - C[i][j]*C[i][j];
+      field_type invariant2Minus3 = secondInvariantOfC/(a*a) - 3;
+      strainEnergy = mooneyrivlin_10 * invariant1Minus3 +
+                        mooneyrivlin_01 * invariant2Minus3  +
+                        mooneyrivlin_20 * invariant1Minus3 * invariant1Minus3 +
+                        mooneyrivlin_02 * invariant2Minus3  * invariant2Minus3  +
+                        mooneyrivlin_11 * invariant1Minus3 * invariant2Minus3  +
+                        mooneyrivlin_30 * invariant1Minus3 * invariant1Minus3 * invariant1Minus3 +
+                        mooneyrivlin_21 * invariant1Minus3 * invariant1Minus3 * invariant2Minus3  +
+                        mooneyrivlin_12 * invariant1Minus3 * invariant2Minus3  * invariant2Minus3  +
+                        mooneyrivlin_03 * invariant2Minus3  * invariant2Minus3  * invariant2Minus3;
       if (mooneyrivlin_energy == "log") {
         using std::log;
         field_type logDetF = log(detF);

test/CMakeLists.txt

@@ -8,5 +8,8 @@ dune_add_test(SOURCES materialtest
 add_dune_ug_flags(materialtest)
 add_dune_adolc_flags(materialtest)
 
+dune_add_test(SOURCES mooneyrivlintest
+              CMAKE_GUARD ADOLC_FOUND)
+
 dune_add_test(SOURCES localhyperdualstiffnesstest
               CMAKE_GUARD ADOLC_FOUND)

test/mooneyrivlintest.cc

+#include <config.h>
+
+// Includes for the ADOL-C automatic differentiation library
+// Need to come before (almost) all others.
+#include <adolc/adouble.h>
+#include <dune/fufem/utilities/adolcnamespaceinjections.hh>
+
+#include <dune/common/parametertree.hh>
+#include <dune/common/parallel/mpihelper.hh>
+
+#include <dune/elasticity/assemblers/localadolcstiffness.hh>
+#include <dune/elasticity/assemblers/feassembler.hh>
+#include <dune/elasticity/materials/localintegralenergy.hh>
+#include <dune/elasticity/materials/mooneyrivlindensity.hh>
+
+#include <dune/functions/functionspacebases/interpolate.hh>
+#include <dune/functions/functionspacebases/lagrangebasis.hh>
+#include <dune/functions/functionspacebases/powerbasis.hh>
+
+#include <dune/grid/utility/structuredgridfactory.hh>
+#include <dune/grid/uggrid.hh>
+
+// grid dimension
+const int dim = 3;
+
+//order of integration
+const int order = 1;
+
+using namespace Dune;
+
+//differentiation method: ADOL-C
+typedef adouble ValueType;
+
+using namespace Dune;
+
+typedef BlockVector<FieldVector<double,dim> > VectorType;
+typedef BCRSMatrix<FieldMatrix<double, dim, dim> > MatrixType;
+
+template <class Basis>
+int assembleAndCompare (const Basis basis, const ParameterTree parameters, const VectorType x,
+                        double expectedEnergy, double expectedGradientTwoNorm, double expectedGradientInfinityNorm, double expectedMatrixFrobeniusNorm) {
+  using LocalView = typename Basis::LocalView;
+
+  std::string mooneyrivlin_energy = parameters.get<std::string>("mooneyrivlin_energy");
+
+  auto elasticDensity = std::make_shared<Elasticity::MooneyRivlinDensity<dim,ValueType>>(parameters);
+  auto elasticEnergy = std::make_shared<Elasticity::LocalIntegralEnergy<LocalView, adouble>>(elasticDensity);
+  auto localADOLCStiffness = std::make_shared<Elasticity::LocalADOLCStiffness<LocalView>>(elasticEnergy);
+  Elasticity::FEAssembler<Basis,VectorType> assembler(basis, localADOLCStiffness);
+
+  //////////////////////////////////////////////////////////////
+  //  Compute the energy, assemble and compare!
+  //////////////////////////////////////////////////////////////
+  
+  VectorType gradient;
+  MatrixType hessianMatrix;
+  double energy = assembler.computeEnergy(x);
+
+  if ( std::abs(energy - expectedEnergy)/expectedEnergy > 1e-3)
+  {
+      std::cerr << std::setprecision(9);
+      std::cerr << "In " << mooneyrivlin_energy << ": Calculated energy is " << energy << " but '" << expectedEnergy << "' was expected!" << std::endl;
+      return 1;
+  }
+
+  assembler.assembleGradientAndHessian(x, gradient, hessianMatrix, true);
+
+  double gradientTwoNorm = gradient.two_norm();
+  double gradientInfinityNorm = gradient.infinity_norm();
+  double matrixFrobeniusNorm = hessianMatrix.frobenius_norm();
+  if ( std::abs(gradientTwoNorm - expectedGradientTwoNorm)/expectedGradientTwoNorm > 1e-4 || std::abs(gradientInfinityNorm - expectedGradientInfinityNorm)/expectedGradientInfinityNorm > 1e-8)
+  {
+      std::cerr << std::setprecision(9);
+      std::cerr << "In " << mooneyrivlin_energy << ": Calculated gradient max norm is " << gradientInfinityNorm << " but '" << expectedGradientInfinityNorm << "' was expected!" << std::endl;
+      std::cerr << "In " << mooneyrivlin_energy << ": Calculated gradient norm is " << gradientTwoNorm << " but '" << expectedGradientTwoNorm << "' was expected!" << std::endl;
+      return 1;
+  }
+
+  if ( std::abs(matrixFrobeniusNorm - expectedMatrixFrobeniusNorm)/expectedMatrixFrobeniusNorm > 1e-4)
+  {
+      std::cerr << std::setprecision(9);
+      std::cerr << "In " << mooneyrivlin_energy << ": Calculated matrix norm is " << matrixFrobeniusNorm << " but '" << expectedMatrixFrobeniusNorm << "' was expected!" << std::endl;
+      return 1;
+  }
+
+  return 0;
+}
+
+int main (int argc, char *argv[])
+{
+
+  Dune::MPIHelper& mpiHelper = MPIHelper::instance(argc, argv);
+
+  // ///////////////////////////////////////
+  //    Create the grid
+  // ///////////////////////////////////////
+  using GridType = UGGrid<dim>;
+  auto grid = StructuredGridFactory<GridType>::createCubeGrid({0,0,0}, {1,1,1}, {3,3,3});
+  int refine = 3;
+  while(refine > 0) {
+    for (auto&& e : elements(grid->leafGridView())){
+      bool refine = false;
+      for (int i = 0; i < e.geometry().corners(); i++) {
+          refine = refine || (e.geometry().corner(i)[0] > 0.99);
+      }
+      grid->mark(refine ? 1 : 0, e);
+    }
+    grid->adapt();
+    refine--;
+  }
+  grid->loadBalance();
+
+  using GridView = GridType::LeafGridView;
+  GridView gridView = grid->leafGridView();
+
+  /////////////////////////////////////////////////////////////////////////
+  //  Create a power basis
+  /////////////////////////////////////////////////////////////////////////
+
+  using namespace Functions::BasisFactory;
+  auto basis = makeBasis(
+    gridView,
+    power<dim>(
+      lagrange<order>(),
+      blockedInterleaved()
+  ));
+
+  using PowerBasis = decltype(basis);
+
+  //////////////////////////////////////////////////////////////
+  //  Prepare the iterate x where we want to assemble
+  //////////////////////////////////////////////////////////////
+  VectorType x(basis.size());
+  Functions::interpolate(basis, x, [](FieldVector<double,dim> x){
+    FieldVector<double,dim> y = {0.5*x[0]*x[0], 0.5*x[1] + 0.2*x[0], 4*std::abs(std::sin(x[2]))};
+    return y;
+  });
+
+  ParameterTree parameters;
+  parameters["mooneyrivlin_10"] = "1.67e+6";
+  parameters["mooneyrivlin_01"] = "1.94e+6";
+  parameters["mooneyrivlin_20"] = "2.42e+6";
+  parameters["mooneyrivlin_02"] = "6.52e+6";
+  parameters["mooneyrivlin_11"] = "-7.34e+6";
+  parameters["mooneyrivlin_30"] = "0";
+  parameters["mooneyrivlin_03"] = "0";
+  parameters["mooneyrivlin_21"] = "0";
+  parameters["mooneyrivlin_12"] = "0";
+  parameters["mooneyrivlin_k"] = "75e+6";
+
+  parameters["mooneyrivlin_energy"] = "square";
+  double expectedEnergy = 177758377;
+  double expectedGradientTwoNorm = 2.13791422e+09;
+  double expectedGradientInfinityNorm = 599188506;
+  double expectedMatrixFrobeniusNorm = 1.62529884e+11;
+  int testSquare = assembleAndCompare(basis, parameters, x, expectedEnergy, expectedGradientTwoNorm, expectedGradientInfinityNorm, expectedMatrixFrobeniusNorm);
+
+  parameters["mooneyrivlin_energy"] = "log";
+  expectedEnergy = 181180273;
+  expectedGradientTwoNorm = 2.29399362e+09;
+  expectedGradientInfinityNorm = 648551554;
+  expectedMatrixFrobeniusNorm = 1.67152642e+11;
+  int testLog = assembleAndCompare(basis, parameters, x, expectedEnergy, expectedGradientTwoNorm, expectedGradientInfinityNorm, expectedMatrixFrobeniusNorm);
+
+  parameters["mooneyrivlin_energy"] = "ciarlet";
+  parameters["mooneyrivlin_a"] = "2.42e+6";
+  parameters["mooneyrivlin_b"] = "6.52e+6";
+  parameters["mooneyrivlin_c"] = "-7.34e+6";
+  expectedEnergy = 83069427.2;
+  expectedGradientTwoNorm = 163210957;
+  expectedGradientInfinityNorm = 49284369.2;
+  expectedMatrixFrobeniusNorm = 5.92126562e+09;
+  int testCiarlet = assembleAndCompare(basis, parameters, x, expectedEnergy, expectedGradientTwoNorm, expectedGradientInfinityNorm, expectedMatrixFrobeniusNorm);
+
+  return testCiarlet + testLog + testSquare;
+}