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 ClassifyMin 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
lambda1 = 0.0
gamma = 1.0/4.0
gamma = 'infinity'
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'
# yName = 'curvature'
# --- Define Parameter this function/quantity depends on:
# Options: mu1 ,lambda1, rho1 , alpha, beta, theta, gamma
# xName = 'theta'
# xName = 'gamma'
# xName = 'lambda1'
xName = 'theta'
# --- define Interval of x-values:
xmin = 0.05
xmax = 0.2
# xmin = 0.01
# xmax = 3.0
numPoints = 5
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 xName != 'gamma' and (gamma=='0' or gamma=='infinity') and lambda1 ==0.0
#
# print('Situation of Lemma1.4')
# q12 = 0.0
# q1 = (1.0/6.0)*harmonicMean(mu_1, beta, theta)
# q2 = (1.0/6.0)*arithmeticMean(mu_1, beta, theta)
# b1 = prestrain_b1(rho_1, beta, alpha,theta)
# b2 = prestrain_b2(rho_1, beta, alpha,theta)
# b3 = 0.0
#
#
#
# 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')
# if gamma == '0':
# Y_Values.append(q2)
# if gamma == 'infinity':
# Y_Values.append(q1)
# elif yName =='q12':
# print('q12 used')
# Y_Values.append(q12)
# elif yName =='b1':
# print('b1 used')
# Y_Values.append(b1)
# elif yName =='b2':
# print('b2 used')
# Y_Values.append(b2)
# elif yName =='b3':
# print('b3 used')
# Y_Values.append(b3)
# elif yName == 'angle' or yName =='type' or yName =='curvature':
# G, angle, Type, curvature = classifyMin_anaVec(alphas,betas,thetas, muGammas, mu1, rho1)
# if yName =='angle':
# print('angle used')
# Y_Values.append(angle)
# if yName =='type':
# print('angle used')
# Y_Values.append(type)
# if yName =='kappa':
# print('angle used')
# Y_Values.append(curvature)
#
#
if yName =='q3' or yName == 'mu_gamma': # if only q3 is needed .. compute 2D problem (much faster)
# fist replace old parameters in parset
SetParametersComputeMuGamma(beta,theta,gamma,mu1,rho1, InputFilePath)
# replace the sought after x value in the parset
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 : # Run full Cell-Problem...
# fist replace old parameters in parset
SetParametersCellProblem(alpha,beta,theta,gamma,mu1,rho1,lambda1,InputFilePath)
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, curvature = 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:', curvature)
if yName =='angle':
print('angle used')
Y_Values.append(angle)
if yName =='type':
print('type used')
Y_Values.append(type)
if yName =='curvature':
print('curvature used')
Y_Values.append(curvature)
# ------------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()
# #---------------------------------------------------------------