experiment folder, wood bilayer setup from Reuggeberg Paper

experiment/wood-bilayer/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
+
+
+
+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=7
+kappa=np.zeros(number)
+for n in range(0,number):
+    #   Read from Date
+    print(str(n))
+    DataPath = './experiment/wood-bilayer/results/'+str(n)
+    QFilePath = DataPath + '/QMatrix.txt'
+    BFilePath = DataPath + '/BMatrix.txt'
+    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)
+            else:
+                E[i,j]=energy(-r[i,j],theta[i,j],Q,B)
+            
+    # 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)
+    plt.colorbar(pcm, extend='max', ticks=colorbarticks, pad=0.1)
+    plt.show()
+
+f = open("./experiment/wood-bilayer/results/kappa_simulation.txt", "w")
+f.write(str(kappa))         
+f.close()   
+
+

experiment/wood-bilayer/cellsolver.parset.wood

+# --- Parameter File as Input for 'Cell-Problem'
+# NOTE: define variables without whitespaces in between! i.e. : gamma=1.0 instead of gamma = 1.0
+# since otherwise these cant be read from other Files!
+# --------------------------------------------------------
+
+# Path for results and logfile
+outputPath=./experiment/wood-bilayer/results/6
+
+# Path for material description
+geometryFunctionPath =experiment/wood-bilayer/
+
+
+# --- DEBUG (Output) Option:
+#print_debug = true  #(default=false)
+
+
+
+
+#############################################
+#  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
+#----------------------------------------------------
+
+numLevels=4 4
+#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 = 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
+
+
+#############################################
+#  Material / Prestrain parameters and ratios
+#############################################
+
+# --- Choose material definition:
+materialFunction = wood_european_beech
+
+
+
+# --- Choose scale ratio gamma:
+gamma=1.0
+
+
+#############################################
+#  Assembly options
+#############################################
+#set_IntegralZero = true            #(default = false)
+#set_oneBasisFunction_Zero = true   #(default = false)
+
+#arbitraryLocalIndex = 7            #(default = 0)
+#arbitraryElementNumber = 3         #(default = 0)
+#############################################
+
+
+#############################################
+#  Solver Type: #1: CG - SOLVER , #2: GMRES - SOLVER, #3: QR - SOLVER (default), #4: UMFPACK - SOLVER
+#############################################
+Solvertype = 2        # recommended to use iterative solver (e.g GMRES) for finer grid-levels
+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:
+write_materialFunctions = true    # VTK indicator function for material/prestrain definition
+#write_prestrainFunctions = true  # VTK norm of B (currently not implemented)
+
+# --- Write Correctos to VTK-File:  
+write_VTK = true
+
+# --- (Optional output) L2Error, integral mean: 
+#write_L2Error = true                   
+#write_IntegralMean = true              
+
+# --- check orthogonality (75) from paper: 
+write_checkOrthogonality = true         
+
+# --- Write corrector-coefficients to log-File:
+#write_corrector_phi1 = true
+#write_corrector_phi2 = true
+#write_corrector_phi3 = true
+
+
+# --- Print Condition number of matrix (can be expensive):
+#print_conditionNumber= true  #(default=false)
+
+# --- write effective quantities to Matlab-folder for symbolic minimization:
+write_toMATLAB = true  # writes effective quantities to .txt-files QMatrix.txt and BMatrix.txt

experiment/wood-bilayer/elasticity_toolbox.py

