Skip to content
Snippets Groups Projects
Commit 46900f16 authored by Klaus Böhnlein's avatar Klaus Böhnlein
Browse files

Update PhaseDiagram.py

parent 5203d9dd
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
src/ClassifyMin.py
+ 27
− 71
View file @ 46900f16
......@@ -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)
......
src/PhaseDiagram.py
+ 37
− 29
View file @ 46900f16
......@@ -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")
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment