diff --git a/microstructure_testsuite/plot.py b/JNLS_plotScripts/SN_energyLandscape_Fig10.py
similarity index 100%
rename from microstructure_testsuite/plot.py
rename to JNLS_plotScripts/SN_energyLandscape_Fig10.py
diff --git a/microstructure_testsuite/helper_functions.py b/microstructure_testsuite/deprecated/helper_functions.py
similarity index 100%
rename from microstructure_testsuite/helper_functions.py
rename to microstructure_testsuite/deprecated/helper_functions.py
diff --git a/microstructure_testsuite/microstructure_testsuite.py b/microstructure_testsuite/deprecated/microstructure_testsuite_DEPRECATED.py
similarity index 100%
rename from microstructure_testsuite/microstructure_testsuite.py
rename to microstructure_testsuite/deprecated/microstructure_testsuite_DEPRECATED.py
diff --git a/microstructure_testsuite/microstructure-testsuite.py b/microstructure_testsuite/microstructure-testsuite.py
new file mode 100644
index 0000000000000000000000000000000000000000..1f534078e117f8876eff46589c7a334dbdbac9e5
--- /dev/null
+++ b/microstructure_testsuite/microstructure-testsuite.py
@@ -0,0 +1,556 @@
+import os
+import sys
+import numpy as np
+import matplotlib as mpl
+import matplotlib.pyplot as plt
+import matplotlib.colors as colors
+from matplotlib.ticker import LogLocator
+import codecs
+import re
+import json
+import math
+import subprocess
+import fileinput
+from matplotlib.ticker import MultipleLocator,FormatStrFormatter,MaxNLocator
+from scipy.optimize import minimize_scalar
+import importlib
+
+# --------------------- Helper functions----------------------------------------------------------
+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
+
+
+
+ #--- Select specific experiment [x, y] with date from results_x/y
+def get_Q_B(basePath, index,perforation=False, perforatedLayer='upper'):
+ # results_index[0]/index[1]/...
+ #DataPath = './experiment/wood-bilayer_PLOS/results_' + str(data[0]) + '/' +str(data[1])
+ # DataPath = './results_' + str(index[0]) + '/' +str(index[1])
+ # DataPath = '.' + dirname + orientation + gridLevel + '/results_' + str(index[0]) + '/' +str(index[1])
+ # DataPath = dirname + orientation + gridLevel + '/results_' + str(index[0]) + '/' +str(index[1])
+ if perforation:
+ DataPath = basePath + '/results_' + perforatedLayer + '_' + str(index[0]) + '/' +str(index[1])
+ else :
+ DataPath = basePath + '/results_' + str(index[0]) + '/' +str(index[1])
+
+ QFilePath = DataPath + '/QMatrix.txt'
+ BFilePath = DataPath + '/BMatrix.txt'
+ # Read Q and B
+ Q, B = ReadEffectiveQuantities(QFilePath,BFilePath)
+ Q=0.5*(np.transpose(Q)+Q) # symmetrize
+ B=np.transpose([B])
+ return (Q,B)
+
+def get_local_minimizer_on_axes(Q,B):
+ invoke_function=lambda kappa: energy(kappa,0,Q,B)
+ result_0 = minimize_scalar(invoke_function, method="golden")
+ invoke_function=lambda kappa: energy(kappa,np.pi/2,Q,B)
+ result_90 = minimize_scalar(invoke_function, method="golden")
+ return np.array([[result_0.x,0],[result_90.x,np.pi/2]])
+
+
+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()
+
+
+# -------------------------- Matplotlib - parameters.--------------------------------------------------------------
+plt.style.use("seaborn")
+mpl.rcParams['text.usetex'] = True
+mpl.rcParams["font.family"] = "serif"
+mpl.rcParams["font.size"] = "8"
+mpl.rcParams['xtick.bottom'] = True
+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'] = 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
+mpl.rcParams.update({"axes.grid" : False}) # Add grid
+mpl.rcParams['axes.labelpad'] = 3
+mpl.rcParams['grid.linewidth'] = 0.25
+mpl.rcParams['grid.alpha'] = 0.9 # 0.75
+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.
+textwidth = 6.26894 #textwidth in inch
+width = textwidth * 0.5
+height = width/1.618 # The golden ratio.
+# ------------------------------------------------------------------------------------------------------------------
+
+
+# ------------------------------------------------------------------------------------------------------------------------------------------------------------------------
+# ------------------------------------------------------------------------------------------------------------------------------------------------------------------------
+computeEffectiveQuantities = True
+plotEnergyLandscape = True
+
+
+
+# 1. Experiment 2. ExperimentPathExtension 3. MicroGridLevels, 4. Executable 5.(additional) parameterArray: [["ParameterName", ParameterValue], ...]
+scenarios = [
+ ["cylindrical_1", "", [3,3], "/micro-problem" , [["mu_phase1", [90]]] ],
+ # ["cylindrical_1", "", [3,3], "/micro-problem" , [["mu_phase1", [90,100,110]]] ],
+ ["parametrized_laminate", "", [3,3], "/micro-problem" , [] ],
+ ]
+
+
+print('sys.argv[0]', sys.argv[0])
+print('sys.argv[1]', sys.argv[1])
+print('sys.argv[2]', sys.argv[2])
+
+
+CONTAINER = int(sys.argv[1])
+slurm_array_task_id = int(sys.argv[2])
+
+parameterFile = scenarios[slurm_array_task_id][0]
+pathExtension = scenarios[slurm_array_task_id][1]
+gridLevels = scenarios[slurm_array_task_id][2]
+executable = scenarios[slurm_array_task_id][3]
+parameterArray = scenarios[slurm_array_task_id][4]
+
+#Test
+# parameterArray = [[" -rho ", [4.0]]]
+
+
+print('pathExtension: ', pathExtension)
+
+
+
+#Get current working directory:
+# cwd = os.getcwd()
+modulePath = os.path.abspath(os.path.join(os.getcwd(), os.pardir))
+print("Parent Directory is: ", modulePath)
+
+
+#Path for parameterFile
+if CONTAINER:
+ #--- Barnard/CONTAINER version
+ pythonPath = "/dune/dune-microstructure/experiment" + pathExtension
+ resultPath = "outputs" + "_" + scenarios[slurm_array_task_id][0]
+ instrumentedPath = resultPath + "/instrumented"
+ executablePath = "/dune/dune-microstructure/build-cmake/src"
+else :
+ #--- Local version
+ pythonPath = modulePath + "/experiment" + pathExtension
+ resultPath = modulePath + "/outputs" + "_" + scenarios[slurm_array_task_id][0]
+ instrumentedPath = resultPath + "/instrumented"
+ executablePath = modulePath + '/build-cmake/src'
+
+try:
+ os.mkdir(resultPath)
+ os.mkdir(instrumentedPath)
+except OSError as error:
+ print(error)
+
+print('pythonPath: ', pythonPath)
+
+
+### Access ParameterFile variables.
+sys.path.insert(0, pythonPath)
+__import__(parameterFile, fromlist=['parameterSet'])
+imported_module = importlib.import_module(parameterFile)
+parameterSet = getattr(imported_module, "parameterSet")
+
+executable = executablePath + executable
+
+resultPathBase = resultPath
+
+#----------------------------------------------------------------------------------------------------------------------
+if computeEffectiveQuantities:
+ print('Computing effective quantities... ')
+ for i in range(0, len(parameterArray)):
+ for v in range(0,len(parameterArray[i][1])):
+ print("------------------")
+ print("New Loop")
+ print("------------------")
+ print('i:', i)
+ print('v:', v)
+ print('len(parameterArray[i][1]):', len(parameterArray[i][1]))
+ # print('parameterArray[i][v]:', parameterArray[i][v])
+
+ print('parameterName:' , parameterArray[i][0])
+ print('parameterValue: ', parameterArray[i][1][v])
+ # print('np.shape(parameterArray)[1]:', np.shape(parameterArray)[1])
+ # print('np.shape(parameterArray[i][1]):', np.shape(parameterArray[i][1]))
+ # print('np.shape(parameterArray[i][1])[1]:', np.shape(parameterArray[i][1])[1])
+ parameterName = parameterArray[i][0]
+ parameterValue = parameterArray[i][1][v]
+
+
+ resultPath = resultPathBase + "/" + parameterName
+ try:
+ os.mkdir(resultPath)
+ except OSError as error:
+ print(error)
+
+ # resultPath = resultPath + "/" + str
+
+
+
+ ## Alternative version: actually changing the parset file..
+ # Change Parameters:
+ print('...change input parameter:' , parameterName)
+ SetParameterMaterialFunction(pythonPath + "/" + scenarios[slurm_array_task_id][0], parameterName , parameterValue)
+
+
+ baseName = scenarios[slurm_array_task_id][0] + "_" + parameterArray[i][0] + "_" + str(parameterArray[i][1][v])
+
+ #--- Compute
+ processList = []
+ for microGridLevel in range(gridLevels[0],gridLevels[1]+1):
+ # if conforming_DiscreteJacobian:
+ # conformity = "_conforming"
+ # else:
+ # conformity = "_nonconforming"
+
+ # LOGFILE = resultPath + "/" + parameterFile + conformity + "_macroGridLevel" + str(macroGridLevel) + parameterName + ".log"
+ LOGFILE = resultPath + "/" + parameterFile + "_microGridLevel" + str(microGridLevel) + "_" + parameterName + "_" + str(parameterValue) + ".log"
+
+ print('LOGFILE:', LOGFILE)
+ print('executable:', executable)
+ print('parameterFile:', parameterFile)
+ print('resultPath:', resultPath)
+
+ # testArray = [ " -rho " + str(8.0), " -beta " + str(0.10) ]
+
+ # start_time = time.time()
+ p = subprocess.Popen(executable + " " + pythonPath + " " + parameterFile
+ + " -macroGridLevel " + str(microGridLevel)
+ + " -resultPath " + str(resultPath)
+ + " -baseName " + str(baseName)
+ # + " -" + parameterName + " " + str(parameterValue)
+ + " | tee " + LOGFILE, shell=True)
+
+ p.wait() # wait
+ # print("--- %s seconds ---" % (time.time() - start_time))
+ print('------FINISHED PROGRAM on macroGridLevel:' + str(microGridLevel))
+ processList.append(p)
+
+ # Wait for all simulation subprocesses before proceeding to the error measurement step
+ exit_codes = [p.wait() for p in processList]
+
+
+# Set Path for effective quantitites.
+QFilePath = resultPath + "/QMatrix.txt"
+BFilePath = resultPath + "/BMatrix.txt"
+
+
+
+#----------------------------------------------------------------------------------------------------------------------
+if plotEnergyLandscape:
+ print('Plot Energy landscape ...')
+ for i in range(0, len(parameterArray)):
+ for v in range(0,len(parameterArray[i][1])):
+
+
+ fig, (ax1) = plt.subplots( nrows=1,figsize=(3,3.5),subplot_kw=dict(projection='polar'), gridspec_kw={'hspace': 0.4})
+
+ parameterName = parameterArray[i][0]
+ parameterValue = parameterArray[i][1][v]
+
+ resultPath = resultPathBase + "/" + parameterName
+ try:
+ os.mkdir(resultPath)
+ except OSError as error:
+ print(error)
+
+ # Set Path for effective quantitites.
+ QFilePath = resultPath + "/QMatrix.txt"
+ BFilePath = resultPath + "/BMatrix.txt"
+
+ # Read the effective quantities from .txt-File
+ Q, B = ReadEffectiveQuantities(QFilePath,BFilePath)
+ print('Qhom: \n', Q)
+ print('Beff: \n', B)
+
+ # Compute local and global minimizer (by Sampling method)
+ kappa=0
+ kappa_pos=0
+ kappa_neg=0
+ #
+ 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]
+ # 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]
+ #
+
+ print('E.min(): ',E.min())
+ print('E.max(): ',E.max())
+ E=E-E.min()
+
+ # Emax=1000
+ Emax = E.max()
+ levs=np.geomspace(1,Emax,200)
+
+
+ pcm=ax1.contourf(theta, r, E, levs, norm=colors.PowerNorm(gamma=0.4), cmap='brg', zorder=0)
+ # ax1.set_title(r"ratio $r=0.12$",pad=20)
+ ax1.set_title(r"Title",pad=20)
+
+
+ print('kappamin_pos: ', kappamin_pos)
+ print('alphamin_pos: ', alphamin_pos)
+ print('kappamin_neg: ', kappamin_neg)
+ print('alphamin_neg: ', alphamin_neg)
+
+ ax1.plot(alphamin_neg, kappamin_neg,
+ markerfacecolor='green',
+ markeredgecolor='white', # marker edgecolor
+ marker='D', # each marker will be rendered as a circle
+ markersize=6, # marker size
+ markeredgewidth=1, # marker edge width
+ linewidth=0, # line width
+ zorder=3,
+ alpha=1)
+
+ ax1.plot(alphamin, kappamin,
+ markerfacecolor='blue',
+ markeredgecolor='white', # marker edgecolor
+ marker='s', # each marker will be rendered as a circle
+ markersize=6, # marker size
+ markeredgewidth=1, # marker edge width
+ linewidth=0,
+ zorder=3,
+ alpha=1, # Change opacity
+ label = r"$\kappa_{0^\circ}(m^{-1})$")
+
+
+ print('alphamin: ' , alphamin )
+ print('kappamin: ' , kappamin )
+ print('Global Minimizer (alphamin, kappamin):' , '('+ str(np.round(np.rad2deg(alphamin),decimals=1)) + ',' + str(np.round(kappamin,decimals=1))+ ')' )
+ print('Local Minimizer (alphamin, kappamin):' , '('+ str(np.round(360-np.rad2deg(alphamin_neg),decimals=1)) + ',' + str(np.round(-kappamin_neg,decimals=1))+ ')' )
+ # ax1.annotate(alphamin, (alphamin,kappamin))
+
+ #--- annotate global minimizer
+ ax1.annotate( '('+ str(np.round(np.rad2deg(alphamin),decimals=1)) + '°' + ',' + str(np.round(kappamin,decimals=1))+ ')' ,
+ # ax1.annotate( r'$(\theta_{m}, \kappa_m)=$' + '('+ str(np.round(np.rad2deg(alphamin),decimals=1)) + '°' + ',' + str(np.round(kappamin,decimals=1))+ ')' ,
+ xy = (alphamin,kappamin),
+ xytext=(0.10, np.round(kappamin,decimals=1)-0.25),
+ color='ivory',
+ fontsize='8')
+
+
+ #--- annotate local minimizer
+ ax1.annotate( '('+ str(np.round(360-np.rad2deg(alphamin_neg),decimals=1)) + '°' + ',' + str(np.round(-kappamin_neg,decimals=1))+ ')' ,
+ xy = (alphamin_neg,kappamin_neg),
+ xytext=(-1.6, np.round(kappamin_neg,decimals=1)+0.9),
+ color='ivory',
+ # zorder=5,
+ fontsize='8')
+
+ # rotate radius axis labels.
+ ax1.set_rlabel_position(135)
+
+ colorbarticks=np.geomspace(0.00001,Emax,10)
+ cbar = plt.colorbar(pcm, ax=[ax1], extend='max', ticks=colorbarticks, pad=0.1, orientation="horizontal")
+ cbar.ax.tick_params(labelsize=8)
+
+ anglelabel=["0°","45°", "90°", "135°","180°","135°","90°","45°"]
+ ax1.set_xticks(np.array([.0,1/4,2/4,3/4,1,5/4,6/4,7/4])*np.pi)
+ ax1.set_xticklabels(anglelabel)
+ ax1.xaxis.grid(True, color='white')
+ ax1.set_yticks(np.array([1,2,3,4]))
+ ax1.set_yticklabels(["1","2","3","4"], zorder=20, color="white")
+ ax1.yaxis.set_tick_params(color='white')
+ ax1.yaxis.grid(True, color='white', alpha=0.75)
+ # Save Figure as pdf.
+ fig.savefig(resultPath + '/energy_landscape.png', dpi=300)
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+# ----------------------------------------------------------------------------------------
+# ------------------------------ EXPERIMENT 1 --------------------------------------------
+# ----------------------------------------------------------------------------------------
+# ----------------------------------------------------------------------------------------
+### PARAMETERS FROM SIMULATION
+
+# [r, h, omega_flat, omega_target, theta, Kappa_experimental]
+# 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)
+
+# materialFunctionParameter=[
+# [ # Dataset Ratio r = 0.12
+# [0.12, 0.0047, 17.32986047, 14.70179844, 0, 1.140351217],
+# [0.12, 0.0047, 17.32986047, 13.6246, 0, 1.691038688],
+# [0.12, 0.0047, 17.32986047, 12.42994508, 0, 2.243918105],
+# [0.12, 0.0047, 17.32986047, 11.69773413, 0, 2.595732726],
+# [0.12, 0.0047, 17.32986047, 11.14159987, 0, 2.945361006],
+# [0.12, 0.0047, 17.32986047, 9.500670278, 0, 4.001528043],
+# [0.12, 0.0047, 17.32986047, 9.005046347, 0, 4.312080261]
+# ],
+# [ # Dataset Ratio r = 0.17
+# [0.17, 0.0049, 17.28772791 , 14.75453569, 0, 1.02915975],
+# [0.17, 0.0049, 17.28772791 , 13.71227639, 0, 1.573720805],
+# [0.17, 0.0049, 17.28772791 , 12.54975012, 0, 2.407706364],
+# [0.17, 0.0049, 17.28772791 , 11.83455959, 0, 2.790518802],
+# [0.17, 0.0049, 17.28772791 , 11.29089521, 0, 3.173814476],
+# [0.17, 0.0049, 17.28772791 , 9.620608917, 0, 4.187433094],
+# [0.17, 0.0049, 17.28772791 , 9.101671742, 0, 4.511739072]
+# ],
+# [ # Dataset Ratio r = 0.22
+# [0.22, 0.0053, 17.17547062, 14.72680026, 0, 1.058078122],
+# [0.22, 0.0053, 17.17547062, 13.64338887, 0, 1.544624544],
+# [0.22, 0.0053, 17.17547062, 12.41305478, 0, 2.317033799],
+# [0.22, 0.0053, 17.17547062, 11.66482931, 0, 2.686043143],
+# [0.22, 0.0053, 17.17547062, 11.09781471, 0, 2.967694189],
+# [0.22, 0.0053, 17.17547062, 9.435795985, 0, 3.913528418],
+# [0.22, 0.0053, 17.17547062, 8.959564147, 0, 4.262750825]
+# ],
+# [ # Dataset Ratio r = 0.34
+# [0.34, 0.0063, 17.14061081 , 14.98380876, 0, 0.789078472],
+# [0.34, 0.0063, 17.14061081 , 13.97154915, 0, 1.1299263],
+# [0.34, 0.0063, 17.14061081 , 12.77309253, 0, 1.738136936],
+# [0.34, 0.0063, 17.14061081 , 12.00959929, 0, 2.159520896],
+# [0.34, 0.0063, 17.14061081 , 11.42001731, 0, 2.370047499],
+# [0.34, 0.0063, 17.14061081 , 9.561447179, 0, 3.088299431],
+# [0.34, 0.0063, 17.14061081 , 8.964704969, 0, 3.18097558]
+# ],
+# [ # Dataset Ratio r = 0.43
+# [0.43, 0.0073, 17.07559686 , 15.11316339, 0, 0.577989364],
+# [0.43, 0.0073, 17.07559686 , 14.17997082, 0, 0.829007544],
+# [0.43, 0.0073, 17.07559686 , 13.05739844, 0, 1.094211707],
+# [0.43, 0.0073, 17.07559686 , 12.32309209, 0, 1.325332511],
+# [0.43, 0.0073, 17.07559686 , 11.74608518, 0, 1.400455154],
+# [0.43, 0.0073, 17.07559686 , 9.812372466, 0, 1.832325697],
+# [0.43, 0.0073, 17.07559686 , 9.10519385 , 0, 2.047483977]
+# ],
+# [ # Dataset Ratio r = 0.49
+# [0.49, 0.008, 17.01520754, 15.30614414, 0, 0.357615902],
+# [0.49, 0.008, 17.01520754, 14.49463867, 0, 0.376287785],
+# [0.49, 0.008, 17.01520754, 13.46629742, 0, 0.851008627],
+# [0.49, 0.008, 17.01520754, 12.78388234, 0, 0.904475291],
+# [0.49, 0.008, 17.01520754, 12.23057715, 0, 1.039744708],
+# [0.49, 0.008, 17.01520754, 10.21852839, 0, 1.346405241],
+# [0.49, 0.008, 17.01520754, 9.341730605, 0, 1.566568558]
+# ]
+# ]
+
+
+# n=len(materialFunctionParameter)
+# m=len(materialFunctionParameter[0])
+# kappa_0=np.zeros([n,m])
+# energy_0=np.zeros([n,m])
+# kappa_90=np.zeros([n,m])
+# energy_90=np.zeros([n,m])
+# energy_global=np.zeros([n,m])
+# kappa_global_min=np.zeros([n,m])
+# energy_origin=np.zeros([n,m])
+# for i in range(0,n):
+# for j in range(0,m):
+# Q, B=get_Q_B(basePath, [i,j])
+# minimizers=get_local_minimizer_on_axes(Q,B)
+# kappa_0[i,j]=minimizers[0,0]
+# energy_0[i,j]=energy(kappa_0[i,j],0,Q,B)
+# kappa_90[i,j]=minimizers[1,0]
+# energy_90[i,j]=energy(kappa_90[i,j],np.pi/2,Q,B)
+# if energy_0[i,j]<energy_90[i,j]:
+# kappa_global_min[i,j]=kappa_0[i,j]
+# energy_global[i,j]=energy_0[i,j]
+# else:
+# kappa_global_min[i,j]=kappa_90[i,j]
+# energy_global[i,j]=energy_90[i,j]
+# energy_origin[i,j]=energy(0,0,Q,B)
+
+
+# # Error plot
+# for i in range(0,n):
+# plt.plot(np.array(materialFunctionParameter[i])[:,3], np.abs(kappa_0[i,:]-np.array(materialFunctionParameter[i])[:,5])/np.array(materialFunctionParameter[i])[:,5], label=str(i))
+
+# plt.legend()
+
+
+# # -------------------------------------------------------------------------------
+# ### Plot diagrams: comparison with experimental data
+
+# print('Create plots for comparison with experimental data ...')
+