diff --git a/experiment/rotation-test/PolarPlotLocalEnergy.py b/experiment/rotation-test/PolarPlotLocalEnergy.py
new file mode 100644
index 0000000000000000000000000000000000000000..67429fda67c053e403a735c2573004e2df30b7a7
--- /dev/null
+++ b/experiment/rotation-test/PolarPlotLocalEnergy.py
@@ -0,0 +1,94 @@
+#!/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=2
+kappa=np.zeros(number)
+for n in range(0,number):
+ # Read from Date
+ print(str(n))
+ Path='./experiment/rotation-test/results/'
+ DataPath = Path+str(n)
+ QFilePath = DataPath + '/QMatrix.txt'
+ BFilePath = DataPath + '/BMatrix.txt'
+ Q, B = ReadEffectiveQuantities(QFilePath,BFilePath)
+ # Q=0.5*(np.transpose(Q)+Q) # symmetrize
+ B=np.transpose([B])
+ #
+
+ N=500
+ length=5
+ r, theta = np.meshgrid(np.linspace(0,length,N),np.radians(np.linspace(0, 360, N)))
+ E=np.zeros(np.shape(r))
+ for i in range(0,N):
+ for j in range(0,N):
+ if theta[i,j]<np.pi:
+ E[i,j]=energy(r[i,j],theta[i,j],Q,B)
+ else:
+ E[i,j]=energy(-r[i,j],theta[i,j],Q,B)
+
+ # Compute Minimizer
+ [imin,jmin]=np.unravel_index(E.argmin(),(N,N))
+ kappamin=r[imin,jmin]
+ alphamin=theta[imin,jmin]
+ kappa[n]=kappamin
+ fig, ax = plt.subplots(figsize=(6,6),subplot_kw=dict(projection='polar'))
+ levs=np.geomspace(E.min(),E.max(),400)
+ pcm=ax.contourf(theta, r, E, levs, norm=colors.PowerNorm(gamma=0.2), cmap='brg')
+ ax.set_xticks([0,np.pi/2])
+ ax.set_yticks([kappamin])
+ colorbarticks=np.linspace(E.min(),E.max(),6)
+ plt.colorbar(pcm, extend='max', ticks=colorbarticks, pad=0.1)
+ plt.show()
+
+f = open("./experiment/wood-bilayer/results/kappa_simulation.txt", "w")
+f.write(str(kappa))
+f.close()
+
+
diff --git a/experiment/rotation-test/cellsolver.parset b/experiment/rotation-test/cellsolver.parset
new file mode 100644
index 0000000000000000000000000000000000000000..e674549614e510e493c929fe6d48145b78df7f03
--- /dev/null
+++ b/experiment/rotation-test/cellsolver.parset
@@ -0,0 +1,96 @@
+# --- Parameter File as Input for 'Cell-Problem'
+# NOTE: define variables without whitespaces in between! i.e. : gamma=1.0 instead of gamma = 1.0
+# since otherwise these cant be read from other Files!
+# --------------------------------------------------------
+
+# Path for results and logfile
+outputPath=./experiment/rotation-test/results/0
+
+# Path for material description
+geometryFunctionPath =experiment/rotation-test/
+
+
+# --- DEBUG (Output) Option:
+#print_debug = true #(default=false)
+
+
+
+
+#############################################
+# Grid parameters
+#############################################
+#----------------------------------------------------
+## numLevels : Number of Levels on which solution is computed. starting with a 2x2x2 cube mesh.
+## {start,finish} computes on all grid from 2^(start) to 2^finish refinement
+#----------------------------------------------------
+
+numLevels=4 4
+#numLevels = 1 1 # computes all levels from first to second entry
+#numLevels = 2 2 # computes all levels from first to second entry
+#numLevels = 1 3 # computes all levels from first to second entry
+#numLevels = 4 4 # computes all levels from first to second entry
+#numLevels = 5 5 # computes all levels from first to second entry
+#numLevels = 6 6 # computes all levels from first to second entry
+#numLevels = 1 6
+
+
+#############################################
+# Material / Prestrain parameters and ratios
+#############################################
+
+# --- Choose material definition:
+materialFunction = isotrop_orthotrop_rotation
+
+
+
+# --- Choose scale ratio gamma:
+gamma=1.0
+
+
+#############################################
+# Assembly options
+#############################################
+#set_IntegralZero = true #(default = false)
+#set_oneBasisFunction_Zero = true #(default = false)
+
+#arbitraryLocalIndex = 7 #(default = 0)
+#arbitraryElementNumber = 3 #(default = 0)
+#############################################
+
+
+#############################################
+# Solver Type: #1: CG - SOLVER , #2: GMRES - SOLVER, #3: QR - SOLVER (default), #4: UMFPACK - SOLVER
+#############################################
+Solvertype = 2 # recommended to use iterative solver (e.g GMRES) for finer grid-levels
+Solver_verbosity = 0 #(default = 2) degree of information for solver output
+
+
+#############################################
+# Write/Output options #(default=false)
+#############################################
+
+# --- (Optional output) write Material / prestrain / Corrector functions to .vtk-Files:
+write_materialFunctions = true # VTK indicator function for material/prestrain definition
+#write_prestrainFunctions = true # VTK norm of B (currently not implemented)
+
+# --- Write Correctos to VTK-File:
+write_VTK = true
+
+# --- (Optional output) L2Error, integral mean:
+#write_L2Error = true
+#write_IntegralMean = true
+
+# --- check orthogonality (75) from paper:
+write_checkOrthogonality = true
+
+# --- Write corrector-coefficients to log-File:
+#write_corrector_phi1 = true
+#write_corrector_phi2 = true
+#write_corrector_phi3 = true
+
+
+# --- Print Condition number of matrix (can be expensive):
+#print_conditionNumber= true #(default=false)
+
+# --- write effective quantities to Matlab-folder for symbolic minimization:
+write_toMATLAB = true # writes effective quantities to .txt-files QMatrix.txt and BMatrix.txt
diff --git a/experiment/rotation-test/elasticity_toolbox.py b/experiment/rotation-test/elasticity_toolbox.py
new file mode 100644
index 0000000000000000000000000000000000000000..8e61952612c0714a5b430a41660775fc0e2c23b5
--- /dev/null
+++ b/experiment/rotation-test/elasticity_toolbox.py
@@ -0,0 +1,123 @@
+import math
+import numpy as np
+
+
+def strain_to_voigt(strain_matrix):
+ # Ensure the input matrix is a 3x3 strain matrix
+ if strain_matrix.shape != (3, 3):
+ raise ValueError("Input matrix should be a 3x3 strain matrix.")
+
+ # Extract the components from the 3x3 strain matrix
+ ε_xx = strain_matrix[0, 0]
+ ε_yy = strain_matrix[1, 1]
+ ε_zz = strain_matrix[2, 2]
+ γ_yz = .5*(strain_matrix[1, 2]+strain_matrix[2,1])
+ γ_xz = .5*(strain_matrix[0, 2]+strain_matrix[0,2])
+ γ_xy = .5*(strain_matrix[0, 1]+strain_matrix[0,1])
+
+ # Create the Voigt notation vector
+ voigt_notation = np.array([ε_xx, ε_yy, ε_zz, γ_yz, γ_xz, γ_xy])
+
+ return voigt_notation
+
+def voigt_to_strain(voigt_notation):
+ # Ensure the input vector has 6 elements
+ if len(voigt_notation) != 6:
+ raise ValueError("Input vector should have 6 elements in Voigt notation.")
+
+ # Extract the components from the Voigt notation vector
+ ε_xx = voigt_notation[0]
+ ε_yy = voigt_notation[1]
+ ε_zz = voigt_notation[2]
+ γ_yz = voigt_notation[3]
+ γ_xz = voigt_notation[4]
+ γ_xy = voigt_notation[5]
+
+ # Create the 3x3 strain matrix
+ strain_matrix = np.array([[ε_xx, γ_xy, γ_xz],
+ [γ_xy, ε_yy, γ_yz],
+ [γ_xz, γ_yz, ε_zz]])
+
+ return strain_matrix
+
+
+def rotation_matrix(ax, angle):
+ cos_theta = np.cos(angle)
+ sin_theta = np.sin(angle)
+ if ax==0:
+ Q=np.array([[0, 0, 1],
+ [0,1,0],
+ [-1,0,0]
+ ])
+ elif ax==1:
+ Q=np.array([[1, 0, 0],
+ [0,0,1],
+ [0,-1,0]
+ ])
+ else:
+ Q=np.array([[1, 0, 0],
+ [0,1,0],
+ [0,0,1]
+ ])
+
+ R = np.array([[cos_theta, -sin_theta, 0],
+ [sin_theta, cos_theta, 0],
+ [0, 0, 1]])
+ return np.dot(np.dot(Q.T, R),Q)
+
+def rotation_matrix_compliance(ax,theta):
+ R=rotation_matrix(ax,theta)
+ Q_xx=R[0,0]
+ Q_xy=R[0,1]
+ Q_xz=R[0,2]
+ Q_yx=R[1,0]
+ Q_yy=R[1,1]
+ Q_yz=R[1,2]
+ Q_zx=R[2,0]
+ Q_zy=R[2,1]
+ Q_zz=R[2,2]
+ return np.array([
+ [Q_xx**2, Q_xy**2, Q_xz**2, Q_xy*Q_xz, Q_xx*Q_xz, Q_xx*Q_xy],
+ [Q_yx**2, Q_yy**2, Q_yz**2, Q_yy*Q_yz, Q_yx*Q_yz, Q_yx*Q_yy],
+ [Q_zx**2, Q_zy**2, Q_zz**2, Q_zy*Q_zz, Q_zx*Q_zz, Q_zx*Q_zy],
+ [2*Q_yx*Q_zx, 2*Q_yy*Q_zy, 2*Q_yz*Q_zz, Q_yy*Q_zz + Q_yz*Q_zy, Q_yx*Q_zz + Q_yz*Q_zx, Q_yx*Q_zy + Q_yy*Q_zx],
+ [2*Q_xx*Q_zx, 2*Q_xy*Q_zy, 2*Q_xz*Q_zz, Q_xy*Q_zz + Q_xz*Q_zy, Q_xx*Q_zz + Q_xz*Q_zx, Q_xx*Q_zy + Q_xy*Q_zx],
+ [2*Q_xx*Q_yx, 2*Q_xy*Q_yy, 2*Q_xz*Q_yz, Q_xy*Q_yz + Q_xz*Q_yy, Q_xx*Q_yz + Q_xz*Q_yx, Q_xx*Q_yy + Q_xy*Q_yx]
+ ])
+
+def rotate_strain(eps, ax, theta):
+ B=voigt_to_strain(np.matmul(rotation_matrix_epsilon(theta,ax),strain_to_voigt(eps)))
+
+import numpy as np
+
+def voigt_to_tensor(voigt_matrix):
+ tensor = np.zeros((6, 6))
+
+ tensor[0, 0] = voigt_matrix[0]
+ tensor[0, 1] = tensor[1, 0] = voigt_matrix[1]
+ tensor[0, 2] = tensor[2, 0] = voigt_matrix[2]
+ tensor[0, 3] = tensor[3, 0] = voigt_matrix[3]
+ tensor[0, 4] = tensor[4, 0] = voigt_matrix[4]
+ tensor[0, 5] = tensor[5, 0] = voigt_matrix[5]
+
+ tensor[1, 1] = voigt_matrix[6]
+ tensor[1, 2] = tensor[2, 1] = voigt_matrix[7]
+ tensor[1, 3] = tensor[3, 1] = voigt_matrix[8]
+ tensor[1, 4] = tensor[4, 1] = voigt_matrix[9]
+ tensor[1, 5] = tensor[5, 1] = voigt_matrix[10]
+
+ tensor[2, 2] = voigt_matrix[11]
+ tensor[2, 3] = tensor[3, 2] = voigt_matrix[12]
+ tensor[2, 4] = tensor[4, 2] = voigt_matrix[13]
+ tensor[2, 5] = tensor[5, 2] = voigt_matrix[14]
+
+ tensor[3, 3] = voigt_matrix[15]
+ tensor[3, 4] = tensor[4, 3] = voigt_matrix[16]
+ tensor[3, 5] = tensor[5, 3] = voigt_matrix[17]
+
+ tensor[4, 4] = voigt_matrix[18]
+ tensor[4, 5] = tensor[5, 4] = voigt_matrix[19]
+
+ tensor[5, 5] = voigt_matrix[20]
+
+ return tensor
diff --git a/experiment/rotation-test/isotrop_orthotrop_rotation.py b/experiment/rotation-test/isotrop_orthotrop_rotation.py
new file mode 100644
index 0000000000000000000000000000000000000000..d4809134ba6dce598d4a59ef62758959b193e2cd
--- /dev/null
+++ b/experiment/rotation-test/isotrop_orthotrop_rotation.py
@@ -0,0 +1,62 @@
+import math
+#from python_matrix_operations import *
+import ctypes
+import os
+import sys
+import numpy as np
+import elasticity_toolbox as elast
+
+# Komposit beschreibt ein Material mit
+# oberer Schicht: orthotrop (Holz), in der Ebene gedreht
+# orientierung y_1-direction: L, y_2-direction: T, x_3-direction: R
+compliance_wood=np.array([[ 7.70543562e-05, -1.56584619e-05, -1.96144058e-05,
+ 0.00000000e+00, 0.00000000e+00, 0.00000000e+00],
+ [-1.56584619e-05, 1.84913679e-03, -5.14793170e-04,
+ 0.00000000e+00, 0.00000000e+00, 0.00000000e+00],
+ [-1.96144058e-05, -5.14793170e-04, 5.92975602e-04,
+ 0.00000000e+00, 0.00000000e+00, 0.00000000e+00],
+ [ 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
+ 2.25174559e-03, 0.00000000e+00, 0.00000000e+00],
+ [ 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
+ 0.00000000e+00, 7.95374719e-04, 0.00000000e+00],
+ [ 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
+ 0.00000000e+00, 0.00000000e+00, 1.19183673e-03]])
+prestrain_wood=np.array([[-0.050821460301886785, 0, 0],
+ [0, -2.134501332679245, 0],
+ [0, 0, -0.8824453561509432]])
+# rotation angle
+param_theta = 0
+
+# untere Schicht: isotrop, ohne prestrain
+
+# --- define geometry
+
+def indicatorFunction(x):
+ factor=1
+ if (x[2]>=0):
+ return 1 #Phase1
+ else :
+ return 2 #Phase2
+
+# --- Number of material phases
+Phases=2
+
+# --- PHASE 1
+phase1_type="general_anisotropic"
+# Drehung um theta um Achse 2 = x_3-Achse
+N=elast.rotation_matrix_compliance(2,param_theta)
+materialParameters_phase1 = np.dot(np.dot(N,compliance_wood),N.T)
+materialParameters_phase1 = (0.5*(materialParameters_phase1.T+materialParameters_phase1)).tolist()
+# rotation of strain
+def prestrain_phase1(x):
+ return elast.voigt_to_strain(np.dot(elast.rotation_matrix_compliance(2,param_theta),np.dot(elast.strain_to_voigt(np.array(prestrain_wood)),N.T))).tolist()
+
+# --- PHASE 2
+phase2_type="isotropic"
+# extract E and nu from wood
+E = 1/compliance_wood[0,0]
+nu = -compliance_wood[0,1]*E
+# [mu, lambda]
+materialParameters_phase2 = [E/(2*(1+nu)), (E*nu)/((1+nu)*(1-2*nu))]
+def prestrain_phase2(x):
+ return [[0,0,0],[0,0,0],[0,0,0]]
diff --git a/experiment/rotation-test/rotation_test.py b/experiment/rotation-test/rotation_test.py
new file mode 100644
index 0000000000000000000000000000000000000000..340617078264671212be03fa925f2f3367744eae
--- /dev/null
+++ b/experiment/rotation-test/rotation_test.py
@@ -0,0 +1,141 @@
+import subprocess
+import re
+import os
+import numpy as np
+import matplotlib.pyplot as plt
+import math
+import fileinput
+import time
+import matplotlib.ticker as tickers
+import matplotlib as mpl
+from matplotlib.ticker import MultipleLocator,FormatStrFormatter,MaxNLocator
+import codecs
+import sys
+import threading
+
+# Schreibe input datei für Parameter
+def SetParametersCellProblem(ParameterSet, ParsetFilePath, outputPath):
+ print('----set Parameters -----')
+ with open(ParsetFilePath, 'r') as file:
+ filedata = file.read()
+ filedata = re.sub('(?m)^materialFunction\s?=.*','materialFunction = '+str(ParameterSet.materialFunction),filedata)
+ filedata = re.sub('(?m)^gamma\s?=.*','gamma='+str(ParameterSet.gamma),filedata)
+ filedata = re.sub('(?m)^numLevels\s?=.*','numLevels='+str(ParameterSet.numLevels)+' '+str(ParameterSet.numLevels) ,filedata)
+ filedata = re.sub('(?m)^outputPath\s?=\s?.*','outputPath='+str(outputPath),filedata)
+ f = open(ParsetFilePath,'w')
+ f.write(filedata)
+ f.close()
+
+
+# Ändere Parameter der MaterialFunction
+def SetParameterMaterialFunction(materialFunction, parameterName, parameterValue):
+ with open(Path+"/"+materialFunction+'.py', 'r') as file:
+ filedata = file.read()
+ filedata = re.sub('(?m)^'+str(parameterName)+'\s?=.*',str(parameterName)+' = '+str(parameterValue),filedata)
+ f = open(Path+"/"+materialFunction+'.py','w')
+ f.write(filedata)
+ f.close()
+
+# Rufe Programm zum Lösen des Cell-Problems auf
+def run_CellProblem(executable, parset,write_LOG):
+ print('----- RUN Cell-Problem ----')
+ processList = []
+ LOGFILE = "Cell-Problem_output.log"
+ print('LOGFILE:',LOGFILE)
+ print('executable:',executable)
+ if write_LOG:
+ p = subprocess.Popen(executable + parset
+ + " | tee " + LOGFILE, shell=True)
+
+ else:
+ p = subprocess.Popen(executable + parset
+ + " | tee " + LOGFILE, shell=True)
+
+ p.wait() # wait
+ processList.append(p)
+ exit_codes = [p.wait() for p in processList]
+
+ return
+
+# Read effective quantites
+def ReadEffectiveQuantities(QFilePath = os.path.dirname(os.getcwd()) + '/outputs/QMatrix.txt', BFilePath = os.path.dirname(os.getcwd())+ '/outputs/BMatrix.txt'):
+ # Read Output Matrices (effective quantities)
+ # From Cell-Problem output Files : ../outputs/Qmatrix.txt , ../outputs/Bmatrix.txt
+ # -- Read Matrix Qhom
+ X = []
+ # with codecs.open(path + '/outputs/QMatrix.txt', encoding='utf-8-sig') as f:
+ with codecs.open(QFilePath, encoding='utf-8-sig') as f:
+ for line in f:
+ s = line.split()
+ X.append([float(s[i]) for i in range(len(s))])
+ Q = np.array([[X[0][2], X[1][2], X[2][2]],
+ [X[3][2], X[4][2], X[5][2]],
+ [X[6][2], X[7][2], X[8][2]] ])
+
+ # -- Read Beff (as Vector)
+ X = []
+ # with codecs.open(path + '/outputs/BMatrix.txt', encoding='utf-8-sig') as f:
+ with codecs.open(BFilePath, encoding='utf-8-sig') as f:
+ for line in f:
+ s = line.split()
+ X.append([float(s[i]) for i in range(len(s))])
+ B = np.array([X[0][2], X[1][2], X[2][2]])
+ return Q, B
+
+# Function for evaluating the energy in terms of kappa, alpha and Q, B
+def eval_energy(kappa,alpha,Q,B) :
+ G=kappa*np.array([[np.cos(alpha)**2],[np.sin(alpha)**2],[np.sqrt(2)*np.cos(alpha)*np.sin(alpha)]])-B
+ return np.matmul(np.transpose(G),np.matmul(Q,G))[0,0]
+
+#-------------------------------------------------------------------------------------------------------
+########################
+#### SET PARAMETERS ####
+########################
+# ----- Setup Paths -----
+Path = "./experiment/rotation-test"
+# parset = ' ./experiment/wood-bilayer/cellsolver.parset.wood'
+ParsetFile = Path + '/cellsolver.parset'
+executable = 'build-cmake/src/Cell-Problem'
+write_LOG = True # writes Cell-Problem output-LOG in "Cell-Problem_output.log"
+
+# ---------------------------------
+# Setup Experiment
+# ---------------------------------
+outputPath = Path + '/results/'
+
+# ----- Define Input parameters --------------------
+class ParameterSet:
+ pass
+
+ParameterSet.materialFunction = "isotrop_orthotrop_rotation"
+ParameterSet.gamma=1.0
+ParameterSet.numLevels=4
+
+# ----- Define Parameters for Material Function --------------------
+# Liste mit Drehwinkel theta
+materialFunctionParameter=[0,np.pi/12, 2*np.pi/12,3*np.pi/12, 4*np.pi/12, 5*np.pi/12]
+
+# ------ Loops through Parameters for Material Function -----------
+for i in range(0,np.shape(materialFunctionParameter)[0]):
+ print("------------------")
+ print("New Loop")
+ print("------------------")
+ # Check output directory
+ path = outputPath + str(i)
+ isExist = os.path.exists(path)
+ if not isExist:
+ # Create a new directory because it does not exist
+ os.makedirs(path)
+ print("The new directory " + path + " is created!")
+ SetParameterMaterialFunction(ParameterSet.materialFunction, "param_theta",materialFunctionParameter[i])
+ SetParametersCellProblem(ParameterSet, ParsetFile, path)
+ #Run Cell-Problem
+ thread = threading.Thread(target=run_CellProblem(executable, " ./"+ParsetFile, write_LOG))
+ thread.start()
+ # ---------------------------------------------------
+ # wait here for the result to be available before continuing
+ thread.join()
+ f = open(path+"/parameter.txt", "w")
+ f.write("theta = "+str(materialFunctionParameter[i])+"\n")
+ f.close()
+ #
diff --git a/experiment/wood-bilayer/wood_european_beech.py b/experiment/wood-bilayer/wood_european_beech.py
index bcaf914a2ba135ed887444828e89fcd9026f2aba..456cb0bc32f45d9ee75544674dcace081cd5248f 100644
--- a/experiment/wood-bilayer/wood_european_beech.py
+++ b/experiment/wood-bilayer/wood_european_beech.py
@@ -138,11 +138,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]]
-# --- PHASE 3
+# --- PHASE 3 = Phase 1 gedreht
# y_1-direction: L
# y_2-direction: R
# x_3-direction: T
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)