%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Given boundary condition on domain (0,1)*(0,1),
% compute the minimal surface using div( Du/ ds ) = 0,
%     where ds = sqrt( 1 + (u_x)^2 + (u_y)^2 )
%
% We use the fnite difference method on a uniform grid.
% Details can be found in the PDF description file.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function minimalsurface
  % Set grid size to be (M+1)*(M+1)
  M = 20;
  h = 1.0/M;
  [x,y] = meshgrid(0:h:1,0:h:1);

  % Example 1
  V = reshape((  (x).^2+(y).^2   )',(M+1)*(M+1),1);
  [SOL,iter] = Newton(V,M,h);
  figure(1);
  surf(reshape(SOL,M+1,M+1)');
  title(['Minimal Surface Example 1: Newton iteration = ' num2str(iter)]);

  % Example 1
  V = reshape((  (x-.5).^2+(y-.5).^2   )',(M+1)*(M+1),1);
  [SOL,iter] = Newton(V,M,h);
  figure(2);
  surf(reshape(SOL,M+1,M+1)');
  title(['Minimal Surface Example 2: Newton iteration = ' num2str(iter)]);

  % Example 3
  V = reshape((  cos(pi*(x-.5))   )',(M+1)*(M+1),1);
  [SOL,iter] = Newton(V,M,h);
  figure(3);
  surf(reshape(SOL,M+1,M+1)');
  title(['Minimal Surface Example 3: Newton iteration = ' num2str(iter)]);

  % Example 4
  V = reshape((  cos(pi*(x-.5))-cos(pi*(y-.5))   )',(M+1)*(M+1),1);
  [SOL,iter] = Newton(V,M,h);
  figure(4);
  surf(reshape(SOL,M+1,M+1)');
  title(['Minimal Surface Example 4: Newton iteration = ' num2str(iter)]);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Newton's method for solving a
% nonlinear system of equations F(V) = 0.
% The iteration formula is
%  Vnew = V - marchingstepsize * (DF(V))^{-1} * F(V)
% where DF(V) is the Jacobian.
%
% The marchingstepsize = 2^(-m),
% where m >=0 is the smallest integer such that
%  norm (Vnew) < ( 1 - alpha * 2^(-m) ) * norm(V)
%  (this guarantees reducing of norm(F(V)) in the iteration)
% Here we set alpha = 1e-4.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [Vnew, iter] = Newton(V,M,h)
  % Maximum iteration is set to be 1000
  for iter=1:1000
    FV = F(V,M,h);
    DirectionVector = ComputeDirectionVector(V,FV,M,h);
    Vnew = V - DirectionVector;

    marchingstepsize = 1.0;
    while (norm(F(Vnew,M,h))>(1-(1e-4)*marchingstepsize)*norm(FV))
      marchingstepsize = .5 * marchingstepsize;
      Vnew = V - marchingstepsize * DirectionVector;
    end

    % Stoping criteria is || Vnew - V || < 1e-6
    if (norm(Vnew-V)<1e-6)
      break;
    end
    V = Vnew;
  end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Compute F(V)
% In actual computation, we reshape V
% to an (M+1)*(M+1) matrix for simplicity.
%
% v and fv are the reshaped matrix form
% for V and FV, respectively.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [FV] = F(V,M,h)
  v = reshape(V,M+1,M+1)';
  ds = Ds(v,M,h);
  fv = Dxf(Dxb(v,M,h)./ds,M,h) + Dyf(Dyb(v,M,h)./ds,M,h);
  % set F(V,M) = 0 on boundary nodes
  fv([1,end],:) = 0;
  fv(:,[1,end]) = 0;
  % reshape it back to a vector
  FV = reshape(fv',(M+1)*(M+1),1);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Assemble matrix DFV = DF(V) and
% compute DirectionVector = (DFV)^(-1) * FV
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [DirectionVector] = ComputeDirectionVector(V,FV,M,h)
  DFV = -speye((M+1)*(M+1),(M+1)*(M+1));

  % Set several intermediate variables
  vv = reshape(V,M+1,M+1)';
  ds = 1./(Ds(vv,M,h)).^3/h^2;
  dxdy = Dxb(vv,M,h).*Dyb(vv,M,h).*ds;
  dxdx = ds.*(1+Dxb(vv,M,h).^2)-dxdy;
  dydy = ds.*(1+Dyb(vv,M,h).^2)-dxdy;

  % Set the value of each row for DFV
  % Details can be found in the PDF description file
  for m=2:M
    for l=2:M
      row = (l-1)*(M+1) + m;
      col = row + [-M-1, -M, -1, 1, M, M+1];
      value = [dxdx(l,m), dxdy(l,m+1), dydy(l,m), ...
	     dydy(l,m+1), dxdy(l+1,m), dxdx(l+1,m) ];

      DFV(row,col) = value;
      DFV(row,row) = -sum(value);
    end
  end

  % Compute DirectionVector = (DFV)^{-1} * FV
  DirectionVector = DFV \ FV;


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Finally, there are several utility functions
% for computing various finite difference schemes.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Backward difference in x direction
function [dv] = Dxb(vv,M,h)
  dv = [zeros(M+1,1), diff(vv,1,2)/h];

% Forward difference in x direction
function [dv] = Dxf(vv,M,h)
  dv = [diff(vv,1,2)/h, zeros(M+1,1)];

% Backward difference in y direction
function [dv] = Dyb(vv,M,h)
  dv = [zeros(1,M+1); diff(vv,1,1)/h];

% Forward difference in y direction
function [dv] = Dyf(vv,M,h)
  dv = [diff(vv,1,1)/h; zeros(1,M+1)];

% unit surface area = sqrt( 1 + Dxb^2 + Dyb^2 )
function [dv] = Ds(vv,M,h)
  dv = sqrt(1 + (Dxb(vv,M,h)).^2 + (Dyb(vv,M,h)).^2);