Add perforated example with square-shaped holes

Plot-Scripts/PerforatedBilayer_square_parameterDependence.py

+import numpy as np
+import matplotlib.pyplot as plt
+import math
+import os
+import subprocess
+import fileinput
+import re
+import sys
+import matplotlib as mpl
+from mpl_toolkits.mplot3d import Axes3D
+import matplotlib.cm as cm
+import matplotlib.ticker as ticker
+from matplotlib.ticker import MultipleLocator,FormatStrFormatter,MaxNLocator
+import seaborn as sns
+import matplotlib.colors as mcolors
+import codecs
+import json
+# ----------------------------------------------------
+# --- Define Plot style:
+# plt.style.use("seaborn-white")
+# plt.style.use("seaborn-pastel")
+# plt.style.use("seaborn-colorblind")
+plt.style.use("seaborn")
+mpl.rcParams['text.usetex'] = True
+mpl.rcParams["font.family"] = "serif"
+mpl.rcParams["font.size"] = "14"
+mpl.rcParams['xtick.bottom'] = True
+mpl.rcParams['xtick.major.size'] = 3
+mpl.rcParams['xtick.minor.size'] = 1.5
+mpl.rcParams['xtick.major.width'] = 0.75
+mpl.rcParams['ytick.left'] = True
+mpl.rcParams['ytick.major.size'] = 3
+mpl.rcParams['ytick.minor.size'] = 1.5
+mpl.rcParams['ytick.major.width'] = 0.75
+
+#Adjust grid:
+mpl.rcParams.update({"axes.grid" : True}) # Add grid
+mpl.rcParams['axes.labelpad'] = 3
+mpl.rcParams['grid.linewidth'] = 0.25
+mpl.rcParams['grid.alpha'] = 0.9 # 0.75
+mpl.rcParams['grid.linestyle'] = '-'
+mpl.rcParams['grid.color']   = 'gray'#'black'
+mpl.rcParams['text.latex.preamble'] = r'\usepackage{amsfonts}' # Makes Use of \mathbb possible.
+# ----------------------------------------------------------------------------------------
+width = 5.79
+height = width / 1.618 # The golden ratio.
+
+
+
+dataset_numbers = [0, 1, 2, 3, 4 ,5]
+
+for dataset_number in dataset_numbers:
+
+    fig, ax = plt.subplots(figsize=(width,height))
+
+
+    Path_KappaUpper = './experiment/perforated-bilayer_square/results_upper_' + str(dataset_number) + '/kappa_simulation.txt' 
+    Path_KappaLower = './experiment/perforated-bilayer_square/results_lower_' + str(dataset_number) + '/kappa_simulation.txt' 
+    
+
+    kappa_sim_upper = open(Path_KappaUpper).read().split(",") 
+    kappa_sim_lower = open(Path_KappaLower).read().split(",") 
+    kappa_sim_upper = [float(i) for i in kappa_sim_upper]
+    kappa_sim_lower = [float(i) for i in kappa_sim_lower]
+
+
+    input = [0 , 0.05, 0.10, 0.15, 0.20, 0.25, 0.30] # Design parameter (volume ratio)
+
+
+    # --------------- Plot Lines + Scatter -----------------------
+    line_1 = ax.plot(input, kappa_sim_upper,                    # data
+                #  color='forestgreen',              # linecolor
+                marker='D',                         # each marker will be rendered as a circle
+                markersize=5,                       # marker size
+                #   markerfacecolor='darkorange',      # marker facecolor
+                markeredgecolor='black',            # marker edgecolor
+                markeredgewidth=0.75,                  # marker edge width
+                # linestyle='dashdot',              # line style will be dash line
+                linewidth=1.5,                      # line width
+                zorder=3,
+                label = r"$\kappa_{sim}$ (upper=passive perf.)")
+
+    line_2 = ax.plot(input, kappa_sim_lower,                     # data
+                # color='orangered',                # linecolor
+                marker='o',                         # each marker will be rendered as a circle
+                markersize=5,                       # marker size
+                #  markerfacecolor='cornflowerblue',   # marker facecolor
+                markeredgecolor='black',            # marker edgecolor
+                markeredgewidth=0.75,                  # marker edge width
+                # linestyle='--',                   # line style will be dash line
+                linewidth=1.5,                      # line width
+                zorder=3,
+                alpha=0.8,                           # Change opacity
+                label = r"$\kappa_{sim}$ (lower=active perf.)")
+
+
+    # --------------- Set Axes  -----------------------
+    # Plot - Title
+    match dataset_number:
+        case 0:
+            ax.set_title(r"ratio $r = 0.12$") 
+        case 1: 
+            ax.set_title(r"ratio $r = 0.17$") 
+        case 2:
+            ax.set_title(r"ratio $r = 0.22$") 
+        case 3: 
+            ax.set_title(r"ratio $r = 0.34$") 
+        case 4:
+            ax.set_title(r"ratio $r = 0.43$") 
+        case 5: 
+            ax.set_title(r"ratio $r = 0.49$") 
+    
+    ax.set_xlabel(r"Volume ratio ", labelpad=4)
+    ax.set_ylabel(r"Curvature $\kappa$", labelpad=4)
+    plt.tight_layout()
+
+
+    # --- Set Legend
+    legend = ax.legend()
+
+    frame = legend.get_frame()
+    frame.set_edgecolor('black')
+
+
+    # --- Adjust left/right spacing:
+    # plt.subplots_adjust(right=0.81)
+    # plt.subplots_adjust(left=0.11)
+
+    # ---------- Output Figure as pdf:
+    fig.set_size_inches(width, height)
+    fig.savefig('PerforatedBilayer_square_dependence_' + str(dataset_number) + '.pdf')
+    # plt.show()
\ No newline at end of file

experiment/perforated-bilayer_square/PolarPlotLocalEnergy.py

