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 time
from ClassifyMin import *
# from scipy.io import loadmat #not Needed anymore?
import codecs
import sys
def ReadEffectiveQuantities(QFilePath = os.path.dirname(os.getcwd()) + '/outputs/QMatrix.txt', BFilePath = os.path.dirname(os.getcwd())+ '/outputs/BMatrix.txt'):
# Read Output Matrices (effective quantities)
# From Cell-Problem output Files : ../outputs/Qmatrix.txt , ../outputs/Bmatrix.txt
# -- Read Matrix Qhom
X = []
with codecs.open(path + '/outputs/QMatrix.txt', encoding='utf-8-sig') as f:
for line in f:
s = line.split()
X.append([float(s[i]) for i in range(len(s))])
Q = np.array([[X[0][2], X[1][2], X[2][2]],
[X[3][2], X[4][2], X[5][2]],
[X[6][2], X[7][2], X[8][2]] ])
# -- Read Beff (as Vector)
X = []
with codecs.open(path + '/outputs/BMatrix.txt', encoding='utf-8-sig') as f:
for line in f:
s = line.split()
X.append([float(s[i]) for i in range(len(s))])
B = np.array([X[0][2], X[1][2], X[2][2]])
return Q, B
def RunCellProblem(alpha,beta,theta,gamma,mu1,rho1, InputFilePath = os.path.dirname(os.getcwd()) +"/inputs/computeMuGamma.parset"):
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
# Optional: Check Time
# t = time.time()
subprocess.run(['./build-cmake/src/Cell-Problem', './inputs/cellsolver.parset'],
capture_output=True, text=True)
# print('elapsed time:', time.time() - t)
# --------------------------------------------------------------------------------------
def Compare_Classification(alpha,beta,theta,gamma,mu1,rho1, InputFilePath = os.path.dirname(os.getcwd()) +"/inputs/computeMuGamma.parset"):
# ---------------------------------------------------------------
# Comparison of the analytical Classification 'ClassifyMin'
# and the symbolic Minimizatio + Classification 'symMinimization'
# ----------------------------------------------------------------
comparison_successful = True
eps = 1e-8
# 1. Substitute Input-Parameters for the Cell-Problem
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()
# 2. --- Run Cell-Problem
print('Run Cell-Problem ...')
# Optional: Check Time
# t = time.time()
subprocess.run(['./build-cmake/src/Cell-Problem', './inputs/cellsolver.parset'],
capture_output=True, text=True)
# 3. --- 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
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)
# 4. --- Read the effective quantities (output values of the Cell-Problem)
# Read Output Matrices (effective quantities)
print('Read effective quantities...')
Q, B = ReadEffectiveQuantities()
print('Q:', Q)
print('B:', B)
q1 = Q[0][0]
q2 = Q[1][1]
q3 = Q[2][2]
q12 = Q[0][1]
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)
# --- Check Assumptions:
# Assumption of Classification-Lemma1.6: [b3 == 0] & [Q orthotropic]
# Check if b3 is close to zero
assert (b3 < eps), "ClassifyMin only defined for b3 == 0 !"
# Check if Q is orthotropic i.e. q13 = q31 = q23 = q32 == 0
assert(Q[2][0] < eps and Q[0][2] < eps and Q[1][2] < eps and Q[2][1] < eps), "Q is not orthotropic !"
# 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)
# 6. Compare
# print('DifferenceMatrix:', G_ana - G )
# print('MinimizerError (FrobeniusNorm):', np.linalg.norm(G_ana - G , 'fro') )
if np.linalg.norm(G_ana - G , 'fro') > eps :
comparison_successful = False
print('Difference between Minimizers is too large !')
if type != type_ana :
comparison_successful = False
print('Difference in type !')
if abs(angle-angle_ana) > eps :
comparison_successful = False
print('Difference in angle is too large!')
if abs(kappa-kappa_ana) > eps :
comparison_successful = False
print('Difference in curvature is too large!')
if comparison_successful:
print('Comparison successful')
return comparison_successful
# ----------------------------------------------------------------------------------------------------------------------------
# ----- Setup Paths -----
InputFile = "/inputs/cellsolver.parset"
OutputFile = "/outputs/output.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)
#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)
#
# #-------------------- PLOT OPTION -------------------------------------------
#
# Gamma_Values = np.linspace(0.01, 2.5, num=6) # 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]
#
#
#
# # Make Plot
# if make_Plot:
# plt.figure()
# plt.title(r'$\mu_\gamma(\gamma)$-Plot')
# plt.plot(Gamma_Values, mu_gamma)
# plt.scatter(Gamma_Values, mu_gamma)
# # 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.xlabel("$\gamma$")
# plt.ylabel("$\mu_\gamma$")
# plt.legend()
# plt.show()
#
# # ------------------------------------------------------------------------
print('---- Input parameters: -----')
mu1 = 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('----------------------------')
# print('RunCellProblem...')
# RunCellProblem(alpha,beta,theta,gamma,mu1,rho1,InputFilePath)
# TEST read Matrix file
# Test = loadmat(path + '/outputs/QMatrix.mat')
print('Compare_Classification...')
Compare_Classification(alpha,beta,theta,gamma,mu1,rho1,InputFilePath)
# ------------- 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)