From 5203d9dde69146092cdcc03f4f0a23b8ab1ddc03 Mon Sep 17 00:00:00 2001
From: Klaus <klaus.boehnlein@tu-dresden.de>
Date: Mon, 13 Sep 2021 20:11:02 +0200
Subject: [PATCH] Add option to plot angle over mu_gamma
---
src/Plotq3-Angle.py | 246 ++++++++++++++++++++++++++++++++++++++++++++
1 file changed, 246 insertions(+)
create mode 100644 src/Plotq3-Angle.py
diff --git a/src/Plotq3-Angle.py b/src/Plotq3-Angle.py
new file mode 100644
index 00000000..fe6b7df6
--- /dev/null
+++ b/src/Plotq3-Angle.py
@@ -0,0 +1,246 @@
+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 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
+# from subprocess import Popen, PIPE
+#import sys
+
+# unabhängig von alpha...
+def GetMuGamma(beta,theta,gamma,mu1,rho1, InputFilePath = os.path.dirname(os.getcwd()) +"/inputs/computeMuGamma.parset"):
+ # ------------------------------------ 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
+
+
+
+
+InputFile = "/inputs/computeMuGamma.parset"
+OutputFile = "/outputs/outputMuGamma.txt"
+# path = os.getcwd()
+# InputFilePath = os.getcwd()+InputFile
+# OutputFilePath = os.getcwd()+OutputFile
+# --------- 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)
+
+
+#1. Define Inputs Gamma-Array..
+#2. for(i=0; i<length(array)) ..compute Q_hom, B_eff from Cell-Problem
+#3
+
+# matrix = np.loadtxt(path + 'Qmatrix.txt', usecols=range(3))
+# print(matrix)
+
+# Use Shell Commands:
+# subprocess.run('ls', shell=True)
+
+
+#---------------------------------------------------------------
+
+# -------------------------- Input Parameters --------------------
+mu1 = 10.0
+rho1 = 1.0
+alpha = 10 #1.05263158
+beta = 40.0 #5.0
+# theta = 1.0/4.0
+theta = 1.0/8.0 # 0.5
+
+# InterestingParameterSet :
+# mu1 = 10.0
+# rho1 = 1.0
+# alpha = 10
+# beta = 40.0
+# theta = 1.0/8.0
+
+
+#set gamma either to 1. '0' 2. 'infinity' or 3. a numerical positive value
+# gamma = '0'
+# gamma = 'infinity'
+# gamma = 0.5
+
+print('---- Input parameters: -----')
+print('mu1: ', mu1)
+print('rho1: ', rho1)
+print('alpha: ', alpha)
+print('beta: ', beta)
+print('theta: ', theta)
+# print('gamma:', gamma)
+print('----------------------------')
+# ----------------------------------------------------------------
+
+
+
+
+
+Gamma_Values = np.linspace(0.01, 5, num=15) # TODO variable Input Parameters...alpha,beta...
+print('(Input) Gamma_Values:', Gamma_Values)
+# mu_gamma = []
+
+
+#
+# # --- Options
+# RUN = True
+# # RUN = False
+# # make_Plot = False
+make_Plot = True # vll besser : Plot_muGamma
+#
+# if RUN:
+# for gamma in Gamma_Values:
+# print("Run Cell-Problem for Gamma = ", gamma)
+# # print('gamma='+str(gamma))
+# with open(InputFilePath, 'r') as file:
+# filedata = file.read()
+# filedata = re.sub('(?m)^gamma=.*','gamma='+str(gamma),filedata)
+# f = open(InputFilePath,'w')
+# f.write(filedata)
+# f.close()
+# # --- Run Cell-Problem
+# 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_gammaValue = float(s[0])
+# print("mu_gamma:", mu_gammaValue)
+# mu_gamma.append(mu_gammaValue)
+# # ------------end of for-loop -----------------
+# print("(Output) Values of mu_gamma: ", mu_gamma)
+# # ----------------end of if-statement -------------
+#
+# # mu_gamma=[2.06099, 1.90567, 1.905]
+# # mu_gamma=[2.08306, 1.905, 1.90482, 1.90479, 1.90478, 1.90477]
+#
+# ##Gamma_Values = np.linspace(0.01, 20, num=12) :
+# #mu_gamma= [2.08306, 1.91108, 1.90648, 1.90554, 1.90521, 1.90505, 1.90496, 1.90491, 1.90487, 1.90485, 1.90483, 1.90482]
+#
+# ##Gamma_Values = np.linspace(0.01, 2.5, num=12)
+# # mu_gamma=[2.08306, 2.01137, 1.96113, 1.93772, 1.92592, 1.91937, 1.91541, 1.91286, 1.91112, 1.90988, 1.90897, 1.90828]
+#
+# # Gamma_Values = np.linspace(0.01, 2.5, num=6)
+# # mu_gamma=[2.08306, 1.95497, 1.92287, 1.91375, 1.9101, 1.90828]
+
+
+GetMuGammaVec = np.vectorize(GetMuGamma)
+muGammas = GetMuGammaVec(beta,theta,Gamma_Values,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)
+
+# Make Plot
+if make_Plot:
+ 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()
+
+#
+# # ------------- RUN Matlab symbolic minimization program
+# eng = matlab.engine.start_matlab()
+# # s = eng.genpath(path + '/Matlab-Programs')
+# s = eng.genpath(path)
+# eng.addpath(s, nargout=0)
+# # print('current Matlab folder:', eng.pwd(nargout=1))
+# eng.cd('Matlab-Programs', nargout=0) #switch to Matlab-Programs folder
+# # print('current Matlab folder:', eng.pwd(nargout=1))
+# Inp = False
+# print('Run symbolic Minimization...')
+# G, angle, type, kappa = eng.symMinimization(Inp,Inp,Inp,Inp, nargout=4) #Name of program:symMinimization
+# # G, angle, type, kappa = eng.symMinimization(Inp,Inp,Inp,Inp,path + "/outputs", nargout=4) #Optional: add Path
+# G = np.asarray(G) #cast Matlab Outout to numpy array
+#
+# # --- print Output ---
+# print('Minimizer G:')
+# print(G)
+# print('angle:', angle)
+# print('type:', type )
+# print('curvature:', kappa)
--
GitLab