+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+"""
+Created on Wed Jul  6 13:17:28 2022
+
+@author: stefan
+"""
+import numpy as np
+import matplotlib.pyplot as plt
+import matplotlib.colors as colors
+import codecs
+import re
+import json
+
+def energy(kappa,alpha,Q,B)  :
+    G=kappa*np.array([[np.cos(alpha)**2],[np.sin(alpha)**2],[np.sqrt(2)*np.cos(alpha)*np.sin(alpha)]])-B
+    return np.matmul(np.transpose(G),np.matmul(Q,G))[0,0]
+
+def xytokappaalpha(x,y):
+   
+    if y>0:
+        return [np.sqrt(x**2+y**2), np.abs(np.arctan2(y,x))]
+    else:
+        return [-np.sqrt(x**2+y**2), np.abs(np.arctan2(y,x))]
+
+# Read effective quantites
+def ReadEffectiveQuantities(QFilePath, BFilePath):
+    # 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
+
+# Number of experiments / folders
+number=7
+show_plot = False
+
+#--- Choose wether to perforate upper (passive) or lower (active) layer
+perforatedLayer = 'upper'
+# perforatedLayer = 'lower'
+
+dataset_numbers = [0, 1, 2, 3, 4, 5]
+# dataset_numbers = [0]
+for dataset_number in dataset_numbers:
+
+
+    kappa=np.zeros(number)
+    for n in range(0,number):
+        #   Read from Date
+        print(str(n))
+        DataPath = './experiment/perforated-bilayer_square/results_'  +  perforatedLayer + '_' + str(dataset_number) + '/' +str(n)
+        QFilePath = DataPath + '/QMatrix.txt'
+        BFilePath = DataPath + '/BMatrix.txt'
+        ParameterPath = DataPath + '/parameter.txt'
+
+        # Read Thickness from parameter file (needed for energy scaling)
+        with open(ParameterPath , 'r') as file:
+            parameterFile  = file.read()
+        thickness = float(re.findall(r'(?m)h = (\d?\d?\d?\.?\d+[Ee]?[+\-]?\d?\d?)',parameterFile)[0])
+
+        Q, B = ReadEffectiveQuantities(QFilePath,BFilePath)
+        # Q=0.5*(np.transpose(Q)+Q) # symmetrize
+        B=np.transpose([B])
+        # 
+        
+        N=500
+        length=5
+        r, theta = np.meshgrid(np.linspace(0,length,N),np.radians(np.linspace(0, 360, N)))
+        E=np.zeros(np.shape(r))
+        for i in range(0,N): 
+            for j in range(0,N):     
+                if theta[i,j]<np.pi:
+                    E[i,j]=energy(r[i,j],theta[i,j],Q,B)  * (thickness**2)
+                else:
+                    E[i,j]=energy(-r[i,j],theta[i,j],Q,B) * (thickness**2)
+                
+        # Compute Minimizer
+        [imin,jmin]=np.unravel_index(E.argmin(),(N,N))
+        kappamin=r[imin,jmin]
+        alphamin=theta[imin,jmin]
+        kappa[n]=kappamin
+        fig, ax = plt.subplots(figsize=(6,6),subplot_kw=dict(projection='polar'))
+        levs=np.geomspace(E.min(),E.max(),400)
+        pcm=ax.contourf(theta, r, E, levs, norm=colors.PowerNorm(gamma=0.2), cmap='brg')
+        ax.set_xticks([0,np.pi/2])
+        ax.set_yticks([kappamin])
+        colorbarticks=np.linspace(E.min(),E.max(),6)
+        cbar = plt.colorbar(pcm, extend='max', ticks=colorbarticks, pad=0.1)
+        cbar.ax.tick_params(labelsize=6)
+        if (show_plot):
+            plt.show()
+        # Save Figure as .pdf
+        width = 5.79 
+        height = width / 1.618 # The golden ratio.
+        fig.set_size_inches(width, height)
+        fig.savefig('PerforatedBilayer_dataset_' +str(dataset_number) + '_exp' +str(n) + '.pdf')
+
+
+    f = open("./experiment/perforated-bilayer_square/results_" +  perforatedLayer + '_' + str(dataset_number) +  "/kappa_simulation.txt", "w")
+    f.write(str(kappa.tolist())[1:-1])     
+    f.close()   
\ No newline at end of file

