diff --git a/src/dune-microstructure.cc b/src/dune-microstructure.cc
index 1a670643ad589c8834288ff070c2892f4aa4a909..639d22e39a5e0dbb0d2a9a50fbcd452039c25c45 100644
--- a/src/dune-microstructure.cc
+++ b/src/dune-microstructure.cc
@@ -44,7 +44,6 @@
 #include <dune/common/fvector.hh>
 #include <dune/common/fmatrix.hh>
 
-
 #include <dune/microstructure/prestrain_material_geometry.hh>
 #include <dune/microstructure/matrix_operations.hh>
 #include <dune/microstructure/vtk_filler.hh>    //TEST
@@ -211,6 +210,7 @@ void computeElementStiffnessMatrix(const LocalView& localView,
             defGradientV[l][2] = (1.0/gamma)*gradients[j][2];           //X3
               
             double energyDensity = linearizedStVenantKirchhoffDensity(mu(quadPos), lambda(quadPos), defGradientU, defGradientV);
+//             double energyDensity = generalizedDensity(mu(quadPos), lambda(quadPos), defGradientU, defGradientV);  // also works..
 
             size_t col = localView.tree().child(k).localIndex(i);                       
             size_t row = localView.tree().child(l).localIndex(j);
@@ -510,8 +510,8 @@ auto energy(const Basis& basis,
 //   using float_50 = boost::multiprecision::cpp_dec_float_50;
 //   float_50 higherPrecEnergy = 0.0;
   double energy = 0.0;
-  constexpr int dim = 3;
-  constexpr int dimWorld = 3;
+  constexpr int dim = Basis::LocalView::Element::dimension;
+  constexpr int dimWorld = dim;
 
   auto localView = basis.localView();
 
@@ -578,8 +578,9 @@ void setOneBaseFunctionToZero(const Basis& basis,
                               )
 {
   std::cout << "Setting one Basis-function to zero" << std::endl;
-  constexpr int dim = 3;
 
+  constexpr int dim = Basis::LocalView::Element::dimension;
+  
   unsigned int arbitraryLeafIndex =  parameterSet.template get<unsigned int>("arbitraryLeafIndex", 0);
   unsigned int arbitraryElementNumber =  parameterSet.template get<unsigned int>("arbitraryElementNumber", 0);
   //Determine 3 global indices (one for each component of an arbitrary local FE-function)
@@ -645,7 +646,8 @@ auto integralMean(const Basis& basis,
                   )
 {
   // computes integral mean of given LocalFunction
-  constexpr int dim = 3;
+
+  constexpr int dim = Basis::LocalView::Element::dimension;
 
   auto GVFunction = Functions::makeDiscreteGlobalBasisFunction<FieldVector<double,dim>>(basis,coeffVector);
   auto localfun = localFunction(GVFunction);
@@ -686,7 +688,8 @@ auto subtractIntegralMean(const Basis& basis,
 {
   // Substract correct Integral mean from each associated component function
   auto IM = integralMean(basis, coeffVector);
-  constexpr int dim = 3;
+
+  constexpr int dim = Basis::LocalView::Element::dimension;
 
   for(size_t k=0; k<dim; k++)
   {
@@ -727,7 +730,7 @@ double computeL2SymError(const Basis& basis,
 
   auto localView = basis.localView();
   
-  //   constexpr int dim = 3;
+
   constexpr int dim = Basis::LocalView::Element::dimension;             //TODO TEST 
   constexpr int dimWorld = 3;                                           // Hier auch möglich? 
   
@@ -802,8 +805,9 @@ double computeL2Norm(const Basis& basis,
 {
   // IMPLEMENTATION with makeDiscreteGlobalBasisFunction
   double error = 0.0;
-  constexpr int dim = 3;
-  constexpr int dimWorld = 3;
+
+  constexpr int dim = Basis::LocalView::Element::dimension;
+  constexpr int dimWorld = dim;
   auto localView = basis.localView();
   auto GVFunc = Functions::makeDiscreteGlobalBasisFunction<FieldVector<double,dim>>(basis,coeffVector);
   auto localfun = localFunction(GVFunc);
@@ -853,6 +857,7 @@ int main(int argc, char *argv[])
   std::fstream log;
   log.open(outputPath ,std::ios::out);
 
+  
   constexpr int dim = 3;
   constexpr int dimWorld = 3;
 
@@ -1065,6 +1070,13 @@ int main(int argc, char *argv[])
     Func2Tensor x3G_1neg = [x3G_1] (const Domain& x) {return -1.0*x3G_1(x);};
     Func2Tensor x3G_2neg = [x3G_2] (const Domain& x) {return -1.0*x3G_2(x);};
     Func2Tensor x3G_3neg = [x3G_3] (const Domain& x) {return -1.0*x3G_3(x);};
+    
+    
+    //TEST 
+//     auto QuadraticForm = [] (const double mu, const double lambda, const MatrixRT& M) {
+// 
+//                                 return lambda*std::pow(trace(M),2) + 2*mu*pow(norm(sym(M)),2);
+//                         };
 
 
     // TEST - energy method ///
@@ -1416,12 +1428,14 @@ int main(int argc, char *argv[])
     log << "computed q1: " << q1_c << std::endl;
     log << "computed q2: " << q2_c << std::endl;
     log << "computed q3: " << q3_c << std::endl;
+    log << "mu_gamma=" << q3_c << std::endl;           // added for Python-Script
     log << "computed b1: " << Beff[0] << std::endl;
     log << "computed b2: " << Beff[1] << std::endl;
     log << "computed b3: " << Beff[2] << std::endl;
     log << "computed b1_hat: " << B_hat[0] << std::endl;
     log << "computed b2_hat: " << B_hat[1] << std::endl;
     log << "computed b3_hat: " << B_hat[2] << std::endl;
+    
 
     //////////////////////////////////////////////////////////////
     // Define Analytic Solutions
@@ -1650,7 +1664,7 @@ int main(int argc, char *argv[])
               << center("L2SymError",tableWidth)     << " | "
               << center("Order",tableWidth)     << " | "
               << center("L2SymNorm",tableWidth)     << " | "
-              << center("L2SNorm_ana",tableWidth)     << " | "
+              << center("L2SymNorm_ana",tableWidth)     << " | "
               << center("L2Norm",tableWidth)  << " | " << "\n";
     std::cout << std::string(tableWidth*6 + 3*6, '-') << "\n";
     log       << center("Levels",tableWidth)       << " | "