Add Mooney-Rivlin-Energy to finite-strain-elasticity

Add option to do calculations using a Mooney-Rivlin Material.
There are mainly three different formulas found in the literature, all
are implemented here and can be switched using the mooneyrivlin_energy
flag in the finite-strain-elasticity.parset.

dune/elasticity/materials/mooneyrivlinenergy.hh

+#ifndef DUNE_ELASTICITY_MATERIALS_MOONEYRIVLINENERGY_HH
+#define DUNE_ELASTICITY_MATERIALS_MOONEYRIVLINENERGY_HH
+
+#include <dune/common/fmatrix.hh>
+#include <dune/common/fmatrixev.hh>
+
+#include <dune/geometry/quadraturerules.hh>
+
+#include <dune/elasticity/assemblers/localfestiffness.hh>
+
+namespace Dune {
+
+template<class GridView, class LocalFiniteElement, class field_type=double>
+class MooneyRivlinEnergy
+  : public LocalFEStiffness<GridView,
+                            LocalFiniteElement,
+                            std::vector<Dune::FieldVector<field_type,GridView::dimension> > >
+{
+  // grid types
+  typedef typename GridView::Grid::ctype DT;
+  typedef typename GridView::template Codim<0>::Entity Entity;
+
+  // some other sizes
+  enum {gridDim=GridView::dimension};
+  enum {dim=GridView::dimension};
+
+public:
+
+  /** \brief Constructor with a set of material parameters
+   * \param parameters The material parameters
+   */
+  MooneyRivlinEnergy(const Dune::ParameterTree& parameters)
+  {
+    // Mooneyrivlin constants
+    mooneyrivlin_a = parameters.get<double>("mooneyrivlin_a");
+    mooneyrivlin_b = parameters.get<double>("mooneyrivlin_b");
+    mooneyrivlin_c = parameters.get<double>("mooneyrivlin_c");
+    
+    mooneyrivlin_10 = parameters.get<double>("mooneyrivlin_10");
+    mooneyrivlin_01 = parameters.get<double>("mooneyrivlin_01");
+    mooneyrivlin_20 = parameters.get<double>("mooneyrivlin_20");
+    mooneyrivlin_02 = parameters.get<double>("mooneyrivlin_02");
+    mooneyrivlin_11 = parameters.get<double>("mooneyrivlin_11");
+    mooneyrivlin_30 = parameters.get<double>("mooneyrivlin_30");
+    mooneyrivlin_03 = parameters.get<double>("mooneyrivlin_03");
+    mooneyrivlin_21 = parameters.get<double>("mooneyrivlin_21");
+    mooneyrivlin_12 = parameters.get<double>("mooneyrivlin_12");
+    mooneyrivlin_k = parameters.get<double>("mooneyrivlin_k");
+    
+    mooneyrivlin_energy = parameters.get<std::string>("mooneyrivlin_energy");
+
+  }
+
+  /** \brief Assemble the energy for a single element */
+  field_type energy (const Entity& e,
+                     const LocalFiniteElement& localFiniteElement,
+                     const std::vector<Dune::FieldVector<field_type, gridDim> >& localConfiguration) const;
+
+  /** \brief Lame constants */
+  field_type mooneyrivlin_a,
+              mooneyrivlin_b,
+              mooneyrivlin_c,
+              mooneyrivlin_10,
+              mooneyrivlin_01,
+              mooneyrivlin_20,
+              mooneyrivlin_02,
+              mooneyrivlin_11, 
+              mooneyrivlin_30, 
+              mooneyrivlin_21, 
+              mooneyrivlin_12, 
+              mooneyrivlin_03,
+              mooneyrivlin_k;
+  std::string mooneyrivlin_energy;
+};
+
+template <class GridView, class LocalFiniteElement, class field_type>
+field_type
+MooneyRivlinEnergy<GridView, LocalFiniteElement, field_type>::
+energy(const Entity& element,
+       const LocalFiniteElement& localFiniteElement,
+       const std::vector<Dune::FieldVector<field_type, gridDim> >& localConfiguration) const
+{
+
+  assert(element.type() == localFiniteElement.type());
+  typedef typename GridView::template Codim<0>::Entity::Geometry Geometry;
+
+  field_type strainEnergyCiarlet = 0;
+  field_type strainEnergy = 0;
+  field_type strainEnergyWithLog = 0;
+  field_type strainEnergyWithSquare = 0;
+
+  // store gradients of shape functions and base functions
+  std::vector<Dune::FieldMatrix<DT,1,gridDim> > referenceGradients(localFiniteElement.size());
+  std::vector<Dune::FieldVector<DT,gridDim> > gradients(localFiniteElement.size());
+
+  int quadOrder = (element.type().isSimplex()) ? localFiniteElement.localBasis().order()
+                                               : localFiniteElement.localBasis().order() * gridDim;
+
+  const Dune::QuadratureRule<DT, gridDim>& quad
+      = Dune::QuadratureRules<DT, gridDim>::rule(element.type(), quadOrder);
+
+  for (size_t pt=0; pt<quad.size(); pt++)
+  {
+    // Local position of the quadrature point
+    const Dune::FieldVector<DT,gridDim>& quadPos = quad[pt].position();
+
+    const DT integrationElement = element.geometry().integrationElement(quadPos);
+
+    const typename Geometry::JacobianInverseTransposed& jacobianInverseTransposed = element.geometry().jacobianInverseTransposed(quadPos);
+
+    DT weight = quad[pt].weight() * integrationElement;
+
+    // get gradients of shape functions
+    localFiniteElement.localBasis().evaluateJacobian(quadPos, referenceGradients);
+
+    // compute gradients of base functions
+    for (size_t i=0; i<gradients.size(); ++i)
+      jacobianInverseTransposed.mv(referenceGradients[i][0], gradients[i]);
+
+    Dune::FieldMatrix<field_type,gridDim,gridDim> derivative(0);
+    for (size_t i=0; i<gradients.size(); i++)
+      for (int j=0; j<gridDim; j++)
+        derivative[j].axpy(localConfiguration[i][j], gradients[i]);
+
+
+    Dune::FieldMatrix<field_type,gridDim,gridDim> C(0);
+
+    for (int i=0; i<gridDim; i++) {
+      for (int j=0; j<gridDim; j++) {
+        for (int k=0; k<gridDim; k++)
+          C[i][j] += derivative[k][i] * derivative[k][j];
+      }
+    }
+
+    //////////////////////////////////////////////////////////
+    //  Eigenvalues of FTF
+    //////////////////////////////////////////////////////////
+
+    // eigenvalues of C, i.e., squared singular values \sigma_i of F
+    Dune::FieldVector<field_type, dim> sigmaSquared;
+    FMatrixHelp::eigenValues(C, sigmaSquared);
+
+    field_type normFSquared = derivative.frobenius_norm2();
+    field_type detF = derivative.determinant();
+
+    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;
+    field_type detFMinus1 = detF - 1; 
+
+    // Add the local energy density
+    strainEnergyCiarlet += weight * (mooneyrivlin_a*normFSquared + mooneyrivlin_b*normFinvSquared*detF + mooneyrivlin_c*detF*detF - ((dim-1)*mooneyrivlin_a + mooneyrivlin_b + 2*mooneyrivlin_c)*std::log(detF));
+    strainEnergy = 0;
+    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;
+    field_type logDetF = std::log(detF);
+    strainEnergyWithLog += weight * ( strainEnergy + 0.5 * mooneyrivlin_k* logDetF * logDetF );
+    strainEnergyWithSquare += weight * ( strainEnergy + mooneyrivlin_k* detFMinus1 * detFMinus1);
+  }
+
+  std::cout << std::scientific;
+  if (mooneyrivlin_energy == "log") {
+    return strainEnergyWithLog;
+  } else if (mooneyrivlin_energy == "square") {
+    return strainEnergyWithSquare;
+  } else {
+    return strainEnergyCiarlet;
+  }
+  std::cout << std::fixed;
+}
+
+}  // namespace Dune
+
+#endif   //#ifndef DUNE_ELASTICITY_MATERIALS_MOONEYRIVLINENERGY_HH