+import math
+import numpy as np
+
+ 
+def strain_to_voigt(strain_matrix):
+    # Ensure the input matrix is a 3x3 strain matrix
+    if strain_matrix.shape != (3, 3):
+        raise ValueError("Input matrix should be a 3x3 strain matrix.")
+
+    # Extract the components from the 3x3 strain matrix
+    ε_xx = strain_matrix[0, 0]
+    ε_yy = strain_matrix[1, 1]
+    ε_zz = strain_matrix[2, 2]
+    γ_yz = .5*(strain_matrix[1, 2]+strain_matrix[2,1])
+    γ_xz = .5*(strain_matrix[0, 2]+strain_matrix[0,2])
+    γ_xy = .5*(strain_matrix[0, 1]+strain_matrix[0,1])
+
+    # Create the Voigt notation vector
+    voigt_notation = np.array([ε_xx, ε_yy, ε_zz, γ_yz, γ_xz, γ_xy])
+
+    return voigt_notation
+
+def voigt_to_strain(voigt_notation):
+    # Ensure the input vector has 6 elements
+    if len(voigt_notation) != 6:
+        raise ValueError("Input vector should have 6 elements in Voigt notation.")
+
+    # Extract the components from the Voigt notation vector
+    ε_xx = voigt_notation[0]
+    ε_yy = voigt_notation[1]
+    ε_zz = voigt_notation[2]
+    γ_yz = voigt_notation[3]
+    γ_xz = voigt_notation[4]
+    γ_xy = voigt_notation[5]
+
+    # Create the 3x3 strain matrix
+    strain_matrix = np.array([[ε_xx, γ_xy, γ_xz],
+                              [γ_xy, ε_yy, γ_yz],
+                              [γ_xz, γ_yz, ε_zz]])
+
+    return strain_matrix
+
+
+def rotation_matrix(ax, angle):
+    cos_theta = np.cos(angle)
+    sin_theta = np.sin(angle)
+    if ax==0:
+        Q=np.array([[0, 0, 1],
+                    [0,1,0],
+                    [-1,0,0]
+                    ])
+    elif ax==1:
+        Q=np.array([[1, 0, 0],
+                    [0,0,1],
+                    [0,-1,0]
+                    ])
+    else:
+        Q=np.array([[1, 0, 0],
+                    [0,1,0],
+                    [0,0,1]
+                    ])
+        
+    R = np.array([[cos_theta, -sin_theta, 0],
+                  [sin_theta, cos_theta, 0],
+                  [0, 0, 1]])
+    return np.dot(np.dot(Q.T, R),Q)
+
+def rotation_matrix_compliance(ax,theta):
+    R=rotation_matrix(ax,theta)
+    Q_xx=R[0,0]
+    Q_xy=R[0,1]
+    Q_xz=R[0,2]
+    Q_yx=R[1,0]
+    Q_yy=R[1,1]
+    Q_yz=R[1,2]
+    Q_zx=R[2,0]
+    Q_zy=R[2,1]
+    Q_zz=R[2,2]
+    return np.array([
+        [Q_xx**2, Q_xy**2, Q_xz**2, Q_xy*Q_xz, Q_xx*Q_xz, Q_xx*Q_xy],
+        [Q_yx**2, Q_yy**2, Q_yz**2, Q_yy*Q_yz, Q_yx*Q_yz, Q_yx*Q_yy],
+        [Q_zx**2, Q_zy**2, Q_zz**2, Q_zy*Q_zz, Q_zx*Q_zz, Q_zx*Q_zy],
+        [2*Q_yx*Q_zx, 2*Q_yy*Q_zy, 2*Q_yz*Q_zz, Q_yy*Q_zz + Q_yz*Q_zy, Q_yx*Q_zz + Q_yz*Q_zx, Q_yx*Q_zy + Q_yy*Q_zx],
+        [2*Q_xx*Q_zx, 2*Q_xy*Q_zy, 2*Q_xz*Q_zz, Q_xy*Q_zz + Q_xz*Q_zy, Q_xx*Q_zz + Q_xz*Q_zx, Q_xx*Q_zy + Q_xy*Q_zx],
+        [2*Q_xx*Q_yx, 2*Q_xy*Q_yy, 2*Q_xz*Q_yz, Q_xy*Q_yz + Q_xz*Q_yy, Q_xx*Q_yz + Q_xz*Q_yx, Q_xx*Q_yy + Q_xy*Q_yx]
+    ])
+
+def rotate_strain(eps, ax, theta):
+    B=voigt_to_strain(np.matmul(rotation_matrix_epsilon(theta,ax),strain_to_voigt(eps)))
+
+import numpy as np
+
+def voigt_to_tensor(voigt_matrix):
+    tensor = np.zeros((6, 6))
+
+    tensor[0, 0] = voigt_matrix[0]
+    tensor[0, 1] = tensor[1, 0] = voigt_matrix[1]
+    tensor[0, 2] = tensor[2, 0] = voigt_matrix[2]
+    tensor[0, 3] = tensor[3, 0] = voigt_matrix[3]
+    tensor[0, 4] = tensor[4, 0] = voigt_matrix[4]
+    tensor[0, 5] = tensor[5, 0] = voigt_matrix[5]
+
+    tensor[1, 1] = voigt_matrix[6]
+    tensor[1, 2] = tensor[2, 1] = voigt_matrix[7]
+    tensor[1, 3] = tensor[3, 1] = voigt_matrix[8]
+    tensor[1, 4] = tensor[4, 1] = voigt_matrix[9]
+    tensor[1, 5] = tensor[5, 1] = voigt_matrix[10]
+
+    tensor[2, 2] = voigt_matrix[11]
+    tensor[2, 3] = tensor[3, 2] = voigt_matrix[12]
+    tensor[2, 4] = tensor[4, 2] = voigt_matrix[13]
+    tensor[2, 5] = tensor[5, 2] = voigt_matrix[14]
+
+    tensor[3, 3] = voigt_matrix[15]
+    tensor[3, 4] = tensor[4, 3] = voigt_matrix[16]
+    tensor[3, 5] = tensor[5, 3] = voigt_matrix[17]
+
+    tensor[4, 4] = voigt_matrix[18]
+    tensor[4, 5] = tensor[5, 4] = voigt_matrix[19]
+
+    tensor[5, 5] = voigt_matrix[20]
+
+    return tensor

