Add option for 2D PhaseDiagram

src/PhaseDiagram.py

@@ -33,10 +33,15 @@ theta = 1.0/8.0
 # print_Cases = True
 # print_Output = True
 
+# make_3D_plot = True
+# make_2D_plot = False
+# make_3D_PhaseDiagram = True
+# make_2D_PhaseDiagram = False
+
 make_3D_plot = False
-make_2D_plot = False
-make_3D_PhaseDiagram = True
-make_2D_PhaseDiagram = False
+make_2D_plot = True
+make_3D_PhaseDiagram = False
+make_2D_PhaseDiagram = True
 
 
 # --- Define q1, q2 , mu_gamma, q12 ---
@@ -45,10 +50,14 @@ make_2D_PhaseDiagram = False
 q1 = harmonicMean(mu_1, beta, theta)
 q2 = arithmeticMean(mu_1, beta, theta)
 
-# TEST
+# TEST / TODO read from Cell-Output
 q12 = 0.0  # (analytical example)
-# --- Set mu_gamma  to value or read from Cell-Output
+
+# TODO read from Cell-Output
+# --- Set mu_gamma  to value
 mu_gamma = q1   # TODO read from Cell-Output
+
+
 b1 = prestrain_b1(rho_1, beta, alpha, theta)
 b2 = prestrain_b2(rho_1, beta, alpha, theta)
 
@@ -75,12 +84,13 @@ print('----------------------------')
 
 # ---------------------- MAKE PLOT / Write to VTK------------------------------------------------------------------------------
 
+SamplePoints_3D = 100 # Number of sample points in each direction
+SamplePoints_2D = 800 # Number of sample points in each direction
+
 if make_3D_PhaseDiagram:
-    alphas_ = np.linspace(-20, 20, 80)
-    # betas_  = np.linspace(0.01,40.01,1)
-    #fix to one value:
-    betas_ = 2.0;
-    thetas_ = np.linspace(0.01,0.99,80)
+    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')
@@ -96,46 +106,58 @@ if make_3D_PhaseDiagram:
     print('Written to VTK-File:', VTKOutputName )
 
 if make_2D_PhaseDiagram:
-    alphas_ = np.linspace(-20, 20, 80)
+    alphas_ = np.linspace(-20, 20, SamplePoints_2D)
     # betas_  = np.linspace(0.01,40.01,1)
     #fix to one value:
     betas_ = 2.0;
-    beta = np.array([betas_])
-    thetas_ = np.linspace(0.01,0.99,80)
+    thetas_ = np.linspace(0.01,0.99,SamplePoints_2D)
     # print('type of alphas', type(alphas_))
     # print('Test:', type(np.array([mu_gamma])) )
-    alphas, thetas = np.meshgrid(alphas_, thetas_, indexing='ij')
+    alphas, betas, thetas = np.meshgrid(alphas_, betas_, thetas_, indexing='ij')
     classifyMin_geoVec = np.vectorize(classifyMin_geo)
-    G, angles, Types, curvature = classifyMin_geoVec(alphas,beta,thetas,mu_gamma,q12, mu_1, rho_1)
+    G, angles, Types, curvature = classifyMin_geoVec(alphas,betas,thetas,mu_gamma,q12, mu_1, rho_1)
     # print('size of G:', G.shape)
     # print('G:', G)
     # Out = classifyMin_geoVec(alphas,betas,thetas)
     # T = Out[2]
     # --- Write to VTK
-    # VTKOutputName = "./PhaseDiagram2D"
-    # gridToVTK(VTKOutputName , alphas, thetas, beta, pointData = {'Type': Types, 'angles': angles, 'curvature': curvature} )
-    # print('Written to VTK-File:', VTKOutputName )
+    VTKOutputName = "./PhaseDiagram2D"
+    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')
-    # pnt3d=ax.scatter(x,y,z,c=z.flat)
     colors = cm.plasma(Types)
-    # pnt3d=ax.scatter(x,y,z,c=T.flat)
     if make_2D_plot: pnt3d=ax.scatter(alphas,thetas,c=Types.flat)
+    # if make_2D_plot: plt.scatter(alphas,thetas,c=Types.flat)
     if make_3D_plot: pnt3d=ax.scatter(alphas,betas,thetas,c=Types.flat)
-    # pnt3d=ax.scatter(x,y,z,c=colors.flat)
-    cbar=plt.colorbar(pnt3d)
-    cbar.set_label("Values (units)")
+    # cbar=plt.colorbar(pnt3d)
+    # cbar.set_label("Values (units)")
     ax.set_xlabel('alpha')
     ax.set_ylabel('beta')
-    ax.set_zlabel('theta')
+    if make_3D_plot: ax.set_zlabel('theta')
     plt.show()
     # plt.savefig('common_labels.png', dpi=300)
-
-    # print('TEST 3D Scatter')
     # print('T:', T)
     # print('Type 1 occured here:', np.where(T == 1))
     # print('Type 2 occured here:', np.where(T == 2))
+
+
+# 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()