Skip to content
Snippets Groups Projects

Compare revisions

Changes are shown as if the source revision was being merged into the target revision. Learn more about comparing revisions.

Source

Select target project
No results found

Target

Select target project
  • s7603593/dune-microstructure
  • s7603593/dune-microstructure-backup
2 results
Show changes
Commits on Source (74)
Hide whitespace changes
Inline Side-by-side
Showing
with 1030 additions and 0 deletions
.gitignore 0 → 100644
View file @ be6bdf92
# ignore all build folders
/build*/
build/
# ignore backup files
*~
# ignore python files
#*.pyc
#ignore kdevelop files
*.kdev4
CMakeLists.txt 0 → 100644
View file @ be6bdf92
# We require CMake version 3.1 to prevent issues
# with dune_enable_all_packages and older CMake versions.
cmake_minimum_required(VERSION 3.1)
project(dune-microstructure CXX)
if(NOT (dune-common_DIR OR dune-common_ROOT OR
"${CMAKE_PREFIX_PATH}" MATCHES ".*dune-common.*"))
string(REPLACE ${CMAKE_PROJECT_NAME} dune-common dune-common_DIR
${PROJECT_BINARY_DIR})
endif()
#find dune-common and set the module path
find_package(dune-common REQUIRED)
list(APPEND CMAKE_MODULE_PATH "${PROJECT_SOURCE_DIR}/cmake/modules"
${dune-common_MODULE_PATH})
#include the dune macros
include(DuneMacros)
# start a dune project with information from dune.module
dune_project()
dune_enable_all_packages()
add_subdirectory(src)
add_subdirectory(dune)
add_subdirectory(doc)
add_subdirectory(cmake/modules)
# finalize the dune project, e.g. generating config.h etc.
finalize_dune_project(GENERATE_CONFIG_H_CMAKE)
Matlab-Programs/Minimization.mlx 0 → 100644
View file @ be6bdf92
File added
Matlab-Programs/Minimization_Script.m 0 → 100644
View file @ be6bdf92
clear all
clc
% status = system('mkdir mynew')
% command = './build-cmake/src/dune-microstructure ./inputs/cellsolver.parset';
% system(['set PATH=' '/home/klaus/Desktop/DUNE/dune-microstructure' ';' command ]);
% --- Change PolynomialDisplayStyle ----
% sympref('PolynomialDisplayStyle','ascend');
% sympref('AbbreviateOutput',false);
syms f_plus(v1,v2,q1,q2,q3,b1,b2,b3)
assume( q1 > 0)
assume( q2 > 0)
assume( q3 > 0)
assume( q3 >= q1)
assume( q2 >= q3)
v = [v1; v2];
%should be sqrt(2) instead of 2
fprintf('Functions to be minimized')
f_plus(v1,v2,q1,q2,q3,b1,b2,b3) = q1*v1^4 + q2*v2^4+2*q3*v1^2*v2^2-2*(q1*b1*v1^2+ q2*b2*v2^2+sqrt(2)*q3*b3*v1*v2);
f_minus(v1,v2,q1,q2,q3,b1,b2,b3) = q1*v1^4 + q2*v2^4+2*q3*v1^2*v2^2+2*(q1*b1*v1^2+ q2*b2*v2^2+sqrt(2)*q3*b3*v1*v2);
% ---- Fix parameters
% 1. Import effective quantities from CellSolver-Code:
Qmat = spconvert(load('QMatrix.txt'));
Qmat = full(Qmat)
Bmat = spconvert(load('BMatrix.txt'));
Bmat = full(Bmat)
% Substitute effective quantitites
f_plus = subs(f_plus,q1,Qmat(1,1));
f_plus = subs(f_plus,q3,Qmat(3,3));
f_plus = subs(f_plus,q2,Qmat(2,2));
f_plus = subs(f_plus,b1,Bmat(1));
f_plus = subs(f_plus,b2,Bmat(2));
f_plus = subs(f_plus,b3,Bmat(3));
% f_plus = subs(f_plus,b3,0);
%
f_minus = subs(f_minus,q1,Qmat(1,1));
f_minus = subs(f_minus,q3,Qmat(3,3));
f_minus = subs(f_minus,q2,Qmat(2,2));
f_minus = subs(f_minus,b1,Bmat(1));
f_minus = subs(f_minus,b2,Bmat(2));
f_minus = subs(f_minus,b3,Bmat(3));
% % f_minus = subs(f_minus,b3,0);
% 2. Substitute specific values:
%%%Compare with 'ClassifyMin-Code'
% % Compare with 'ClassifyMin-Code'
% mu_1 = 1;
% rho_1 = 1;
% % --- type 1 Situation:
% % beta = 2;
% % alpha = 2;
% % theta = 1/4;
% % --- type 2 Situation:
% % beta = 3.0714;
% % alpha = -20;
% % theta = 0.3;
% % --- type 3 Situation:
% % beta = 2.2857;
% % alpha = -20;
% % theta = 0.3;
%
% % interesting:
% alpha = 18.3673;
% beta = 8.57143;
% theta= 0.913265;
%
% % alpha = 2.85714;
% % beta = 7;
% % theta= 0.3;
%
%
%
% set_mu_gamma = 'q1';
% % set_mu_gamma = 'q2';
% print_output = false;
%
% q1i = mu_1.*(beta./(theta+(1-theta).*beta))
% q2i = mu_1.*((1-theta)+theta.*beta)
% q3i = q1i
% % b1i = (mu_1*rho_1/4).*(beta./(theta+(1-theta).*beta)).*(1-theta.*(1+alpha))
% % b2i = mu_1.*(rho_1/8).*(1-theta.*(1+beta.*alpha))
% b3i = 0
%
% %TEST (new)
% b1i = (3*rho_1/2).*beta.*(1-theta.*(1+alpha));
% b2i = (3*rho_1/(4*((1-theta)+theta.*beta))).*(1-theta.*(1+beta.*alpha));
%
% f_plus = subs(f_plus,q1,q1i);
% f_plus = subs(f_plus,q3,q3i);
% f_plus = subs(f_plus,q2,q2i);
% f_plus = subs(f_plus,b1,b1i);
% f_plus = subs(f_plus,b2,b2i);
% f_plus = subs(f_plus,b3,b3i)
%
% f_minus = subs(f_minus,q1,q1i);
% f_minus = subs(f_minus,q3,q3i);
% f_minus = subs(f_minus,q2,q2i);
% f_minus = subs(f_minus,b1,b1i);
% f_minus = subs(f_minus,b2,b2i);
% f_minus = subs(f_minus,b3,b3i)
%
%
%
% [A,angle,V] = classifyMIN(mu_1,rho_1,alpha,beta,theta,set_mu_gamma,print_output)
% Substitute random values...
% rndm1 = randi([1 20],1,1);
% rndm2 = randi([1 20],1,1);
% rndm3 = randi([1 20],1,1);
% f_plus = subs(f_plus,b1,rndm1);
% f_plus = subs(f_plus,b2,rndm2);
% f_plus = subs(f_plus,b3,rndm3);
% Compute the Gradients
df_plusx = diff(f_plus,v1);
df_plusy = diff(f_plus,v2);
df_minusx = diff(f_minus,v1);
df_minusy = diff(f_minus,v2);
% Setup Equations Grad(f) = 0
eq1 = df_plusx == 0;
eq2 = df_plusy == 0;
eqns = [eq1, eq2];
eq3 = df_minusx == 0;
eq4 = df_minusy == 0;
eqns_minus = [eq3, eq4];
% Symbolically Solve Equations:
% More robust (works even for values b_3 ~ 1e-08 ):
S = solve(eqns,v1,v2,'MaxDegree' , 5);
S_minus = solve(eqns_minus,v1,v2,'MaxDegree' , 5);
%Tests:
% S = solve(eqns,v1,v2,'MaxDegree' , 5, 'Real', true);
% S_minus = solve(eqns_minus,v1,v2,'MaxDegree' , 5, 'Real', true);
% S = solve(eqns,v1,v2,'MaxDegree' , 5, 'IgnoreAnalyticConstraints',true);
% S = solve(eqns,v1,v2,'MaxDegree' , 5, 'IgnoreAnalyticConstraints',true, 'Real', true);
% S = solve(eqns)
A = S.v1;
B = S.v2;
A_minus = S_minus.v1;
B_minus = S_minus.v2;
% A = simplify(A);
% B = simplify(B)
%---------- TEST if Grad(f) = 0 ---------------------
% fprintf('Testing equation grad(f) = 0 with stationary points')
%
% for i = 1:size(A,1)
% fprintf('Testing %d.point (f_plus): ',i )
% [ double(subs(subs(df_plusx,v1,A(i)),v2,B(i))), double(subs(subs(df_plusy,v1,A(i)),v2,B(i))) ]
% end
% for i = 1:size(A_minus,1)
% fprintf('Testing %d.point (f_minus): ',i )
% [double(subs(subs(df_minusx,v1,A_minus(i)),v2,B_minus(i))), double(subs(subs(df_minusy,v1,A_minus(i)),v2,B_minus(i)))]
% end
% ------------------------------------
fprintf('stationary points of f_plus:')
A = double(A); %safe symbolic values
B = double(B);
fprintf('stationary points of f_minus:')
A_minus = double(A_minus);
B_minus = double(B_minus);
% Extract only Real-Solutions
fprintf('real stationary points of f_plus:')
tmp1 = A(imag(A)==0 & imag(B) == 0);
tmp2 = B(imag(A)==0 & imag(B) == 0);
A = tmp1;
B = tmp2;
% A(abs(imag(A)) <1e-3 & abs(imag(B)) <1e-3 )
SP_Plus = [A,B]
fprintf('real stationary points of f_minus:')
tmp1 = A_minus(imag(A_minus)==0 & imag(B_minus) == 0);
tmp2 = B_minus(imag(A_minus)==0 & imag(B_minus) == 0);
A_minus = tmp1;
B_minus = tmp2;
% A_minus(abs(imag(A_minus)) <1e-3 & abs(imag(B_minus)) <1e-3 )
SP_Minus = [A_minus,B_minus]
% Determine global Minimizer from stationary points:
fprintf('function values at stationary points (f_plus):')
T = arrayfun(@(v1,v2) double(f_plus(v1,v2,q1,q2,q3,b1,b2,b3)),A,B,'UniformOutput', false)
T = cell2mat(T);
% Min_plus = min(T, [], 'all')
[Min_plus,MinIdx_plus] = min(T, [], 'all', 'linear');
fprintf('function values at stationary points (f_minus):')
T_minus = arrayfun(@(v1,v2) double(f_minus(v1,v2,q1,q2,q3,b1,b2,b3)),A_minus,B_minus,'UniformOutput', false)
T_minus = cell2mat(T_minus);
% Min_minus = min(T_minus, [], 'all')
[Min_minus,MinIdx_minus] = min(T_minus, [], 'all', 'linear');
[globalMinimizerValue,GlobalIdx] = min([Min_plus,Min_minus]);
if GlobalIdx == 1 %Min_plus
GlobalMinimizer = SP_Plus(MinIdx_plus,:);
sign = 1.0;
elseif GlobalIdx == 2 %Min_minus
GlobalMinimizer = SP_Minus(MinIdx_minus,:);
sign = -1.0;
end
fprintf('Global Minimizer:(%d,%d)', GlobalMinimizer(1),GlobalMinimizer(2) )
fprintf('Global Minimizer Value : %d', globalMinimizerValue )
% Plot functions
fsurf(@(x,y) f_plus(x,y,q1,q2,q3,b1,b2,b3))
hold on
plot3(A,B,T,'g*')
%Plot GlobalMinimizer:
hold on
plot3(GlobalMinimizer(1),GlobalMinimizer(2),globalMinimizerValue, 'o', 'Color','c')
% view(90,0)
% view(2)
figure
fsurf(@(x,y) f_minus(x,y,q1,q2,q3,b1,b2,b3))
hold on
plot3(A_minus,B_minus,T_minus,'g*')
hold on
plot3(GlobalMinimizer(1),GlobalMinimizer(2),globalMinimizerValue, 'o', 'Color','c')
%Write to txt-File
fileID = fopen('txt.txt','w');
fprintf(fileID,'%s' , latex(S.v1));
fclose(fileID);
%%%Compare with 'ClassifyMin-Code'
fprintf('----------------compare with ClassifyMIN----------------')
fprintf('Output Minimizer Matrix from symbolic Minimization')
sign*GlobalMinimizer'*GlobalMinimizer %sign correct? should do this with symbolic Values! TODO
% GlobalMinimizer'*GlobalMinimizer
%check with higher Precision:
% vpa(GlobalMinimizer'*GlobalMinimizer)
% %
% % %Output from Classify Min:
% % [A,angle,type,K] = classifyMIN(mu_1,rho_1,alpha,beta,theta,set_mu_gamma,print_output);
% % fprintf('Output Minimizer Matrix from ClassifyMIN')
% %
% % % [A(1) sign*sqrt(A(1)*A(2)) ; sign*sqrt(A(1)*A(2)) A(2)] %sign correct?
% % [A(1) sqrt(A(1)*A(2)) ; sqrt(A(1)*A(2)) A(2)] %sign correct?
% %
% % %check with higher Precision:
% % % vpa([A(1) sqrt(A(1)*A(2)) ; sqrt(A(1)*A(2)) A(2)])
% %
% %
% % e = [sqrt(A(1)), sqrt(A(2))]; %TODO .. this might be complex?!
% %
% % norm = sqrt((A(1)+A(2)));
% %
% % e = e./norm;
% %
% % K*(e'*e)
% %
% % % e'*e
% % % K
% %
% % fprintf('angle: %d', angle)
% % fprintf('Type: %d', type)
%% Compare with "Task2"
% fprintf('----------------compare with Task2----------------')
%
% B = [b1i b3i; b3i b2i];
% x = 0:0.01:2*pi;
%
% y1 = arrayfun(@(alpha)compute_F(alpha,B,q1i,q2i,q3i),x,'UniformOutput', false);
% y1 = cell2mat(y1);
%
%
% figure
% plot(x,y1,'b')
% hold on
%
% fun = @(a) compute_F(a,B,q1i,q2i,q3i);
% [alphaMin,f] = fminunc(fun,0)
% [alphaMin,f] = fminunc(fun,3) % Different initial value
% plot(alphaMin,f,'*')
%
% %compute radius
% rMin = compute_r(alphaMin,B,q1i,q2i,q3i)
%
% %compute Minimizer:
% v_alpha = [cos(alphaMin);sin(alphaMin)];
%
%
%
% G = rMin.*(v_alpha*v_alpha')
%%Determine Minimizer Type (in general case)
% % T = [GlobalMinimizer' e1']
% % det(T) % also works?
% %
% %
% % % symbolically :
% %
% % if GlobalIdx == 1 %Min_plus
% % A_sym = S.v1;
% % B_sym = S.v2
% % Index = MinIdx_plus;
% % elseif GlobalIdx == 2 %Min_minus
% % A_sym = S_minus.v1;
% % B_sym = S_minus.v2;
% % Index = MinIdx_minus;
% % end
% %
% % % Check Determinant symbolically?!?!
% %
% % g_sym = [A_sym(Index) B_sym(Index)]
% % G_sym = g_sym'*g_sym
% %
% % e1 = [1 0];
% % e2 = [0 1];
% %
% % % check alignment with e1
% % % if ....
% % det([g_sym' e1'])
% % % ... bending in e1 direction
% % % check alignment with e2
% % % if..
% % det([g_sym' e2'])
% % double(det([g_sym' e2']))
% % % bending in e2 direction
% % %Else
% % %....
Matlab-Programs/Minimization_abstract.mlx 0 → 100755
View file @ be6bdf92
File added
Matlab-Programs/PhaseDiagrams.mlx 0 → 100644
View file @ be6bdf92
File added
Matlab-Programs/PhaseDiagramsScript.m 0 → 100644
View file @ be6bdf92
clear all
clc
% INPUT Parameters
mu_1 = 1;
rho_1 = 1;
theta= 0.3;
alpha = 2;
beta= 2;
fprintf('===========================================================================================================================');
fprintf(' Possible cases for global minimizers: (I) , (II) , (III) , (IV) ')
fprintf('============================================== INPUT PARAMETERS ===========================================================');
fprintf('mu_1: %d rho:1: %d theta: %d', mu_1,rho_1,theta );
fprintf('===========================================================================================================================');
% choose u_gamma = q3 to be either q1, q2 or something in between
% set_mu_gamma = 'q1';
set_mu_gamma = 'q2'; %(hyperbolic-case)
% set_mu_gamma = 'm'; % mean of q1,q2
% print_output = true;
print_output = false;
% Test for fixed value
% [A,angle,type] = classifyMIN(mu_1,rho_1,alpha,beta,theta,set_mu_gamma,print_output);
% PLOT
x = linspace(-20,20,45); %~alpha
y = linspace(1.5,40,45); %~beta
z = linspace(0.05,0.95,45); %~theta
[X1,Y1] = meshgrid(x,y);
[A_2D,angle_2D,V_2D] = arrayfun(@(a,b)classifyMIN(mu_1,rho_1,a,b,theta,set_mu_gamma,print_output),X1,Y1,'UniformOutput', false);
[X,Y,Z] = meshgrid(x,y,z);
[A,angle_3D,V_3D] = arrayfun(@(a,b,theta)classifyMIN(mu_1,rho_1,a,b,theta,set_mu_gamma,print_output),X,Y,Z,'UniformOutput', false);
color_3D = cell2mat(V_3D);
color_2D = cell2mat(V_2D);
angle_2D = cell2mat(angle_2D);
angle_3D = cell2mat(angle_3D);
A = cell2mat(A);
A_2D = cell2mat(A_2D);
X1 = reshape(X1,[],1);
Y1 = reshape(Y1,[],1);
X = reshape(X,[],1);
Y = reshape(Y,[],1);
Z = reshape(Z,[],1);
color_2D = reshape(color_2D,[],1);
color_3D = reshape(color_3D,[],1);
angle_2D = reshape(angle_2D,[],1);
angle_3D = reshape(angle_3D,[],1);
% Structure result depending on Type/Color
% V_2D = reshape(V_2D,[],1);
V_2DT1 = (color_2D == 1);
V_2DT2 = (color_2D == 2);
V_2DT3 = (color_2D == 3);
V_2DT1 = reshape(V_2DT1,[],1);
V_2DT2 = reshape(V_2DT2,[],1);
V_2DT3 = reshape(V_2DT3,[],1);
X1T1 = V_2DT1.*X1;
Y1T1 = V_2DT1.*Y1;
X1T2 = V_2DT2.*X1;
Y1T2 = V_2DT2.*Y1;
X1T3 = V_2DT3.*X1;
Y1T3 = V_2DT3.*Y1;
%%% 2D - Plot (fixed Theta)
cm_2D = redblue(3);
scatter(X1T1,Y1T1,[], 'MarkerFaceColor',[0 0 1 ], 'MarkerEdgeColor','black');
colormap(cm_2D)
hold on
scatter(X1T2,Y1T2,[], 'MarkerFaceColor',[1 0 0], 'MarkerEdgeColor','black');
hold on
% scatter(X1T3,Y1T3,[],'MarkerFaceColor',[0 0 1 ], 'MarkerEdgeColor','black');
scatter(X1T3,Y1T3,[], angle_2D, 'filled','MarkerEdgeColor','black');
colormap(cm_2D)
legend('Type 1', 'Type 2', 'Type 3')
title('Fixed Theta Plot')
xlabel('alpha')
ylabel('beta')
hold off
%%% 3D - Plot
V_3DT1 = (color_3D == 1);
V_3DT2 = (color_3D == 2);
V_3DT3 = (color_3D == 3);
V_3DT1 = reshape(V_3DT1,[],1);
V_3DT2 = reshape(V_3DT2,[],1);
V_3DT3 = reshape(V_3DT3,[],1);
XT1 = V_3DT1.*X;
YT1 = V_3DT1.*Y;
ZT1 = V_3DT1.*Z;
XT2 = V_3DT2.*X;
YT2 = V_3DT2.*Y;
ZT2 = V_3DT2.*Z;
XT3 = V_3DT3.*X;
YT3 = V_3DT3.*Y;
ZT3 = V_3DT3.*Z;
cm = redblue(90);
figure
%fixed Color
% scatter3(XT1,YT1,ZT1, [], 'MarkerFaceColor',[0.75 0 0],'MarkerEdgeColor', 'none');
%variing Color
scatter3(YT1,ZT1,XT1, [], 'MarkerFaceColor',[0 0 1 ], 'MarkerEdgeColor','black');
colormap(cm);
hold on
scatter3(YT2,ZT2,XT2, [], 'MarkerFaceColor',[1 0 0],'MarkerEdgeColor', 'none');
hold on
% scatter3(XT3,YT3,ZT3, [],'MarkerFaceColor',[0 0 1 ],'MarkerEdgeColor', 'none');
scatter3(YT3,ZT3,XT3, [], angle_3D, 'filled');
legend('Type 1', 'Type 2', 'Type 3')
title('Classification of Minimizers, theta:')
% xlabel('alpha')
% ylabel('beta')
% zlabel('theta')
xlabel('beta')
ylabel('theta')
zlabel('alpha')
Matlab-Programs/SamplePlot2.mlx 0 → 100644
View file @ be6bdf92
File added
Matlab-Programs/Task2.mlx 0 → 100755
View file @ be6bdf92
File added
Matlab-Programs/TestCompute.m 0 → 100644
View file @ be6bdf92
mu_1 = 10;
rho_1 = 1;
theta = 0.05;
% alpha = -0.5;
% alpha = 0.5;
% alpha = 1;
alpha = 18;
beta = 0.5; %egal welcher Wert Beta ... alha = 1 & Theta = 0.5 -> type 2 ?
beta = 2;
% Interessant:
alpha = 10;
% alpha = 18;
theta = 0.05;
% unabhänhig von beta?
beta = 2;
% Compute components of B_eff
% b1 = (mu_1*rho_1/4).*(beta./(theta+(1-theta).*beta)).*(1-theta.*(1+alpha));
% b2 = mu_1.*(rho_1/8).*(1-theta.*(1+beta.*alpha));
%TEST (new)
b1 = (3*rho_1/2).*beta.*(1-theta.*(1+alpha));
b2 = (3*rho_1/(4*((1-theta)+theta.*beta))).*(1-theta.*(1+beta.*alpha));
mu_h = @(beta,theta) mu_1.*(beta./(theta+(1-theta).*beta)); % harmonic mean
mu_bar = @(b,t) mu_1.*((1-theta)+theta.*beta); % mu_bar
% q1 q2 q3..
q1 = mu_h(beta,theta);
q2 = mu_bar(beta,theta);
fprintf('--------------------')
fprintf(' alpha:%d' , alpha)
fprintf(' beta:%d ' , beta)
fprintf(' theta:%d ', theta )
fprintf('-------------------- \n')
fprintf('q1*b1^2:')
q1*b1^2
fprintf('q2*b2^2:')
q2*b2^2
fprintf('Test')
(9/4)*(beta/(1+beta))*((1/2)-alpha+(1/2)*alpha^2)
\ No newline at end of file
Matlab-Programs/classifyMIN.m 0 → 100755
View file @ be6bdf92
function [A, angle, type, kappa] = classifyMIN (mu_1,rho_1,a,b,t,set_mu_gamma,print_output)
% returns
% A : Matrix of basis coefficients [a1,a2,a3]
%
% type :
% Type of minimizer 1 = (I) , 2 = (II) , 3 = (III) , 4 = (IV)
%
type = 0; % either 1,2,3,4
mu_h = @(b,t) mu_1.*(b./(t+(1-t).*b)); % harmonic mean
mu_bar = @(b,t) mu_1.*((1-t)+t.*b); % mu_bar
if (set_mu_gamma == 'q1')
mu_gamma = @(b,t) mu_h(b,t);
end
if (set_mu_gamma == 'q2')
mu_gamma = @(b,t) mu_bar(b,t);
end
if (set_mu_gamma == 'm')
mu_gamma = @(b,t) 0.5*(mu_h(b,t) + mu_bar(b,t));
end
% q1 q2 q3..
q1 = mu_h(b,t);
q2 = mu_bar(b,t);
q3 = mu_gamma(b,t);
% values for q1,q2,q3 should be positiv
% assert((q1 > 0 ) & (q2 > 0 ) & (q3 > 0), 'At least one of q1,q2 or q3 is not positive' )
% Compute components of B_eff (old)
% b1 = (mu_1*rho_1/4).*(b./(t+(1-t).*b)).*(1-t.*(1+a));
% b2 = mu_1.*(rho_1/8).*(1-t.*(1+b.*a));
%TEST (new)
b1 = (3*rho_1/2).*b.*(1-t.*(1+a));
b2 = (3*rho_1./(4.*((1-t)+t.*b))).*(1-t.*(1+b.*a));
% H = [q1 q3; q3 q2];
%check condition of H first
% fprintf('condition number of Matrix H: %d \n', cond(H));
CaseCount = 0; %check if more than one case ever happens
epsilon = 1.e-18;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% PARABOLIC CASE %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% if abs(det(A)) < epsilon * min(abs(det(A)),0)
if abs(q1*q2-q3^2) < epsilon
fprintf('determinant equal zero (parabolic case)')
fprintf('SHOULD NOT HAPPEN')
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% ELLIPTIC CASE %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if (q1*q2-q3^2 > epsilon)
% fprintf('determinant greater than zero (elliptic case)');
% a1_star = (q2(b,t).*q1(b,t).*b1(a,b,t) - q3(b,t).*q2(b,t).*b2(a,b,t))./(q1(b,t).*q2(b,t)-q3(b,t).^2);
% a2_star = (q1(b,t).*q2(b,t).*b2(a,b,t) - q3(b,t).*q1(b,t).*b1(a,b,t))./(q1(b,t).*q2(b,t)-q3(b,t).^2);
a1_star = (q2.*q1.*b1 - q3.*q2.*b2)./(q1.*q2-q3.^2);
a2_star = (q1.*q2.*b2 - q3.*q1.*b1)./(q1.*q2-q3.^2);
prod = a1_star*a2_star;
if(prod > epsilon) % (E1) inside Lamba %
% (a1_star,a2_star) is unique minimizer lies inside Lambda
% therefore Minimizer not aligned with axes
% fprintf('\n elliptic-case: (E1)');
a1 = a1_star;
a2 = a2_star;
type = 3;
CaseCount = CaseCount + 1;
end
% Make distinction between boundary & outside (prod < 0 )
if(abs(prod) < epsilon) % (E2) on boundary of Lambda
% fprintf('\n elliptic-case: (E2)');
% Uniqueness of gloal minimizer if lies on boundary if prod = 0
% ----- %
% global minimizer lies on the boundary of Lambda depending on
% condition:
if (q2*b2^2 < q1*b1^2) % global Minimizer given by (b1,0)
a1 = b1;
a2 = 0*b1;
type = 1; % Minimizer aligned with x1-axis
CaseCount = CaseCount + 1;
end
if (q2*b2^2 > q1*b1^2)% global Minimizer given by (0,b2)
a1 = 0*b1;
a2 = b2;
type = 2; % Minimizer aligned with x2-axis
CaseCount = CaseCount + 1;
end
if abs(q1*b1^2-q2*b2^2) < epsilon * min(q2*b2^2,q1*b1^2) % q1b1^2 = q2b2^2
% two Minimizers ..pick one
a1 = b1;
a2 = 0*b1;
type = 4;
CaseCount = CaseCount + 1;
end
end
if((prod) < -1*epsilon) %Outside of Lambda
if (q2*b2^2 < q1*b1^2) % global Minimizer given by (b1,0)
a1 = b1;
a2 = 0*b1;
type = 1; % Minimizer aligned with x1-axis
CaseCount = CaseCount + 1;
end
if (q2*b2^2 > q1*b1^2)% global Minimizer given by (0,b2)
a1 = 0*b1;
a2 = b2;
type = 2; % Minimizer aligned with x2-axis
CaseCount = CaseCount + 1;
end
end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% HYPERBOLIC CASE %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if (q1*q2-q3^2 < -1*epsilon)
% fprintf('determinant less than zero (hyperbolic case)');
if (q2*b2^2 < q1*b1^2) % global Minimizer given by (b1,0)
a1 = b1;
a2 = 0*b1;
type = 1; % Minimizer aligned with x1-axis
CaseCount = CaseCount + 1;
end
if (q2*b2^2 > q1*b1^2)% global Minimizer given by (0,b2)
a1 = 0*b1;
a2 = b2;
type = 2; % Minimizer aligned with x2-axis
CaseCount = CaseCount + 1;
end
if abs(q1*b1^2-q2*b2^2) < epsilon * min(q2*b2^2,q1*b1^2) % q1b1^2 = q2b2^2
% two Minimizers ..pick one
a1 = b1;
a2 = 0*b1;
type = 4;
CaseCount = CaseCount + 1;
end
% CAN NOT BE TYPE 3!!
end
% Compute a3 from a1 % a2
a3 = sqrt(2*a1*a2); % WRONG ?
% compute angle between [sqrt(a1) , sqrt(a2)] and e1:
% angle = atan2(sqrt(a2),sqrt(a1));
v = [sqrt(a1), sqrt(a2)];
e1 = [1 0];
if (type == 3 )
% angle = atan2(a2,a1);
% angle = acos(v*e1/(sqrt(a1+a2)));
else
angle = 0;
end
% angle = atan2(norm(cross(a,b)), dot(a,b))
% Try to compute angle of Matrices
E1 = [1 0;0 0 ];
% V = [a1 sqrt(a1*a2)/sqrt(2); sqrt(a1*a2)/sqrt(2) a2]; % This ?? NO
V = [a1 sqrt(a1*a2); sqrt(a1*a2) a2]; % This!
% CHECK
% U = trace(V'*E1)/(sqrt(trace(V'*V))*sqrt(trace(E1'*E1)))
% if abs(U) > 1
% fprintf('value greater 1')
% end
% angle = atan2(sqrt(abs(a2)), sqrt(abs(a1))); % does this make sense?
% angle = acos(trace(V'*E1)/(sqrt(trace(V'*V))*sqrt(trace(E1'*E1)))); %angle in radians
% angle = acos(trace(V'*E1)/(sqrt(trace(V'*V))*sqrt(trace(E1'*E1))))/pi * 180; % angle in degrees
% CHECK does case (0,-b) ever happen? Yes
% if (q2*b2^2 > q1*b1^2 && b2 < 0)
% fprintf('point lies on negative half of y-axis');
% end
% CHECK if Minimizer ever lies inside lambda ... YES
% if (a1 ~= 0 && a2 ~= 0)
% fprintf('minimizer lies inside lambda');
% end
% Alternative atan2d(y,x): returns angle in degrees
% CHECK
e = [sqrt(a1), sqrt(a2)]; % Might be complex...
norm = sqrt(a1+a2); % Might be complex...
e = e./norm;
angle = atan2(e(2), e(1)); %TEST
if (imag(e(1))~= 0 || imag(e(2))~=0)
fprintf('complex vector e');
end
% if (imag(norm)~=0)
% fprintf('complex norm'); % happens a lot ..
% end
% do it this way to avoid complex numbers:
e = [sqrt((a1/(a1+a2))), sqrt((a2/(a1+a2)))]; % always positive under sqrt here .. basically takes absolute value here
angle = atan2(e(2), e(1));
if (imag(e(1))~= 0 || imag(e(2))~=0)
fprintf('complex vector e');
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% TEST -------------------
% the magnitude kappa of the vector does not matter for the angle ..
% In order to also get negative angles just use:
% atan2(a2,a1)
% angle = atan2(a2,a1); %% FEHLER : muss [sqrt(a1) sqrt(a2)] betrachten..
% angle = atan2(sqrt(a2),sqrt(a1)); % Inputs must be real...
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% angle2 = acos(e(1)/(sqrt(e(1)^2 + e(2)^2) ));
% if angle ~= angle2
% interesting = 1;
% end
% if (angle ~= 0 | angle ~= 3.1416)
% fprintf('angle not zero or pi.. plot type')
% type
% angle
% end
%compute Kappa?
% fprintf('Output Kappa:')
% k = sqrt(abs(a1) + abs(a2)); % ?
kappa = (a1 + a2);
% e = [sqrt(a1), sqrt(a2)];
% norm = sqrt((a1+a2));
% e = e./norm;
% K.*(e'*e);
if (CaseCount > 1)
fprintf('=============================== ERROR========================================= \n')
fprintf('Error: More than one Case happened!')
end
% Coefficients of minimizer
if(print_output)
fprintf(' \n ')
fprintf('=============================== OUTPUT ========================================= \n')
fprintf('Minimizer is of Type: %d \n' , type);
fprintf('Coefficients a1,a2,a3 given by : %d, %d, %d \n', a1, a2, a3);
fprintf('================================================================================ \n')
end
A = [a1, a2, a3];
end
Matlab-Programs/compute_F.m 0 → 100755
View file @ be6bdf92
function F = compute_F(alpha,B,q1,q2,q3)
% r = compute_r(alpha,B,q1,q2,q3);
v = [cos(alpha);sin(alpha)];
b1 = B(1,1);
b2 = B(2,2);
b3 = B(1,2);
%compute Q(v_alpha x v_alpha)
Q = q1.*v(1).^4 + q2.*v(2).^4 + 2.*q3.*v(1).^2.*v(2).^2;
%
% TP = v*v';
% L = stVenant(TP,mu,lambda);
tmp1 = q1.*(v(1).^2+b1).^2 + q2.*(v(2).^2+b2).^2 + q3.*(sqrt(2)*v(1)*v(2)+b3).^2;
tmp2 = q1.*b1.^2 + q2.*b2.^2+ q3.*b3.^2;
L = 0.5*(tmp1-Q-tmp2) ; %Polarization identity
r = L./Q;
% F = r.^2.*Q - 2.*r.*trace(L'*B)
F = (r.^2).*Q - 2.*r.*L;
end
Matlab-Programs/plotBeff.mlx 0 → 100644
View file @ be6bdf92
File added
Matlab-Programs/redblue.m 0 → 100755
View file @ be6bdf92
function c = redblue(m)
%REDBLUE Shades of red and blue color map
% REDBLUE(M), is an M-by-3 matrix that defines a colormap.
% The colors begin with bright blue, range through shades of
% blue to white, and then through shades of red to bright red.
% REDBLUE, by itself, is the same length as the current figure's
% colormap. If no figure exists, MATLAB creates one.
%
% For example, to reset the colormap of the current figure:
%
% colormap(redblue)
%
% See also HSV, GRAY, HOT, BONE, COPPER, PINK, FLAG,
% COLORMAP, RGBPLOT.
% Adam Auton, 9th October 2009
if nargin < 1, m = size(get(gcf,'colormap'),1); end
if (mod(m,2) == 0)
% From [0 0 1] to [1 1 1], then [1 1 1] to [1 0 0];
m1 = m*0.5;
r = (0:m1-1)'/max(m1-1,1);
g = r;
r = [r; ones(m1,1)];
g = [g; flipud(g)];
b = flipud(r);
else
% From [0 0 1] to [1 1 1] to [1 0 0];
m1 = floor(m*0.5);
r = (0:m1-1)'/max(m1,1);
g = r;
r = [r; ones(m1+1,1)];
g = [g; 1; flipud(g)];
b = flipud(r);
end
c = [r g b];
Matlab-Programs/resources/addons_core.xml 0 → 100755
View file @ be6bdf92
<?xml version="1.0" encoding="UTF-8"?>
<addonCore>
<label>Red Blue Colormap</label>
<version>1.0.0.0</version>
<type>zip</type>
<identifier>e5698820-4a80-11e4-9553-005056977bd0</identifier>
<createdBy name="Adam Auton"/>
<image>resources/screenshot.png</image>
</addonCore>
Matlab-Programs/resources/matlab_path_entries.xml 0 → 100755
View file @ be6bdf92
<?xml version="1.0" encoding="UTF-8"?>
<paths>
<path>.</path>
</paths>
Matlab-Programs/resources/metadata.xml 0 → 100755
View file @ be6bdf92
<?xml version="1.0" encoding="UTF-8"?>
<addOn>
<identifier>e5698820-4a80-11e4-9553-005056977bd0</identifier>
<displayType>Function</displayType>
<translatedDisplayType>
<en_US>Function</en_US>
<ja_JP>関数</ja_JP>
<ko_KR>함수</ko_KR>
<zh_CN>函数</zh_CN>
</translatedDisplayType>
<name>Red Blue Colormap</name>
<author>Adam Auton</author>
<version>1.0.0.0</version>
<downloadUrl>https://www.mathworks.com/matlabcentral/mlc-downloads/downloads/submissions/25536/versions/1/download/zip?src=addons_ml_desktop_install&amp;profile_id=14169257&amp;license=40758619&amp;release_family=R2020b</downloadUrl>
<licenseUrl>https://addons.mathworks.com/registry/v1/e5698820-4a80-11e4-9553-005056977bd0/1.0.0.0/-/license</licenseUrl>
<previewImageUrl>https://www.mathworks.com/responsive_image/160/120/0/0/0/cache/matlabcentral/mlc-downloads/downloads/submissions/25536/versions/1/screenshot.png</previewImageUrl>
<releaseNotesUrl>https://www.mathworks.com/add-ons/e5698820-4a80-11e4-9553-005056977bd0/d1034754-ea44-08f8-b658-22b590dfae7e/releaseNotes</releaseNotesUrl>
<installationFolder>Functions</installationFolder>
</addOn>
Matlab-Programs/resources/previewImage.png 0 → 100755
View file @ be6bdf92
Matlab-Programs/resources/previewImage.png

18.8 KiB

Matlab-Programs/resources/redblue.zip 0 → 100755
View file @ be6bdf92
File added
Matlab-Programs/resources/screenshot.png 0 → 100755
View file @ be6bdf92
Matlab-Programs/resources/screenshot.png

17 KiB