experiment/wood-bilayer/results/auswertung.py

+import matplotlib.pyplot as plt
+omega=[
+    14.72680026, 13.64338887, 12.41305478, 11.66482931, 11.09781471, 9.435795985, 8.959564147]
+kappa_sim=[1.20240481, 1.73346693, 2.3246493,  2.68537074, 2.95591182, 3.74749499,
+ 3.97795591]
+kappa_exp=[
+     1.058078122, 1.544624544, 2.317033799, 2.686043143, 2.967694189, 3.913528418, 4.262750825]
+plt.scatter(omega,kappa_sim, label="simulation")
+plt.scatter(omega, kappa_exp, marker='x', label="experiment")
+plt.xlabel="target humidity"
+plt.ylabel="curvature"
+plt.show()

experiment/wood-bilayer/wood_european_beech.py

+import math
+#from python_matrix_operations import *
+import ctypes
+import os
+import sys
+import numpy as np
+import elasticity_toolbox as elast
+
+#---------------------------------------------------------------
+# 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
+
+# Parameters of the model
+# -- ratio between upper layer / thickness, thickness [m]
+param_r = 0.22
+param_h = 0.0053
+# -- moisture content in the flat state and target state [%]
+param_omega_flat = 17.17547062
+param_omega_target = 8.959564147
+param_theta = 0
+#
+#
+#
+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 properties
+#E_R=1500
+E_R = properties_coefficients[0,0]+properties_coefficients[0,1]*omega
+# E_T=450
+E_T = properties_coefficients[1,0]+properties_coefficients[1,1]*omega
+#E_L=12000
+E_L = properties_coefficients[2,0]+properties_coefficients[2,1]*omega
+# G_RT=400
+G_RT = properties_coefficients[3,0]+properties_coefficients[3,1]*omega
+# G_RL=1200
+G_LR = properties_coefficients[4,0]+properties_coefficients[4,1]*omega
+# G_TL=800
+G_LT  = properties_coefficients[5,0]+properties_coefficients[5,1]*omega
+# nu_TR=0.28
+nu_TR  = properties_coefficients[6,0]+properties_coefficients[6,1]*omega
+# nu_LR=0.2125
+nu_LR  = properties_coefficients[7,0]+properties_coefficients[7,1]*omega
+# nu_LT=0.175
+nu_LT  = properties_coefficients[8,0]+properties_coefficients[8,1]*omega
+#
+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
+
+def indicatorFunction(x):
+    factor=1
+    if (x[2]>=(0.5-param_r)):
+        return 1  #Phase1
+    else :
+        return 2   #Phase2
+
+# --- Number of material phases
+Phases=3
+
+# --- 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]
+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):
+    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
+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
+#phase2_type="orthotropic"
+#materialParameters_phase2 = [E_R, E_L, E_T,G_RL, G_TL, G_RT,nu_LR, nu_TR,nu_LT]
+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]]
+
+# --- PHASE 3
+# y_1-direction: L
+# y_2-direction: R
+# x_3-direction: T
+phase3_type="general_anisotropic"
+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()

