Skip to content
Snippets Groups Projects
Commit bf9c6df3 authored by Klaus Böhnlein's avatar Klaus Böhnlein
Browse files

add energy scaling to PlotLocalEnergy

parent 8318c0b1
No related branches found
No related tags found
No related merge requests found
...@@ -9,7 +9,7 @@ import numpy as np ...@@ -9,7 +9,7 @@ import numpy as np
import matplotlib.pyplot as plt import matplotlib.pyplot as plt
import matplotlib.colors as colors import matplotlib.colors as colors
import codecs import codecs
import re
def energy(kappa,alpha,Q,B) : def energy(kappa,alpha,Q,B) :
...@@ -49,7 +49,10 @@ def ReadEffectiveQuantities(QFilePath, BFilePath): ...@@ -49,7 +49,10 @@ def ReadEffectiveQuantities(QFilePath, BFilePath):
return Q, B return Q, B
# Number of experiments / folders # Number of experiments / folders
number=4 number=7
show_plot = False
kappa=np.zeros(number) kappa=np.zeros(number)
for n in range(0,number): for n in range(0,number):
# Read from Date # Read from Date
...@@ -58,21 +61,28 @@ for n in range(0,number): ...@@ -58,21 +61,28 @@ for n in range(0,number):
DataPath = './experiment/perforated-bilayer/results/'+str(n) DataPath = './experiment/perforated-bilayer/results/'+str(n)
QFilePath = DataPath + '/QMatrix.txt' QFilePath = DataPath + '/QMatrix.txt'
BFilePath = DataPath + '/BMatrix.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, B = ReadEffectiveQuantities(QFilePath,BFilePath)
# Q=0.5*(np.transpose(Q)+Q) # symmetrize # Q=0.5*(np.transpose(Q)+Q) # symmetrize
B=np.transpose([B]) B=np.transpose([B])
# #
N=300 N=500
length=12 length=12
r, theta = np.meshgrid(np.linspace(0,length,N),np.radians(np.linspace(0, 360, N))) r, theta = np.meshgrid(np.linspace(0,length,N),np.radians(np.linspace(0, 360, N)))
E=np.zeros(np.shape(r)) E=np.zeros(np.shape(r))
for i in range(0,N): for i in range(0,N):
for j in range(0,N): for j in range(0,N):
if theta[i,j]<np.pi: if theta[i,j]<np.pi:
E[i,j]=energy(r[i,j],theta[i,j],Q,B) E[i,j]=energy(r[i,j],theta[i,j],Q,B) * (thickness**2)
else: else:
E[i,j]=energy(-r[i,j],theta[i,j],Q,B) E[i,j]=energy(-r[i,j],theta[i,j],Q,B) * (thickness**2)
# Compute Minimizer # Compute Minimizer
[imin,jmin]=np.unravel_index(E.argmin(),(N,N)) [imin,jmin]=np.unravel_index(E.argmin(),(N,N))
...@@ -85,8 +95,16 @@ for n in range(0,number): ...@@ -85,8 +95,16 @@ for n in range(0,number):
ax.set_xticks([0,np.pi/2]) ax.set_xticks([0,np.pi/2])
ax.set_yticks([kappamin]) ax.set_yticks([kappamin])
colorbarticks=np.linspace(E.min(),E.max(),6) colorbarticks=np.linspace(E.min(),E.max(),6)
plt.colorbar(pcm, extend='max', ticks=colorbarticks, pad=0.1) cbar = plt.colorbar(pcm, extend='max', ticks=colorbarticks, pad=0.1)
plt.show() 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('Plot_PerforatedBilayer.pdf')
f = open("./experiment/perforated-bilayer/results/kappa_simulation.txt", "w") f = open("./experiment/perforated-bilayer/results/kappa_simulation.txt", "w")
f.write(str(kappa)) f.write(str(kappa))
......
...@@ -107,8 +107,8 @@ def eval_energy(kappa,alpha,Q,B) : ...@@ -107,8 +107,8 @@ def eval_energy(kappa,alpha,Q,B) :
# pythonPath = '/home/klaus/Desktop/Dune_release/dune-microstructure/experiment/wood-bilayer' # pythonPath = '/home/klaus/Desktop/Dune_release/dune-microstructure/experiment/wood-bilayer'
path = os.getcwd() + '/experiment/perforated-bilayer/results/' path = os.getcwd() + '/experiment/perforated-bilayer/results/'
pythonPath = os.getcwd() + '/experiment/perforated-bilayer' pythonPath = os.getcwd() + '/experiment/perforated-bilayer'
# pythonModule = "perforated_wood_upper" pythonModule = "perforated_wood_upper"
pythonModule = "perforated_wood_lower" # pythonModule = "perforated_wood_lower"
executable = os.getcwd() + '/build-cmake/src/Cell-Problem' executable = os.getcwd() + '/build-cmake/src/Cell-Problem'
# --------------------------------- # ---------------------------------
# Setup Experiment # Setup Experiment
...@@ -127,11 +127,12 @@ gamma = 1.0 ...@@ -127,11 +127,12 @@ gamma = 1.0
#Experiment: Perforate "active" bilayer phase #Experiment: Perforate "active" bilayer phase
materialFunctionParameter=[ materialFunctionParameter=[
[0.22, 0.0053, 17.17547062, 8.959564147, 0.0, 0.0 ], [0.22, 0.0053, 17.17547062, 8.959564147, 0.0, 0.0 ],
[0.22, 0.0053, 17.17547062, 8.959564147, 0.0, 0.05 ],
[0.22, 0.0053, 17.17547062, 8.959564147, 0.0, 0.1 ], [0.22, 0.0053, 17.17547062, 8.959564147, 0.0, 0.1 ],
[0.22, 0.0053, 17.17547062, 8.959564147, 0.0, 0.15 ],
[0.22, 0.0053, 17.17547062, 8.959564147, 0.0, 0.2 ], [0.22, 0.0053, 17.17547062, 8.959564147, 0.0, 0.2 ],
[0.22, 0.0053, 17.17547062, 8.959564147, 0.0, 0.25 ],
[0.22, 0.0053, 17.17547062, 8.959564147, 0.0, 0.3 ] [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 ----------- # ------ Loops through Parameters for Material Function -----------
......
...@@ -9,7 +9,7 @@ import numpy as np ...@@ -9,7 +9,7 @@ import numpy as np
import matplotlib.pyplot as plt import matplotlib.pyplot as plt
import matplotlib.colors as colors import matplotlib.colors as colors
import codecs import codecs
import re
def energy(kappa,alpha,Q,B) : def energy(kappa,alpha,Q,B) :
...@@ -55,7 +55,8 @@ def ReadEffectiveQuantities(QFilePath, BFilePath): ...@@ -55,7 +55,8 @@ def ReadEffectiveQuantities(QFilePath, BFilePath):
number=7 number=7
show_plot = False show_plot = False
dataset_numbers = [0, 1, 2, 3, 4, 5] # dataset_numbers = [0, 1, 2, 3, 4, 5]
dataset_numbers = [0]
for dataset_number in dataset_numbers: for dataset_number in dataset_numbers:
...@@ -66,6 +67,13 @@ for dataset_number in dataset_numbers: ...@@ -66,6 +67,13 @@ for dataset_number in dataset_numbers:
DataPath = './experiment/wood-bilayer/results_'+ str(dataset_number) + '/' +str(n) DataPath = './experiment/wood-bilayer/results_'+ str(dataset_number) + '/' +str(n)
QFilePath = DataPath + '/QMatrix.txt' QFilePath = DataPath + '/QMatrix.txt'
BFilePath = DataPath + '/BMatrix.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, B = ReadEffectiveQuantities(QFilePath,BFilePath)
# Q=0.5*(np.transpose(Q)+Q) # symmetrize # Q=0.5*(np.transpose(Q)+Q) # symmetrize
B=np.transpose([B]) B=np.transpose([B])
...@@ -78,9 +86,9 @@ for dataset_number in dataset_numbers: ...@@ -78,9 +86,9 @@ for dataset_number in dataset_numbers:
for i in range(0,N): for i in range(0,N):
for j in range(0,N): for j in range(0,N):
if theta[i,j]<np.pi: if theta[i,j]<np.pi:
E[i,j]=energy(r[i,j],theta[i,j],Q,B) E[i,j]=energy(r[i,j],theta[i,j],Q,B) * (thickness**2)
else: else:
E[i,j]=energy(-r[i,j],theta[i,j],Q,B) E[i,j]=energy(-r[i,j],theta[i,j],Q,B) * (thickness**2)
# Compute Minimizer # Compute Minimizer
[imin,jmin]=np.unravel_index(E.argmin(),(N,N)) [imin,jmin]=np.unravel_index(E.argmin(),(N,N))
...@@ -93,11 +101,12 @@ for dataset_number in dataset_numbers: ...@@ -93,11 +101,12 @@ for dataset_number in dataset_numbers:
ax.set_xticks([0,np.pi/2]) ax.set_xticks([0,np.pi/2])
ax.set_yticks([kappamin]) ax.set_yticks([kappamin])
colorbarticks=np.linspace(E.min(),E.max(),6) colorbarticks=np.linspace(E.min(),E.max(),6)
plt.colorbar(pcm, extend='max', ticks=colorbarticks, pad=0.1) cbar = plt.colorbar(pcm, extend='max', ticks=colorbarticks, pad=0.1)
cbar.ax.tick_params(labelsize=6)
if (show_plot): if (show_plot):
plt.show() plt.show()
# Save Figure as .pdf # Save Figure as .pdf
width = 5.79 width = 5.79
height = width / 1.618 # The golden ratio. height = width / 1.618 # The golden ratio.
fig.set_size_inches(width, height) fig.set_size_inches(width, height)
fig.savefig('Plot_dataset_' +str(dataset_number) + '_exp' +str(n) + '.pdf') fig.savefig('Plot_dataset_' +str(dataset_number) + '_exp' +str(n) + '.pdf')
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment