%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% curvature.m - level set framework for motion of a curve via curvature
%
%   Benjamin Ong, Simon Fraser University
%   May 5, 2002
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

clear;
close all;

% Geometry
disp('Initializing Mesh');
hmesh = 0.1;
fem.geom=geomcsg({rect2(0,1,0,1)});
fem.mesh=meshinit(fem,'hmax',hmesh);
%fem.mesh=meshrefine(fem);
%fem.mesh=meshrefine(fem);
figure(1);
meshplot(fem);
axis([0 1 0 1]);
hold on;
fem.mesh = meshinit(fem,'jiggle','on');
meshplot(fem,'edgecolor',[1 0 0],'boundcolor',[0 1 0]);
fem.mesh = meshinit(fem,'jiggle','on');
meshplot(fem,'edgecolor',[1 1 0],'boundcolor',[1 0 1]);
pause(0);

fem.mesh=meshrefine(fem);

% initializing phi
disp('Initializing Initial Contour');
x = fem.mesh.p(1,:);
y = fem.mesh.p(2,:);

[m,n] = size(x);
phi = zeros(m,n);

%    p1            p2
%    ----------------
%    |   p7          | p5
%    |   ------------
%    |  |           
%    |  |p8           p6
%    |  -------------     
%    |               |
%    -----------------
%    p4               p3

p1 = [0.2 0.8]; p2 = [0.8 0.8]; p3 = [0.8 0.2];
p4 = [0.2 0.2]; p5 = [0.8,0.6]; p6 = [0.8,0.4];
%p7 = [0.25,0.6]; p8 = [0.25,0.4];
p7 = [0.5,0.6]; p8 = [0.5,0.4];
%p7 = [0.8 0.6]; p8 = [0.8 0.4];

temp = - sqrt((x-p1(1)).^2+(y-p1(2)).^2);
phi = phi + (x<=p1(1)).*(y>=p1(2)).*temp;
temp = - sqrt((x-p2(1)).^2+(y-p2(2)).^2);
phi = phi + (x>=p2(1)).*(y>=p2(2)).*temp;
temp = - sqrt((x-p3(1)).^2+(y-p3(2)).^2);
phi = phi + (x>=p3(1)).*(y<=p3(2)).*temp;
temp = - sqrt((x-p4(1)).^2+(y-p4(2)).^2);
phi = phi + (x<=p4(1)).*(y<=p4(2)).*temp;
temp = x - p4(1);
phi = phi + (x<=p4(1)).*(y>=p4(2)).*(y<=p1(2)).*temp;
temp = p1(2) - y;
phi = phi + (y>=p1(2)).*(x>=p1(1)).*(x<=p2(1)).*temp;
temp = y - p4(2);
phi = phi + (y<=p4(2)).*(x>=p4(1)).*(x<=p3(1)).*temp;
temp = p2(1) - x;
phi = phi + (x>=p2(1)).*(y<=p2(2)).*(y>=p5(2)).*temp;
temp = p6(1) - x;
phi = phi + (x>=p6(1)).*(y<=p6(2)).*(y>=p3(2)).*temp;
temp1 = -sqrt((x-p5(1)).^2+(y-p5(2)).^2);
temp2 = -sqrt((x-p6(1)).^2+(y-p6(2)).^2);
temp = max(temp1,temp2);
phi = phi + (x>=p6(1)).*(y<=p5(2)).*(y>=p6(2)).*temp;
temp1 = y - p5(2);
temp2 = p6(2) - y;
temp3 = p7(1) - x;
temp4 = max(temp1,temp2);
temp = max(temp3,temp4);
phi = phi + (x>=p7(1)).*(x<=p5(1)).*(y<=p7(2)).*(y>=p8(2)).*temp;
temp1 = p1(2)-y;
temp2 = y - p7(2);
temp3 = x - p1(1);
temp4 = p5(1)-x;
temp5 = min(temp1,temp2);
temp6 = min(temp3,temp4);
temp = min(temp5,temp6);
phi = phi + (x>=p7(1)).*(x<=p5(1)).*(y>=p7(2)).*(y<=p1(2)).*temp;
temp1 = p8(2) - y;
temp2 = y - p4(2);
temp3 = x - p4(1);
temp4 = p6(1) - x;
temp5 = min(temp1,temp2);
temp6 = min(temp3,temp4);
temp = min(temp5,temp6);
phi = phi + (x>=p8(1)).*(x<=p6(1)).*(y>=p4(2)).*(y<=p8(2)).*temp;
temp1 = x - p1(1);
temp2 = p7(1) - x;
temp3 = p1(2) - y;
temp4 = y - p4(2);
temp5 = min(temp1,temp2);
temp6 = min(temp3,temp4);
temp = min(temp5,temp6);
phi = phi + (x>=p1(1)).*(x<=p7(1)).*(y<=p7(2)).*(y>=p8(2)).*temp;
temp1 = p1(2)-y;
temp2 = x-p1(1);
temp3 = sqrt((x-p7(1)).^2+(y-p7(2)).^2);
temp4 = min(temp1,temp2);
temp = min(temp3,temp4);
phi = phi + (x>=p1(1)).*(x<=p7(1)).*(y<=p1(2)).*(y>=p7(2)).*temp;
temp1 = x - p4(1);
temp2 = y - p4(2);
temp3 = sqrt((x-p8(1)).^2+(y-p8(2)).^2);
temp4 = min(temp1,temp2);
temp = min(temp3,temp4);
phi = phi + (x>=p4(1)).*(x<=p8(1)).*(y<=p8(2)).*(y>=p4(2)).*temp;