experiment/wood-bilayer/wood_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()
+
+# 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 -----
+Path = "./experiment/wood-bilayer"
+# parset = ' ./experiment/wood-bilayer/cellsolver.parset.wood'
+ParsetFile = Path + '/cellsolver.parset.wood'
+executable = 'build-cmake/src/Cell-Problem'
+write_LOG = True   # writes Cell-Problem output-LOG in "Cell-Problem_output.log"
+
+# ---------------------------------
+# Setup Experiment
+# ---------------------------------
+outputPath = Path + '/results/'
+
+# ----- Define Input parameters  --------------------
+class ParameterSet:
+    pass
+
+ParameterSet.materialFunction  = "wood_european_beech"
+ParameterSet.gamma=1.0
+ParameterSet.numLevels=4
+
+# ----- Define Parameters for Material Function  --------------------
+# [r, h, omega_flat, omega_target, theta, experimental_kappa]
+materialFunctionParameter=[
+   [0.22, 0.0053, 17.17547062, 14.72680026, 0, 1.058078122],
+   [0.22, 0.0053, 17.17547062, 13.64338887, 0, 1.544624544],
+   [0.22, 0.0053, 17.17547062, 12.41305478, 0, 2.317033799],
+   [0.22, 0.0053, 17.17547062, 11.66482931, 0, 2.686043143],
+   [0.22, 0.0053, 17.17547062, 11.09781471, 0, 2.967694189],
+   [0.22, 0.0053, 17.17547062, 9.435795985, 0, 3.913528418],
+   [0.22, 0.0053,  17.17547062, 8.959564147, 0, 4.262750825]
+]
+# materialFunctionParameter=[
+#      [0.22, 0.0053,  17.17547062, 8.959564147, 0, 4.262750825],
+#    [0.22, 0.0053,  17.17547062, 8.959564147, np.pi/2, 4.262750825]
+# ]
+
+# ------ Loops through Parameters for Material Function -----------
+for i in range(0,np.shape(materialFunctionParameter)[0]):
+    print("------------------")
+    print("New Loop")
+    print("------------------")
+   # Check output directory
+    path = outputPath + str(i)
+    isExist = os.path.exists(path)
+    if not isExist:
+        # Create a new directory because it does not exist
+        os.makedirs(path)
+        print("The new directory " + path + " is created!")
+    SetParameterMaterialFunction(ParameterSet.materialFunction, "param_r",materialFunctionParameter[i][0])
+    SetParameterMaterialFunction(ParameterSet.materialFunction, "param_h",materialFunctionParameter[i][1])
+    SetParameterMaterialFunction(ParameterSet.materialFunction, "param_omega_flat",materialFunctionParameter[i][2])
+    SetParameterMaterialFunction(ParameterSet.materialFunction, "param_omega_target",materialFunctionParameter[i][3])
+    SetParameterMaterialFunction(ParameterSet.materialFunction, "param_theta",materialFunctionParameter[i][4])    
+    SetParametersCellProblem(ParameterSet, ParsetFile, path)
+    #Run Cell-Problem
+    thread = threading.Thread(target=run_CellProblem(executable, " ./"+ParsetFile, write_LOG))
+    thread.start()
+    # ---------------------------------------------------
+    # wait here for the result to be available before continuing
+    thread.join()
+    f = open(path+"/parameter.txt", "w")
+    f.write("r = "+str(materialFunctionParameter[i][0])+"\n")
+    f.write("h = "+str(materialFunctionParameter[i][1])+"\n")
+    f.write("omega_flat = "+str(materialFunctionParameter[i][2])+"\n")        
+    f.write("omega_target = "+str(materialFunctionParameter[i][3])+"\n")         
+    f.close()   
+    #