function mit18086_levelset_front(method)
%MIT18086_LEVELSET_FRONT
%    Computes the movement of fronts under a given velocity field
%    using the level set method. The iteration scheme is first order
%    accurate, and a reinitialization iteration is used.
%    The front velocity can be chosen to be constant or proportional
%    to the negative curvature.

% 03/2008 by Benjamin Seibold
% Feel free to modify for teaching and learning.
%-----------------------------------------------------------------------
n = 90;       % number of space gridpoints
dt = 8e-3;    % time step
tf = 2e-0;    % final time
nr = 1;       % number of reinitialization steps
nsteps = 20;  % number of steps with graphic output
%-----------------------------------------------------------------------
if nargin<1, method = 1; end
nt = ceil(tf/dt); dt = tf/nt;
x = linspace(0,1,n)'; h = x(2)-x(1);
[X,Y] = meshgrid(x); ax = [min(x) max(x) min(x) max(x)];
%-----------------------------------------------------------------------
% initial conditions
P = sqrt((X-.45).^2+((Y-.8)*3).^2)-.4;         % place an ellipse
P = max(P,.07-sqrt((X-.6).^2+(Y-.78).^2));     % cut a hole out
P = min(P,max(max(.7-X,X-.87),max(.2-Y,Y-.6))); % add a rectangle
P = min(P,sqrt((X-.3).^2+(Y-.25).^2)-.195);    % add a circle
%-----------------------------------------------------------------------
for it = 1:nt
    switch method
        case 1, F = 5e-2;              % movement with constant velocity
        case 2, F = -curvature(P,h)*5e-3;     % movement under curvature
    end
    P = P-dt*FabsgradP(P,h,F);                        % level set update
    for ir = 1:nr                               % reinitialization steps
        P = P-dt*FabsgradP(P,h,P./sqrt(P.^2+(2*h)^2),1);
    end
    if it==1|floor(nsteps*it/nt)>floor(nsteps*(it-1)/nt) % visualization
        clf
        subplot(1,2,1), contourf(x,x,-P,[0 0],'k-')
        axis equal, axis(ax)
        title(sprintf('geometry at t=%0.2f',it*dt))
        subplot(1,2,2), surf(x,x,-P,'EdgeAlpha',.2), hold on
        patch([0 1 1 0],[0 0 1 1],[0 0 0 0],'k','FaceAlpha',.5)
        hold off, axp = [min(0,min(min(-P))) max(0,max(max(-P)))];
        axis([ax axp]), title('level set function')
        drawnow
    end
end

%=======================================================================

function dP = FabsgradP(P,h,F,c)
if nargin<4, c = 0; if nargin<3, F = 1; end, end
DxP = diff(P)/h;   DxmP = DxP([1 1:end],:); DxpP = DxP([1:end end],:);
DyP = diff(P')'/h; DymP = DyP(:,[1 1:end]); DypP = DyP(:,[1:end end]);
Np = sqrt(max(DxmP,0).^2+min(DxpP,0).^2+max(DymP,0).^2+min(DypP,0).^2);
Nm = sqrt(min(DxmP,0).^2+max(DxpP,0).^2+min(DymP,0).^2+max(DypP,0).^2);
dP = max(F,0).*(Np-c)+min(F,0).*(Nm-c);

function F = curvature(P,h)
% computes curvature by central differences
Pxx = diff(P([1 1:end end],:),2)/h^2;
Pyy = diff(P(:,[1 1:end end])',2)'/h^2;
Px = (P(3:end,:)-P(1:end-2,:))/(2*h); Px = Px([1 1:end end],:);
Py = (P(:,3:end)-P(:,1:end-2))/(2*h); Py = Py(:,[1 1:end end]);
Pxy = (Px(:,3:end)-Px(:,1:end-2))/(2*h); Pxy = Pxy(:,[1 1:end end]);
F = (Pxx.*Py.^2-2*Px.*Py.*Pxy+Pyy.*Px.^2)./(Px.^2+Py.^2).^1.5;
F = min(max(F,-1/h),1/h);