experiment/perforated-bilayer_square/perfBilayer_test.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
+import threading
+
+# # Schreibe input datei für Parameter
+# def SetParametersCellProblem(ParameterSet, ParsetFilePath, outputPath):
+#     print('----set Parameters -----')
+#     with open(ParsetFilePath, 'r') as file:
+#         filedata = file.read()
+#         filedata = re.sub('(?m)^materialFunction\s?=.*','materialFunction = '+str(ParameterSet.materialFunction),filedata)
+#         filedata = re.sub('(?m)^gamma\s?=.*','gamma='+str(ParameterSet.gamma),filedata)
+#         filedata = re.sub('(?m)^numLevels\s?=.*','numLevels='+str(ParameterSet.numLevels)+' '+str(ParameterSet.numLevels) ,filedata)
+#         filedata = re.sub('(?m)^outputPath\s?=\s?.*','outputPath='+str(outputPath),filedata)
+#         f = open(ParsetFilePath,'w')
+#         f.write(filedata)
+#         f.close()
+
+
+# # Ändere Parameter der MaterialFunction
+# def SetParameterMaterialFunction(materialFunction, parameterName, parameterValue):
+#     with open(Path+"/"+materialFunction+'.py', 'r') as file:
+#         filedata = file.read()
+#         filedata = re.sub('(?m)^'+str(parameterName)+'\s?=.*',str(parameterName)+' = '+str(parameterValue),filedata)
+#         f = open(Path+"/"+materialFunction+'.py','w')
+#         f.write(filedata)
+#         f.close()
+
+def SetParameterMaterialFunction(inputFunction, parameterName, parameterValue):
+    with open(inputFunction+'.py', 'r') as file:
+        filedata = file.read()
+        filedata = re.sub('(?m)^'+str(parameterName)+'\s?=.*',str(parameterName)+' = '+str(parameterValue),filedata)
+        f = open(inputFunction+'.py','w')
+        f.write(filedata)
+        f.close()
+
+
+
+# # Rufe Programm zum Lösen des Cell-Problems auf
+# def run_CellProblem(executable, parset,write_LOG):
+#     print('----- RUN Cell-Problem ----')
+#     processList = []
+#     LOGFILE = "Cell-Problem_output.log"
+#     print('LOGFILE:',LOGFILE)
+#     print('executable:',executable)
+#     if write_LOG:
+#         p = subprocess.Popen(executable + parset
+#                                         + " | tee " + LOGFILE, shell=True)
+
+#     else:
+#         p = subprocess.Popen(executable + parset
+#                                         + " | tee " + LOGFILE, shell=True)
+
+#     p.wait() # wait
+#     processList.append(p)
+#     exit_codes = [p.wait() for p in processList]
+
+#     return
+
+# Read effective quantites
+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 for evaluating the energy in terms of kappa, alpha and Q, B
+def eval_energy(kappa,alpha,Q,B)  :
+    G=kappa*np.array([[np.cos(alpha)**2],[np.sin(alpha)**2],[np.sqrt(2)*np.cos(alpha)*np.sin(alpha)]])-B
+    return np.matmul(np.transpose(G),np.matmul(Q,G))[0,0]
+
+#-------------------------------------------------------------------------------------------------------
+########################
+#### SET PARAMETERS ####
+########################
+# ----- Setup Paths -----
+# write_LOG = True   # writes Cell-Problem output-LOG in "Cell-Problem_output.log"
+# path='/home/klaus/Desktop/Dune_release/dune-microstructure/experiment/wood-bilayer/results/'  
+# pythonPath = '/home/klaus/Desktop/Dune_release/dune-microstructure/experiment/wood-bilayer'
+
+#--- Choose wether to perforate upper (passive) or lower (active) layer
+# perforatedLayer = 'upper'
+perforatedLayer = 'lower'
+
+# dataset_numbers = [0, 1, 2, 3, 4, 5]
+dataset_numbers = [0]
+
+for dataset_number in dataset_numbers:
+    print("------------------")
+    print(str(dataset_number) + "th data set")
+    print("------------------")
+
+
+    # path = os.getcwd() + '/experiment/perforated-bilayer/results_' +  perforatedLayer + '/'
+    path = os.getcwd() + '/experiment/perforated-bilayer_square/results_' +  perforatedLayer + '_' + str(dataset_number) + '/'
+    pythonPath = os.getcwd() + '/experiment/perforated-bilayer_square'
+    pythonModule = "perforated_wood_" + perforatedLayer
+    executable = os.getcwd() + '/build-cmake/src/Cell-Problem'
+    # ---------------------------------
+    # Setup Experiment
+    # ---------------------------------
+    gamma = 1.0
+
+    # ----- Define Parameters for Material Function  --------------------
+    # [r, h, omega_flat, omega_target, theta, experimental_kappa]
+    # r = (thickness upper layer)/(thickness)
+    # h = thickness [meter]
+    # omega_flat = moisture content in the flat state before drying [%]
+    # omega_target = moisture content in the target state [%]
+    # theta = rotation angle (not implemented and used)
+    # beta = design parameter for perforation = ratio of Volume of (cylindrical) perforation to Volume of active/passive layer 
+
+    #Experiment: Perforate "active"/"passive" bilayer phase 
+
+
+    materialFunctionParameter=[
+    [  # Dataset Ratio r = 0.12
+    [0.12, 0.0047, 17.17547062, 8.959564147, 0.0, 0.0 ],
+    [0.12, 0.0047, 17.17547062, 8.959564147, 0.0, 0.05],
+    [0.12, 0.0047, 17.17547062, 8.959564147, 0.0, 0.1 ],
+    [0.12, 0.0047, 17.17547062, 8.959564147, 0.0, 0.15],
+    [0.12, 0.0047, 17.17547062, 8.959564147, 0.0, 0.2 ],
+    [0.12, 0.0047, 17.17547062, 8.959564147, 0.0, 0.25],
+    [0.12, 0.0047, 17.17547062, 8.959564147, 0.0, 0.3 ]
+    ],
+    [  # Dataset Ratio r = 0.17
+    [0.17, 0.0049, 17.17547062, 8.959564147, 0.0, 0.0 ],
+    [0.17, 0.0049, 17.17547062, 8.959564147, 0.0, 0.05],
+    [0.17, 0.0049, 17.17547062, 8.959564147, 0.0, 0.1 ],
+    [0.17, 0.0049, 17.17547062, 8.959564147, 0.0, 0.15],
+    [0.17, 0.0049, 17.17547062, 8.959564147, 0.0, 0.2 ],
+    [0.17, 0.0049, 17.17547062, 8.959564147, 0.0, 0.25],
+    [0.17, 0.0049, 17.17547062, 8.959564147, 0.0, 0.3 ]
+    ],
+    [ # Dataset Ratio r = 0.22
+    [0.22, 0.0053,  17.17547062, 8.959564147, 0.0, 0.0 ],
+    [0.22, 0.0053,  17.17547062, 8.959564147, 0.0, 0.05 ],
+    [0.22, 0.0053,  17.17547062, 8.959564147, 0.0, 0.1 ],
+    [0.22, 0.0053,  17.17547062, 8.959564147, 0.0, 0.15 ],
+    [0.22, 0.0053,  17.17547062, 8.959564147, 0.0, 0.2 ],
+    [0.22, 0.0053,  17.17547062, 8.959564147, 0.0, 0.25 ],
+    [0.22, 0.0053,  17.17547062, 8.959564147, 0.0, 0.3 ]
+    ],
+    [  # Dataset Ratio r = 0.34
+    [0.34, 0.0063, 17.17547062, 8.959564147, 0.0, 0.0 ],
+    [0.34, 0.0063, 17.17547062, 8.959564147, 0.0, 0.05],
+    [0.34, 0.0063, 17.17547062, 8.959564147, 0.0, 0.1 ],
+    [0.34, 0.0063, 17.17547062, 8.959564147, 0.0, 0.15],
+    [0.34, 0.0063, 17.17547062, 8.959564147, 0.0, 0.2 ],
+    [0.34, 0.0063, 17.17547062, 8.959564147, 0.0, 0.25],
+    [0.34, 0.0063, 17.17547062, 8.959564147, 0.0, 0.3 ]
+    ],
+    [  # Dataset Ratio r = 0.43
+    [0.43, 0.0073, 17.17547062, 8.959564147, 0.0, 0.0 ],
+    [0.43, 0.0073, 17.17547062, 8.959564147, 0.0, 0.05],
+    [0.43, 0.0073, 17.17547062, 8.959564147, 0.0, 0.1 ],
+    [0.43, 0.0073, 17.17547062, 8.959564147, 0.0, 0.15],
+    [0.43, 0.0073, 17.17547062, 8.959564147, 0.0, 0.2 ],
+    [0.43, 0.0073, 17.17547062, 8.959564147, 0.0, 0.25],
+    [0.43, 0.0073, 17.17547062, 8.959564147, 0.0, 0.3 ]
+    ],
+    [  # Dataset Ratio r = 0.49
+    [0.49, 0.008,  17.17547062, 8.959564147, 0.0, 0.0 ],
+    [0.49, 0.008,  17.17547062, 8.959564147, 0.0, 0.05],
+    [0.49, 0.008,  17.17547062, 8.959564147, 0.0, 0.1 ],
+    [0.49, 0.008,  17.17547062, 8.959564147, 0.0, 0.15],
+    [0.49, 0.008,  17.17547062, 8.959564147, 0.0, 0.2 ],
+    [0.49, 0.008,  17.17547062, 8.959564147, 0.0, 0.25],
+    [0.49, 0.008,  17.17547062, 8.959564147, 0.0, 0.3 ]
+    ]
+    ]
+    
+
+
+    #--- Different moisture values for different thicknesses:
+
+    # materialFunctionParameter=[
+    # [  # Dataset Ratio r = 0.12
+    # [0.12, 0.0047, 17.32986047, 14.70179844, 0.0, 0.0 ],
+    # [0.12, 0.0047, 17.32986047, 13.6246,     0.0, 0.05],
+    # [0.12, 0.0047, 17.32986047, 12.42994508, 0.0, 0.1 ],
+    # [0.12, 0.0047, 17.32986047, 11.69773413, 0.0, 0.15],
+    # [0.12, 0.0047, 17.32986047, 11.14159987, 0.0, 0.2 ],
+    # [0.12, 0.0047, 17.32986047, 9.500670278, 0.0, 0.25],
+    # [0.12, 0.0047, 17.32986047, 9.005046347, 0.0, 0.3 ]
+    # ],
+    # [  # Dataset Ratio r = 0.17
+    # [0.17, 0.0049, 17.28772791 , 14.75453569, 0.0, 0.0 ],
+    # [0.17, 0.0049, 17.28772791 , 13.71227639, 0.0, 0.05],
+    # [0.17, 0.0049, 17.28772791 , 12.54975012, 0.0, 0.1 ],
+    # [0.17, 0.0049, 17.28772791 , 11.83455959, 0.0, 0.15],
+    # [0.17, 0.0049, 17.28772791 , 11.29089521, 0.0, 0.2 ],
+    # [0.17, 0.0049, 17.28772791 , 9.620608917, 0.0, 0.25],
+    # [0.17, 0.0049, 17.28772791 , 9.101671742, 0.0, 0.3 ]
+    # ],
+    # [ # Dataset Ratio r = 0.22
+    # [0.22, 0.0053,  17.17547062, 8.959564147, 0.0, 0.0 ],
+    # [0.22, 0.0053,  17.17547062, 8.959564147, 0.0, 0.05 ],
+    # [0.22, 0.0053,  17.17547062, 8.959564147, 0.0, 0.1 ],
+    # [0.22, 0.0053,  17.17547062, 8.959564147, 0.0, 0.15 ],
+    # [0.22, 0.0053,  17.17547062, 8.959564147, 0.0, 0.2 ],
+    # [0.22, 0.0053,  17.17547062, 8.959564147, 0.0, 0.25 ],
+    # [0.22, 0.0053,  17.17547062, 8.959564147, 0.0, 0.3 ]
+    # ],
+    # [  # Dataset Ratio r = 0.34
+    # [0.34, 0.0063, 17.14061081 , 14.98380876, 0.0, 0.0 ],
+    # [0.34, 0.0063, 17.14061081 , 13.97154915, 0.0, 0.05],
+    # [0.34, 0.0063, 17.14061081 , 12.77309253, 0.0, 0.1 ],
+    # [0.34, 0.0063, 17.14061081 , 12.00959929, 0.0, 0.15],
+    # [0.34, 0.0063, 17.14061081 , 11.42001731, 0.0, 0.2 ],
+    # [0.34, 0.0063, 17.14061081 , 9.561447179, 0.0, 0.25],
+    # [0.34, 0.0063, 17.14061081 , 8.964704969, 0.0, 0.3 ]
+    # ],
+    # [  # Dataset Ratio r = 0.43
+    # [0.43, 0.0073, 17.07559686 , 15.11316339, 0.0, 0.0 ],
+    # [0.43, 0.0073, 17.07559686 , 14.17997082, 0.0, 0.05],
+    # [0.43, 0.0073, 17.07559686 , 13.05739844, 0.0, 0.1 ],
+    # [0.43, 0.0073, 17.07559686 , 12.32309209, 0.0, 0.15],
+    # [0.43, 0.0073, 17.07559686 , 11.74608518, 0.0, 0.2 ],
+    # [0.43, 0.0073, 17.07559686 , 9.812372466, 0.0, 0.25],
+    # [0.43, 0.0073, 17.07559686 , 9.10519385 , 0.0, 0.3 ]
+    # ],
+    # [  # Dataset Ratio r = 0.49
+    # [0.49, 0.008,  17.01520754, 15.30614414, 0.0, 0.0 ],
+    # [0.49, 0.008,  17.01520754, 14.49463867, 0.0, 0.05],
+    # [0.49, 0.008,  17.01520754, 13.46629742, 0.0, 0.1 ],
+    # [0.49, 0.008,  17.01520754, 12.78388234, 0.0, 0.15],
+    # [0.49, 0.008,  17.01520754, 12.23057715, 0.0, 0.2 ],
+    # [0.49, 0.008,  17.01520754, 10.21852839, 0.0, 0.25],
+    # [0.49, 0.008,  17.01520754, 9.341730605, 0.0, 0.3 ]
+    # ]
+    # ]
+
+
+
+    # ------ Loops through Parameters for Material Function -----------
+    for i in range(0,np.shape(materialFunctionParameter)[1]):
+        print("------------------")
+        print("New Loop")
+        print("------------------")
+    # Check output directory
+        outputPath = path + str(i)
+        isExist = os.path.exists(outputPath)
+        if not isExist:
+            # Create a new directory because it does not exist
+            os.makedirs(outputPath)
+            print("The new directory " + outputPath + " is created!")
+
+        # thread = threading.Thread(target=run_CellProblem(executable, pythonModule, pythonPath, LOGFILE))
+        # thread.start()
+
+        #TODO: apperently its not possible to pass a variable via subprocess and "calculate" another input value inside the python file.
+        #      Therefore we use this instead.
+        SetParameterMaterialFunction(pythonPath + "/" + pythonModule, "param_r",materialFunctionParameter[dataset_number][i][0])
+        SetParameterMaterialFunction(pythonPath + "/" + pythonModule, "param_h",materialFunctionParameter[dataset_number][i][1])
+        SetParameterMaterialFunction(pythonPath + "/" + pythonModule, "param_omega_flat",materialFunctionParameter[dataset_number][i][2])
+        SetParameterMaterialFunction(pythonPath + "/" + pythonModule, "param_omega_target",materialFunctionParameter[dataset_number][i][3])
+        SetParameterMaterialFunction(pythonPath + "/" + pythonModule, "param_theta",materialFunctionParameter[dataset_number][i][4])   
+        SetParameterMaterialFunction(pythonPath + "/" + pythonModule, "param_beta",materialFunctionParameter[dataset_number][i][5])    
+        # SetParameterMaterialFunction(pythonPath + "/" + pythonModule, "perfDepth",materialFunctionParameter[dataset_number][i][0])   
+        LOGFILE = outputPath + "/" + pythonModule + "_output" + "_" + str(i) + ".log"
+
+        processList = []
+        p = subprocess.Popen(executable + " " + pythonPath + " " + pythonModule
+                                        + " -outputPath " + outputPath
+                                        + " -gamma " + str(gamma) 
+                                        + " | tee " + LOGFILE, shell=True)
+
+        # p = subprocess.Popen(executable + " " + pythonPath + " " + pythonModule
+        #                                 + " -outputPath " + outputPath
+        #                                 + " -gamma " + str(gamma) 
+        #                                 + " -param_r " + str(materialFunctionParameter[i][0])
+        #                                 + " -param_h " + str(materialFunctionParameter[i][1])
+        #                                 + " -param_omega_flat " + str(materialFunctionParameter[i][2])
+        #                                 + " -param_omega_target " + str(materialFunctionParameter[i][3])
+        #                                 + " -phase2_angle " + str(materialFunctionParameter[i][4])
+        #                                 + " | tee " + LOGFILE, shell=True)
+
+        p.wait() # wait
+        processList.append(p)
+        exit_codes = [p.wait() for p in processList]
+        # ---------------------------------------------------
+        # wait here for the result to be available before continuing
+        # thread.join()
+        f = open(outputPath+"/parameter.txt", "w")
+        f.write("r = "+str(materialFunctionParameter[dataset_number][i][0])+"\n")
+        f.write("h = "+str(materialFunctionParameter[dataset_number][i][1])+"\n")
+        f.write("omega_flat = "+str(materialFunctionParameter[dataset_number][i][2])+"\n")        
+        f.write("omega_target = "+str(materialFunctionParameter[dataset_number][i][3])+"\n")     
+        f.write("param_beta = "+str(materialFunctionParameter[dataset_number][i][5])+"\n")         
+        f.close()   
+        #

