diff --git a/dune/microstructure/prestrain_material_geometry.hh b/dune/microstructure/prestrain_material_geometry.hh
index 8e429bf0683bdfa3e3db803d3388329f78d2cd39..f86071c976f51522f7160f07076899939a3d6d9c 100644
--- a/dune/microstructure/prestrain_material_geometry.hh
+++ b/dune/microstructure/prestrain_material_geometry.hh
@@ -46,12 +46,27 @@ public:
 
 
     	if (imp == "homogeneous"){    
-		    double mu     = parameters.get<double>("mu", 10);
+		    double mu1     = parameters.get<double>("mu1", 10);
 
-		    auto muTerm = [mu] (const Domain& z) {return mu;};
+		    auto muTerm = [mu1] (const Domain& z) {return mu1;};
 		    
 			return muTerm;
 		}
+		else if (imp == "isotropic_bilayer")
+    	{
+            double mu1     = parameters.get<double>("mu1",10.0);
+            double beta = parameters.get<double>("beta",2.0); 
+		    double mu2 = beta*mu1;
+    
+ 		   auto muTerm = [mu1, mu2, theta] (const Domain& z) {
+                    if (z[2]>0)                                                   
+                        return mu1;
+                    else
+                        return mu2;
+                    };
+		    
+           return muTerm;
+		}
         else if (imp == "analytical_Example"){	
 		    double mu1     = parameters.get<double>("mu1",10.0);
             double beta = parameters.get<double>("beta",2.0); 
@@ -391,8 +406,21 @@ public:
 		    
 			return lambdaTerm;
 		}
-		
-		
+        else if (imp == "isotropic_bilayer")
+    	{
+            double lambda1 = parameters.get<double>("lambda1",0.0);
+            double beta = parameters.get<double>("beta",2.0); 
+            double lambda2 = beta*lambda1;
+    
+ 		   auto lambdaTerm = [lambda1,lambda2, theta] (const Domain& z) {
+                    if (z[2]>0)                                                   
+                        return lambda1;
+                    else
+                        return lambda2;
+                    };
+           return lambdaTerm;
+		}
+
         else if (imp == "analytical_Example"){
             double lambda1 = parameters.get<double>("lambda1",0.0);
             double beta = parameters.get<double>("beta",2.0); 
@@ -405,8 +433,7 @@ public:
                         return lambda1;
                     else
                         return lambda2;
-                    };
-		    
+                    }; 
            return lambdaTerm;
 		}
 		else if (imp == "circle_fiber"){
@@ -771,13 +798,28 @@ public:
         double alpha = parameters.get<double>("alpha", 2.0);
         double p2 = alpha*p1;
         
-    	if (imp == "isotropic_bilayer")
+//     	if (imp == "isotropic_bilayer")
+//     	{
+//             Func2Tensor B = [p1,p2,theta] (const Domain& x)                 //  ISOTROPIC PRESSURE
+//             {             
+//                 if (abs(x[0]) > (theta/2)  && x[2] > 0)
+//                     return MatrixRT{{p1, 0.0 , 0.0}, {0.0, p1, 0.0}, {0.0, 0.0, p1}};
+//                 else if (abs(x[0]) < (theta/2)  && x[2] < 0)
+//                     return MatrixRT{{p2, 0.0 , 0.0}, {0.0, p2, 0.0}, {0.0, 0.0, p2}};
+//                 else
+//                     return MatrixRT{{0.0, 0.0 , 0.0}, {0.0, 0.0, 0.0}, {0.0, 0.0, 0.0}};
+// 
+//             };
+//             std::cout <<" Prestrain Type: isotropic_bilayer " << std::endl;
+//             return B;
+//         }
+        if (imp == "isotropic_bilayer")
     	{
             Func2Tensor B = [p1,p2,theta] (const Domain& x)                 //  ISOTROPIC PRESSURE
             {             
-                if (abs(x[0]) > (theta/2)  && x[2] > 0)
+                if (x[2] > 0)
                     return MatrixRT{{p1, 0.0 , 0.0}, {0.0, p1, 0.0}, {0.0, 0.0, p1}};
-                else if (abs(x[0]) < (theta/2)  && x[2] < 0)
+                else if (x[2] < 0)
                     return MatrixRT{{p2, 0.0 , 0.0}, {0.0, p2, 0.0}, {0.0, 0.0, p2}};
                 else
                     return MatrixRT{{0.0, 0.0 , 0.0}, {0.0, 0.0, 0.0}, {0.0, 0.0, 0.0}};
@@ -786,6 +828,7 @@ public:
             std::cout <<" Prestrain Type: isotropic_bilayer " << std::endl;
             return B;
         }
+        
         else if (imp == "analytical_Example")
         {
             Func2Tensor B = [p1,theta] (const Domain& x)                // Bilayer with one rectangular Fiber & ISOTROPIC PRESSURE   
diff --git a/src/Cell-Problem.cc b/src/Cell-Problem.cc
index 21f0e885b0856c53f662fcaedd496fee31f12e62..dbc2cfd525c614b7137e331113e986d0742a8796 100644
--- a/src/Cell-Problem.cc
+++ b/src/Cell-Problem.cc
@@ -876,6 +876,7 @@ int main(int argc, char *argv[])
   ///////////////////////////////////
   // Get Material Parameters
   ///////////////////////////////////
+  std::string imp = parameterSet.get<std::string>("material_implementation", "analytical_Example");
   double beta = parameterSet.get<double>("beta",2.0); 
   double mu1 = parameterSet.get<double>("mu1",10.0);;
   double mu2 = beta*mu1;
@@ -1430,12 +1431,6 @@ int main(int argc, char *argv[])
     std::cout<< "q1 : " << q1 << std::endl;
     std::cout<< "q2 : " << q2 << std::endl;
     std::cout<< "q3 : " << q3 << std::endl;
-//     std::cout<< "b1hat_c: " << B_hat[0] << std::endl;
-//     std::cout<< "b2hat_c: " << B_hat[1] << std::endl;
-//     std::cout<< "b3hat_c: " << B_hat[2] << std::endl;
-//     std::cout<< "b1eff_c: " << Beff[0] << std::endl;
-//     std::cout<< "b2eff_c: " << Beff[1] << std::endl;
-//     std::cout<< "b3eff_c: " << Beff[2] << std::endl;
     printvector(std::cout, B_hat, "B_hat", "--");
     printvector(std::cout, Beff, "Beff", "--");
     
@@ -1504,12 +1499,25 @@ int main(int argc, char *argv[])
     log << "q2_ana : "     << q2_ana << std::endl;
     log << "q3 should be between q1 and q2"  << std::endl;
     
-    Storage_Quantities.push_back(std::abs(q1_ana - q1));
-    Storage_Quantities.push_back(std::abs(q2_ana - q2));
-    Storage_Quantities.push_back( q3 );
-    Storage_Quantities.push_back(std::abs(b1_eff_ana - Beff[0]));
-    Storage_Quantities.push_back(std::abs(b2_eff_ana - Beff[1]));
-    Storage_Quantities.push_back(std::abs(b3_eff_ana - Beff[2]));
+    if (imp == "analytical_Example")   // print Errors only for analytical_Example
+    {
+        Storage_Quantities.push_back(std::abs(q1_ana - q1));
+        Storage_Quantities.push_back(std::abs(q2_ana - q2));
+        Storage_Quantities.push_back(q3);
+        Storage_Quantities.push_back(std::abs(b1_eff_ana - Beff[0]));
+        Storage_Quantities.push_back(std::abs(b2_eff_ana - Beff[1]));
+        Storage_Quantities.push_back(std::abs(b3_eff_ana - Beff[2]));
+    }
+    else
+    {
+        Storage_Quantities.push_back(q1);
+        Storage_Quantities.push_back(q2);
+        Storage_Quantities.push_back(q3);
+        Storage_Quantities.push_back(Beff[0]);
+        Storage_Quantities.push_back(Beff[1]);
+        Storage_Quantities.push_back(Beff[2]);
+    }
+
 
     auto symPhi_1_analytic = [mu1, mu2, theta, muTerm] (const Domain& x) {
                         return MatrixRT{{  (((mu1*mu2)/((theta*mu1 +(1-theta)*mu2)*muTerm(x))) - 1)*x[2] , 0.0, 0.0}, {0.0, 0.0, 0.0}, {0.0, 0.0, 0.0}};
@@ -1672,60 +1680,88 @@ int main(int argc, char *argv[])
   
     
   
-  
 
     /////////////////////////////////////////
     // Print Storage
     /////////////////////////////////////////
     int tableWidth = 12;
-    std::cout << center("Levels",tableWidth)       << " | "
-              << center("L2SymError",tableWidth)     << " | "
-              << center("Order",tableWidth)     << " | "
-              << center("L2SymNorm",tableWidth)     << " | "
-              << center("L2SymNorm_ana",tableWidth)     << " | "
-              << center("L2Norm",tableWidth)  << " | " << "\n";
-    std::cout << std::string(tableWidth*6 + 3*6, '-') << "\n";
-    log       << center("Levels",tableWidth)       << " | "
-              << center("L2SymError",tableWidth)     << " | "
-              << center("Order",tableWidth)     << " | "
-              << center("L2SymNorm",tableWidth)     << " | "
-              << center("L2SNorm_ana",tableWidth)     << " | "
-              << center("L2Norm",tableWidth)  << " | " << "\n";
-    log       << std::string(tableWidth*6 + 3*6, '-') << "\n";
-    
-    
-    int StorageCount = 0;
-
-    for(auto& v: Storage_Error) 
+    if (imp == "analytical_Example")   // print Errors only for analytical_Example
     {
-        std::visit([tableWidth](auto&& arg){std::cout << center(prd(arg,5,1),tableWidth)      << " | ";}, v);      // Anzahl-Nachkommastellen
-        std::visit([tableWidth, &log](auto&& arg){log << center(prd(arg,5,1),tableWidth)      << " & ";}, v);
+        std::cout << std::string(tableWidth*7 + 3*7, '-') << "\n";
+        std::cout << center("Levels",tableWidth)       << " | "
+                << center("L2SymError",tableWidth)     << " | "
+                << center("Order",tableWidth)     << " | "
+                << center("L2SymNorm",tableWidth)     << " | "
+                << center("L2SymNorm_ana",tableWidth)     << " | "
+                << center("L2Norm",tableWidth)  << " | " << "\n";
+        std::cout << std::string(tableWidth*6 + 3*6, '-') << "\n";
+        log     << std::string(tableWidth*6 + 3*6, '-') << "\n";
+        log     << center("Levels",tableWidth)       << " | "
+                << center("L2SymError",tableWidth)     << " | "
+                << center("Order",tableWidth)     << " | "
+                << center("L2SymNorm",tableWidth)     << " | "
+                << center("L2SNorm_ana",tableWidth)     << " | "
+                << center("L2Norm",tableWidth)  << " | " << "\n";
+        log     << std::string(tableWidth*6 + 3*6, '-') << "\n";
+        
+        int StorageCount = 0;
 
-        StorageCount++;
-        if(StorageCount % 6 == 0 )
+        for(auto& v: Storage_Error) 
         {
-            std::cout << std::endl;
-            log << std::endl;
+            std::visit([tableWidth](auto&& arg){std::cout << center(prd(arg,5,1),tableWidth)      << " | ";}, v);      // Anzahl-Nachkommastellen
+            std::visit([tableWidth, &log](auto&& arg){log << center(prd(arg,5,1),tableWidth)      << " & ";}, v);
+            StorageCount++;
+            if(StorageCount % 6 == 0 )
+            {
+                std::cout << std::endl;
+                log << std::endl;
+            }
         }
     }
-
+    
     //////////////// OUTPUT QUANTITIES TABLE //////////////
+    if (imp == "analytical_Example")   // print Errors only for analytical_Example
+    {
+    std::cout << std::string(tableWidth*7 + 3*7, '-') << "\n";
     std::cout << center("Levels ",tableWidth)       << " | "
               << center("|q1_ana-q1|",tableWidth)       << " | "
               << center("|q2_ana-q2|",tableWidth)       << " | "
-              << center(" q3_c",tableWidth)           << " | "
+              << center("q3",tableWidth)           << " | "
               << center("|b1_ana-b1|",tableWidth)       << " | "
               << center("|b2_ana-b2|",tableWidth)       << " | "
               << center("|b3_ana-b3|",tableWidth)       << " | " << "\n";
     std::cout << std::string(tableWidth*7 + 3*7, '-') << "\n";
+    log       << std::string(tableWidth*7 + 3*7, '-') << "\n";
     log       << center("Levels ",tableWidth)       << " | "
               << center("|q1_ana-q1|",tableWidth)       << " | "
               << center("|q2_ana-q2|",tableWidth)       << " | "
-              << center(" q3",tableWidth)           << " | "
+              << center("q3",tableWidth)           << " | "
               << center("|b1_ana-b1|",tableWidth)       << " | "
               << center("|b2_ana-b2|",tableWidth)       << " | "
               << center("|b3_ana-b3|",tableWidth)       << " | " << "\n";
     log       << std::string(tableWidth*7 + 3*7, '-') << "\n";
+    }
+    else
+    {
+    std::cout << center("Levels ",tableWidth)       << " | "
+              << center("q1",tableWidth)       << " | "
+              << center("q2",tableWidth)       << " | "
+              << center("q3",tableWidth)           << " | "
+              << center("b1",tableWidth)       << " | "
+              << center("b2",tableWidth)       << " | "
+              << center("b3",tableWidth)       << " | " << "\n";
+    std::cout << std::string(tableWidth*7 + 3*7, '-') << "\n";
+    log       << std::string(tableWidth*7 + 3*7, '-') << "\n";   
+    log       << center("Levels ",tableWidth)       << " | "
+              << center("q1",tableWidth)       << " | "
+              << center("q2",tableWidth)       << " | "
+              << center("q3",tableWidth)           << " | "
+              << center("b1",tableWidth)       << " | "
+              << center("b2",tableWidth)       << " | "
+              << center("b3",tableWidth)       << " | " << "\n";
+    log       << std::string(tableWidth*7 + 3*7, '-') << "\n";   
+        
+    }
     
     int StorageCount2 = 0;
     for(auto& v: Storage_Quantities) 
@@ -1739,6 +1775,8 @@ int main(int argc, char *argv[])
             log << std::endl;
         }
     }
+    std::cout << std::string(tableWidth*7 + 3*7, '-') << "\n";
+    log       << std::string(tableWidth*7 + 3*7, '-') << "\n";  
 
     log.close();
 }