From f68601ec024cff93e4834cdeb97364d10000d86a Mon Sep 17 00:00:00 2001
From: Klaus <klaus.boehnlein@tu-dresden.de>
Date: Tue, 10 Oct 2023 18:40:05 +0200
Subject: [PATCH] Test/Add rotation in wood_test example
---
.../isotrop_orthotrop_rotation.py | 11 +++--
experiment/rotation-test/rotation_test.py | 1 -
.../wood-bilayer/PolarPlotLocalEnergy.py | 4 +-
.../wood-bilayer/wood_european_beech.py | 42 +++++++++++-------
experiment/wood-bilayer/wood_test.py | 44 ++++++++++++++++---
5 files changed, 72 insertions(+), 30 deletions(-)
diff --git a/experiment/rotation-test/isotrop_orthotrop_rotation.py b/experiment/rotation-test/isotrop_orthotrop_rotation.py
index 6a5b26f6..2ea848b9 100644
--- a/experiment/rotation-test/isotrop_orthotrop_rotation.py
+++ b/experiment/rotation-test/isotrop_orthotrop_rotation.py
@@ -44,7 +44,7 @@ prestrain_wood=np.array([[-0.050821460301886785, 0, 0],
[0, 0, -0.8824453561509432]])
# rotation angle
# param_theta = '1.5707963267948966'
-param_theta = '0.5'
+
# untere Schicht: isotrop, ohne prestrain
@@ -71,10 +71,9 @@ parameterSet.phase1_type="general_anisotropic"
# Drehung um theta um Achse 2 = x_3-Achse
parameterSet.phase1_axis = 2
-# phase1_angle = 0.0
-# phase1_angle = float("1.5707963267948966")
-# phase1_angle = float(param_theta)
-parameterSet.phase1_angle = 0.5
+
+parameterSet.phase1_angle = 0.0
+# parameterSet.phase1_angle = param_theta
materialParameters_phase1 = compliance_wood
@@ -146,7 +145,7 @@ parameterSet.subsamplingRefinement = 2
#parameterSet.write_IntegralMean = 1
# --- check orthogonality (75) from paper:
-parameterSet.write_checkOrthogonality = 1
+parameterSet.write_checkOrthogonality = 0
# --- Write corrector-coefficients to log-File:
#parameterSet.write_corrector_phi1 = 1
diff --git a/experiment/rotation-test/rotation_test.py b/experiment/rotation-test/rotation_test.py
index ed34f34a..431e5f40 100644
--- a/experiment/rotation-test/rotation_test.py
+++ b/experiment/rotation-test/rotation_test.py
@@ -164,7 +164,6 @@ for i in range(0,np.shape(materialFunctionParameter)[0]):
p = subprocess.Popen(executable + " " + pythonPath + " " + pythonModule
+ " -outputPath " + outputPath
+ " -gamma " + str(gamma)
- + " -param_theta " + str(materialFunctionParameter[i])
+ " -phase1_angle " + str(materialFunctionParameter[i])
+ " | tee " + LOGFILE, shell=True)
diff --git a/experiment/wood-bilayer/PolarPlotLocalEnergy.py b/experiment/wood-bilayer/PolarPlotLocalEnergy.py
index 8a7cc06c..6e623b15 100644
--- a/experiment/wood-bilayer/PolarPlotLocalEnergy.py
+++ b/experiment/wood-bilayer/PolarPlotLocalEnergy.py
@@ -63,8 +63,8 @@ for n in range(0,number):
B=np.transpose([B])
#
- N=500
- length=5
+ 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):
diff --git a/experiment/wood-bilayer/wood_european_beech.py b/experiment/wood-bilayer/wood_european_beech.py
index c8c326ab..135ca8b2 100644
--- a/experiment/wood-bilayer/wood_european_beech.py
+++ b/experiment/wood-bilayer/wood_european_beech.py
@@ -49,8 +49,9 @@ param_h = 0.0053
param_omega_flat = 17.17547062
# -- moisture content in the target state [%]
param_omega_target = 8.959564147
-# -- Drehwinkel (nicht implementiert)
+# -- Drehwinkel
param_theta = 0
+
#
#
#
@@ -151,18 +152,27 @@ 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]]
-# --- 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()
+#Rotation um 2. Achse (= L)
+parameterSet.phase2_axis = 1
+# phase2_angle = param_theta
+# -- Drehwinkel
+parameterSet.phase2_angle = param_theta
+
+
+
+# # --- 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:
@@ -177,7 +187,7 @@ parameterSet.gamma=1.0
## 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= '4 4' # computes all levels from first to second entry
+parameterSet.numLevels= '3 3' # computes all levels from first to second entry
#############################################
@@ -191,7 +201,7 @@ parameterSet.set_oneBasisFunction_Zero = 1 #(default = false)
#############################################
# Solver Options, Type: #1: CG - SOLVER , #2: GMRES - SOLVER, #3: QR - SOLVER (default), #4: UMFPACK - SOLVER
#############################################
-parameterSet.Solvertype = 2 # recommended to use iterative solver (e.g GMRES) for finer grid-levels
+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
@@ -216,7 +226,7 @@ parameterSet.subsamplingRefinement = 2
#parameterSet.write_IntegralMean = 1
# --- check orthogonality (75) from paper:
-parameterSet.write_checkOrthogonality = 1
+parameterSet.write_checkOrthogonality = 0
# --- Write corrector-coefficients to log-File:
#parameterSet.write_corrector_phi1 = 1
diff --git a/experiment/wood-bilayer/wood_test.py b/experiment/wood-bilayer/wood_test.py
index 6f1d5df7..5844cf4c 100644
--- a/experiment/wood-bilayer/wood_test.py
+++ b/experiment/wood-bilayer/wood_test.py
@@ -36,6 +36,15 @@ import threading
# 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):
@@ -113,6 +122,8 @@ gamma = 1.0
# omega_target = moisture content in the target state [%]
# theta = rotation angle (not implemented and used)
# experimental_kappa = curvature measure in experiment
+
+#First Experiment:
materialFunctionParameter=[
[0.22, 0.0053, 17.17547062, 14.72680026, 0, 1.058078122],
[0.22, 0.0053, 17.17547062, 13.64338887, 0, 1.544624544],
@@ -122,6 +133,14 @@ materialFunctionParameter=[
[0.22, 0.0053, 17.17547062, 9.435795985, 0, 3.913528418],
[0.22, 0.0053, 17.17547062, 8.959564147, 0, 4.262750825]
]
+
+#Second Experiment: Rotate "active" bilayer phase
+# materialFunctionParameter=[
+# [0.22, 0.0053, 17.17547062, 8.959564147, 0, 4.262750825],
+# [0.22, 0.0053, 17.17547062, 8.959564147, (np.pi/4.0), 4.262750825],
+# [0.22, 0.0053, 17.17547062, 8.959564147, (np.pi/2.0), 4.262750825]
+# ]
+
# ------ Loops through Parameters for Material Function -----------
for i in range(0,np.shape(materialFunctionParameter)[0]):
print("------------------")
@@ -137,18 +156,33 @@ for i in range(0,np.shape(materialFunctionParameter)[0]):
# 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])
+
LOGFILE = outputPath + "/" + pythonModule + "_output" + "_" + str(i) + ".log"
+
processList = []
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])
- + " -param_theta " + str(materialFunctionParameter[i][4])
+ " | 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]
--
GitLab