From 6f6cde2728c528e0f552bfb3aa7341e344dbe2c3 Mon Sep 17 00:00:00 2001
From: Stefan Neukamm <stefan.neukamm@tu-dresden.de>
Date: Wed, 12 Jul 2023 07:31:15 +0200
Subject: [PATCH] experiment folder, wood bilayer setup from Reuggeberg Paper
---
.../wood-bilayer/PolarPlotLocalEnergy.py | 93 ++++++++++
.../wood-bilayer/cellsolver.parset.wood | 96 +++++++++++
experiment/wood-bilayer/elasticity_toolbox.py | 123 ++++++++++++++
...Experimental_data_beech-beech-bilayers.ods | Bin 0 -> 18897 bytes
experiment/wood-bilayer/results/auswertung.py | 12 ++
.../wood-bilayer/wood_european_beech.py | 151 +++++++++++++++++
experiment/wood-bilayer/wood_test.py | 160 ++++++++++++++++++
7 files changed, 635 insertions(+)
create mode 100644 experiment/wood-bilayer/PolarPlotLocalEnergy.py
create mode 100644 experiment/wood-bilayer/cellsolver.parset.wood
create mode 100644 experiment/wood-bilayer/elasticity_toolbox.py
create mode 100644 experiment/wood-bilayer/results/Experimental_data_beech-beech-bilayers.ods
create mode 100644 experiment/wood-bilayer/results/auswertung.py
create mode 100644 experiment/wood-bilayer/wood_european_beech.py
create mode 100644 experiment/wood-bilayer/wood_test.py
diff --git a/experiment/wood-bilayer/PolarPlotLocalEnergy.py b/experiment/wood-bilayer/PolarPlotLocalEnergy.py
new file mode 100644
index 00000000..fae4dbb6
--- /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 00000000..aee5f271
--- /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 00000000..8e619526
--- /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
GIT binary patch
literal 18897
zcmWIWW@Zs#VBlb2h_bp9)|zA4rN_X)0Kyy$3=FxMxv3?U1*wSz1v#0?i6xo&dHQ8}
zDSG*d#hJx=`30$YDf!8zxv6<2dc_4rsfj7Y8L6oysAe)C0eJ=n2Iu^|w9NF<BCu)2
zM*5k#iRr1u`c9c8xrqhE`nWaAV9}hIoSd4IT9jClUxY`e6c(K&`T02oiFv6xc=Tbj
zv7jhFy(qP~I4O}-8;eU4OG*jolfdCp{I+4!m6nsASdvPr8w>IaN(*vR^GflA3O3t7
zS`cpJ;NSo!C#0BX-~lD4#JtS3)Z!AoqLj3=lh5Wg8wj+#w`om1cf&!Sm9<FY)XuE=
z2Ta={SoJ^661umqZ*o}Rl2xzfUwORF&G>%bmf3F;WlnT&T5&MMBw17Slt;$am)cgl
z^>!}3ZFTF$Hf{M?vLCaWm(JJPRi3A189%L6*PJKk{N_a(6|+9u`0r6Vc-5lNgXh?m
zJ`PT~UZ>2jPimwr&pqNy3lA)Ly<kG+cLm$-yH+a~Y~o7X7PPm0o5I9Khh9#fQ!ORU
zGQ0e|j~Y+=9B6rP;oC>Y6ZT$iUSyN-LA&q#je8wD%bvzP*t_9fO!NJksTo!8lB-HD
zAN_b=c}l_@qj&zVc3W5efBsMJ{C~+u>$jA-<T#vLoc?RV^DUp`L*3U0c(ZeC`JQo!
zm63s=nVEqBlK8O&A2>IrmL%#`<mU7S+2$Ws5V`YPTW8kQx^0h_ac+`b>Ac({Hfm+G
zOT&&W5;_M8>*|vKT-<W%LT~w<yN~^z+%DE_7B&yd2wkAk6?j1Os?Um54*BcPKa|Qn
zy*oy}zp3bqNpAMKn22KcrtXl)7|CNPjdpXg4IHkQ25#@!vQnY<v1-VH4&{fQAKuCJ
zcE!akxTdkDQ&dblsJL79=&pxZQH}LGZrEqfI3KiPee~@~XVaMfl`S%|`aJ8J$N40Y
z9Z&Y|Z&mS1c$Bp^z36|J;(HP0=EQk7>we|8Y_yh7vWXX+x2Py<`}wG6b{}K5lpf)H
z&BAN*?u6r~39rI-c1!RoJPq=Tm}S?q!mg(1m70y4K~#85q4ONGr+zX#^&fX9il$yn
zo}yRgwDtMl-&I@VQxA$J-j}%h+fht^UX)_wjAtK?>$k@{KDf4Ab3Y>}$orm!T(V_k
zU?^h37vx+F3=G93l{u-!pg_MHm79OtOyJ+Xc*PY5HZTTj+|CuR?rh!~yL-Bs#4Wzx
z5f_A1g_pWCE_u6Tf1P{Rq_1*ESI$lHx)5^wd0+YYyMZEayI0R{v5ZzydHQmxwdLp7
z_9cw`nNsn){@2vJox3~y^n1Muk&6vFPI5UyleGL-MY}CwKD9C?rO87eDK#g@W!Hm<
zpq(CtPEVXzx#vyX5)iwkR4naQjGTvM@bq_4D{~qz<ZZXUefHdrj*G{?d-PX!-O&hH
zbUk9l#$c&}QoWtM$6qmB;jWm}xL@0Cm77yzdGgu}9f9aj@tby0YbOflD$V=#&sfFd
zSgY0rs~;ODXXfZe&E0hIV^%?^wIWx4=(QWG-tjJe`F!0Rzh|E<jTWnxtQES(e|0gd
z$)n`IJ13s$uJHd?VRDV9-yuKhnXSg%GiFVa>W3Kw-L`bOOy2AnxojiH;d3h#jvBL6
zMI`2XpM2-}QuH}<>PL~M-+xQZX?>hsAa+fPZ*xL?sFGf_&7;{SXMP)Rch56YS|NAn
z%nIEsv01)`5{HhsO9elhBoyWn;Iw-g>*6>i+gVyl728&1TrXbr_4xCvXTmOqSD5Y#
z<=L@nZr&tO!Ns2$cjlK$x9pTY^WQHZVF43A+e0hYtq<B>XZP`}xRCItfobEKFr6#Y
zgx6fN&8_RZcW}Ya=XY<{ZthOeSu7PD{e#JI^TpF?ed|Px{i~|Vn;0HWT$@*MeWPHw
z+_K>ENE=IKKL4{E>_@mg4?Ov#biPuv>{IpIn@fLBb7GWCo6B#(dilc-%>|sP+o!+#
zZt?ijbT&QlU4Cb`9IH0;R5T7yE^=SGM$u>$lgCjp?yRD|cLlxcR~E%h)b%UmvHH8u
zv>|!U?x|fW{2$-V(#~5^5clur^v@4d*O&V~2>;CY_37<XpL>G^T|JccB$io+Mr`Yw
zCH95s(2VE<#ySN$1>3$EMf~ju5c!tyQ((bCuJd~)c(n%a{pI-AmV1(IYL)YaFq6J-
zQ(}*Z^DO_!cHF=#NSR~7>!}R?@9im=t!DE=Z%)~=gAClstnW_mlGqs5-8lVv?(JV6
z&z|s2=e)Ld`U%tATl($)c9`9K<#G6vQv89pR&l3$*Baf=ZIG#+e&LGE?4%3anLbln
zOlEXi<%m^g`hU>+e*Eovjw|c`SRbtVZnxmgl&c|mPxt;hQmdN$>AZ!_oq0c=8L1`~
zsk{_7FsPlivP<WF(XtH+$=k%=bM#){ps?P8TRU<`vxxNd<tzVZ2DLC8oo86{`l`UW
zCkt+H-{b46EvoCkvG&h|7t0S{^q1HC^-^&0e~x=LD;BW-VinUanvnZA<Yf5DiwU`{
z&F42>h*n*%V9=z^@h$R*p}*nhurt9<#givpo6+>J$S3-#(v@nT|8LS?ShyFtIX^$S
zbAA1<@IPKc>CduMFWM&w#TT!;wBEU_K4w~zVPW8nWtSdkZ+WKsW8(fJa)t(Xi(Bq!
zCyLGvs5|B!l;Cspi$j^omEFaA0{5_Ae12u2NL@fg?vssHe|VoiblkQ!y8TOqzzjM0
z%U%2a|EWA)Rdpt8kD9dVH`OAopf`ui+ot^d^JK}M1M2IAe2;rLH7)qOvf9V#-%%yo
z67`!ZPiJ$zF<0DS@wY4Ck<8@RCt~;)J*lgH9pB8LR68fUZE?Hlhhr0%`Dzk3ERhZ3
zp4KUR^-%5YWdG>%tJ{{o-C*9!exa=Y@7u+f%%qp(mSx>?D(>ZfVmtBha{lFUyj^pD
zWtwdXUEX$G_KL-o&D*Qyv1m(G9CbJt?qj3DzAw;p-~H5+N8U1Q<$wKUqrJ_<y)Eqv
zj<+li`<)kYSi9zP&dF0s6<?d|zq-6QylQJN|2^fekqZTr1rKv5Pn>qW_0Xc%mqf$8
z=1=_nv{Cuh_Y=E{s`N$ATg9x|-5g~pF>&|H{~G02H@CjLdHLtRoQ0PpCHX$TXaAA@
z@5!DEuWo$&@oRa?@}>ox23J>24qI}+PO)}z&E43mE6UeyV&<H-w6lnHvr}&7ivzoE
zvCDeLy!-d*P>r+V=K80zCdlmFvuoEb2j7KLQ+S1>W?s)+V_&U*m%}*n_Y;re;B)`B
z_Z_{LvHi2~ucu39tXIn388BnjM7tMESsyFJx4wTX?!Y!<?zS0c-?MaTp7;K}FXcJY
zmOqOB=P&9CUHJ2Fzphi3@imuaQ>EfsB1AXN%6AAc7rMS-r;g)2vCX~zN)FCHXY+rj
zxo!FMef?+FzmL{8Dfl@#*n%hU$G<~=rumgm4ZOSm`0BMyg(0B^^3l4_Kd*AqvkuZ-
z^zNPe-TVKe{;`89t_40}0*!184EAdHDlTpY28QJPyb@5eydpQ}?cC_%$J<1Y-~X*?
z@@~syy9tN>K3Zqa@A>OO+V2a#`u=Ax2&EP)bxKY8^z`=k1|g3d+SUD=s#DIIoIE1w
zZO>58^k8QL!{2we+piWCg`4Guu3xWLaZXR|X;s+y)nCtEJ!)V7%lnz$zS{r0%fBWC
zo;`KC()RYcQ>D|-zA9z^Ki}bHY4y&`Q(|@6yY@}W-Mnv2_qHtk>ql4R=<mLnx>xi1
z#y3xwr|SJKE1!S<R#;iM`1aWCJLmYlfB*h-^~~#M&wg26?{m03GTg4@vd4Xs&p$s2
z=t@k!Jp1LI7hP}GX}zzP_qMrHwz^q$$M?Eh_0|h=zdBv5`8A*4zHZLlX{(i1U&{?&
zlm4`Tal7l~#Oj6Zmrv)v&)bm7*jgBoyMTNB^wg-6do<r|+{z-xX03ZK=NKc0YFV@V
z`dI(jMZxB=OAc2gNOQlLy!x%&bk$!~%yXv2O^t2TmASFy&EZW^|37^);kfTs5|}Jq
znKAw5qPoO@SFGD+RWUw2VlZQ~<FmBeI@J?D-`o;2ZJ~$tQ%m{0_Fe1ZB`r<w{<ZTH
zwmj}HZ0Wx~Ezo>*d+o!83u}3w9$uESDqu%MlZwWKtDjl#XvE&UGW+z#CiU;N8|U;I
zuiPn_b9T|ab<eGh3m5zfc=O4V_k4a!`0tj@nF(#HSDkt)QyyBnEo**!kjCGV6V8=?
zHhf>J`gHdrsru?K1yNOdV`HcDRy7<9sw-+gn0jx&)xvkPxvCTnrmugL*fB5W&THXo
z{?!`OBChREYI;*`y7FNsr{=EP;%gda*EQ(>mY1*dS#G<{eg63mW!FNax2(K#(O<22
zZ`DQF*6KBTxW6SB6fFGtt(7Th*YBz?C7;h%@9I9ePk(ZOW6YlZjR*Y-tYz&Y9C{|`
zc16xSd-~?4+W&dxd)xZeci(qbkBR*gaWLjBd*P>z5zPi4@≻Paj|1Z7()^*Wb^^
zdEXx%ed2kw@qFGXQy-=oeQ#^8Zu9T(y?Qp`%KNTmua;l`Vtnb<?NzgD;_r#Re#CZa
z&ziTA9fIF)ls$dredX`ocdzzotKI+mHGgvd{k<#RS^l_s#W!B+{om(Dx2=9?w>m8{
z{=9+UulfG#Qo}wk+_?AKJ%$GN*u3!9-4RBA4Q@9u>djuKb^GL^eS7jwJ(~VzSIbLY
z(Y9sDz14sA&1#R1d)oH-=Bn%SPgMHOwwP)*lRLZbvYP)=-zPyQ9&cHBIaBzS>~i1v
zc|KFE+@d~TT6^Kj_wywO)mN^#8~*U{?M=IOvbG;CTFPqjZ`Zxt=(MlPR)-xm^7bx2
zyY0y9cU^~fJ#N?LFgm4nTf{tdoA$@rof_Ltzq{$GD!zp~=UM3PNU{EE+uygGelVW&
z-nQy_z~xORcfIys>KDJo`M3VM*FPV>IHq8z`|NFQWZC=sZ)foZo-Vs-J$H@n?h8pF
zrq_0^xb`_GGG@AW!LIJ4n@uNVg{`(}<}YYdYreqp>D~id{=n~7i|n*#Ra$@h_3re)
z=W=(x-s1P{`CQy_w{D-_>6^LR{zp27H%{-H`dxfZukx2A_x>KSYxwr9;{M;=fv3`*
z7R_AhBePF)m!D(6@}AOnP4{eXIqlCswMbiW*F~1KT1{4(B|_i$BD$D;*<OG3k-hRg
z#-PTeWy;Z_=(M|!*|@gHGMC>ts`9LueP46eA@kYa<xeHq{tRchb<+H`_Rr<zm;ai-
zUv~bs)rDW_RYz+s9glbXbbMVy$<mqr&z{C!m7e--b=S4MtTkqC`|B?3`%<)j_4V#o
zX}h~`_V2!+pR;pjjBmubqtPZ=_ZI|y6bsk7Ds8grNwmqU{;VknivE}_te<^h`RohN
zA3Obvz5e;+s&%*WZcmx~<HgI__BFBT=T-MUj?sC&Xm`aY-X-bn)rrUX<{sEr(I%5K
zC)(r1USprgDLpOf!YXsN9;)gS6{$3GypxqWQ!v7D;&S07GajB<Sgr9XrD?|1Thi$p
znw0$2JF*#jHMT6U^?&%7|Nd9YU>^?sJ;(0;SN6Z0d7)wIrA*a=#WOynG%fjfD5LY;
z!V{A@gPby-?s4Y7IZs{qNw3JIj4u;h=lM)XbF*G;8y4{7zJZu!j&s4Z<`WGo?-=Fn
z4fbF*`&81jrOEpA12(gVFIhMrOE^cZZc!J0GPmSR<dy{|ggGx+KA$A5InmkdQ%aLc
zi=tQIlO=LLg{xBa<~)mwTEi$^VYU9vTGLyLb+ewVWGOtSFi%8GoKM{*lE3l#^ivjh
zws=)+lDl*)lK-K~wss|dHIIT!9%uF-=K8!iN%Q9`g!E@hFP!kK_V^^ujWHbA8<y2R
zd0FM=X!Jm)NHp=uNtVjfTR+K|KNX&*W7F{CO!Tt}t2vbY)jTRwml&T=l>RAP^>_nk
zfN>j7v+%rg($>noPlew0ez%#K8Gb73_=YMaziamA3Vuvq82s(T<gfjqLZ7@-Lgsxa
zIml=(_v$!1d-cCNHHPozE%EYjUDxv6%lJ-|&!NejOy3htAA2y{e=4y(qI&t%qw0+F
zlQ}nq8R?sxi%_+RoN!>v)E}YKCaMcB=@GoA^+ZPMXYj9-OqI=H$-FlXT-I&kO6!@j
zLv?qHy79+J+AB@0JR~ZN9DO(sE7Tq4>av-s;NdzYm`PXTQ%NFI%Wby>?m7MQXZs&>
zwG7GBJm1ScdF=#Y`-dl4EO)Hm;%?XHQL#(T=oDxFmQYUrXOo?eJ=yE<?ZS!4oRgF+
zYTLv_DtE0o>iN9Ig^g3W-fR9E)9{|hpP!t#{4?eKlIcR*cTD(DQn|^%Cgt%9m&)^(
z*PfI*w>AAm>X#FfIg55jOW$>#xLjC8)S~vw<o-r`^YuO62Uc9t<13a_>pQ;gY@f|c
zg@-dgM=h9}@Z=<mVZLm*m@X)FhFDZ*e@t|Ga*`!7b;f0xozo0HrNl|@*|$w>_dfei
z>R(UY;W~Xmt?t3gJ)4@UZ%cmpxAeqi-L|ZihvzYWNqRAF`R7yDsylivG%5K{pAnU;
z;Ce_ekH==F!uy#Qn-&ItKQVuHxSO6`z*@^6F~(^+!Y97oY*{Xp!os}xt)*UoM}?N7
zR~EP8vcr-;ThygDZIIVo^EaVg=Tk_NjQAz4)B?xqL&85(3zF_+KHOj}?DnsvZl-p6
zk8RRd!%rptQ)dPQDY+XyJjs&zQsTLgtcJ_HcRPzBRLoC=E^JvI?Qh`LY<X<G`kTsK
zs-<s#FaP~Q^{Z`tUR~rvwYf^g%{`gR6xt5$I@q{e=!m5G1D6T@tito$CY(NaJd1fH
z<GOaCc|H>Yr)l=znLF=8iL+Q4>tg2V9l<9XTMl@&-JLf1gZh_K37cfYCoKkrXHDP4
zL!Xq-IlIpD5in<ZrdyT&>%e4Ap<Bl9FLb|R_Mg5lTXFAEyV4GexiTMk&pnuaMcMz+
z<oWX&*z|1w>G17*qt41FGc)Y@*L4$bZriE-%lXm$HUDORG`5%jH}l&Y*1AQ4a}@qN
zMf3e{C^fqg{Iv4?zX=!LGQVZD>@?$L_Sf{N=rw;fQC;N7nvEwUPCa8(l{0)+YjLWD
zC#q^%6<_1=1L_MW+HV(moy)(y?Yf}HS=OSth7VV=7(U$j?ZlPJBPX9^9AER8X(xBd
zla*@wO68;iUMSZy`>T3Hzn(jtZFiIThmfWv+n&w1E5lr$v|MNH1+Mh{Dfc!XVKbFE
z<z^S^SK;Ltmi<I|X3O%Drl^d(0=4Mpw%%`9ELDnb%(VM3C06BOoRePv$vx@zw_X32
zFOM#MY&rRy#S^*yiOYqS=+4l7Heq!br@xwqrtrNNO?SAgKBY{Y@*z&*N!e?U;JY1u
z6=L@mwv?V?_J1~IE$4eS=XHnF{&ddWWYD90eq-mJ^p(6ecdnf{ZoujPa>~>XCo<=K
zSW>BQ$H$%H*2gJ&7cwN+?yOD!G;gWjz23=eJDV>)v90d&`{~=ivnxj9<<7?*@-se!
zG>L4UCYkeU>W7dfmF+%W-FMXHeF*Wp>2-W*V~+64Y3%cIY#ATL@`d-t=1W#=@^CCZ
zHb?Nr(Ho8SYm=WjeSF+g^FVWN-&LNSp0*D)|8DvowDA5{fl9CChg+YX_%8EtPx{O2
z8+hH{P5WATe$w*yQ70B(`hH;k?D)FNi`OY#m)qo<IB_<!&^)gRUww+7Ptcyg>92aC
zWxGrJ<hM7DsC^1)vgm0{;QMJ}_i*K(O(Ic;>QC(b%jusgW5X@Y#QN##ZO-~fEs>co
z*isVaeb@P_W#;1ZBiMA}auF3B;o}z1=Bbx92Ngci{9N~crb7KEhRNcZ8V^0if9<{E
z;cE5%X0gwV`(L)KK4B}t&ar5M^kz^FI8%Az{Q<v6EUGq<{H=4PyR35#*?dTm+`~FK
zZzKC&;d!SFS3gKHIsWfw!$kGd50-8>P~D(!E&Wqi{mg@_AJ*6P+pO&6xBO{nX7^!9
zlgz&<=jE1{SO(aBs+@nC*=z1Bsb?45v)JZsYCd#L`3X-@<@tGw{p%iYUU$LxA^#eS
zHXk<Od0qnLVV)-^tatIQSY>dqnB&D9r+2NM6<UsGc78bK^F(+~%krWoKQ51_J?`%{
z`)nfrt?NmWdYvHj??_yvbYs0`-4oBr!RtQMf7AFhrD@HT7tNEI{eQoi68U7;gp<|N
ztfIpNq-8I=+WPNWz@RDgNpkl2>Tfc(pQY!Qh3<TkczTLMo`Fw=mg6jEwl|VRn+p#G
za|*ddoRulE)O~m|Oho2$Sa!}5fx=Z^Y@a;ht9_*hiWWZix3b%sh5l`+i!7+~S$FyI
z6t=5cy}jRZ?rpwecU)_R#_=6b4@@o-ykR>>azmxV%K5s{GlX6+N85U@^(fSvz#gu!
zlIzX-52s$L2Nr%f_IXXs4U^Q{VPALgE;s$o75!twjP2V~Z}!T^UH^Ri<929w>+S!2
zAAbJ_k8hrJ_9%MC%fKK#8DBRMG$LP|T2hjkmku7`d>eH;@3EOc?fY<+s`47|UKfu&
z1$#W-Ub`zXZSuC?n-(nD=-cv$L*?n|x?0bhd_rrtgyvkl!TM{Gi%?~niH(i1^V464
zt&-Bc1tux&FP!z9_c>2RgW}WA($5R4;&+$WmDKM%cC&}=L7Rk%U_t2Q1ykgd&VMP<
zyUXyi=M`t?13e+nrKdevE`>IB8-Dm}^8IUux#_mOo71M)-?!jinzTJGer?$E)!!_;
zD)&~FmWKMzf7h<`=(I!g-mf~ccd|_MUfq5@{o>2e$m3=CGp)A${vDbT<yS7>{`YnH
z`ttXi{)R5*QWBoM`EG9Vn`0OCc6)QJp1#7e^7FlOm5V0*;;CXfA(EW_KOi^h?}VbK
zCDPZw<~!bFwP~8UsripM?^PxKU3Yip#L4j5Gg!qx*;QA1RgX8m+aqSNO^azuWblSl
z@k>hj)Ov1ADU=91Z18w#h?DVwmr`~n_n)q0Jp02^)k@+2(-*fj%U?^dn9fmu<Lx2x
z^AR&k_pNv8sVR#k<C_-zb@32Mo)Yr%@5l01o2u?V{;emNMZRdy3jdzuc6GjlQ(2{$
zM&Aua1rgqyh3R|lv#&TM`{!0Y_m8~}7gd+M5I*#r^<kzg+rIDHlR|GF_Ulr4?r|sU
zYHfd2w^w3_;)nDj8LI-dQnxVHD(ReDw5-oey8T$gnS&Yc%P(59H&t`T{8h^cyH@MM
z)?yv2Dr8Ym{K0`;_>ktBtIc}WX*}V|GbSelr)0a$&@q@+{JHd%R@IzmGdBIUIBUoH
za&62k%P4l{#-5)pBG;tPp7wjG^V|MmH%C6tomabFSLz%P(PJ0b)HCzEg2C7Lu>HGv
z{;7v-TkzJS$?>hS^F;6IVaxAk{k^p^AyStm)nUmuo6a9B;misPHMMdjK8i7TX-{PR
zvarTuN}lb(+b(uSFCESPhbO6|{&5HwKNBJmW&3`I;58>J*^MXLE<IM6_o&5#H8LnP
zF@bx{{mk2|BR3l`t~?ofqas7qHsI~14?FK3-}oW!QG1DJ>cZcmjnA#?HA=UA*lRlV
zp4)USCd<E3XD-d(&39`q*AxB;3r@}dos+N9Bw@DqIl~LRPmbO7BE}J~t=f;hio0$j
zmZE2RW%87Q-)U}EJn9CY58Qkf==392TC^+Rs8->mBgvt<I%0b2nmRgR-DTO+Hs<X=
zv#H5z({V1li5GV6SURVpH~7D!(b}^?2g|%Kc0UeDJK}!#{=9_dut@DWdUN|Hv<UC@
z56s@w>iSNV>vsz)r}C-^F7C^Oqce{vd(Y6DbmDo-rui%N@2)(fUOTb=U$9eLnMATs
z^A2%|r(R*Kj_>q46dJdsWN9*9db8xz`BJwj_V&jXE=%~hvi|q+oGZMX^S8-c%2&)^
zxJB{kHR<yU|1h)ouf9=uh<WOSIsW3I61UcEQhXy)lAGsR()Q``yWLxspUi%$vfy}0
z_|&|E_b<K|mOH@66mmo6>i(ZTar(d5o8IjQjZN+EViLDzVPI%dfE7BhS>=lxS>A@J
zGBAKJX!0;5qck@uFEKNxSic0u)GNqK?+oyB=jD<D&Aae=dbk8JFfd$YU|{&q!OXzG
z5V$h=Ap--`-vFNwR|W<~RxVLdF-ajAO&Jw^F?CBdNj@cMaceCJbyaOm9dlh{D@$!f
zOGABgQ!^b)FFo4;a~oF!`#@_aZyWbu8~<caQxzv`^8j0IM_X5KC!0_=s{~&g7iT9I
z7jJJ@XHPc|9}hQgFJE`>5TBr^Fdz3|zksM9pYX5<uZVQd*rM?0l;F4`zr@n8#N61#
z)TpHFsPyurK<A`TpUhCl>?ohq=<vcg|AP3Cr1*sF#HixLkm8i6l8o4j^zh2G=$eeM
z>a6JGq?Dx8tklfhtdyjjjMSX$?9`mHw1TqK!kXNIlI-H@((Kfd!ra<|)av4#^78VS
zyym#V*5tg#^s=_>+KFYgEroTH^BSjCHFlJ@Pp|A<SUYiP-Q*RGc`+@8agBLNZDr}L
zH6@ekaylBTrZ?qIZYrA8S~az!W@c;gjP{DT9i_9os~7fGE}K}}*4Eb6Gr4Qh>?u91
zvnF*-pFX{L(&EOcD_du->zul%d&ZLP*{j=UZ|I!2Y3kf1eG69fF5Eh6(W=Rdw@h5N
zYs$*KvzKn1y<*Gk)w`BYZ(llR%BH!U3+63axoFnbg}s}XOy0C~=9ZPScP^c@XXVu0
zt7k7?zI@4=t*bWe+PY@RhV>hFZCbf`^X5%kcW&Fcee=#eyS8oGy>s)fUAtziJve9M
z(RrIsEZ%%*(Y8}dcN|~3^UUgP2iNR4vTFC~<$Ery+IMl=-XnVt>|Jy4^3Fr2wx77Z
z>)gY=7atv5GvnySxra9|KDcew(H+Z<>|Aql>%ucTmYm+R=E$Bc7x%6_e_-9|Lp!e=
zTyy!zri;gRT|Kew%E{d~k8ZqmV(aZQd+uD^_wekF#~1e=J9g~InafAc-8gyf(#Z>F
zPhGlo<;;;=mrmZedE?NP$A_-HJa*&hv74_>Uwd%&=Ht_MUZ1%4;mVyScOTw6_u%!p
zM;~uKes$x;=X)RjKD}|`(e2BR?%jNO@6w06C*C|b_v-P@FOSZBd2#dS$Ggv8Jb(J;
z^P5jEUVr}a>Fv|+AD?{r`1;qUCqF;G{`>Rom(QO+ef#<O*Y__!|NQ*+@z2jM|Ni~^
zv)=X~1A~%?r;B4q#jUq<c|(H3_r2Jgf3N!cxih=ZZ7!Quel@c;T!61veBrqij-{)a
zyAlOCbF3J>ld_pLUun-0SQM(|ACRq(@%Mpt$(k2OpC8vc-}1|_?}D3b^xB20^5@oC
zJzx=9?OyfC&~V<DCx5QrR_EHfH~Z49dwas}pNwLkwE6qZlkeXY|BFw5|9>6}|KsdS
z-+w-x*SP-8iyhBvU!0XUxT-JC{O9sL{kK*7SsqON`{V4%()e||lRy91o&32dTU-BM
z!OBgS9yCk-ukulP8#b*=hjWz(Tf>LVER&at$wV*z^Eu+*rAsl<j8Prmx&O}0dZ^ZN
z!O^JZ-*)p9`<PWtvkpHst)B0s@znP07hj_jUzl@lRuq|T$zEo<K6ux$OKPi)M4jex
zH;6hsJJJ5>MdIH>y+0oFPw$%3-}g?ms?jH9+LkE}&omCK-J<ip_vRb@=xdtxJ14p^
z`2K#=GJWx>PxIqDe!ptC$nUx41-sADJrxzt{|SFyx_0KPMU_Pw2WnZT8trh+*i>{)
ze)i&NnG<h>1!x`oZTvTAee^-utLEz4zwLT+a7tr*t!(a&@4-K(wl{O_dN=irT7}1o
zM@w2?8vMN~^@Z8<sM+GI|9Aa!J|+B}etK;|%@g}qyJc2oJpL88Va2QcsnTz?jIRfo
zZc*N%#uj?i>eiyd6|0%qa%1oN`?f9Eb?-v`)V2953$w2Dom?-yx14+b-SZjiBw794
zz1hX?w3NS0T+?~V{EpkbH|6^(h3iePFm}(DV6g~)q<p9#dOC+rZB8xEs+lqJg8$y8
z_O}_gUyeA~Rr}68d`3R^Jfr#MbEk36?6Md9C!5;;%Xok6{Ii=F1HK(}l*;=jcJKL-
zZx?wP1pRVrFN@D~0GZ2{af{u{&qU~ldG}0-kj({(t_NBxe+InpNoCcK-Kl#ZWaZ3%
z%UnZ`I=NL&)Yf5J`Z_xB;d}Poj_Dr`MHtSBDsAtY61S^g?L<wDR8$F@j>z`+GZ$EG
zGQ3v1b?J?`V6`ow*X}O8bp7AI=f11^mvj8R`ApYeEq>h(myFeyL^q$(zCY{K^AfAW
zH_k}cJWCP}&=yUZbn@l@XTEoyTrtdRV%jk=O}#Bw@cCoE6LOuJ%x~F*g)<CU*Cu?8
zT4;aa&Q!6?jo%^~v*fb2f4i}_ahlw^3;)-y__e=6+JMKQ{zlS;w})!Yznq<YApT$L
zk(J9mq7NioF#LPw?+fknd0!abn|=!Ue|`m%8pl_L+r}H8w>`R7<Q!rdT>p8}cfUI(
zcXG?uz0;L$$yp%#Yu?H{^+OZ?^DTeE*UbAXI_BWC-0lyymN%MQStX+P#3mT><oe1U
zyQFksb9dG=AHPjU+g#6FxUt~iw?#5KUn68CSMc~{JkPt`@9?^5<Dd7xxz0#v2K|jb
zs&}x!t@n`P38R<Z8v>Q@&r0%-yqMsxo@ZBiw%}9cr|kh5^Q|;wR@F$a&QZ4texvv+
zXpV9Gy>hR2id@<KQ}<lyb}jhlxA~RJ%ru6UGv9r-?oGVM%KC90)2d~An2Uua*x#+5
zn!nrIVTvJ_F=v=+-QT}QlV_ZoB|oEQ0%x9Bvh04zP1Yg~Av0G82~LWWJRKf)?9pbG
zCE1$6KLu>8f7u52+HRXQS9`ai^rq9-**Onwo%>pQbG4v=X!~)NBQrYH{4|dRHE&+C
zN#No39SSZL*?u<X)#e`v(Y|N5s-r~1a#i<_X-8faz3ey@&E@p}yDszeoeuv6Wx~td
zrLLNlXL?H?czHTGH-RI|RPXTf-dFP&icGU+wMQRZx9+f6lj~i!vsssys~PM#eW*99
zc(?hDs(n^GYcynDH2Rlr3%GSy<!1lg(u-S}`KRt%R4*}4EkF80$i)ldr@EJJ+rnAE
z=DFP9!8yUx_M%<kZ|6)eUMX|-vsw4^K+B+<M>>&rbz%*wZJsUhVOJJXXP6_NlJ{n^
zhxy72!K<}ZH)nA#soe1(u-Y$gyLF}io|S0_)+eo48_hZA_0govn_Cp8Z1~&iaA^Iy
zyr1!%0f}NpzwRD=<CE||)9^s)&mH=X3AUW`|B5}By4O5dt13U~&a%T1+*_9%(9yqh
zDt_m~s*gn%k68rQYP8Er?>n|IrFe(=bJwd3Q~yUtZr9n$A=o#cSDVvV$8hT0XdZJf
zwN?E~9|unq-|9VeZ*3DRZvjjBnQwUvi)ETVCsuoY7D;p8vvLFPrcZCzFZ|VDy!5T|
zgg3fI$?rcL*;8_3V;DnphwX%gQ8Rw7D6f1kK6kd<*AE*#=T=J^xkW1dQgJ%cuuo6Q
zY{A>+rrmY*HtQXhs^}i8HEpZwk@#C@QTAd_Mx&3t@Kn=t^Nm9qiekJ~k8B8NTEB;1
z{BgYfLss!WOL*e`aKBiu!5Ol*xbEU=!I~&1rh@-0is46HE!0ogd^xtJ<=y*8wz|6a
z?Z1uMXG+YEh}-lzgx@sOL;7>t$ApWlV#fEpLZx;^8uPqmR@}-JyJ>ax0?FO;5)RHQ
z|LnqZ@Lt#dE&Pd&>zuOBD}8lLDD^Q8+THU?&OB&8pW9b~m21-<`7HJBdun~}nhMi|
z%{vY_Tu{8dah6Z*`ah-|Ro!zq_pM2@v=u(FJgy=6^0E0(CP%F)J8!#9Uco7P)#@`P
z64U%rTmC6TzV%*r)A{mpk(B-qMVH$u;v4;sf9v2(=>8SvUV5P86Bmd3&WT*hvyP`7
zK3&Q(Y2%B%M^j#I>i+&>wLsPrHxmcG#b=@mUfgZFb3T2+&vgn>E5tkE*Pqazr2gAq
z{P^rw=JD=J^;dfypCuwKFhlF``M{^Y77C_%iMh9&7Q4Dc{d&5eP3T`8*(q1O&NWnh
z59M(BFh9NQoWzHNySq51rgA%Py2R7Mx-zA3*`m5q!#~@-BaIqAMa)h)mw4fQ^W>uH
zUHQ8Ym2=;2y_J}_c#hYm$Za<1vS;MvBZUf=1o34*edMTL&XwgeEo8On$+kIDmsTDw
zHeLH|%4F4VTil|L-Y#6W>S9K+jXP(CNyMxP@|JoBC3OUQtn&A}RV!R+KH7Ay&Hh$_
zsE*;A$CnN*zIe^>K|$o63-8+(9Bi<zG_NShO}H4Ttjosrt2EB*`eE;~d9NNW<0`Ru
zzkJtG{*dFx=fu3&(zv+y(Dq9LH}207uwc|)_gQV<>^HKywhuNsepu_sQO3uzWBIvx
zX7UlYGJdUXv<Tkl68@>y#yf9b;l#6B(-P`e6->W2@7{cdgUiMK$K}~rtZTiPb7SAB
zj=5{HubHRx9TKWowswQ{_p@*ER<~}@d2@eD@zla^#<{1Q3*Tu=ZF={d)3R*h0^8%P
z+vUQaS1M{Woceb3=F*_cw!gS+Ze_j8Z4_TO{lT~JJ-*IUCCl<+ChF%t;NpDdwCAg6
zbX`2ZuJ%8cgwSt?&Wc=kuex-4*}7NNb?dYzROyReJ!hNyQ=t07_4&*f)*rkb)$znY
z|2?}>5BK(Gn)ADjLl6FSEto3K+;Po3AZJo}Wj3Gawii=0PyK!OrL)7)YOy)H(29a7
z|D8D39_6aaH0OOMWU<S_;mQnV&$<7W^?qFXD@V?0`g`YeafQu-nM~OS`d;Y>g}h?_
zln@YYW%KUg;!3HYUbbyJ7hL2Qd-*3z)hXANW9D7<4PX20xK>;^P`O3zo|*Of^%Yw4
zwD{If{x7sTWIj`4OlHZE1KqPv#uUEY_|sBP)SP9H{?lW^E36_E*1ny*pzy6g<S&^W
z#bP;3H}o9V9yqt@rRrx*QLnh)=U0fx-;-T6r7=)vmqfbOrV>A$+;4^lXWGVWs7=4f
z_WSB>mWSofmn5H+a%lSg)U4{XoPgNg$8z$)G8>jR-z|+k(EH)|w207kGZXywv8}l=
z`T3W9O}1{Q;&|Mi#k$T}*|nn5@bu1!cYU`q<o(=zgY(+9eOFfbdz$^-#aLN<PSH7O
zspO{FtcTuJZ#B5)&k$`I6&c@rPv_0b3c0d8wWDlKPnZ_Y=RU$Nw&Rg^X;g#S?Bwbn
z0Vn@H^YV1#3e97RU%v6k<;aza7oOzV|K3t2pQ2LL+*JKYY@*US$F9bik{Q2>Q=%_t
zD<n0Egq*2i=Xka`=|{+;uHUsskL_YQeyVQGO|z{&$w%7a?q~b<uUh?8BYnk&9J{c^
zi?~}>Rrs_#@#PL({q9=%;*(cHJ72w8rXBrx$?6@2Y+5GArG*$YwdN+z4-R=9@Ndg6
zaqg?^kMx*}lcvPZFnRetq+<0!wR5{)$Iml)>Avs(vWCm`_oZ#ueR9{k@oB#MjhFpB
z$3@D0W_^j@^l@1Yd)-IlJ;kZ|Kfm7<@Hrg*r2bXw^!jf<%VK`;PgWECJZbLosvWj9
zZx7CR#iJcEseA6>nOW*vUhg?8Rod>hdSA-*H;2@Ab!zjNy;6Jh=7tG7ySPNR$?H9i
z%z6pD7uH?l$u_mx<@)?<##MHvUDqmtQ~64dNSrH)p39-Q=iz*9<>=XBw=Xw`1Uolq
z6)3%WY0o+@sGzm#b7!Hh&#Lsp%eDN(_G~(DyZdB{Cfi>TH<!Q^IoFrj=bxG9F8Q;6
z)7l%t@vAm!Y90D-`b+CjrgH79kbQT~oJ;sOt9`4C*sWESbxU0tP1aqQU;QFZ_e^hj
zvpv_6*9SlOe193bbH06>O8PQ|pQqEl)rTzY;#d29{^(8deLq!9{=U8*;QIXV+Ba+u
zH$>P@&b-Ax!Mf|}njoEgzwgDZ2fwF$52`&D7&^m~@q!rR#i=QaA1rBVQngJgWpdrj
zsjl+$V_&pN;Q1r73SKH_a0t|0D)n;{`qrFYf4X=7?x}3aO!JF;(}P!;d;H+4KJRx#
z*UopbWO>zvm-QEIrY$nD7OLM~eN17|s#8xK6OSL<ylHLw-KMng5~&QEl{T;DMTlpt
zTy-V5%Ve5&LP&hCN!ag&&5Z%l2O>@PJ!MuCT6uL!n6aL&%&JSTK5K1xw!?hIWghly
z6OjkZGCS4<<|dU|+24Dsv~}&Cl;<nu@>V5;pUHH1qcH7FOw3$qM*f1{rfj_pC2zUd
zckd}U>(acueOvk4UlofiZ_nkupA^RShWS9$+z1gq#;Z^HOG<AY$mDHq?=rsnmgV&w
z3BgP;RR@lIhWe??5?&hAHoHb}pUf{~oov<l@5;9?{9V6=-|qhM`0(aqOEx#^JW%=l
zO=z!={`#oci1-AtcUo)Idz#jE@5*r(H(Nfx<F|*E*=A>n*9QwyVz`;jf9#0a&;IX&
z)bAa;-d#?L@XTr3`EAd%KN%jiGOHzStPoq)sB=z&_pHc))=7(|HA%6iEVfBWo;myI
zo4$(R)*>I)>H9l#5Atrjwd>~zD~p|dZx>F_*K;Yp-gRNi{i)y1C;P2!t6$!FQ1RP0
z&XC~e2YUZB&0BhR=F-c%=Qg?=Q!YC=k+*Q!u6cag-;TfXR=O2$a(3mtGGDi@|NQIQ
zqGn93-mOrbpy2V~+-&Qr^D$<6I)!^@@EX0dkviAEBb8-ut?t1&Q)Vk2>hpFyG_n2d
z=C>c0H607`)zR&7Inn8|)NQl$|I`O_cJ58M_qx~Pu*u@<k)oSeLsy7qsp)vOhlHd}
z?=tBP|0)zx%d7alatWvO#96BxBa)RbhZ$Zn(Vubuly}*@uH{t=R`hijuRFBzQAAU`
z*1ewA1S6*@A-AL^bn#CJ+Q_-^$>A!`%mA&m+Y~aJT$_%~=boHo_<bAO$-9s4G|O$=
z@qLcVn=Enf6`wX4{kbR-8O9VfMe=+{3H#@DN(!B#Pw#AE<Y4<Z%UtH#d*+Lm4zKaJ
zb?s&Igc(Pq7M|nrjL^}VBA(wfe`Z(QnrjF4Zb;rx|0FEyn}@@M-M@Y3bM#Gp=Iquu
zQ&odIJXG>}$Cr5Swxr`X{ZGF6Z_yR&5-Vgr!TLDsq0F_uTp5qlw#=>Jy|sE{yv#f0
zj*k}AlMflSHZIm|_N-fB(;|E3Nv7@8?8^_Ed^CJ}#frGz&QY2duAsT^S6$K@i|Tsz
zciPk1M2*iX=)7lDd?-7cSLs=royps)*9{*`b3XQB?MDmi>^rSjziG1CRBt`O61iBg
zNxDs;;rnFK8_(bWkrkZlyRs*H&R_9Vs~YhmyO>T|Y8_3_W@>Ru5nu5n@6XGOd0S>%
z-)EJY_;KQ_EVui-zk=TWm~x}I(sj1%UB!iJ%Z+@uuhd_^!ZChU=DLt>AEmhq?AK>B
zN1HC4w|m874_99P<oO50F9$Ac`0Ds1YDP7?-5%C=GKVJpd$9FkB1c=H^%e1lu3SoW
z2`e}BXU|^k{d1ADtxvIhW}f*`NCaP4G0pCw@B#1oU$;g3_DHGh5HOBiaZ5lX_o!XY
zwAWkaxUxL4`2AL`N^M7LH{-Ux*Q>+%dK0J5%vD@)$!D46Ebj#mlC`G^K6Tce^T2v9
z*RM(6uU7OQyS(j1w#J8aPoZQVJBw8ue$uVM-3v26?RRwZy~Sk7xIJjA{L>y^GXcH+
z0KdhTyA>5A^OKIo6>$lFU9gc&?ZGadWiPqzweg2BrX*a{_|`u)@9y!L3Lh79`L3Q9
z7R<9=rZ3Qa{Xw^j>|7B)k62H5$?BmtV`HX&h`*T~=k0jj%f+WJPu0A)tX9}J>ve~9
z$jrroT4s!T5BQ`R^RC&yTeaNVxU1H6Sw&o<oz~@v1~UYe&-qQ;_MH8;(~-Q3@Agdl
zlbiZ><<zZHrJ1Lkn14|%$LIcBhrrtN#*xAw=WM=e%5_O(qusr&o7V*T<<)+A+*S3o
zE=PLPHA&llPt(?Xon>aT``VM9)mt_gH*TBpIPa~k&k43hKQ5VAo*U(c2@lys_?$kh
zRV;XPS3u_2TAwT1qZAE9UcGqn;kIwj^2LXbl-=C%YR$&Qho#T1bFTVrZI!c!=gG<n
zt$AI`A7P_@Jq)}U|M@55!cBEAFW-D+{om72d`{JRzAx=^3|j5G^dk7*9<~ZKpR?o6
znmlW>qmk_DKD&N>zY>#IvGD$a%M#w2SKmeTq|K16UUrVZwd+G(@w3LF0xi`BpH;~g
zCzW^4<<c(EOgUAsQ1Zj?m&)F|`lqdBoN6v4t2EX06GMQjd+}bA8GK%=&*)}U-kDq<
z&2XCclXNg|n_}J{DaNVd_v;zM*J*w&*1R%>pW*767bkYVX{^ZzT5tDLOT{{6#;er1
zkE%-6{aCxjI;8yg%fgUWeumHu$qa8~o9(nW$nr?|^3A^Z^o2ypt%j_JkM6Z9hZR&X
zmpeUMvT$dq&U?iPZ{~;9r?&q)<?L|uw%C9BfQ^Umy}xtJcCp94i?5|b^R*WGeTm&z
zrWca8WB#wuPfW!%E7zTN*#F;p+eE{v`p#{g)}PNknko71=kCR)_j{Kc-}&(&zwT`J
zjPq&7r*pmjHoba(vUIE7rXTyxHT`F<obz6LugkyK*&jNNe&1?vBy!`C(~u1GV8xxL
z8Iu?EyQec*7Zs@fT{yuj$*%c`b1O%s?V-Rh)w?Wik|xTTH}@Vl(LGkr<Wbe>@~Gc$
zL6F*nYD*2L)Qah6I|Am3E$n~fliIUW;pW?IE7V?{kZ8SpsPxyn2`f}<cKIA&C`s$x
z!S4{X!FX2F(>PZ@sl}TYROc^!x5sMQtRp%1H|#!Sqn5nzlfm@OPrh&6`R_-3V`Thv
z=dXP~1srN$?>w{N-?QWPUQS1@zG#X~Pr1k$)u?Y>9C%u>S=y8Rmv2k<Qkhj}p7^M&
z$ok;(dGlulZo!Y|wtW9OGjV75viO<mKAL+toTSY&dHU8gzX_~xu3O&p@3iE-#l=rQ
zoK(8;DX$=VRa}GUoV(vQMTYwuTx_j<oNL%vSXImO+U$nmTHOy;67siQ4Bzjo(p$SM
z#j~_h@}kp?GXM4Y{x34PxA?qF**Aw{Mo!J5<1ehH?f3cEwPx}1J+3b{rD*j0kVt;H
zDKuu)uH~<}q#yaDZvA#-0axViHIC{Lx@z|33+Dck@aB_dSu1kGH*ERIq^}Q^ete49
z)37!Df^b{5*@rLDxf31r*JT&7uBp@8EL`!yHdFqQs<r>ug5bG<+mCHIzWAW__FYqR
zLv|mC<7dn*y|HK=dv}%aJ)fM-RSFBP`t10T9)G+(PjSn!j*!&_2Mxk@sx9Bu_cDC8
zl2z~I+$%e7P1HJf^L9k1{^Q9KQ4zwis;cQf&Mv#YZfW7o>)g4=pGe*BUgnZ^a`w-Z
z+n%g*7rZI${~f+$TJ`g3mrX9zuD-XyGVSuirShVIr&nLDo#}1fxAcXyO(yFU`*-09
zD$jbpW_;WBx_<?$_jk7RAI>jtOx!#5n4ahlZo!Ll7J11lFX{HHGp)Y;n={DLPL$`1
zF#E$}N_*T_#<FX_>SWF4a&Pc(vz-3O=WFTexka_5dnLk3re~)a#XOn4bE*_y*tO-Y
z)AcKO=Uuz%*0Shy%CG4wBBbhUFNJfPy!Lp!Yy#`Vmrqu2h}ojAaY~9&k=6Xsn%2mq
z1dFh`Zx3f4b`0()Nc*$yUrSHfk_n=BeWebq%<)|^MR(=7CAnNiFC{O}EY$i{bIKsq
zTD@HD=hXVPuD#Pj&Mi7+W%0(ZSZL~N+s=Fqou97Ezbu~f?^$~JRnpV+>#s`egeyvz
zFU0v*mzb9RdzM%5s_MW(<HlHX7Pq?r8<y7HdH4NA30o1vuhPR8&NJ*i`=HglG55+#
zHkn<D+YTG>922w8%y<ln+YP}-=dAU-e{sDCZ%iL|-m$a;)z1a;rt3#FF4K4){>pvP
zu?_uQ`}SPib!G11#iBKZ%P;-6nRZ8fW$?zQdRGJ%zX^Wr|9oXi)Rf!1OFzwzT<G^i
ztDfbUOV-U#t+Nj9{_4Mcs(bK>JrPA#u{+=GSZ1>{(02R6O+FQoD;<x2e7GcHub7Pe
zx8pNk)s(K`na0@lG|iFUd7Jn8SotrJpE#en?aM#u;qv6nt~a8Ej=z?RU)kF8UDUOB
zb<u&cxiZ>nI+stccXD|neBgH_(^oZlyBb%?O|2?B3zu<u9<nrC?+|LTcfr1ELBHQs
z6)2qjucJGkSHIV#sMPY#`?!1DH^15*>wKEe>?&U<cKPchpL2(wzv)`Gow@SCypQ@{
zHiyeR{nTaa=w$jOXHIe5{OL#LY&_?+c<0n<WzBM&!qMxTHub&=J}_y2g8bf0{(pCG
zn{PT4>yuD(tjeKA^nb29YjDf_78V&}j>Ce1KQbb%LQcd^SaP3x^O1~Cs|%-#X}H}y
z$<pHSvv*}siRiuikHxyL%Lse#@lJojyJEr-&AG=^J8qOD)xG=Cc~N-98=bQI`RkO5
zkER9mIk_-x@YVU4y|i0q=c8Q-jMJl>o<3K#_CFeITM&{SxJD;<jaXv)1YL)vcBY#w
zkF1F6-}^n+LpaY}%G&AkpQsXrrtV{p+qF`jT`O_0Jj}6^N$bfWxqnw3ex|Mnyv}+>
zOL$!Yr}*)oD`jRs-g2y>=D|G~&qtq*U085mM&!+n6Aa9u%c=rqnf7d&{`JJUQg-1l
zcfWm1@2!32U~<D-ujZxF%(#o*%eXggObMBF*XTsRTVAeJe5@>{_a+y9TYs$XmQc0p
z>|&#smz*jW)+;YJT_&p0Fokio56@xUh4Nf`V_M^<>jwK=>G#`soUh>ML+Sp4{)3Yz
zo?r`2n>JxtRa)A*XErr#Q}z}vn<Z;)yk_H@b3D&Q?^?-DmoUG%|BvzlhyBNO3lIAG
z?D6IE?{`uAwE6XQx9GF&3vM~MTz72c_cFN3UKG$}5Z%b#kzmqT9rJ$m&RiFs_#ZJz
z;r<JHMc$~DO*cKYefriVjThv-UY#+py?5>D-HnX@I<((Rd3k*PjQ2k8H&vba(0={F
zlU1p&x1G)wIkfFi-f!d8zgWYZRPXEl*t@3kV*aY#a|)MnmGU0nbko{;caxdI{UrzT
zc$@dOv*pbDy=zWMdb)b)S?RS0mQI%Y{=wS1n#tKAsl0ji&AjbS1>1$!%$@6;qMo-%
zyXW<$Y3Kavj#+_9-<!#4F>H(5cC(&y5x3b|dF9o$hHdBMFKlyidEQ=kRzT+YG`G1Y
zcAP2FUYgXtgs&&7G&W%Q=PDCj9+Tw8CzcEs>Q~0EU6;nYtD?HMvm~s?V)yPIp94nD
z4xD?ReEI5r?RUo7+~s=`Ug)m3%WUzUk*lfL`Hb!SjvW@)-ufNPdBw7z_Weoev-Ljy
z-(H#TnV8t-5oWhT{<iSPxsf56n~q0E*)+2limz3@yvVyx*HXuHUA|Dl^A)ifv%<W3
z9KO9cwOvtq^}IDHQ+8J#;udum=YPu6bH`xM^xLtjitCEFH|QGc_9*Dh^=m8-TWN6a
zG4t_DGrA&f-j=!87;-vblSgo<Q>&t;kG$~mg2nGQxqq3#HM=&_*nMAwYHJzCyP(3C
zNiR2ND!rc=6_HrCeuwSG0^g-`ESEP#`s|rs*1z<^DT$EzTpIF!3tij=ul)}Wy>#cT
zTFbU%qeClAejBa(6Ij-3p)0NS<+Xo*`rP=9Hz#{o9-Qt|_<Dhs@3$!`S--gtp87ja
zFL}|*)j@K{{PutP&{whka<SWGvDK@uWNm%VAjW+|LV;b0tyeoB___L4_6@x2wtY>$
zm~dssf}OYO`?q|U!Te)IV6sHht}HI^?t;xxxomkKL%yfJ`0nJ9{vbTz@k+TQty#PG
z?N5t}{W4S5cH>49yN-S7;@^~C)?a5it9z&Sd$H=5x$hS)KCiI&_T0yhT$fHyIJ)ua
zbHmqv+yuXVcwX%ibad?^kpQm?4vI=KGo=n`I(0Yv{rYI_BvxVV824qxk*B8nwI*|i
zuD<2{k@0!mAq~MqQN_mEzLK{yEibVZ95+2@$Dtm+`>CPwrRbegW{b~j5AQgvu;fIE
zM#k&J+<;7>fSapqCm8Mu-92Mby}hB*POq)v^I5%)yfC-^#Hh8NWkQGKGVR@(UJ)gN
z9n<X+Ii=ory{-*ST(?7Lvexq^rBwm5>wa;j^xpGokN2B(q;=+AZFSq7j|_AQ&;Beh
z)cJX_Mb_f|i+k2lKd#mP-Q9hMqhWD%k>ro|`OKD?TQ9LXrRx9w@$%~VW1H8XoA&Gc
zxqiv?e5aZHeFefIT4HaW%iUkzD4KWf-mbQ4uKliC7+XxajpeWYbbh?;Xpo8dQTe5(
zBpOvJj>?>HocucaLC7x6R}Z^2mbV|+o%Ptm)Vn8&@rdE;n0QmQ_9D+rh4O>1-?3gk
zC3VztS*FlS2SFiz{Rd)^y926g?unNv%(+v(PEld5&P`kS`wp#l0{WEFX1bU}D4hNL
zq{H00sOW>I&hy7_Y()LPHU#$irk=VwX$JGH_?1t?4!Uu_oRy?-lv~QL`{?wo9F^G@
z4lM{-l@PkZD&oh{Q_>S^IUb}jO=>^&R3J!eL6B0+`sq)2lrv3373>mR9!+`tXwK=6
zOSF_~nl>0+E4g}^<Mrwb4l0*3=Q2NUk8rZ6QtnUCDqDKkOV9rI%1vUkpC1vw+_X(B
zq4sIq6VA!$DFts$zB|NwMJB4|sjcSl{$}vz&AX2GRx5m;G1ZHC^;`^ZY`wGSk${#*
z>Y~pNzwpMSMJ#REFy+9-(hZw;t*kKU{8`X+O4t5S;Mwb{M&|Ro&mC!={CcC5_lD&+
zmis@w=E!^Y_*$_ehOcWf)P5cR^xR?p`~2{glTs3OWi}ivGG6xRTC#XszR>2xcDC$?
zdv8_U*vS{|dt=eA-$7o_mKg|qIh_CR5bMMBkN2IDN!^&bY(x2zX9cep$lYwve`2w|
zIm)54WYLUnt=NE`wVWMSjx4AuDh*!bc70EQqyIX;ik$(v7G|M~54!%Gw#a?I#M7#j
z^^33iW?XvnwExQO?qy~!?w3}!vKL=j+9$K6t@zJ;t*c*-Og+VJA{rb~Ft3?kDa7Is
z@7;saKla?HdK$E$-Y;0?sEdZ0M__!9R98!+z4;pU_e)-6)b0_!-XFl=_{Z*q?-A9X
zJBw^Y{f`EJ455^|TL0(ue<sOAHI>Ql!>54`Uhs7Fb6Mw<(1iO;0Z|4924B|@M_o@p
zH~n1L*$frAIcq2T<{efLaQklAdUTu7mdH>IKe35bO<x$U@bC%A3Gf6K{Z7*MiwT{W
zq-kB(_T83k<ITC%;d2$#UT;j4|1&qVVeNIX-EFn<`(F8PyS>hAuT!?d#5ruCmhCUi
zawGaQ_IjClhS;ijuKmIo@Z0QrgUbJFrww14&E{NlQ}nIP9QC|6PgsoX7KeqbTJKoE
zC-8Q|(m7v_-4J`yW1yBBe)TWcF{>5cm3!{6`joBX>+0hb=HL6ymmxRQcRB-au%gB7
zKBq@5v-)ixd|v%oL`h?AAH$I&LR0s;C-0VVuZrqc7Y%ZYvpTT(>DxIQbZVj=2hLV+
z^1d3)89Re<3ENfn8m;&0V$Q;`GE1f>duwM_oSpufDe!0RdT~3M;{VO3cBy1YeQLPd
zEO2`@-;oXRL1vHc_uv1=1X}+#<KDdhLq-M$9Tw2~H%2BA2HfYWfZ_*%bD^iK1bCzB
zLa1P10G*tI;D=y2CIsUg72Jl&AWULlV3>%-Fw}EZaGL}=Sp{L(X)GoopR<D798?!E
zGh!YSQ&5zTbnXgnlTcj*I^_ZVz?ovCb69YjBZ2T10|P@T78ikh208PDSW{MGF(nOt
zP77{RQ2q1}izx+2=eFQB2h~rKOqijJIL8HI3hXG7El7uyfYyy80>__;fdO826s4rW
z5;06C@>v@k2yF}u41LTD47t!FP0)2AZ*&H&2SsRPWJNxe16?2T%6HH<XoS99+$fu(
z(KRBkUPblbLrDe(<Yn?Oow!!FBCHIPV_?9w))i(7ESQi>4^X*`Fy)amj#4_no0Scu
ROprm4A(e%JA<PxT0{|P#;k^I=
literal 0
HcmV?d00001
diff --git a/experiment/wood-bilayer/results/auswertung.py b/experiment/wood-bilayer/results/auswertung.py
new file mode 100644
index 00000000..12e9fc7a
--- /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 00000000..bcaf914a
--- /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 00000000..56a48fd0
--- /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()
+ #
--
GitLab