diff --git a/dune/elasticity/materials/mooneyrivlindensity.hh b/dune/elasticity/materials/mooneyrivlindensity.hh
index a95f758187dbd34a53b80f92b0f0f6172ea36229..e067923871b8e0be80c3ae078e55a734b2397f74 100644
--- a/dune/elasticity/materials/mooneyrivlindensity.hh
+++ b/dune/elasticity/materials/mooneyrivlindensity.hh
@@ -19,23 +19,26 @@ public:
   */
   MooneyRivlinDensity(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");
+    // 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 Evaluation with the deformation gradient
@@ -84,11 +87,14 @@ public:
     // \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
-    auto strainEnergyCiarlet = (mooneyrivlin_a*normFSquared + mooneyrivlin_b*normFinvSquared*detF + mooneyrivlin_c*detF*detF - ((dim-1)*mooneyrivlin_a + mooneyrivlin_b + 2*mooneyrivlin_c)*std::log(detF));
-    auto strainEnergy = mooneyrivlin_10 * trCTildeMinus3 +
+    field_type strainEnergy = 0;
+
+    if (mooneyrivlin_energy == "ciarlet")
+      return 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 = mooneyrivlin_10 * trCTildeMinus3 +
                         mooneyrivlin_01 * c2TildeMinus3  +
                         mooneyrivlin_20 * trCTildeMinus3 * trCTildeMinus3 +
                         mooneyrivlin_02 * c2TildeMinus3  * c2TildeMinus3  +
@@ -97,20 +103,14 @@ public:
                         mooneyrivlin_21 * trCTildeMinus3 * trCTildeMinus3 * c2TildeMinus3  +
                         mooneyrivlin_12 * trCTildeMinus3 * c2TildeMinus3  * c2TildeMinus3  +
                         mooneyrivlin_03 * c2TildeMinus3  * c2TildeMinus3  * c2TildeMinus3;
-    field_type logDetF = std::log(detF);
-    auto strainEnergyWithLog =  ( strainEnergy + 0.5 * mooneyrivlin_k* logDetF * logDetF );
-    auto strainEnergyWithSquare =  ( strainEnergy + mooneyrivlin_k* detFMinus1 * detFMinus1);
-
-
-    std::cout << std::scientific;
-    if (mooneyrivlin_energy == "log") {
-      return strainEnergyWithLog;
-    } else if (mooneyrivlin_energy == "square") {
-      return strainEnergyWithSquare;
+      if (mooneyrivlin_energy == "log") {
+        field_type logDetF = std::log(detF);
+        return strainEnergy + 0.5 * mooneyrivlin_k* logDetF * logDetF;
+      } else if (mooneyrivlin_energy == "square") {
+        field_type detFMinus1 = detF - 1;
+        return strainEnergy + mooneyrivlin_k* detFMinus1 * detFMinus1;
+      }
     }
-    std::cout << std::fixed;
-
-    return strainEnergyCiarlet;
   }
 
   /** \brief Lame constants */