diff --git a/dune.module b/dune.module
index 074b1c049229dd47967848355e35feea6c39a848..cbf9d40aef30499b5e199ea6d34cefd3b675fd8e 100644
--- a/dune.module
+++ b/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:
diff --git a/dune/microstructure/prestrain_material_geometry.hh b/dune/microstructure/prestrain_material_geometry.hh
index f96ac4878668bf8977f07afbd6e7b9dc36eff85e..d39bdbe88a2f1a21bdf6632bcdfb843dc0c1c235 100644
--- a/dune/microstructure/prestrain_material_geometry.hh
+++ b/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
{
diff --git a/geometries/__pycache__/two_phase_material.cpython-38.pyc b/geometries/__pycache__/two_phase_material.cpython-38.pyc
new file mode 100644
index 0000000000000000000000000000000000000000..e7624af2ace8ebb1cb9aa0900aead6aba77a25ba
Binary files /dev/null and b/geometries/__pycache__/two_phase_material.cpython-38.pyc differ
diff --git a/geometries/__pycache__/two_phase_material_1.cpython-38.pyc b/geometries/__pycache__/two_phase_material_1.cpython-38.pyc
new file mode 100644
index 0000000000000000000000000000000000000000..f1216220340d68114f39bfa32fc9fcac3375aee3
Binary files /dev/null and b/geometries/__pycache__/two_phase_material_1.cpython-38.pyc differ
diff --git a/geometries/__pycache__/two_phase_material_2.cpython-38.pyc b/geometries/__pycache__/two_phase_material_2.cpython-38.pyc
new file mode 100644
index 0000000000000000000000000000000000000000..de9026c11a1592fe0b109938ac6617277d4588af
Binary files /dev/null and b/geometries/__pycache__/two_phase_material_2.cpython-38.pyc differ
diff --git a/geometries/__pycache__/two_phase_material_3.cpython-38.pyc b/geometries/__pycache__/two_phase_material_3.cpython-38.pyc
new file mode 100644
index 0000000000000000000000000000000000000000..82d3ba4f63541c8c5816860928b1a814eaa272ba
Binary files /dev/null and b/geometries/__pycache__/two_phase_material_3.cpython-38.pyc differ
diff --git a/geometries/two_phase_material_1.py b/geometries/two_phase_material_1.py
new file mode 100644
index 0000000000000000000000000000000000000000..8c3418f48755232a24c143d57ffb1d3f842746cf
--- /dev/null
+++ b/geometries/two_phase_material_1.py
@@ -0,0 +1,13 @@
+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
diff --git a/geometries/two_phase_material_2.py b/geometries/two_phase_material_2.py
new file mode 100644
index 0000000000000000000000000000000000000000..b643b490f72d73dd77f9a72bd2d55a1e5fca9101
--- /dev/null
+++ b/geometries/two_phase_material_2.py
@@ -0,0 +1,13 @@
+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
diff --git a/geometries/two_phase_material_3.py b/geometries/two_phase_material_3.py
new file mode 100644
index 0000000000000000000000000000000000000000..fea31cb1e031f7a2850503c19e0ff22a42d49b0a
--- /dev/null
+++ b/geometries/two_phase_material_3.py
@@ -0,0 +1,13 @@
+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
diff --git a/inputs/cellsolver.parset b/inputs/cellsolver.parset
index 71dc1e54751e09a138e9d63845203c23ac82f321..1b3fb6f8bba60146ac1c2dc3a4d1a92b641d2665 100644
--- a/inputs/cellsolver.parset
+++ b/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
diff --git a/microstructure_testsuite/helper_functions.py b/microstructure_testsuite/helper_functions.py
new file mode 100644
index 0000000000000000000000000000000000000000..1e3e3c1ea65700ba2339b9841761ed5b5b7104a0
--- /dev/null
+++ b/microstructure_testsuite/helper_functions.py
@@ -0,0 +1,133 @@
+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()
+
+#-------------------------------------------------------------------------------------------------------
+#-------------------------------------------------------------------------------------------------------
diff --git a/microstructure_testsuite/microstructure_testsuite.py b/microstructure_testsuite/microstructure_testsuite.py
index 4cf31a637d831f249c534c4936bf413f9b1449df..000594440ff6f898f522bcafb3af2b3684f0bf26 100644
--- a/microstructure_testsuite/microstructure_testsuite.py
+++ b/microstructure_testsuite/microstructure_testsuite.py
@@ -11,109 +11,58 @@ 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!"
- # 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
+gamma=1.0 #scale ratio
+
+# ---------------------------------------------------
+
-#-------------------------------------------------------------------------------------------------------
-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= "analytical_Example"
+ #material_prestrain_imp= "parametrized_Laminate" #(Example used in the Paper)
+ 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
@@ -126,8 +75,8 @@ 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
@@ -135,22 +84,42 @@ 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 +140,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])
diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt
index abee48768ad5ffe02a85ea78fb487a08a5ae1ebd..d1e55346dd42b2d001b7657188e28011fd4f7ce0 100644
--- a/src/CMakeLists.txt
+++ b/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)
diff --git a/src/Cell-Problem.cc b/src/Cell-Problem.cc
index 3813a62be3e2aee8d2740682ed797c8e06f37174..348371e34edccaf7bb7e63853165062ea5dfd332 100644
--- a/src/Cell-Problem.cc
+++ b/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