function [p,e,t] = initmesh(g, NVArgs)

%INITMESH Build an initial PDE triangular mesh.
%
%   [P,E,T]=INITMESH(G) returns a triangular mesh of the geometry
%   representation given by G. The geometry can be a decomposed geometry
%   matrix DECSG or a geometry representation defined in a MATLAB-file, see
%   PDEGEOM for details. The outputs P, E, and T are the mesh data; node
%   points, edge associativity and triangular elements T. The target element
%   size is determined from the shape of the geometry.
%
%   In the node point matrix P, the first and second rows contain x- and
%   y-coordinates of the points in the mesh.
%
%   In the edge association matrix E, the first and second rows contain indices
%   of the starting and ending point, the third and fourth rows contain
%   the starting and ending parameter values, the fifth row contains
%   the boundary edge segment number, and the sixth and seventh row
%   contain the left- and right-hand side subdomain numbers.
%
%   In the triangle matrix T, the first three rows contain indices to
%   the corner points, given in counter clockwise order, and the fourth
%   row contains the subdomain number.
%
%   [P,E,T]=INITMESH(G,'Property',Value,...) can take property-value pairs to
%   be used by the meshing algorithm. Valid property/value pairs include:
%
%   Property        Value/{Default}         Description
%   -----------------------------------------------------------
%   Hmax            numeric {estimate}      Maximum edge size
%   Hgrad           numeric {1.3}           Triangle growth rate
%   Box             on|{off}                Preserve bounding box
%   Init            on|{off}                Boundary triangulation
%   Jiggle          on|off|{mean}|minimum   Call JIGGLEMESH
%   JiggleIter      numeric {10}            Maximum iterations
%   MesherVersion   {preR2013a}|R2013a      Meshing algorithm used
%
%   The Hmax parameter controls the size of the triangles in the
%   mesh.  INITMESH creates a mesh where no triangle side exceeds
%   Hmax.
%
%   Both the Box and Init parameters are related to the way the
%   mesh algorithm works. By turning on Box you can get a good idea
%   of how the mesh generation algorithm works within the bounding box.
%   By turning on Init you can see the initial triangulation of
%   the boundaries.
%
%   The Jiggle and JiggleIter parameters are used to control whether
%   jiggling of the mesh should be attempted. See JIGGLEMESH for
%   details.
%
%   The MesherVersion parameter selects the specific method used to
%   produce a mesh. Option R2013a selects a Delaunay-based algorithm that is
%   faster and more robust at dealing with non-standard geometries. The
%   default option, preR2013a, selects the method used in previous versions
%   of the toolbox.
%
%   See also DECSG, PDEGEOM, JIGGLEMESH, REFINEMESH, CREATEPDE, pde.PDEModel/generateMesh

%   Copyright 1994-2023 The MathWorks, Inc.

% Analyze input arguments

arguments
    g {mustBeA(g, ["double", "char", "string", "function_handle", "pde.AnalyticGeometry"])}
    NVArgs.Hmax (1,1) {mustBeNonnegative} = 0
    NVArgs.Hgrad (1,1) {mustBeInRange(NVArgs.Hgrad, 1, 2, "exclusive")} = 1.3
    NVArgs.Box {mustBeNonempty, mustBeMember(NVArgs.Box, ["off", "on"])} = "off"
    NVArgs.Init {mustBeNonempty, mustBeMember(NVArgs.Init, ["off", "on"])} = "off"
    NVArgs.Jiggle {mustBeNonempty, mustBeMember(NVArgs.Jiggle, ["mean", "minimum", "on", "off"])} = "mean"
    NVArgs.JiggleIter (1,1) {mustBeInteger, mustBeNonnegative} = 10
    NVArgs.MesherVersion {mustBeNonempty, mustBeMember(NVArgs.MesherVersion, ["preR2013a", "R2013a"])} = "preR2013a"
    NVArgs.MaxRefineSteps (1,1) {mustBeNumeric, mustBeNonnegativeIntOrInf} = Inf % undocumented
    NVArgs.ConstrainedDelaunayRefinement (1,1) {mustBeNumericOrLogical, mustBeMember(NVArgs.ConstrainedDelaunayRefinement, [0 1])} = false % undocumented
