Update plot routines

src/PhaseDiagram.py

@@ -115,11 +115,11 @@ mu1 = 1.0
 rho1 = 1.0
 alpha = 2.0
 beta = 2.0
-beta = 10.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 = 'infinity'
 # gamma = 0.5
 # gamma = 0.25
 # gamma = 1.0
@@ -159,8 +159,9 @@ print('----------------------------')
 
 print('type of gamma:', type(gamma))
 # # #
-Gamma_Values = ['0', 'infinity']
-# Gamma_Values = ['infinity']
+# Gamma_Values = ['0', 'infinity']
+Gamma_Values = ['infinity']
+Gamma_Values = ['0']
 print('(Input) Gamma_Values:', Gamma_Values)
 
 for gamma in Gamma_Values:
@@ -180,13 +181,13 @@ for gamma in Gamma_Values:
     generalCase = False
 
     # make_3D_plot = True
-    make_3D_PhaseDiagram = True
+    # make_3D_PhaseDiagram = True
     make_2D_plot = False
     make_2D_PhaseDiagram = False
     make_3D_plot = False
-    # make_3D_PhaseDiagram = False
+    make_3D_PhaseDiagram = False
     # make_2D_plot = True
-    # make_2D_PhaseDiagram = True
+    make_2D_PhaseDiagram = True
     #
 
     # --- Define effective quantities: q1, q2 , q3 = mu_gamma, q12 ---
@@ -232,7 +233,7 @@ for gamma in Gamma_Values:
     # 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 = 200 # 4000 # Number of sample points in each direction
 
     if make_3D_PhaseDiagram:
         alphas_ = np.linspace(-20, 20, SamplePoints_3D)
@@ -307,7 +308,11 @@ for gamma in Gamma_Values:
         # betas_  = np.linspace(0.01,40.01,1)
         #fix to one value:
         # betas_ = 2.0;
-        betas_ = 10.0;
+        # betas_ = 10.0;
+        betas_ = 5.0;
+
+
+
 
         # TEST
         # alphas_ = np.linspace(-8, 8, SamplePoints_2D)

src/PhaseDiagram_intermediateGamma.py

+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
+
+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
+rho1 = 1.0
+alpha = 2.0
+beta = 2.0
+# beta = 10.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
+
+
+
+
+
+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...
+
+Gamma_Values = [0.05, 0.25, 0.5, 0.75, 1.0, 1.25, 1.5]
+# Gamma_Values = ['infinity']
+print('(Input) Gamma_Values:', Gamma_Values)
+# #
+for gamma in Gamma_Values:
+
+
+        # 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_2D = 200 # 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(-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)
+
+        #bEFORE PrestrainChange:
+        # alphas_ = np.linspace(8, 10, SamplePoints_2D)
+        # thetas_ = np.linspace(0.05,0.16,SamplePoints_2D)
+
+        alphas_ = np.linspace(-2, 1, SamplePoints_2D)
+        thetas_ = np.linspace(0.4,0.6,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;
+        alphas, betas, thetas = np.meshgrid(alphas_, betas_, 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!!!
+            # 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)")
+        plt.axvline(x = 8, color = 'b', linestyle = ':', label='$q_1$')
+        plt.axhline(y = 0.083333333, color = 'b', linestyle = ':', label='$q_1$')
+
+        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_)
+        # 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()