Plot-Energy_Axial.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 == -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 == -1:
        return r"$-\pi/4$"
    elif N == 2:
        return r"$\pi/2$"
    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: -----')

# q1=1;
# q2=2;
# q12=1/2;
# q3=((4*q1*q2)**0.5-q12)/2;
# # H=[2*q1,q12+2*q3;q12+2*q3,2*q2];
#
# H = np.array([[2*q1, q12+2*q3], [q12+2*q3,2*q2] ])
# A = np.array([[q1,1/2*q12], [1/2*q12,q2] ])
# abar = np.array([q12+2*q3, 2*q2])
# abar = (1.0/math.sqrt((q12+2*q3)**2+(2*q2)**2))*abar
#
# print('abar:',abar)
#
# b = np.linalg.lstsq(A, abar)[0]
# print('b',b)
#
#
# # print('abar:',np.shape(abar))
# # print('np.transpose(abar):',np.shape(np.transpose(abar)))
# sstar = (1/(q1+q2))*abar.dot(A.dot(b))
# # sstar = (1/(q1+q2))*abar.dot(tmp)
# print('sstar', sstar)
# abarperp= np.array([abar[1],-abar[0]])
# print('abarperp:',abarperp)


# -------------------------- Input Parameters --------------------

mu1 = 1.0
rho1 = 1.0
# alpha = 2.0
# theta = 1.0/4.0
beta = 2.0

#set gamma either to 1. '0' 2. 'infinity' or 3. a numerical positive value
gamma = '0'
gamma = 'infinity'


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('----------------------------')
# ----------------------------------------------------------------
print('----------------------------')

# ----------------------------------------------------------------





# Curve parametrised by \theta_rho = alpha in parameter space
N=100;
theta_rho = np.linspace(1, 3, num=N)
print('theta_rho:', theta_rho)


theta_values = []


for t in theta_rho:

        s = (1.0/10.0)*t+0.1
        theta_values.append(s)





theta_rho = np.array(theta_rho)
theta_values = np.array(theta_values)

betas_ = 2.0


# alphas, betas, thetas = np.meshgrid(theta_rho, betas_, theta_values, indexing='ij')
#
#
# harmonicMeanVec = np.vectorize(harmonicMean)
# arithmeticMeanVec = np.vectorize(arithmeticMean)
# prestrain_b1Vec = np.vectorize(prestrain_b1)
# prestrain_b2Vec = np.vectorize(prestrain_b2)
#
# GetMuGammaVec = np.vectorize(GetMuGamma)
# muGammas = GetMuGammaVec(betas,thetas,gamma,mu1,rho1,InputFilePath ,OutputFilePath )
#
# q1_vec = harmonicMeanVec(mu1, betas, thetas)
# q2_vec = arithmeticMeanVec(mu1, betas, thetas)
#
# b1_vec = prestrain_b1Vec(rho1, betas, alphas, thetas)
# b2_vec = prestrain_b2Vec(rho1, betas, alphas, thetas)

# special case: q12 == 0!! .. braucht eigentlich nur b1 & b2 ...

# print('type b1_values:', type(b1_values))



# print('size(q1)',q1.shape)


energy_axial1 = []
energy_axial2 = []

# for b1 in b1_values:
for i in range(len(theta_rho)):
    print('index i:', i)

    print('theta_rho[i]',theta_rho[i])
    print('theta_values[i]',theta_values[i])

    q1 = (1.0/6.0)*harmonicMean(mu1, beta, theta_values[i])
    q2 = (1.0/6.0)*arithmeticMean(mu1, beta, theta_values[i])
    q12 = 0.0
    q3 = GetMuGamma(beta, theta_values[i],gamma,mu1,rho1,InputFilePath ,OutputFilePath )
    b1 = prestrain_b1(rho1,beta, theta_rho[i],theta_values[i] )
    b2 = prestrain_b2(rho1,beta, theta_rho[i],theta_values[i] )


    # q2_vec = arithmeticMean(mu1, betas, thetas)
    #
    # b1_vec = prestrain_b1Vec(rho1, betas, alphas, thetas)
    # b2_vec = prestrain_b2Vec(rho1, betas, alphas, thetas)
    print('q1[i]',q1)
    print('q2[i]',q2)
    print('q3[i]',q3)
    print('b1[i]',b1)
    print('b2[i]',b2)
    # print('q1[i]',q1[0][i])
    # print('q2[i]',q2[i])
    # print('b1[i]',b1[i])
    # print('b2[i]',b2[i])
    #compute axial energy #1 ...

    a_axial1 = np.array([b1,0])
    a_axial2 = np.array([0,b2])
    b = np.array([b1,b2])

    H = np.array([[2*q1, q12+2*q3], [q12+2*q3,2*q2] ])
    A = np.array([[q1,1/2*q12], [1/2*q12,q2] ])


    tmp = H.dot(a_axial1)

    print('H',H)
    print('A',A)
    print('b',b)
    print('a_axial1',a_axial1)
    print('tmp',tmp)

    tmp = (1/2)*a_axial1.dot(tmp)
    print('tmp',tmp)

    tmp2 = A.dot(b)
    print('tmp2',tmp2)
    tmp2 = 2*a_axial1.dot(tmp2)

    print('tmp2',tmp2)
    energy_1 = tmp - tmp2
    print('energy_1',energy_1)


    energy_axial1.append(energy_1)


    tmp = H.dot(a_axial2)

    print('H',H)
    print('A',A)
    print('b',b)
    print('a_axial2',a_axial2)
    print('tmp',tmp)

    tmp = (1/2)*a_axial2.dot(tmp)
    print('tmp',tmp)

    tmp2 = A.dot(b)
    print('tmp2',tmp2)
    tmp2 = 2*a_axial2.dot(tmp2)

    print('tmp2',tmp2)
    energy_2 = tmp - tmp2
    print('energy_2',energy_2)


    energy_axial2.append(energy_2)