end

Hmax        = NVArgs.Hmax;
Hgrad       = NVArgs.Hgrad;
Box         = NVArgs.Box;
Init        = NVArgs.Init;
Jiggle      = NVArgs.Jiggle;
JiggleIter  = NVArgs.JiggleIter;
useDt       = strcmpi(NVArgs.MesherVersion, 'R2013a');
maxRefineSteps = NVArgs.MaxRefineSteps; % undocumented
constrainedDelaunayRefine = NVArgs.ConstrainedDelaunayRefinement; % undocumented

if isa(g, 'pde.AnalyticGeometry')
    g = g.geom;
end

% Start by creating a coarse mesh

% Sort out mesh geometry
[p,e,Hmin,factor] = pdemgeom(g);
chkTopology(p, e);

% Check the geometry for intersecting edges
delTri = pdeDelTri(p');
cm = getEdgeConstraints(e, false);
warnState(1) = warning('off','MATLAB:DelaunayTri:ConsSplitPtWarnId');
warnState(2) = warning('off', 'MATLAB:DelaunayTri:ConsConsSplitWarnId');
warnState(3) = warning('off', 'MATLAB:DelaunayTri:DupConsWarnId');
warnState(4) = warning('off', 'MATLAB:DelaunayTri:LoopConsWarnId');
delTri.setConstraints(cm);
warning(warnState);
numCmOld = size(cm,1);
numCmNew = size(delTri.getConstraints(),1);

if(numCmOld ~= numCmNew)
    checkSelfIntersections(g)
end

% Design bounding box
x               = p(1,:);
[xmin,xmax]     = bounds(x);
xdiff           = xmax-xmin;
x1              = xmin-xdiff;
x2              = xmax+xdiff;
y               = p(2,:);
[ymin,ymax]     = bounds(y);
ydiff           = ymax-ymin;
y1              = ymin-ydiff;
y2              = ymax+ydiff;
p               = [x1 x2 x2 x1; y1 y1 y2 y2];
e(1:2,:)        = e(1:2,:)+4;
t               = [1 1; 2 3; 3 4; 0 0];

if ~Hmax
    Hmax = 0.1*max(xdiff,ydiff);
end

% Display warning message for large number of triangles
if Hmax~=inf
    ntg = floor(2*(xmax-xmin)*(ymax-ymin)/Hmax^2);
    if ntg>20000
        warning(message('pdelib:initmesh:LargeMesh', ntg));
    end
end

% determine scale and tolerance of model
small = 10000*eps;
scale = max(xmax-xmin,ymax-ymin);
tol   = small*scale;
tol2  = small*scale^2;

% Triangulate edges
if useDt
    % Create a new triangulation without constraints
    delTri = pdeDelTri(p');
    delTri.pdevoron(x, y);
    p = delTri.getPoints();
    t = delTri.getTriangles();
    c = delTri.getCircumcenters();
else
    [p,t,c] = pdevoron(p,t,[],Hmax,x,y,tol,Hmax,Hgrad);
end

% Make sure that triangulation respects boundaries
[p,e,t,c] = pderespe(g,p,e,t,c,Hmax,tol,tol2,Hmax,Hgrad);

% Do we want Hmax==inf with bounding box ?
if Hmax==inf && strcmp(Box,'on')
    return
end

% Determine internal triangles
[k,t] = pdeintrn(p,e,t);
it = t(:,k);

% Do we want Hmax==inf?
if Hmax==inf
    % Remove outside points
    [p,e,t] = pdermpnt(p,e,it);
    return
end

% Check for narrow geometries
h = pdehloc(p,e,it,Hmax,Hmin,Hgrad);

% Distribute edge points
[x,y,e,h] = pdedistr(g,p,e,it,Hmax,Hgrad,tol,h,factor);

% Triangulate edges
if ~useDt
    try
        [p,t,c,h] = pdevoron(p,t,[],h,x,y,tol,Hmax,Hgrad);
    catch err
        pdevoronFailureMessage(err);
    end
    % Make sure that triangulation respects boundaries
    [p,e,t,c,h] = pderespe(g,p,e,t,c,h,tol,tol2,Hmax,Hgrad);
else
    delTri = pdeDelTri(p');
    [~,~,c] = delTri.pdevoron(x, y);
    cm = getEdgeConstraints(e, false);
    delTri.setConstraints(cm);
    p = delTri.getPoints();
    t = delTri.getTriangles();
end

% Determine internal triangles
[k, t] = pdeRegionCalc(e, t);

it = t(:,k);
ic = c(:,k);

% Do we want initial mesh with bounding box?
if strcmp(Init,'on') && strcmp(Box,'on')
    return
end

% Do we want initial mesh?
if strcmp(Init,'on')
    % Remove outside points
    [p,e,t] = pdermpnt(p,e,it);
    return
end

if useDt
    delTri = pdeDelTri(p');
    cm = getEdgeConstraints(e, ~constrainedDelaunayRefine);
    delTri.setConstraints(cm);
    p = delTri.getPoints();
    t = delTri.getTriangles();
    c = delTri.getCircumcenters();
    hEdge = calcEdgeH(p, e);
    hInterpol = delTri.buildSizeInterp(hEdge,Hmax,Hgrad);
    k = delTri.getInOutStatus();
    it = t(:,k);
    ic = c(:,k);
end

t(:,k) = [];
c(:,k) = [];
ne = size(t,2);
t = [t,it];
c = [c,ic];
% Determine large triangles inside geometry and refine
numRefine = 0;
while true
    it1 = it(1,:);
    it2 = it(2,:);
    it3 = it(3,:);

    % Triangle side lengths
    l1 = (p(1,it2)-p(1,it1)).^2 + (p(2,it2)-p(2,it1)).^2;
    l2 = (p(1,it3)-p(1,it2)).^2 + (p(2,it3)-p(2,it2)).^2;
    l3 = (p(1,it1)-p(1,it3)).^2 + (p(2,it1)-p(2,it3)).^2;

    % Triangle size
    maxTriLen = sqrt(max(l1,max(l2,l3)));

    % Compute angles
    a = pdehloc(p,e,it);

    % Finished?
    if useDt
        if size(h,2)>1
            hit1 = hInterpol(p(:,it1)')';
            hit2 = hInterpol(p(:,it2)')';
            hit3 = hInterpol(p(:,it3)')';
            hh = prod([hit1;hit2;hit3]).^(1/3);
        else
            hh = h;
        end
        enclosingTri = delTri.getEnclosingTriangles(ic(1:2,:)');
        allTri = delTri.getTriangles();
        net = pdeFilterInsertPts(p, e, allTri, ic, enclosingTri);
        k = find(maxTriLen>=1.2*hh & min(abs(a-pi/2))>pi/8 & net);
    else
        hh = h;
        if size(h,2)>1
            hh = prod([h(it1);h(it2);h(it3)]).^(1/3);
        end
        k  = find(maxTriLen>=1.2*hh & min(abs(a-pi/2))>pi/8);
    end

    if isempty(k)
        break
    end

    % Remaining triangles
    ic = ic(:,k);
    maxTriLen = maxTriLen(k);

    % Choose center of circumscribed circle of largest triangles
    % Ensure that we do not select points inside previous circles
    [~,si] = sort(maxTriLen);

    sj = filterTrianglesByCircumcenter(ic, si, tol); % mex version

    x = ic(1,sj);
    y = ic(2,sj);

    % Something is going wrong - ignore!
    i = (x>x2 | x<x1 | y>y2 | y<y1);
    x(i) = [];
    y(i) = [];
    if isempty(x)
        break
    end

    % Triangulate
    if ~useDt
        try
            [pNew,t1,c1,h]=pdevoron(p,t,c,h,x,y,tol,Hmax,Hgrad);
        catch err
            pdevoronFailureMessage(err);
        end
        npNew = size(pNew,2) - size(p,2);
        p = pNew;

        % Check if outer triangles changed (internal boundaries may not be
        % respected).
        if any(any(t(:,1:ne)~=t1(:,1:ne)))
            % Make sure that triangulation respects boundaries
            [p,e,t,c,h] = pderespe(g,p,e,t1,c1,h,tol,tol2,Hmax,Hgrad);

            % Determine internal triangles
            k = pdeintrn(p,e,t);

            it = t(:,k);
            ic = c(:,k);
            t(:,k) = [];
            c(:,k) = [];
            ne = size(t,2);
            t  = [t,it]; %#ok<AGROW>
            c  = [c,ic]; %#ok<AGROW>
        else
            t  = t1;
            c  = c1;
            it = t(:,ne+1:size(t,2));
            ic = c(:,ne+1:size(t,2));
        end
    else
        [p,t1,c1,npNew] = delTri.pdevoron(x, y);
        if constrainedDelaunayRefine
            k = pdeRegionCalc(e, t);
        else
            k = delTri.getInOutStatus();
        end
        it = t1(:,k);
        ic = c1(:,k);
        t  = t1;
        c  = c1;
    end

    numRefine = numRefine + 1;

    if (numRefine>maxRefineSteps || npNew==0)
        break
    end

end

if ~useDt
    % Make sure that triangulation respects internal boundaries
    [p,e,t,~,~] = pderespe(g,p,e,t,c,h,tol,tol2,Hmax,Hgrad);
else
    % add all constraints
    cm = getEdgeConstraints(e, false);
    delTri.setConstraints(cm);
    p = delTri.getPoints();
    t = delTri.getTriangles();
end

% Select internal triangles
[k, t] = pdeRegionCalc(e, t);


if strcmp(Box,'on')
    return
end

% Remove outside points
t = t(:,k);
[p,e,t] = pdermpnt(p,e,t);

if ~strcmp(Jiggle,"off")
    if strcmp(Jiggle,"on")
        Opt = "off";
    else
        Opt = Jiggle;
    end
    p = jigglemesh(p,e,t,"Iter",JiggleIter,"Opt",Opt);
end
end % END initmesh

function c = getEdgeConstraints(e, excludeInner)
if excludeInner
    eOuter = (e(6,:)==0 | e(7,:)==0);
    c = [e(1,eOuter)' e(2,eOuter)'];
else
    c = [e(1,:)' e(2,:)'];
end
end % END getEdgeConstraints

function h = calcEdgeH(p, e)
v1 = e(1,:);
v2 = e(2,:);
x1 = p(:,v1);
x2 = p(:,v2);
dx = x2 - x1;
hEdge = sqrt(dx(1,:).*dx(1,:) + dx(2,:).*dx(2,:));
h = zeros(1, size(p,2));
h(v1) = hEdge;
h(v2) = hEdge;
end % END calcEdgeH

function chkTopology(p, e)
np = size(p,2);
T = sparse(e(1,:), e(2,:), 1, np, np);
T = T+T';
singleEdgeNodes = find(sum(T)==1);
if ~isempty(singleEdgeNodes)
    xyInterleaved = [p(1,singleEdgeNodes); p(2,singleEdgeNodes)];
    xyLocs = [newline sprintf('(%g,%g)\n',xyInterleaved)];
    error(message('pdelib:initmesh:InvalidGeometry', xyLocs))
end
end % END chkTopology

function pdevoronFailureMessage(err)
if strcmp(err.identifier,'pdelib:pdevoron:GeomError')
    warnState = warning('backtrace', 'off');
    warning(message('pdelib:initmesh:TryNewMesher'));
    warning(warnState);
end
rethrow(err);
end % END pdevoronFailureMessage

function mustBeNonnegativeIntOrInf(MaxRefineSteps)
if MaxRefineSteps~=Inf
    mustBeNonnegative(MaxRefineSteps)
    mustBeInteger(MaxRefineSteps)
end
end % END mustBeNonnegativeIntOrInf