experiment/perforated-bilayer_square/perforated_wood_lower.py

+import math
+#from python_matrix_operations import *
+import ctypes
+import os
+import sys
+import numpy as np
+# import elasticity_toolbox as elast
+
+class ParameterSet(dict):
+    def __init__(self, *args, **kwargs):
+        super(ParameterSet, self).__init__(*args, **kwargs)
+        self.__dict__ = self
+
+parameterSet = ParameterSet()
+#---------------------------------------------------------------
+#############################################
+#  Paths
+#############################################
+# Path for results and logfile
+parameterSet.outputPath='/home/klaus/Desktop/Dune_release/dune-microstructure/experiment/perforated-bilayer_square/results'
+parameterSet.baseName= 'perforated_wood_lower'   #(needed for Output-Filename)
+
+# Path for material description
+# parameterSet.geometryFunctionPath =experiment/wood-bilayer/
+
+#---------------------------------------------------------------
+# Wooden bilayer, https://www.ncbi.nlm.nih.gov/pmc/articles/PMC6191116/#pone.0205607.ref015
+#--- define indicator function
+# x[0] : y1-component -1/2 to 1/2
+# x[1] : y2-component -1/2 to 1/2
+# x[2] : x3-component range -1/2 to 1/2
+#--- define indicator function
+def indicatorFunction(x):
+    pRadius = math.sqrt((param_beta*(1.0-param_r))/(4.0*perfDepth))
+    if (x[2]>=(0.5-param_r)):
+        return 1      #Phase1
+    else :  
+        if( (max(abs(x[0]), abs(x[1])) < pRadius) and (x[2] <= (-0.5+perfDepth)) ): #inside perforation    
+            return 3  #Phase3
+        else:  
+            return 2  #Phase2
+    
+# def indicatorFunction(x):
+#     factor=1
+#     pRadius = 0.25
+#     if (x[2]>=(0.5-param_r) and np.sqrt(x[0]**2 + x[1]**2) < pRadius):
+#         return 3
+#     elif((x[2]>=(0.5-param_r))):  
+#         return 1  #Phase1
+#     else :
+#         return 2      #Phase2
+
+# # --- Number of material phases
+# parameterSet.Phases=3
+
+# def indicatorFunction(x):
+#     factor=1
+#     pRadius = 1
+#     # if (np.sqrt(x[0]*x[0] + x[1]*x[1]) < pRadius):
+#     if ((x[0] < 0 and math.sqrt(pow(x[1],2) + pow(x[2],2) ) < pRadius/4.0)  or ( 0 < x[0] and math.sqrt(pow(x[1],2) + pow(x[2],2) ) > pRadius/4.0)):
+#         return 1
+#     else :
+#         return 2      #Phase2
+
+# --- Number of material phases
+parameterSet.Phases=3
+
+
+# Parameters of the model
+# -- (thickness upper layer) / (thickness)
+# param_r = 0.22
+param_r = 0.12
+# -- thickness [meter]
+param_h = 0.0047
+# -- moisture content in the flat state [%]
+param_omega_flat = 17.17547062
+# -- moisture content in the target state [%]
+param_omega_target = 8.959564147
+# -- Drehwinkel
+param_theta = 0.0
+
+# Design Parameter ratio between perforaton volume and volume of perforated layer
+param_beta = 0.3
+# Depth of perforation
+perfDepth = (1.0-param_r) * (3.0/4.0)
+
+#
+#
+#
+# -- increment of the moisture content
+delta_omega=param_omega_target-param_omega_flat
+# moisture content for material law
+omega=param_omega_target
+
+# --- Material properties from https://www.ncbi.nlm.nih.gov/pmc/articles/PMC6191116/#pone.0205607.ref015
+# --- for European beech, moisture content omega = 15%
+# --- L=direction orthogonal to layering and fibres = orthogonal to wood stem cross-section
+# --- T=tangential zu layering
+# --- R=orthogonal zu layering
+# --- in MPa
+# --- Properties are defined by affine function in dependence of moisture content omega via property = b_0+b_1 \omega
+# --- coefficients of affine function are contained in the following array 
+# --- data taken from http://dx.doi.org/10.1016/j.cma.2014.10.031
+
+properties_coefficients=np.array([
+    # [b_0, b_1]    
+    [2565.6,-59.7], # E_R [MPa]
+    [885.4, -23.4], # E_T [MPa]
+    [17136.7,-282.4], # E_L [MPa]
+    [667.8, -15.19], # G_RT [MPa]
+    [1482, -15.26], # G_RL [MPa]
+    [1100, -17.72], # G_TL [MPa]
+    [0.2933, -0.001012], # nu_TR [1]
+    [0.383, -0.008722], # nu_LR [1]
+    [0.3368, -0.009071] # nu_LT [1]
+    ])
+# Compute actual material properties
+E_R = properties_coefficients[0,0]+properties_coefficients[0,1]*omega
+E_T = properties_coefficients[1,0]+properties_coefficients[1,1]*omega
+E_L = properties_coefficients[2,0]+properties_coefficients[2,1]*omega
+G_RT = properties_coefficients[3,0]+properties_coefficients[3,1]*omega
+G_LR = properties_coefficients[4,0]+properties_coefficients[4,1]*omega
+G_LT  = properties_coefficients[5,0]+properties_coefficients[5,1]*omega
+nu_TR  = properties_coefficients[6,0]+properties_coefficients[6,1]*omega
+nu_LR  = properties_coefficients[7,0]+properties_coefficients[7,1]*omega
+nu_LT  = properties_coefficients[8,0]+properties_coefficients[8,1]*omega
+# Compute the remaining Poisson ratios
+nu_TL=nu_LT*E_T/E_L
+nu_RT=nu_TR*E_R/E_T
+nu_RL=nu_LR*E_R/E_L
+#
+# --- differential swelling strain
+# --- relation to swelling strain eps: eps=alpha* delta_omega with delta_omega = change of water content in %
+alpha_L=0.00011 # PLOS paper
+alpha_R=0.00191 # PLOS paper
+alpha_T=0.00462 # PLOS paper
+# Umrechnen
+#alpha_L=(1-1/(1+delta_omega*alpha_L))/delta_omega
+#alpha_R=(1-1/(1+delta_omega*alpha_R))/delta_omega
+#alpha_T=(1-1/(1+delta_omega*alpha_T))/delta_omega
+# --- define geometry
+
+
+
+# --- PHASE 1
+# y_1-direction: L
+# y_2-direction: T
+# x_3-direction: R
+# phase1_type="orthotropic"
+# materialParameters_phase1 = [E_L,E_T,E_R,G_TL,G_RT,G_RL,nu_LT,nu_LR,nu_TR]
+parameterSet.phase1_type="general_anisotropic"
+[E_1,E_2,E_3]=[E_L,E_T,E_R]
+[nu_12,nu_13,nu_23]=[nu_LT,nu_LR,nu_TR]
+[nu_21,nu_31,nu_32]=[nu_TL,nu_RL,nu_RT]
+[G_12,G_31,G_23]=[G_LT,G_LR,G_RT]
+compliance_S=np.array([[1/E_1,         -nu_21/E_2,     -nu_31/E_3,   0.0,         0.0,  0.0],
+                       [-nu_12/E_1,      1/E_2,        -nu_32/E_3,   0.0,         0.0,  0.0],
+                       [-nu_13/E_1,     -nu_23/E_2,         1/E_3,   0.0,         0.0,  0.0],
+                       [0.0,                0.0,              0.0,   1/G_23,         0.0,  0.0],
+                       [0.0,                0.0,              0.0,   0.0,         1/G_31,  0.0],
+                       [0.0,                0.0,              0.0,   0.0,         0.0,  1/G_12]]);
+materialParameters_phase1 = compliance_S
+
+def prestrain_phase1(x):
+    # hB=delta_omega * alpha with delta_omega increment of moisture content and alpha swelling factor.
+    return [[1/param_h*delta_omega*alpha_L, 0, 0], [0,1/param_h*delta_omega*alpha_T,0], [0,0,1/param_h*delta_omega*alpha_R]]
+
+# --- PHASE 2
+# y_1-direction: R
+# y_2-direction: L
+# x_3-direction: T
+parameterSet.phase2_type="general_anisotropic"
+[E_1,E_2,E_3]=[E_R,E_L,E_T]
+[nu_12,nu_13,nu_23]=[nu_RL,nu_RT,nu_LT]
+[nu_21,nu_31,nu_32]=[nu_LR,nu_TR,nu_TL]
+[G_12,G_31,G_23]=[G_LR,G_RT,G_LT]
+compliance_S=np.array([[1/E_1,         -nu_21/E_2,     -nu_31/E_3,   0.0,         0.0,  0.0],
+                       [-nu_12/E_1,      1/E_2,        -nu_32/E_3,   0.0,         0.0,  0.0],
+                       [-nu_13/E_1,     -nu_23/E_2,         1/E_3,   0.0,         0.0,  0.0],
+                       [0.0,                0.0,              0.0,   1/G_23,         0.0,  0.0],
+                       [0.0,                0.0,              0.0,   0.0,         1/G_31,  0.0],
+                       [0.0,                0.0,              0.0,   0.0,         0.0,  1/G_12]]);
+materialParameters_phase2 = compliance_S
+def prestrain_phase2(x):
+    return [[1/param_h*delta_omega*alpha_R, 0, 0], [0,1/param_h*delta_omega*alpha_L,0], [0,0,1/param_h*delta_omega*alpha_T]]
+
+#Rotation um 2. Achse (= L) 
+parameterSet.phase2_axis = 1
+# phase2_angle = param_theta
+# -- Drehwinkel
+parameterSet.phase2_angle = param_theta
+
+
+# --- PHASE 3 
+parameterSet.phase3_type="isotropic"
+epsilon = 1e-8
+materialParameters_phase3 = [epsilon, epsilon]
+
+def prestrain_phase3(x):
+    return [[0, 0, 0], [0,0,0], [0,0,0]]
+
+# # --- PHASE 3 = Phase 1 gedreht
+# # y_1-direction: L
+# # y_2-direction: R
+# # x_3-direction: T
+# parameterSet.phase3_type="general_anisotropic"
+# # Drehung um theta um Achse 2 = x_3-Achse
+# N=elast.rotation_matrix_compliance(2,param_theta)
+# materialParameters_phase3 = np.dot(np.dot(N,materialParameters_phase1),N.T)
+# materialParameters_phase3 = 0.5*(materialParameters_phase3.T+materialParameters_phase3)
+# # rotation of strain
+# def prestrain_phase3(x):
+#     return elast.voigt_to_strain(np.dot(elast.rotation_matrix_compliance(2,param_theta),np.dot(elast.strain_to_voigt(np.array(prestrain_phase1(x))),N.T))).tolist()
+
+
+
+# --- Choose scale ratio gamma:
+parameterSet.gamma=1.0
+
+
+
+
+#############################################
+#  Grid parameters
+#############################################
+## numLevels : Number of Levels on which solution is computed. starting with a 2x2x2 cube mesh.
+## {start,finish} computes on all grid from 2^(start) to 2^finish refinement
+#----------------------------------------------------
+parameterSet.numLevels= '3 3'      # computes all levels from first to second entry
+# parameterSet.numLevels= '4 4' 
+
+#############################################
+#  Assembly options
+#############################################
+parameterSet.set_IntegralZero = 1            #(default = false)
+parameterSet.set_oneBasisFunction_Zero = 1   #(default = false)
+#parameterSet.arbitraryLocalIndex = 7            #(default = 0)
+#parameterSet.arbitraryElementNumber = 3         #(default = 0)
+
+#############################################
+#  Solver Options, Type: #1: CG - SOLVER , #2: GMRES - SOLVER, #3: QR - SOLVER (default), #4: UMFPACK - SOLVER
+#############################################
+parameterSet.Solvertype = 3        # recommended to use iterative solver (e.g GMRES) for finer grid-levels
+parameterSet.Solver_verbosity = 0  #(default = 2)  degree of information for solver output
+
+
+#############################################
+#  Write/Output options      #(default=false)
+#############################################
+# --- (Optional output) write Material / prestrain / Corrector functions to .vtk-Files:
+parameterSet.write_materialFunctions = 1   # VTK indicator function for material/prestrain definition
+#parameterSet.write_prestrainFunctions = 1  # VTK norm of B (currently not implemented)
+
+# --- (Additional debug output)
+parameterSet.print_debug = 0  #(default=false)
+
+# --- Write Correctos to VTK-File:  
+parameterSet.write_VTK = 1
+
+# The grid can be refined several times for a higher resolution in the VTK-file.
+parameterSet.subsamplingRefinement = 2
+
+# --- (Optional output) L2Error, integral mean: 
+#parameterSet.write_L2Error = 1
+#parameterSet.write_IntegralMean = 1      
+
+# --- check orthogonality (75) from paper: 
+parameterSet.write_checkOrthogonality = 0
+
+# --- Write corrector-coefficients to log-File:
+#parameterSet.write_corrector_phi1 = 1
+#parameterSet.write_corrector_phi2 = 1
+#parameterSet.write_corrector_phi3 = 1
+
+# --- Print Condition number of matrix (can be expensive):
+#parameterSet.print_conditionNumber= 1  #(default=false)
+
+# --- write effective quantities to Matlab-folder for symbolic minimization:
+parameterSet.write_toMATLAB = 1  # writes effective quantities to .txt-files QMatrix.txt and BMatrix.txt

