function [p,t,stats] = mesh2d(node,edge,hdata,options)

%  MESH2D: 2D unstructured triangular mesh generation.
%
% A 2D unstructured triangular mesh is generated based on a piecewise-
% linear geometry input. An iterative method is implemented to optimise 
% mesh quality. General multiply connected domains and element size 
% functions can be specified.
%
% Returns the final coordinates of the nodes P, and their triangulation T
% (with a counter-clockwise node ordering).
%
%  SHORT SYNTAX:
%
%  [p,t] = mesh2d(node);
%
% NODE defines the geometry nodes as an Nx2 array:
%
%  node  = [x1 y1; x2 y2; etc], geometry nodes specified in consectutive
%                               order, such that NODE(2,:) is joined with
%                               NODE(1,:) etc.
%
% An element size function is automatically generated based on the 
% complexity of the geometry. Generally this produces meshes with the 
% fewest number of triangles.
%
%  LONG SYNTAX:
%
%  [p,t] = mesh2d(node,edge,hdata,options);
%
% Blank arguments can be passed using the empty placeholder "[]".
%
% EDGE defines the connectivity between the points in NODE as a list of
% edges:
%
%   edge = [n1 n2; n2 n3; etc], connectivity between nodes to form
%                               geometry edges. If EDGE is specified it is
%                               not required that NODE be consectutive.
%
% HDATA is a structure containing user defined element size information. 
% HDATA can include the following fields:
%
%  hdata.hmax  = h0;                   Max allowable global element size.
%  hdata.edgeh = [e1,h1; e2,h2; etc];  Element size on specified geometry 
%                                      edges.
%  hdata.fun   = 'fun' or @fun;        User defined size function.
%  hdata.args  = {arg1, arg2, etc};    Additional arguments for HDATA.FUN.
%
% Calls to user specified functions must accept vectorised input of the 
% form H = FUN(X,Y,ARGS{:}), where X,Y are the xy coordinates where the
% element size will be evaluated and ARGS are optional additional arguments 
% as passed by HDATA.ARGS.
%
% An automatic size function is always generated to ensure that the
% geometry is adequately resolved. The overall size function is the minimum
% of the user specified and automatic functions.
%
% OPTIONS is a structure array that allows some of the "tuning" parameters
% used in the solver to be modified:
%
%   options.mlim   : The convergence tolerance. The maximum percentage 
%                    change in edge length per iteration must be less than 
%                    MLIM { 0.05, 5% }. 
%   options.maxit  : The maximum allowable number of iterations { 20 }.
%   options.dhmax  : The maximum allowable (relative) gradient in the size 
%                    function { 0.3 }.
%   options.output : Displays the mesh and the mesh statistics upon
%                    completion { TRUE }.
%
% EXAMPLE:
%
%   meshdemo                  % Will run the standard demos
%   mesh_collection(n)        % Will run some additional demos
%
% See also REFINE, SMOOTHMESH, DELAUNAYN

% STATS is an undocumented output used in debugging. Returns the algorithm 
% statistics usually printed to screen as a structure.

% Mesh2d is a delaunay based algorithm with a "Laplacian-like" smoothing
% operation built into the mesh generation process. 
% 
% An unbalanced quadtree decomposition is used to evaluate the element size 
% distribution required to resolve the geometry. The quadtree is 
% triangulated and used as a backgorund mesh to store the element size 
% data.  
%
% The main method attempts to optimise the node location and mesh topology 
% through an iterative process. In each step a constrained delaunay 
% triangulation is generated with a series of "Laplacian-like" smoothing 
% operations used to improve triangle quality. Nodes are added or removed 
% from the mesh to ensure the required element size distribution is 
% approximated.  
%
% The optimisation process generally returns well shaped meshes with no
% small angles and smooth element size variations. Mesh2d shares some 
% similarities with the Distmesh code: 
%
%   [1] P.-O. Persson, G. Strang, A Simple Mesh Generator in MATLAB.
%       SIAM Review, Volume 46 (2), pp. 329-345, June 2004
%
%   Darren Engwirda : 2005-07
%   Email           : d_engwirda@hotmail.com
%   Last updated    : 20/07/2007 with MATLAB 7.0
%
% Mesh2d is Copyright (C) 2007 Darren Engwirda. See "copyright.m" for full 
% details.
%
% Please email me any un-meshable geometries, meshing benchmarks or
% suggestions!

