.gitignore

+# ignore all build folders
+/build*/
+build/
+# ignore backup files
+*~
+# ignore python files
+*.pyc
+
+

.kdev4/dune-microstructure.kdev4

+[Buildset]
+BuildItems=@Variant(\x00\x00\x00\t\x00\x00\x00\x00\x01\x00\x00\x00\x0b\x00\x00\x00\x00\x01\x00\x00\x00&\x00d\x00u\x00n\x00e\x00-\x00m\x00i\x00c\x00r\x00o\x00s\x00t\x00r\x00u\x00c\x00t\x00u\x00r\x00e)
+
+[CMake]
+Build Directory Count=1
+Current Build Directory Index-Host System=0
+
+[CMake][CMake Build Directory 0]
+Build Directory Path=/home/klaus/Desktop/DUNE2.8/dune-microstructure/build-cmake
+Build Type=\n
+CMake Binary=/usr/bin/cmake
+CMake Executable=/usr/bin/cmake
+Environment Profile=
+Extra Arguments=
+Install Directory=
+Runtime=Host System
+
+[Launch]
+Launch Configurations=Launch Configuration 0
+
+[Launch][Launch Configuration 0]
+Configured Launch Modes=execute
+Configured Launchers=nativeAppLauncher
+Name=dune-microstructure
+Type=Native Application
+
+[Launch][Launch Configuration 0][Data]
+Arguments=
+Dependencies=@Variant(\x00\x00\x00\t\x00\x00\x00\x00\x01\x00\x00\x00\x0b\x00\x00\x00\x00\x03\x00\x00\x00&\x00d\x00u\x00n\x00e\x00-\x00m\x00i\x00c\x00r\x00o\x00s\x00t\x00r\x00u\x00c\x00t\x00u\x00r\x00e\x00\x00\x00\x06\x00s\x00r\x00c\x00\x00\x00&\x00d\x00u\x00n\x00e\x00-\x00m\x00i\x00c\x00r\x00o\x00s\x00t\x00r\x00u\x00c\x00t\x00u\x00r\x00e)
+Dependency Action=Build
+EnvironmentGroup=
+Executable=file:///home/klaus/Desktop/DUNE2.8/dune-microstructure
+External Terminal=konsole --noclose --workdir %workdir -e %exe
+Project Target=dune-microstructure,src,dune-microstructure
+Use External Terminal=false
+Working Directory=file:///home/klaus/Desktop/DUNE2.8/dune-microstructure/build-cmake/src
+isExecutable=false
+
+[Project]
+VersionControlSupport=kdevgit

CMakeLists.txt

+# 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/BMatrix.txt

+1 1 -0.375507808365302675
+1 2 -0.59999999999997744
+1 3 6.94558662207737461e-10

Matlab-Programs/Minimization_Script.m

+clear all
+clc
+
+% --- 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(q3 >= q1)
+assume(q2 >= q3)
+
+
+v = [v1; v2];
+
+%should be sqrt(2) instead of 2!
+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
+% f_plus = subs(f_plus,b3,0)  % set b3
+% f_plus = subs(f_plus,q3,q1)
+
+% f_plus = subs(f_plus,q1,40);
+% f_plus = subs(f_plus,q3,63.9);
+% f_plus = subs(f_plus,q2,34.9);
+% f_plus = subs(f_plus,b1,4);
+% f_plus = subs(f_plus,b2,2.4);
+% f_plus = subs(f_plus,b3,2.4);
+
+% f_plus = subs(f_plus,q1,40);
+% f_plus = subs(f_plus,q3,20);
+% f_plus = subs(f_plus,q2,50);
+% f_plus = subs(f_plus,b1,5);
+% f_plus = subs(f_plus,b2,6);
+% f_plus = subs(f_plus,b3,7);
+
+f_plus = subs(f_plus,q1,40);
+f_plus = subs(f_plus,q3,63.9);
+f_plus = subs(f_plus,q2,70);
+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);
+
+
+f_minus = subs(f_minus,q1,40);
+f_minus = subs(f_minus,q3,63.9);
+f_minus = subs(f_minus,q2,70);
+f_minus = subs(f_minus,b1,rndm1);
+f_minus = subs(f_minus,b2,rndm2);
+f_minus = subs(f_minus,b3,rndm3);
+
+
+% Compute Gradient 
+df_plusx = diff(f_plus,v1);
+df_plusy = diff(f_plus,v2);
+
+df_minusx = diff(f_minus,v1);
+df_minusy = diff(f_minus,v2);
+
+
+eq1 = df_plusx == 0;
+eq2 = df_plusy == 0;
+eqns = [eq1, eq2]
+
+eq3 = df_minusx == 0;
+eq4 = df_minusy == 0;
+eqns_minus = [eq3, eq4]
+
+
+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);
+% 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)
+% S_minus = solve(eqns_minus,v1,v2);
+% S_minus.v1
+% S_minus.v2
+
+%---------- TEST ---------------------
+% r = subs(subs(df_plusx,v1,A(3)),v2,B(3));
+fprintf('Testing equation grad(f) = 0  with stationary points')
+double(subs(subs(df_plusx,v1,A(1)),v2,B(1)))
+double(subs(subs(df_plusx,v1,A(2)),v2,B(2)))
+double(subs(subs(df_plusx,v1,A(3)),v2,B(3)))
+% ------------------------------------
+
+
+
+
+
+fprintf('print stationary points of f_plus:')
+double(A)
+double(B)
+fprintf('print stationary points of f_minus:')
+double(A_minus)
+double(B_minus)
+
+
+% determine global Minimizer from stationary points:
+fprintf('function values at stationary points:')
+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')
+
+
+% determine global Minimizer from stationary points:
+fprintf('function values at stationary points:')
+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')
+
+globalMinimizerValue = min(Min_plus,Min_minus)
+
+% Plot function
+fsurf(@(x,y) f_plus(x,y,q1,q2,q3,b1,b2,b3))
+hold on 
+plot3(A,B,T,'g*')
+% 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,'g*')
+
+
+
+
+
+%Write to txt-File
+fileID = fopen('txt.txt','w');
+fprintf(fileID,'%s' , latex(S.v1));
+fclose(fileID);
+
+
+
+
+
+
+
+
+

Matlab-Programs/PhaseDiagramsScript.m

+
+
+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/QMatrix.txt

+1 1 1.90510815738902228
+1 2 2.60132454317226442e-26
+1 3 -4.93750948541723109e-11
+2 1 -5.99537621698159769e-27
+2 2 2.08333333333341164
+2 3 -2.82386576263450663e-18
+3 1 -2.16574525495542866e-10
+3 2 -2.82386567902393238e-18
+3 3 1.92269764555762257

Matlab-Programs/TestCompute.m

+
+
+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

+
+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
+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);
+
+
+
+
+% compute angle between [sqrt(a1) , sqrt(a2)] and e1:
+% angle = atan2(sqrt(a2),sqrt(a1));
+if (type == 3 )
+   angle = atan2(a2,a1);
+else
+   angle = 0;
+end
+% angle = atan2(norm(cross(a,b)), dot(a,b))
+
+
+%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

+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/redblue.m

+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

+<?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

+<?xml version="1.0" encoding="UTF-8"?>
+<paths>
+  <path>.</path>
+</paths>

Matlab-Programs/resources/metadata.xml

+<?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>