experiment/perforated-bilayer_square/perforated_wood_upper.py

+import math
+#from python_matrix_operations import *
+import ctypes
+import os
+import sys
+import numpy as np
+# import elasticity_toolbox as elast
+
+class ParameterSet(dict):
+    def __init__(self, *args, **kwargs):
+        super(ParameterSet, self).__init__(*args, **kwargs)
+        self.__dict__ = self
+
+parameterSet = ParameterSet()
+#---------------------------------------------------------------
+#############################################
+#  Paths
+#############################################
+# Path for results and logfile
+parameterSet.outputPath='/home/klaus/Desktop/Dune_release/dune-microstructure/experiment/perforated-bilayer_square/results'
+parameterSet.baseName= 'perforated_wood_upper'   #(needed for Output-Filename)
+
+# Path for material description
+# parameterSet.geometryFunctionPath =experiment/wood-bilayer/
+
+#---------------------------------------------------------------
+# Wooden bilayer, https://www.ncbi.nlm.nih.gov/pmc/articles/PMC6191116/#pone.0205607.ref015
+#--- define indicator function
+# x[0] : y1-component -1/2 to 1/2
+# x[1] : y2-component -1/2 to 1/2
+# x[2] : x3-component range -1/2 to 1/2
+#--- define indicator function
+def indicatorFunction(x):
+    pRadius = math.sqrt((param_beta*param_r)/(4.0*perfDepth))
+    if (x[2]>=(0.5-param_r)):
+        if( (max(abs(x[0]), abs(x[1])) < pRadius) and (x[2] >= (0.5-perfDepth)) ): #inside perforation    
+            return 3  #Phase3
+        else:  
+            return 1  #Phase1
+    else :
+        return 2      #Phase2
+    
+# def indicatorFunction(x):
+#     factor=1
+#     pRadius = 0.25
+#     if (x[2]>=(0.5-param_r) and np.sqrt(x[0]**2 + x[1]**2) < pRadius):
+#         return 3
+#     elif((x[2]>=(0.5-param_r))):  
+#         return 1  #Phase1
+#     else :
+#         return 2      #Phase2
+
+# # --- Number of material phases
+# parameterSet.Phases=3
+
+# def indicatorFunction(x):
+#     factor=1
+#     pRadius = 1
+#     # if (np.sqrt(x[0]*x[0] + x[1]*x[1]) < pRadius):
+#     if ((x[0] < 0 and math.sqrt(pow(x[1],2) + pow(x[2],2) ) < pRadius/4.0)  or ( 0 < x[0] and math.sqrt(pow(x[1],2) + pow(x[2],2) ) > pRadius/4.0)):
+#         return 1
+#     else :
+#         return 2      #Phase2
+
+# --- Number of material phases
+parameterSet.Phases=3
+
+
+# Parameters of the model
+# -- (thickness upper layer) / (thickness)
+# param_r = 0.22
+param_r = 0.12
+# -- thickness [meter]
+param_h = 0.0047
+# -- moisture content in the flat state [%]
+param_omega_flat = 17.17547062
+# -- moisture content in the target state [%]
+param_omega_target = 8.959564147
+# -- Drehwinkel
+param_theta = 0.0
+
+# Design Parameter ratio between perforaton (cylindrical) volume and volume of perforated layer
+param_beta = 0.3
+# Depth of perforation
+perfDepth = param_r * (3.0/4.0)
+
+#
+#
+#
+# -- increment of the moisture content
+delta_omega=param_omega_target-param_omega_flat
+# moisture content for material law
+omega=param_omega_target
+
+# --- Material properties from https://www.ncbi.nlm.nih.gov/pmc/articles/PMC6191116/#pone.0205607.ref015
+# --- for European beech, moisture content omega = 15%
+# --- L=direction orthogonal to layering and fibres = orthogonal to wood stem cross-section
+# --- T=tangential zu layering
+# --- R=orthogonal zu layering
+# --- in MPa
+# --- Properties are defined by affine function in dependence of moisture content omega via property = b_0+b_1 \omega
+# --- coefficients of affine function are contained in the following array 
+# --- data taken from http://dx.doi.org/10.1016/j.cma.2014.10.031
+
+properties_coefficients=np.array([
+    # [b_0, b_1]    
+    [2565.6,-59.7], # E_R [MPa]
+    [885.4, -23.4], # E_T [MPa]
+    [17136.7,-282.4], # E_L [MPa]
+    [667.8, -15.19], # G_RT [MPa]
+    [1482, -15.26], # G_RL [MPa]
+    [1100, -17.72], # G_TL [MPa]
+    [0.2933, -0.001012], # nu_TR [1]
+    [0.383, -0.008722], # nu_LR [1]
+    [0.3368, -0.009071] # nu_LT [1]
+    ])
+# Compute actual material properties
+E_R = properties_coefficients[0,0]+properties_coefficients[0,1]*omega
+E_T = properties_coefficients[1,0]+properties_coefficients[1,1]*omega
+E_L = properties_coefficients[2,0]+properties_coefficients[2,1]*omega
+G_RT = properties_coefficients[3,0]+properties_coefficients[3,1]*omega
+G_LR = properties_coefficients[4,0]+properties_coefficients[4,1]*omega
+G_LT  = properties_coefficients[5,0]+properties_coefficients[5,1]*omega
+nu_TR  = properties_coefficients[6,0]+properties_coefficients[6,1]*omega
+nu_LR  = properties_coefficients[7,0]+properties_coefficients[7,1]*omega
+nu_LT  = properties_coefficients[8,0]+properties_coefficients[8,1]*omega
+# Compute the remaining Poisson ratios
+nu_TL=nu_LT*E_T/E_L
+nu_RT=nu_TR*E_R/E_T
+nu_RL=nu_LR*E_R/E_L
+#
+# --- differential swelling strain
+# --- relation to swelling strain eps: eps=alpha* delta_omega with delta_omega = change of water content in %
+alpha_L=0.00011 # PLOS paper
+alpha_R=0.00191 # PLOS paper
+alpha_T=0.00462 # PLOS paper
+# Umrechnen
+#alpha_L=(1-1/(1+delta_omega*alpha_L))/delta_omega
+#alpha_R=(1-1/(1+delta_omega*alpha_R))/delta_omega
+#alpha_T=(1-1/(1+delta_omega*alpha_T))/delta_omega
+# --- define geometry
+
+
+
+# --- PHASE 1
+# y_1-direction: L
+# y_2-direction: T
+# x_3-direction: R
+# phase1_type="orthotropic"
+# materialParameters_phase1 = [E_L,E_T,E_R,G_TL,G_RT,G_RL,nu_LT,nu_LR,nu_TR]
+parameterSet.phase1_type="general_anisotropic"
+[E_1,E_2,E_3]=[E_L,E_T,E_R]
+[nu_12,nu_13,nu_23]=[nu_LT,nu_LR,nu_TR]
+[nu_21,nu_31,nu_32]=[nu_TL,nu_RL,nu_RT]
+[G_12,G_31,G_23]=[G_LT,G_LR,G_RT]
+compliance_S=np.array([[1/E_1,         -nu_21/E_2,     -nu_31/E_3,   0.0,         0.0,  0.0],
+                       [-nu_12/E_1,      1/E_2,        -nu_32/E_3,   0.0,         0.0,  0.0],
+                       [-nu_13/E_1,     -nu_23/E_2,         1/E_3,   0.0,         0.0,  0.0],
+                       [0.0,                0.0,              0.0,   1/G_23,         0.0,  0.0],
+                       [0.0,                0.0,              0.0,   0.0,         1/G_31,  0.0],
+                       [0.0,                0.0,              0.0,   0.0,         0.0,  1/G_12]]);
+materialParameters_phase1 = compliance_S
+
+def prestrain_phase1(x):
+    # hB=delta_omega * alpha with delta_omega increment of moisture content and alpha swelling factor.
+    return [[1/param_h*delta_omega*alpha_L, 0, 0], [0,1/param_h*delta_omega*alpha_T,0], [0,0,1/param_h*delta_omega*alpha_R]]
+
+# --- PHASE 2
+# y_1-direction: R
+# y_2-direction: L
+# x_3-direction: T
+parameterSet.phase2_type="general_anisotropic"
+[E_1,E_2,E_3]=[E_R,E_L,E_T]
+[nu_12,nu_13,nu_23]=[nu_RL,nu_RT,nu_LT]
+[nu_21,nu_31,nu_32]=[nu_LR,nu_TR,nu_TL]
+[G_12,G_31,G_23]=[G_LR,G_RT,G_LT]
+compliance_S=np.array([[1/E_1,         -nu_21/E_2,     -nu_31/E_3,   0.0,         0.0,  0.0],
+                       [-nu_12/E_1,      1/E_2,        -nu_32/E_3,   0.0,         0.0,  0.0],
+                       [-nu_13/E_1,     -nu_23/E_2,         1/E_3,   0.0,         0.0,  0.0],
+                       [0.0,                0.0,              0.0,   1/G_23,         0.0,  0.0],
+                       [0.0,                0.0,              0.0,   0.0,         1/G_31,  0.0],
+                       [0.0,                0.0,              0.0,   0.0,         0.0,  1/G_12]]);
+materialParameters_phase2 = compliance_S
+def prestrain_phase2(x):
+    return [[1/param_h*delta_omega*alpha_R, 0, 0], [0,1/param_h*delta_omega*alpha_L,0], [0,0,1/param_h*delta_omega*alpha_T]]
+
+#Rotation um 2. Achse (= L) 
+parameterSet.phase2_axis = 1
+# phase2_angle = param_theta
+# -- Drehwinkel
+parameterSet.phase2_angle = param_theta
+
+
+# --- PHASE 3 
+parameterSet.phase3_type="isotropic"
+epsilon = 1e-8
+materialParameters_phase3 = [epsilon, epsilon]
+
+def prestrain_phase3(x):
+    return [[0, 0, 0], [0,0,0], [0,0,0]]
+
+# # --- PHASE 3 = Phase 1 gedreht
+# # y_1-direction: L
+# # y_2-direction: R
+# # x_3-direction: T
+# parameterSet.phase3_type="general_anisotropic"
+# # Drehung um theta um Achse 2 = x_3-Achse
+# N=elast.rotation_matrix_compliance(2,param_theta)
+# materialParameters_phase3 = np.dot(np.dot(N,materialParameters_phase1),N.T)
+# materialParameters_phase3 = 0.5*(materialParameters_phase3.T+materialParameters_phase3)
+# # rotation of strain
+# def prestrain_phase3(x):
+#     return elast.voigt_to_strain(np.dot(elast.rotation_matrix_compliance(2,param_theta),np.dot(elast.strain_to_voigt(np.array(prestrain_phase1(x))),N.T))).tolist()
+
+
+
+# --- Choose scale ratio gamma:
+parameterSet.gamma=1.0
+
+
+
+
+#############################################
+#  Grid parameters
+#############################################
+## numLevels : Number of Levels on which solution is computed. starting with a 2x2x2 cube mesh.
+## {start,finish} computes on all grid from 2^(start) to 2^finish refinement
+#----------------------------------------------------
+parameterSet.numLevels= '3 3'      # computes all levels from first to second entry
+# parameterSet.numLevels= '4 4' 
+
+#############################################
+#  Assembly options
+#############################################
+parameterSet.set_IntegralZero = 1            #(default = false)
+parameterSet.set_oneBasisFunction_Zero = 1   #(default = false)
+#parameterSet.arbitraryLocalIndex = 7            #(default = 0)
+#parameterSet.arbitraryElementNumber = 3         #(default = 0)
+
+#############################################
+#  Solver Options, Type: #1: CG - SOLVER , #2: GMRES - SOLVER, #3: QR - SOLVER (default), #4: UMFPACK - SOLVER
+#############################################
+parameterSet.Solvertype = 3        # recommended to use iterative solver (e.g GMRES) for finer grid-levels
+parameterSet.Solver_verbosity = 0  #(default = 2)  degree of information for solver output
+
+
+#############################################
+#  Write/Output options      #(default=false)
+#############################################
+# --- (Optional output) write Material / prestrain / Corrector functions to .vtk-Files:
+parameterSet.write_materialFunctions = 1   # VTK indicator function for material/prestrain definition
+#parameterSet.write_prestrainFunctions = 1  # VTK norm of B (currently not implemented)
+
+# --- (Additional debug output)
+parameterSet.print_debug = 0  #(default=false)
+
+# --- Write Correctos to VTK-File:  
+parameterSet.write_VTK = 1
+
+# The grid can be refined several times for a higher resolution in the VTK-file.
+parameterSet.subsamplingRefinement = 0
+
+# --- (Optional output) L2Error, integral mean: 
+#parameterSet.write_L2Error = 1
+#parameterSet.write_IntegralMean = 1      
+
+# --- check orthogonality (75) from paper: 
+parameterSet.write_checkOrthogonality = 0
+
+# --- Write corrector-coefficients to log-File:
+#parameterSet.write_corrector_phi1 = 1
+#parameterSet.write_corrector_phi2 = 1
+#parameterSet.write_corrector_phi3 = 1
+
+# --- Print Condition number of matrix (can be expensive):
+#parameterSet.print_conditionNumber= 1  #(default=false)
+
+# --- write effective quantities to Matlab-folder for symbolic minimization:
+parameterSet.write_toMATLAB = 1  # writes effective quantities to .txt-files QMatrix.txt and BMatrix.txt