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

Add classifyMin with Matrix Input & generalize PhaseDiagram to case of q12 != 0

parent 26ec8c07
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
src/ClassifyMin.py
+ 11
− 0
View file @ 2cdcb851
......@@ -75,6 +75,17 @@ def classifyMin_ana(alpha,beta,theta,q3,mu_1,rho_1,print_Cases=False, print_Outp
# 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 )
......
src/HelperFunctions.py
+ 7
− 2
View file @ 2cdcb851
......@@ -94,13 +94,18 @@ def RunCellProblem(alpha,beta,theta,gamma,mu1,rho1, InputFilePath = os.path.dirn
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)
# print('Q:', Q)
# print('B:', B)
return Q, B
......
src/PhaseDiagram.py
+ 154
− 135
View file @ 2cdcb851
......@@ -118,9 +118,9 @@ 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 = 'infinity'
# gamma = 0.5
# gamma = 0.25
gamma = 0.25
# gamma = 1.0
print('---- Input parameters: -----')
......@@ -133,150 +133,169 @@ 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
# gamma_min = 0.05
# gamma_max = 2.0
#TODO
# generalCase = False #Read Output from Cell-Problem instead of using Lemma1.4 (special case)
gamma_min = 1
gamma_max = 1
Gamma_Values = np.linspace(gamma_min, gamma_max, num=1)
# 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
# Gamma_Values = np.linspace(gamma_min, gamma_max, num=13) # TODO variable Input Parameters...alpha,beta...
print('(Input) Gamma_Values:', Gamma_Values)
for gamma in Gamma_Values:
# --- 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
# 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
#TODO
generalCase = True #Read Output from Cell-Problem instead of using Lemma1.4 (special case)
# 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
# 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 = 300 # 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)
# thetas_ = np.linspace(0.01,0.99,SamplePoints_2D)
alphas_ = np.linspace(8, 12, SamplePoints_2D)
thetas_ = np.linspace(0.05,0.2,SamplePoints_2D)
# betas_ = np.linspace(0.01,40.01,1)
#fix to one value:
betas_ = 2.0;
alphas, betas, thetas = np.meshgrid(alphas_, betas_, thetas_, indexing='ij')
# if generalCase: #TODO
# classifyMinVec = np.vectorize(classifyMin)
# GetCellOutputVec = np.vectorize(GetCellOutput)
# Q, B = GetCellOutputVec(alpha,betas,thetas,gamma,mu1,rho1,InputFilePath ,OutputFilePath )
#
# print('type of Q:', type(Q))
# print('Q:', Q)
# --- 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)
#
# else:
classifyMin_anaVec = np.vectorize(classifyMin_ana)
GetMuGammaVec = np.vectorize(GetMuGamma)
muGammas = GetMuGammaVec(betas,thetas,gamma,mu1,rho1,InputFilePath ,OutputFilePath )
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('---- 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 = 300 # Number of sample points in each direction
SamplePoints_2D = 5 # 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)
# thetas_ = np.linspace(0.01,0.99,SamplePoints_2D)
alphas_ = np.linspace(9, 10, SamplePoints_2D)
thetas_ = np.linspace(0.075,0.14,SamplePoints_2D)
# alphas_ = np.linspace(8, 12, SamplePoints_2D)
# thetas_ = np.linspace(0.05,0.2,SamplePoints_2D)
# betas_ = np.linspace(0.01,40.01,1)
#fix to one value:
betas_ = 2.0;
alphas, betas, thetas = np.meshgrid(alphas_, betas_, thetas_, indexing='ij')
if generalCase: #TODO
classifyMin_matVec = np.vectorize(classifyMin_mat)
GetCellOutputVec = np.vectorize(GetCellOutput)
Q, B = GetCellOutputVec(alphas,betas,thetas,gamma,mu1,rho1,InputFilePath ,OutputFilePath )
print('type of Q:', type(Q))
print('Q:', Q)
G, angles, Types, curvature = classifyMin_matVec(Q,B)
else:
classifyMin_anaVec = np.vectorize(classifyMin_ana)
GetMuGammaVec = np.vectorize(GetMuGamma)
muGammas = GetMuGammaVec(betas,thetas,gamma,mu1,rho1,InputFilePath ,OutputFilePath )
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_)
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment