modify material_neukamm to 3Phases

dune/microstructure/prestrain_material_geometry.hh

@@ -53,7 +53,8 @@ public:
 
     }
     else if (imp == "material_neukamm"){    
-      std::array<double,2> mu = parameters.get<std::array<double,2>>("mu", {1.0,3.0});
+//       std::array<double,2> mu = parameters.get<std::array<double,2>>("mu", {1.0,3.0});
+      std::array<double,3> mu = parameters.get<std::array<double,3>>("mu", {1.0,3.0,2.0});
 
       //             Python::Module module = Python::import(parameters.get<std::string>("material_neukamm"));
       Python::Module module = Python::import("material_neukamm");
@@ -64,8 +65,10 @@ public:
 		      //              std::cout << "Test:" << indicatorFunction(z) << std::endl;
 		      if (indicatorFunction(z) == 1) 
 			return mu[0];
-		      else 
+		      else if (indicatorFunction(z) == 2) 
 			return mu[1];
+              else 
+            return mu[2];
 		    };
       return muTerm;
     }
@@ -491,7 +494,8 @@ public:
       return lambdaTerm;
     }
     if (imp == "material_neukamm"){    
-      std::array<double,2> lambda = parameters.get<std::array<double,2>>("lambda", {1.0,3.0});
+//       std::array<double,2> lambda = parameters.get<std::array<double,2>>("lambda", {1.0,3.0});
+      std::array<double,3> lambda = parameters.get<std::array<double,3>>("lambda", {1.0,3.0,2.0});
 
       Python::Module module = Python::import("material_neukamm");
       auto indicatorFunction = Python::make_function<double>(module.get("f"));
@@ -501,8 +505,10 @@ public:
 			  //              std::cout << "Test:" << indicatorFunction(z) << std::endl;
 			  if (indicatorFunction(z) == 1) 
 			    return lambda[0];
-			  else 
+			  else if (indicatorFunction(z) == 2) 
 			    return lambda[1];
+              else
+                return lambda[2];
 			};
       return lambdaTerm;
     }
@@ -995,22 +1001,33 @@ public:
       return B;
     }
     else if (imp == "material_neukamm"){    
-      std::array<double,2> rho = parameters.get<std::array<double,2>>("rho", {1.0,3.0});
+//       std::array<double,2> rho = parameters.get<std::array<double,2>>("rho", {1.0,3.0});
 
       //             Python::Module module = Python::import(parameters.get<std::string>("two_phase_material_1"));
       Python::Module module = Python::import("material_neukamm");
       auto indicatorFunction = Python::make_function<double>(module.get("f"));
-
-    
-      Func2Tensor B = [rho,indicatorFunction] (const Domain& x) 
+      
+      auto B1 = Python::make_function<MatrixRT>(module.get("b1"));
+      auto B2 = Python::make_function<MatrixRT>(module.get("b2"));
+      auto B3 = Python::make_function<MatrixRT>(module.get("b3"));
+      
+//       Func2Tensor B = [rho,indicatorFunction] (const Domain& x) 
+      Func2Tensor B = [indicatorFunction,B1,B2,B3] (const Domain& x) 
 		      {
 			//              std::cout << "Test:" << indicatorFunction(z) << std::endl;
 			if (indicatorFunction(x) == 1) 
-			  return MatrixRT{{rho[0], 0.0 , 0.0}, {0.0, rho[0], 0.0}, {0.0, 0.0, rho[0]}};
-			else 
-			  return MatrixRT{{rho[1], 0.0 , 0.0}, {0.0, rho[1], 0.0}, {0.0, 0.0, rho[1]}};
+            {
+// 			  return MatrixRT{{rho[0], 0.0 , 0.0}, {0.0, rho[0], 0.0}, {0.0, 0.0, rho[0]}};
+//               printmatrix(std::cout, B1(x), "Matrix B1(x)", "--");
+              return B1(x);
+            }
+			else if (indicatorFunction(x) == 2) 
+// 			  return MatrixRT{{rho[1], 0.0 , 0.0}, {0.0, rho[1], 0.0}, {0.0, 0.0, rho[1]}};
+              return B2(x);
+            else 
+              return B3(x);
 		      };
-      std::cout <<" Prestrain Type: two_phase_material_1" << std::endl;
+      std::cout <<" Prestrain Type: material_neukamm" << std::endl;
       return B;
     }
     else if (imp == "two_phase_material_2"){    

geometries/material_neukamm.py

@@ -12,6 +12,16 @@ def f(x):
     if ((abs(x[0]) < theta/2) and x[2]<-1/2+theta):
         return 1    #Phase1
     elif ((abs(x[1]) < factor*theta/2) and x[2]>1/2-theta):
-        return 1    #Phase1
+        return 2    #Phase2
     else :
-        return 0    #Phase2
+        return 0    #Phase3
+
+
+def b1(x):
+    return [[1, 4, 5], [-5, 8, 9], [-5, 8, 0]]
+
+def b2(x):
+    return [[9, 1, 9], [1, 9, 1], [1, 1, 9]]
+
+def b3(x):
+    return [[0,0,0], [0, 0, 0], [0, 8, 0]]

inputs/cellsolver.parset

@@ -10,7 +10,7 @@
 #############################################
 ### Remove/Comment this when running via Python-Script:
 outputPath=/home/stefan/DUNE/dune-microstructure/outputs
-
+#outputPath=/home/klaus/Desktop/DUNE/dune-microstructure/outputs
 
 
 #############################################
@@ -64,8 +64,8 @@ lambda1=1.0
 
 
 # better: pass material parameters as a vector
-mu=80 60
-lambda=80 25
+mu=80 60 20
+lambda=80 25 10
 rho=1.0 0
 
 
@@ -80,6 +80,7 @@ theta=0.25
 
 #--- choose composite-Implementation:
 geometryFunctionPath = /home/stefan/DUNE/dune-microstructure/geometries
+#geometryFunctionPath = /home/klaus/Desktop/DUNE/dune-microstructure/geometries
 material_prestrain_imp= "material_neukamm"   #(Python-indicator-function with same name determines material phases)
 #material_prestrain_imp= "two_phase_material_2"  #(Python-indicator-function with same name determines material phases)
 #material_prestrain_imp= "two_phase_material_3"  #(Python-indicator-function with same name determines material phases)