%hdata = [];

debug    = false;                                                         
t_tri    = 0;
t_inpoly = 0;
t_edge   = 0;
t_sparse = 0;
t_search = 0;
t_smooth = 0;
n_tri    = 0;

ts = cputime;
wbar = waitbar(0,'Checking geometry');

% Error checks
if nargin<4
   options = [];
   if nargin<3
      hdata = [];
      if nargin<2
         edge = [];
         if nargin<1
            error('Wrong number of inputs');
         end
      end
   end
elseif nargin>4
   error('Wrong number of inputs');
end
if nargout>3
   error('Wrong number of outputs');
end

% Get user options
[mlim,maxit,dhmax,output] = getoptions(options);

% Check geometry
[node,edge,hdata] = checkgeometry(node,edge,hdata);                        % Error checking for geometry
edgexy = [node(edge(:,1),:), node(edge(:,2),:)];                           % Geometry as edge endpoints [x1,y1,x2,y2]

% QUADTREE DECOMPOSITION
%  PH    : Background mesh nodes
%  TH    : Background mesh triangles
%  HH    : Size function value at PH
tic

[ph,th,hh] = quadtree(edgexy,hdata,dhmax,wbar); 
t_quad = toc;

waitbar(0,wbar,'Initialising mesh');

% INITIALISE MESH
%  P     : Initial nodes
%  T     : Initial triangulation
%  WNDX  : Closest edge for each node as indices into EDGEXY
%  TNDX  : Enclosing triangle for each node as indices into TH
%  FIX   : Indices of FIXED nodes in P
tic
[p,t,wndx,tndx,fix,fixed] = initmesh(ph,th,node,edge,edgexy); 
t_init = toc;

% MAIN LOOP
retri = false;
subit = 2;
for iter = 1:maxit

   tic
   if retri
      t = MyDelaunayn(p);                                                  % Delaunay triangulation via QHULL
      n_tri = n_tri+1;
   end
   t_tri = t_tri+toc;
   
   tic
   pc = (p(t(:,1),:)+p(t(:,2),:)+p(t(:,3),:))/3;                           % Centroids
   t = t(inpoly(pc,node,edge),:);                                          % Take triangles with internal centroids
   t_inpoly = t_inpoly+toc;
   
   tic
   [e,bnd] = getedges(t,size(p,1));                                        % Unique edges and boundary nodes
   bnd(fix) = false;                                                       % Don't flag fixed nodes
   bnd = find(bnd);
   t_edge = t_edge+toc;
   
   tic
   nume = size(e,1);
   S = sparse(e(:),[1:nume,1:nume],1,size(p,1),nume);                      % Sparse node-to-edge connectivity matrix
   N = max(sum(S,2),eps);                                                  % Num edge attachments per node
   t_sparse = t_sparse+toc;
   
   tic
   tndx = mytsearch(ph(:,1),ph(:,2),th,p(:,1),p(:,2),tndx);                % Find enclosing triangle in background mesh for nodes
   h = tinterp(ph,th,hh,p,tndx);                                           % Size function at nodes via linear interpolation
   h = 0.5*(h(e(:,1))+h(e(:,2)));                                          % from the background mesh. Average to edge midpoints.
   t_search = t_search+toc;

   % Inner smoothing iterations
   tic
   L = max(sqrt(sum((p(e(:,1),:)-p(e(:,2),:)).^2,2)),eps);                 % Edge length
   done = false;
   subit = max(subit,iter);                                                % Increment sub-iters with outer iters to aid convergence
   for subiter = 1:subit   
      
      r = L./h;                                                            % Ratio of actual to desired edge length
      
      pm = 0.5*[r,r].*(p(e(:,1),:)+p(e(:,2),:));                           % Edge midpoints, size function weighted

      W = max(S*r,eps);                                                    % Size function weighting
      p = (S*pm)./[W,W];                                                   % Weighted Laplacian-like smoothing
      p(fix,:) = fixed;                                                    % Don't update fixed nodes
      
      [p,wndx] = project2poly(p,bnd,edgexy,wndx,2*(S*L)./N);               % Project BND nodes onto the closest geometry edge
      
      Lnew = max(sqrt(sum((p(e(:,1),:)-p(e(:,2),:)).^2,2)),eps);           % Edge length
      move = norm((Lnew-L)./L,'inf');                                      % Percentage change
      L = Lnew;

      if move<mlim                                                         % Test convergence
         done = true;
         break
      end
      
   end
   t_smooth = t_smooth+toc;
   
   msg = ['Generating mesh (Iteration ',num2str(iter),')'];                % Show progress
   waitbar(mlim/max(move,eps),wbar,msg);

   r = L./h;
   if done && max(r)<3                                                     % Main loop convergence
      break
   end

   % Nodal density control
   retri = false;
   if iter<maxit
      test = find(r<=2/3);                                                 % Remove both nodes for edges with R<2/3
      if ~isempty(test) || any(N<2)
         prob = false(size(p,1),1);                                        % True for nodes to be removed
         prob(e(test,:)) = true;                                           % Edges with R<0.5
         prob(N<2) = true;                                                 % Remove nodes with less than 2 edge connections
         prob(fix) = false;                                                % Don't remove fixed nodes
         pnew = p(~prob,:);                                                % New node list

         tmp_wndx = wndx(~prob);
         tmp_tndx = tndx(~prob);
         
         j = zeros(size(p,1),1);                                           % Re-index FIX to keep consistent
         j(~prob) = 1;
         j = cumsum(j);
         fix = j(fix);
         
         retri = true;
      else
         pnew = p;                                                  
         tmp_wndx = wndx;
         tmp_tndx = tndx;
      end
      test = find(r>=2);                                                   % Add node at midpoint for edges with R>=2
      if ~isempty(test)
         p = [pnew; 0.5*(p(e(test,1),:)+p(e(test,2),:))];
         wndx = [tmp_wndx; zeros(length(test),1)];
         tndx = [tmp_tndx; zeros(length(test),1)];
         retri = true;
      else
         p = pnew;                                             
         wndx = tmp_wndx;
         tndx = tmp_tndx;
      end
   end

