src/CellScript.py

+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 HelperFunctions import *
+# from scipy.io import loadmat #not Needed anymore?
+import codecs
+import sys
+
+
+# ----------------------------------------------------------------------------------------------------------------------------
+
+
+
+
+
+
+# ----- 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)
+
+
+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)
+
+print('Read effective quantities...')
+Q, B = ReadEffectiveQuantities()
+print('Q:', Q)
+print('B:', B)
+
+
+
+# Compare symbolicMinimization with Classification 'ClassifyMin' :
+# print('Compare_Classification...')
+# Compare_Classification(alpha,beta,theta,gamma,mu1,rho1,InputFilePath)
+
+
+
+# ------------- 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)

src/ClassifyMin.py

+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 sys
+# print(sys.executable)
+# from subprocess import Popen, PIPE
+
+
+# --------------------------------------------------
+# 'classifyMin' classifies Minimizers by utilizing the result of
+# Lemma1.6
+#
+#
+#
+#
+# 'classifyMin_ana': (Special-Case : Lemma1.4)
+# ..additionally assumes Poisson-ratio=0 => q12==0
+#
+#
+#
+# Output : MinimizingMatrix, Angle, Type, Curvature
+
+
+
+
+
+def harmonicMean(mu_1, beta, theta):
+    return mu_1*(beta/(theta+(1-theta)*beta))
+
+
+def arithmeticMean(mu_1, beta, theta):
+    return mu_1*((1-theta)+theta*beta)
+
+
+def prestrain_b1(rho_1, beta, alpha, theta):
+    return (3.0*rho_1/2.0)*beta*(1-(theta*(1+alpha)))
+
+
+def prestrain_b2(rho_1, beta, alpha, theta):
+    return (3.0*rho_1/(4.0*((1.0-theta) + theta*beta)))*(1-theta*(1+beta*alpha))
+
+
+# Define function to be minimized
+def f(a1, a2, q1, q2, q3, q12, b1, b2):
+    A = np.array([[q1, q3 + q12/2.0], [q3 + q12/2.0, q2]])
+    B = np.array([-2.0*q1*b1-q12*b2, -2.0*q2*b2-q12*b1])
+    a = np.array([a1, a2])
+    tmp = np.dot(A, a)
+    tmp2 = np.dot(a, tmp)
+    tmpB = np.dot(B, a)
+    return tmp2 + tmpB + q1*(b1**2) + q2*(b2**2) + 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):
+
+    # 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)
+    # print('q2: ', q2)
+    b1 = prestrain_b1(rho_1, beta, alpha,theta)
+    b2 = prestrain_b2(rho_1, beta, alpha,theta)
+    return classifyMin(q1, q2, q3, q12,  b1, b2,  print_Cases, print_Output)
+
+
+
+# Matrix Version that just gets matrices Q & B
+def classifyMin_mat(Q,B,print_Cases=False, print_Output=False):
+    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]
+    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 )
+# where
+# R1 : unique local (global) minimizer which is not axial
+# R2 : continuum of local (global) minimizers which are not axial
+# R3 : one or two local (global) minimizers which are axial
+# Partition given by
+# 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?
+    # Assumption of Classification-Lemma1.6:
+    #  1. [b3 == 0]
+    #  2. Q is orthotropic i.e. q13 = q31 = q23 = q32 == 0
+
+    # TODO: check if Q is orthotropic here - assert()
+
+
+    if print_Output: print("Run ClassifyMin...")
+    CaseCount = 0
+    epsilon = sys.float_info.epsilon #Machine epsilon
+
+    B = np.array([-2.0*q1*b1-q12*b2, -2.0*q2*b2-q12*b1])
+    A = np.array([[q1, q3 + q12/2.0], [q3 + q12/2.0, q2]])
+    # print('Matrix A:',  A)
+    # print('Matrix B:', B)
+    # print('Matrix rank of A:', np.linalg.matrix_rank(A))
+    # print('shape of A:', A.shape)
+    # print('shape of B:', B.shape)
+    # print('Matrix [A,B]:', np.c_[A, B])
+    # print('Matrix rank of [A,B]:', np.linalg.matrix_rank(np.c_[A, B]))
+    # print('shape of [A,B]:', C.shape)
+    #
+    # x = np.linalg.solve(A, B)     # works only if A is not singular!!
+    # print("Solution LGS:", x)
+    # # print("sym Matrix", sym.Matrix(([A],[B])))
+    #
+    # # Test
+    # C = np.array([[1, 0], [0, 0]])
+    # d = np.array([5, 0])
+    # y = np.linalg.lstsq(C, d)[0]
+    # print("Solution LGS:", y)
+    # T = np.c_[C, d]
+    # print('T:', T)
+    # Trref = sym.Matrix(T).rref()[0]
+    # Out = np.array(Trref, dtype=float)
+    # print('rref:', Out)
+
+    determinant = q1*q2 - (q3**2 + 2*q3*q12 + q12**2)
+    if print_Cases: print("determinant:", determinant)
+
+    # Define values for axial-Solutions (b1*,0) & (0,b2*)
+    b1_star = (2.0*q1*b1 + b2*q12)/(2*q1)
+    b2_star = (2.0*q2*b2 + b1*q12)/(2*q2)
+
+    # ------------------------------------ Parabolic Case -----------------------------------
+    if abs(determinant) < epsilon:
+        if print_Cases: print('P : parabolic case (determinant equal zero)')
+        # if print_Cases: print('P : parabolic case (determinant equal zero)')
+
+        # check if B is in range of A
+        # check if rank(A) == rank([A,B])
+        # OK this way? (or use Sympy?)
+        if np.linalg.matrix_rank(A) == np.linalg.matrix_rank(np.c_[A, B]):
+            if print_Cases: print('P1 (B is in the range of A)')
+
+            if (q12+q3)/2.0 <= -1.0*epsilon:
+                print('should not happen(not compatible with det = 0)')
+
+            if (abs(B[0]) < epsilon and abs(B[1]) < epsilon) and (q12+q3)/2.0 >= epsilon:
+                if print_Cases: print('P1.1')
+                a1 = 0.0
+                a2 = 0.0
+                type = 3
+                CaseCount += 1
+
+            if (abs(B[0]) >= epsilon or abs(B[1]) >= epsilon) and (q12+q3)/2.0 >= epsilon:
+                # Continuum of minimizers located along a line of negative slope in Lambda
+                if print_Cases: print('P1.2 (Continuum of minimizers located along a line of negative slope in Lambda) ')
+
+                # Just solve Aa* = b (alternatively using SymPy ?)
+                # we know that A is singular (det A = 0 ) here..
+                # & we know that there are either infinitely many solutions or a unique solution ...
+                # ---- determine one via Least Squares
+                # "If A is square and of full rank, then x (but for round-off error)
+                # is the “exact” solution of the equation. Else, x minimizes the
+                # Euclidean 2-norm || b-Ax ||. If there are multiple minimizing solutions,
+                # the one with the smallest 2-norm  is returned. ""
+                a = np.linalg.lstsq(A, B)[0]                            # TODO check is this Ok ?
+                print("Solution LGS: a =", a)
+                a1 = a[0]
+                a2 = a[1]
+
+                type = 2
+                CaseCount += 1
+        else:
+            if print_Cases: print('P2 (B is not in the range of A)')
+            # local Minimizers occur on the boundary of Lambda...
+            # There are at most two local minima and they are either
+            # (b1_star, 0) or (0, b2_star)
+            # could also outsource this to another function..
+            if f(b1_star, 0, q1, q2, q3, q12, b1, b2) < f(0, b2_star, q1, q2, q3, q12, b1, b2):
+                a1 = b1_star
+                a2 = 0.0
+                type = 3  # 1
+                CaseCount += 1
+            if f(b1_star, 0, q1, q2, q3, q12, b1, b2) > f(0, b2_star, q1, q2, q3, q12, b1, b2):
+                a1 = 0
+                a2 = b2_star
+                type = 3  # 2
+                CaseCount += 1
+            # TODO Problem: angle depends on how you choose... THE angle is not defined for this case
+            if f(b1_star, 0, q1, q2, q3, q12, b1, b2) == f(0, b2_star, q1, q2, q3, q12, b1, b2):
+                # Two Minimizers pick one
+                a1 = b1_star
+                a2 = 0.0
+                type = 3  # 4
+                CaseCount += 1
+
+    # ------------------------------------ Elliptic Case -----------------------------------
+    if determinant >= epsilon:
+        if print_Cases: print('E : elliptic case (determinant greater zero)')
+        a1_star = -(2*(b1*(q12**2) + 2*b1*q3*q12 - 4*b1*q1*q2 + 4*b2*q2*q3)) / \
+            (4*(q3**2) + 4*q3*q12 + (q12**2) - 4*q1*q2)
+        a2_star = -(2*(b2*(q12**2) + 2*b2*q3*q12 + 4*b1*q1*q3 - 4*b2*q1*q2)) / \
+            (4*(q3**2) + 4*q3*q12 + (q12**2) - 4*q1*q2)
+        prod = a1_star*a2_star
+
+        if prod >= epsilon:
+            if print_Cases: print('(E1) - inside Lambda ')
+            a1 = a1_star
+            a2 = a2_star
+            type = 1  # non-axial Minimizer
+            CaseCount += 1
+        if abs(prod) < epsilon:   # same as prod = 0 ? or better use <=epsilon ?
+            if print_Cases: print('(E2) - on the boundary of Lambda ')
+            a1 = a1_star
+            a2 = a2_star
+            type = 3  # could check which axis: if a1_star or a2_star close to zero.. ?
+            CaseCount += 1
+
+        if prod <= -1.0*epsilon:
+            if print_Cases: print('(E3) - Outside Lambda ')
+            if f(b1_star, 0, q1, q2, q3, q12, b1, b2) < f(0, b2_star, q1, q2, q3, q12, b1, b2):
+                a1 = b1_star
+                a2 = 0.0
+                type = 3  # 1
+                CaseCount += 1
+            if f(b1_star, 0, q1, q2, q3, q12, b1, b2) > f(0, b2_star, q1, q2, q3, q12, b1, b2):
+                a1 = 0
+                a2 = b2_star
+                type = 3  # 2
+                CaseCount += 1
+
+            # TODO Problem: angle depends on how you choose... THE angle is not defined for this case
+            if f(b1_star, 0, q1, q2, q3, q12, b1, b2) == f(0, b2_star, q1, q2, q3, q12, b1, b2):
+                # Two Minimizers pick one
+                a1 = b1_star
+                a2 = 0.0
+                type = 3  # 4
+                CaseCount += 1
+
+    # ------------------------------------ Hyperbolic Case -----------------------------------
+    if determinant <= -1.0*epsilon:
+        if print_Cases: print('H : hyperbolic case (determinant smaller zero)')
+        # One or two minimizers wich are axial
+        type = 3 # (always type 3)
+        if f(b1_star, 0, q1, q2, q3, q12, b1, b2) < f(0, b2_star, q1, q2, q3, q12, b1, b2):
+            a1 = b1_star
+            a2 = 0.0
+            # type = 3  # 1
+            CaseCount += 1
+        if f(b1_star, 0, q1, q2, q3, q12, b1, b2) > f(0, b2_star, q1, q2, q3, q12, b1, b2):
+            a1 = 0
+            a2 = b2_star
+            # type = 3  # 2
+            CaseCount += 1
+        # TODO can add this case to first or second ..
+        if f(b1_star, 0, q1, q2, q3, q12, b1, b2) == f(0, b2_star, q1, q2, q3, q12, b1, b2):
+            # Two Minimizers pick one
+            a1 = b1_star
+            a2 = 0.0
+            # type = 3  # 4
+            CaseCount += 1
+    # ---------------------------------------------------------------------------------------
+
+    if (CaseCount > 1):
+        print('Error: More than one Case happened!')
+
+    # compute a3
+    # a3 = math.sqrt(2.0*a1*a2)   # never needed?
+
+    # compute the angle <(e,e_1) where Minimizer = kappa* (e (x) e)
+    e = [math.sqrt((a1/(a1+a2))), math.sqrt((a2/(a1+a2)))]
+    angle = math.atan2(e[1], e[0])
+
+    # compute kappa
+    kappa = (a1 + a2)
+
+    # Minimizer G
+    Minimizer = np.array([[a1, math.sqrt(a1*a2)], [math.sqrt(a1*a2), a2]],dtype=object)
+    # Minimizer = np.array([[a1, math.sqrt(a1*a2)], [math.sqrt(a1*a2), a2]])
+
+    if print_Output:
+        print('--- Output ClassifyMin ---')
+        print("Minimizing Matrix G:")
+        print(Minimizer)
+        print("angle = ", angle)
+        print("type: ", type)
+        print("kappa = ", kappa)
+
+    return Minimizer, angle, type, kappa
+# ------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
+
+
+# ---------------------------------------------- Main ---------------------
+
+# --- Input Parameters ----
+# mu_1 = 1.0
+# rho_1 = 1.0
+# alpha = 9.0
+# beta = 2.0
+# theta = 1.0/8.0
+
+# # define q1, q2 , mu_gamma, q12
+# # 1. read from Cell-Output
+# # 2. define values from analytic formulas (expect for mu_gamma)
+# q1 = (1.0/6.0)*harmonicMean(mu_1, beta, theta)
+# q2 = (1.0/6.0)*arithmeticMean(mu_1, beta, theta)
+# # TEST
+# q12 = 0.0  # (analytical example)
+# # q12 = 12.0  # (analytical example)
+
+# set mu_gamma  to value or read from Cell-Output
+# mu_gamma = q1   # TODO read from Cell-Output
+# b1 = prestrain_b1(rho_1, beta, alpha, theta)
+# b2 = prestrain_b2(rho_1, beta, alpha, theta)
+
+
+# print('---- Input parameters: -----')
+# print('mu_1: ', mu_1)
+# print('rho_1: ', rho_1)
+# print('alpha: ', alpha)
+# print('beta: ', beta)
+# print('theta: ', theta)
+# print("q1: ", q1)
+# print("q2: ", q2)
+# print("mu_gamma: ", mu_gamma)
+# print("q12: ", q12)
+# print("b1: ", b1)
+# print("b2: ", b2)
+# print('----------------------------')
+# # print("machine epsilon", sys.float_info.epsilon)
+#
+#
+# # ------- Options --------
+# print_Cases = True
+# print_Output = True
+
+# G, angle, type, kappa = classifyMin(q1, q2, mu_gamma, q12,  b1, b2, print_Cases, print_Output)
+#
+# G, angle, type, kappa = classifyMin_ana(alpha, beta, theta, mu_gamma, q12, print_Cases, print_Output)
+#
+# Out = classifyMin_ana(alpha, beta, theta, mu_gamma, q12, print_Cases, print_Output)
+
+# print('TEST:')
+# Out = classifyMin_ana(alpha, beta, theta)
+
+# print('Out[0]', Out[0])
+# print('Out[1]', Out[1])
+# print('Out[2]', Out[2])
+# print('Out[3]', Out[3])
+
+
+# #supress certain Outout..
+# _,_,T,_ = classifyMin_ana(alpha, beta, theta, mu_gamma, q12, print_Cases, print_Output)
+# print('Output only type..:', T)
+
+# Test = f(1,2 ,q1,q2,mu_gamma,q12,b1,b2)
+# print("Test", Test)
+
+# -----------------------------------------------------------------------------------------------------------------------------------------------------------------

