PARMG cleanup, write out deformationfunction to a file, Continue the trustregionsolver with a smaller trustregion radius in case of an exception

• Clean up PARMG compiler-flags
• Write out the deformation to a file; this can then be read by dune-gfe
• Continue the trustregionsolver with a smaller trustregion radius in case of an exception @oliver.sander_at_tu-dresden.de: Can you have a look?

dune/elasticity/materials/mooneyrivlinenergy.hh

@@ -31,24 +31,28 @@ public:
    */
   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");
-    
+    // Type of mooneyrivlin energy definition
     mooneyrivlin_energy = parameters.get<std::string>("mooneyrivlin_energy");
 
+    // Mooneyrivlin constants
+    if (mooneyrivlin_energy == "ciarlet") {
+      mooneyrivlin_a = parameters.get<double>("mooneyrivlin_a");
+      mooneyrivlin_b = parameters.get<double>("mooneyrivlin_b");
+      mooneyrivlin_c = parameters.get<double>("mooneyrivlin_c");
+    } else if (mooneyrivlin_energy == "log" or mooneyrivlin_energy == "square") {
+      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");
+    } else {
+      DUNE_THROW(Exception, "Error: Selected mooneyrivlin implementation (" << mooneyrivlin_energy << ") not available!");
+    }
   }
 
   /** \brief Assemble the energy for a single element */
@@ -84,10 +88,7 @@ energy(const Entity& element,
   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());
@@ -162,9 +163,10 @@ energy(const Entity& element,
     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 + 
+    if (mooneyrivlin_energy == "ciarlet") {
+      strainEnergy += weight * (mooneyrivlin_a*normFSquared + mooneyrivlin_b*normFinvSquared*detF + mooneyrivlin_c*detF*detF - ((dim-1)*mooneyrivlin_a + mooneyrivlin_b + 2*mooneyrivlin_c)*std::log(detF));
+    } else { 
+      strainEnergy += weight * (mooneyrivlin_10 * trCTildeMinus3 + 
                   mooneyrivlin_01 * c2TildeMinus3 + 
                   mooneyrivlin_20 * trCTildeMinus3 * trCTildeMinus3 + 
                   mooneyrivlin_02 * c2TildeMinus3 * c2TildeMinus3 + 
@@ -172,21 +174,17 @@ energy(const Entity& element,
                   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;
+                  mooneyrivlin_03 * c2TildeMinus3 * c2TildeMinus3 * c2TildeMinus3);
+      if (mooneyrivlin_energy == "log") {
+        field_type logDetF = std::log(detF);
+        strainEnergy += weight * ( 0.5 * mooneyrivlin_k* logDetF * logDetF );
+      } else if (mooneyrivlin_energy == "square") {
+        strainEnergy += weight * (mooneyrivlin_k* detFMinus1 * detFMinus1);
+      }
+    }
   }
-  std::cout << std::fixed;
 
-  return strainEnergyCiarlet;
+  return strainEnergy;
 }
 
 }  // namespace Dune