From 78ce2695149ef23adb5ce19c9c5a7410344861d5 Mon Sep 17 00:00:00 2001
From: Klaus <klaus.boehnlein@tu-dresden.de>
Date: Sat, 18 Sep 2021 00:12:16 +0200
Subject: [PATCH] Fix error in generalized form of phase Diagram
---
Matlab-Programs/classifyMIN.m | 2 +-
src/PhaseDiagram.py | 269 +++++++++++++++++-----------------
2 files changed, 136 insertions(+), 135 deletions(-)
diff --git a/Matlab-Programs/classifyMIN.m b/Matlab-Programs/classifyMIN.m
index 3f36a2c4..d73a7692 100755
--- a/Matlab-Programs/classifyMIN.m
+++ b/Matlab-Programs/classifyMIN.m
@@ -54,8 +54,8 @@ b2 = (3*rho_1./(4.*((1-t)+t.*b))).*(1-t.*(1+b.*a));
CaseCount = 0; %check if more than one case ever happens
+epsilon = eps;
-epsilon = 1.e-18;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% PARABOLIC CASE %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% if abs(det(A)) < epsilon * min(abs(det(A)),0)
diff --git a/src/PhaseDiagram.py b/src/PhaseDiagram.py
index 3d68b230..51ccbc0b 100644
--- a/src/PhaseDiagram.py
+++ b/src/PhaseDiagram.py
@@ -137,14 +137,14 @@ print('----------------------------')
# gamma_min = 0.05
# gamma_max = 2.0
-gamma_min = 1
-gamma_max = 1
-Gamma_Values = np.linspace(gamma_min, gamma_max, num=1)
-
-# 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:
+# gamma_min = 1
+# gamma_max = 1
+# Gamma_Values = np.linspace(gamma_min, gamma_max, num=1)
+#
+# # 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:
# muGamma = GetMuGamma(beta,theta,gamma,mu1,rho1,InputFilePath)
@@ -155,147 +155,148 @@ for gamma in Gamma_Values:
# 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
-
+ #TODO
+generalCase = True #Read Output from Cell-Problem instead of using Lemma1.4 (special case)
- # --- 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
+# 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_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')
+# 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 = 1 # 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, otypes=[np.ndarray, np.ndarray])
+ Q, B = GetCellOutputVec(alphas,betas,thetas,gamma,mu1,rho1,InputFilePath ,OutputFilePath )
+
+ # output = output.astype(object)
+
+ # print('type of Q:', type(Q))
+ # print('Q:', Q)
+ G, angles, Types, curvature = classifyMin_matVec(Q,B)
+
+ else:
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!!!
+ 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
- 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))
+ # 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_)
--
GitLab