function [p,t] = meshpoly(node,edge,qtree,p,options)
%*****************************************************************************80
%
%% meshpoly(): Core meshing routine called by mesh2d and meshfaces.
%
% Discussion:
%
% 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:
%
% Do not call this routine directly, use mesh2d or meshfaces instead!
%
% Licensing:
%
% This code is distributed under the GNU LGPL license.
%
% Modified:
%
% 13 March 2013
%
% Author:
%
% Darren Engwirda
%
% Reference:
%
% [1] P.-O. Persson, G. Strang, A Simple Mesh Generator in MATLAB.
% SIAM Review, Volume 46 (2), pp. 329-345, June 2004
%
% Input:
%
% NODE : Nx2 array of geometry XY co-ordinates
% EDGE : Mx2 array of connections between NODE, defining geometry
% edges
% QTREE : Quadtree data structure, defining background mesh and element
% size function
% P : Qx2 array of potential boundary nodes
% OPTIONS : Meshing options data structure
% WBAR : Handle to progress bar
%
% Output:
%
% P : Nx2 array of triangle nodes
% T : Mx3 array of triangles as indices into P
%
shortedge = 0.75;
longedge = 1.5;
smalltri = 0.25;
largetri = 4.0;
qlimit = 0.5;
dt = 0.2;
stats = struct('t_init',0.0,'t_tri',0.0,'t_inpoly',0.0,'t_edge',0.0, ...
't_sparse',0.0,'t_search',0.0,'t_smooth',0.0,'t_density',0.0, ...
'n_tri',0);
% Initialise mesh
% P : Initial nodes
% T : Initial triangulation
% TNDX : Enclosing triangle for each node as indices into TH
% FIX : Indices of FIXED nodes in P
tic
if options.stats
fprintf('Initialising Mesh \n');
end
[p,fix,tndx] = initmesh(p,qtree.p,qtree.t,qtree.h,node,edge);
stats.t_init = toc;
% Main loop
if options.stats
fprintf('Iteration Convergence (%%)\n');
end
for iter = 1:options.maxit
[p,i,j] = unique(p,'rows'); % Ensure unique node list
fix = j(fix);
tndx = tndx(i);
tic
[p,t] = cdt(p,node,edge); % Constrained Delaunay triangulation
stats.n_tri = stats.n_tri+1;
stats.t_tri = stats.t_tri+toc;
tic
e = getedges(t,size(p,1)); % Unique edges
stats.t_edge = stats.t_edge+toc;
% Sparse node-to-edge connectivity matrix
tic
nume = size(e,1);
S = sparse(e(:),[1:nume,1:nume],[ones(nume,1); -ones(nume,1)],size(p,1),nume);
stats.t_sparse = stats.t_sparse+toc;
tic
%
% Find enclosing triangle in background mesh for nodes
%
opt = 4;
if ( opt == 0 )
tndx = mytsearch ( qtree.p(:,1), qtree.p(:,2), qtree.t, p(:,1), p(:,2), tndx );
elseif ( opt == 1 )
tndx = tsearch ( qtree.p(:,1), qtree.p(:,2), qtree.t, p(:,1), p(:,2) );
elseif ( opt == 2 )
tndx = tsearch_mex ( qtree.p(:,1), qtree.p(:,2), qtree.t, p(:,1), p(:,2) );
elseif ( opt == 3 )
dtri = DelaunayTri ( qtree.p );
tndx = pointLocation ( dtri, p );
elseif ( opt == 4 )
tndx = tsearchn ( qtree.p, qtree.t, p );
end
hn = tinterp(qtree.p,qtree.t,qtree.h,p,tndx); % Size function at nodes via linear interpolation
h = 0.5*(hn(e(:,1))+hn(e(:,2))); % from the background mesh. Average to edge midpoints.
stats.t_search = stats.t_search+toc;
edgev = p(e(:,1),:)-p(e(:,2),:);
L = max(sqrt(sum((edgev).^2,2)),eps); % Edge lengths
% Inner smoothing sub-iterations
time = cputime;
move = 1.0;
done = false;
for subiter = 1:(iter-1)
moveold = move;
% Spring based smoothing
L0 = h*sqrt(sum(L.^2)/sum(h.^2));
F = max(L0./L-1.0,-0.1);
F = S*(edgev.*[F,F]);
F(fix,:) = 0.0;
p = p+dt*F;
% Measure convergence
edgev = p(e(:,1),:)-p(e(:,2),:);
L0 = max(sqrt(sum((edgev).^2,2)),eps); % Edge lengths
move = norm((L0-L)./L,'inf'); % Percentage change
L = L0;
if movelargetri*Ah; % Large triangles
r = longest(p,t)./tricentre(t,hn);
k = (r>longedge) & (quality(p,t)0,:);
fix = (1:size(p,1))';
% Initial nodes taken as fixed boundary nodes + internal nodes from
% quadtree.
[i,j] = inpoly(ph,node,edge);
p = [p; ph(i&~j,:)];
tndx = zeros(size(p,1),1);
end % initmesh()
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function e = getedges(t,n)
% Get the unique edges and boundary nodes in a triangulation.
e = sortrows( sort([t(:,[1,2]); t(:,[1,3]); t(:,[2,3])],2) );
idx = all(diff(e,1)==0,2); % Find shared edges
idx = [idx;false]|[false;idx]; % True for all shared edges
bnd = e(~idx,:); % Boundary edges
e = e(idx,:); % Internal edges
e = [bnd; e(1:2:end-1,:)]; % Unique edges and bnd edges for tri's
end % getedges()
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function fc = tricentre(t,f)
% Interpolate nodal F to the centroid of the triangles T.
fc = (f(t(:,1),:)+f(t(:,2),:)+f(t(:,3),:))/3.0;
end % tricentre()
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function d = longest(p,t)
% Return the length of the longest edge in each triangle.
d1 = sum((p(t(:,2),:)-p(t(:,1),:)).^2,2);
d2 = sum((p(t(:,3),:)-p(t(:,2),:)).^2,2);
d3 = sum((p(t(:,1),:)-p(t(:,3),:)).^2,2);
d = sqrt(max([d1,d2,d3],[],2));
end % longest()