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 import matplotlib.ticker as ticker import matplotlib as mpl from matplotlib.ticker import MultipleLocator,FormatStrFormatter,MaxNLocator import pandas as pd # from matplotlib import rc # rc('text', usetex=True) # Use LaTeX font # # import seaborn as sns # sns.set(color_codes=True) def format_func(value, tick_number): # # find number of multiples of pi/2 # N = int(np.round(2 * value / np.pi)) # if N == 0: # return "0" # elif N == 1: # return r"$\pi/2$" # elif N == 2: # return r"$\pi$" # elif N % 2 > 0: # return r"${0}\pi/2$".format(N) # else: # return r"${0}\pi$".format(N // 2) # find number of multiples of pi/2 N = int(np.round(4 * value / np.pi)) if N == 0: return "0" elif N == 1: return r"$\pi/4$" elif N == 2: return r"$\pi/2$" elif N % 2 > 0: return r"${0}\pi/2$".format(N) else: return r"${0}\pi$".format(N // 2) def find_nearest(array, value): array = np.asarray(array) idx = (np.abs(array - value)).argmin() return array[idx] def find_nearestIdx(array, value): array = np.asarray(array) idx = (np.abs(array - value)).argmin() return idx 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) print('---- Input parameters: -----') alpha = 10 mu1 = 1.0 rho1 = 1.0 beta = 2.0 #5.0 theta = 1.0/8.0 # alpha = -0.5 beta = 40.0 theta= 1/8.0 # # INTERESTING! from pi/2: alpha = -0.5 beta = 40.0 theta= 1/8.0 # # # # INTERESTING! from pi/2: # alpha = -0.2 # beta = 25.0 # theta= 1/2 # INTERESTING!: # alpha = -0.5 # beta = 5.0 # theta= 1/30 # INTERESTING!: # alpha = -0.25 # beta = 10.0 # theta= 3/4 # # INTERESTING!: alpha = -0.25 beta = 10.0 theta= 1/8 # # INTERESTING!: # alpha = -0.25 # beta = 5.0 # theta= 1/8 # # # INTERESTING!: alpha = -0.5 beta = 10.0 theta= 1/8 alpha_1 = -1.0 alpha_2 = -0.75 alpha_3 = -0.70 # alpha_1 = -1 # alpha_2 = -0.5 # alpha_3 = -0.25 angles_1 = [] angles_2 = [] angles_3 = [] beta = 2.0 theta= 0.25 # beta = 10.0 # theta= 0.5 print('mu1: ', mu1) print('rho1: ', rho1) print('alpha_1: ', alpha_1) print('alpha_2: ', alpha_2) print('alpha_3: ', alpha_3) print('beta: ', beta) print('theta: ', theta) # print('gamma:', gamma) print('----------------------------') # ---------------------------------------------------------------- gamma_min = 0.01 gamma_max = 1.5 # Gamma_Values = np.linspace(gamma_min, gamma_max, num=100) # TODO variable Input Parameters...alpha,beta... Gamma_Values = np.linspace(gamma_min, gamma_max, num=40) # TODO variable Input Parameters...alpha,beta... print('(Input) Gamma_Values:', Gamma_Values) # mu_gamma = [] # Gamma_Values = '0' # Get values for mu_Gamma GetMuGammaVec = np.vectorize(GetMuGamma) muGammas = GetMuGammaVec(beta,theta,Gamma_Values,mu1,rho1, InputFilePath ,OutputFilePath ) print('muGammas:', muGammas) q12 = 0.0 q1 = (1.0/6.0)*harmonicMean(mu1, beta, theta) q2 = (1.0/6.0)*arithmeticMean(mu1, beta, theta) print('q1: ', q1) print('q2: ', q2) b1 = prestrain_b1(rho1, beta, alpha,theta) b2 = prestrain_b2(rho1, beta, alpha,theta) q3_star = math.sqrt(q1*q2) print('q3_star:', q3_star) # TODO these have to be compatible with input parameters!!! # compute certain ParameterValues that this makes sense # b1 = q3_star # b2 = q1 print('b1: ', b1) print('b2: ', b2) # return classifyMin(q1, q2, q3, q12, b1, b2, print_Cases, print_Output) # classifyMin_anaVec = np.vectorize(classifyMin_ana) # G, angles, Types, curvature = classifyMin_anaVec(alpha, beta, theta, muGammas, mu1, rho1) classifyMin_anaVec = np.vectorize(classifyMin_ana) G, angles_1, Types, curvature = classifyMin_anaVec(alpha_1, beta, theta, muGammas, mu1, rho1) G, angles_2, Types, curvature = classifyMin_anaVec(alpha_2, beta, theta, muGammas, mu1, rho1) G, angles_3, Types, curvature = classifyMin_anaVec(alpha_3, beta, theta, muGammas, mu1, rho1) # _,angles,_,_ = classifyMin_anaVec(alpha, beta, theta, muGammas, mu1, rho1) print('angles_1:', angles_1) print('angles_2:', angles_2) print('angles_3:', angles_3) idx = find_nearestIdx(muGammas, q3_star) print('GammaValue Idx closest to q_3^*', idx) gammaClose = Gamma_Values[idx] print('GammaValue(Idx) with mu_gamma closest to q_3^*', gammaClose) determinantVec = np.vectorize(determinant) detValues = determinantVec(q1,q2,muGammas,q12) print('detValues:', detValues) detZeroidx = find_nearestIdx(detValues, 0) print('idx where det nearest to zero', idx) gammaClose = Gamma_Values[detZeroidx] print('gammaClose:', gammaClose) # --- Convert to numpy array Gamma_Values = np.array(Gamma_Values) angles_1 = np.array(angles_1) angles_2 = np.array(angles_2) angles_3 = np.array(angles_3) # ---------------- Create Plot ------------------- # plt.figure() mpl.rcParams['text.usetex'] = True mpl.rcParams["font.family"] = "serif" mpl.rcParams["font.size"] = "9" width = 6.28 height = width / 1.618 height = width / 2.5 fig = plt.figure() fig,ax = plt.subplots(nrows=1,ncols=3,figsize=(width,height)) # more than one plot # fig,ax = plt.subplots(nrows=1,ncols=3,figsize=(width,height),sharey=True) # Share Y-axis fig = plt.figure() gs = fig.add_gridspec(1,3, hspace=0.2, wspace=0.1) ax = gs.subplots(sharey=True) # ax = plt.axes((0.15,0.21 ,0.75,0.75)) # ax = plt.axes((0.15,0.21 ,0.8,0.75)) # ax.tick_params(axis='x',which='major', direction='out',pad=5) # ax.tick_params(axis='y',which='major', length=3, width=1, direction='out',pad=3) # ax.xaxis.set_major_locator(MultipleLocator(0.1)) # ax.xaxis.set_minor_locator(MultipleLocator(0.05)) ax[0].yaxis.set_major_locator(plt.MultipleLocator(np.pi / 8)) ax[0].yaxis.set_minor_locator(plt.MultipleLocator(np.pi / 16)) ax[0].yaxis.set_major_formatter(plt.FuncFormatter(format_func)) ax[1].yaxis.set_major_locator(plt.MultipleLocator(np.pi / 8)) ax[1].yaxis.set_minor_locator(plt.MultipleLocator(np.pi / 16)) ax[1].yaxis.set_major_formatter(plt.FuncFormatter(format_func)) ax[2].yaxis.set_major_locator(plt.MultipleLocator(np.pi / 8)) ax[2].yaxis.set_minor_locator(plt.MultipleLocator(np.pi / 16)) ax[2].yaxis.set_major_formatter(plt.FuncFormatter(format_func)) ax[0].grid(True,which='major',axis='both',alpha=0.3) ax[1].grid(True,which='major',axis='both',alpha=0.3) ax[2].grid(True,which='major',axis='both',alpha=0.3) # ax.grid(True,which='major',axis='both',alpha=0.3) # # plt.title(r'angle$-\mu_\gamma(\gamma)$-Plot') # # plt.title(r'angle$-\gamma$-Plot') # plt.plot(Gamma_Values, angles) # plt.scatter(Gamma_Values, angles) # plt.plot(muGammas, angles) # plt.scatter(muGammas, angles) # # plt.axis([0, 6, 0, 20]) # # plt.axhline(y = 1.90476, color = 'b', linestyle = ':', label='$q_1$') # # plt.axhline(y = 2.08333, color = 'r', linestyle = 'dashed', label='$q_2$') # plt.axvline(x = q1, color = 'b', linestyle = ':', label='$q_1$') # plt.axvline(x = q2, color = 'r', linestyle = 'dashed', label='$q_2$') # # plt.axvline(x = q3_star, color = 'r', linestyle = 'dashed', label='$\gamma^*$') ax[0].plot(Gamma_Values, angles_1, 'royalblue', zorder=3, ) ax[1].plot(Gamma_Values, angles_2, 'royalblue', zorder=3, ) ax[2].plot(Gamma_Values, angles_3, 'royalblue', zorder=3, ) ax[0].set_xlabel(r"$\gamma$") ax[0].set_ylabel(r"angle $\alpha$") ax[0].xaxis.set_minor_locator(MultipleLocator(0.25)) ax[0].xaxis.set_major_locator(MultipleLocator(0.5)) ax[1].set_xlabel(r"$\gamma$") # ax[1].set_ylabel(r"angle $\alpha$") ax[1].xaxis.set_minor_locator(MultipleLocator(0.25)) ax[1].xaxis.set_major_locator(MultipleLocator(0.5)) ax[2].set_xlabel(r"$\gamma$") # ax[2].set_ylabel(r"angle $\alpha$") ax[2].xaxis.set_minor_locator(MultipleLocator(0.25)) ax[2].xaxis.set_major_locator(MultipleLocator(0.5)) # Labels to use in the legend for each line line_labels = [r"$\theta_\mu = 1.0$", r"$\theta_\mu = 2.0$", r"$\theta_\mu = 5.0$", r"$\theta_\mu = 10.0$"] # ax.plot(Gamma_Values, angles, 'royalblue', zorder=3, ) # ax.set_xlabel(r"$\gamma$") # ax.set_ylabel(r"angle $\alpha$") # plt.xlabel("$q_3$") # plt.xlabel("$\gamma$") # plt.ylabel("angle") # ax.grid(True) # ax.yaxis.set_major_locator(plt.MultipleLocator(np.pi / 2)) # ax.yaxis.set_minor_locator(plt.MultipleLocator(np.pi / 12)) # ax.yaxis.set_major_locator(plt.MultipleLocator(np.pi / 4)) # ax.yaxis.set_minor_locator(plt.MultipleLocator(np.pi / 12)) # ax.yaxis.set_major_formatter(plt.FuncFormatter(format_func)) # ax.yaxis.set_major_formatter(ticker.FormatStrFormatter('%g $\pi$')) # ax.yaxis.set_major_locator(ticker.MultipleLocator(base=0.25)) # ax.yaxis.set_major_formatter(ticker.FuncFormatter( # lambda val,pos: '{:.0g}$\pi$'.format(2*val/np.pi) if val !=0 else '0')) # ax.yaxis.set_major_locator(ticker.MultipleLocator(base=0.5*np.pi)) # ---------------------------- show pi values ------------------------------------ # ax.axvline(x = q1, color = 'b', linestyle = ':', label='$q_1$') # ax.axvline(x = q2, color = 'r', linestyle = 'dashed', label='$q_2$') # ax.legend() # # ax.set(xlim=(1.750, 1.880), ylim=(0, math.pi/2.0)) # ax.set(xlim=(1.760, 1.880), ylim=(-0.1, np.pi/4.0)) # # ax.set_yticks([0, np.pi/4 ,np.pi/2]) # # labels = ['$0$', r'$\pi/4$', r'$\pi/2$'] # ax.set_yticks([0, np.pi/8, np.pi/4 ]) # labels = ['$0$',r'$\pi/8$', r'$\pi/4$'] # ax.set_yticklabels(labels) # --------------------------------------------------------------- # ax.set(xlim=(1.750, 1.880), ylim=(0, math.pi/2.0)) # ax.set(xlim=(1.760, 1.880), ylim=(-0.1, np.pi/4.0)) # ax.set(xlim=(1.760, 1.880), ylim=(-0.1, np.pi/4.0)) # ax.set_yticks([0, np.pi/4 ,np.pi/2]) # labels = ['$0$', r'$\pi/4$', r'$\pi/2$'] # OLD : # ax.set_yticks([0, np.pi/8, np.pi/4 ]) # labels = ['$0$',r'$\pi/8$', r'$\pi/4$'] labels = ['$0$',r'$\pi/8$', r'$\pi/4$' ,r'$3\pi/8$' , r'$\pi/2$'] ax[0].set_yticks([0, np.pi/8, np.pi/4, 3*np.pi/8 , np.pi/2, ]) ax[1].set_yticks([0, np.pi/8, np.pi/4, 3*np.pi/8 , np.pi/2 ]) ax[2].set_yticks([0, np.pi/8, np.pi/4, 3*np.pi/8 , np.pi/2 ]) ax[0].set_yticklabels(labels) ax[1].set_yticklabels(labels) ax[2].set_yticklabels(labels) for i in range(3): ax[i].set_ylim([0-0.1, np.pi/2+0.1]) # Plot Gamma Value that is closest to q3_star l1 = ax[0].axvline(x = gammaClose, color = 'midnightblue', linestyle = 'dashed', linewidth=1, label='$\gamma^*$') l2 = ax[1].axvline(x = gammaClose, color = 'midnightblue', linestyle = 'dashed', linewidth=1, label='$\gamma^*$') l3 = ax[2].axvline(x = gammaClose, color = 'midnightblue', linestyle = 'dashed', linewidth=1, label='$\gamma^*$') # color elliptic/hyperbolic region # ax.axvspan(gamma_min, gammaClose, color='red', alpha=0.2) # ax.axvspan(gammaClose, gamma_max, color='green', alpha=0.2) # ax[0].legend(loc='upper right') # ax[2].legend(l3, # labels=['$\gamma^*$'], # loc='upper left', # bbox_to_anchor=(1,1)) # ax[2].legend( # loc='upper left', # bbox_to_anchor=(1,1)) # fig.legend([l1], # The line objects # # label='$\gamma^*$', # The labels for each line # loc="center right", # Position of legend # borderaxespad=0.15 # Small spacing around legend box # # title="Legend Title" # Title for the legend # ) line_labels = [r"$\gamma^*$"] fig.legend([l1], [r"$\gamma^*$"], # bbox_to_anchor=[0.5, 0.92], bbox_to_anchor=[0.5, 0.94], loc='center', ncol=3) # plt.subplots_adjust(wspace=0.4, hspace=0.0) # plt.tight_layout() # Adjust the scaling factor to fit your legend text completely outside the plot # (smaller value results in more space being made for the legend) # plt.subplots_adjust(right=0.9) plt.subplots_adjust(bottom=0.2) fig.set_size_inches(width, height) fig.savefig('Plot-Angle-Gamma.pdf') plt.show() # plt.figure() # plt.title(r'angle$-\mu_\gamma(\gamma)$-Plot') # plt.plot(muGammas, angles) # plt.scatter(muGammas, angles) # # plt.axis([0, 6, 0, 20]) # # plt.axhline(y = 1.90476, color = 'b', linestyle = ':', label='$q_1$') # # plt.axhline(y = 2.08333, color = 'r', linestyle = 'dashed', label='$q_2$') # plt.axvline(x = 1.90476, color = 'b', linestyle = ':', label='$q_1$') # plt.axvline(x = 2.08333, color = 'r', linestyle = 'dashed', label='$q_2$') # plt.xlabel("$\mu_\gamma$") # plt.ylabel("angle") # plt.legend() # plt.show() #