end
close(wbar)

if (iter==maxit) && ( ~done || (max(r)>3) )
   disp('WARNING: Maximum number of iterations reached. Solution did not converge!')
   disp('Please email the geometry and settings to d_engwirda@hotmail.com')
end

% Ensure consistent, CCW ordered triangulation
[p,t] = fixmesh(p,t);                                                      

% Element quality
q = quality(p,t);

% Method statistics
stats = struct('Iterations',iter,'Time',cputime-ts,'Triangles',size(t,1), ...
                  'Nodes',size(p,1),'Mean_quality',mean(q),'Min_quality',min(q));

if output
   figure('Name','Mesh')
   patch('faces',t,'vertices',p,'facecolor','None','edgecolor','b')
   hold on
   patch('faces',edge,'vertices',node,'facecolor','None','edgecolor','r')
   axis equal off; hold off
   disp(stats);
end
if debug
   disp(struct('t_quad',t_quad,'t_init',t_init,'t_tri',t_tri,'t_inpoly', ...
                  t_inpoly,'t_edge',t_edge,'t_sparse',t_sparse,'t_search', ...
                     t_search,'t_smooth',t_smooth,'n_tri',n_tri,'mean_size', ...
                        mean(r),'min_size',min(r),'max_size',max(r)));
end


%% SUB-FUNCTIONS
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [p,t,wndx,tndx,fix,fixed] = initmesh(p,t,node,edge,edgexy)

% Initialise the nodes, triangulation and data structures for the mesh 
% based on the quadtree data and geometry.

% Impose the geometry constraints on the quadtree mesh (approximately!)
pc = ( p(t(:,1),:)+p(t(:,2),:)+p(t(:,3),:) )/3;
ti = inpoly(pc,node,edge);
ok = false(size(p,1),1);
ok(t(ti,:)) = true;

% At this stage some nodes have been accepted that are not inside the 
% geometry. This is done because the quadtree triangulation will overlap 
% the geometry in some cases, so nodes invloved in the overlap are accepted 
% to get a reasonable distribution near the edges.

% Get an upper bound on the distance between the polygon and the points to
% be projected. The points cannot lie much more than an edge length from
% the polygon based on the way the quadtree works.
e = [t(:,[1,2]); t(:,[2,3]); t(:,[3,1])];
len = sum((p(e(:,2),:)-p(e(:,1),:)).^2,2);
lim = zeros(size(p,1),1);
for k = 1:size(e,1)
   lenk = len(k);
   if lenk>lim(e(k,1)), lim(e(k,1)) = lenk; end
   if lenk>lim(e(k,2)), lim(e(k,2)) = lenk; end
