diff --git a/src/ClassifyMin.py b/src/ClassifyMin.py
index 1f78c5ef9c9151c4cfee49b3a455d8e73882f7ae..e2d08c74f37047b61c26f367482b9339ddda038d 100644
--- a/src/ClassifyMin.py
+++ b/src/ClassifyMin.py
@@ -75,6 +75,17 @@ def classifyMin_ana(alpha,beta,theta,q3,mu_1,rho_1,print_Cases=False, print_Outp
+# Matrix Version that just gets matrices Q & B
+def classifyMin_mat(Q,B,print_Cases=False, print_Output=False):
+ q1 = Q[0][0]
+ q2 = Q[1][1]
+ q3 = Q[2][2]
+ q12 = Q[0][1]
+ b1 = B[0]
+ b2 = B[1]
+ b3 = B[2]
+ return classifyMin(q1, q2, q3, q12, b1, b2, print_Cases, print_Output)
+
# --------------------------------------------------------------------
# Classify Type of minimizer 1 = R1 , 2 = R2 , 3 = R3 # before : destinction between which axis.. (4Types )
diff --git a/src/HelperFunctions.py b/src/HelperFunctions.py
index b819eaf89ac5da9620d8db27921cec6b1e3f7960..7007d5de2ee64059746a5a4470d0f0ed55391ef1 100644
--- a/src/HelperFunctions.py
+++ b/src/HelperFunctions.py
@@ -94,13 +94,18 @@ def RunCellProblem(alpha,beta,theta,gamma,mu1,rho1, InputFilePath = os.path.dirn
+
+
+
+
+
def GetCellOutput(alpha,beta,theta,gamma,mu1,rho1, InputFilePath = os.path.dirname(os.getcwd()) +"/inputs/computeMuGamma.parset",
OutputFilePath = os.path.dirname(os.getcwd()) + "/outputs/outputMuGamma.txt" ):
RunCellProblem(alpha,beta,theta,gamma,mu1,rho1, InputFilePath)
print('Read effective quantities...')
Q, B = ReadEffectiveQuantities()
- print('Q:', Q)
- print('B:', B)
+ # print('Q:', Q)
+ # print('B:', B)
return Q, B
diff --git a/src/PhaseDiagram.py b/src/PhaseDiagram.py
index 35081eb7b376b71c0878a81987ce96f05d6546d5..3d68b230ee06de24c4d64733163e86a4db08e0ad 100644
--- a/src/PhaseDiagram.py
+++ b/src/PhaseDiagram.py
@@ -118,9 +118,9 @@ beta = 2.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 = 'infinity'
# gamma = 0.5
-# gamma = 0.25
+gamma = 0.25
# gamma = 1.0
print('---- Input parameters: -----')
@@ -133,150 +133,169 @@ print('gamma:', gamma)
print('----------------------------')
# ----------------------------------------------------------------
-# 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
+# gamma_min = 0.05
+# gamma_max = 2.0
- #TODO
-# generalCase = False #Read Output from Cell-Problem instead of using Lemma1.4 (special case)
+gamma_min = 1
+gamma_max = 1
+Gamma_Values = np.linspace(gamma_min, gamma_max, num=1)
-# 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
+# Gamma_Values = np.linspace(gamma_min, gamma_max, num=13) # TODO variable Input Parameters...alpha,beta...
+print('(Input) Gamma_Values:', Gamma_Values)
+for gamma in Gamma_Values:
-# --- 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
+ # 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)
+ # 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
-# 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_2D = 10 # Number of sample points in each direction
-
-
-
-if make_3D_PhaseDiagram:
- alphas_ = np.linspace(-20, 20, SamplePoints_3D)
- betas_ = np.linspace(0.01,40.01,SamplePoints_3D)
- thetas_ = np.linspace(0.01,0.99,SamplePoints_3D)
- # 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!!!
- # print('size of G:', G.shape)
- # print('G:', G)
- # 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)
- # thetas_ = np.linspace(0.01,0.99,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;
- alphas, betas, thetas = np.meshgrid(alphas_, betas_, thetas_, indexing='ij')
-
-
- # if generalCase: #TODO
- # classifyMinVec = np.vectorize(classifyMin)
- # GetCellOutputVec = np.vectorize(GetCellOutput)
- # Q, B = GetCellOutputVec(alpha,betas,thetas,gamma,mu1,rho1,InputFilePath ,OutputFilePath )
- #
- # print('type of Q:', type(Q))
- # print('Q:', Q)
+
+ # --- 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)
#
- # 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!!!
- # 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"
- 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)
- # 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)
- # cbar=plt.colorbar(pnt3d)
- # cbar.set_label("Values (units)")
- ax.set_xlabel('alpha')
- ax.set_ylabel('beta')
- if make_3D_plot: ax.set_zlabel('theta')
- 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('---- 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_2D = 5 # Number of sample points in each direction
+
+
+
+ if make_3D_PhaseDiagram:
+ alphas_ = np.linspace(-20, 20, SamplePoints_3D)
+ betas_ = np.linspace(0.01,40.01,SamplePoints_3D)
+ thetas_ = np.linspace(0.01,0.99,SamplePoints_3D)
+ # 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!!!
+ # print('size of G:', G.shape)
+ # print('G:', G)
+ # 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)
+ # thetas_ = np.linspace(0.01,0.99,SamplePoints_2D)
+ alphas_ = np.linspace(9, 10, SamplePoints_2D)
+ thetas_ = np.linspace(0.075,0.14,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;
+ alphas, betas, thetas = np.meshgrid(alphas_, betas_, thetas_, indexing='ij')
+
+
+ if generalCase: #TODO
+ classifyMin_matVec = np.vectorize(classifyMin_mat)
+ GetCellOutputVec = np.vectorize(GetCellOutput)
+ Q, B = GetCellOutputVec(alphas,betas,thetas,gamma,mu1,rho1,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!!!
+ # 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"
+ 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)
+ # 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)
+ # cbar=plt.colorbar(pnt3d)
+ # cbar.set_label("Values (units)")
+ ax.set_xlabel('alpha')
+ ax.set_ylabel('beta')
+ if make_3D_plot: ax.set_zlabel('theta')
+ 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_)