diff --git a/src/ClassifyMin.py b/src/ClassifyMin.py
index 87e9f935fb8773b2261ece1826eb6c727f4b6ce4..a0b090961cd16a98904890d10234fb5bd83f3bc8 100644
--- a/src/ClassifyMin.py
+++ b/src/ClassifyMin.py
@@ -8,192 +8,366 @@ import fileinput
import re
import matlab.engine
import sys
+# print(sys.executable)
# from subprocess import Popen, PIPE
-
-def harmonicMean(mu_1,beta,theta):
+def harmonicMean(mu_1, beta, theta):
return mu_1*(beta/(theta+(1-theta)*beta))
-def arithmeticMean(mu_1,beta,theta):
+def arithmeticMean(mu_1, beta, theta):
return mu_1*((1-theta)+theta*beta)
-def prestrain_b1(rho_1,beta,alpha):
- return (3.0*rho_1/2.0)*beta*(1-(theta*(1+alpha)));
+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):
- return (3.0*rho_1/(4.0*((1.0-theta)+ theta*beta)))*(1-theta*(1+beta*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):
+# 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]])
- print(A)
B = np.array([-2.0*q1*b1-q12*b2, -2.0*q2*b2-q12*b1])
- print(B)
a = np.array([a1, a2])
- print(a)
- tmp = np.dot(A,a)
- print(tmp)
- tmp2 = np.dot(a,tmp)
- print(tmp2)
- tmpB = np.dot(B,a)
- print(tmpB)
+ # print(A)
+ # print(B)
+ # print(a)
+ tmp = np.dot(A, a)
+ # print(tmp)
+ tmp2 = np.dot(a, tmp)
+ # print(tmp2)
+ tmpB = np.dot(B, a)
+ # print(tmpB)
# print(q1*(b1**2))
# print(q2*(b2**2))
# print(q12*b1*b2)
return tmp2 + tmpB + q1*(b1**2) + q2*(b2**2) + q12*b1*b2
-# Classify Type of minimizer 1 = (I) , 2 = (II) , 3 = (III) , 4 = (IV)
-def classifyMin(q1, q2, q3, q12, b1, b2, print_output = "True"):
- print("Run ClassifyMin...")
+ # ---- Alternative Version using alpha,beta,theta ,mu_1,rho_1
+
+def classifyMin_geo(alpha,beta,theta,q3,q12,mu_1,rho_1,print_Cases=False, print_Output=False):
+ q1 = harmonicMean(mu_1, beta, theta)
+ q2 = arithmeticMean(mu_1, beta, theta)
+ 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)
+
+# TEST
+# def classifyMin_geo(alpha,beta,theta,q3,q12,mu_1, rho_1, print_Cases=False, print_Output=False):
+# mu_1 = 1.0
+# rho_1 = 1.0
+# q12 = 0.0
+# q1 = harmonicMean(mu_1, beta, theta)
+# q2 = arithmeticMean(mu_1, beta, theta)
+# b1 = prestrain_b1(rho_1, beta, alpha,theta)
+# b2 = prestrain_b2(rho_1, beta, alpha,theta)
+# q3 = q1
+#
+#
+# return classifyMin(q1, q2, q3, q12, b1, b2)
+
+
+# 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
+ epsilon = sys.float_info.epsilon #Machine epsilon
- #TEST
- # A = np.array([[1, 2], [3, 4]])
- # kappa = 1.0
- # angle = 1.57
- # type = 2
+ 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)
- print("determinant:", determinant)
+ if print_Cases: print("determinant:", determinant) # TODO ..add everywhere if print_Cases:
+ # 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:
- print('parabolic case (determinant equal zero)')
+ 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 (TODO)
+ # 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) ')
+
+ # TODO - what to set for a1, a2 ?
+
+ # 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?...
+ 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:
- print('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)
+ 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
- b1_star = (2.0*q1*b1 + b2*q12)/(2*q1)
- b2_star = (2.0*q2*b2 + b1*q12)/(2*q2)
-
if prod >= epsilon:
- print('(E1) - inside Lambda ')
+ 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 # non-axial Minimizer
+ type = 3 # could check which axis if a1_star or a2_star close to zero.. ?
CaseCount += 1
- if abs(prod) < epsilon:
- print('(E2) - on the boundary of Lambda ')
+
# if q2*b2**2 < q1*b1**2: # Needs to be updated to Mixed Term!!! just define f as a function and check value?!
- #check
- if f(b1_star,0,q1,q2,mu_gamma,q12,b1,b2) < f(0, b2_star, q1,q2,mu_gamma,q12,b1,b2):
+ # check
+ # 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 = 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 = 2
+ # CaseCount += 1
+ # 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 = 4
+ # 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 = 1
+ type = 3 # 1
CaseCount += 1
- if f(b1_star,0,q1,q2,mu_gamma,q12,b1,b2) > f(0, b2_star, q1,q2,mu_gamma,q12,b1,b2):
+ 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 = 2
+ type = 3 # 2
CaseCount += 1
- if f(b1_star,0,q1,q2,mu_gamma,q12,b1,b2) = f(0, b2_star, q1,q2,mu_gamma,q12,b1,b2):
+ # TODO ...does type4 happen here?
+ # TODO Problem? angle depends on how you choose?...
+ 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 = 4
- CaseCount += 1
-
- if prod < -1.0*epsilon:
- print('(E3) - Outside Lambda ')
- if f(b1_star,0,q1,q2,mu_gamma,q12,b1,b2) < f(0, b2_star, q1,q2,mu_gamma,q12,b1,b2):
- a1 = b1_star
- a2 = 0.0
- type = 1
- CaseCount += 1
- if f(b1_star,0,q1,q2,mu_gamma,q12,b1,b2) > f(0, b2_star, q1,q2,mu_gamma,q12,b1,b2):
- a1 = 0
- a2 = b2_star
- type = 2
+ type = 3 # 4
CaseCount += 1
- #TODO ...does type4 happen here?
+ # ------------------------------------ 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
+ 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) # ever needed?
- return A, angle ,type ,kappa
+ # 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)
+ #Test
+ 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)
-# ---------------------------------------------- Main ---------------------
-print('Running Python Code')
+ return Minimizer, angle, type, kappa
+# ------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
-# --- Input Parameters ----
-
-mu_1 = 1.0
-rho_1 = 1.0
-alpha = 14.0
-beta = 20.0
-theta = 1.0/4.0
+# ---------------------------------------------- Main ---------------------
-# define q1, q2 , mu_gamma, q12
-# 1. read from Cell-Output
-# 2. define values from analytic formulas (expect for mu_gamma)
-q1 = harmonicMean(mu_1,beta,theta)
-q2 = arithmeticMean(mu_1,beta,theta)
-q12 = 0.0 # (analytical example)
+# --- 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 = harmonicMean(mu_1, beta, theta)
+# q2 = 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)
-b2 = prestrain_b2(rho_1,beta,alpha)
-
-
-print("q1: ", q1)
-print("q2: ", q2)
-print("mu_gamma : ", mu_gamma)
-print("b1: ", b1)
-print("b2: ", b2)
-
-print("machine epsilon", sys.float_info.epsilon)
-
-## ------- options --------
-# print_output = false;
-
-
-
-
-
-
-
-A, angle, type, kappa = classifyMin(q1, q2, mu_gamma, q12, b1, b2)
-
-print("Matrix A")
-print(A)
-print("angle",angle)
-print("type",type)
-print("kappa", kappa)
-
-
-Test = f(1,2 ,q1,q2,mu_gamma,q12,b1,b2)
-print("Test", Test)
-
-
-
-# Run several times / make Array :
-# Gamma_Values = np.linspace(0.01, 2.5, num=12)
-# print(Gamma_Values)
-
-
-# Make Plot
-# plt.figure()
-# plt.title(r'$\mu_\gamma(\gamma)$-Plot') # USE MATHEMATICAL EXPRESSIONS IN TITLE
-# plt.plot(Gamma_Values, mu_gamma)
-# plt.scatter(Gamma_Values, mu_gamma)
-# plt.xlabel("$\gamma$")
-# plt.ylabel("$\mu_\gamma$")
-# plt.legend()
-# plt.show()
+# 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_geo(alpha, beta, theta, mu_gamma, q12, print_Cases, print_Output)
+#
+# Out = classifyMin_geo(alpha, beta, theta, mu_gamma, q12, print_Cases, print_Output)
+
+# print('TEST:')
+# Out = classifyMin_geo(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_geo(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)
+
+# -----------------------------------------------------------------------------------------------------------------------------------------------------------------
diff --git a/src/PhaseDiagram.py b/src/PhaseDiagram.py
new file mode 100644
index 0000000000000000000000000000000000000000..9d8de8d40a546e43ae7dfe5ebe54f1e1021a6ce0
--- /dev/null
+++ b/src/PhaseDiagram.py
@@ -0,0 +1,141 @@
+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 mpl_toolkits.mplot3d import Axes3D
+import matplotlib.cm as cm
+from vtk.util import numpy_support
+from pyevtk.hl import gridToVTK
+# print(sys.executable)
+# from subprocess import Popen, PIPE
+
+# Not Needed:
+# def ndm(*args):
+# return [x[(None,)*i+(slice(None),)+(None,)*(len(args)-i-1)] for i, x in enumerate(args)]
+
+# ---------------------------------------------- Main ---------------------s
+# --- Input Parameters ----
+mu_1 = 1.0
+rho_1 = 1.0
+alpha = 9.0
+beta = 2.0
+theta = 1.0/8.0
+# print('Type of theta', type(theta))
+
+# ------- Options --------
+# print_Cases = True
+# print_Output = True
+
+make_3D_plot = False
+make_2D_plot = False
+make_3D_PhaseDiagram = True
+make_2D_PhaseDiagram = False
+
+
+# --- Define q1, q2 , mu_gamma, q12 ---
+# 1. read from Cell-Output
+# 2. define values from analytic formulas (expect for mu_gamma)
+q1 = harmonicMean(mu_1, beta, theta)
+q2 = arithmeticMean(mu_1, beta, theta)
+
+# TEST
+q12 = 0.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)
+
+
+# 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------------------------------------------------------------------------------
+
+if make_3D_PhaseDiagram:
+ alphas_ = np.linspace(-20, 20, 80)
+ # betas_ = np.linspace(0.01,40.01,1)
+ #fix to one value:
+ betas_ = 2.0;
+ thetas_ = np.linspace(0.01,0.99,80)
+ # print('type of alphas', type(alphas_))
+ # print('Test:', type(np.array([mu_gamma])) )
+ alphas, betas, thetas = np.meshgrid(alphas_, betas_, thetas_, indexing='ij')
+ classifyMin_geoVec = np.vectorize(classifyMin_geo)
+ G, angles, Types, curvature = classifyMin_geoVec(alphas,betas,thetas,mu_gamma,q12, mu_1, rho_1)
+ # print('size of G:', G.shape)
+ # print('G:', G)
+ # Out = classifyMin_geoVec(alphas,betas,thetas)
+ # T = Out[2]
+ # --- Write to VTK
+ VTKOutputName = "./PhaseDiagram3D"
+ 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, 80)
+ # betas_ = np.linspace(0.01,40.01,1)
+ #fix to one value:
+ betas_ = 2.0;
+ beta = np.array([betas_])
+ thetas_ = np.linspace(0.01,0.99,80)
+ # print('type of alphas', type(alphas_))
+ # print('Test:', type(np.array([mu_gamma])) )
+ alphas, thetas = np.meshgrid(alphas_, thetas_, indexing='ij')
+ classifyMin_geoVec = np.vectorize(classifyMin_geo)
+ G, angles, Types, curvature = classifyMin_geoVec(alphas,beta,thetas,mu_gamma,q12, mu_1, rho_1)
+ # print('size of G:', G.shape)
+ # print('G:', G)
+ # Out = classifyMin_geoVec(alphas,betas,thetas)
+ # T = Out[2]
+ # --- Write to VTK
+ # VTKOutputName = "./PhaseDiagram2D"
+ # gridToVTK(VTKOutputName , alphas, thetas, beta, 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')
+ # pnt3d=ax.scatter(x,y,z,c=z.flat)
+ colors = cm.plasma(Types)
+ # pnt3d=ax.scatter(x,y,z,c=T.flat)
+ 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)
+ # pnt3d=ax.scatter(x,y,z,c=colors.flat)
+ cbar=plt.colorbar(pnt3d)
+ cbar.set_label("Values (units)")
+ ax.set_xlabel('alpha')
+ ax.set_ylabel('beta')
+ ax.set_zlabel('theta')
+ plt.show()
+ # plt.savefig('common_labels.png', dpi=300)
+
+ # print('TEST 3D Scatter')
+ # print('T:', T)
+ # print('Type 1 occured here:', np.where(T == 1))
+ # print('Type 2 occured here:', np.where(T == 2))