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

Fix error in generalized form of phase Diagram

parent 2cdcb851
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
Matlab-Programs/classifyMIN.m
+ 1
− 1
View file @ 78ce2695
...@@ -54,8 +54,8 @@ b2 = (3*rho_1./(4.*((1-t)+t.*b))).*(1-t.*(1+b.*a)); ...@@ -54,8 +54,8 @@ b2 = (3*rho_1./(4.*((1-t)+t.*b))).*(1-t.*(1+b.*a));
CaseCount = 0; %check if more than one case ever happens CaseCount = 0; %check if more than one case ever happens
epsilon = eps;
epsilon = 1.e-18;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% PARABOLIC CASE %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% PARABOLIC CASE %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% if abs(det(A)) < epsilon * min(abs(det(A)),0) % if abs(det(A)) < epsilon * min(abs(det(A)),0)
......
src/PhaseDiagram.py
+ 135
− 134
View file @ 78ce2695
...@@ -137,14 +137,14 @@ print('----------------------------') ...@@ -137,14 +137,14 @@ print('----------------------------')
# gamma_min = 0.05 # gamma_min = 0.05
# gamma_max = 2.0 # gamma_max = 2.0
gamma_min = 1 # gamma_min = 1
gamma_max = 1 # gamma_max = 1
Gamma_Values = np.linspace(gamma_min, gamma_max, num=1) # Gamma_Values = np.linspace(gamma_min, gamma_max, num=1)
#
# Gamma_Values = np.linspace(gamma_min, gamma_max, num=13) # TODO variable Input Parameters...alpha,beta... # # Gamma_Values = np.linspace(gamma_min, gamma_max, num=13) # TODO variable Input Parameters...alpha,beta...
print('(Input) Gamma_Values:', Gamma_Values) # print('(Input) Gamma_Values:', Gamma_Values)
#
for gamma in Gamma_Values: # for gamma in Gamma_Values:
# muGamma = GetMuGamma(beta,theta,gamma,mu1,rho1,InputFilePath) # muGamma = GetMuGamma(beta,theta,gamma,mu1,rho1,InputFilePath)
...@@ -155,147 +155,148 @@ for gamma in Gamma_Values: ...@@ -155,147 +155,148 @@ for gamma in Gamma_Values:
# print_Cases = True # print_Cases = True
# print_Output = True # print_Output = True
#TODO #TODO
generalCase = True #Read Output from Cell-Problem instead of using Lemma1.4 (special case) 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
# --- Define effective quantities: q1, q2 , q3 = mu_gamma, q12 --- # make_3D_plot = True
# q1 = harmonicMean(mu1, beta, theta) # make_3D_PhaseDiagram = True
# q2 = arithmeticMean(mu1, beta, theta) make_2D_plot = False
# --- Set q12 # make_2D_PhaseDiagram = False
# q12 = 0.0 # (analytical example) # TEST / TODO read from Cell-Output 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------------------------------------------------------------------------------ # b1 = prestrain_b1(rho1, beta, alpha, theta)
# b2 = prestrain_b2(rho1, beta, alpha, theta)
# SamplePoints_3D = 10 # Number of sample points in each direction #
# SamplePoints_2D = 10 # Number of sample points in each direction # print('---- Input parameters: -----')
SamplePoints_3D = 300 # Number of sample points in each direction # print('mu1: ', mu1)
SamplePoints_2D = 5 # Number of sample points in each direction # print('rho1: ', rho1)
# print('alpha: ', alpha)
# print('beta: ', beta)
# print('theta: ', theta)
if make_3D_PhaseDiagram: # print("q1: ", q1)
alphas_ = np.linspace(-20, 20, SamplePoints_3D) # print("q2: ", q2)
betas_ = np.linspace(0.01,40.01,SamplePoints_3D) # print("mu_gamma: ", mu_gamma)
thetas_ = np.linspace(0.01,0.99,SamplePoints_3D) # print("q12: ", q12)
# print('type of alphas', type(alphas_)) # print("b1: ", b1)
# print('Test:', type(np.array([mu_gamma])) ) # print("b2: ", b2)
alphas, betas, thetas = np.meshgrid(alphas_, betas_, thetas_, indexing='ij') # 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 = 1 # 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, otypes=[np.ndarray, np.ndarray])
Q, B = GetCellOutputVec(alphas,betas,thetas,gamma,mu1,rho1,InputFilePath ,OutputFilePath )
# output = output.astype(object)
# print('type of Q:', type(Q))
# print('Q:', Q)
G, angles, Types, curvature = classifyMin_matVec(Q,B)
else:
classifyMin_anaVec = np.vectorize(classifyMin_ana) classifyMin_anaVec = np.vectorize(classifyMin_ana)
# Get MuGamma values ...
GetMuGammaVec = np.vectorize(GetMuGamma) GetMuGammaVec = np.vectorize(GetMuGamma)
muGammas = GetMuGammaVec(betas, thetas, gamma, mu1, rho1) muGammas = GetMuGammaVec(betas,thetas,gamma,mu1,rho1,InputFilePath ,OutputFilePath )
# Classify Minimizers.... G, angles, Types, curvature = classifyMin_anaVec(alphas,betas,thetas, muGammas, 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('size of G:', G.shape)
# print('G:', G) # print('G:', G)
# print('Types:', Types)
# Out = classifyMin_anaVec(alphas,betas,thetas) # Out = classifyMin_anaVec(alphas,betas,thetas)
# T = Out[2] # T = Out[2]
# --- Write to VTK # --- Write to VTK
GammaString = str(gamma) # VTKOutputName = + path + "./PhaseDiagram2DNEW"
VTKOutputName = "outputs/PhaseDiagram3D" + "Gamma" + GammaString GammaString = str(gamma)
gridToVTK(VTKOutputName , alphas, betas, thetas, pointData = {'Type': Types, 'angles': angles, 'curvature': curvature} ) VTKOutputName = "outputs/PhaseDiagram2D" + "Gamma_" + GammaString
print('Written to VTK-File:', VTKOutputName ) 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) # --- Make 3D Scatter plot
alphas_ = np.linspace(9, 10, SamplePoints_2D) if(make_3D_plot or make_2D_plot):
thetas_ = np.linspace(0.075,0.14,SamplePoints_2D) fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
# alphas_ = np.linspace(8, 12, SamplePoints_2D) colors = cm.plasma(Types)
# thetas_ = np.linspace(0.05,0.2,SamplePoints_2D) # if make_2D_plot: pnt3d=ax.scatter(alphas,thetas,c=Types.flat)
# betas_ = np.linspace(0.01,40.01,1) # if make_3D_plot: pnt3d=ax.scatter(alphas,betas,thetas,c=Types.flat)
#fix to one value: if make_2D_plot: pnt3d=ax.scatter(alphas,thetas,c=angles.flat)
betas_ = 2.0; if make_3D_plot: pnt3d=ax.scatter(alphas,betas,thetas,c=angles.flat)
alphas, betas, thetas = np.meshgrid(alphas_, betas_, thetas_, indexing='ij') # cbar=plt.colorbar(pnt3d)
# cbar.set_label("Values (units)")
ax.set_xlabel('alpha')
if generalCase: #TODO ax.set_ylabel('beta')
classifyMin_matVec = np.vectorize(classifyMin_mat) if make_3D_plot: ax.set_zlabel('theta')
GetCellOutputVec = np.vectorize(GetCellOutput) plt.show()
Q, B = GetCellOutputVec(alphas,betas,thetas,gamma,mu1,rho1,InputFilePath ,OutputFilePath ) # plt.savefig('common_labels.png', dpi=300)
# print('T:', T)
# print('Type 1 occured here:', np.where(T == 1))
print('type of Q:', type(Q)) # print('Type 2 occured here:', np.where(T == 2))
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_) # 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