diff --git a/Matlab-Programs/symMinimization.m b/Matlab-Programs/symMinimization.m
index 74ad4aca2241b62bb84fc7eeaa82983016eb8617..f714a3203b0aa45b5f71fee8b565afcb66cd56ec 100644
--- a/Matlab-Programs/symMinimization.m
+++ b/Matlab-Programs/symMinimization.m
@@ -28,6 +28,11 @@ f_minus(v1,v2,q1,q2,q3,q12,b1,b2,b3) = q1*v1^4 + q2*v2^4+2*q3*v1^2*v2^2+2*(q1*b1
% ---- Fix parameters
+
+% Epsilon used:
+epsilon = 1e-8;
+
+
if ~exist('InputPath','var')
% third parameter does not exist, so default it to something
absPath = "/home/klaus/Desktop/DUNE/dune-microstructure/outputs";
@@ -52,7 +57,7 @@ Bmat = full(Bmat);
-% --- TODO CHECK: assert if Q is orthotropic ???
+% --- TODO CHECK: assert if Q is orthotropic ??? check ifq13=q31=q23=q32= 0 ?
if print_Input
@@ -66,15 +71,16 @@ if print_Input
% Test: issymmetric(Qmat) does not work for float matrices?
% symmetric part 0.5*(Qmat+Qmat')
% anti-symmetric part 0.5*(Qmat-Qmat')
- if norm(0.5*(Qmat-Qmat'),'fro') < 1e-10
+ if norm(0.5*(Qmat-Qmat'),'fro') < 1e-8
fprintf('Qmat (close to) symmetric \n')
+ norm(0.5*(Qmat-Qmat'),'fro') % TEST
else
fprintf('Qmat not symmetric \n')
end
% Check if B_eff is diagonal this is equivalent to b3 == 0
- if abs(Bmat(3)) < 1e-10
+ if abs(Bmat(3)) < 1e-8
fprintf('B_eff is diagonal (b3 == 0) \n')
else
fprintf('B_eff is NOT diagonal (b3 != 0) \n')
@@ -231,11 +237,10 @@ end
% SP_value = T_minus(idx) % not needed?
Matrix = sign*(SP'*SP);
- if norm(double(Matrix-Minimizer),'fro') < eps %check is this sufficient here?
- % both StationaryPoints correspond to the same
- % (Matrix-)Minimizer
+ if norm(double(Matrix-Minimizer),'fro') < 1e-8 %check is this sufficient here?
+% fprintf('both StationaryPoints correspond to the same(Matrix-)Minimizer')
else
- % StationaryPoint corresponds to a different (Matrix-)Minimizer
+% fprintf('StationaryPoint corresponds to a different (Matrix-)Minimizer')
MinimizerCount = MinimizerCount + 1;
end
end
@@ -289,7 +294,7 @@ G = double(Minimizer);
if MinimizerCount == 1
% fprintf('Unique Minimzier')
% Check if Minimizer is axial or non-axial:
- if (abs(angle-pi/2) < 1e-9 || abs(angle) < 1e-9) % axial Minimizer
+ if (abs(angle-pi/2) < 1e-8 || abs(angle) < 1e-8) % axial Minimizer
Type = 3;
else % non-axial Minimizer
Type = 1;
@@ -297,7 +302,7 @@ if MinimizerCount == 1
elseif MinimizerCount == 2
% fprintf('Two Minimziers')
% Check if Minimizer is axial or non-axial:
- if (abs(angle-pi/2) < 1e-9 || abs(angle) < 1e-9) % axial Minimizer
+ if (abs(angle-pi/2) < 1e-8 || abs(angle) < 1e-8) % axial Minimizer
Type = 3;
else % non-axial Minimizer
fprintf('ERROR: Two non-axial Minimizers cannot happen!')
@@ -305,7 +310,7 @@ elseif MinimizerCount == 2
else
% fprintf('1-Parameter family of Minimziers')
% Check if Minimizer is axial or non-axial:
- if (abs(angle-pi/2) < 1e-9 || abs(angle) < 1e-9) % axial Minimizer
+ if (abs(angle-pi/2) < 1e-8 || abs(angle) < 1e-8) % axial Minimizer
% fprintf('ERROR: axial Minimizers cannot happen for 1-Parameter Family!')
else % non-axial Minimizer
Type = 2;
diff --git a/src/CellScript.py b/src/CellScript.py
index 19cec12ef2c51c49374e10a2926e755cbe31c6f4..b7afd08b49733e1ce4035a65eac2d74a07088711 100644
--- a/src/CellScript.py
+++ b/src/CellScript.py
@@ -9,168 +9,12 @@ import re
import matlab.engine
import time
from ClassifyMin import *
+from HelperFunctions 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 -----
@@ -188,97 +32,6 @@ 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
@@ -298,23 +51,17 @@ print('----------------------------')
-
+#
# print('RunCellProblem...')
# RunCellProblem(alpha,beta,theta,gamma,mu1,rho1,InputFilePath)
-
-# TEST read Matrix file
-# Test = loadmat(path + '/outputs/QMatrix.mat')
-
-
-
-
+# TEST:
print('Compare_Classification...')
Compare_Classification(alpha,beta,theta,gamma,mu1,rho1,InputFilePath)
-# ------------- RUN Matlab symbolic minimization program
+# ------------- RUN Matlab symbolic minimization program 'symMinimization'
eng = matlab.engine.start_matlab()
# s = eng.genpath(path + '/Matlab-Programs')
s = eng.genpath(path)
@@ -323,11 +70,12 @@ eng.addpath(s, nargout=0)
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...')
-G, angle, type, kappa = eng.symMinimization(Inp,Inp,Inp,Inp, nargout=4) #Name of program:symMinimization
+#Arguments: symMinization(print_Input,print_statPoint,print_Output,make_FunctionPlot, InputPath)
+G, angle, type, kappa = eng.symMinimization(Inp_T,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)