src/finite-strain-elasticity-bending.parset

@@ -62,6 +62,22 @@ energy = stvenantkirchhoff
 mu = 2.7191e+6
 lambda = 4.4364e+6
 
+mooneyrivlin_a = 2.7191e+6
+mooneyrivlin_b = 2.7191e+6
+mooneyrivlin_c = 2.7191e+6 
+mooneyrivlin_10 = -7.28e+5
+mooneyrivlin_01 = 9.17e+5
+mooneyrivlin_20 = 1.23e+5
+mooneyrivlin_02 = 8.15e+5
+mooneyrivlin_11 = -5.14e+5
+mooneyrivlin_30 = 0
+mooneyrivlin_21 = 0
+mooneyrivlin_12 = 0
+mooneyrivlin_03 = 0
+
+mooneyrivlin_k = 1e+6	# volume-preserving parameter
+mooneyrivlin_energy = log # log, square or ciarlet; different ways to compute the Mooney-Rivlin-Energy
+
 []
 
 #############################################

src/finite-strain-elasticity.cc

@@ -38,6 +38,7 @@
 #include <dune/elasticity/materials/stvenantkirchhoffenergy.hh>
 #include <dune/elasticity/materials/henckyenergy.hh>
 #include <dune/elasticity/materials/exphenckyenergy.hh>
+#include <dune/elasticity/materials/mooneyrivlinenergy.hh>
 #include <dune/elasticity/materials/neohookeenergy.hh>
 #include <dune/elasticity/materials/neumannenergy.hh>
 #include <dune/elasticity/materials/sumenergy.hh>
@@ -285,11 +286,14 @@ int main (int argc, char *argv[]) try
       elasticEnergy = std::make_shared<HenckyEnergy<GridView,
                                                     FEBasis::LocalView::Tree::FiniteElement,
                                                     adouble> >(materialParameters);
-
     if (parameterSet.get<std::string>("energy") == "exphencky")
       elasticEnergy = std::make_shared<ExpHenckyEnergy<GridView,
                                                        FEBasis::LocalView::Tree::FiniteElement,
                                                        adouble> >(materialParameters);
+    if (parameterSet.get<std::string>("energy") == "mooneyrivlin")
+      elasticEnergy = std::make_shared<MooneyRivlinEnergy<GridView,
+                                                               FEBasis::LocalView::Tree::FiniteElement,
+                                                               adouble> >(materialParameters);
 
     if(!elasticEnergy)
       DUNE_THROW(Exception, "Error: Selected energy not available!");