Add Python-Interface from dune-fufem and geometries-folder for definition of indicator-functions

dune.module

@@ -7,6 +7,6 @@ Module: dune-microstructure
 Version: 1.0
 Maintainer: klaus.boehnlein@tu-dresden.de
 # Required build dependencies
-Depends: dune-common dune-istl dune-geometry dune-localfunctions dune-uggrid dune-grid dune-typetree dune-functions 
+Depends: dune-common dune-istl dune-geometry dune-localfunctions dune-uggrid dune-grid dune-typetree dune-functions dune-fufem
 # Optional build dependencies
 #Suggests:

dune/microstructure/prestrain_material_geometry.hh

@@ -6,6 +6,8 @@
 #include <dune/grid/yaspgrid.hh>
 #include <dune/microstructure/matrix_operations.hh>
 
+#include <dune/fufem/dunepython.hh>
+
 using namespace Dune;
 using namespace MatrixOperations;
 using std::pow;
@@ -16,8 +18,6 @@ using std::cos;
 
 
 
-
-
 template <int dim> 
 class IsotropicMaterialImp 
 {
@@ -44,11 +44,64 @@ public:
 
 
     	if (imp == "homogeneous"){    
+
 		    double mu1     = parameters.get<double>("mu1", 10);
 
 		    auto muTerm = [mu1] (const Domain& z) {return mu1;};
 		    
 			return muTerm;
+
+		}
+    	else if (imp == "two_phase_material_1"){    
+            std::array<double,2> mu = parameters.get<std::array<double,2>>("mu", {1.0,3.0});
+
+//             Python::Module module = Python::import(parameters.get<std::string>("two_phase_material_1"));
+            Python::Module module = Python::import("two_phase_material_1");
+            auto indicatorFunction = Python::make_function<double>(module.get("f"));
+
+            auto muTerm = [mu,indicatorFunction] (const Domain& z) 
+            {
+//              std::cout << "Test:" << indicatorFunction(z) << std::endl;
+                if (indicatorFunction(z) == 1) 
+                    return mu[0];
+                else 
+                    return mu[1];
+            };
+			return muTerm;
+		}
+    	else if (imp == "two_phase_material_2"){    
+            std::array<double,2> mu = parameters.get<std::array<double,2>>("mu", {1.0,3.0});
+
+//             Python::Module module = Python::import(parameters.get<std::string>("two_phase_material_1"));
+            Python::Module module = Python::import("two_phase_material_2");
+            auto indicatorFunction = Python::make_function<double>(module.get("f"));
+
+            auto muTerm = [mu,indicatorFunction] (const Domain& z) 
+            {
+//              std::cout << "Test:" << indicatorFunction(z) << std::endl;
+                if (indicatorFunction(z) == 1) 
+                    return mu[0];
+                else 
+                    return mu[1];
+            };
+			return muTerm;
+		}
+    	else if (imp == "two_phase_material_3"){    
+            std::array<double,2> mu = parameters.get<std::array<double,2>>("mu", {1.0,3.0});
+
+//             Python::Module module = Python::import(parameters.get<std::string>("two_phase_material_1"));
+            Python::Module module = Python::import("two_phase_material_3");
+            auto indicatorFunction = Python::make_function<double>(module.get("f"));
+
+            auto muTerm = [mu,indicatorFunction] (const Domain& z) 
+            {
+//              std::cout << "Test:" << indicatorFunction(z) << std::endl;
+                if (indicatorFunction(z) == 1) 
+                    return mu[0];
+                else 
+                    return mu[1];
+            };
+			return muTerm;
 		}
 		else if (imp == "isotropic_bilayer")
     	{
@@ -413,11 +466,62 @@ public:
         double width = parameters.get<double>("width", 1.0);
         double height = parameters.get<double>("height", 1.0);
 
-		if (imp == "homogenenous"){    
-		    double lambda     = parameters.get<double>("lambda",384.615);
+		if (imp == "homogeneous"){    
+		    double lambda1 = parameters.get<double>("lambda1",0.0);
 
-		    auto lambdaTerm = [lambda] (const Domain& z) {return lambda;};
+		    auto lambdaTerm = [lambda1] (const Domain& z) {return lambda1;};
 		    
+			return lambdaTerm;
+		}
+    	if (imp == "two_phase_material_1"){    
+            std::array<double,2> lambda = parameters.get<std::array<double,2>>("lambda", {1.0,3.0});
+
+//             Python::Module module = Python::import(parameters.get<std::string>("two_phase_material_1"));
+            Python::Module module = Python::import("two_phase_material_1");
+            auto indicatorFunction = Python::make_function<double>(module.get("f"));
+
+            auto lambdaTerm = [lambda,indicatorFunction] (const Domain& z) 
+            {
+//              std::cout << "Test:" << indicatorFunction(z) << std::endl;
+                if (indicatorFunction(z) == 1) 
+                    return lambda[0];
+                else 
+                    return lambda[1];
+            };
+			return lambdaTerm;
+		}
+    	if (imp == "two_phase_material_2"){    
+            std::array<double,2> lambda = parameters.get<std::array<double,2>>("lambda", {1.0,3.0});
+
+//             Python::Module module = Python::import(parameters.get<std::string>("two_phase_material_1"));
+            Python::Module module = Python::import("two_phase_material_2");
+            auto indicatorFunction = Python::make_function<double>(module.get("f"));
+
+            auto lambdaTerm = [lambda,indicatorFunction] (const Domain& z) 
+            {
+//              std::cout << "Test:" << indicatorFunction(z) << std::endl;
+                if (indicatorFunction(z) == 1) 
+                    return lambda[0];
+                else 
+                    return lambda[1];
+            };
+			return lambdaTerm;
+		}
+    	if (imp == "two_phase_material_3"){    
+            std::array<double,2> lambda = parameters.get<std::array<double,2>>("lambda", {1.0,3.0});
+
+//             Python::Module module = Python::import(parameters.get<std::string>("two_phase_material_1"));
+            Python::Module module = Python::import("two_phase_material_3");
+            auto indicatorFunction = Python::make_function<double>(module.get("f"));
+
+            auto lambdaTerm = [lambda,indicatorFunction] (const Domain& z) 
+            {
+//              std::cout << "Test:" << indicatorFunction(z) << std::endl;
+                if (indicatorFunction(z) == 1) 
+                    return lambda[0];
+                else 
+                    return lambda[1];
+            };
 			return lambdaTerm;
 		}
         else if (imp == "isotropic_bilayer")
@@ -812,10 +916,8 @@ protected:
 
 public:
 
-//     PrestrainImp(double p1, double p2, double theta, double width) : p1(p1), p2(p2), theta(theta), width(width){}
 
     
-    
 //     Func2Tensor getPrestrain(std::string imp)
     Func2Tensor getPrestrain(ParameterTree parameters)
     {
@@ -830,22 +932,74 @@ public:
         double alpha = parameters.get<double>("alpha", 2.0);
         double p2 = alpha*p1;
         
-//     	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")
+
+        
+        if (imp == "homogeneous"){    
+            Func2Tensor B = [p1] (const Domain& x)                 //  ISOTROPIC PRESSURE
+            {             
+                return MatrixRT{{p1, 0.0 , 0.0}, {0.0, p1, 0.0}, {0.0, 0.0, p1}};
+            };
+            std::cout <<" Prestrain Type: homogeneous" << std::endl;
+            return B;
+		}
+    	else if (imp == "two_phase_material_1"){    
+            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("two_phase_material_1");
+            auto indicatorFunction = Python::make_function<double>(module.get("f"));
+
+    
+            Func2Tensor B = [rho,indicatorFunction] (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]}};
+            };
+            std::cout <<" Prestrain Type: two_phase_material_1" << std::endl;
+			return B;
+		}
+    	else if (imp == "two_phase_material_2"){    
+            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("two_phase_material_2");
+            auto indicatorFunction = Python::make_function<double>(module.get("f"));
+
+    
+            Func2Tensor B = [rho,indicatorFunction] (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]}};
+            };
+            std::cout <<" Prestrain Type: two_phase_material_2" << std::endl;
+			return B;
+		}
+    	else if (imp == "two_phase_material_3"){    
+            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("two_phase_material_3");
+            auto indicatorFunction = Python::make_function<double>(module.get("f"));
+
+    
+            Func2Tensor B = [rho,indicatorFunction] (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]}};
+            };
+            std::cout <<" Prestrain Type: two_phase_material_3" << std::endl;
+			return B;
+		}
+        else if (imp == "isotropic_bilayer")
     	{
             Func2Tensor B = [p1,p2,theta] (const Domain& x)                 //  ISOTROPIC PRESSURE
             {             

geometries/two_phase_material_1.py

+import math
+
+
+#Indicator function that determines both phases
+# x[0] : x-component
+# x[1] : y-component
+# x[2] : z-component
+def f(x):
+    # --- replace with your definition of indicatorFunction:
+    if (abs(x[0]) > 0.25):
+        return 1    #Phase1
+    else :
+        return 0    #Phase2

geometries/two_phase_material_2.py

+import math
+
+
+#Indicator function that determines both phases
+# x[0] : x-component
+# x[1] : y-component
+# x[2] : z-component
+def f(x):
+    # --- replace with your definition of indicatorFunction:
+    if (abs(x[1]) > 0.25):
+        return 1    #Phase1
+    else :
+        return 0    #Phase2

geometries/two_phase_material_3.py

+import math
+
+
+#Indicator function that determines both phases
+# x[0] : x-component
+# x[1] : y-component
+# x[2] : z-component
+def f(x):
+    # --- replace with your definition of indicatorFunction:
+    if (abs(x[2]) > 0.25):
+        return 1    #Phase1
+    else :
+        return 0    #Phase2

inputs/cellsolver.parset

@@ -11,8 +11,13 @@
 #outputPath = "../../outputs/output.txt"
 #outputPath = "/home/klaus/Desktop/DUNE/dune-microstructure/outputs/output.txt"
 #outputPath = "/home/klaus/Desktop/DUNE/dune-microstructure/outputs"
+
+
 ### Remove/Comment this when running via Python-Script:
-outputPath = "../../outputs"
+#outputPath = "../../outputs"
+
+
+
 
 #############################################
 #  Cell Domain
@@ -28,50 +33,50 @@ cellDomain=1
 ## {start,finish} computes on all grid from 2^(start) to 2^finish refinement
 #----------------------------------------------------
 
+numLevels=2 2
 #numLevels =  1 1   # computes all levels from first to second entry
 #numLevels =  2 2   # computes all levels from first to second entry
 #numLevels =  1 3   # computes all levels from first to second entry
-numLevels = 3 3     # computes all levels from first to second entry
 #numLevels = 4 4    # computes all levels from first to second entry
 #numLevels = 5 5    # computes all levels from first to second entry
 #numLevels = 6 6    # computes all levels from first to second entry
 #numLevels = 1 6
 
+
+
 ########################################################################################
 
 # --- Choose scale ratio gamma:
-gamma=1.0  #(default = 1.0)
-#gamma=50.0
-#gamma=0.01
-#gamma=3.0
-#gamma=0.5
-#gamma = 0.05
+gamma=0.45
+
 
 #############################################
 #  Material / Prestrain parameters and ratios
 #############################################
 
 #-- stiffness-ratio (ratio between material parameters mu1 & mu2 .... beta = 1.0 corresponds to homogeneous case)
-beta = 2.0
-
+beta=3.0
 
 $--- strength of prestrain
-rho1 = 1.0
+rho1=1.0
 
 #--- prestrain-contrast (ratio between prestrain parameters rho1 & rho2)
-#alpha = 5.0
-alpha = 2.0
+alpha=8.0
 
 
 #--- Lame-Parameters
 mu1=1.0
+lambda1=1.0
+
 
-lambda1=0.0
-#lambda1=5.0
+# better: pass material parameters as a vector
+mu=1.0 3.0
+lambda=1.0 3.0
+rho=1.0 8.0
 
 
 # ---volume fraction  (default value = 1.0/4.0)
-theta = 0.25
+theta=0.25
 #theta = 0.3
 #theta = 0.75
 #theta = 0.125
@@ -80,7 +85,18 @@ theta = 0.25
 
 
 #--- choose composite-Implementation:
-material_prestrain_imp= "parametrized_Laminate"
+
+geometryFunctionPath = /home/klaus/Desktop/DUNE/dune-microstructure/geometries
+#geometryFunction = two_phase_material_1
+material_prestrain_imp= "two_phase_material_1"   #(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)
+
+
+
+
+#material_prestrain_imp= "parametrized_Laminate"
+#material_prestrain_imp= "homogeneous"
 #material_prestrain_imp= "analytical_Example"
 #material_prestrain_imp="isotropic_bilayer"
 #material_prestrain_imp= "circle_fiber"    #TEST
@@ -91,8 +107,8 @@ material_prestrain_imp= "parametrized_Laminate"
 
 
 # --- (Optional output) write Material / prestrain / Corrector functions to .vtk-Files (default=false):
-#write_materialFunctions = true
-#write_prestrainFunctions = true  # VTK norm of B ,
+write_materialFunctions = true
+write_prestrainFunctions = true  # VTK norm of B ,
 #write_VTK = true
 
 

microstructure_testsuite/helper_functions.py

+import subprocess
+import re
+import os
+import numpy as np
+import matplotlib.pyplot as plt
+import math
+import fileinput
+import time
+import matplotlib.ticker as tickers
+import matplotlib as mpl
+from matplotlib.ticker import MultipleLocator,FormatStrFormatter,MaxNLocator
+import codecs
+import sys
+
+
+#-------------------------------------------------------------------------------------------------------
+def run_CellProblem(executable, parset, gridLevel, gamma, mu1,lambda1, rho1, alpha, beta, theta, material_prestrain_imp, outputPath, write_materialFunctions, write_prestrainFunctions, write_VTK,write_L2Error ,write_IntegralMean, write_LOG ):
+    print('----- RUN Cell-Problem ----')
+    processList = []
+    LOGFILE = "Cell-Problem_output.log"
+    print('LOGFILE:',LOGFILE)
+    print('executable:',executable)
+    print('parset:',parset)
+
+    # test = {5,5}
+    # test = '5 5'
+    # levels = [gridLevel, gridLevel]
+    # print(listToString(gridLevel))
+        # start_time = time.time()
+    if write_LOG:
+        p = subprocess.Popen(executable + parset
+                                        # + " -numLevels " + str(str(gridLevel) + "  " + str(gridLevel))
+                                        # + " -numLevels " + str(4) + " " + str(4)
+                                        # + " -numLevels " + {listToString(gridLevel)}
+                                        # + " -numLevels " + " " + test
+                                        # + " -gamma " + str(gamma)
+                                        # + " -mu1 " + str(mu1)
+                                        # + " -lambda1 " + str(lambda1)
+                                        # + " -alpha " + str(alpha)
+                                        # + " -beta " + str(beta)
+                                        # + " -theta " + str(theta)
+                                        + " -material_prestrain_imp " + str(material_prestrain_imp)
+                                        + " -write_materialFunctions " + str(write_materialFunctions)
+                                        + " -write_prestrainFunctions " + str(write_prestrainFunctions)
+                                        + " -write_VTK " + str(write_VTK)
+                                        + " -write_L2Error " + str(write_L2Error)
+                                        + " -write_IntegralMean" + str(write_IntegralMean)
+                                        + " -outputPath " + str(outputPath)
+                                        + " | tee " + LOGFILE, shell=True)
+
+    else:
+        p = subprocess.Popen(executable + parset
+                                        # + " -numLevels " + str(gridLevel) + " " + str(gridLevel)
+                                        # + " -gamma " + str(gamma)
+                                        # + " -mu1 " + str(mu1)
+                                        # + " -lambda1 " + str(lambda1)
+                                        # + " -alpha " + str(alpha)
+                                        # + " -beta " + str(beta)
+                                        # + " -theta " + str(theta)
+                                        + " -material_prestrain_imp " + str(material_prestrain_imp)
+                                        + " -write_materialFunctions " + str(write_materialFunctions)
+                                        + " -write_prestrainFunctions " + str(write_prestrainFunctions)
+                                        + " -write_VTK " + str(write_VTK)
+                                        + " -write_L2Error " + str(write_L2Error)
+                                        + " -write_IntegralMean" + str(write_IntegralMean)
+                                        + " -outputPath " + str(outputPath), stdout=subprocess.DEVNULL, stderr=subprocess.DEVNULL, shell=True )  # surpress Logging
+
+    p.wait() # wait
+
+    # print("--- %s seconds ---" % (time.time() - start_time))
+    # print('------FINISHED PROGRAM on level:' + str(gridLevel))
+    processList.append(p)
+    ###Wait for all simulation subprocesses before proceeding
+    exit_codes = [p.wait() for p in processList]
+
+    return
+
+#-------------------------------------------------------------------------------------------------------
+def ReadEffectiveQuantities(QFilePath = os.path.dirname(os.getcwd()) + '/outputs/QMatrix.txt', BFilePath = os.path.dirname(os.getcwd())+ '/outputs/BMatrix.txt'):
+    # Read Output Matrices (effective quantities)
+    # From Cell-Problem output Files : ../outputs/Qmatrix.txt , ../outputs/Bmatrix.txt
+    # -- Read Matrix Qhom
+    X = []
+    # with codecs.open(path + '/outputs/QMatrix.txt', encoding='utf-8-sig') as f:
+    with codecs.open(QFilePath, encoding='utf-8-sig') as f:
+        for line in f:
+            s = line.split()
+            X.append([float(s[i]) for i in range(len(s))])
+    Q = np.array([[X[0][2], X[1][2], X[2][2]],
+                  [X[3][2], X[4][2], X[5][2]],
+                  [X[6][2], X[7][2], X[8][2]] ])
+
+    # -- Read Beff (as Vector)
+    X = []
+    # with codecs.open(path + '/outputs/BMatrix.txt', encoding='utf-8-sig') as f:
+    with codecs.open(BFilePath, encoding='utf-8-sig') as f:
+        for line in f:
+            s = line.split()
+            X.append([float(s[i]) for i in range(len(s))])
+    B = np.array([X[0][2], X[1][2], X[2][2]])
+    return Q, B
+
+
+# Function to convert a list to string
+def listToString(s):
+    # initialize an empty string
+    str1 = " "
+    # return string
+    # return (str1.join(str(s)))
+    return (str1.join(str(e) for e in s))
+
+
+def SetParametersCellProblem(mu_,lambda_,rho_,alpha,beta,theta,gamma,gridLevel, ParsetFilePath = os.path.dirname(os.getcwd()) +"/inputs/cellsolver.parset"):
+    print('----set Parameters -----')
+    with open(ParsetFilePath, 'r') as file:
+        filedata = file.read()
+    filedata = re.sub('(?m)^numLevels\s?=.*','numLevels='+str(gridLevel)+' '+str(gridLevel) ,filedata)
+    filedata = re.sub('(?m)^alpha\s?=.*','alpha='+str(alpha),filedata)
+    filedata = re.sub('(?m)^beta\s?=.*','beta='+str(beta),filedata)
+    filedata = re.sub('(?m)^theta\s?=.*','theta='+str(theta),filedata)
+    filedata = re.sub('(?m)^gamma\s?=.*','gamma='+str(gamma),filedata)
+    filedata = re.sub('(?m)^mu\s?=.*','mu='+listToString(mu_),filedata)
+    filedata = re.sub('(?m)^lambda\s?=.*','lambda='+listToString(lambda_),filedata)
+    filedata = re.sub('(?m)^rho\s?=.*','rho='+listToString(rho_),filedata)
+    filedata = re.sub('(?m)^mu1\s?=.*','mu1='+str(mu_[0]),filedata)
+    filedata = re.sub('(?m)^rho1\s?=.*','rho1='+str(rho_[0]),filedata)
+    filedata = re.sub('(?m)^lambda1\s?=.*','lambda1='+str(lambda_[0]),filedata)
+    f = open(ParsetFilePath,'w')
+    f.write(filedata)
+    f.close()
+
+#-------------------------------------------------------------------------------------------------------
+#-------------------------------------------------------------------------------------------------------

microstructure_testsuite/microstructure_testsuite.py

@@ -11,146 +11,109 @@ import matplotlib as mpl
 from matplotlib.ticker import MultipleLocator,FormatStrFormatter,MaxNLocator
 import codecs
 import sys
-# import sympy as sym
-# from prettytable import PrettyTable
-# import latextable
-# from texttable import Texttable
-# from tabulate import tabulate
+from helper_functions import *
+
 #-------------------------------------------------------------------------------------------------------
-def run_CellProblem(executable, parset, gridLevel, gamma, mu1,lambda1, rho1, alpha, beta, theta, material_prestrain_imp, outputPath, write_materialFunctions, write_prestrainFunctions, write_VTK,write_L2Error ,write_IntegralMean, write_LOG ):
-    print('----- RUN Cell-Problem ----')
-    processList = []
-    LOGFILE = "Cell-Problem_output.log"
-    print('LOGFILE:',LOGFILE)
-    print('executable:',executable)
-    print('parset:',parset)
+########################
+#### SET PARAMETERS ####
+########################
+# ----- Setup Paths -----
+# outputPath = ' ../outputs'
+outputPath = os.path.dirname(os.getcwd()) + '/outputs'
+QFilePath = os.path.dirname(os.getcwd()) + '/outputs/QMatrix.txt'
+BFilePath = os.path.dirname(os.getcwd())+ '/outputs/BMatrix.txt'
+ParsetFilePath = os.path.dirname(os.getcwd()) +"/inputs/cellsolver.parset"
+parset = ' ../inputs/cellsolver.parset'
+executable = ' ../build-cmake/src/Cell-Problem'
 
-        # start_time = time.time()
-    if write_LOG:
-        p = subprocess.Popen(executable + parset
-                                        # + " -numLevels " + str(gridLevel) + " " + str(gridLevel)
-                                        + " -gamma " + str(gamma)
-                                        + " -mu1 " + str(mu1)
-                                        + " -lambda1 " + str(lambda1)
-                                        + " -alpha " + str(alpha)
-                                        + " -beta " + str(beta)
-                                        + " -theta " + str(theta)
-                                        + " -material_prestrain_imp " + str(material_prestrain_imp)
-                                        + " -write_materialFunctions " + str(write_materialFunctions)
-                                        + " -outputPath " + str(outputPath)
-                                        + " | tee " + LOGFILE, shell=True)
+# ----- Define Input parameters  --------------------
+phases = 2 # number of phases
 
-    else:
-        p = subprocess.Popen(executable + parset
-                                        # + " -numLevels " + str(gridLevel) + " " + str(gridLevel)
-                                        + " -gamma " + str(gamma)
-                                        + " -mu1 " + str(mu1)
-                                        + " -lambda1 " + str(lambda1)
-                                        + " -alpha " + str(alpha)
-                                        + " -beta " + str(beta)
-                                        + " -theta " + str(theta)
-                                        + " -material_prestrain_imp " + str(material_prestrain_imp)
-                                        + " -write_materialFunctions " + str(write_materialFunctions)
-                                        + " -outputPath " + str(outputPath), stdout=subprocess.DEVNULL, stderr=subprocess.DEVNULL, shell=True )  # surpress Logging
+# --- Setup Material & Prestrain parameters:
+mu_ = np.array([1.0, 3.0])      # Vector of [mu1, mu2 , ....]
+lambda_ = np.array([1.0, 3.0])  # Vector of [lambda1, lambda2 , ....]
+rho_ = np.array([1.0, 8.0])     # Vector of [rho1, rho2 , ....]
 
-    p.wait() # wait
+assert phases == np.size(mu_)== np.size(lambda_) == np.size(rho_) , "Number of Phases must align with number of material and prestrain parameters!"
+
+
+gamma=1.0 #scale ratio
+
+# ---------------------------------------------------
 
-    # print("--- %s seconds ---" % (time.time() - start_time))
-    # print('------FINISHED PROGRAM on level:' + str(gridLevel))
-    processList.append(p)
-    ###Wait for all simulation subprocesses before proceeding
-    exit_codes = [p.wait() for p in processList]
 
-    return
 
-#-------------------------------------------------------------------------------------------------------
-def ReadEffectiveQuantities(QFilePath = os.path.dirname(os.getcwd()) + '/outputs/QMatrix.txt', BFilePath = os.path.dirname(os.getcwd())+ '/outputs/BMatrix.txt'):
-    # Read Output Matrices (effective quantities)
-    # From Cell-Problem output Files : ../outputs/Qmatrix.txt , ../outputs/Bmatrix.txt
-    # -- Read Matrix Qhom
-    X = []
-    # with codecs.open(path + '/outputs/QMatrix.txt', encoding='utf-8-sig') as f:
-    with codecs.open(QFilePath, encoding='utf-8-sig') as f:
-        for line in f:
-            s = line.split()
-            X.append([float(s[i]) for i in range(len(s))])
-    Q = np.array([[X[0][2], X[1][2], X[2][2]],
-                  [X[3][2], X[4][2], X[5][2]],
-                  [X[6][2], X[7][2], X[8][2]] ])
 
-    # -- Read Beff (as Vector)
-    X = []
-    # with codecs.open(path + '/outputs/BMatrix.txt', encoding='utf-8-sig') as f:
-    with codecs.open(BFilePath, encoding='utf-8-sig') as f:
-        for line in f:
-            s = line.split()
-            X.append([float(s[i]) for i in range(len(s))])
-    B = np.array([X[0][2], X[1][2], X[2][2]])
-    return Q, B
 
-#-------------------------------------------------------------------------------------------------------
-#-------------------------------------------------------------------------------------------------------
-########################
-#### SET PARAMETERS ####
-########################
-# ----- Setup Paths -----
-outputPath = " ../outputs"
-QFilePath = os.path.dirname(os.getcwd()) + '/outputs/QMatrix.txt'
-BFilePath = os.path.dirname(os.getcwd())+ '/outputs/BMatrix.txt'
 
-print('---- Input parameters: -----')
-mu1 = 10.0
-rho1 = 1.0
-alpha = 2.8
-beta = 2.0
-theta = 1.0/4.0
-gamma = 0.75
-lambda1= 10.0
 
-gridLevel = 4
+# --- Choose grid-Level for computation:
+gridLevel = 3
+
 #############################################
-#  Choose preferred Geometry/Prestrain/Material
+#  Choose preferred Geometry/Prestrain/Material             //TODO: Add Option for more Phases
 #############################################
-#material_prestrain_imp= "analytical_Example"
-material_prestrain_imp= "parametrized_Laminate"
-
+if phases == 1:
+    material_prestrain_imp= "homogeneous"
+elif phases == 2:
+    material_prestrain_imp= "two_phase_material_1"  #read as a python-function
+    # material_prestrain_imp= "two_phase_material_2"  #read as a python-function
+    # material_prestrain_imp= "two_phase_material_3"    #read as a python-function
+# elif phases== 3:  #//TODO:
+#     material_prestrain_imp= "three_phase_material"     #read as a python-function
 
 #############################################
 #  Solver Type: #1: CG - SOLVER (default), #2: GMRES - SOLVER, #3: QR - SOLVER
 #############################################
 Solvertype = 1 #(default = 1)
 Solver_verbosity = 0  #(default = 2)  degree of information for solver output
-
 set_IntegralZero = True
 
 #############################################
 #  Output-Options
 #############################################
-write_materialFunctions = False
-write_prestrainFunctions = False
+write_materialFunctions = True
+write_prestrainFunctions = True
 write_VTK = False
-# write_L2Error = True
-# write_IntegralMean = True
+
 write_L2Error = False
 write_IntegralMean = False
 
 write_LOG = False   # writes Cell-Problem output-LOG in "Cell-Problem_output.log"
-
+# write_LOG = True   # writes Cell-Problem output-LOG in "Cell-Problem_output.log"
+
+#---- Some of the old Material definitions use the following input parameters: (not needed when using "two_phase_material_x"):   --------
+mu1 = mu_[0]
+mu2 = mu_[1]
+lambda1 = lambda_[0]
+lambda2 = lambda_[1]
+rho1 = rho_[0]
+rho2 = rho_[1]
+beta = mu_[1]/mu_[0]   #ratio between material parameters
+alpha= rho_[1]/rho_[0] #prestrain-contrast
+theta = 1.0/4.0        #(volume fraction) needed for some of the older material definitions. Doesn't matter when using a indicatorFunction
+# ---------------------------------------------------
+print('---- Input parameters: -----')
 print('mu1: ', mu1)
+print('mu2: ', mu2)
+print('lambda1: ', lambda1)
+print('lambda2: ', lambda2)
 print('rho1: ', rho1)
+print('rho2: ', rho2)
 print('alpha: ', alpha)
 print('beta: ', beta)
 print('theta: ', theta)
 print('gamma:', gamma)
 print('material_prestrain_imp: ', material_prestrain_imp)
 print('gridLevel: ', gridLevel)
-
-parset = ' ../inputs/cellsolver.parset'
-executable = ' ../build-cmake/src/Cell-Problem'
 print('---------------------------------------------------------')
 ###########################################################################################################
 
+#Set Parameters
+SetParametersCellProblem(mu_,lambda_,rho_,alpha,beta,theta,gamma,gridLevel, ParsetFilePath)
 
+
+#Run Cell-Problem
 run_CellProblem(executable,
                 parset,
                 gridLevel,
@@ -171,367 +134,18 @@ run_CellProblem(executable,
                 write_LOG
                 )
 
+print('---------------------------------------------------------')
+#Read effective quantities
 print('Read effective quantities...')
 Q, B = ReadEffectiveQuantities(QFilePath,BFilePath)
 # Q, B = ReadEffectiveQuantities()
 print('Q:', Q)
+print('---------------------------------------------------------')
 print('B:', B)
 
 
 
-
-
-
-
-
-
-
-
-
-
-
-
-
-#
-#
-# for kappa in [4]:
-#
-#     print('--- RUN COMPUTATION FOR KAPPA=',kappa)
-#
-#     # --- map R2 into R2
-#     if domainDim == 2 and targetDim ==2:
-#         print('domainDim == 2 and targetDim == 2')
-#
-#         if kappa == 0 or kappa == 1:
-#             initialIterate = "initialIterateR2_R2"
-#         else:
-#             initialIterate = "initialIterateR2_R2_deg"+str(kappa)
-#         # initialIterate = "initialIterateR2_R2"
-#         # initialIterate = "initialIterateR2_R2_deg2"
-#         # initialIterate = "initialIterateR2_R2_deg3"
-#         # initialIterate = "initialIterateR2_R2_deg4"
-#         # initialIterate = "initialIterateR2_R2_deg5"
-#         # initialIterate = "initialIterateR2_R2_deg8"
-#
-#         executable = '../build-cmake/src/harmonicmaps_2d'
-#         parset = ' ../parsets/harmonicmaps_EOC_2d.parset'
-#
-#     # --- map R2 into R3
-#     if domainDim == 2 and targetDim ==3:
-#         print('domainDim == 2 and targetDim == 3')
-#
-#         # initialIterate = "initialIterateR2_R3"
-#         # initialIterate = "initialIterateR2_R3_smoothed"
-#         if kappa == 0 or kappa == 1:
-#             # initialIterate = "initialIterateR2_R3"
-#             initialIterate = "initialIterateR2_R3_smoothed"
-#         else:
-#             initialIterate = "initialIterateR2_R2_deg"+str(kappa)
-#
-#         # initialIterate = "initialIterateR2_R3_deg2"
-#         # initialIterate = "initialIterateR2_R3_deg3"
-#         # initialIterate = "initialIterateR2_R3_deg4"
-#         # initialIterate = "initialIterateR2_R3_deg5"
-#         executable = '../build-cmake/src/harmonicmaps_2d'
-#         parset = ' ../parsets/harmonicmaps_EOC_2d.parset'
-#
-#      # --- map R3 into R3
-#     if domainDim == 3 and targetDim ==3:
-#         print('domainDim == 3 and targetDim == 3')
-#
-#         if kappa == 0 or kappa == 1:
-#             initialIterate = "initialIterateR3_R3"
-#         else:
-#             initialIterate = "initialIterateR3_R3_deg"+str(kappa)
-#
-#
-#         # initialIterate = "initialIterateR3_R3_pertubed"
-#         # initialIterate = "initialIterateR2_R2"
-#         # initialIterate = "initialIterateR3_R3_deg2"
-#         # initialIterate = "initialIterateR3_R3_deg3"
-#         # initialIterate = "initialIterateR3_R3_deg4"
-#         # initialIterate = "initialIterateR3_R3_deg5"
-#         # initialIterate = "initialIterateR3_R3_deg8"
-#         # initialIterate = "initialIterateR3_R3_deg12"
-#         # initialIterate = "initialIterateR3_R3_deg3_Pert"
-#         # initialIterate = "initialIterateR3_R3_deg3_newOrigin"
-#
-#
-#         executable = '../build-cmake/src/harmonicmaps_3d'
-#         parset = ' ../parsets/harmonicmaps_EOC_3d.parset'
-#
-#
-#
-#
-#
-#
-#
-#     print(' InitialIterate used: ', initialIterate)
-#
-#     path = os.getcwd()
-#     print("Path: ", path)
-#
-#     for order in [1]:
-
-        # minLevel =  5 + 1 - order;
-        # maxLevel = 10 + 1 - order;
-        # minLevel =  2 + 1 - order;
-        # maxLevel = 8 + 1 - order;
-
-        #--- Setup LatexTable
-        # table_1 = Texttable()
-        # # table_1.maketitle = "Test"
-        # table_1.set_cols_align(["l", "r", "c"])
-        # table_1.set_cols_valign(["t", "m", "b"])
-        # table_1.add_rows([["Name", "Age", "Nickname"],
-        #                  ["Mr\nXavier\nHuon", 32, "Xav'"],
-        #                  ["Mr\nBaptiste\nClement", 1, "Baby"],
-        #                  ["Mme\nLouise\nBourgeau", 28, "Lou\n \nLoue"]])
-        #
-        # # table_1.add_rows("\multicolumn{3}{c}{ Projection-based-FE of order: 1, $\kappa = 2$}")
-        # print('-- Example 1: Basic --')
-        # print('Texttable Output:')
-        # print(table_1.draw())
-        # print('\nLatextable Output:')
-        # print(latextable.draw_latex(table_1, caption="An example table.", label="table:example_table"))
-        #
-
-
-
-        # rows = [['Rocket', 'Organisation', 'LEO Payload (Tonnes)', 'Maiden Flight'],
-        #         ['Saturn V', 'NASA', '140', '1967'],
-        #         ['Space Shuttle', 'NASA', '24.4', '1981'],
-        #         ['Falcon 9 FT-Expended', 'SpaceX', '22.8', '2017'],
-        #         ['Ariane 5 ECA', 'ESA', '21', '2002']]
-        #
-        # table = Texttable()
-        # table.set_cols_align(["c"] * 4)
-        # table.set_deco(Texttable.HEADER | Texttable.VLINES)
-        # table.add_rows(rows)
-        #
-        # print('Tabulate Table:')
-        # print(tabulate(rows, headers='firstrow'))
-        #
-        # # print('\nTexttable Table:')
-        # # print(table.draw())
-        #
-        # print('\nTabulate Latex:')
-        # print(tabulate(rows, headers='firstrow', tablefmt='latex'))
-
-        # print('\nTexttable Latex:')
-        # print(latextable.draw_latex(table, caption="A comparison of rocket features."))
-
-
-
-#
-#         #--- Setup Output-Table:
-#         x = PrettyTable()
-#         x.title = interpolationMethod + "-FE" + ' of order ' + str(order) + " into R" + str(targetDim)
-#         # x.field_names = ["k", "#Triang/#DOF", "# TR-Steps","$E^{Di}(0)$", "$E^{Di}$" , "$EOC_h^{L^2}$", "$EOC_h^{H^1}$", "$EOC_h^{E}$", "$\lnorm{u_h-u_{2h}$", "$\lnorm{u_{2h}-u_{4h}$", "$\seminorm{u_h-u_{2h}}$", "$\seminorm{u_{2h}-u_{4h}}$" ]
-#         x.field_names = ["k", "#Triang/#DOF", "# TR-Steps","$E^{Di}(0)$", "$E^{Di}$", "$\delta_1[u_h]$" , "$EOC_h^{L^2}$", "$EOC_h^{H^1}$", "$EOC_h^{E}$", "$\lnorm{u_h-u_{2h}$", "$\lnorm{u_{2h}-u_{4h}$", "$\seminorm{u_h-u_{2h}}$", "$\seminorm{u_{2h}-u_{4h}}$" , "wall-time" ]
-#         x.align["k"] = "l" # Left align city names
-#         x.padding_width = 1 # One space between column edges and contents (default)
-#
-#         rows = []
-#         rows.append(["k", "#Triang/#DOF", "# TR-Steps","$E^{Di}(0)$", "$E^{Di}$", "$\delta_1[u_h]$" , "$EOC_h^{L^2}$", "$EOC_h^{H^1}$", "$EOC_h^{E}$" , "wall-time" ])
-#         # rows.append(["k", "#Triang/#DOF", "# TR-Steps","wall-time","$E^{Di}(0)$", "$E^{Di}$" , "$EOC_h^{L^2}$", "$EOC_h^{H^1}$", "$EOC_h^{E}$" ])
-#         print('rows:', rows)
-#
-#         if compute_HarmonicMaps:
-#
-#             computeHarmonicMap(minLevel,maxLevel,domainDim,targetDim,order,interpolationMethod,maxTrustRegionSteps,initialIterate,randomInitialIterate,readConfiguration,configurationFile,pertubation,pertubationRadius,epsilon,tolerance,executable,parset  )
-#
-#         # --------------------
-#
-#
-#         if compute_Error:
-#             print("-------------------------------------------------------")
-#             print("Now measuring errors with discretizationErrorMode:" + str(discretizationErrorMode))
-#         # subprocess.call(["echo", "Now measuring errors"])
-#
-#             # EOC_l2_list = ["-","-"]
-#
-#             energyList = []
-#             initial_energyList = []
-#             stepList = []
-#             ndofList = []
-#             timeList = []
-#             # EOC_l2_list = []
-#
-#         #
-#             for numLevels in range(minLevel,maxLevel+1):
-#
-#                 print("NUMBER OF LEVELS:", numLevels)
-#                 # Measure the discretization errors against the solution on the finest grid
-#                 # LOGFILE = "./compute-disc-error_" + str(order) + "_" + str(numLevels) + ".log"
-#
-#
-#                 #subprocess.Popen(["../build-cmake/src/compute-disc-error", "compute-disc-error-skyrmions-hexagon.parset",
-#                                                 #"-order", str(order),
-#                                                 #"-numLevels", str(numLevels),
-#                                                 #"-numReferenceLevels", str(maxLevel),
-#                                                 #"-simulationData", "harmonicmaps-result-" + str(order) + "-" + str(numLevels) + ".data",
-#                                                 #"-referenceData", "harmonicmaps-result-" + str(order) + "-" + str(maxLevel) + ".data"])
-#
-#
-#                 # ------------  Get Energy / nDofs/Steps etc. -----------------------------------------------------------------------------
-#                 LOGFILE_comp = "./harmonicmapsR"+ str(domainDim) +"_intoR"+ str(targetDim) +  "_deg" + str(kappa) + "_"+ interpolationMethod + "_" + str(order) + "_" + str(numLevels) + ".log"
-#                 # Read Energy Values:
-#                 with open(LOGFILE_comp, 'r') as file:
-#                     output  = file.read()
-#
-#                 try:
-#                     # tmp_energy = re.search(r'(?m)energy: (-?\d\d?\.\d+[Ee]?[+\-]?\d\d?)',output).group(1)     # muss nun nichtmehr am Zeilenanfang stehen! :)
-#                     # tmp_energy = re.findall(r'(?m)energy: (-?\d\d?\.\d+[Ee]?[+\-]?\d\d?)',output).group(1)
-#                     # tmp_energy = re.findall(r'(?m)energy: (-?\d\d?\.\d+[Ee]?[+\-]?\d\d?)',output)
-#                     # tmp_energy = re.findall(r'(?m)energy: (-?\d\d?\d?\d?\.\d+[Ee]?[+\-]?\d\d?)',output)
-#                     tmp_energy = re.findall(r'(?m)energy: (-?\d?\d?\d?\.?\d+[Ee]?[+\-]?\d?\d?)',output)
-#                     tmp_step = re.findall(r'(?m)Step Number: (\d+)',output)
-#                     tmp_dofs = re.findall(r'(?m)powerBasis: (\d+)',output)
-#                     # tmp_time = re.findall(r'(?m)Time: (-?\d\d?\d?\.\d+[Ee]?[+\-]?\d\d?)',output)
-#                     # tmp_time = re.findall(r'(?m)Time: (-?\d+\.\d+[Ee]?[+\-]?\d\d?)',output)
-#                     tmp_time = re.findall(r'(?m)Time: (-?\d+\.\d+\d?)',output)
-#
-#                 except AttributeError:
-#                     # tmp_energy = re.search(r'(?m)energy: (-?\d\d?\.\d+[Ee]?[+\-]?\d\d?)',output)    # muss nun nichtmehr am Zeilenanfang stehen! :)
-#                     tmp_energy = re.search(r'(?m)energy: (-?\d?\d?\d?\.?\d+[Ee]?[+\-]?\d?\d?)',output)
-#                     tmp_step = re.findall(r'(?m)Step Number: (\d+)',output)
-#                     tmp_dofs = re.findall(r'(?m)powerBasis: (\d+)',output)
-#                     # tmp_time = re.findall(r'(?m)Time: (-?\d\d?\d?\d?\d?\.\d+[Ee]?[+\-]?\d\d?)',output)
-#                     # tmp_time = re.findall(r'(?m)Time: (-?\d+\.\d+[Ee]?[+\-]?\d\d?)',output)
-#                     tmp_time = re.findall(r'(?m)Time: (-?\d+\.\d+\d?)',output)
-#
-#                 if maxTrustRegionSteps>0 :
-#                     print('tmp_step:',tmp_step)
-#                     print('tmp_step last:', int(tmp_step[-1]))
-#
-#                     print('tmp_dofs:',tmp_dofs)
-#                     print('tmp_dofs last:', int(tmp_dofs[-1]))
-#
-#                     print('type tmp_energy:', type(tmp_energy))
-#                     print('tmp_energy:', tmp_energy)
-#                     print('tmp_energy last:', tmp_energy[-1])
-#
-#                     print('tmp_time:', tmp_time)
-#
-#                     energy = float(tmp_energy[-1])   #[1] since otherwise it recognizes "2" from L2error...
-#                     energyList.append(energy)
-#
-#                     initialEnergy = float(tmp_energy[0])
-#                     initial_energyList.append(initialEnergy)
-#
-#                     steps = int(tmp_step[-1])+1   #starts with zero therefore +1
-#                     stepList.append(steps)
-#
-#                     ndofs = int(tmp_dofs[-1])
-#                     ndofList.append(ndofs)
-#
-#                     time = float(tmp_time[-1])
-#                     timeList.append(time)
-#                 # -----------------------------------------------------------------------------------------------
-#
-#                 if numLevels >= (minLevel+2) and discretizationErrorMode=='discrete':
-#
-#                     [EOC_l2, EOC_h1, EOC_energy, L2_fine , L2_coarse , H1_fine , H1_coarse, constraintError] \
-#                         = computeError_discrete(numLevels,minLevel,domainDim,targetDim,order,interpolationMethod,discretizationErrorMode,maxTrustRegionSteps,energyList,ndofList,stepList,initial_energyList)
-#
-#                     # x.add_row([numLevels, currentDofs ,  currentSteps,currentInitialEnergy,  currentEnergy, EOC_l2, EOC_h1, EOC_energy , format(l2_fine, "5.3e"), format(l2_coarse, "5.3e"), format(h1_fine, "5.3e"), format(h1_coarse, "5.3e") ])
-#                     # x.add_row([numLevels, "-" , "-", "energy", EOC_l2, EOC_h1, "EOC_h^{E}" , format(l2_fine, "5.3e"), format(l2_coarse, "5.3e"), format(h1_fine, "5.3e"), format(h1_coarse, "5.3e") ])
-#
-#                 elif numLevels >= (minLevel+1)  and discretizationErrorMode=='analytical':
-#
-#                     [EOC_l2, EOC_h1, EOC_energy,L2_fine , L2_coarse , H1_fine , H1_coarse, constraintError] \
-#                         = computeError_analytical(numLevels,domainDim,targetDim,order,interpolationMethod,discretizationErrorMode,maxTrustRegionSteps,energyList,ndofList,stepList,initial_energyList,referenceSolution,numReferenceLevels)
-#
-#                 else:
-#                     EOC_l2 = "-"
-#                     EOC_h1 = "-"
-#                     EOC_energy = "-"
-#                     L2_fine = "-"
-#                     L2_coarse = "-"
-#                     H1_fine = "-"
-#                     H1_coarse = "-"
-#                     constraintError = "-"
-#
-#                     # if maxTrustRegionSteps > 0 :
-#                     #     print('energyList', energyList)
-#                     #     # currentEnergy = energyList[numLevels-1]
-#                     #     currentEnergy = energyList[-1]              #  better like this?
-#                     #     print('currentEnergy', currentEnergy)
-#                     #
-#                     #         # currentDofs = ndofList[numLevels-1]
-#                     #     currentDofs = ndofList[-1]
-#                     #     print('currentDofs', currentDofs)
-#                     #
-#                     #     # currentSteps = stepList[numLevels-1]
-#                     #     currentSteps = stepList[-1]
-#                     #     print('currentSteps', currentSteps)
-#                     #
-#                     #     print('initial_energyList', initial_energyList)
-#                     #     currentInitialEnergy = initial_energyList[-1]
-#                     #     print('currentInitialEnergy', currentInitialEnergy)
-#                     #
-#                     # else:
-#                     #     currentEnergy = "-"
-#                     #     currentDofs = "-"
-#                     #     currentSteps = "-"
-#                     #     currentInitialEnergy = "-"
-#
-#                     # x.add_row([numLevels, currentDofs ,  currentSteps, currentInitialEnergy,  currentEnergy, "-", "-" , "-" , "-", "-", "-", "-"])
-#                     # print('Energy List:', energyList)
-#
-#
-#                 if maxTrustRegionSteps > 0 :
-#                     print('energyList', energyList)
-#                     # currentEnergy = energyList[numLevels-1]
-#                     currentEnergy = energyList[-1]              #  better like this?
-#                     print('currentEnergy', currentEnergy)
-#
-#                         # currentDofs = ndofList[numLevels-1]
-#                     currentDofs = ndofList[-1]
-#                     print('currentDofs', currentDofs)
-#
-#                     # currentSteps = stepList[numLevels-1]
-#                     currentSteps = stepList[-1]
-#                     print('currentSteps', currentSteps)
-#
-#                     print('initial_energyList', initial_energyList)
-#                     currentInitialEnergy = initial_energyList[-1]
-#                     print('currentInitialEnergy', currentInitialEnergy)
-#
-#
-#
-#
-#                 else:
-#                     currentEnergy = "-"
-#                     currentDofs = "-"
-#                     currentSteps = "-"
-#                     currentInitialEnergy = "-"
-#
-#                 print('timeList:', timeList)
-#                 currentTime = timeList[-1]
-#                 print('currentSteps:', currentSteps)
-#
-#                 x.add_row([numLevels, currentDofs, currentSteps, currentInitialEnergy, currentEnergy,constraintError, EOC_l2, EOC_h1, EOC_energy,L2_fine , L2_coarse , H1_fine , H1_coarse , currentTime])
-#                 rows.append([numLevels, currentDofs, str(currentSteps), currentInitialEnergy, currentEnergy,constraintError, EOC_l2, EOC_h1, EOC_energy , currentTime])
-#
-#     ## Add extra column:
-#     # print(*EOC_l2_list)
-#     # x.add_column('EOC', EOC_l2_list)
-#
-#     # Write Table to Text-File:
-#     tablefile = open("EOC-table_R"+str(domainDim)+"intoR"+ str(targetDim)+ "_deg" + str(kappa) + "_"  + interpolationMethod+ "_order" +str(order) + "tol" + str(tolerance) + ".txt", "w")
-#     tablefile.write(str(x))
-#     tablefile.write('\n')
-#     tablefile.write(str(tabulate(rows, headers='firstrow', tablefmt='latex')))
-#
-#     tablefile.close()
-#     # print Table
-#     print(x)
-#
-#
-#
-# #--- Try Tabulate:
-#
-# print(tabulate(rows, headers='firstrow', tablefmt='latex'))
-# ##########  end of kappa-loop ##########################
+print('---------------------------------------------------------')
+print('access entries')
+print('Q[1][1]',Q[1][1])
+print('B[2]',B[2])

src/CMakeLists.txt

@@ -12,10 +12,10 @@ set(programs Cell-Problem
 foreach(_program ${programs})
   add_executable(${_program} ${_program}.cc)
   target_link_dune_default_libraries(${_program})
-  #add_dune_pythonlibs_flags(${_program})
-  #target_link_libraries(${_program} tinyxml2)
+  add_dune_pythonlibs_flags(${_program})
+  target_link_libraries(${_program} tinyxml2)
 #  target_compile_options(${_program} PRIVATE "-fpermissive")
 endforeach()
 
 
-set(CMAKE_BUILD_TYPE Debug)
+#set(CMAKE_BUILD_TYPE Debug)

src/Cell-Problem.cc

@@ -51,6 +51,9 @@
 
 #include <dune/functions/functionspacebases/hierarchicvectorwrapper.hh>
 
+#include <dune/fufem/dunepython.hh>
+#include <python2.7/Python.h>
+
 // #include <dune/fufem/discretizationerror.hh>
 // #include <boost/multiprecision/cpp_dec_float.hpp>
 // #include <boost/any.hpp>
@@ -878,6 +881,24 @@ int main(int argc, char *argv[])
 //   parameterSet.report(log); // short Alternativ
   
   
+    std::string geometryFunctionPath = parameterSet.get<std::string>("geometryFunctionPath");
+    //Start Python interpreter
+    Python::start();
+    Python::Reference main = Python::import("__main__");
+    Python::run("import math");
+
+    //"sys.path.append('/home/klaus/Desktop/DUNE/dune-gfe/problems')"
+    Python::runStream()
+        << std::endl << "import sys"
+        << std::endl << "sys.path.append('" << geometryFunctionPath << "')"
+        << std::endl;
+//         
+//     // Use python-function for initialIterate
+//     // Read initial iterate into a Python function
+//     Python::Module module = Python::import(parameterSet.get<std::string>("geometryFunction"));
+//     auto pythonInitialIterate = Python::make_function<double>(module.get("f"));
+  
+  
   constexpr int dim = 3;
   constexpr int dimWorld = 3;
 
@@ -949,6 +970,9 @@ int main(int argc, char *argv[])
   std::array<int,2> numLevels = parameterSet.get<std::array<int,2>>("numLevels", {1,3});
 
   int levelCounter = 0;
+   
+  
+//   FieldVector<double,2> mu = parameterSet.get<FieldVector<double,2>>("mu", {1.0,3.0});
 
   ///////////////////////////////////
   // Create Data Storage