From f785b70a1d99a2f79f0a4e5c6b502542c86d0583 Mon Sep 17 00:00:00 2001
From: Klaus <klaus.boehnlein@tu-dresden.de>
Date: Tue, 14 Sep 2021 21:38:47 +0200
Subject: [PATCH] Add helper functions to Run the C++ programs, Read Output
from Files & compareClassifications
---
src/Cell-Problem.cc | 13 +-
src/CellScript.py | 341 ++++++++++++++++++++++++++++++++---------
src/ClassifyMin.py | 14 +-
src/HelperFunctions.py | 176 +++++++++++++++++++++
4 files changed, 461 insertions(+), 83 deletions(-)
create mode 100644 src/HelperFunctions.py
diff --git a/src/Cell-Problem.cc b/src/Cell-Problem.cc
index 2657bbc0..bc152182 100644
--- a/src/Cell-Problem.cc
+++ b/src/Cell-Problem.cc
@@ -1448,13 +1448,14 @@ int main(int argc, char *argv[])
std::cout << "Theta: " << theta << std::endl;
std::cout << "Gamma: " << gamma << std::endl;
std::cout << "Number of Elements in each direction: " << nElements << std::endl;
- log << "q1: " << q1 << std::endl;
- log << "q2: " << q2 << std::endl;
- log << "q3: " << q3 << std::endl;
+ log << "q1=" << q1 << std::endl;
+ log << "q2=" << q2 << std::endl;
+ log << "q3=" << q3 << std::endl;
+ log << "q12=" << Q[0][1] << std::endl;
- log << "effective b1: " << Beff[0] << std::endl;
- log << "effective b2: " << Beff[1] << std::endl;
- log << "effective b3: " << Beff[2] << std::endl;
+ log << "b1=" << Beff[0] << std::endl;
+ log << "b2=" << Beff[1] << std::endl;
+ log << "b3=" << Beff[2] << std::endl;
log << "b1_hat: " << B_hat[0] << std::endl;
log << "b2_hat: " << B_hat[1] << std::endl;
log << "b3_hat: " << B_hat[2] << std::endl;
diff --git a/src/CellScript.py b/src/CellScript.py
index ca11ec26..19cec12e 100644
--- a/src/CellScript.py
+++ b/src/CellScript.py
@@ -8,18 +8,174 @@ import fileinput
import re
import matlab.engine
import time
-# from subprocess import Popen, PIPE
-#import sys
+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"
-# 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)
@@ -42,79 +198,120 @@ print("Path: ", path)
# 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()
+#
+# # ------------------------------------------------------------------------
-#---------------------------------------------------------------
-Gamma_Values = np.linspace(0.01, 2.5, num=6) # TODO variable Input Parameters...alpha,beta...
-print('(Input) Gamma_Values:', Gamma_Values)
-mu_gamma = []
+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('----------------------------')
-# --- Options
-RUN = True
-# RUN = False
-# make_Plot = False
-make_Plot = True
-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 -> much faster!!!
-
- # subprocess.run(['./build-cmake/src/Cell-Problem_muGamma', './inputs/cellsolver.parset'],
- # capture_output=True, text=True)
- print('elapsed time:', time.time() - t)
- #Extract mu_gamma from Output-File
- with open(OutputFilePath, 'r') as file:
- output = file.read()
- tmp = re.search(r'(?m)^mu_gamma=.*',output).group()
- 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('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
diff --git a/src/ClassifyMin.py b/src/ClassifyMin.py
index 0c84fdab..9f7e7d22 100644
--- a/src/ClassifyMin.py
+++ b/src/ClassifyMin.py
@@ -40,9 +40,13 @@ def f(a1, a2, q1, q2, q3, q12, b1, b2):
# ---- Alternative Version using alpha,beta,theta ,mu_1,rho_1
def classifyMin_ana(alpha,beta,theta,q3,mu_1,rho_1,print_Cases=False, print_Output=False):
- # TODO: assert(q12 == 0!)?
+ # Assumption of Classification-Lemma1.6:
+ # 1. [b3 == 0]
+ # 2. Q is orthotropic i.e. q13 = q31 = q23 = q32 == 0
+ # 3. This additionally assumes that Poisson-Ratio = 0 => q12 == 0
q12 = 0.0
+
q1 = (1.0/6.0)*harmonicMean(mu_1, beta, theta)
q2 = (1.0/6.0)*arithmeticMean(mu_1, beta, theta)
# print('q1: ', q1)
@@ -51,8 +55,8 @@ def classifyMin_ana(alpha,beta,theta,q3,mu_1,rho_1,print_Cases=False, print_Outp
b2 = prestrain_b2(rho_1, beta, alpha,theta)
return classifyMin(q1, q2, q3, q12, b1, b2, print_Cases, print_Output)
-
-# Classify Type of minimizer 1 = R1 , 2 = R2 , 3 = R3 # before : destinction between which axis.. (4Types )
+# --------------------------------------------------------------------
+# Classify Type of minimizer 1 = R1 , 2 = R2 , 3 = R3 # before : destinction between which axis.. (4Types )
# where
# R1 : unique local (global) minimizer which is not axial
# R2 : continuum of local (global) minimizers which are not axial
@@ -61,8 +65,8 @@ def classifyMin_ana(alpha,beta,theta,q3,mu_1,rho_1,print_Cases=False, print_Outp
# R1 = E1
# R2 = P1.2
# R3 = E2 U E3 U P1.1 U P2 U H
-# -------------------------------------------------------------
-def classifyMin(q1, q2, q3, q12, b1, b2, print_Cases=False, print_Output=False): #ClassifyMin_hom?
+# -------------------------------------------------------------------
+def classifyMin(q1, q2, q3, q12, b1, b2, print_Cases=False, print_Output=False): #ClassifyMin_hom?
if print_Output: print("Run ClassifyMin...")
CaseCount = 0
epsilon = sys.float_info.epsilon #Machine epsilon
diff --git a/src/HelperFunctions.py b/src/HelperFunctions.py
new file mode 100644
index 00000000..426c87a2
--- /dev/null
+++ b/src/HelperFunctions.py
@@ -0,0 +1,176 @@
+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:
+ with codecs.open(QFilePath, 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:
+ with codecs.open(BFilePath, 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
+
+
+
+
+# ----------------------------------------------------------------------------------------------------------------------------
--
GitLab