.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/PlotTest.m

+[t,y] = ode45(@cluster,[0:0.01:1],[1 2 3]);
+
+
+% figure(1)
+% plot(t,y(:,3));                 % plot of z(t) versus time
+% figure(2)
+% plot(t,y(:,1)); 
+% figure(3)
+% plot(y(:,1),y(:,3));            % plot of z versus x
+% figure(4)
+% plot3(y(:,1),y(:,2),y(:,3));    % 3D plot of trajectory
+% figure(5)
+% plot(y(:,1),y(:,2));            % plot of z versus x
+% figure(6)
+% plot(y(:,3),y(:,1));
+
+
+[x1,y1,z1] = meshgrid(-2:0.2:2,-2:0.2:2,-2:0.2:2);
+u = zeros(size(x1));
+v = zeros(size(y1));
+w = zeros(size(z1));
+
+t=0;
+for i = 1:numel(x1)
+    Yprime = cluster(t,[x1(i); y1(i); z1(i)]);
+    u(i) = Yprime(1);
+    v(i) = Yprime(2);
+    w(i) = Yprime(3);
+end
+
+
+quiver3(x1,y1,z1,u,v,w); figure(gcf)
\ No newline at end of file

Matlab-Programs/VanDePolOscillator.m

+
+
+
+mu = 100;
+
+
+F = @(t,y) [y(2); mu*(1-y(1)^2)*y(2)-y(1)];
+
+
+y0 = [0  1]';
+
+opts = odeset('stats','on')
+
+tspan= (0:1/36:1)*2*pi;
+
+tic
+[t,y] = ode15s(F,[0 460],y0,opts)    % stiff solver
+toc
+
+
+
+% plot(t,y(:,1), '.')
+% axis square
+% axis(1.2*[-1 1 -1 1])
+
+% plot phase-plane 
+plot(y(:,1),y(:,2),'.-')
\ No newline at end of file

Matlab-Programs/classifyMIN.m

+
+function [A, angle, type] = 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));
+
+
+
+% H = [q1 q3; q3 q2];         
+%check condition of H first
+% fprintf('condition number of Matrix H: %d \n', cond(H));
+
+
+
+
+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;
+    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
+        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
+        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;
+        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
+        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
+        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
+    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
+    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;
+    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? 
+k = sqrt(abs(a1) + abs(a2));  % ? 
+
+
+% 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/cluster.m

+function f = cluster(t,y)
+  %BD
+  a=1;
+  b=1.2;
+  %equilibrium values
+  %c1 equilbrium value
+    f = zeros(size(y));
+    f(1) = -50*a*y(1)-b*y(1)+15*a*y(1)*y(2)+20*a*y(2)*y(3)+y(2)*b+9*a*y(2)^2+6*a*y(1)^2-60*a*y(2)-80*a*y(3)+24*a*y(2)*y(3)+16*a*y(3)^2;
+    f(2) = 10*a*y(1) - a*y(1)*y(2) -4*a*y(1)*y(3) -2*a*y(1)^2 -b*y(2) +3*a*y(2)^2 -10*a*y(2) +4*a*y(2)*y(3) +b*y(3);
+    f(3) = -2*a*y(1)*y(2) - 3*a*y(2)^2 -4*a*y(2)*y(3) +10*a*y(2)-b*y(3);
+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/compute_r.m

+function 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);
+
+
+% TP = v*v';
+% L = stVenant(TP,mu,lambda);
+
+
+
+%compute Q(v_alpha x v_alpha)
+Q = q1.*v(1).^4 + q2.*v(2).^4 + 2.*q3.*v(1).^2.*v(2).^2;
+
+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 = trace(L'*B)./ Q;
+r = L./Q;
+
+
+end
+

Matlab-Programs/example1.m

+
+
+
+
+
+
+F = @(t,y) 2*y;
+
+t0 = 0;
+h = 1;
+tfinal = 3;
+y0 = 10;
+
+ode1(F,t0,h,tfinal,y0);
\ No newline at end of file

Matlab-Programs/exp.txt

+     x       exp(x)
+  0.00   1.00000000
+  0.10   1.10517092
+  0.20   1.22140276
+  0.30   1.34985881
+  0.40   1.49182470
+  0.50   1.64872127
+  0.60   1.82211880
+  0.70   2.01375271
+  0.80   2.22554093
+  0.90   2.45960311
+  1.00   2.71828183