src/HelperFunctions.py

+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 SetParametersCellProblem(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)^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)^gamma=.*','gamma='+str(gamma),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()
+
+#TODO combine these...
+
+def SetParametersComputeMuGamma(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)^beta=.*','beta='+str(beta),filedata)
+    filedata = re.sub('(?m)^theta=.*','theta='+str(theta),filedata)
+    filedata = re.sub('(?m)^gamma=.*','gamma='+str(gamma),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()
+
+
+
+
+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()
+        SetParametersCellProblem(alpha,beta,theta,gamma,mu1,rho1, InputFilePath)
+        # --- 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 GetCellOutput(alpha,beta,theta,gamma,mu1,rho1, InputFilePath = os.path.dirname(os.getcwd()) +"/inputs/computeMuGamma.parset",
+                OutputFilePath = os.path.dirname(os.getcwd()) + "/outputs/outputMuGamma.txt" ):
+        RunCellProblem(alpha,beta,theta,gamma,mu1,rho1, InputFilePath)
+        print('Read effective quantities...')
+        Q, B = ReadEffectiveQuantities()
+        # print('Q:', Q)
+        # print('B:', B)
+        return Q, B
+
+
+# unabhängig von alpha...
+def GetMuGamma(beta,theta,gamma,mu1,rho1, InputFilePath = os.path.dirname(os.getcwd()) +"/inputs/computeMuGamma.parset",
+                OutputFilePath = os.path.dirname(os.getcwd()) + "/outputs/outputMuGamma.txt" ):
+    # ------------------------------------ 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)
+
+        SetParametersComputeMuGamma(beta,theta,gamma,mu1,rho1, InputFilePath)
+        # 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'
+        # 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(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
+        # 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')
+        else:
+            print('Comparison unsuccessful')
+
+        return comparison_successful
+
+
+
+
+# ----------------------------------------------------------------------------------------------------------------------------