diff --git a/Plot-Scripts/WoodBilayer_ExperimentComparisonError_localminimizer.py b/Plot-Scripts/WoodBilayer_ExperimentComparisonError_localminimizer.py
new file mode 100644
index 0000000000000000000000000000000000000000..11750cb1e69648d1cb66cd5ced1ee30f3e12ce4a
--- /dev/null
+++ b/Plot-Scripts/WoodBilayer_ExperimentComparisonError_localminimizer.py
@@ -0,0 +1,153 @@
+import numpy as np
+import matplotlib.pyplot as plt
+import math
+import os
+import subprocess
+import fileinput
+import re
+import sys
+import matplotlib as mpl
+from mpl_toolkits.mplot3d import Axes3D
+import matplotlib.cm as cm
+import matplotlib.ticker as ticker
+from matplotlib.ticker import MultipleLocator,FormatStrFormatter,MaxNLocator
+import seaborn as sns
+import matplotlib.colors as mcolors
+# ----------------------------------------------------
+# --- Define Plot style:
+# plt.style.use("seaborn-white")
+# plt.style.use("seaborn-pastel")
+# plt.style.use("seaborn-colorblind")
+plt.style.use("seaborn")
+mpl.rcParams['text.usetex'] = True
+mpl.rcParams["font.family"] = "serif"
+mpl.rcParams["font.size"] = "14"
+mpl.rcParams['xtick.bottom'] = True
+mpl.rcParams['xtick.major.size'] = 3
+mpl.rcParams['xtick.minor.size'] = 1.5
+mpl.rcParams['xtick.major.width'] = 0.75
+mpl.rcParams['ytick.left'] = True
+mpl.rcParams['ytick.major.size'] = 3
+mpl.rcParams['ytick.minor.size'] = 1.5
+mpl.rcParams['ytick.major.width'] = 0.75
+
+#Adjust grid:
+mpl.rcParams.update({"axes.grid" : True}) # Add grid
+mpl.rcParams['axes.labelpad'] = 3
+mpl.rcParams['grid.linewidth'] = 0.25
+mpl.rcParams['grid.alpha'] = 0.9 # 0.75
+mpl.rcParams['grid.linestyle'] = '-'
+mpl.rcParams['grid.color'] = 'gray'#'black'
+mpl.rcParams['text.latex.preamble'] = r'\usepackage{amsfonts}' # Makes Use of \mathbb possible.
+# ----------------------------------------------------------------------------------------
+width = 5.79
+height = width / 1.618 # The golden ratio.
+fig, ax = plt.subplots(figsize=(width,height))
+fig.subplots_adjust(left=.15, bottom=.16, right=.95, top=.92)
+
+# ax.tick_params(axis='x',which='major', direction='out',pad=3)
+
+ax.xaxis.set_major_locator(MultipleLocator(1.0))
+ax.xaxis.set_minor_locator(MultipleLocator(0.5))
+ax.yaxis.set_major_locator(MultipleLocator(0.5))
+
+data = [ # Dataset Ratio r = 0.49
+ [15.30614414,14.49463867,13.46629742,12.78388234,12.23057715,10.21852839,9.341730605], # input (moisture-content in %)
+ [0.60150376, 0.87218045, 1.23308271, 1.5037594 , 1.71428571, 2.46616541,2.79699248], # curvature kappa from Simulation
+ [0.4008016032064128, 0.5911823647294588, 0.8416833667334669, 1.002004008016032, 1.1322645290581161, 1.6132264529058116, 1.813627254509018], # curvature kappa of local minimizerfrom Simulation
+ [0.357615902,0.376287785,0.851008627,0.904475291,1.039744708,1.346405241,1.566568558] # curvature kappa from Experiment
+ ]
+
+# relative_error = (np.array(data[dataset_number][1]) - np.array(dataset[dataset_number][2])) / np.array(dataset[dataset_number][2])
+#print('relative_error:', relative_error)
+
+#--------------- Plot Lines + Scatter -----------------------
+line_1 = ax.plot(np.array(data[0]), np.array(data[1]), # data
+ # color='forestgreen', # linecolor
+ marker='D', # each marker will be rendered as a circle
+ markersize=5, # marker size
+ # markerfacecolor='darkorange', # marker facecolor
+ markeredgecolor='black', # marker edgecolor
+ markeredgewidth=0.75, # marker edge width
+ # linestyle='dashdot', # line style will be dash line
+ linewidth=1.5, # line width
+ zorder=3,
+ label = r"$\kappa_{1,sim}$")
+
+line_2 = ax.plot(np.array(data[0]), np.array(data[2]), # data
+ color='red', # linecolor
+ marker='s', # each marker will be rendered as a circle
+ markersize=5, # marker size
+ # markerfacecolor='cornflowerblue', # marker facecolor
+ markeredgecolor='black', # marker edgecolor
+ markeredgewidth=0.75, # marker edge width
+ # linestyle='--', # line style will be dash line
+ linewidth=1.5, # line width
+ zorder=3,
+ alpha=0.8, # Change opacity
+ label = r"$\kappa_{2,sim}$")
+
+line_3 = ax.plot(np.array(data[0]), np.array(data[3]), # data
+ # color='orangered', # linecolor
+ marker='o', # each marker will be rendered as a circle
+ markersize=5, # marker size
+ # markerfacecolor='cornflowerblue', # marker facecolor
+ markeredgecolor='black', # marker edgecolor
+ markeredgewidth=0.75, # marker edge width
+ # linestyle='--', # line style will be dash line
+ linewidth=1.5, # line width
+ zorder=3,
+ alpha=0.8, # Change opacity
+ label = r"$\kappa_{exp}$")
+
+ # --- Plot order line
+ # x = np.linspace(0.01,1/2,100)
+ # y = CC_L2[0]*x**2
+ # OrderLine = ax.plot(x,y,linestyle='--', label=r"$\mathcal{O}(h)$")
+
+
+
+ # Fix_value = 7.674124
+ # l3 = plt.axhline(y = Fix_value, color = 'black', linewidth=0.75, linestyle = 'dashed')
+ # --------------- Set Axes -----------------------
+ # ax.set_title(r"ratio $r = 0.22$") # Plot - Title
+
+# Plot - Titl
+ax.set_title(r"ratio $r = 0.49$")
+
+# plt.xscale('log') # Use Logarithmic-Scale
+# plt.yscale('log')
+ax.set_xlabel(r"Wood moisture content $\omega (\%)$", labelpad=4)
+ax.set_ylabel(r"curvature $\kappa$", labelpad=4)
+plt.tight_layout()
+
+# # --- Set Line labels
+# line_labels = [r"$CC_{L_2}$",r"$CC_{H_1}$", r"$\mathcal{O}(h)$"]
+
+# --- Set Legend
+legend = ax.legend()
+# legend = fig.legend([line_1 , line_2, OrderLine],
+# labels = line_labels,
+# bbox_to_anchor=[0.97, 0.50],
+# # bbox_to_anchor=[0.97, 0.53],
+# # loc='center',
+# ncol=1, # Number of columns used for legend
+# # borderaxespad=0.15, # Small spacing around legend box
+# frameon=True,
+# prop={'size': 10})
+
+
+frame = legend.get_frame()
+frame.set_edgecolor('black')
+
+
+# --- Adjust left/right spacing:
+# plt.subplots_adjust(right=0.81)
+# plt.subplots_adjust(left=0.11)
+
+# ---------- Output Figure as pdf:
+fig.set_size_inches(width, height)
+fig.savefig('WoodBilayer_expComparison_local_global.pdf')
+plt.show()
+
+
diff --git a/experiment/perforated-bilayer/PolarPlotLocalEnergy.py b/experiment/perforated-bilayer/PolarPlotLocalEnergy.py
index acce25b0ab988909fc5d87b45eae677b5d6b52d6..2c17e8e02c1016c7cf5839e4cf833c9c03e0c302 100644
--- a/experiment/perforated-bilayer/PolarPlotLocalEnergy.py
+++ b/experiment/perforated-bilayer/PolarPlotLocalEnergy.py
@@ -57,7 +57,6 @@ show_plot = False
perforatedLayer = 'lower'
dataset_numbers = [0, 1, 2, 3, 4, 5]
-# dataset_numbers = [0]
for dataset_number in dataset_numbers:
@@ -110,7 +109,7 @@ for dataset_number in dataset_numbers:
width = 5.79
height = width / 1.618 # The golden ratio.
fig.set_size_inches(width, height)
- fig.savefig('PerforatedBilayer_dataset_' +str(dataset_number) + '_exp' +str(n) + '.pdf')
+ fig.savefig('./experiment/perforated-bilayer/PerforatedBilayer_dataset_' +str(dataset_number) + '_exp' +str(n) + '.pdf')
f = open("./experiment/perforated-bilayer/results_" + perforatedLayer + '_' + str(dataset_number) + "/kappa_simulation.txt", "w")
diff --git a/experiment/perforated-bilayer/perfBilayer_test.py b/experiment/perforated-bilayer/perfBilayer_test.py
index 01d1c57ea4a6373e6d37ae7ee1bd48a2424c47a6..2a4417200032f97804c94b4c7eaf40af6a294078 100644
--- a/experiment/perforated-bilayer/perfBilayer_test.py
+++ b/experiment/perforated-bilayer/perfBilayer_test.py
@@ -129,7 +129,7 @@ for dataset_number in dataset_numbers:
gamma = 1.0
# ----- Define Parameters for Material Function --------------------
- # [r, h, omega_flat, omega_target, theta, experimental_kappa]
+ # [r, h, omega_flat, omega_target, theta, beta]
# r = (thickness upper layer)/(thickness)
# h = thickness [meter]
# omega_flat = moisture content in the flat state before drying [%]
diff --git a/experiment/wood-bilayer_PLOS/Auswertung.py b/experiment/wood-bilayer_PLOS/Auswertung.py
new file mode 100644
index 0000000000000000000000000000000000000000..9bd21562b7099b48f9b71e9fc1aa17bf0d884f96
--- /dev/null
+++ b/experiment/wood-bilayer_PLOS/Auswertung.py
@@ -0,0 +1,136 @@
+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+"""
+Created on Wed Jul 6 13:17:28 2022
+
+@author: stefan
+"""
+import numpy as np
+import matplotlib.pyplot as plt
+import matplotlib.colors as colors
+import codecs
+import re
+import json
+
+def energy(kappa,alpha,Q,B) :
+ G=kappa*np.array([[np.cos(alpha)**2],[np.sin(alpha)**2],[np.sqrt(2)*np.cos(alpha)*np.sin(alpha)]])-B
+ return np.matmul(np.transpose(G),np.matmul(Q,G))[0,0]
+
+def xytokappaalpha(x,y):
+
+ if y>0:
+ return [np.sqrt(x**2+y**2), np.abs(np.arctan2(y,x))]
+ else:
+ return [-np.sqrt(x**2+y**2), np.abs(np.arctan2(y,x))]
+
+# Read effective quantites
+def ReadEffectiveQuantities(QFilePath, BFilePath):
+ # Read Output Matrices (effective quantities)
+ # From Cell-Problem output Files : ../outputs/Qmatrix.txt , ../outputs/Bmatrix.txt
+ # -- Read Matrix Qhom
+ X = []
+ # with codecs.open(path + '/outputs/QMatrix.txt', encoding='utf-8-sig') as f:
+ with codecs.open(QFilePath, encoding='utf-8-sig') as f:
+ for line in f:
+ s = line.split()
+ X.append([float(s[i]) for i in range(len(s))])
+ Q = np.array([[X[0][2], X[1][2], X[2][2]],
+ [X[3][2], X[4][2], X[5][2]],
+ [X[6][2], X[7][2], X[8][2]] ])
+
+ # -- Read Beff (as Vector)
+ X = []
+ # with codecs.open(path + '/outputs/BMatrix.txt', encoding='utf-8-sig') as f:
+ with codecs.open(BFilePath, encoding='utf-8-sig') as f:
+ for line in f:
+ s = line.split()
+ X.append([float(s[i]) for i in range(len(s))])
+ B = np.array([X[0][2], X[1][2], X[2][2]])
+ return Q, B
+
+# Number of experiments / folders
+number=7
+show_plot = False
+save_plot = False
+
+#--- Choose wether to perforate upper (passive) or lower (active) layer
+# perforatedLayer = 'upper'
+perforatedLayer = 'lower'
+
+dataset_indices=[[0,6],[1,6],[2,6],[3,6],[4,6],[5,6]]
+#dataset_indices=[[5,0]]
+for dataset_index in dataset_indices:
+ dataset_number=dataset_index[0]
+ n=dataset_index[1]
+ kappa=np.zeros(number)
+ kappa_pos=np.zeros(number)
+ kappa_neg=np.zeros(number)
+ # Read from Date
+ print(str(n))
+ DataPath = './experiment/wood-bilayer_PLOS/results_' + perforatedLayer + '_' + str(dataset_number) + '/' +str(n)
+ QFilePath = DataPath + '/QMatrix.txt'
+ BFilePath = DataPath + '/BMatrix.txt'
+ ParameterPath = DataPath + '/parameter.txt'
+
+ # Read Thickness from parameter file (needed for energy scaling)
+ with open(ParameterPath , 'r') as file:
+ parameterFile = file.read()
+ thickness = float(re.findall(r'(?m)h = (\d?\d?\d?\.?\d+[Ee]?[+\-]?\d?\d?)',parameterFile)[0])
+
+ Q, B = ReadEffectiveQuantities(QFilePath,BFilePath)
+ # Q=0.5*(np.transpose(Q)+Q) # symmetrize
+ B=np.transpose([B])
+ #
+
+ N=500
+ length=5
+ r, theta = np.meshgrid(np.linspace(0,length,N),np.radians(np.linspace(0, 360, N)))
+ E=np.zeros(np.shape(r))
+ for i in range(0,N):
+ for j in range(0,N):
+ if theta[i,j]<np.pi:
+ E[i,j]=energy(r[i,j],theta[i,j],Q,B) * (thickness**2)
+ else:
+ E[i,j]=energy(-r[i,j],theta[i,j],Q,B) * (thickness**2)
+
+ # Compute Minimizer
+ [imin,jmin]=np.unravel_index(E.argmin(),(N,N))
+ kappamin=r[imin,jmin]
+ alphamin=theta[imin,jmin]
+ Emin=E[imin,jmin]
+ kappa[n]=kappamin
+ # Positiv curvature region
+ N_mid=int(N/2)
+ [imin,jmin]=np.unravel_index(E[:N_mid,:].argmin(),(N_mid,N))
+ kappamin_pos=r[imin,jmin]
+ alphamin_pos=theta[imin,jmin]
+ Emin_pos=E[imin,jmin]
+ kappa_pos[n]=kappamin_pos
+ # Negative curvature region
+ [imin,jmin]=np.unravel_index(E[N_mid:,:].argmin(),(N_mid,N))
+ kappamin_neg=r[imin+N_mid,jmin]
+ alphamin_neg=theta[imin+N_mid,jmin]
+ Emin_neg=E[imin+N_mid,jmin]
+ kappa_neg[n]=kappamin_neg
+ #
+ fig, ax = plt.subplots(figsize=(6,6),subplot_kw=dict(projection='polar'))
+ levs=np.geomspace(E.min(),E.max(),400)
+ pcm=ax.contourf(theta, r, E, levs, norm=colors.PowerNorm(gamma=0.2), cmap='brg')
+ ax.set_xticks([0,np.pi/2])
+ ax.set_yticks([kappamin])
+ colorbarticks=np.linspace(E.min(),E.max(),6)
+ cbar = plt.colorbar(pcm, extend='max', ticks=colorbarticks, pad=0.1)
+ cbar.ax.tick_params(labelsize=6)
+ # We compute the relative energy gap between the minimal energy and the energy of the flat state. (E(0)-Emin)/Emin
+ energygap=(energy(0,0,Q,B)*(thickness**2)-E.min())/(energy(0,0,Q,B)*(thickness**2))
+ print("Minimum global, positiv, negativ: ", Emin, Emin_pos, Emin_neg)
+ print("Kappa global, positiv, negativ: ", kappamin, kappamin_pos, kappamin_neg)
+ print("Energy gap:", energygap)
+ if (show_plot):
+ plt.show()
+ if (save_plot):
+ # Save Figure as .pdf
+ width = 5.79
+ height = width / 1.618 # The golden ratio.
+ fig.set_size_inches(width, height)
+ fig.savefig('./experiment/wood-bilayer_PLOS/wood-bilayer_PLOS_dataset_' +str(dataset_number) + '_exp_auswertung' +str(n) + '.pdf')
diff --git a/experiment/wood-bilayer_PLOS/PolarPlotLocalEnergy.py b/experiment/wood-bilayer_PLOS/PolarPlotLocalEnergy.py
new file mode 100644
index 0000000000000000000000000000000000000000..7c3feaac02477d6f6ee29abfb36937587aa94aa6
--- /dev/null
+++ b/experiment/wood-bilayer_PLOS/PolarPlotLocalEnergy.py
@@ -0,0 +1,131 @@
+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+"""
+Created on Wed Jul 6 13:17:28 2022
+
+@author: stefan
+"""
+import numpy as np
+import matplotlib.pyplot as plt
+import matplotlib.colors as colors
+import codecs
+import re
+import json
+
+def energy(kappa,alpha,Q,B) :
+ G=kappa*np.array([[np.cos(alpha)**2],[np.sin(alpha)**2],[np.sqrt(2)*np.cos(alpha)*np.sin(alpha)]])-B
+ return np.matmul(np.transpose(G),np.matmul(Q,G))[0,0]
+
+def xytokappaalpha(x,y):
+
+ if y>0:
+ return [np.sqrt(x**2+y**2), np.abs(np.arctan2(y,x))]
+ else:
+ return [-np.sqrt(x**2+y**2), np.abs(np.arctan2(y,x))]
+
+# Read effective quantites
+def ReadEffectiveQuantities(QFilePath, BFilePath):
+ # Read Output Matrices (effective quantities)
+ # From Cell-Problem output Files : ../outputs/Qmatrix.txt , ../outputs/Bmatrix.txt
+ # -- Read Matrix Qhom
+ X = []
+ # with codecs.open(path + '/outputs/QMatrix.txt', encoding='utf-8-sig') as f:
+ with codecs.open(QFilePath, encoding='utf-8-sig') as f:
+ for line in f:
+ s = line.split()
+ X.append([float(s[i]) for i in range(len(s))])
+ Q = np.array([[X[0][2], X[1][2], X[2][2]],
+ [X[3][2], X[4][2], X[5][2]],
+ [X[6][2], X[7][2], X[8][2]] ])
+
+ # -- Read Beff (as Vector)
+ X = []
+ # with codecs.open(path + '/outputs/BMatrix.txt', encoding='utf-8-sig') as f:
+ with codecs.open(BFilePath, encoding='utf-8-sig') as f:
+ for line in f:
+ s = line.split()
+ X.append([float(s[i]) for i in range(len(s))])
+ B = np.array([X[0][2], X[1][2], X[2][2]])
+ return Q, B
+
+# Number of experiments / folders
+show_plot = False
+
+#--- Choose wether to perforate upper (passive) or lower (active) layer
+
+dataset_numbers = [0, 1, 2, 3, 4, 5]
+number=7
+for dataset_number in dataset_numbers:
+ kappa=np.zeros(number)
+ kappa_pos=np.zeros(number)
+ kappa_neg=np.zeros(number)
+ for n in range(0,number):
+ # Read from Date
+ print(str(n))
+ DataPath = './experiment/wood-bilayer_PLOS/results_' + str(dataset_number) + '/' +str(n)
+ QFilePath = DataPath + '/QMatrix.txt'
+ BFilePath = DataPath + '/BMatrix.txt'
+ ParameterPath = DataPath + '/parameter.txt'
+
+ # Read Thickness from parameter file (needed for energy scaling)
+ with open(ParameterPath , 'r') as file:
+ parameterFile = file.read()
+ thickness = float(re.findall(r'(?m)h = (\d?\d?\d?\.?\d+[Ee]?[+\-]?\d?\d?)',parameterFile)[0])
+
+ Q, B = ReadEffectiveQuantities(QFilePath,BFilePath)
+ # Q=0.5*(np.transpose(Q)+Q) # symmetrize
+ B=np.transpose([B])
+ #
+
+ N=500
+ length=5
+ r, theta = np.meshgrid(np.linspace(0,length,N),np.radians(np.linspace(0, 360, N)))
+ E=np.zeros(np.shape(r))
+ for i in range(0,N):
+ for j in range(0,N):
+ if theta[i,j]<np.pi:
+ E[i,j]=energy(r[i,j],theta[i,j],Q,B) * (thickness**2)
+ else:
+ E[i,j]=energy(-r[i,j],theta[i,j],Q,B) * (thickness**2)
+
+ # Compute Minimizer
+ [imin,jmin]=np.unravel_index(E.argmin(),(N,N))
+ kappamin=r[imin,jmin]
+ alphamin=theta[imin,jmin]
+ kappa[n]=kappamin
+ # Positiv curvature region
+ N_mid=int(N/2)
+ [imin,jmin]=np.unravel_index(E[:N_mid,:].argmin(),(N_mid,N))
+ kappamin_pos=r[imin,jmin]
+ alphamin_pos=theta[imin,jmin]
+ Emin_pos=E[imin,jmin]
+ kappa_pos[n]=kappamin_pos
+ # Negative curvature region
+ [imin,jmin]=np.unravel_index(E[N_mid:,:].argmin(),(N_mid,N))
+ kappamin_neg=r[imin+N_mid,jmin]
+ alphamin_neg=theta[imin+N_mid,jmin]
+ Emin_neg=E[imin+N_mid,jmin]
+ kappa_neg[n]=kappamin_neg
+ #
+ fig, ax = plt.subplots(figsize=(6,6),subplot_kw=dict(projection='polar'))
+ levs=np.geomspace(E.min(),E.max(),400)
+ pcm=ax.contourf(theta, r, E, levs, norm=colors.PowerNorm(gamma=0.2), cmap='brg')
+ ax.set_xticks([0,np.pi/2])
+ ax.set_yticks([kappamin])
+ colorbarticks=np.linspace(E.min(),E.max(),6)
+ cbar = plt.colorbar(pcm, extend='max', ticks=colorbarticks, pad=0.1)
+ cbar.ax.tick_params(labelsize=6)
+ if (show_plot):
+ plt.show()
+ # Save Figure as .pdf
+ width = 5.79
+ height = width / 1.618 # The golden ratio.
+ fig.set_size_inches(width, height)
+ fig.savefig('./experiment/wood-bilayer_PLOS/wood-bilayer_PLOS_dataset_' +str(dataset_number) + '_exp' +str(n) + '.pdf')
+
+
+ f = open("./experiment/wood-bilayer_PLOS/results_" + str(dataset_number) + "/kappa_simulation.txt", "w")
+ f.write(str(kappa.tolist())[1:-1]+"\n")
+ f.write(str(kappa_pos.tolist())[1:-1]+"\n")
+ f.write(str(kappa_neg.tolist())[1:-1]+"\n")
+ f.close()
\ No newline at end of file
diff --git a/experiment/wood-bilayer_PLOS/perforated_wood_lower.py b/experiment/wood-bilayer_PLOS/perforated_wood_lower.py
new file mode 100644
index 0000000000000000000000000000000000000000..3627cf878b980eec37461757cfa50b5dd6cf7c6c
--- /dev/null
+++ b/experiment/wood-bilayer_PLOS/perforated_wood_lower.py
@@ -0,0 +1,283 @@
+import math
+#from python_matrix_operations import *
+import ctypes
+import os
+import sys
+import numpy as np
+# import elasticity_toolbox as elast
+
+class ParameterSet(dict):
+ def __init__(self, *args, **kwargs):
+ super(ParameterSet, self).__init__(*args, **kwargs)
+ self.__dict__ = self
+
+parameterSet = ParameterSet()
+#---------------------------------------------------------------
+#############################################
+# Paths
+#############################################
+# Path for results and logfile
+parameterSet.outputPath='/home/klaus/Desktop/Dune_release/dune-microstructure/experiment/perforated-bilayer/results'
+parameterSet.baseName= 'perforated_wood_lower' #(needed for Output-Filename)
+
+# Path for material description
+# parameterSet.geometryFunctionPath =experiment/wood-bilayer/
+
+#---------------------------------------------------------------
+# Wooden bilayer, https://www.ncbi.nlm.nih.gov/pmc/articles/PMC6191116/#pone.0205607.ref015
+#--- define indicator function
+# x[0] : y1-component -1/2 to 1/2
+# x[1] : y2-component -1/2 to 1/2
+# x[2] : x3-component range -1/2 to 1/2
+#--- define indicator function
+def indicatorFunction(x):
+ pRadius = math.sqrt((param_beta*param_r)/(np.pi*perfDepth)) # perforation radius
+ if (x[2]>=(0.5-param_r)):
+ # if(np.sqrt(x[0]**2 + x[1]**2) < pRadius): #inside perforation
+ return 1 #Phase1
+ else :
+ if(((x[0]**2 + x[1]**2) < pRadius**2) and (x[2] <= (-0.5+perfDepth))): #inside perforation
+ return 3 #Phase3
+ else:
+ return 2 #Phase2
+
+# def indicatorFunction(x):
+# factor=1
+# pRadius = 0.25
+# if (x[2]>=(0.5-param_r) and np.sqrt(x[0]**2 + x[1]**2) < pRadius):
+# return 3
+# elif((x[2]>=(0.5-param_r))):
+# return 1 #Phase1
+# else :
+# return 2 #Phase2
+
+# # --- Number of material phases
+# parameterSet.Phases=3
+
+# def indicatorFunction(x):
+# factor=1
+# pRadius = 1
+# # if (np.sqrt(x[0]*x[0] + x[1]*x[1]) < pRadius):
+# if ((x[0] < 0 and math.sqrt(pow(x[1],2) + pow(x[2],2) ) < pRadius/4.0) or ( 0 < x[0] and math.sqrt(pow(x[1],2) + pow(x[2],2) ) > pRadius/4.0)):
+# return 1
+# else :
+# return 2 #Phase2
+
+# --- Number of material phases
+parameterSet.Phases=3
+
+
+# Parameters of the model
+# -- (thickness upper layer) / (thickness)
+# param_r = 0.22
+param_r = 0.49
+# -- thickness [meter]
+param_h = 0.008
+# -- moisture content in the flat state [%]
+param_omega_flat = 17.01520754
+# -- moisture content in the target state [%]
+param_omega_target = 9.341730605
+# -- Drehwinkel
+param_theta = 0.0
+
+# Design Parameter ratio between perforaton (cylindrical) volume and volume of upper layer
+param_beta = 0
+
+# Depth of perforation
+# perfDepth = 0.12
+perfDepth = (1.0-param_r)
+# perfDepth = (1.0-param_r) * (2.0/3.0)
+
+#
+#
+#
+# -- increment of the moisture content
+delta_omega=param_omega_target-param_omega_flat
+# moisture content for material law
+omega=param_omega_target
+
+# --- Material properties from https://www.ncbi.nlm.nih.gov/pmc/articles/PMC6191116/#pone.0205607.ref015
+# --- for European beech, moisture content omega = 15%
+# --- L=direction orthogonal to layering and fibres = orthogonal to wood stem cross-section
+# --- T=tangential zu layering
+# --- R=orthogonal zu layering
+# --- in MPa
+# --- Properties are defined by affine function in dependence of moisture content omega via property = b_0+b_1 \omega
+# --- coefficients of affine function are contained in the following array
+# --- data taken from http://dx.doi.org/10.1016/j.cma.2014.10.031
+
+properties_coefficients=np.array([
+ # [b_0, b_1]
+ [2565.6,-59.7], # E_R [MPa]
+ [885.4, -23.4], # E_T [MPa]
+ [17136.7,-282.4], # E_L [MPa]
+ [667.8, -15.19], # G_RT [MPa]
+ [1482, -15.26], # G_RL [MPa]
+ [1100, -17.72], # G_TL [MPa]
+ [0.2933, -0.001012], # nu_TR [1]
+ [0.383, -0.008722], # nu_LR [1]
+ [0.3368, -0.009071] # nu_LT [1]
+ ])
+# Compute actual material properties
+E_R = properties_coefficients[0,0]+properties_coefficients[0,1]*omega
+E_T = properties_coefficients[1,0]+properties_coefficients[1,1]*omega
+E_L = properties_coefficients[2,0]+properties_coefficients[2,1]*omega
+G_RT = properties_coefficients[3,0]+properties_coefficients[3,1]*omega
+G_LR = properties_coefficients[4,0]+properties_coefficients[4,1]*omega
+G_LT = properties_coefficients[5,0]+properties_coefficients[5,1]*omega
+nu_TR = properties_coefficients[6,0]+properties_coefficients[6,1]*omega
+nu_LR = properties_coefficients[7,0]+properties_coefficients[7,1]*omega
+nu_LT = properties_coefficients[8,0]+properties_coefficients[8,1]*omega
+# Compute the remaining Poisson ratios
+nu_TL=nu_LT*E_T/E_L
+nu_RT=nu_TR*E_R/E_T
+nu_RL=nu_LR*E_R/E_L
+#
+# --- differential swelling strain
+# --- relation to swelling strain eps: eps=alpha* delta_omega with delta_omega = change of water content in %
+alpha_L=0.00011 # PLOS paper
+alpha_R=0.00191 # PLOS paper
+alpha_T=0.00462 # PLOS paper
+# Umrechnen
+#alpha_L=(1-1/(1+delta_omega*alpha_L))/delta_omega
+#alpha_R=(1-1/(1+delta_omega*alpha_R))/delta_omega
+#alpha_T=(1-1/(1+delta_omega*alpha_T))/delta_omega
+# --- define geometry
+
+
+
+# --- PHASE 1
+# y_1-direction: L
+# y_2-direction: T
+# x_3-direction: R
+# phase1_type="orthotropic"
+# materialParameters_phase1 = [E_L,E_T,E_R,G_TL,G_RT,G_RL,nu_LT,nu_LR,nu_TR]
+parameterSet.phase1_type="general_anisotropic"
+[E_1,E_2,E_3]=[E_L,E_T,E_R]
+[nu_12,nu_13,nu_23]=[nu_LT,nu_LR,nu_TR]
+[nu_21,nu_31,nu_32]=[nu_TL,nu_RL,nu_RT]
+[G_12,G_31,G_23]=[G_LT,G_LR,G_RT]
+compliance_S=np.array([[1/E_1, -nu_21/E_2, -nu_31/E_3, 0.0, 0.0, 0.0],
+ [-nu_12/E_1, 1/E_2, -nu_32/E_3, 0.0, 0.0, 0.0],
+ [-nu_13/E_1, -nu_23/E_2, 1/E_3, 0.0, 0.0, 0.0],
+ [0.0, 0.0, 0.0, 1/G_23, 0.0, 0.0],
+ [0.0, 0.0, 0.0, 0.0, 1/G_31, 0.0],
+ [0.0, 0.0, 0.0, 0.0, 0.0, 1/G_12]]);
+materialParameters_phase1 = compliance_S
+
+def prestrain_phase1(x):
+ # hB=delta_omega * alpha with delta_omega increment of moisture content and alpha swelling factor.
+ return [[1/param_h*delta_omega*alpha_L, 0, 0], [0,1/param_h*delta_omega*alpha_T,0], [0,0,1/param_h*delta_omega*alpha_R]]
+
+# --- PHASE 2
+# y_1-direction: R
+# y_2-direction: L
+# x_3-direction: T
+parameterSet.phase2_type="general_anisotropic"
+[E_1,E_2,E_3]=[E_R,E_L,E_T]
+[nu_12,nu_13,nu_23]=[nu_RL,nu_RT,nu_LT]
+[nu_21,nu_31,nu_32]=[nu_LR,nu_TR,nu_TL]
+[G_12,G_31,G_23]=[G_LR,G_RT,G_LT]
+compliance_S=np.array([[1/E_1, -nu_21/E_2, -nu_31/E_3, 0.0, 0.0, 0.0],
+ [-nu_12/E_1, 1/E_2, -nu_32/E_3, 0.0, 0.0, 0.0],
+ [-nu_13/E_1, -nu_23/E_2, 1/E_3, 0.0, 0.0, 0.0],
+ [0.0, 0.0, 0.0, 1/G_23, 0.0, 0.0],
+ [0.0, 0.0, 0.0, 0.0, 1/G_31, 0.0],
+ [0.0, 0.0, 0.0, 0.0, 0.0, 1/G_12]]);
+materialParameters_phase2 = compliance_S
+def prestrain_phase2(x):
+ return [[1/param_h*delta_omega*alpha_R, 0, 0], [0,1/param_h*delta_omega*alpha_L,0], [0,0,1/param_h*delta_omega*alpha_T]]
+
+#Rotation um 2. Achse (= L)
+parameterSet.phase2_axis = 1
+# phase2_angle = param_theta
+# -- Drehwinkel
+parameterSet.phase2_angle = param_theta
+
+
+# --- PHASE 3
+parameterSet.phase3_type="isotropic"
+epsilon = 1e-8
+materialParameters_phase3 = [epsilon, epsilon]
+
+def prestrain_phase3(x):
+ return [[0, 0, 0], [0,0,0], [0,0,0]]
+
+# # --- PHASE 3 = Phase 1 gedreht
+# # y_1-direction: L
+# # y_2-direction: R
+# # x_3-direction: T
+# parameterSet.phase3_type="general_anisotropic"
+# # Drehung um theta um Achse 2 = x_3-Achse
+# N=elast.rotation_matrix_compliance(2,param_theta)
+# materialParameters_phase3 = np.dot(np.dot(N,materialParameters_phase1),N.T)
+# materialParameters_phase3 = 0.5*(materialParameters_phase3.T+materialParameters_phase3)
+# # rotation of strain
+# def prestrain_phase3(x):
+# return elast.voigt_to_strain(np.dot(elast.rotation_matrix_compliance(2,param_theta),np.dot(elast.strain_to_voigt(np.array(prestrain_phase1(x))),N.T))).tolist()
+
+
+
+# --- Choose scale ratio gamma:
+parameterSet.gamma=1.0
+
+
+
+
+#############################################
+# Grid parameters
+#############################################
+## numLevels : Number of Levels on which solution is computed. starting with a 2x2x2 cube mesh.
+## {start,finish} computes on all grid from 2^(start) to 2^finish refinement
+#----------------------------------------------------
+parameterSet.numLevels= '3 3' # computes all levels from first to second entry
+# parameterSet.numLevels= '4 4'
+
+#############################################
+# Assembly options
+#############################################
+parameterSet.set_IntegralZero = 1 #(default = false)
+parameterSet.set_oneBasisFunction_Zero = 1 #(default = false)
+#parameterSet.arbitraryLocalIndex = 7 #(default = 0)
+#parameterSet.arbitraryElementNumber = 3 #(default = 0)
+
+#############################################
+# Solver Options, Type: #1: CG - SOLVER , #2: GMRES - SOLVER, #3: QR - SOLVER (default), #4: UMFPACK - SOLVER
+#############################################
+parameterSet.Solvertype = 3 # recommended to use iterative solver (e.g GMRES) for finer grid-levels
+parameterSet.Solver_verbosity = 0 #(default = 2) degree of information for solver output
+
+
+#############################################
+# Write/Output options #(default=false)
+#############################################
+# --- (Optional output) write Material / prestrain / Corrector functions to .vtk-Files:
+parameterSet.write_materialFunctions = 1 # VTK indicator function for material/prestrain definition
+#parameterSet.write_prestrainFunctions = 1 # VTK norm of B (currently not implemented)
+
+# --- (Additional debug output)
+parameterSet.print_debug = 0 #(default=false)
+
+# --- Write Correctos to VTK-File:
+parameterSet.write_VTK = 1
+
+# The grid can be refined several times for a higher resolution in the VTK-file.
+parameterSet.subsamplingRefinement = 2
+
+# --- (Optional output) L2Error, integral mean:
+#parameterSet.write_L2Error = 1
+#parameterSet.write_IntegralMean = 1
+
+# --- check orthogonality (75) from paper:
+parameterSet.write_checkOrthogonality = 0
+
+# --- Write corrector-coefficients to log-File:
+#parameterSet.write_corrector_phi1 = 1
+#parameterSet.write_corrector_phi2 = 1
+#parameterSet.write_corrector_phi3 = 1
+
+# --- Print Condition number of matrix (can be expensive):
+#parameterSet.print_conditionNumber= 1 #(default=false)
+
+# --- write effective quantities to Matlab-folder for symbolic minimization:
+parameterSet.write_toMATLAB = 1 # writes effective quantities to .txt-files QMatrix.txt and BMatrix.txt
diff --git a/experiment/wood-bilayer_PLOS/perforated_wood_upper.py b/experiment/wood-bilayer_PLOS/perforated_wood_upper.py
new file mode 100644
index 0000000000000000000000000000000000000000..06be5119e77f13d37f0950bce810822c59fded37
--- /dev/null
+++ b/experiment/wood-bilayer_PLOS/perforated_wood_upper.py
@@ -0,0 +1,280 @@
+import math
+#from python_matrix_operations import *
+import ctypes
+import os
+import sys
+import numpy as np
+# import elasticity_toolbox as elast
+
+class ParameterSet(dict):
+ def __init__(self, *args, **kwargs):
+ super(ParameterSet, self).__init__(*args, **kwargs)
+ self.__dict__ = self
+
+parameterSet = ParameterSet()
+#---------------------------------------------------------------
+#############################################
+# Paths
+#############################################
+# Path for results and logfile
+parameterSet.outputPath='/home/klaus/Desktop/Dune_release/dune-microstructure/experiment/perforated-bilayer/results'
+parameterSet.baseName= 'perforated_wood_upper' #(needed for Output-Filename)
+
+# Path for material description
+# parameterSet.geometryFunctionPath =experiment/wood-bilayer/
+
+#---------------------------------------------------------------
+# Wooden bilayer, https://www.ncbi.nlm.nih.gov/pmc/articles/PMC6191116/#pone.0205607.ref015
+#--- define indicator function
+# x[0] : y1-component -1/2 to 1/2
+# x[1] : y2-component -1/2 to 1/2
+# x[2] : x3-component range -1/2 to 1/2
+#--- define indicator function
+def indicatorFunction(x):
+ pRadius = math.sqrt((param_beta*param_r)/(np.pi*perfDepth)) # perforation radius
+ if (x[2]>=(0.5-param_r)):
+ if(((x[0]**2 + x[1]**2) < pRadius**2) and (x[2] >= (0.5-perfDepth))): #inside perforation
+ return 3 #Phase3
+ else:
+ return 1 #Phase1
+ else :
+ return 2 #Phase2
+
+# def indicatorFunction(x):
+# factor=1
+# pRadius = 0.25
+# if (x[2]>=(0.5-param_r) and np.sqrt(x[0]**2 + x[1]**2) < pRadius):
+# return 3
+# elif((x[2]>=(0.5-param_r))):
+# return 1 #Phase1
+# else :
+# return 2 #Phase2
+
+# # --- Number of material phases
+# parameterSet.Phases=3
+
+# def indicatorFunction(x):
+# factor=1
+# pRadius = 1
+# # if (np.sqrt(x[0]*x[0] + x[1]*x[1]) < pRadius):
+# if ((x[0] < 0 and math.sqrt(pow(x[1],2) + pow(x[2],2) ) < pRadius/4.0) or ( 0 < x[0] and math.sqrt(pow(x[1],2) + pow(x[2],2) ) > pRadius/4.0)):
+# return 1
+# else :
+# return 2 #Phase2
+
+# --- Number of material phases
+parameterSet.Phases=3
+
+
+# Parameters of the model
+# -- (thickness upper layer) / (thickness)
+# param_r = 0.22
+param_r = 0.49
+# -- thickness [meter]
+param_h = 0.008
+# -- moisture content in the flat state [%]
+param_omega_flat = 17.17547062
+# -- moisture content in the target state [%]
+param_omega_target = 8.959564147
+# -- Drehwinkel
+param_theta = 0.0
+
+# Design Parameter ratio between perforaton (cylindrical) volume and volume of upper layer
+param_beta = 0.3
+# Depth of perforation
+# perfDepth = 0.12
+perfDepth = param_r
+# perfDepth = param_r * (2.0/3.0)
+#
+#
+#
+# -- increment of the moisture content
+delta_omega=param_omega_target-param_omega_flat
+# moisture content for material law
+omega=param_omega_target
+
+# --- Material properties from https://www.ncbi.nlm.nih.gov/pmc/articles/PMC6191116/#pone.0205607.ref015
+# --- for European beech, moisture content omega = 15%
+# --- L=direction orthogonal to layering and fibres = orthogonal to wood stem cross-section
+# --- T=tangential zu layering
+# --- R=orthogonal zu layering
+# --- in MPa
+# --- Properties are defined by affine function in dependence of moisture content omega via property = b_0+b_1 \omega
+# --- coefficients of affine function are contained in the following array
+# --- data taken from http://dx.doi.org/10.1016/j.cma.2014.10.031
+
+properties_coefficients=np.array([
+ # [b_0, b_1]
+ [2565.6,-59.7], # E_R [MPa]
+ [885.4, -23.4], # E_T [MPa]
+ [17136.7,-282.4], # E_L [MPa]
+ [667.8, -15.19], # G_RT [MPa]
+ [1482, -15.26], # G_RL [MPa]
+ [1100, -17.72], # G_TL [MPa]
+ [0.2933, -0.001012], # nu_TR [1]
+ [0.383, -0.008722], # nu_LR [1]
+ [0.3368, -0.009071] # nu_LT [1]
+ ])
+# Compute actual material properties
+E_R = properties_coefficients[0,0]+properties_coefficients[0,1]*omega
+E_T = properties_coefficients[1,0]+properties_coefficients[1,1]*omega
+E_L = properties_coefficients[2,0]+properties_coefficients[2,1]*omega
+G_RT = properties_coefficients[3,0]+properties_coefficients[3,1]*omega
+G_LR = properties_coefficients[4,0]+properties_coefficients[4,1]*omega
+G_LT = properties_coefficients[5,0]+properties_coefficients[5,1]*omega
+nu_TR = properties_coefficients[6,0]+properties_coefficients[6,1]*omega
+nu_LR = properties_coefficients[7,0]+properties_coefficients[7,1]*omega
+nu_LT = properties_coefficients[8,0]+properties_coefficients[8,1]*omega
+# Compute the remaining Poisson ratios
+nu_TL=nu_LT*E_T/E_L
+nu_RT=nu_TR*E_R/E_T
+nu_RL=nu_LR*E_R/E_L
+#
+# --- differential swelling strain
+# --- relation to swelling strain eps: eps=alpha* delta_omega with delta_omega = change of water content in %
+alpha_L=0.00011 # PLOS paper
+alpha_R=0.00191 # PLOS paper
+alpha_T=0.00462 # PLOS paper
+# Umrechnen
+#alpha_L=(1-1/(1+delta_omega*alpha_L))/delta_omega
+#alpha_R=(1-1/(1+delta_omega*alpha_R))/delta_omega
+#alpha_T=(1-1/(1+delta_omega*alpha_T))/delta_omega
+# --- define geometry
+
+
+
+# --- PHASE 1
+# y_1-direction: L
+# y_2-direction: T
+# x_3-direction: R
+# phase1_type="orthotropic"
+# materialParameters_phase1 = [E_L,E_T,E_R,G_TL,G_RT,G_RL,nu_LT,nu_LR,nu_TR]
+parameterSet.phase1_type="general_anisotropic"
+[E_1,E_2,E_3]=[E_L,E_T,E_R]
+[nu_12,nu_13,nu_23]=[nu_LT,nu_LR,nu_TR]
+[nu_21,nu_31,nu_32]=[nu_TL,nu_RL,nu_RT]
+[G_12,G_31,G_23]=[G_LT,G_LR,G_RT]
+compliance_S=np.array([[1/E_1, -nu_21/E_2, -nu_31/E_3, 0.0, 0.0, 0.0],
+ [-nu_12/E_1, 1/E_2, -nu_32/E_3, 0.0, 0.0, 0.0],
+ [-nu_13/E_1, -nu_23/E_2, 1/E_3, 0.0, 0.0, 0.0],
+ [0.0, 0.0, 0.0, 1/G_23, 0.0, 0.0],
+ [0.0, 0.0, 0.0, 0.0, 1/G_31, 0.0],
+ [0.0, 0.0, 0.0, 0.0, 0.0, 1/G_12]]);
+materialParameters_phase1 = compliance_S
+
+def prestrain_phase1(x):
+ # hB=delta_omega * alpha with delta_omega increment of moisture content and alpha swelling factor.
+ return [[1/param_h*delta_omega*alpha_L, 0, 0], [0,1/param_h*delta_omega*alpha_T,0], [0,0,1/param_h*delta_omega*alpha_R]]
+
+# --- PHASE 2
+# y_1-direction: R
+# y_2-direction: L
+# x_3-direction: T
+parameterSet.phase2_type="general_anisotropic"
+[E_1,E_2,E_3]=[E_R,E_L,E_T]
+[nu_12,nu_13,nu_23]=[nu_RL,nu_RT,nu_LT]
+[nu_21,nu_31,nu_32]=[nu_LR,nu_TR,nu_TL]
+[G_12,G_31,G_23]=[G_LR,G_RT,G_LT]
+compliance_S=np.array([[1/E_1, -nu_21/E_2, -nu_31/E_3, 0.0, 0.0, 0.0],
+ [-nu_12/E_1, 1/E_2, -nu_32/E_3, 0.0, 0.0, 0.0],
+ [-nu_13/E_1, -nu_23/E_2, 1/E_3, 0.0, 0.0, 0.0],
+ [0.0, 0.0, 0.0, 1/G_23, 0.0, 0.0],
+ [0.0, 0.0, 0.0, 0.0, 1/G_31, 0.0],
+ [0.0, 0.0, 0.0, 0.0, 0.0, 1/G_12]]);
+materialParameters_phase2 = compliance_S
+def prestrain_phase2(x):
+ return [[1/param_h*delta_omega*alpha_R, 0, 0], [0,1/param_h*delta_omega*alpha_L,0], [0,0,1/param_h*delta_omega*alpha_T]]
+
+#Rotation um 2. Achse (= L)
+parameterSet.phase2_axis = 1
+# phase2_angle = param_theta
+# -- Drehwinkel
+parameterSet.phase2_angle = param_theta
+
+
+# --- PHASE 3
+parameterSet.phase3_type="isotropic"
+epsilon = 1e-8
+materialParameters_phase3 = [epsilon, epsilon]
+
+def prestrain_phase3(x):
+ return [[0, 0, 0], [0,0,0], [0,0,0]]
+
+# # --- PHASE 3 = Phase 1 gedreht
+# # y_1-direction: L
+# # y_2-direction: R
+# # x_3-direction: T
+# parameterSet.phase3_type="general_anisotropic"
+# # Drehung um theta um Achse 2 = x_3-Achse
+# N=elast.rotation_matrix_compliance(2,param_theta)
+# materialParameters_phase3 = np.dot(np.dot(N,materialParameters_phase1),N.T)
+# materialParameters_phase3 = 0.5*(materialParameters_phase3.T+materialParameters_phase3)
+# # rotation of strain
+# def prestrain_phase3(x):
+# return elast.voigt_to_strain(np.dot(elast.rotation_matrix_compliance(2,param_theta),np.dot(elast.strain_to_voigt(np.array(prestrain_phase1(x))),N.T))).tolist()
+
+
+
+# --- Choose scale ratio gamma:
+parameterSet.gamma=1.0
+
+
+
+
+#############################################
+# Grid parameters
+#############################################
+## numLevels : Number of Levels on which solution is computed. starting with a 2x2x2 cube mesh.
+## {start,finish} computes on all grid from 2^(start) to 2^finish refinement
+#----------------------------------------------------
+parameterSet.numLevels= '3 3' # computes all levels from first to second entry
+# parameterSet.numLevels= '4 4'
+
+#############################################
+# Assembly options
+#############################################
+parameterSet.set_IntegralZero = 1 #(default = false)
+parameterSet.set_oneBasisFunction_Zero = 1 #(default = false)
+#parameterSet.arbitraryLocalIndex = 7 #(default = 0)
+#parameterSet.arbitraryElementNumber = 3 #(default = 0)
+
+#############################################
+# Solver Options, Type: #1: CG - SOLVER , #2: GMRES - SOLVER, #3: QR - SOLVER (default), #4: UMFPACK - SOLVER
+#############################################
+parameterSet.Solvertype = 3 # recommended to use iterative solver (e.g GMRES) for finer grid-levels
+parameterSet.Solver_verbosity = 0 #(default = 2) degree of information for solver output
+
+
+#############################################
+# Write/Output options #(default=false)
+#############################################
+# --- (Optional output) write Material / prestrain / Corrector functions to .vtk-Files:
+parameterSet.write_materialFunctions = 1 # VTK indicator function for material/prestrain definition
+#parameterSet.write_prestrainFunctions = 1 # VTK norm of B (currently not implemented)
+
+# --- (Additional debug output)
+parameterSet.print_debug = 0 #(default=false)
+
+# --- Write Correctos to VTK-File:
+parameterSet.write_VTK = 1
+
+# The grid can be refined several times for a higher resolution in the VTK-file.
+parameterSet.subsamplingRefinement = 2
+
+# --- (Optional output) L2Error, integral mean:
+#parameterSet.write_L2Error = 1
+#parameterSet.write_IntegralMean = 1
+
+# --- check orthogonality (75) from paper:
+parameterSet.write_checkOrthogonality = 0
+
+# --- Write corrector-coefficients to log-File:
+#parameterSet.write_corrector_phi1 = 1
+#parameterSet.write_corrector_phi2 = 1
+#parameterSet.write_corrector_phi3 = 1
+
+# --- Print Condition number of matrix (can be expensive):
+#parameterSet.print_conditionNumber= 1 #(default=false)
+
+# --- write effective quantities to Matlab-folder for symbolic minimization:
+parameterSet.write_toMATLAB = 1 # writes effective quantities to .txt-files QMatrix.txt and BMatrix.txt
diff --git a/experiment/wood-bilayer_PLOS/wood_bilayer_test.py b/experiment/wood-bilayer_PLOS/wood_bilayer_test.py
new file mode 100644
index 0000000000000000000000000000000000000000..e319d8a4784a6afcf5c3a77e052684516a8ee78a
--- /dev/null
+++ b/experiment/wood-bilayer_PLOS/wood_bilayer_test.py
@@ -0,0 +1,257 @@
+import subprocess
+import re
+import os
+import numpy as np
+import matplotlib.pyplot as plt
+import math
+import fileinput
+import time
+import matplotlib.ticker as tickers
+import matplotlib as mpl
+from matplotlib.ticker import MultipleLocator,FormatStrFormatter,MaxNLocator
+import codecs
+import sys
+import threading
+
+# # Schreibe input datei für Parameter
+# def SetParametersCellProblem(ParameterSet, ParsetFilePath, outputPath):
+# print('----set Parameters -----')
+# with open(ParsetFilePath, 'r') as file:
+# filedata = file.read()
+# filedata = re.sub('(?m)^materialFunction\s?=.*','materialFunction = '+str(ParameterSet.materialFunction),filedata)
+# filedata = re.sub('(?m)^gamma\s?=.*','gamma='+str(ParameterSet.gamma),filedata)
+# filedata = re.sub('(?m)^numLevels\s?=.*','numLevels='+str(ParameterSet.numLevels)+' '+str(ParameterSet.numLevels) ,filedata)
+# filedata = re.sub('(?m)^outputPath\s?=\s?.*','outputPath='+str(outputPath),filedata)
+# f = open(ParsetFilePath,'w')
+# f.write(filedata)
+# f.close()
+
+
+# # Ändere Parameter der MaterialFunction
+# def SetParameterMaterialFunction(materialFunction, parameterName, parameterValue):
+# with open(Path+"/"+materialFunction+'.py', 'r') as file:
+# filedata = file.read()
+# filedata = re.sub('(?m)^'+str(parameterName)+'\s?=.*',str(parameterName)+' = '+str(parameterValue),filedata)
+# f = open(Path+"/"+materialFunction+'.py','w')
+# f.write(filedata)
+# f.close()
+
+def SetParameterMaterialFunction(inputFunction, parameterName, parameterValue):
+ with open(inputFunction+'.py', 'r') as file:
+ filedata = file.read()
+ filedata = re.sub('(?m)^'+str(parameterName)+'\s?=.*',str(parameterName)+' = '+str(parameterValue),filedata)
+ f = open(inputFunction+'.py','w')
+ f.write(filedata)
+ f.close()
+
+
+
+# # Rufe Programm zum Lösen des Cell-Problems auf
+# def run_CellProblem(executable, parset,write_LOG):
+# print('----- RUN Cell-Problem ----')
+# processList = []
+# LOGFILE = "Cell-Problem_output.log"
+# print('LOGFILE:',LOGFILE)
+# print('executable:',executable)
+# if write_LOG:
+# p = subprocess.Popen(executable + parset
+# + " | tee " + LOGFILE, shell=True)
+
+# else:
+# p = subprocess.Popen(executable + parset
+# + " | tee " + LOGFILE, shell=True)
+
+# p.wait() # wait
+# processList.append(p)
+# exit_codes = [p.wait() for p in processList]
+
+# return
+
+# Read effective quantites
+def ReadEffectiveQuantities(QFilePath = os.path.dirname(os.getcwd()) + '/outputs/QMatrix.txt', BFilePath = os.path.dirname(os.getcwd())+ '/outputs/BMatrix.txt'):
+ # Read Output Matrices (effective quantities)
+ # From Cell-Problem output Files : ../outputs/Qmatrix.txt , ../outputs/Bmatrix.txt
+ # -- Read Matrix Qhom
+ X = []
+ # with codecs.open(path + '/outputs/QMatrix.txt', encoding='utf-8-sig') as f:
+ with codecs.open(QFilePath, encoding='utf-8-sig') as f:
+ for line in f:
+ s = line.split()
+ X.append([float(s[i]) for i in range(len(s))])
+ Q = np.array([[X[0][2], X[1][2], X[2][2]],
+ [X[3][2], X[4][2], X[5][2]],
+ [X[6][2], X[7][2], X[8][2]] ])
+
+ # -- Read Beff (as Vector)
+ X = []
+ # with codecs.open(path + '/outputs/BMatrix.txt', encoding='utf-8-sig') as f:
+ with codecs.open(BFilePath, encoding='utf-8-sig') as f:
+ for line in f:
+ s = line.split()
+ X.append([float(s[i]) for i in range(len(s))])
+ B = np.array([X[0][2], X[1][2], X[2][2]])
+ return Q, B
+
+# Function for evaluating the energy in terms of kappa, alpha and Q, B
+def eval_energy(kappa,alpha,Q,B) :
+ G=kappa*np.array([[np.cos(alpha)**2],[np.sin(alpha)**2],[np.sqrt(2)*np.cos(alpha)*np.sin(alpha)]])-B
+ return np.matmul(np.transpose(G),np.matmul(Q,G))[0,0]
+
+#-------------------------------------------------------------------------------------------------------
+########################
+#### SET PARAMETERS ####
+########################
+# ----- Setup Paths -----
+# write_LOG = True # writes Cell-Problem output-LOG in "Cell-Problem_output.log"
+# path='/home/klaus/Desktop/Dune_release/dune-microstructure/experiment/wood-bilayer/results/'
+# pythonPath = '/home/klaus/Desktop/Dune_release/dune-microstructure/experiment/wood-bilayer'
+
+#--- Choose wether to perforate upper (passive) or lower (active) layer
+# perforatedLayer = 'upper'
+perforatedLayer = 'lower'
+
+dataset_numbers = [0, 1, 2, 3, 4, 5]
+
+for dataset_number in dataset_numbers:
+ print("------------------")
+ print(str(dataset_number) + "th data set")
+ print("------------------")
+
+
+ # path = os.getcwd() + '/experiment/perforated-bilayer/results_' + perforatedLayer + '/'
+ path = os.getcwd() + '/experiment/wood-bilayer_PLOS/results_' + str(dataset_number) + '/'
+ pythonPath = os.getcwd() + '/experiment/wood-bilayer_PLOS'
+ pythonModule = "perforated_wood_" + perforatedLayer
+ executable = os.getcwd() + '/build-cmake/src/Cell-Problem'
+ # ---------------------------------
+ # Setup Experiment
+ # ---------------------------------
+ gamma = 1.0
+
+ # ----- Define Parameters for Material Function --------------------
+ # [r, h, omega_flat, omega_target, theta, beta]
+ # r = (thickness upper layer)/(thickness)
+ # h = thickness [meter]
+ # omega_flat = moisture content in the flat state before drying [%]
+ # omega_target = moisture content in the target state [%]
+ # theta = rotation angle (not implemented and used)
+ # beta = design parameter for perforation = ratio of Volume of (cylindrical) perforation to Volume of active/passive layer
+
+ #--- Different moisture values for different thicknesses:
+
+ materialFunctionParameter=[
+ [ # Dataset Ratio r = 0.12
+ [0.12, 0.0047, 17.32986047, 14.70179844, 0.0, 0.0 ],
+ [0.12, 0.0047, 17.32986047, 13.6246, 0.0, 0],
+ [0.12, 0.0047, 17.32986047, 12.42994508, 0.0, 0 ],
+ [0.12, 0.0047, 17.32986047, 11.69773413, 0.0, 0],
+ [0.12, 0.0047, 17.32986047, 11.14159987, 0.0, 0],
+ [0.12, 0.0047, 17.32986047, 9.500670278, 0.0, 0],
+ [0.12, 0.0047, 17.32986047, 9.005046347, 0.0, 0]
+ ],
+ [ # Dataset Ratio r = 0.17
+ [0.17, 0.0049, 17.28772791 , 14.75453569, 0.0, 0 ],
+ [0.17, 0.0049, 17.28772791 , 13.71227639, 0.0, 0],
+ [0.17, 0.0049, 17.28772791 , 12.54975012, 0.0, 0 ],
+ [0.17, 0.0049, 17.28772791 , 11.83455959, 0.0, 0],
+ [0.17, 0.0049, 17.28772791 , 11.29089521, 0.0, 0 ],
+ [0.17, 0.0049, 17.28772791 , 9.620608917, 0.0, 0],
+ [0.17, 0.0049, 17.28772791 , 9.101671742, 0.0, 0 ]
+ ],
+ [ # Dataset Ratio r = 0.22
+ [0.22, 0.0053, 17.17547062, 8.959564147, 0.0, 0 ],
+ [0.22, 0.0053, 17.17547062, 8.959564147, 0.0, 0 ],
+ [0.22, 0.0053, 17.17547062, 8.959564147, 0.0, 0 ],
+ [0.22, 0.0053, 17.17547062, 8.959564147, 0.0, 0 ],
+ [0.22, 0.0053, 17.17547062, 8.959564147, 0.0, 0 ],
+ [0.22, 0.0053, 17.17547062, 8.959564147, 0.0, 0 ],
+ [0.22, 0.0053, 17.17547062, 8.959564147, 0.0, 0 ]
+ ],
+ [ # Dataset Ratio r = 0.34
+ [0.34, 0.0063, 17.14061081 , 14.98380876, 0.0, 0 ],
+ [0.34, 0.0063, 17.14061081 , 13.97154915, 0.0, 0],
+ [0.34, 0.0063, 17.14061081 , 12.77309253, 0.0, 0 ],
+ [0.34, 0.0063, 17.14061081 , 12.00959929, 0.0, 0],
+ [0.34, 0.0063, 17.14061081 , 11.42001731, 0.0, 0 ],
+ [0.34, 0.0063, 17.14061081 , 9.561447179, 0.0, 0],
+ [0.34, 0.0063, 17.14061081 , 8.964704969, 0.0, 0 ]
+ ],
+ [ # Dataset Ratio r = 0.43
+ [0.43, 0.0073, 17.07559686 , 15.11316339, 0.0, 0 ],
+ [0.43, 0.0073, 17.07559686 , 14.17997082, 0.0, 0],
+ [0.43, 0.0073, 17.07559686 , 13.05739844, 0.0, 0 ],
+ [0.43, 0.0073, 17.07559686 , 12.32309209, 0.0, 0],
+ [0.43, 0.0073, 17.07559686 , 11.74608518, 0.0, 0 ],
+ [0.43, 0.0073, 17.07559686 , 9.812372466, 0.0, 0],
+ [0.43, 0.0073, 17.07559686 , 9.10519385 , 0.0, 0 ]
+ ],
+ [ # Dataset Ratio r = 0.49
+ [0.49, 0.008, 17.01520754, 15.30614414, 0.0, 0 ],
+ [0.49, 0.008, 17.01520754, 14.49463867, 0.0, 0],
+ [0.49, 0.008, 17.01520754, 13.46629742, 0.0, 0 ],
+ [0.49, 0.008, 17.01520754, 12.78388234, 0.0, 0],
+ [0.49, 0.008, 17.01520754, 12.23057715, 0.0, 0 ],
+ [0.49, 0.008, 17.01520754, 10.21852839, 0.0, 0],
+ [0.49, 0.008, 17.01520754, 9.341730605, 0.0, 0 ]
+ ]
+ ]
+
+
+
+ # ------ Loops through Parameters for Material Function -----------
+ for i in range(0,np.shape(materialFunctionParameter)[1]):
+ print("------------------")
+ print("New Loop")
+ print("------------------")
+ # Check output directory
+ outputPath = path + str(i)
+ isExist = os.path.exists(outputPath)
+ if not isExist:
+ # Create a new directory because it does not exist
+ os.makedirs(outputPath)
+ print("The new directory " + outputPath + " is created!")
+
+ # thread = threading.Thread(target=run_CellProblem(executable, pythonModule, pythonPath, LOGFILE))
+ # thread.start()
+
+ #TODO: apperently its not possible to pass a variable via subprocess and "calculate" another input value inside the python file.
+ # Therefore we use this instead.
+ SetParameterMaterialFunction(pythonPath + "/" + pythonModule, "param_r",materialFunctionParameter[dataset_number][i][0])
+ SetParameterMaterialFunction(pythonPath + "/" + pythonModule, "param_h",materialFunctionParameter[dataset_number][i][1])
+ SetParameterMaterialFunction(pythonPath + "/" + pythonModule, "param_omega_flat",materialFunctionParameter[dataset_number][i][2])
+ SetParameterMaterialFunction(pythonPath + "/" + pythonModule, "param_omega_target",materialFunctionParameter[dataset_number][i][3])
+ SetParameterMaterialFunction(pythonPath + "/" + pythonModule, "param_theta",materialFunctionParameter[dataset_number][i][4])
+ SetParameterMaterialFunction(pythonPath + "/" + pythonModule, "param_beta",materialFunctionParameter[dataset_number][i][5])
+
+ LOGFILE = outputPath + "/" + pythonModule + "_output" + "_" + str(i) + ".log"
+
+ processList = []
+ p = subprocess.Popen(executable + " " + pythonPath + " " + pythonModule
+ + " -outputPath " + outputPath
+ + " -gamma " + str(gamma)
+ + " | tee " + LOGFILE, shell=True)
+
+ # p = subprocess.Popen(executable + " " + pythonPath + " " + pythonModule
+ # + " -outputPath " + outputPath
+ # + " -gamma " + str(gamma)
+ # + " -param_r " + str(materialFunctionParameter[i][0])
+ # + " -param_h " + str(materialFunctionParameter[i][1])
+ # + " -param_omega_flat " + str(materialFunctionParameter[i][2])
+ # + " -param_omega_target " + str(materialFunctionParameter[i][3])
+ # + " -phase2_angle " + str(materialFunctionParameter[i][4])
+ # + " | tee " + LOGFILE, shell=True)
+
+ p.wait() # wait
+ processList.append(p)
+ exit_codes = [p.wait() for p in processList]
+ # ---------------------------------------------------
+ # wait here for the result to be available before continuing
+ # thread.join()
+ f = open(outputPath+"/parameter.txt", "w")
+ f.write("r = "+str(materialFunctionParameter[dataset_number][i][0])+"\n")
+ f.write("h = "+str(materialFunctionParameter[dataset_number][i][1])+"\n")
+ f.write("omega_flat = "+str(materialFunctionParameter[dataset_number][i][2])+"\n")
+ f.write("omega_target = "+str(materialFunctionParameter[dataset_number][i][3])+"\n")
+ f.write("param_beta = "+str(materialFunctionParameter[dataset_number][i][5])+"\n")
+ f.close()
+ #