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: -----') # mu1 = 10.0 # rho1 = 1.0 # alpha = 2.56140350877193 #2.56140350877193, 4.0852130325814535 # beta = 2.0 #5.0 # theta = 1.0/4.0 # theta = 1.0/8.0 # 0.5 # theta = 0.075 # 0.5 # mu1 = 10.0 # rho1 = 1.0 alpha = 10.0 # beta = 40.0 # theta = 0.125 # mu1 = 10.0 mu1 = 1.0 rho1 = 1.0 # alpha = 2.0 beta = 2.0 #5.0 # theta = 1.0/4.0 theta = 1.0/8.0 # # mu1 = 10.0 # rho1 = 1.0 # alpha = 10.0 # beta = 2.0 #5.0 # theta = 1.0/8.0 # # # mu1 = 10.0 # rho1 = 1.0 # # alpha = 10.02021333 # alpha = 10.0 # beta = 2.0 # theta = 0.125 # mu1 = 10.0 # rho1 = 1.0 # # alpha = 10.02021333 # alpha = 9.0 # beta = 2.0 # theta = 0.075 # # mu1 = 10.0 # rho1 = 1.0 # alpha = 4.0 # beta = 10.0 # theta = 1.0/4.0 # alpha = 10 #1.05263158 # beta = 40.0 #5.0 # # theta = 1.0/4.0 # theta = 1.0/8.0 # 0.5 # # # alpha = 2.0 # beta = 2.0 #5.0 # theta = 1/4.0 # rho1 = 1.0 # mu1 = 10.0 # InterestingParameterSet : # mu1 = 10.0 # rho1 = 1.0 # alpha = 10 # beta = 40.0 # theta = 0.124242 # gamma = 0.75 #set gamma either to 1. '0' 2. 'infinity' or 3. a numerical positive value # gamma = '0' # gamma = 'infinity' # gamma = 0.5 # gamma = 0.25 print('mu1: ', mu1) print('rho1: ', rho1) print('alpha: ', alpha) print('beta: ', beta) print('theta: ', theta) # print('gamma:', gamma) print('----------------------------') # ---------------------------------------------------------------- gamma_min = 0.01 gamma_max = 1.0 Gamma_Values = np.linspace(gamma_min, gamma_max, num=20) # 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, Types, curvature = classifyMin_anaVec(alpha, beta, theta, muGammas, mu1, rho1) # _,angles,_,_ = classifyMin_anaVec(alpha, beta, theta, muGammas, mu1, rho1) print('angles:', angles) 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 = np.array(angles) # ---------------- Create Plot ------------------- # plt.figure() mpl.rcParams['text.usetex'] = True mpl.rcParams["font.family"] = "serif" mpl.rcParams["font.size"] = "9" width = 6.28 *0.5 height = width / 1.618 fig = plt.figure() # 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.yaxis.set_major_locator(plt.MultipleLocator(np.pi / 8)) ax.yaxis.set_minor_locator(plt.MultipleLocator(np.pi / 16)) ax.yaxis.set_major_formatter(plt.FuncFormatter(format_func)) 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^*$') # f,ax=plt.subplots(1) # ax.plot(muGammas, angles) # ax.scatter(muGammas, angles) ax.plot(Gamma_Values, angles, 'royalblue', zorder=3, ) # ax.scatter(Gamma_Values, angles) # ax.set_xlabel(r"$q_3(\gamma)$") ax.set_xlabel(r"$\gamma$") # ax.set_ylabel(r"angle $\angle$") 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$'] ax.set_yticklabels(labels) # Plot Gamma Value that is closest to q3_star ax.axvline(x = gammaClose, color = 'midnightblue', linestyle = 'dashed', 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.legend(loc='upper right') # plt.xlabel("$q_3(\gamma)$") # plt.xlabel("$\gamma$") # plt.ylabel("angle") # plt.legend(loc='upper center') 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() #