Clean up mooneyrivlinenergy: Now the program only takes the values it actually...

Clean up mooneyrivlinenergy: Now the program only takes the values it actually needs from the parameter-file.

dune/elasticity/materials/mooneyrivlinenergy.hh

@@ -31,11 +31,15 @@ public:
    */
   MooneyRivlinEnergy(const Dune::ParameterTree& parameters)
   {
+    // 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");
@@ -46,9 +50,9 @@ public:
       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");
-
+    } 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;
+                  mooneyrivlin_03 * c2TildeMinus3 * c2TildeMinus3 * c2TildeMinus3);
       if (mooneyrivlin_energy == "log") {
-    return strainEnergyWithLog;
+        field_type logDetF = std::log(detF);
+        strainEnergy += weight * ( 0.5 * mooneyrivlin_k* logDetF * logDetF );
       } else if (mooneyrivlin_energy == "square") {
-    return strainEnergyWithSquare;
+        strainEnergy += weight * (mooneyrivlin_k* detFMinus1 * detFMinus1);
+      }
+    }
   }
-  std::cout << std::fixed;
 
-  return strainEnergyCiarlet;
+  return strainEnergy;
 }
 
 }  // namespace Dune