end
lim = sqrt(lim);

% Find overlapping nodes
[e,bnd] = getedges(t(ti,:),size(p,1));
bnd = find(bnd);

% Project overlapping points onto the polygon
[p,wndx] = project2poly(p,bnd,edgexy,[],2*lim);

% Set the geometry nodes as "fixed" throughout the iteration. Replace the
% closest nodes in the quadtree mesh with the fixed nodes.
fixed = node;
fix = dsearch(p(:,1),p(:,2),t,fixed(:,1),fixed(:,2));

% Ensure FIX is unique
[i,i] = unique(fix);
j = true(size(fix,1),1);
j(i) = false;

% For any non-unique FIX do a slow search to find the closest unique node
% in the quadtree mesh
if any(j)
   
   kk = true(size(p,1),1);                % Nodes ok to select from
   kk(fix) = false;
   
   jj = find(j);
   for k = 1:length(jj)
      [ii,ii] = min( (p(kk,1)-fixed(jj(k),1)).^2+(p(kk,2)-fixed(jj(k),2)).^2 );
      ikk = find(kk);
      fix(jj(k)) = ikk(ii);
      kk(ikk(ii)) = false;
   end
   
end
p(fix,:) = fixed;

% New mesh
ok(fix) = true;
p = p(ok,:);
t = t(ti,:);

% Re-index to keep indexing consistent
wndx = wndx(ok);
j = zeros(length(ok),1);
j(ok) = 1;
j = cumsum(j);
t = j(t);
fix = j(fix);
tndx = zeros(size(p,1),1);


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [e,bnd] = getedges(t,n)

% Get the unique edges and boundary nodes in a triangulation.

e = [t(:,[1,2]); t(:,[1,3]); t(:,[2,3])];                                  % Non-unique edges

swap = e(:,2)<e(:,1);                                                      % Ensure e(:,1) contains the lower value
e(swap,:) = e(swap,[2,1]);

e = sortrows(e);
idx = all(diff(e,1)==0,2);                                                 % Find shared edges
idx = [idx;false]|[false;idx];                                             % True for all shared edges
bnde = e(~idx,:);                                                          % Boundary edges
e = e(idx,:);                                                              % Internal edges
e = [bnde; e(1:2:end-1,:)];                                                % Unique edges and bnd edges for tri's

bnd = false(n,1);                                                          % True for boundary nodes
bnd(bnde) = true;


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [mlim,maxit,dhmax,output] = getoptions(options)

% Extract the user defined options

if ~isempty(options)
   if ~isstruct(options)
      error('OPTIONS must be a structure array');
   end
   if numel(options)~=1
      error('Options cannot be an array of structures');
   end
   fields = fieldnames(options);
   names = {'mlim','maxit','dhmax','output'};
   for k = 1:length(fields)
      if strcmp(fields{k},names)
         error('Invalid field in OPTIONS');
      end
   end
   if isfield(options,'mlim')                                              % Movement tolerance
      mlim = checkposscalar(options.mlim,'options.mlim');
   else
      mlim = 0.05;
   end
   if isfield(options,'maxit')                                             % Maximum iterations
      maxit = round(checkposscalar(options.maxit,'options.maxit'));
   else
      maxit = 20;
   end
   if isfield(options,'dhmax')                                             % Size function gradient limit
      dhmax = checkposscalar(options.dhmax,'options.dhmax');
   else
      dhmax = 0.3;
   end
   if isfield(options,'output')                                            % Output on/off
      output = checklogicalscalar(options.output,'options.output');
   else
      output = true;
   end
else                                                                       % Default values
   mlim = 0.05;
   maxit = 20;
   dhmax = 0.3;
   output = true;
end


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function var = checkposscalar(var,name)

% Helper function to check if var is a positive scalar.

if var<0 || any(size(var)>1)
   error([name,' must be a positive scalar']);
end


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function var = checklogicalscalar(var,name)

% Helper function to check if var is a logical scalar.

if ~islogical(var) || any(size(var)>1)
   error([name,' must be a logical scalar']);
end
