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 subprocess import Popen, PIPE
+#import sys
+
+
+
+
+
+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)
+
+
+#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)
+
+
+#---------------------------------------------------------------
+
+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
+
+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()
+
+
+# ------------- 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)

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
+
+
+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):
+    # TODO: assert(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)
+
+
+# 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?
+    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/PhaseDiagram.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
+from ClassifyMin import *
+# from CellScript import *
+from mpl_toolkits.mplot3d import Axes3D
+import matplotlib.cm as cm
+from vtk.util import numpy_support
+from pyevtk.hl import gridToVTK
+import time
+# print(sys.executable)
+
+# --------------------------------------------------------------------
+# START :
+# INPUT (Parameters):   alpha, beta, theta, gamma, mu1, rho1
+#
+# -Option 1 : (Case lambda = 0 => q12 = 0)
+#   compute q1,q2,b1,b2 from Formula
+#       Option 1.1 :
+#           set mu_gamma = 'q1' or 'q2' (extreme regimes: gamma \in {0,\infty})
+#       Option 1.2 :
+#           compute mu_gamma with 'Compute_MuGamma' (2D problem much faster then Cell-Problem)
+# -Option 2 :
+#   compute Q_hom & B_eff by running 'Cell-Problem'
+#
+# -> CLASSIFY ...
+#
+# OUTPUT: Minimizer G, angle , type, curvature
+# -----------------------------------------------------------------------
+
+
+# unabhängig von alpha...
+def GetMuGamma(beta,theta,gamma,mu1,rho1, InputFilePath = os.path.dirname(os.getcwd()) +"/inputs/computeMuGamma.parset"):
+    # ------------------------------------ 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)
+        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
+
+
+
+
+
+# --- SETUPS PATHS
+# InputFile  = "/inputs/cellsolver.parset"
+# OutputFile = "/outputs/output.txt"
+
+InputFile  = "/inputs/computeMuGamma.parset"
+OutputFile = "/outputs/outputMuGamma.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)
+
+
+# -------------------------- Input Parameters --------------------
+mu1 = 10.0               # TODO : here must be the same values as in the Parset for computeMuGamma
+rho1 = 1.0
+alpha = 2.0
+beta = 2.0
+theta = 1.0/4.0
+#set gamma either to 1. '0' 2. 'infinity' or 3. a numerical positive value
+gamma = '0'
+# gamma = 'infinity'
+# gamma = 0.5
+
+print('---- Input parameters: -----')
+print('mu1: ', mu1)
+print('rho1: ', rho1)
+print('alpha: ', alpha)
+print('beta: ', beta)
+print('theta: ', theta)
+print('gamma:', gamma)
+print('----------------------------')
+# ----------------------------------------------------------------
+
+# muGamma = GetMuGamma(beta,theta,gamma,mu1,rho1,InputFilePath)
+# # muGamma = GetMuGamma(beta,theta,gamma,mu1,rho1)
+# print('Test MuGamma:', muGamma)
+
+# ------- Options --------
+# print_Cases = True
+# print_Output = True
+
+
+make_3D_plot = True
+make_3D_PhaseDiagram = True
+make_2D_plot = False
+make_2D_PhaseDiagram = False
+# make_3D_plot = False
+# make_3D_PhaseDiagram = False
+# make_2D_plot = True
+# make_2D_PhaseDiagram = True
+
+
+# --- Define effective quantities: q1, q2 , q3 = mu_gamma, q12 ---
+# q1 = harmonicMean(mu1, beta, theta)
+# q2 = arithmeticMean(mu1, beta, theta)
+# --- Set q12
+# q12 = 0.0  # (analytical example)              # TEST / TODO read from Cell-Output
+
+
+
+
+
+# b1 = prestrain_b1(rho1, beta, alpha, theta)
+# b2 = prestrain_b2(rho1, beta, alpha, theta)
+#
+# print('---- Input parameters: -----')
+# print('mu1: ', mu1)
+# print('rho1: ', rho1)
+# 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)
+
+# G, angle, type, kappa = classifyMin(q1, q2, mu_gamma, q12,  b1, b2, print_Cases, print_Output)
+# Test = f(1,2 ,q1,q2,mu_gamma,q12,b1,b2)
+# print("Test", Test)
+
+# ---------------------- MAKE PLOT / Write to VTK------------------------------------------------------------------------------
+
+SamplePoints_3D = 10 # Number of sample points in each direction
+SamplePoints_2D = 10 # Number of sample points in each direction
+SamplePoints_3D = 20 # Number of sample points in each direction
+SamplePoints_2D = 10 # Number of sample points in each direction
+
+
+
+if make_3D_PhaseDiagram:
+    alphas_ = np.linspace(-20, 20, SamplePoints_3D)
+    betas_  = np.linspace(0.01,40.01,SamplePoints_3D)
+    thetas_ = np.linspace(0.01,0.99,SamplePoints_3D)
+    # print('type of alphas', type(alphas_))
+    # print('Test:', type(np.array([mu_gamma])) )
+    alphas, betas, thetas = np.meshgrid(alphas_, betas_, thetas_, indexing='ij')
+    classifyMin_anaVec = np.vectorize(classifyMin_ana)
+
+    # Get MuGamma values ...
+    GetMuGammaVec = np.vectorize(GetMuGamma)
+    muGammas = GetMuGammaVec(betas, thetas, gamma, mu1, rho1)
+    # Classify Minimizers....
+    G, angles, Types, curvature = classifyMin_anaVec(alphas, betas, thetas, muGammas,  mu1, rho1)   # Sets q12 to zero!!!
+    # print('size of G:', G.shape)
+    # print('G:', G)
+    # Out = classifyMin_anaVec(alphas,betas,thetas)
+    # T = Out[2]
+    # --- Write to VTK
+    GammaString = str(gamma)
+
+    VTKOutputName = "outputs/PhaseDiagram3D" + "Gamma" + GammaString
+    gridToVTK(VTKOutputName , alphas, betas, thetas, pointData = {'Type': Types, 'angles': angles, 'curvature': curvature} )
+    print('Written to VTK-File:', VTKOutputName )
+
+if make_2D_PhaseDiagram:
+    alphas_ = np.linspace(-20, 20, SamplePoints_2D)
+    # betas_  = np.linspace(0.01,40.01,1)
+    #fix to one value:
+    betas_ = 2.0;
+    thetas_ = np.linspace(0.01,0.99,SamplePoints_2D)
+    # print('type of alphas', type(alphas_))
+    # print('Test:', type(np.array([mu_gamma])) )
+    alphas, betas, thetas = np.meshgrid(alphas_, betas_, thetas_, indexing='ij')
+    classifyMin_anaVec = np.vectorize(classifyMin_ana)
+
+
+    GetMuGammaVec = np.vectorize(GetMuGamma)
+    muGammas = GetMuGammaVec(betas,thetas,gamma,mu1,rho1)
+    G, angles, Types, curvature = classifyMin_anaVec(alphas,betas,thetas, muGammas,  mu1, rho1)    # Sets q12 to zero!!!
+    # print('size of G:', G.shape)
+    # print('G:', G)
+    print('Types:', Types)
+    # Out = classifyMin_anaVec(alphas,betas,thetas)
+    # T = Out[2]
+    # --- Write to VTK
+    # VTKOutputName = + path + "./PhaseDiagram2DNEW"
+    GammaString = str(gamma)
+
+    VTKOutputName = "outputs/PhaseDiagram2D" + "Gamma_" + GammaString
+    gridToVTK(VTKOutputName , alphas, betas, thetas, pointData = {'Type': Types, 'angles': angles, 'curvature': curvature} )
+    print('Written to VTK-File:', VTKOutputName )
+
+
+# --- Make 3D Scatter plot
+if(make_3D_plot or make_2D_plot):
+    fig = plt.figure()
+    ax = fig.add_subplot(111, projection='3d')
+    colors = cm.plasma(Types)
+    # if make_2D_plot: pnt3d=ax.scatter(alphas,thetas,c=Types.flat)
+    # if make_3D_plot: pnt3d=ax.scatter(alphas,betas,thetas,c=Types.flat)
+    if make_2D_plot: pnt3d=ax.scatter(alphas,thetas,c=angles.flat)
+    if make_3D_plot: pnt3d=ax.scatter(alphas,betas,thetas,c=angles.flat)
+    # cbar=plt.colorbar(pnt3d)
+    # cbar.set_label("Values (units)")
+    ax.set_xlabel('alpha')
+    ax.set_ylabel('beta')
+    if make_3D_plot: ax.set_zlabel('theta')
+    plt.show()
+    # plt.savefig('common_labels.png', dpi=300)
+    # print('T:', T)
+    # print('Type 1 occured here:', np.where(T == 1))
+    # print('Type 2 occured here:', np.where(T == 2))
+
+
+print(alphas_)
+print(betas_)
+# ALTERNATIVE
+# colors = ("red", "green", "blue")
+# groups = ("Type 1", "Type2", "Type3")
+#
+# # Create plot
+# fig = plt.figure()
+# ax = fig.add_subplot(1, 1, 1)
+#
+# for data, color, group in zip(Types, colors, groups):
+#     # x, y = data
+#     ax.scatter(alphas, thetas, alpha=0.8, c=color, edgecolors='none', label=group)
+#
+# plt.title('Matplot scatter plot')
+# plt.legend(loc=2)
+# plt.show()