diff --git a/experiment/wood-bilayer/PolarPlotLocalEnergy.py b/experiment/wood-bilayer/PolarPlotLocalEnergy.py
new file mode 100644
index 0000000000000000000000000000000000000000..fae4dbb6d3c091afa7c5e21c22cbd6a068da1a52
--- /dev/null
+++ b/experiment/wood-bilayer/PolarPlotLocalEnergy.py
@@ -0,0 +1,93 @@
+#!/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()
+
+
diff --git a/experiment/wood-bilayer/cellsolver.parset.wood b/experiment/wood-bilayer/cellsolver.parset.wood
new file mode 100644
index 0000000000000000000000000000000000000000..aee5f271338618094934f67f7b2ca61840d955a5
--- /dev/null
+++ b/experiment/wood-bilayer/cellsolver.parset.wood
@@ -0,0 +1,96 @@
+# --- 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
diff --git a/experiment/wood-bilayer/elasticity_toolbox.py b/experiment/wood-bilayer/elasticity_toolbox.py
new file mode 100644
index 0000000000000000000000000000000000000000..8e61952612c0714a5b430a41660775fc0e2c23b5
--- /dev/null
+++ b/experiment/wood-bilayer/elasticity_toolbox.py
@@ -0,0 +1,123 @@
+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
diff --git a/experiment/wood-bilayer/results/Experimental_data_beech-beech-bilayers.ods b/experiment/wood-bilayer/results/Experimental_data_beech-beech-bilayers.ods
new file mode 100644
index 0000000000000000000000000000000000000000..1ad6db2b60d5b7bd9c2c322f26661a56dbfea4fe
Binary files /dev/null and b/experiment/wood-bilayer/results/Experimental_data_beech-beech-bilayers.ods differ
diff --git a/experiment/wood-bilayer/results/auswertung.py b/experiment/wood-bilayer/results/auswertung.py
new file mode 100644
index 0000000000000000000000000000000000000000..12e9fc7a2c88da4429613778cf6044c325abe1eb
--- /dev/null
+++ b/experiment/wood-bilayer/results/auswertung.py
@@ -0,0 +1,12 @@
+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()
diff --git a/experiment/wood-bilayer/wood_european_beech.py b/experiment/wood-bilayer/wood_european_beech.py
new file mode 100644
index 0000000000000000000000000000000000000000..bcaf914a2ba135ed887444828e89fcd9026f2aba
--- /dev/null
+++ b/experiment/wood-bilayer/wood_european_beech.py
@@ -0,0 +1,151 @@
+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()
diff --git a/experiment/wood-bilayer/wood_test.py b/experiment/wood-bilayer/wood_test.py
new file mode 100644
index 0000000000000000000000000000000000000000..56a48fd0d4126f0f8c38739f820ef1121bb6e2a2
--- /dev/null
+++ b/experiment/wood-bilayer/wood_test.py
@@ -0,0 +1,160 @@
+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()
+ #