diff --git a/src/HelperFunctions.py b/src/HelperFunctions.py index 426c87a2e7710f5b7d557a6bc5bfe18f7347ce3e..578fcfdbb7115d783b3ebe31babbf0a95db4a4fd 100644 --- a/src/HelperFunctions.py +++ b/src/HelperFunctions.py @@ -62,6 +62,63 @@ def RunCellProblem(alpha,beta,theta,gamma,mu1,rho1, InputFilePath = os.path.dirn + +# 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 + + + + + + def Compare_Classification(alpha,beta,theta,gamma,mu1,rho1, InputFilePath = os.path.dirname(os.getcwd()) +"/inputs/computeMuGamma.parset"): # --------------------------------------------------------------- # Comparison of the analytical Classification 'ClassifyMin' @@ -93,7 +150,7 @@ def Compare_Classification(alpha,beta,theta,gamma,mu1,rho1, InputFilePath = os.p # 3. --- Run Matlab symbolic minimization program: 'symMinimization' eng = matlab.engine.start_matlab() # s = eng.genpath(path + '/Matlab-Programs') - s = eng.genpath(path) + s = eng.genpath(os.path.dirname(os.getcwd())) eng.addpath(s, nargout=0) # print('current Matlab folder:', eng.pwd(nargout=1)) eng.cd('Matlab-Programs', nargout=0) #switch to Matlab-Programs folder @@ -123,13 +180,13 @@ def Compare_Classification(alpha,beta,theta,gamma,mu1,rho1, InputFilePath = os.p b1 = B[0] b2 = B[1] b3 = B[2] - # print("q1:", q1) - # print("q2:", q2) - # print("q3:", q3) - # print("q12:", q12) - # print("b1:", b1) - # print("b2:", b2) - # print("b3:", b3) + print("q1:", q1) + print("q2:", q2) + print("q3:", q3) + print("q12:", q12) + print("b1:", b1) + print("b2:", b2) + print("b3:", b3) # --- Check Assumptions: # Assumption of Classification-Lemma1.6: [b3 == 0] & [Q orthotropic] @@ -141,11 +198,11 @@ def Compare_Classification(alpha,beta,theta,gamma,mu1,rho1, InputFilePath = os.p # 5. --- Get output from the analytical Classification 'ClassifyMin' G_ana, angle_ana, type_ana, kappa_ana = classifyMin(q1, q2, q3, q12, b1, b2) - # print('Minimizer G_ana:') - # print(G_ana) - # print('angle_ana:', angle_ana) - # print('type_ana:', type_ana ) - # print('curvature_ana:', kappa_ana) + print('Minimizer G_ana:') + print(G_ana) + print('angle_ana:', angle_ana) + print('type_ana:', type_ana ) + print('curvature_ana:', kappa_ana) # 6. Compare # print('DifferenceMatrix:', G_ana - G ) @@ -167,6 +224,8 @@ def Compare_Classification(alpha,beta,theta,gamma,mu1,rho1, InputFilePath = os.p if comparison_successful: print('Comparison successful') + else: + print('Comparison unsuccessful') return comparison_successful diff --git a/src/makePlot.py b/src/makePlot.py new file mode 100644 index 0000000000000000000000000000000000000000..b73422667bafad7cd4b39c25ce46d72b30a9d74c --- /dev/null +++ b/src/makePlot.py @@ -0,0 +1,236 @@ +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 +from HelperFunctions import * +# from subprocess import Popen, PIPE +#import sys + +###################### makePlot.py ######################### +# Generalized Plot-Script giving the option to define +# quantity of interest and the parameter it depends on +# to create a plot +# +# Input: Define y & x for "x-y plot" as Strings +# - Run the 'Cell-Problem' for the different Parameter-Points +# (alternatively run 'Compute_MuGamma' if quantity of interest +# is q3=muGamma for a significant Speedup) + +########################################################### + +# TODO +# - Fallunterscheidung (Speedup) falls gesuchter value mu_gamma = q3 +# - Also Add option to plot Minimization Output + + +# ----- Setup Paths ----- +InputFile = "/inputs/cellsolver.parset" +OutputFile = "/outputs/output.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) + +#--------------------------------------------------------------- + +print('---- Input parameters: -----') +mu1 = 10.0 +lambda1 = 10.0 +rho1 = 1.0 +alpha = 2.8 +beta = 2.0 +theta = 1.0/4.0 +gamma = 0.75 + +print('mu1: ', mu1) +print('rho1: ', rho1) +print('alpha: ', alpha) +print('beta: ', beta) +print('theta: ', theta) +print('gamma:', gamma) +print('----------------------------') + + +# TODO? : Ask User for Input ... +# function = input("Enter value you want to plot (y-value):\n") +# print(f'You entered {function}') +# parameter = input("Enter Parameter this value depends on (x-value) :\n") +# print(f'You entered {parameter}') + +# Add Option to change NumberOfElements used for computation of Cell-Problem + + +# --- Define Quantity of interest: +# Options: 'q1', 'q2', 'q3', 'q12' ,'q21', 'q31', 'q13' , 'q23', 'q32' , 'b1', 'b2' ,'b3' +# TODO: EXTRA (MInimization Output) 'Minimizer (norm?)' 'angle', 'type', 'curvature' +yName = 'q12' +# yName = 'b1' +yName = 'q3' +yName = 'angle' + +# --- Define Parameter this function/quantity depends on: +# Options: mu1 ,lambda1, rho1 , alpha, beta, theta, gamma +# xName = 'theta' +xName = 'gamma' +# xName = 'lambda1' + +# --- define Interval of x-values: +# xmin = 0.0 +# xmax = 10.0 +xmin = 0.01 +xmax = 3.0 +numPoints = 3 +X_Values = np.linspace(xmin, xmax, num=numPoints) +print(X_Values) + + +Y_Values = [] + + + +# --- Options +RUN = True +# RUN = False +# make_Plot = False +make_Plot = True + +if RUN: + for x in X_Values: + + if yName =='q3' or yName == 'mu_gamma': + with open(path+"/inputs/computeMuGamma.parset", 'r') as file: + filedata = file.read() + filedata = re.sub('(?m)^'+xName+'=.*',xName+'='+str(x),filedata) + f = open(path+"/inputs/computeMuGamma.parset",'w') + f.write(filedata) + f.close() + #Run Compute_MuGamma since its much faster + print("Run Compute_MuGamma for " + xName + " =" , x) + subprocess.run(['./build-cmake/src/Compute_MuGamma', './inputs/computeMuGamma.parset'], + capture_output=True, text=True) + #Extract mu_gamma from Output-File + with open(path + '/outputs/outputMuGamma.txt', '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) + Y_Values.append(float(s[0])) + else : + with open(InputFilePath, 'r') as file: + filedata = file.read() + filedata = re.sub('(?m)^'+xName+'=.*',xName+'='+str(x),filedata) + f = open(InputFilePath,'w') + f.write(filedata) + f.close() + # Run Cell-Problem + print("Run Cell-Problem for " + xName + " =" , x) + subprocess.run(['./build-cmake/src/Cell-Problem', './inputs/cellsolver.parset'], + capture_output=True, text=True) + #Extract quantity from Output-File + Q, B = ReadEffectiveQuantities() + if yName == 'q1': # TODO: Better use dictionary?... + print('q1 used') + Y_Values.append(Q[0][0]) + elif yName =='q2': + print('q2 used') + Y_Values.append(Q[1][1]) + # elif yName =='q3': + # print('q3 used') + # Y_Values.append(Q[0][0]) + elif yName =='q12': + print('q12 used') + Y_Values.append(Q[0][1]) + elif yName =='q21': + print('q21 used') + Y_Values.append(Q[1][0]) + elif yName =='q13': + print('q13 used') + Y_Values.append(Q[0][2]) + elif yName =='q31': + print('q31 used') + Y_Values.append(Q[2][0]) + elif yName =='q23': + print('q23 used') + Y_Values.append(Q[1][2]) + elif yName =='q32': + print('q32 used') + Y_Values.append(Q[2][1]) + elif yName =='b1': + print('b1 used') + Y_Values.append(B[0]) + elif yName =='b2': + print('b2 used') + Y_Values.append(B[1]) + elif yName =='b3': + print('b3 used') + Y_Values.append(B[2]) + elif yName == 'angle' or yName =='type' or yName =='curvature': + # ------------- Run Matlab symbolic minimization program 'symMinimization' + 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 + Inp_T = True + print('Run symbolic Minimization...') + #Arguments: symMinization(print_Input,print_statPoint,print_Output,make_FunctionPlot, InputPath) + 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) + if yName =='angle': + print('angle used') + Y_Values.append(angle) + if yName =='type': + print('angle used') + Y_Values.append(angle) + if yName =='kappa': + print('angle used') + Y_Values.append(angle) + # ------------end of for-loop ----------------- + print("(Output) Values of " + yName + ": ", Y_Values) +# ----------------end of if-statement ------------- + + + + + + + + +# ---------------- Create Plot ------------------- +plt.figure() +# plt.title(r''+ yName + '-Plot') +plt.plot(X_Values, Y_Values) +plt.scatter(X_Values, Y_Values) +# plt.axis([0, 6, 0, 20]) +plt.xlabel(xName) +plt.ylabel(yName) + + +# plt.axhline(y = 1.90476, color = 'b', linestyle = ':', label='$q_1$') +# plt.axhline(y = 2.08333, color = 'r', linestyle = 'dashed', label='$q_2$') +# plt.legend() +plt.show() +# #---------------------------------------------------------------