Add Script for AngleCurv-Gamma Subplots

src/Plot-AngleCurvature-GammaV2_SubPlots.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
+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
+
+angles_1 = []
+angles_2 = []
+angles_3 = []
+
+beta = 2.0
+theta= 0.25
+
+
+
+
+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...
+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_1 = classifyMin_anaVec(alpha_1, beta, theta, muGammas,  mu1, rho1)
+G, angles_2, Types, curvature_2 = classifyMin_anaVec(alpha_2, beta, theta, muGammas,  mu1, rho1)
+G, angles_3, Types, curvature_3 = 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)
+
+print('curvature_1:', curvature_1)
+print('curvature_2:', curvature_2)
+print('curvature_3:', curvature_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)
+
+curvature_1 = np.array(curvature_1)
+curvature_2 = np.array(curvature_2)
+curvature_3 = np.array(curvature_3)
+
+# ---------------- Create Plot -------------------
+# plt.figure()
+
+# Styling
+# plt.style.use("seaborn-darkgrid")
+# plt.rcParams["font.family"] = "Avenir"
+# plt.rcParams["font.size"] = 16
+
+
+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(figsize=(width,height))
+
+# fig,ax = plt.subplots(nrows=2,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.tight_layout()
+#
+#
+# fig = plt.figure()
+
+gs = fig.add_gridspec(nrows=2,ncols=3, hspace=0.15, wspace=0.1)
+# ax = gs.subplots(sharey=True)
+
+# Create Three Axes Objects
+
+
+ax4 = fig.add_subplot(gs[1, 0])
+ax5 = fig.add_subplot(gs[1, 1],sharey=ax4)
+ax6 = fig.add_subplot(gs[1, 2],sharey=ax4)
+plt.setp(ax5.get_yticklabels(), visible=False)
+plt.setp(ax6.get_yticklabels(), visible=False)
+
+ax1 = fig.add_subplot(gs[0, 0],sharex=ax4)
+ax2 = fig.add_subplot(gs[0, 1],sharey=ax1)
+ax3 = fig.add_subplot(gs[0, 2],sharey=ax1)
+plt.setp(ax1.get_xticklabels(), visible=False)
+plt.setp(ax2.get_xticklabels(), visible=False)
+plt.setp(ax3.get_xticklabels(), visible=False)
+plt.setp(ax2.get_yticklabels(), visible=False)
+plt.setp(ax3.get_yticklabels(), visible=False)
+
+# 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,0].yaxis.set_major_locator(plt.MultipleLocator(np.pi / 8))
+# ax[0,0].yaxis.set_minor_locator(plt.MultipleLocator(np.pi / 16))
+# ax[0,0].yaxis.set_major_formatter(plt.FuncFormatter(format_func))
+# ax[0,1].yaxis.set_major_locator(plt.MultipleLocator(np.pi / 8))
+# ax[0,1].yaxis.set_minor_locator(plt.MultipleLocator(np.pi / 16))
+# ax[0,1].yaxis.set_major_formatter(plt.FuncFormatter(format_func))
+# ax[0,2].yaxis.set_major_locator(plt.MultipleLocator(np.pi / 8))
+# ax[0,2].yaxis.set_minor_locator(plt.MultipleLocator(np.pi / 16))
+# ax[0,2].yaxis.set_major_formatter(plt.FuncFormatter(format_func))
+#
+# ax[0,0].grid(True,which='major',axis='both',alpha=0.3)
+# ax[0,1].grid(True,which='major',axis='both',alpha=0.3)
+# ax[0,2].grid(True,which='major',axis='both',alpha=0.3)
+
+ax1.yaxis.set_major_locator(plt.MultipleLocator(np.pi / 8))
+ax1.yaxis.set_minor_locator(plt.MultipleLocator(np.pi / 16))
+ax1.yaxis.set_major_formatter(plt.FuncFormatter(format_func))
+ax2.yaxis.set_major_locator(plt.MultipleLocator(np.pi / 8))
+ax2.yaxis.set_minor_locator(plt.MultipleLocator(np.pi / 16))
+ax2.yaxis.set_major_formatter(plt.FuncFormatter(format_func))
+ax3.yaxis.set_major_locator(plt.MultipleLocator(np.pi / 8))
+ax3.yaxis.set_minor_locator(plt.MultipleLocator(np.pi / 16))
+ax3.yaxis.set_major_formatter(plt.FuncFormatter(format_func))
+
+ax1.grid(True,which='major',axis='both',alpha=0.3)
+ax2.grid(True,which='major',axis='both',alpha=0.3)
+ax3.grid(True,which='major',axis='both',alpha=0.3)
+
+
+ax1.plot(Gamma_Values, angles_1, 'royalblue', zorder=3, )
+ax2.plot(Gamma_Values, angles_2, 'royalblue', zorder=3, )
+ax3.plot(Gamma_Values, angles_3, 'royalblue', zorder=3, )
+
+# ax1.set_xlabel(r"$\gamma$")
+ax1.set_ylabel(r"angle  $\alpha$")
+ax1.xaxis.set_minor_locator(MultipleLocator(0.25))
+ax1.xaxis.set_major_locator(MultipleLocator(0.5))
+
+# ax2.set_xlabel(r"$\gamma$")
+ax2.xaxis.set_minor_locator(MultipleLocator(0.25))
+ax2.xaxis.set_major_locator(MultipleLocator(0.5))
+
+# ax3.set_xlabel(r"$\gamma$")
+# ax[2].set_ylabel(r"angle  $\alpha$")
+ax3.xaxis.set_minor_locator(MultipleLocator(0.25))
+ax3.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$"]
+labels = ['$0$',r'$\pi/8$', r'$\pi/4$' ,r'$3\pi/8$' , r'$\pi/2$']
+ax1.set_yticks([0, np.pi/8, np.pi/4, 3*np.pi/8 , np.pi/2, ])
+ax2.set_yticks([0, np.pi/8, np.pi/4, 3*np.pi/8 , np.pi/2 ])
+ax3.set_yticks([0, np.pi/8, np.pi/4, 3*np.pi/8 , np.pi/2 ])
+
+ax1.set_yticklabels(labels)
+ax2.set_yticklabels(labels)
+ax3.set_yticklabels(labels)
+
+ax1.set_ylim([0-0.1, np.pi/2+0.1])
+ax2.set_ylim([0-0.1, np.pi/2+0.1])
+ax3.set_ylim([0-0.1, np.pi/2+0.1])
+
+# for i in range(3):
+#     ax1[i].set_ylim([0-0.1, np.pi/2+0.1])
+
+# Plot Gamma Value that is closest to q3_star
+l1 = ax1.axvline(x = gammaClose, color = 'midnightblue', linestyle = 'dashed', linewidth=1, label='$\gamma^*$')
+l2 = ax2.axvline(x = gammaClose, color = 'midnightblue', linestyle = 'dashed', linewidth=1, label='$\gamma^*$')
+l3 = ax3.axvline(x = gammaClose, color = 'midnightblue', linestyle = 'dashed', linewidth=1, label='$\gamma^*$')
+
+
+# -------------------------------------------------------------------------------------
+
+# ax[1,0].grid(True,which='major',axis='both',alpha=0.3)
+# ax[1,1].grid(True,which='major',axis='both',alpha=0.3)
+# ax[1,2].grid(True,which='major',axis='both',alpha=0.3)
+
+# ax[1,0].set_xlabel(r"$\gamma$")
+# ax[1,0].set_ylabel(r"curvature $\kappa$")
+# ax[1,0].xaxis.set_minor_locator(MultipleLocator(0.5))
+# ax[1,0].xaxis.set_major_locator(MultipleLocator(1))
+# ax[1,0].yaxis.set_minor_locator(MultipleLocator(0.5))
+# ax[1,0].yaxis.set_major_locator(MultipleLocator(1))
+# ax[1,1].set_xlabel(r"$\gamma$")
+# # ax[1].set_ylabel(r"angle  $\alpha$")
+# ax[1,1].xaxis.set_minor_locator(MultipleLocator(0.5))
+# ax[1,1].xaxis.set_major_locator(MultipleLocator(1))
+# ax[1,2].set_xlabel(r"$\gamma$")
+# # ax[2].set_ylabel(r"angle  $\alpha$")
+# ax[1,2].xaxis.set_minor_locator(MultipleLocator(0.5))
+# ax[1,2].xaxis.set_major_locator(MultipleLocator(1))
+# l4 = ax[1,0].axvline(x = gammaClose, color = 'midnightblue', linestyle = 'dashed', linewidth=1, label='$\gamma^*$', zorder=4)
+# l5 = ax[1,1].axvline(x = gammaClose, color = 'midnightblue', linestyle = 'dashed', linewidth=1, label='$\gamma^*$', zorder=4)
+# l6 = ax[1,2].axvline(x = gammaClose, color = 'midnightblue', linestyle = 'dashed', linewidth=1, label='$\gamma^*$' ,zorder=4)
+
+
+
+ax4.grid(True,which='major',axis='both',alpha=0.3)
+ax5.grid(True,which='major',axis='both',alpha=0.3)
+ax6.grid(True,which='major',axis='both',alpha=0.3)
+ax4.plot(Gamma_Values, curvature_1, 'forestgreen', zorder=3, )
+ax5.plot(Gamma_Values, curvature_2, 'forestgreen', zorder=3, )
+ax6.plot(Gamma_Values, curvature_3, 'forestgreen', zorder=3, )
+# ax2.plot(Gamma_Values, curvature_1, 'forestgreen', zorder=3, )
+# ax2.plot(Gamma_Values, curvature_2, 'forestgreen', zorder=3, )
+# ax2.plot(Gamma_Values, curvature_3, 'forestgreen', zorder=3, )
+ax4.set_xlabel(r"$\gamma$")
+ax4.set_ylabel(r"curvature $\kappa$")
+# ax4.set_ylabel(r"curvature $\kappa$", labelpad=10)
+ax4.xaxis.set_minor_locator(MultipleLocator(0.25))
+ax4.xaxis.set_major_locator(MultipleLocator(0.5))
+# ax4.yaxis.set_minor_locator(MultipleLocator(0.1))
+ax4.yaxis.set_major_locator(MultipleLocator(0.05))
+ax5.set_xlabel(r"$\gamma$")
+# ax[1].set_ylabel(r"angle  $\alpha$")
+ax5.xaxis.set_minor_locator(MultipleLocator(0.25))
+ax5.xaxis.set_major_locator(MultipleLocator(0.5))
+ax6.set_xlabel(r"$\gamma$")
+# ax[2].set_ylabel(r"angle  $\alpha$")
+ax6.xaxis.set_minor_locator(MultipleLocator(0.25))
+ax6.xaxis.set_major_locator(MultipleLocator(0.5))
+l4 = ax4.axvline(x = gammaClose, color = 'midnightblue', linestyle = 'dashed', linewidth=1, label='$\gamma^*$', zorder=4)
+l5 = ax5.axvline(x = gammaClose, color = 'midnightblue', linestyle = 'dashed', linewidth=1, label='$\gamma^*$', zorder=4)
+l6 = ax6.axvline(x = gammaClose, color = 'midnightblue', linestyle = 'dashed', linewidth=1, label='$\gamma^*$' ,zorder=4)
+
+
+
+
+
+#
+#
+
+
+
+## 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.align_ylabels()
+
+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()
+#