diff --git a/Plot-Scripts/WoodBilayer_ExperimentComparison.py b/Plot-Scripts/WoodBilayer_ExperimentComparison.py
index c32095861010dfbc81d1a5806baf6482718df822..c632ce91e8d0272e6d090f279e5d396091294671 100644
--- a/Plot-Scripts/WoodBilayer_ExperimentComparison.py
+++ b/Plot-Scripts/WoodBilayer_ExperimentComparison.py
@@ -233,7 +233,7 @@ for dataset_number in dataset_numbers:
# plt.xscale('log') # Use Logarithmic-Scale
# plt.yscale('log')
ax.set_xlabel(r"Wood moisture content $\omega (\%)$")
- ax.set_ylabel(r"Curvature $\kappa$")
+ ax.set_ylabel(r"Curvature $\kappa(m^{-1})$")
# plt.tight_layout()
# # --- Set Line labels
diff --git a/Plot-Scripts/WoodBilayer_ExperimentComparisonError_localminimizer.py b/Plot-Scripts/WoodBilayer_ExperimentComparisonError_localminimizer.py
index 11750cb1e69648d1cb66cd5ced1ee30f3e12ce4a..9c839ec2ed8c435c896e70326f167a1e3886ed84 100644
--- a/Plot-Scripts/WoodBilayer_ExperimentComparisonError_localminimizer.py
+++ b/Plot-Scripts/WoodBilayer_ExperimentComparisonError_localminimizer.py
@@ -13,23 +13,35 @@ 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["font.size"] = "8"
mpl.rcParams['xtick.bottom'] = True
-mpl.rcParams['xtick.major.size'] = 3
+mpl.rcParams['xtick.major.size'] = 2
mpl.rcParams['xtick.minor.size'] = 1.5
mpl.rcParams['xtick.major.width'] = 0.75
+mpl.rcParams['xtick.labelsize'] = 8
+mpl.rcParams['xtick.major.pad'] = 1
+
mpl.rcParams['ytick.left'] = True
-mpl.rcParams['ytick.major.size'] = 3
+mpl.rcParams['ytick.major.size'] = 2
mpl.rcParams['ytick.minor.size'] = 1.5
mpl.rcParams['ytick.major.width'] = 0.75
+mpl.rcParams['ytick.labelsize'] = 8
+mpl.rcParams['ytick.major.pad'] = 1
+
+mpl.rcParams['axes.titlesize'] = 8
+mpl.rcParams['axes.titlepad'] = 1
+mpl.rcParams['axes.labelsize'] = 8
+
+#Adjust Legend:
+mpl.rcParams['legend.frameon'] = True # Use frame for legend
+# mpl.rcParams['legend.framealpha'] = 0.5
+mpl.rcParams['legend.fontsize'] = 8 # fontsize of legend
+
#Adjust grid:
mpl.rcParams.update({"axes.grid" : True}) # Add grid
@@ -40,8 +52,12 @@ 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.
+# width = 5.79
+# height = width / 1.618 # The golden ratio.
+textwidth = 6.26894 #textwidth in inch
+width = textwidth * 0.5
+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)
@@ -65,24 +81,24 @@ data = [ # Dataset Ratio r = 0.49
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
+ markersize=3.5, # marker size
# markerfacecolor='darkorange', # marker facecolor
markeredgecolor='black', # marker edgecolor
- markeredgewidth=0.75, # marker edge width
+ markeredgewidth=0.5, # marker edge width
# linestyle='dashdot', # line style will be dash line
- linewidth=1.5, # line width
+ linewidth=1, # 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
+ markersize=3.5, # marker size
# markerfacecolor='cornflowerblue', # marker facecolor
markeredgecolor='black', # marker edgecolor
- markeredgewidth=0.75, # marker edge width
+ markeredgewidth=0.5, # marker edge width
# linestyle='--', # line style will be dash line
- linewidth=1.5, # line width
+ linewidth=1, # line width
zorder=3,
alpha=0.8, # Change opacity
label = r"$\kappa_{2,sim}$")
@@ -90,12 +106,12 @@ line_2 = ax.plot(np.array(data[0]), np.array(data[2]), # data
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
+ markersize=3.5, # marker size
# markerfacecolor='cornflowerblue', # marker facecolor
markeredgecolor='black', # marker edgecolor
- markeredgewidth=0.75, # marker edge width
+ markeredgewidth=0.5, # marker edge width
# linestyle='--', # line style will be dash line
- linewidth=1.5, # line width
+ linewidth=1, # line width
zorder=3,
alpha=0.8, # Change opacity
label = r"$\kappa_{exp}$")
@@ -118,7 +134,7 @@ 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)
+ax.set_ylabel(r"Curvature $\kappa$($m^{-1}$)", labelpad=4)
plt.tight_layout()
# # --- Set Line labels
@@ -139,6 +155,7 @@ legend = ax.legend()
frame = legend.get_frame()
frame.set_edgecolor('black')
+frame.set_linewidth(0.5)
# --- Adjust left/right spacing:
diff --git a/Plot-Scripts/WoodBilayer_Experiment_EnergyLandscape_localmin.py b/Plot-Scripts/WoodBilayer_Experiment_EnergyLandscape_localmin.py
new file mode 100644
index 0000000000000000000000000000000000000000..e42b5a04b8b0f5be3be5fb5f980e64a788c99897
--- /dev/null
+++ b/Plot-Scripts/WoodBilayer_Experiment_EnergyLandscape_localmin.py
@@ -0,0 +1,148 @@
+#!/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
+from matplotlib.ticker import LogLocator
+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
+
+#--- Select specific experiment [x, y] with date from results_x/y
+data=[5,6]
+DataPath = './experiment/wood-bilayer_PLOS/results_' + str(data[0]) + '/' +str(data[1])
+#DataPath = './results_' + str(data[0]) + '/' +str(data[1])
+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])
+energyscalingfactor = thickness**2
+# Read Q and B
+Q, B = ReadEffectiveQuantities(QFilePath,BFilePath)
+Q=0.5*(np.transpose(Q)+Q) # symmetrize
+B=np.transpose([B])
+#
+# Compute lokal and global minimizer
+kappa=0
+kappa_pos=0
+kappa_neg=0
+#
+N=500
+length=4
+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) * energyscalingfactor
+ else:
+ E[i,j]=energy(-r[i,j],theta[i,j],Q,B) * energyscalingfactor
+#
+# Compute Minimizer
+[imin,jmin]=np.unravel_index(E.argmin(),(N,N))
+kappamin=r[imin,jmin]
+alphamin=theta[imin,jmin]
+# 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]
+# 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]
+#
+E=E/E.min()
+print(Emin_pos/Emin_neg)
+fig, ax = plt.subplots(figsize=(6,6),subplot_kw=dict(projection='polar'))
+levs=np.geomspace(1,E.max(),1000)
+pcm=ax.contourf(theta, r, E, levs, norm=colors.PowerNorm(gamma=0.4), cmap='brg')
+ax.set_xticks(np.array([.0,1/4,2/4,3/4,1,5/4,6/4,7/4])*np.pi)
+anglelabel=["0°","45°", "90°", "135°","180°","135°","90°","45°"]
+ax.set_xticklabels(anglelabel)
+ax.set_yticks([1,2,3,4])
+#ax.set_yticklabels(["1$m^{-1}$","2$m^{-1}$","3$m^{-1}$","4$m^{-1}$"])
+#
+ax.plot([alphamin_pos,alphamin_pos+np.pi], [kappamin_pos,kappamin_pos],
+ markerfacecolor='red',
+ markeredgecolor='black', # marker edgecolor
+ marker='s', # each marker will be rendered as a circle
+ markersize=5, # marker size
+ markeredgewidth=0.5, # marker edge width
+ linewidth=0,
+ zorder=3,
+ alpha=1, # Change opacity
+ label = r"$\kappa_{2,sim}(m^{-1})$")
+ax.plot(alphamin_neg, kappamin_neg,
+ markerfacecolor='blue',
+ markeredgecolor='black', # marker edgecolor
+ marker='D', # each marker will be rendered as a circle
+ markersize=5, # marker size
+ markeredgewidth=0.5, # marker edge width
+ linewidth=0, # line width
+ zorder=3,
+ alpha=1, # Change opacity
+ label = r"$\kappa_{1,sim}(m^{-1})$")
+colorbarticks=np.linspace(1,15,10)
+cbar = plt.colorbar(pcm, extend='max', ticks=colorbarticks, pad=0.1)
+#bounds = ['0','1/80','1/20','1/5','4/5']
+cbar.ax.tick_params(labelsize=8)
+#cbar.set_ticklabels(bounds)
+fig.legend(loc="upper left")
+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('./wood-bilayer_PLOS_dataset_' +str(data[0]) + '_exp' + str(data[1]) + '.pdf', dpi=300)
\ No newline at end of file
diff --git a/experiment/wood-bilayer_PLOS/Auswertung.py b/experiment/wood-bilayer_PLOS/Auswertung.py
deleted file mode 100644
index 9bd21562b7099b48f9b71e9fc1aa17bf0d884f96..0000000000000000000000000000000000000000
--- a/experiment/wood-bilayer_PLOS/Auswertung.py
+++ /dev/null
@@ -1,136 +0,0 @@
-#!/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
index 7c3feaac02477d6f6ee29abfb36937587aa94aa6..6379e849fdaf9814bb1bfe56dd0bfd2d2c8f34b5 100644
--- a/experiment/wood-bilayer_PLOS/PolarPlotLocalEnergy.py
+++ b/experiment/wood-bilayer_PLOS/PolarPlotLocalEnergy.py
@@ -111,7 +111,7 @@ for dataset_number in dataset_numbers:
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])
+ ax.set_yticks([kappamin_pos,kappamin_neg])
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)