print('theta_values', theta_values)

mpl.rcParams['text.usetex'] = True
mpl.rcParams["font.family"] = "serif"
mpl.rcParams["font.size"] = "9"
width = 6.28 *0.5
# width = 6.28
height = width / 1.618
fig = plt.figure()
ax = plt.axes((0.17,0.21 ,0.75,0.75))
# ax = plt.axes((0.15,0.18,0.8,0.8))
# 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.xaxis.set_major_locator(plt.MultipleLocator(np.pi / 8))
# ax.xaxis.set_minor_locator(plt.MultipleLocator(np.pi / 16))
# ax.xaxis.set_major_locator(plt.MultipleLocator(np.pi / 2))
# ax.xaxis.set_minor_locator(plt.MultipleLocator(np.pi / 4))
# ax.xaxis.set_major_formatter(plt.FuncFormatter(format_func))
ax.grid(True,which='major',axis='both',alpha=0.3)




# ax.plot(theta_rho, theta_values, 'royalblue', zorder=3, )


ax.plot(theta_rho, energy_axial1, 'royalblue', zorder=3, label=r"axialMin1")
ax.plot(theta_rho, energy_axial2, 'forestgreen', zorder=3, label=r"axialMin2")
# ax.plot(-1.0*alphas, kappas, 'red', zorder=3, )




lg = ax.legend(bbox_to_anchor=(0.0, 0.75), loc='upper left')




#
#
#
#
# kappas = []
# alphas = []
# # G.append(float(s[0]))
#
#
#
#
# for t in T :
#
#     abar_current = sstar*abar+t*abarperp;
#     # print('abar_current', abar_current)
#     abar_current[abar_current < 1e-10] = 0
#     # print('abar_current', abar_current)
#
#     # G = np.array([[2*q1, q12+2*q3], [q12+2*q3,2*q2] ])
#     G = [abar_current[0], abar_current[1] , (2*abar_current[0]*abar_current[1])**0.5 ]
#
#     e = [(abar_current[0]/(abar_current[0]+abar_current[1]))**0.5, (abar_current[1]/(abar_current[0]+abar_current[1]))**0.5]
#     kappa = abar_current[0]+abar_current[1]
#     alpha = math.atan2(e[1], e[0])
#
#     print('angle current:', alpha)
#
#     kappas.append(kappa)
#     alphas.append(alpha)
#
#
#
# alphas = np.array(alphas)
# kappas = np.array(kappas)
#
#
# print('kappas:',kappas)
# print('alphas:',alphas)
# print('min alpha:', min(alphas))
# print('min kappa:', min(kappas))
#
# 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.xaxis.set_major_locator(plt.MultipleLocator(np.pi / 8))
# # ax.xaxis.set_minor_locator(plt.MultipleLocator(np.pi / 16))
# ax.xaxis.set_major_locator(plt.MultipleLocator(np.pi / 2))
# ax.xaxis.set_minor_locator(plt.MultipleLocator(np.pi / 4))
# ax.xaxis.set_major_formatter(plt.FuncFormatter(format_func))
# ax.grid(True,which='major',axis='both',alpha=0.3)
#
#
#
#
# ax.plot(alphas, kappas, 'royalblue', zorder=3, )
# ax.plot(-1.0*alphas, kappas, 'red', zorder=3, )



ax.set_xlabel(r"$\theta_\rho$")
# ax.set_ylabel(r"$\theta$")
ax.set_ylabel(r"energy")


# ax.set_xticks([-np.pi/2, -np.pi/4  ,0,  np.pi/4,  np.pi/2 ])
# labels = ['$0$',r'$\pi/8$', r'$\pi/4$' ,r'$3\pi/8$' , r'$\pi/2$']
# ax.set_yticklabels(labels)





# ax.legend(loc='upper right')



fig.set_size_inches(width, height)
fig.savefig('Plot-Energy_Axial.pdf')

plt.show()