diff --git a/experiment/perforated-bilayer/PolarPlotLocalEnergy.py b/experiment/perforated-bilayer/PolarPlotLocalEnergy.py
new file mode 100644
index 0000000000000000000000000000000000000000..4ed2e4aae50f9a3f985db2ed7c174b1af2281e17
--- /dev/null
+++ b/experiment/perforated-bilayer/PolarPlotLocalEnergy.py
@@ -0,0 +1,95 @@
+#!/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 of experiments / folders
+number=4
+kappa=np.zeros(number)
+for n in range(0,number):
+ # Read from Date
+ print(str(n))
+ DataPath = './experiment/perforated-bilayer/results/'+str(n)
+ DataPath = './experiment/perforated-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=300
+ length=12
+ 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/perforated-bilayer/results/kappa_simulation.txt", "w")
+f.write(str(kappa))
+f.close()
+
+
diff --git a/experiment/perforated-bilayer/perfBilayer_test.py b/experiment/perforated-bilayer/perfBilayer_test.py
new file mode 100644
index 0000000000000000000000000000000000000000..c580079027edf990fd4fa9ad5801fe2e66a8e80b
--- /dev/null
+++ b/experiment/perforated-bilayer/perfBilayer_test.py
@@ -0,0 +1,193 @@
+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'
+path = os.getcwd() + '/experiment/perforated-bilayer/results/'
+pythonPath = os.getcwd() + '/experiment/perforated-bilayer'
+# pythonModule = "perforated_wood_upper"
+pythonModule = "perforated_wood_lower"
+executable = os.getcwd() + '/build-cmake/src/Cell-Problem'
+# ---------------------------------
+# Setup Experiment
+# ---------------------------------
+gamma = 1.0
+
+# ----- Define Parameters for Material Function --------------------
+# [r, h, omega_flat, omega_target, theta, experimental_kappa]
+# r = (thickness upper layer)/(thickness)
+# h = thickness [meter]
+# omega_flat = moisture content in the flat state before drying [%]
+# omega_target = moisture content in the target state [%]
+# theta = rotation angle (not implemented and used)
+# beta = design parameter for perforation = ratio of Volume of (cylindrical) perforation to Volume of active/passive layer
+
+#Experiment: Perforate "active" bilayer phase
+materialFunctionParameter=[
+ [0.22, 0.0053, 17.17547062, 8.959564147, 0.0, 0.0 ],
+ [0.22, 0.0053, 17.17547062, 8.959564147, 0.0, 0.1 ],
+ [0.22, 0.0053, 17.17547062, 8.959564147, 0.0, 0.2 ],
+ [0.22, 0.0053, 17.17547062, 8.959564147, 0.0, 0.3 ]
+ # [0.22, 0.0053, 17.17547062, 8.959564147, 0.0, 0.4 ],
+ # [0.22, 0.0053, 17.17547062, 8.959564147, 0.0, 0.5 ]
+]
+
+# ------ Loops through Parameters for Material Function -----------
+for i in range(0,np.shape(materialFunctionParameter)[0]):
+ 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[i][0])
+ SetParameterMaterialFunction(pythonPath + "/" + pythonModule, "param_h",materialFunctionParameter[i][1])
+ SetParameterMaterialFunction(pythonPath + "/" + pythonModule, "param_omega_flat",materialFunctionParameter[i][2])
+ SetParameterMaterialFunction(pythonPath + "/" + pythonModule, "param_omega_target",materialFunctionParameter[i][3])
+ SetParameterMaterialFunction(pythonPath + "/" + pythonModule, "param_theta",materialFunctionParameter[i][4])
+ SetParameterMaterialFunction(pythonPath + "/" + pythonModule, "param_beta",materialFunctionParameter[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[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.write("param_beta = "+str(materialFunctionParameter[i][5])+"\n")
+ f.close()
+ #
diff --git a/experiment/perforated-bilayer/perforated_wood_lower.py b/experiment/perforated-bilayer/perforated_wood_lower.py
new file mode 100644
index 0000000000000000000000000000000000000000..95730419f1f1a8d1c8c5dce00f350b429ca24edb
--- /dev/null
+++ b/experiment/perforated-bilayer/perforated_wood_lower.py
@@ -0,0 +1,278 @@
+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/np.pi)
+ 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): #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.22
+# -- thickness [meter]
+param_h = 0.0053
+# -- 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
+
+#
+#
+#
+# -- 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
+
+
+#############################################
+# 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/perforated-bilayer/perforated_wood_upper.py b/experiment/perforated-bilayer/perforated_wood_upper.py
new file mode 100644
index 0000000000000000000000000000000000000000..9316c2b1dadd3031a71b6c0e1cc073da8042b80b
--- /dev/null
+++ b/experiment/perforated-bilayer/perforated_wood_upper.py
@@ -0,0 +1,278 @@
+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/np.pi)
+ if (x[2]>=(0.5-param_r)):
+ # if(np.sqrt(x[0]**2 + x[1]**2) < pRadius): #inside perforation
+ if((x[0]**2 + x[1]**2) < pRadius**2): #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.22
+# -- thickness [meter]
+param_h = 0.0053
+# -- 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
+
+#
+#
+#
+# -- 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
+
+
+#############################################
+# 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