diff --git a/src/PhaseDiagram_Contour.py b/src/PhaseDiagram_Contour.py
new file mode 100644
index 0000000000000000000000000000000000000000..98138fc9fefe8f55ec3be1b808767a0e67cc54b6
--- /dev/null
+++ b/src/PhaseDiagram_Contour.py
@@ -0,0 +1,481 @@
+import numpy as np
+import matplotlib.pyplot as plt
+import sympy as sym
+import math
+import os
+import subprocess
+import fileinput
+import re
+import matlab.engine
+import sys
+from ClassifyMin import *
+from HelperFunctions import *
+# from CellScript import *
+from mpl_toolkits.mplot3d import Axes3D
+import matplotlib.cm as cm
+from vtk.util import numpy_support
+from pyevtk.hl import gridToVTK
+from matplotlib.ticker import MultipleLocator,FormatStrFormatter,MaxNLocator
+
+import time
+# print(sys.executable)
+
+# --------------------------------------------------------------------
+# START :
+# INPUT (Parameters): alpha, beta, theta, gamma, mu1, rho1
+#
+# -Option 1 : (Case lambda = 0 => q12 = 0)
+# compute q1,q2,b1,b2 from Formula
+# Option 1.1 :
+# set mu_gamma = 'q1' or 'q2' (extreme regimes: gamma \in {0,\infty})
+# Option 1.2 :
+# compute mu_gamma with 'Compute_MuGamma' (2D problem much faster then Cell-Problem)
+# -Option 2 :
+# compute Q_hom & B_eff by running 'Cell-Problem'
+#
+# -> CLASSIFY ...
+#
+# OUTPUT: Minimizer G, angle , type, curvature
+# -----------------------------------------------------------------------
+#
+#
+# def GetMuGamma(beta,theta,gamma,mu1,rho1, InputFilePath = os.path.dirname(os.getcwd()) +"/inputs/computeMuGamma.parset",
+# OutputFilePath = os.path.dirname(os.getcwd()) + "/outputs/outputMuGamma.txt" ):
+# # ------------------------------------ get mu_gamma ------------------------------
+# # ---Scenario 1.1: extreme regimes
+# if gamma == '0':
+# print('extreme regime: gamma = 0')
+# mu_gamma = (1.0/6.0)*arithmeticMean(mu1, beta, theta) # = q2
+# print("mu_gamma:", mu_gamma)
+# elif gamma == 'infinity':
+# print('extreme regime: gamma = infinity')
+# mu_gamma = (1.0/6.0)*harmonicMean(mu1, beta, theta) # = q1
+# print("mu_gamma:", mu_gamma)
+# else:
+# # --- Scenario 1.2: compute mu_gamma with 'Compute_MuGamma' (much faster than running full Cell-Problem)
+# # print("Run computeMuGamma for Gamma = ", gamma)
+# with open(InputFilePath, 'r') as file:
+# filedata = file.read()
+# filedata = re.sub('(?m)^gamma=.*','gamma='+str(gamma),filedata)
+# # filedata = re.sub('(?m)^alpha=.*','alpha='+str(alpha),filedata)
+# filedata = re.sub('(?m)^beta=.*','beta='+str(beta),filedata)
+# filedata = re.sub('(?m)^theta=.*','theta='+str(theta),filedata)
+# filedata = re.sub('(?m)^mu1=.*','mu1='+str(mu1),filedata)
+# filedata = re.sub('(?m)^rho1=.*','rho1='+str(rho1),filedata)
+# f = open(InputFilePath,'w')
+# f.write(filedata)
+# f.close()
+# # --- Run Cell-Problem
+#
+# # Check Time
+# # t = time.time()
+# # subprocess.run(['./build-cmake/src/Cell-Problem', './inputs/cellsolver.parset'],
+# # capture_output=True, text=True)
+# # --- Run Cell-Problem_muGama -> faster
+# # subprocess.run(['./build-cmake/src/Cell-Problem_muGamma', './inputs/cellsolver.parset'],
+# # capture_output=True, text=True)
+# # --- Run Compute_muGamma (2D Problem much much faster)
+#
+# subprocess.run(['./build-cmake/src/Compute_MuGamma', './inputs/computeMuGamma.parset'],
+# capture_output=True, text=True)
+# # print('elapsed time:', time.time() - t)
+#
+# #Extract mu_gamma from Output-File TODO: GENERALIZED THIS FOR QUANTITIES OF INTEREST
+# with open(OutputFilePath, 'r') as file:
+# output = file.read()
+# tmp = re.search(r'(?m)^mu_gamma=.*',output).group() # Not necessary for Intention of Program t output Minimizer etc.....
+# s = re.findall(r"[-+]?\d*\.\d+|\d+", tmp)
+# mu_gamma = float(s[0])
+# # print("mu_gamma:", mu_gammaValue)
+# # --------------------------------------------------------------------------------------
+# return mu_gamma
+#
+
+
+
+# ----------- SETUP PATHS
+# InputFile = "/inputs/cellsolver.parset"
+# OutputFile = "/outputs/output.txt"
+InputFile = "/inputs/computeMuGamma.parset"
+OutputFile = "/outputs/outputMuGamma.txt"
+# --------- Run from src folder:
+path_parent = os.path.dirname(os.getcwd())
+os.chdir(path_parent)
+path = os.getcwd()
+print(path)
+InputFilePath = os.getcwd()+InputFile
+OutputFilePath = os.getcwd()+OutputFile
+print("InputFilepath: ", InputFilePath)
+print("OutputFilepath: ", OutputFilePath)
+print("Path: ", path)
+
+
+# -------------------------- Input Parameters --------------------
+# mu1 = 10.0 # TODO : here must be the same values as in the Parset for computeMuGamma
+mu1 = 1.0
+rho1 = 1.0
+alpha = 2.0
+beta = 2.0
+# beta = 5.0
+theta = 1.0/4.0
+#set gamma either to 1. '0' 2. 'infinity' or 3. a numerical positive value
+gamma = '0'
+# gamma = 'infinity'
+# gamma = 0.5
+# gamma = 0.25
+# gamma = 1.0
+
+# gamma = 5.0
+
+#added
+# lambda1 = 10.0
+lambda1 = 0.0
+
+#Test:
+# rho1 = -1.0
+
+
+
+print('---- Input parameters: -----')
+print('mu1: ', mu1)
+print('rho1: ', rho1)
+print('alpha: ', alpha)
+print('beta: ', beta)
+print('theta: ', theta)
+print('gamma:', gamma)
+
+print('lambda1: ', lambda1)
+print('----------------------------')
+# ----------------------------------------------------------------
+
+#
+# gamma_min = 0.5
+# gamma_max = 1.0
+#
+# # gamma_min = 1
+# # gamma_max = 1
+# Gamma_Values = np.linspace(gamma_min, gamma_max, num=3)
+# # #
+# # # Gamma_Values = np.linspace(gamma_min, gamma_max, num=13) # TODO variable Input Parameters...alpha,beta...
+# print('(Input) Gamma_Values:', Gamma_Values)
+
+print('type of gamma:', type(gamma))
+# # #
+# Gamma_Values = ['0', 'infinity']
+Gamma_Values = ['infinity']
+Gamma_Values = ['0']
+print('(Input) Gamma_Values:', Gamma_Values)
+
+for gamma in Gamma_Values:
+
+ print('Run for gamma = ', gamma)
+ print('type of gamma:', type(gamma))
+ # muGamma = GetMuGamma(beta,theta,gamma,mu1,rho1,InputFilePath)
+ # # muGamma = GetMuGamma(beta,theta,gamma,mu1,rho1)
+ # print('Test MuGamma:', muGamma)
+
+ # ------- Options --------
+ # print_Cases = True
+ # print_Output = True
+
+ #TODO
+ # generalCase = True #Read Output from Cell-Problem instead of using Lemma1.4 (special case)
+ generalCase = False
+
+ # make_3D_plot = True
+ # make_3D_PhaseDiagram = True
+ make_2D_plot = False
+ make_2D_PhaseDiagram = False
+ make_3D_plot = False
+ make_3D_PhaseDiagram = False
+ make_2D_plot = True
+ make_2D_PhaseDiagram = True
+ #
+
+ # --- Define effective quantities: q1, q2 , q3 = mu_gamma, q12 ---
+ # q1 = harmonicMean(mu1, beta, theta)
+ # q2 = arithmeticMean(mu1, beta, theta)
+ # --- Set q12
+ # q12 = 0.0 # (analytical example) # TEST / TODO read from Cell-Output
+
+
+
+
+
+ # b1 = prestrain_b1(rho1, beta, alpha, theta)
+ # b2 = prestrain_b2(rho1, beta, alpha, theta)
+ #
+ # print('---- Input parameters: -----')
+ # print('mu1: ', mu1)
+ # print('rho1: ', rho1)
+ # print('alpha: ', alpha)
+ # print('beta: ', beta)
+ # print('theta: ', theta)
+ # print("q1: ", q1)
+ # print("q2: ", q2)
+ # print("mu_gamma: ", mu_gamma)
+ # print("q12: ", q12)
+ # print("b1: ", b1)
+ # print("b2: ", b2)
+ # print('----------------------------')
+ # print("machine epsilon", sys.float_info.epsilon)
+
+ # G, angle, type, kappa = classifyMin(q1, q2, mu_gamma, q12, b1, b2, print_Cases, print_Output)
+ # Test = f(1,2 ,q1,q2,mu_gamma,q12,b1,b2)
+ # print("Test", Test)
+
+ # ---------------------- MAKE PLOT / Write to VTK------------------------------------------------------------------------------
+
+ # SamplePoints_3D = 10 # Number of sample points in each direction
+ # SamplePoints_2D = 10 # Number of sample points in each direction
+ SamplePoints_3D = 300 # Number of sample points in each direction
+ # SamplePoints_3D = 150 # Number of sample points in each direction
+ # SamplePoints_3D = 100 # Number of sample points in each direction
+ # SamplePoints_3D = 200 # Number of sample points in each direction
+ # SamplePoints_3D = 400 # Number of sample points in each direction
+ # SamplePoints_2D = 7500 # Number of sample points in each direction
+ # SamplePoints_2D = 4000 # 4000 # Number of sample points in each direction
+ # SamplePoints_2D = 400 # 4000 # Number of sample points in each direction
+ # SamplePoints_2D = 500 # 4000 # Number of sample points in each direction
+ # SamplePoints_2D = 100 # 4000 # Number of sample points in each direction
+ # SamplePoints_2D = 2000 # 4000 # Number of sample points in each direction
+ SamplePoints_2D = 1000 # 4000 # Number of sample points in each direction
+
+ if make_3D_PhaseDiagram:
+ alphas_ = np.linspace(-20, 20, SamplePoints_3D)
+ # alphas_ = np.linspace(-10, 10, SamplePoints_3D)
+
+ # betas_ = np.linspace(0.01,40.01,SamplePoints_3D) # Full Range
+ # betas_ = np.linspace(0.01,20.01,SamplePoints_3D) # FULL Range
+
+
+
+ # betas_ = np.linspace(0.01,0.99,SamplePoints_3D) # weird part
+ betas_ = np.linspace(1.01,40.01,SamplePoints_3D) #TEST !!!!! For Beta <1 weird tings happen...
+ thetas_ = np.linspace(0.01,0.99,SamplePoints_3D)
+
+
+ # TEST
+ # alphas_ = np.linspace(-2, 2, SamplePoints_3D)
+ # betas_ = np.linspace(1.01,10.01,SamplePoints_3D)
+ # print('betas:', betas_)
+
+ # TEST :
+ # alphas_ = np.linspace(-40, 40, SamplePoints_3D)
+ # betas_ = np.linspace(0.01,80.01,SamplePoints_3D) # Full Range
+
+ # print('type of alphas', type(alphas_))
+ # print('Test:', type(np.array([mu_gamma])) )
+ alphas, betas, thetas = np.meshgrid(alphas_, betas_, thetas_, indexing='ij')
+ classifyMin_anaVec = np.vectorize(classifyMin_ana)
+
+ # Get MuGamma values ...
+ GetMuGammaVec = np.vectorize(GetMuGamma)
+ muGammas = GetMuGammaVec(betas, thetas, gamma, mu1, rho1)
+ # Classify Minimizers....
+ G, angles, Types, curvature = classifyMin_anaVec(alphas, betas, thetas, muGammas, mu1, rho1) # Sets q12 to zero!!!
+
+ # G, angles, Types, curvature = classifyMin_anaVec(alphas, betas, thetas, muGammas, mu1, rho1, True, True)
+ # print('size of G:', G.shape)
+ # print('G:', G)
+
+ # Option to print angles
+ # print('angles:', angles)
+
+
+ # Out = classifyMin_anaVec(alphas,betas,thetas)
+ # T = Out[2]
+ # --- Write to VTK
+
+ GammaString = str(gamma)
+ VTKOutputName = "outputs/PhaseDiagram3D" + "Gamma" + GammaString
+ gridToVTK(VTKOutputName , alphas, betas, thetas, pointData = {'Type': Types, 'angles': angles, 'curvature': curvature} )
+ print('Written to VTK-File:', VTKOutputName )
+
+ if make_2D_PhaseDiagram:
+ # alphas_ = np.linspace(-20, 20, SamplePoints_2D)
+ # alphas_ = np.linspace(0, 1, SamplePoints_2D)
+ thetas_ = np.linspace(0.01,0.99,SamplePoints_2D)
+ alphas_ = np.linspace(-5, 5, SamplePoints_2D)
+ # alphas_ = np.linspace(-5, 15, SamplePoints_2D)
+ # thetas_ = np.linspace(0.05,0.25,SamplePoints_2D)
+
+
+ # good range:
+ # alphas_ = np.linspace(9, 10, SamplePoints_2D)
+ # thetas_ = np.linspace(0.075,0.14,SamplePoints_2D)
+
+ # range used:
+ # alphas_ = np.linspace(8, 10, SamplePoints_2D)
+ # thetas_ = np.linspace(0.05,0.16,SamplePoints_2D)
+
+ # alphas_ = np.linspace(8, 12, SamplePoints_2D)
+ # thetas_ = np.linspace(0.05,0.2,SamplePoints_2D)
+ # betas_ = np.linspace(0.01,40.01,1)
+ #fix to one value:
+ betas_ = 2.0;
+ # betas_ = 10.0;
+ # betas_ = 5.0;
+ # betas_ = 0.5;
+
+
+ #intermediate Values
+ # alphas_ = np.linspace(-2, 1, SamplePoints_2D)
+ # thetas_ = np.linspace(0.4,0.6,SamplePoints_2D)
+ # betas_ = 10.0;
+
+ # TEST
+ # alphas_ = np.linspace(-8, 8, SamplePoints_2D)
+ # thetas_ = np.linspace(0.01,0.99,SamplePoints_2D)
+ # betas_ = 1.0; #TEST Problem: disvison by zero if alpha = 9, theta = 0.1 !
+ # betas_ = 0.9;
+ # betas_ = 0.5; #TEST!!!
+ # alphas, betas, thetas = np.meshgrid(alphas_, betas_, thetas_, indexing='ij')
+ betas = betas_
+ alphas, thetas = np.meshgrid(alphas_, thetas_, indexing='ij')
+
+ if generalCase:
+ classifyMin_matVec = np.vectorize(classifyMin_mat)
+ GetCellOutputVec = np.vectorize(GetCellOutput, otypes=[np.ndarray, np.ndarray])
+ Q, B = GetCellOutputVec(alphas,betas,thetas,gamma,mu1,rho1,lambda1, InputFilePath ,OutputFilePath )
+
+
+ # print('type of Q:', type(Q))
+ # print('Q:', Q)
+ G, angles, Types, curvature = classifyMin_matVec(Q,B)
+
+ else:
+ classifyMin_anaVec = np.vectorize(classifyMin_ana)
+ GetMuGammaVec = np.vectorize(GetMuGamma)
+ # muGammas = GetMuGammaVec(betas,thetas,gamma,mu1,rho1,InputFilePath ,OutputFilePath )
+ # G, angles, Types, curvature = classifyMin_anaVec(alphas,betas,thetas, muGammas, mu1, rho1) # Sets q12 to zero!!!
+ muGammas = GetMuGammaVec(betas,thetas,gamma,mu1,rho1,InputFilePath ,OutputFilePath )
+ G, angles, Types, curvature = classifyMin_anaVec(alphas,betas,thetas, muGammas, mu1, rho1) # Sets q12 to zero!!!
+
+ # print('size of G:', G.shape)
+ # print('G:', G)
+ # print('Types:', Types)
+ # Out = classifyMin_anaVec(alphas,betas,thetas)
+ # T = Out[2]
+ # --- Write to VTK
+ # VTKOutputName = + path + "./PhaseDiagram2DNEW"
+
+ print('angles:',angles)
+ # GammaString = str(gamma)
+ # VTKOutputName = "outputs/PhaseDiagram2D" + "Gamma_" + GammaString
+ # gridToVTK(VTKOutputName , alphas, betas, thetas, pointData = {'Type': Types, 'angles': angles, 'curvature': curvature} )
+ # print('Written to VTK-File:', VTKOutputName )
+
+
+ # --- Make 3D Scatter plot
+ if(make_3D_plot or make_2D_plot):
+ # fig = plt.figure()
+ # ax = fig.add_subplot(111, projection='3d')
+ # colors = cm.plasma(Types)
+ colors = cm.coolwarm(angles)
+
+
+ width = 6.28
+ height = width / 1.618
+ # height = width / 2.5
+ fig, ax = plt.subplots()
+ # ax = plt.axes((0.15,0.21 ,0.8,0.75))
+
+
+
+
+ # if make_2D_plot: pnt3d=ax.scatter(alphas,thetas,c=Types.flat)
+ # if make_3D_plot: pnt3d=ax.scatter(alphas,betas,thetas,c=Types.flat)
+ #
+ # if make_2D_plot: pnt3d=ax.scatter(alphas,thetas,c=angles.flat)
+ # if make_3D_plot: pnt3d=ax.scatter(alphas,betas,thetas,c=angles.flat)
+
+
+ # pnt=ax.scatter(alphas,thetas,c=angles,cmap='coolwarm')
+ # # ax.colorbar()
+ # CS = ax.contourf(alphas, thetas, angles,6, cmap=plt.cm.coolwarm, linestyle=dashed)
+ # # CS = ax.contour(alphas, thetas, angles,6, colors='k')
+ # ax.clabel(CS, inline=True, fontsize=7.5)
+ # # ax.set_title('Simplest default with labels')
+
+ if gamma =='0':
+ CS = ax.contourf(alphas, thetas, angles, 10, cmap=plt.cm.coolwarm)
+ # CS = ax.contourf(alphas, thetas, angles, 10, cmap='RdBu')
+ CS2 = ax.contour(CS, levels=CS.levels[::2], colors='black',inline=True, linewidths=(0.5,))
+ # ax.clabel(CS2, inline=True, fontsize=9, colors='black')
+ # ax.clabel(CS2, inline=True, inline_spacing=3, rightside_up=True, colors='k', fontsize=8)
+ manual_locations = [
+ (-0.5, 0.3), (-0.7, 0.4), (-0.8, 0.5), (-0.9, 0.6), (-1,0.7)]
+ # ax.clabel(CS2, inline=True, fontsize=6, colors='black', manual=manual_locations)
+ # ax.clabel(CS2, inline=True, fontsize=6, colors='black')
+ # ax.clabel(CS2, CS2.levels, inline=True, fontsize=10)
+ # ax.clabel(CS, fontsize=5, colors='black')
+ # cbar = fig.colorbar(CS,label=r'angle $\alpha$', ticks=[0, np.pi/8, np.pi/4, 3*np.pi/8 , np.pi/2 ])
+ cbar = fig.colorbar(CS, ticks=[0, np.pi/2 ])
+ cbar.ax.set_yticklabels(['$0$', r'$\pi/2$'])
+ cbar.ax.set_title(r'angle $\alpha$')
+
+ if gamma == 'infinity':
+ CS = ax.contourf(alphas, thetas, angles, 10, cmap=plt.cm.coolwarm)
+ # CS = ax.contourf(alphas, thetas, angles, 10, cmap='RdBu')
+ CS2 = ax.contour(CS, levels=CS.levels[::2], colors='black',inline=True, linewidths=(0.5,))
+ # ax.clabel(CS2, inline=True, fontsize=9, colors='black')
+ # ax.clabel(CS2, inline=True, inline_spacing=3, rightside_up=True, colors='k', fontsize=8)
+ manual_locations = [
+ (-0.5, 0.3), (-0.7, 0.4), (-0.8, 0.5), (-0.9, 0.6), (-1,0.7)]
+ ax.clabel(CS2, inline=True, fontsize=6, colors='black', manual=manual_locations)
+
+ # ax.clabel(CS2, CS2.levels, inline=True, fontsize=10)
+ # ax.clabel(CS, fontsize=5, colors='black')
+ # cbar = fig.colorbar(CS,label=r'angle $\alpha$', ticks=[0, np.pi/8, np.pi/4, 3*np.pi/8 , np.pi/2 ])
+ cbar = fig.colorbar(CS, ticks=[0, np.pi/8, np.pi/4, 3*np.pi/8 , np.pi/2 ])
+ cbar.ax.set_yticklabels(['$0$',r'$\pi/8$', r'$\pi/4$' ,r'$3\pi/8$' , r'$\pi/2$'])
+ cbar.ax.set_title(r'angle $\alpha$')
+ # cbar=plt.colorbar(pnt3d)
+ # cbar.set_label("Values (units)")
+ # plt.axvline(x = 8, color = 'b', linestyle = ':', label='$q_1$')
+ # plt.axhline(y = 0.083333333, color = 'b', linestyle = ':', label='$q_1$')
+
+ ax.set_xlabel(r'$\theta_\rho$',fontsize=10)
+ # ax.xaxis.set_minor_locator(MultipleLocator(0.5))
+ ax.yaxis.set_major_locator(MultipleLocator(0.1))
+ ax.xaxis.set_major_locator(MultipleLocator(1))
+ # ax.set_ylabel('beta')
+ ax.set_ylabel(r'$\theta$ ',fontsize=10, rotation=0)
+ # if make_3D_plot: ax.set_zlabel('theta')
+ plt.subplots_adjust(bottom=0.2)
+ plt.grid( linestyle = '--', linewidth = 0.25)
+
+ fig.set_size_inches(width, height)
+ outputName = 'Plot-Contour_Gamma' +str(gamma) + '.pdf'
+ fig.savefig(outputName)
+ # fig.savefig('Plot-Contour.pdf')
+ plt.show()
+ # plt.savefig('common_labels.png', dpi=300)
+ # print('T:', T)
+ # print('Type 1 occured here:', np.where(T == 1))
+ # print('Type 2 occured here:', np.where(T == 2))
+
+
+ # print(alphas_)
+ # print(betas_)
+
+
+
+
+
+# ALTERNATIVE
+# colors = ("red", "green", "blue")
+# groups = ("Type 1", "Type2", "Type3")
+#
+# # Create plot
+# fig = plt.figure()
+# ax = fig.add_subplot(1, 1, 1)
+#
+# for data, color, group in zip(Types, colors, groups):
+# # x, y = data
+# ax.scatter(alphas, thetas, alpha=0.8, c=color, edgecolors='none', label=group)
+#
+# plt.title('Matplot scatter plot')
+# plt.legend(loc=2)
+# plt.show()
diff --git a/src/PhaseDiagram_ContourSubPlots.py b/src/PhaseDiagram_ContourSubPlots.py
new file mode 100644
index 0000000000000000000000000000000000000000..59a1a421f3788f05ee57bbd32e168c58866e4d3b
--- /dev/null
+++ b/src/PhaseDiagram_ContourSubPlots.py
@@ -0,0 +1,529 @@
+import numpy as np
+import matplotlib.pyplot as plt
+import sympy as sym
+import math
+import os
+import subprocess
+import fileinput
+import re
+import matlab.engine
+import sys
+from ClassifyMin import *
+from HelperFunctions import *
+# from CellScript import *
+from mpl_toolkits.mplot3d import Axes3D
+import matplotlib.cm as cm
+from vtk.util import numpy_support
+from pyevtk.hl import gridToVTK
+from matplotlib.ticker import MultipleLocator,FormatStrFormatter,MaxNLocator
+
+from mpl_toolkits.axes_grid1.inset_locator import inset_axes
+
+import time
+# print(sys.executable)
+
+# --------------------------------------------------------------------
+# START :
+# INPUT (Parameters): alpha, beta, theta, gamma, mu1, rho1
+#
+# -Option 1 : (Case lambda = 0 => q12 = 0)
+# compute q1,q2,b1,b2 from Formula
+# Option 1.1 :
+# set mu_gamma = 'q1' or 'q2' (extreme regimes: gamma \in {0,\infty})
+# Option 1.2 :
+# compute mu_gamma with 'Compute_MuGamma' (2D problem much faster then Cell-Problem)
+# -Option 2 :
+# compute Q_hom & B_eff by running 'Cell-Problem'
+#
+# -> CLASSIFY ...
+#
+# OUTPUT: Minimizer G, angle , type, curvature
+# -----------------------------------------------------------------------
+#
+#
+# def GetMuGamma(beta,theta,gamma,mu1,rho1, InputFilePath = os.path.dirname(os.getcwd()) +"/inputs/computeMuGamma.parset",
+# OutputFilePath = os.path.dirname(os.getcwd()) + "/outputs/outputMuGamma.txt" ):
+# # ------------------------------------ get mu_gamma ------------------------------
+# # ---Scenario 1.1: extreme regimes
+# if gamma == '0':
+# print('extreme regime: gamma = 0')
+# mu_gamma = (1.0/6.0)*arithmeticMean(mu1, beta, theta) # = q2
+# print("mu_gamma:", mu_gamma)
+# elif gamma == 'infinity':
+# print('extreme regime: gamma = infinity')
+# mu_gamma = (1.0/6.0)*harmonicMean(mu1, beta, theta) # = q1
+# print("mu_gamma:", mu_gamma)
+# else:
+# # --- Scenario 1.2: compute mu_gamma with 'Compute_MuGamma' (much faster than running full Cell-Problem)
+# # print("Run computeMuGamma for Gamma = ", gamma)
+# with open(InputFilePath, 'r') as file:
+# filedata = file.read()
+# filedata = re.sub('(?m)^gamma=.*','gamma='+str(gamma),filedata)
+# # filedata = re.sub('(?m)^alpha=.*','alpha='+str(alpha),filedata)
+# filedata = re.sub('(?m)^beta=.*','beta='+str(beta),filedata)
+# filedata = re.sub('(?m)^theta=.*','theta='+str(theta),filedata)
+# filedata = re.sub('(?m)^mu1=.*','mu1='+str(mu1),filedata)
+# filedata = re.sub('(?m)^rho1=.*','rho1='+str(rho1),filedata)
+# f = open(InputFilePath,'w')
+# f.write(filedata)
+# f.close()
+# # --- Run Cell-Problem
+#
+# # Check Time
+# # t = time.time()
+# # subprocess.run(['./build-cmake/src/Cell-Problem', './inputs/cellsolver.parset'],
+# # capture_output=True, text=True)
+# # --- Run Cell-Problem_muGama -> faster
+# # subprocess.run(['./build-cmake/src/Cell-Problem_muGamma', './inputs/cellsolver.parset'],
+# # capture_output=True, text=True)
+# # --- Run Compute_muGamma (2D Problem much much faster)
+#
+# subprocess.run(['./build-cmake/src/Compute_MuGamma', './inputs/computeMuGamma.parset'],
+# capture_output=True, text=True)
+# # print('elapsed time:', time.time() - t)
+#
+# #Extract mu_gamma from Output-File TODO: GENERALIZED THIS FOR QUANTITIES OF INTEREST
+# with open(OutputFilePath, 'r') as file:
+# output = file.read()
+# tmp = re.search(r'(?m)^mu_gamma=.*',output).group() # Not necessary for Intention of Program t output Minimizer etc.....
+# s = re.findall(r"[-+]?\d*\.\d+|\d+", tmp)
+# mu_gamma = float(s[0])
+# # print("mu_gamma:", mu_gammaValue)
+# # --------------------------------------------------------------------------------------
+# return mu_gamma
+#
+
+
+
+# ----------- SETUP PATHS
+# InputFile = "/inputs/cellsolver.parset"
+# OutputFile = "/outputs/output.txt"
+InputFile = "/inputs/computeMuGamma.parset"
+OutputFile = "/outputs/outputMuGamma.txt"
+# --------- Run from src folder:
+path_parent = os.path.dirname(os.getcwd())
+os.chdir(path_parent)
+path = os.getcwd()
+print(path)
+InputFilePath = os.getcwd()+InputFile
+OutputFilePath = os.getcwd()+OutputFile
+print("InputFilepath: ", InputFilePath)
+print("OutputFilepath: ", OutputFilePath)
+print("Path: ", path)
+
+
+# -------------------------- Input Parameters --------------------
+# mu1 = 10.0 # TODO : here must be the same values as in the Parset for computeMuGamma
+mu1 = 1.0
+rho1 = 1.0
+alpha = 2.0
+beta = 2.0
+# beta = 5.0
+theta = 1.0/4.0
+#set gamma either to 1. '0' 2. 'infinity' or 3. a numerical positive value
+gamma = '0'
+# gamma = 'infinity'
+# gamma = 0.5
+# gamma = 0.25
+# gamma = 1.0
+
+# gamma = 5.0
+
+#added
+# lambda1 = 10.0
+lambda1 = 0.0
+
+#Test:
+# rho1 = -1.0
+
+
+
+print('---- Input parameters: -----')
+print('mu1: ', mu1)
+print('rho1: ', rho1)
+print('alpha: ', alpha)
+print('beta: ', beta)
+print('theta: ', theta)
+print('gamma:', gamma)
+
+print('lambda1: ', lambda1)
+print('----------------------------')
+# ----------------------------------------------------------------
+
+#
+# gamma_min = 0.5
+# gamma_max = 1.0
+#
+# # gamma_min = 1
+# # gamma_max = 1
+# Gamma_Values = np.linspace(gamma_min, gamma_max, num=3)
+# # #
+# # # Gamma_Values = np.linspace(gamma_min, gamma_max, num=13) # TODO variable Input Parameters...alpha,beta...
+# print('(Input) Gamma_Values:', Gamma_Values)
+
+print('type of gamma:', type(gamma))
+# # #
+Gamma_Values = ['0', 'infinity']
+# Gamma_Values = ['infinity']
+# Gamma_Values = ['0']
+print('(Input) Gamma_Values:', Gamma_Values)
+
+for gamma in Gamma_Values:
+
+ print('Run for gamma = ', gamma)
+ print('type of gamma:', type(gamma))
+ # muGamma = GetMuGamma(beta,theta,gamma,mu1,rho1,InputFilePath)
+ # # muGamma = GetMuGamma(beta,theta,gamma,mu1,rho1)
+ # print('Test MuGamma:', muGamma)
+
+ # ------- Options --------
+ # print_Cases = True
+ # print_Output = True
+
+ #TODO
+ # generalCase = True #Read Output from Cell-Problem instead of using Lemma1.4 (special case)
+ generalCase = False
+
+ # make_3D_plot = True
+ # make_3D_PhaseDiagram = True
+ make_2D_plot = False
+ make_2D_PhaseDiagram = False
+ make_3D_plot = False
+ make_3D_PhaseDiagram = False
+ make_2D_plot = True
+ make_2D_PhaseDiagram = True
+ #
+
+ # --- Define effective quantities: q1, q2 , q3 = mu_gamma, q12 ---
+ # q1 = harmonicMean(mu1, beta, theta)
+ # q2 = arithmeticMean(mu1, beta, theta)
+ # --- Set q12
+ # q12 = 0.0 # (analytical example) # TEST / TODO read from Cell-Output
+
+
+
+
+
+ # b1 = prestrain_b1(rho1, beta, alpha, theta)
+ # b2 = prestrain_b2(rho1, beta, alpha, theta)
+ #
+ # print('---- Input parameters: -----')
+ # print('mu1: ', mu1)
+ # print('rho1: ', rho1)
+ # print('alpha: ', alpha)
+ # print('beta: ', beta)
+ # print('theta: ', theta)
+ # print("q1: ", q1)
+ # print("q2: ", q2)
+ # print("mu_gamma: ", mu_gamma)
+ # print("q12: ", q12)
+ # print("b1: ", b1)
+ # print("b2: ", b2)
+ # print('----------------------------')
+ # print("machine epsilon", sys.float_info.epsilon)
+
+ # G, angle, type, kappa = classifyMin(q1, q2, mu_gamma, q12, b1, b2, print_Cases, print_Output)
+ # Test = f(1,2 ,q1,q2,mu_gamma,q12,b1,b2)
+ # print("Test", Test)
+
+ # ---------------------- MAKE PLOT / Write to VTK------------------------------------------------------------------------------
+
+ # SamplePoints_3D = 10 # Number of sample points in each direction
+ # SamplePoints_2D = 10 # Number of sample points in each direction
+ SamplePoints_3D = 300 # Number of sample points in each direction
+ # SamplePoints_3D = 150 # Number of sample points in each direction
+ # SamplePoints_3D = 100 # Number of sample points in each direction
+ # SamplePoints_3D = 200 # Number of sample points in each direction
+ # SamplePoints_3D = 400 # Number of sample points in each direction
+ # SamplePoints_2D = 7500 # Number of sample points in each direction
+ # SamplePoints_2D = 4000 # 4000 # Number of sample points in each direction
+ # SamplePoints_2D = 400 # 4000 # Number of sample points in each direction
+ # SamplePoints_2D = 500 # 4000 # Number of sample points in each direction
+ # SamplePoints_2D = 100 # 4000 # Number of sample points in each direction
+ SamplePoints_2D = 2000 # 4000 # Number of sample points in each direction
+ # SamplePoints_2D = 1000 # 4000 # Number of sample points in each direction
+
+ if make_3D_PhaseDiagram:
+ alphas_ = np.linspace(-20, 20, SamplePoints_3D)
+ # alphas_ = np.linspace(-10, 10, SamplePoints_3D)
+
+ # betas_ = np.linspace(0.01,40.01,SamplePoints_3D) # Full Range
+ # betas_ = np.linspace(0.01,20.01,SamplePoints_3D) # FULL Range
+
+
+
+ # betas_ = np.linspace(0.01,0.99,SamplePoints_3D) # weird part
+ betas_ = np.linspace(1.01,40.01,SamplePoints_3D) #TEST !!!!! For Beta <1 weird tings happen...
+ thetas_ = np.linspace(0.01,0.99,SamplePoints_3D)
+
+
+ # TEST
+ # alphas_ = np.linspace(-2, 2, SamplePoints_3D)
+ # betas_ = np.linspace(1.01,10.01,SamplePoints_3D)
+ # print('betas:', betas_)
+
+ # TEST :
+ # alphas_ = np.linspace(-40, 40, SamplePoints_3D)
+ # betas_ = np.linspace(0.01,80.01,SamplePoints_3D) # Full Range
+
+ # print('type of alphas', type(alphas_))
+ # print('Test:', type(np.array([mu_gamma])) )
+ alphas, betas, thetas = np.meshgrid(alphas_, betas_, thetas_, indexing='ij')
+ classifyMin_anaVec = np.vectorize(classifyMin_ana)
+
+ # Get MuGamma values ...
+ GetMuGammaVec = np.vectorize(GetMuGamma)
+ muGammas = GetMuGammaVec(betas, thetas, gamma, mu1, rho1)
+ # Classify Minimizers....
+ G, angles, Types, curvature = classifyMin_anaVec(alphas, betas, thetas, muGammas, mu1, rho1) # Sets q12 to zero!!!
+
+ # G, angles, Types, curvature = classifyMin_anaVec(alphas, betas, thetas, muGammas, mu1, rho1, True, True)
+ # print('size of G:', G.shape)
+ # print('G:', G)
+
+ # Option to print angles
+ # print('angles:', angles)
+
+
+ # Out = classifyMin_anaVec(alphas,betas,thetas)
+ # T = Out[2]
+ # --- Write to VTK
+
+ GammaString = str(gamma)
+ VTKOutputName = "outputs/PhaseDiagram3D" + "Gamma" + GammaString
+ gridToVTK(VTKOutputName , alphas, betas, thetas, pointData = {'Type': Types, 'angles': angles, 'curvature': curvature} )
+ print('Written to VTK-File:', VTKOutputName )
+
+ if make_2D_PhaseDiagram:
+ # alphas_ = np.linspace(-20, 20, SamplePoints_2D)
+ # alphas_ = np.linspace(0, 1, SamplePoints_2D)
+ thetas_ = np.linspace(0.01,0.99,SamplePoints_2D)
+ alphas_ = np.linspace(-5, 5, SamplePoints_2D)
+ # alphas_ = np.linspace(-5, 15, SamplePoints_2D)
+ # thetas_ = np.linspace(0.05,0.25,SamplePoints_2D)
+
+
+ # good range:
+ # alphas_ = np.linspace(9, 10, SamplePoints_2D)
+ # thetas_ = np.linspace(0.075,0.14,SamplePoints_2D)
+
+ # range used:
+ # alphas_ = np.linspace(8, 10, SamplePoints_2D)
+ # thetas_ = np.linspace(0.05,0.16,SamplePoints_2D)
+
+ # alphas_ = np.linspace(8, 12, SamplePoints_2D)
+ # thetas_ = np.linspace(0.05,0.2,SamplePoints_2D)
+ # betas_ = np.linspace(0.01,40.01,1)
+ #fix to one value:
+ betas_ = 2.0;
+ # betas_ = 10.0;
+ # betas_ = 5.0;
+ # betas_ = 0.5;
+
+
+ #intermediate Values
+ alphas_ = np.linspace(-2, 1, SamplePoints_2D)
+ # thetas_ = np.linspace(0.4,0.6,SamplePoints_2D)
+ # betas_ = 10.0;
+
+ # TEST
+ # alphas_ = np.linspace(-8, 8, SamplePoints_2D)
+ # thetas_ = np.linspace(0.01,0.99,SamplePoints_2D)
+ # betas_ = 1.0; #TEST Problem: disvison by zero if alpha = 9, theta = 0.1 !
+ # betas_ = 0.9;
+ # betas_ = 0.5; #TEST!!!
+ # alphas, betas, thetas = np.meshgrid(alphas_, betas_, thetas_, indexing='ij')
+ betas = betas_
+ alphas, thetas = np.meshgrid(alphas_, thetas_, indexing='ij')
+
+ if generalCase:
+ classifyMin_matVec = np.vectorize(classifyMin_mat)
+ GetCellOutputVec = np.vectorize(GetCellOutput, otypes=[np.ndarray, np.ndarray])
+ Q, B = GetCellOutputVec(alphas,betas,thetas,gamma,mu1,rho1,lambda1, InputFilePath ,OutputFilePath )
+
+
+ # print('type of Q:', type(Q))
+ # print('Q:', Q)
+ G, angles, Types, curvature = classifyMin_matVec(Q,B)
+
+ else:
+ classifyMin_anaVec = np.vectorize(classifyMin_ana)
+ GetMuGammaVec = np.vectorize(GetMuGamma)
+ # muGammas = GetMuGammaVec(betas,thetas,gamma,mu1,rho1,InputFilePath ,OutputFilePath )
+ # G, angles, Types, curvature = classifyMin_anaVec(alphas,betas,thetas, muGammas, mu1, rho1) # Sets q12 to zero!!!
+ muGammas = GetMuGammaVec(betas,thetas,gamma,mu1,rho1,InputFilePath ,OutputFilePath )
+
+ if gamma == '0':
+ G, angles_0, Types, curvature = classifyMin_anaVec(alphas,betas,thetas, muGammas, mu1, rho1) # Sets q12 to zero!!!
+ if gamma == 'infinity':
+ G, angles_inf, Types, curvature = classifyMin_anaVec(alphas,betas,thetas, muGammas, mu1, rho1) # Sets q12 to zero!!!
+
+ # print('size of G:', G.shape)
+ # print('G:', G)
+ # print('Types:', Types)
+ # Out = classifyMin_anaVec(alphas,betas,thetas)
+ # T = Out[2]
+ # --- Write to VTK
+ # VTKOutputName = + path + "./PhaseDiagram2DNEW"
+
+ # print('angles:',angles)
+ # GammaString = str(gamma)
+ # VTKOutputName = "outputs/PhaseDiagram2D" + "Gamma_" + GammaString
+ # gridToVTK(VTKOutputName , alphas, betas, thetas, pointData = {'Type': Types, 'angles': angles, 'curvature': curvature} )
+ # print('Written to VTK-File:', VTKOutputName )
+
+
+# --- Make 3D Scatter plot
+if(make_3D_plot or make_2D_plot):
+ # fig = plt.figure()
+ # ax = fig.add_subplot(111, projection='3d')
+ # colors = cm.plasma(Types)
+ colors = cm.coolwarm(angles_inf)
+
+
+ width = 6.28
+ height = width / 1.618
+ # height = width / 2.5
+ # fig, ax = plt.subplots()
+ fig,ax = plt.subplots(nrows=1,ncols=2,figsize=(width,height), sharey=True)
+ # ax = plt.axes((0.15,0.21 ,0.8,0.75))
+
+
+
+
+ # if make_2D_plot: pnt3d=ax.scatter(alphas,thetas,c=Types.flat)
+ # if make_3D_plot: pnt3d=ax.scatter(alphas,betas,thetas,c=Types.flat)
+ #
+ # if make_2D_plot: pnt3d=ax.scatter(alphas,thetas,c=angles.flat)
+ # if make_3D_plot: pnt3d=ax.scatter(alphas,betas,thetas,c=angles.flat)
+
+
+ # pnt=ax.scatter(alphas,thetas,c=angles,cmap='coolwarm')
+ # # ax.colorbar()
+ # CS = ax.contourf(alphas, thetas, angles,6, cmap=plt.cm.coolwarm, linestyle=dashed)
+ # # CS = ax.contour(alphas, thetas, angles,6, colors='k')
+ # ax.clabel(CS, inline=True, fontsize=7.5)
+ # # ax.set_title('Simplest default with labels')
+
+
+ CS_0 = ax[0].contourf(alphas, thetas, angles_0, 10, cmap=plt.cm.coolwarm)
+ # CS = ax.contourf(alphas, thetas, angles, 10, cmap='RdBu')
+ CS_02 = ax[0].contour(CS_0, levels=CS_0.levels[::2], colors='black',inline=True, linewidths=(0.5,))
+ # ax.clabel(CS2, inline=True, fontsize=9, colors='black')
+ # ax.clabel(CS2, inline=True, inline_spacing=3, rightside_up=True, colors='k', fontsize=8)
+ # manual_locations = [
+ # (-0.5, 0.3), (-0.7, 0.4), (-0.8, 0.5), (-0.9, 0.6), (-1,0.7)]
+ manual_locations = [
+ (-0.4, 0.2),(-0.6, 0.3), (-0.7, 0.4), (-0.8, 0.5), (-0.9, 0.6), (-1,0.7)]
+ # ax.clabel(CS2, inline=True, fontsize=6, colors='black', manual=manual_locations)
+ # ax.clabel(CS2, inline=True, fontsize=6, colors='black')
+ # ax.clabel(CS2, CS2.levels, inline=True, fontsize=10)
+ # ax.clabel(CS, fontsize=5, colors='black')
+ # cbar = fig.colorbar(CS,label=r'angle $\alpha$', ticks=[0, np.pi/8, np.pi/4, 3*np.pi/8 , np.pi/2 ])
+ # cbar = fig.colorbar(CS_0, ticks=[0, np.pi/2 ])
+ # cbar.ax.set_yticklabels(['$0$', r'$\pi/2$'])
+ # cbar.ax.set_title(r'angle $\alpha$')
+
+
+ CS_1 = ax[1].contourf(alphas, thetas, angles_inf, 10, cmap=plt.cm.coolwarm)
+ # CS = ax.contourf(alphas, thetas, angles, 10, cmap='RdBu')
+ CS_12 = ax[1].contour(CS_1, levels=CS_1.levels[::2], colors='black',inline=True, linewidths=(0.5,))
+ # ax.clabel(CS2, inline=True, fontsize=9, colors='black')
+ # ax.clabel(CS2, inline=True, inline_spacing=3, rightside_up=True, colors='k', fontsize=8)
+ manual_locations = [
+ (-0.5, 0.3), (-0.7, 0.4), (-0.8, 0.5), (-0.9, 0.6), (-1,0.7)]
+ ax[1].clabel(CS_12, inline=True, fontsize=10, colors='black', manual=manual_locations)
+ # ax[1].clabel(CS_12, inline=True, fontsize=8, colors='black')
+
+
+ axins1 = inset_axes(ax[1],
+ width="5%", # width = 5% of parent_bbox width
+ height="100%", # height : 50%
+ loc='lower left',
+ bbox_to_anchor=(1.05, 0., 1, 1),
+ bbox_transform=ax[1].transAxes,
+ borderpad=0,
+ )
+
+ # ax.clabel(CS2, CS2.levels, inline=True, fontsize=10)
+ # ax.clabel(CS, fontsize=5, colors='black')
+ # cbar = fig.colorbar(CS,label=r'angle $\alpha$', ticks=[0, np.pi/8, np.pi/4, 3*np.pi/8 , np.pi/2 ])
+ # cbar = fig.colorbar(CS_1, ticks=[0, np.pi/8, np.pi/4, 3*np.pi/8 , np.pi/2 ])
+
+ # cbar_ax = fig.add_axes([0.85, 0.15, 0.05, 0.7])
+ cbar = fig.colorbar(CS_1, cax=axins1, ticks=[0, np.pi/8, np.pi/4, 3*np.pi/8 , np.pi/2 ])
+ # cbar = fig.colorbar(CS_1, cax=cbar_ax, shrink=0.2, location='right', ticks=[0, np.pi/8, np.pi/4, 3*np.pi/8 , np.pi/2 ])
+ # cbar = fig.colorbar(CS_1, ax=ax[:], shrink=0.8, location='right', ticks=[0, np.pi/8, np.pi/4, 3*np.pi/8 , np.pi/2 ])
+
+ cbar.ax.set_yticklabels(['$0$',r'$\pi/8$', r'$\pi/4$' ,r'$3\pi/8$' , r'$\pi/2$'])
+ cbar.ax.set_title(r'angle $\alpha$')
+ # cbar=plt.colorbar(pnt3d)
+ # cbar.set_label("Values (units)")
+ # plt.axvline(x = 8, color = 'b', linestyle = ':', label='$q_1$')
+ # plt.axhline(y = 0.083333333, color = 'b', linestyle = ':', label='$q_1$')
+
+ ax[0].set_xlabel(r'$\theta_\rho$',fontsize=10)
+ # ax[0].yaxis.set_major_locator(MultipleLocator(0.1))
+ # ax[0].xaxis.set_major_locator(MultipleLocator(1))
+ ax[0].yaxis.set_major_locator(MultipleLocator(0.1))
+ ax[0].xaxis.set_major_locator(MultipleLocator(0.5))
+ ax[0].set_ylabel(r'$\theta$ ',fontsize=10, rotation=0)
+ ax[0].tick_params(axis='x', labelsize=7 )
+ ax[0].tick_params(axis='y', labelsize=7 )
+
+ ax[1].set_xlabel(r'$\theta_\rho$',fontsize=10)
+ # ax.xaxis.set_minor_locator(MultipleLocator(0.5))
+ # ax[1].yaxis.set_major_locator(MultipleLocator(0.1))
+ # ax[1].xaxis.set_major_locator(MultipleLocator(1))
+ ax[1].yaxis.set_major_locator(MultipleLocator(0.1))
+ ax[1].xaxis.set_major_locator(MultipleLocator(0.5))
+ ax[1].tick_params(axis='x', labelsize=7 )
+ ax[1].tick_params(axis='y', labelsize=7 )
+ # ax.set_ylabel('beta')
+ # ax[1].set_ylabel(r'$\theta$ ',fontsize=10, rotation=0)
+ # if make_3D_plot: ax.set_zlabel('theta')
+ # plt.subplots_adjust(bottom=0.2)
+ # plt.subplots_adjust(wspace=0.22, hspace=0.1)
+ plt.subplots_adjust(hspace=0.15, wspace=0.1)
+
+ # fig.subplots_adjust(right=0.75)
+
+
+ ax[0].grid( linestyle = '--', linewidth = 0.25)
+ ax[1].grid( linestyle = '--', linewidth = 0.25)
+
+
+
+ fig.set_size_inches(width, height)
+ outputName = 'Plot-Contour_Gamma' +str(gamma) + '.pdf'
+ fig.savefig(outputName)
+ # fig.savefig('Plot-Contour.pdf')
+ plt.show()
+ # plt.savefig('common_labels.png', dpi=300)
+ # print('T:', T)
+ # print('Type 1 occured here:', np.where(T == 1))
+ # print('Type 2 occured here:', np.where(T == 2))
+
+
+ # print(alphas_)
+ # print(betas_)
+
+
+
+
+
+# ALTERNATIVE
+# colors = ("red", "green", "blue")
+# groups = ("Type 1", "Type2", "Type3")
+#
+# # Create plot
+# fig = plt.figure()
+# ax = fig.add_subplot(1, 1, 1)
+#
+# for data, color, group in zip(Types, colors, groups):
+# # x, y = data
+# ax.scatter(alphas, thetas, alpha=0.8, c=color, edgecolors='none', label=group)
+#
+# plt.title('Matplot scatter plot')
+# plt.legend(loc=2)
+# plt.show()