[p,wndx] = project2poly(p,find(~in&ok),edgexy);     

% Find the closest node in P for each FIXED node and replace P(i,:) with
% FIXED(i,:)
fix = zeros(size(fixed,1),1);
for k = 1:size(ndx,1)
   
   x = fixed(k,1);
   y = fixed(k,2);

   if isnan(ndx(k))
      % Slow search for all p(ok,:)
      [tmp,tmp] = min( (p(ok,1)-x).^2+(p(ok,2)-y).^2 );
      fix(k) = tmp;
   else
      % Search nodes in ndx(k)
      d = inf;
      j = 1;
      while j<=3
         cn = t(ndx(k),j);
         if ok(cn)
            dkj = (p(cn,1)-x)^2+(p(cn,2)-y)^2;
            if dkj<d
               fix(k) = cn;
               d = dkj;
            end
         end
         j = j+1;
      end
   end
   
   if fix(k)==0      
      % Slow search for all p(ok,:)
      [tmp,tmp] = min( (p(ok,1)-x).^2+(p(ok,2)-y).^2 );
      fix(k) = tmp;
   end

end
p(fix,:) = fixed;

% Take internal nodes
p = p(ok,:);

% Re-index to keeps lists consistent
wndx = wndx(ok);
j = zeros(length(ok),1);
j(ok) = 1;
j = cumsum(j);
t = j(t(tin,:));
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 [p,ndxnew] = project2poly(p,bnd,edgexy,ndx)

% Project the points in P(BND) onto the closest edge of the polygon defined
% by the edge segments in EDGEXY. NDX is an optional argument defining
% the edge to project onto: P(BND) is projected onto EDGEXY(NDX(BND)).

% Uses (something like?) a double sweep-line approach to reduce the number
% of edges that are required to be tested in order to determine the closest
% edge for each point. On average only size(EDGEXY)/4 comparisons need to
% be made for each point.

if nargin<4 || isempty(ndx)
   ndx = zeros(size(p,1),1);
end
ndxnew = zeros(size(p,1),1);
todo = true(size(bnd));

% Check NDX first
for k = 1:length(bnd)
   cn = bnd(k);
   if ndx(cn)>0
      j = ndx(cn);

      x1 = edgexy(j,1); x2mx1 = edgexy(j,3)-x1;
      y1 = edgexy(j,2); y2my1 = edgexy(j,4)-y1;

      r = ((p(cn,1)-x1)*x2mx1+(p(cn,2)-y1)*y2my1)/(x2mx1^2+y2my1^2);
      if (r>0) && (r<1)
         todo(k) = false;
         p(cn,1) = x1+r*x2mx1;
         p(cn,2) = y1+r*y2my1;
         ndxnew(cn) = j;
      end

   end
end

% Do a full search for points not already projected
if any(todo)

   bnd = bnd(todo);

   % Choose the direction with the biggest range as the "y-coordinate" for the
   % test. This should ensure that the sorting is done along the best
   % direction for long and skinny problems wrt either the x or y axes.
   dxy = max(p)-min(p);
   if dxy(1)>dxy(2)
      % Flip co-ords if x range is bigger
      p = p(:,[2,1]);
      edgexy  = edgexy(:,[2,1,4,3]);
      flip = true;
   else
      flip = false;
   end

   % Ensure edgexy(:,[1,2]) contains the lower y value
   swap = edgexy(:,4)<edgexy(:,2);
   edgexy(swap,:) = edgexy(swap,[3,4,1,2]);

   % Sort edges
   [ilower,ilower] = sort(edgexy(:,2));                                    % Sort edges by lower y value
   edgexy_lower = edgexy(ilower,:);
   [iupper,iupper] = sort(edgexy(:,4));                                    % Sort edges by upper y value
   edgexy_upper = edgexy(iupper,:);

   % Mean edge y value
   ne = size(edgexy,1);
   ymean = 0.5*( sum(sum(edgexy(:,[2,4]))) )/ne;

   % Loop through points
   for k = 1:length(bnd)

      cn = bnd(k);
      x = p(cn,1);
      y = p(cn,2);
      d = inf;

      if y<ymean

         % Loop through edges bottom up
         for j = 1:ne
            y2 = edgexy_lower(j,4);
            if y2>=(y-d)
               y1 = edgexy_lower(j,2);
               if y1<=(y+d)

                  % Calculate the distance along the normal projection from [x,y] to the jth edge
                  x1 = edgexy_lower(j,1); 
                  x2mx1 = edgexy_lower(j,3)-x1;
                  y2my1 = y2-y1;

                  r = ((x-x1)*x2mx1+(y-y1)*y2my1)/(x2mx1^2+y2my1^2);
                  if r>1                                                   % Limit to wall endpoints
                     r = 1;
                  elseif r<0
                     r = 0;
                  end
                  xn = x1+r*x2mx1;
                  yn = y1+r*y2my1;

                  dj = (xn-x)^2+(yn-y)^2;
                  if ( dj<d^2 )
                     d = sqrt(dj);
                     p(cn,1) = xn;
                     p(cn,2) = yn;
                     ndxnew(cn) = ilower(j);
                  end

               else
                  break
               end
            end
         end

      else

         % Loop through edges top down
         for j = ne:-1:1
            y1 = edgexy_upper(j,2);
            if y1<=(y+d)
               y2 = edgexy_upper(j,4);
               if y2>=(y-d)

                  % Calculate the distance along the normal projection from [x,y] to the jth edge
                  x1 = edgexy_upper(j,1); 
                  x2mx1 = edgexy_upper(j,3)-x1;
                  y2my1 = y2-y1;

                  r = ((x-x1)*x2mx1+(y-y1)*y2my1)/(x2mx1^2+y2my1^2);
                  if r>1                                                   % Limit to wall endpoints
                     r = 1;
                  elseif r<0
                     r = 0;
                  end
                  xn = x1+r*x2mx1;
                  yn = y1+r*y2my1;

                  dj = (xn-x)^2+(yn-y)^2;
                  if ( dj<d^2 )
                     d = sqrt(dj);
                     p(cn,1) = xn;
                     p(cn,2) = yn;
                     ndxnew(cn) = iupper(j);
                  end

               else
                  break
               end
            end
         end

      end

   end

   if flip
      p = p(:,[2,1]);
   end

end


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
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