clear temp p1 p2 p3 p4 p5 p6 p7 p8
clear temp temp1 temp2 temp3 temp4 temp5 temp6

% initial contour
fem.init = phi';

% constants
fem.variables.epsilon = 1e-8;

% PDE coefficients
disp('Setting PDE Coefficients');
clear equ
equ.da={{{'1./((ux.^2+uy.^2).^0.5+1e-8)'}}};
equ.c={{{'1./((ux.^2+uy.^2).^0.5+1e-8)'}}};
equ.al={{{'0';'0'}}};
equ.ga={{{'0';'0'}}};
equ.be={{{'0';'0'}}};
equ.a={{{'0'}}};
equ.f={{{'0'}}};
equ.var={'absux','sqrt(ux.^2+uy.^2)'};
equ.ind=1;
fem.equ = equ;

% Space dimensions
fem.dim={'u'};
fem.sdim={'x','y'};

% Boundary Conditions
fem.border = 0;

% Usage
fem.usage = 1;

% Problem form
fem.form='coefficient';

% Differentiation
fem.diff='off';

% Differentiation simplification
fem.simplify='on';

%Boundary Conditions
clear bnd
bnd.q={{{'0'}}};
bnd.g={{{'0'}}};
bnd.h={{{'1'}}};
bnd.r={{{'0'}}};
bnd.ind=ones(1,4);
bnd.type={'neu'};
fem.bnd=bnd;

% Solution
disp('Solving problem');

fem.sol=femtime(fem,...
	'tlist',  0:1e-3:1e-2,...
	'atol',   1e-8,...
	'rtol',   1e-6,...
	'jacobian','equ',...
	'mass',   'full',...
	'ode',    'ode15s',...
	'odeopt', struct('InitialStep',{[]},'MaxOrder',{5},'MaxStep',{[]}),...
	'out',    'sol',...
	'stop',   'on',...
	'report', 'on',...
	'context','local',...
	'sd',     'off',...
	'Epoint', 'gauss2',...
	'Nullfun','flnullorth',...
	'Tpoint', 'gauss2',...
	'Solcomp',1);

figure(2);

disp('displaying solution');

postmovie(fem,...
	'context','local',...
	'contdata',{'u','cont','on'},...
	'contlevels',[-1 0 1],...
	'contstyle','color',...
	'contlabel','off',...
	'contmaxmin','off',...
	'contbar','on',...
	'contmap','cool',...
	'geom',   'on',...
	'geomcol','bginv',...
	'movie',  'on',...
	'repeat', 1,...
	'fps',    1,...
	'title',  'Contour: u (u)  ',...
	'renderer','zbuffer')

figure(3);

        
% Plot solution
postplot(fem,...
	'context','local',...
	'contdata',{'u','cont','on'},...
	 'contlevels',[0 0],...
	 'contstyle','color',...
	 'contlabel','off',...
	 'contmaxmin','off',...
	 'contbar','off',...
	 'contmap','gray',...
	'geom',   'on',...
	'geomcol','bginv',...
	'view',   [0 90],...
	'title',  'Time=0',...
	'renderer','zbuffer',...
	'solnum', 1,...
	'axisvisible','off')