From daea3744e179d67323bdc543df031c5f92e3fd6d Mon Sep 17 00:00:00 2001
From: Klaus <klaus.boehnlein@tu-dresden.de>
Date: Mon, 13 Sep 2021 20:11:54 +0200
Subject: [PATCH] Update PhaseDiagram.py
---
src/ClassifyMin.py | 98 +++++++++++++--------------------------------
src/PhaseDiagram.py | 66 ++++++++++++++++--------------
2 files changed, 64 insertions(+), 100 deletions(-)
diff --git a/src/ClassifyMin.py b/src/ClassifyMin.py
index a0b09096..0c84fdab 100644
--- a/src/ClassifyMin.py
+++ b/src/ClassifyMin.py
@@ -33,44 +33,24 @@ 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])
- # 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
+# ---- 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!)?
- # ---- 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)
+ 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)
-# 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
@@ -114,7 +94,7 @@ def classifyMin(q1, q2, q3, q12, b1, b2, print_Cases=False, print_Output=False
# print('rref:', Out)
determinant = q1*q2 - (q3**2 + 2*q3*q12 + q12**2)
- if print_Cases: print("determinant:", determinant) # TODO ..add everywhere if print_Cases:
+ 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)
@@ -125,9 +105,8 @@ def classifyMin(q1, q2, q3, q12, b1, b2, print_Cases=False, print_Output=False
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 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)')
@@ -146,8 +125,6 @@ def classifyMin(q1, q2, q3, q12, b1, b2, print_Cases=False, print_Output=False
# 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 ...
@@ -156,7 +133,7 @@ def classifyMin(q1, q2, q3, q12, b1, b2, print_Cases=False, print_Output=False
# 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 ?
+ a = np.linalg.lstsq(A, B)[0] # TODO check is this Ok ?
print("Solution LGS: a =", a)
a1 = a[0]
a2 = a[1]
@@ -179,7 +156,7 @@ def classifyMin(q1, q2, q3, q12, b1, b2, print_Cases=False, print_Output=False
a2 = b2_star
type = 3 # 2
CaseCount += 1
- # TODO Problem? angle depends on how you choose?...
+ # 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
@@ -206,28 +183,9 @@ def classifyMin(q1, q2, q3, q12, b1, b2, print_Cases=False, print_Output=False
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.. ?
+ type = 3 # could check which axis: if a1_star or a2_star close to zero.. ?
CaseCount += 1
- # 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,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):
@@ -240,8 +198,8 @@ def classifyMin(q1, q2, q3, q12, b1, b2, print_Cases=False, print_Output=False
a2 = b2_star
type = 3 # 2
CaseCount += 1
- # TODO ...does type4 happen here?
- # TODO Problem? angle depends on how you choose?...
+
+ # 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
@@ -253,23 +211,23 @@ def classifyMin(q1, q2, q3, q12, b1, b2, print_Cases=False, print_Output=False
if determinant <= -1.0*epsilon:
if print_Cases: print('H : hyperbolic case (determinant smaller zero)')
# One or two minimizers wich are axial
- type = 3
+ 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
+ # 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
+ # 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
+ # type = 3 # 4
CaseCount += 1
# ---------------------------------------------------------------------------------------
@@ -277,7 +235,7 @@ def classifyMin(q1, q2, q3, q12, b1, b2, print_Cases=False, print_Output=False
print('Error: More than one Case happened!')
# compute a3
- a3 = math.sqrt(2.0*a1*a2) # ever needed?
+ # 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)))]
@@ -286,11 +244,10 @@ def classifyMin(q1, q2, q3, q12, b1, b2, print_Cases=False, print_Output=False
# compute kappa
kappa = (a1 + a2)
- #Test
+ # 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:")
@@ -299,7 +256,6 @@ def classifyMin(q1, q2, q3, q12, b1, b2, print_Cases=False, print_Output=False
print("type: ", type)
print("kappa = ", kappa)
-
return Minimizer, angle, type, kappa
# ------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
@@ -316,8 +272,8 @@ def classifyMin(q1, q2, q3, q12, b1, b2, print_Cases=False, print_Output=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)
+# 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)
@@ -350,12 +306,12 @@ def classifyMin(q1, q2, q3, q12, b1, b2, print_Cases=False, print_Output=False
# 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)
+# G, angle, type, kappa = classifyMin_ana(alpha, beta, theta, mu_gamma, q12, print_Cases, print_Output)
#
-# Out = classifyMin_geo(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_geo(alpha, beta, theta)
+# Out = classifyMin_ana(alpha, beta, theta)
# print('Out[0]', Out[0])
# print('Out[1]', Out[1])
@@ -364,7 +320,7 @@ def classifyMin(q1, q2, q3, q12, b1, b2, print_Cases=False, print_Output=False
# #supress certain Outout..
-# _,_,T,_ = classifyMin_geo(alpha, beta, theta, mu_gamma, q12, print_Cases, print_Output)
+# _,_,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)
diff --git a/src/PhaseDiagram.py b/src/PhaseDiagram.py
index 5fec95a9..486380a8 100644
--- a/src/PhaseDiagram.py
+++ b/src/PhaseDiagram.py
@@ -50,8 +50,7 @@ def GetMuGamma(beta,theta,gamma,mu1,rho1, InputFilePath = os.path.dirname(os.get
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)
- # print('gamma='+str(gamma))
+ # print("Run computeMuGamma for Gamma = ", gamma)
with open(InputFilePath, 'r') as file:
filedata = file.read()
filedata = re.sub('(?m)^gamma=.*','gamma='+str(gamma),filedata)
@@ -64,16 +63,19 @@ def GetMuGamma(beta,theta,gamma,mu1,rho1, InputFilePath = os.path.dirname(os.get
f.write(filedata)
f.close()
# --- Run Cell-Problem
- t = time.time()
+
+ # 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)
+ # 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:
@@ -115,8 +117,8 @@ 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
+# gamma = 'infinity'
+# gamma = 0.5
print('---- Input parameters: -----')
print('mu1: ', mu1)
@@ -128,24 +130,23 @@ print('gamma:', gamma)
print('----------------------------')
# ----------------------------------------------------------------
-muGamma = GetMuGamma(beta,theta,gamma,mu1,rho1,InputFilePath)
-# muGamma = GetMuGamma(beta,theta,gamma,mu1,rho1)
-
-print('Test MuGamma:', muGamma)
+# 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
+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 ---
@@ -184,8 +185,8 @@ make_2D_PhaseDiagram = True
SamplePoints_3D = 10 # Number of sample points in each direction
SamplePoints_2D = 10 # Number of sample points in each direction
-SamplePoints_3D = 40 # Number of sample points in each direction
-SamplePoints_2D = 3 # 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
@@ -200,15 +201,17 @@ if make_3D_PhaseDiagram:
# Get MuGamma values ...
GetMuGammaVec = np.vectorize(GetMuGamma)
- muGammas = GetMuGammaVec(betas,thetas,gamma,mu1,rho1)
+ muGammas = GetMuGammaVec(betas, thetas, gamma, mu1, rho1)
# Classify Minimizers....
- G, angles, Types, curvature = classifyMin_anaVec(alphas,betas,thetas,mu_gamma, mu1, rho1) # Sets q12 to zero!!!
+ 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
- VTKOutputName = "outputs/PhaseDiagram3D"
+ 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 )
@@ -221,20 +224,22 @@ if make_2D_PhaseDiagram:
# 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)
- classifyMin_anaVec = np.vectorize(classifyMin_ana)
-
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"
- VTKOutputName = "outputs/PhaseDiagram2D"
+ 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 )
@@ -244,9 +249,10 @@ 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_2D_plot: plt.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=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')
@@ -259,6 +265,8 @@ if(make_3D_plot or make_2D_plot):
# print('Type 2 occured here:', np.where(T == 2))
+print(alphas_)
+print(betas_)
# ALTERNATIVE
# colors = ("red", "green", "blue")
# groups = ("Type 1", "Type2", "Type3")
--
GitLab