A technique to improve quality of elements on the boundary

[krober10nd]

• Using the functions below, one can discretize a polygon so the spacing of its nodes follows a sizing function.
• Along each boundary segment of the polygon, you can place a new node at a distance from the segment's midpoint so a equilateral triangle is formed.
• During mesh generation, these points associated with the boundary are fixed so they cannot move. In other words you can pass these points as fixed. Note that you should remove any existing initialized points in the vicinity of these fixed points. Second note: in areas where the boundary folds in on itself more work is needed to avoid degenerate triangles.

The result is dramatically improved mesh quality along the boundary

function [pout, t, converged] = mesh1d(poly, edgeFuncHandles, minEdgeLength, fixedPoints, boundaryBox, boxNum, varargin)
% Mesh Generator using Distance Functions.
%
% Inputs:
%   poly - Nx2 array of polygon vertices.
%   edgeFuncHandles - Cell array of functions for edge length determination.
%   minEdgeLength - Minimum edge length.
%   fixedPoints - Indices of fixed polygon vertices.
%   boundaryBox - Boundary box for mesh.
%   boxNum - Index of the boundary box for initial edge length.
%
% Outputs:
%   pout - Nx2 array of resampled node positions.
%   t - Edge indices (NTx2).
%   converged - Boolean flag indicating convergence.

% Initialize variables
X = poly(:,1);
Y = poly(:,2);
xd = diff(X);
yd = diff(Y);
distances = sqrt(xd.^2 + yd.^2);
cumulativeDist = [0; cumsum(distances)];
totalDist = cumulativeDist(end);

% Calculate number of points and parameterize
numPoints = totalDist / min(minEdgeLength);
t = linspace(0, max(cumulativeDist), ceil(numPoints));

% Interpolate new points
newX = interp1(cumulativeDist, X, t);
newY = interp1(cumulativeDist, Y, t);

% Prepare for mesh generation
p = [newX', newY']; % Initial point distribution
fixedIndices = unique([1, length(p)]); % Ensure the start and end points are fixed
p(fixedIndices, :) = poly(fixedPoints, :); % Apply fixed points

% Distance function for the domain
fdist = @(x, y) sqrt(min((x - boundaryBox(:,1)).^2 + (y - boundaryBox(:,2)).^2));

% Edge length function (simplified for demonstration)
fh = @(x, y) minEdgeLength;

% Main mesh generation loop
[p, t, converged] = generateMesh(p, fdist, fh, fixedIndices, totalDist, minEdgeLength);

% Map the final points back to the original coordinates
pout = [interp1(cumulativeDist, X, p(:,1)), interp1(cumulativeDist, Y, p(:,2))];

if any(isnan(pout))
    warning('NaNs detected in 1D mesh');
end

end

function [p, t, converged] = generateMesh(p, fdist, fh, fixedIndices, totalDist, minEdgeLength)
    % Initialize variables
    converged = false;
    maxIterations = 100;
    iteration = 0;
    ptol = 1e-3; % Tolerance for point movement
    
    % Start the iteration
    while ~converged && iteration < maxIterations
        iteration = iteration + 1;
        
        % Calculate edge lengths and desired adjustments
        edges = diff(p, 1, 1); % Differences between points
        edgeLengths = sqrt(sum(edges.^2, 2));
        desiredLengths = fh(p(:,1), p(:,2)); % Assuming fh returns a scalar or vector the same length as p
        
        % Calculate movement for each point based on edge length discrepancy
        moveDist = (desiredLengths - edgeLengths) .* (edges ./ edgeLengths);
        moveDist = [0, 0; moveDist]; % No movement for the first point
        
        % Apply movements, excluding fixed points
        isFixed = false(size(p, 1), 1);
        isFixed(fixedIndices) = true;
        p(~isFixed, :) = p(~isFixed, :) + moveDist(~isFixed, :);
        
        % Check for convergence
        maxMove = max(sqrt(sum(moveDist.^2, 2)));
        if maxMove < ptol
            converged = true;
        end
    end
    
    % Generate edge indices (for visualization or further processing)
    t = [(1:size(p, 1)-1)', (2:size(p, 1))'];
end

function [initial_points, layer_of_points] = form_initial_points_near_and_along_boundary(outer_polygon, ...
    bounding_polygon, ...
    minimum_resolution_in_meters, ...
    mesh_sizing_fx)
% Form a layer of points one row along the boundary rim of the domain
% so that when triangulated it forms perfectly equilateral triangles. 

minimum_resolution_in_degrees = minimum_resolution_in_meters / 111e3; 

% ensure the first point is not repeated etc.
outer_polygon = unique(outer_polygon,'rows','stable'); 

%% Resample the outer polygon to match the point spacing following mesh_sizing_fx 
[outer_polygon_rs,~,~]=mesh1d(outer_polygon, ...
                                                   mesh_sizing_fx, ...
                                                   minimum_resolution_in_degrees, ...
                                                   [], ...
                                                   bounding_polygon, ...
                                                   1, ... % box_number (assumed 1). 
                                                   [] ...
                                                   );
layer_of_points = [];

%figure;

for i = 1 : length(outer_polygon_rs)-1
    % 1) Compute the midpoint
    p1 = outer_polygon_rs(i,:);
    p2 = outer_polygon_rs(i+1,:);
    midpoint = (p1 + p2) / 2.0;
    
    % 2) Compute the tangent and then the unit normal vector pointing inwards
    tangent = p2 - p1; % Tangent vector
    tangent_length = norm(tangent); % Length of the tangent vector
    unit_tangent = tangent / tangent_length; % Normalize to get unit tangent
    unit_normal = [-unit_tangent(2), unit_tangent(1)]; % Rotate by 90 degrees to get the unit normal
    
    % 3) Calculate the distance for the equilateral triangle and compute p3
    local_size_in_degrees = mesh_sizing_fx{1}(midpoint) / 111e3; 
    % The distance from the midpoint to the third point to form an equilateral triangle
    % Assuming local_size_in_degrees gives the desired edge length of the equilateral triangle
    height_of_triangle = (sqrt(3) / 2) * local_size_in_degrees; % Height of an equilateral triangle
    p3 = midpoint + height_of_triangle * unit_normal; % New point forming equilateral triangle

    % For debugging plot
%     hold on; 
%     plot([p1(1), p2(1)], [p1(2), p2(2)], 'k-'); % Plot edge
%     plot(midpoint(1), midpoint(2), 'k.'); % Plot midpoint
%     quiver(midpoint(1), midpoint(2), unit_normal(1), unit_normal(2), 0.05, 'r'); % Plot unit normal
%     plot(p3(1), p3(2), 'bs'); % Plot new point

    layer_of_points = [layer_of_points; p3];
end

% Combine initial resampled outer_polygon with layer_of_points 
initial_points = [outer_polygon_rs; layer_of_points];

end

[oceanywalker]

It's excellent work, can you provide a case to show how it works?

[krober10nd]

hi @oceanywalker I haven't had much time but recently I applied to this a project mesh and it worked well. The red points are the initially formed triangles that become locked/fixed. Donations of any amount are always welcome. Cheers.

See where it says "form_initial_points_along_boundary" below. All other functions should be part of the latest release on the Projection branch.

classdef meshgen
    %   MESHGEN: Mesh generation class
    %   Handles input parameters to create a meshgen class object that can be
    %   used to build a msh class.
    %   Copyright (C) 2018  Keith Roberts & William Pringle
    %
    %   This program is free software: you can redistribute it and/or modify
    %   it under the terms of the GNU General Public License as published by
    %   the Free Software Foundation, either version 3 of the License, or
    %   (at your option) any later version.
    %
    %   This program is distributed in the hope that it will be useful,
    %   but WITHOUT ANY WARRANTY; without even the implied warranty of
    %   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
    %   GNU General Public License for more details.
    %
    %   You should have received a copy of the GNU General Public License
    %   along with this program.  If not, see <http://www.gnu.org/licenses/>.
    %
    %   Available options:
    %         ef            % edgefx class
    %         bou           % geodata class
    %         h0            % minimum edge length (optional if bou exists)
    %         bbox          % bounding box [xmin,ymin; xmax,ymax] (manual specification, no bou)
    %         proj          % structure containing the m_map projection info
    %         plot_on       % flag to plot (def: 1) or not (0)
    %         nscreen       % how many it to plot and write temp files (def: 5)
    %         itmax         % maximum number of iterations.
    %         pfix          % fixed node positions (nfix x 2 )
    %         egfix         % edge constraints
    %         outer         % meshing boundary (manual specification, no bou)
    %         inner         % island boundaries (manual specification, no bou)
    %         mainland      % the shoreline boundary (manual specification, no bou)
    %         fixboxes      % a flag that indicates which boxes will use fixed constraints
    %         memory_gb     % memory in GB allowed to use for initial rejector
    %         cleanup       % logical flag or string to trigger cleaning of topology (default is on).
    %         direc_smooth  % logical flag to trigger direct smoothing of mesh in the cleanup
    %         dj_cutoff     % the cutoff area fraction for disjoint portions to delete
    %         qual_tol      % tolerance for the accepted negligible change in quality
    %         enforceWeirs  % whether or not to enforce weirs in meshgen
    %         enforceMin    % whether or not to enfore minimum edgelength for all edgefxs
    %         big_mesh      % set to 1 to remove the bou data from memory
    %      improve_boundary  % run a gradient desc. to improve the boundary conformity
    properties
        fd            % handle to distance function
        fh            % handle to edge function
        h0            % minimum edge length
        edgefx        % edgefx class
        bbox          % bounding box [xmin,ymin; xmax,ymax]
        pfix          % fixed node positions (nfix x 2 )
        egfix         % edge constraints
        fixboxes      % a flag that indicates which boxes will use fixed constraints
        plot_on       % flag to plot (def: 1) or not (0)
        nscreen       % how many it to plot and write temp files (def: 5)
        bou           % geodata class
        ef            % edgefx class
        itmax         % maximum number of iterations.
        outer         % meshing boundary (manual specification, no bou)
        inner         % island boundaries (manual specification, no bou)
        mainland      % the shoreline boundary (manual specification, no bou)
        boubox        % the bbox as a polygon 2-tuple
        inpoly_flip   % used to flip the inpoly test to determine the signed distance.
        memory_gb     % memory in GB allowed to use for initial rejector
        cleanup       % logical flag or string to trigger cleaning of topology (default is on).
        direc_smooth  % logical flag to trigger direct smoothing of mesh in the cleanup
        dj_cutoff     % the cutoff area fraction for disjoint portions to delete
        grd = msh();  % create empty mesh class to return p and t in.
        big_mesh      % release bou data from memory
        qual          % mean, lower 3rd sigma, and the minimum element quality.
        qual_tol      % tolerance for the accepted negligible change in quality
        proj          % structure containing the m_map projection info
        anno          % Approx. Nearest Neighbor search object.
        annData       % datat contained with KD-tree in anno
        Fb            % bathymetry data interpolant
        enforceWeirs  % whether or not to enforce weirs in meshgen
        enforceMin    % whether or not to enfore minimum edgelength for all edgefxs
        improve_boundary  % improve the boundary representation
        high_fidelity     % flag to form pfix and egfix for this domain
    end
    
    
    methods
       
        function obj = plot(obj)
            if ~isempty(obj.pfix) && ~isempty(obj.egfix)
                figure;
                % plot the bounding boxes & shorelines for visual aid
                for box_number = 1 : length(obj.boubox)
                    
                    iboubox = obj.boubox{box_number};
                    hold on; plot(iboubox(:,1),iboubox(:,2),'g-','linewidth',2);
                    
                    touter = obj.outer(box_number);
                    tedges = Get_poly_edges(touter{1});
                    [touter, ~]=filter_constraints(touter{1},tedges,obj.boubox,box_number);
                    hold on; plot(touter(:,1),touter(:,2),'bs');
   
                    tinner = obj.inner(box_number);
                    if ~isempty(tinner{1})
                        tedges = Get_poly_edges(tinner{1});
                        [tinner, ~]=filter_constraints(tinner{1},tedges,obj.boubox,box_number);
                        hold on; plot(tinner(:,1),tinner(:,2),'rs');
                    end
                    
                end
                
                % plot the fixed constraints
                drawedge2(obj.pfix,obj.egfix,[0,0,0]);
                hold on; plot(obj.pfix(:,1),obj.pfix(:,2),'ks','MarkerSize',15);
                axis equal;
                title('Resampled Point & Edge Constraints');
            else
                disp('No constraints to plot!')
            end
        end
        
        % class constructor/default grd generation options
        function obj = meshgen(varargin)
            rng(0); 
            % Check for m_map dir
            M_MAP_EXISTS=0;
            if exist('m_proj','file')==2
                M_MAP_EXISTS=1 ;
            end
            if M_MAP_EXISTS~=1
                error('Where''s m_map? Chief, you need to read the user guide')
            end
            
            % Check for utilties dir
            UTIL_DIR_EXISTS=0 ;
            if exist('inpoly.m','file')
                UTIL_DIR_EXISTS=1;
            end
            if UTIL_DIR_EXISTS~=1
                error('Where''s the utilities directory? Chief, you need to read the user guide')
            end
            
            p = inputParser;
            % unpack options and set default ones, catch errors.
            
            defval = 0; % placeholder value if arg is not passed.
            % add name/value pairs
            addOptional(p,'h0',defval);
            addOptional(p,'bbox',defval);
            addOptional(p,'fh',defval);
            addOptional(p,'pfix',defval);
            addOptional(p,'egfix',defval);
            addOptional(p,'fixboxes',defval);
            addOptional(p,'inner',defval);
            addOptional(p,'outer',defval);
            addOptional(p,'mainland',defval);
            addOptional(p,'bou',defval);
            addOptional(p,'ef',defval);
            addOptional(p,'plot_on',defval);
            addOptional(p,'nscreen',defval);
            addOptional(p,'itmax',defval);
            addOptional(p,'memory_gb',1);
            addOptional(p,'cleanup',1);
            addOptional(p,'direc_smooth',1);
            addOptional(p,'dj_cutoff',0.25);
            addOptional(p,'big_mesh',defval);
            addOptional(p,'proj',defval);
            addOptional(p,'qual_tol',defval);
            addOptional(p,'enforceWeirs',0);
            addOptional(p,'enforceMin',1);
            
            % parse the inputs
            parse(p,varargin{:});
            %if isempty(varargin); return; end
            % store the inputs as a struct
            inp = p.Results;
            % kjr...order these argument so they are processed in a predictable
            % manner. Process the general opts first, then the OceanMesh
            % classes...then basic non-critical options.
            inp = orderfields(inp,{'h0','bbox','enforceWeirs','enforceMin','fh',...
                'inner','outer','mainland',...
                'bou','ef',... %<--OceanMesh classes come after
                'egfix','pfix','fixboxes',...
                'plot_on','nscreen','itmax',...
                'memory_gb','qual_tol','cleanup',...
                'direc_smooth','dj_cutoff',...
                'big_mesh','proj'});
            % get the fieldnames of the edge functions
            fields = fieldnames(inp);
            % loop through and determine which args were passed.
            % also, assign reasonable default values if some options were
            % not assigned.
            for i = 1 : numel(fields)
                type = fields{i};
                switch type
                    % parse aux options first
                    case('h0')
                        % Provide in meters
                        obj.h0 = inp.(fields{i});
                    case('fh')
                        if isa(inp.(fields{i}),'function_handle')
                            obj.fh = inp.(fields{i});
                        end
                        % can't check for errors here yet.
                    case('bbox')
                        obj.bbox= inp.(fields{i});
                        if iscell(obj.bbox)
                            % checking bbox extents
                            ob_min = obj.bbox{1}(:,1);
                            ob_max = obj.bbox{1}(:,2);
                            for ii = 2:length(obj.bbox)
                                if any(obj.bbox{ii}(:,1) < ob_min) || ...
                                        any(obj.bbox{ii}(:,2) > ob_max)
                                    error(['Outer bbox must contain all ' ...
                                        'inner bboxes: inner box #' ...
                                        num2str(ii) ' violates this'])
                                end
                            end
                        end
                        
                        % if user didn't pass anything explicitly for
                        % bounding box make it empty so it can be populated
                        % from ef as a cell-array
                        if obj.bbox(1)==0
                            obj.bbox = [];
                        end
                    case('pfix')
                        obj.pfix= inp.(fields{i});
                        if obj.pfix(1)~=0
                            obj.pfix(:,:) = inp.(fields{i});
                        else
                            obj.pfix = [];
                        end
                        if obj.enforceWeirs
                            for j = 1 : length(obj.bou)
                                if  ~isempty(obj.bou{j}.weirPfix)
                                    obj.pfix = [obj.pfix ; obj.bou{j}.weirPfix];
                                end
                            end
                        end
                    case('egfix')
                        obj.egfix= inp.(fields{i});
                        if ~isempty(obj.egfix) && obj.egfix(1)~=0
                            obj.egfix = inp.(fields{i});
                        else
                            obj.egfix = [];
                        end
                        if obj.enforceWeirs
                            for j = 1 : length(obj.bou)
                                if ~isempty(obj.bou{j}.weirEgfix) && ~isempty(obj.egfix)
                                    obj.egfix = [obj.egfix ; obj.bou{j}.weirEgfix+max(obj.egfix(:))];
                                else
                                    obj.egfix =  obj.bou{j}.weirEgfix;
                                end
                            end
                        end
                        obj.egfix = renumberEdges(obj.egfix);
                    case('fixboxes')
                        obj.fixboxes= inp.(fields{i});
                    case('bou')
                        % got it from user arg
                        if obj.outer~=0, continue; end
                        
                        obj.outer = {} ;
                        obj.inner = {} ;
                        obj.mainland = {} ;
                        
                        obj.bou = inp.(fields{i});
                        
                        % handle when not a cell
                        if ~iscell(obj.bou)
                            boutemp = obj.bou;
                            obj.bou = cell(1);
                            obj.bou{1} = boutemp;
                        end
                        
                        % then the geodata class was provide, unpack
                        for ee = 1:length(obj.bou)
                            try
                                arg = obj.bou{ee} ;
                            catch
                                arg = obj.bou;
                            end
                            if isa(arg,'geodata')
                                obj.high_fidelity{ee} = obj.bou{ee}.high_fidelity;
                                obj.outer{ee} = obj.bou{ee}.outer;
                                obj.inner{ee} = obj.bou{ee}.inner;
                                
                                % save bathy interpolant to meshgen
                                if ~isempty(obj.bou{ee}.Fb)
                                    obj.Fb{ee} = obj.bou{ee}.Fb ;
                                end
                                
                                if ~isempty(obj.inner{ee}) && ...
                                        obj.inner{ee}(1)~= 0
                                    obj.outer{ee} = [obj.outer{ee};
                                        obj.inner{ee}];
                                end
                                obj.mainland{ee} = obj.bou{ee}.mainland;
                                obj.boubox{ee} = obj.bou{ee}.boubox;
                                obj.inpoly_flip{ee} = obj.bou{ee}.inpoly_flip;
                                if obj.big_mesh
                                    % release gdat's
                                    obj.bou{ee}.mainland= [];
                                    obj.bou{ee}.outer= [];
                                    if ~isempty(obj.bou{ee}.inner)
                                        obj.bou{ee}.inner= [];
                                    end
                                end
                            end
                        end
                        
                    case('ef')
                        tmp = inp.(fields{i});
                        if isa(tmp, 'function_handle')
                            error('Please specify your edge function handle through the name/value pair fh');
                        end
                        obj.ef = tmp;
                        
                        % handle when not a cell
                        if ~iscell(obj.ef)
                            eftemp = obj.ef;
                            obj.ef = cell(1);
                            obj.ef{1} = eftemp;
                        end
                        
                        % Gather boxes from ef class.
                        for ee = 1 : length(obj.ef)
                            if isa(obj.ef{ee},'edgefx')
                                obj.bbox{ee} = obj.ef{ee}.bbox;
                            end
                        end
                        
                        % checking bbox extents
                        if iscell(obj.bbox)
                            ob_min = obj.bbox{1}(:,1);
                            ob_max = obj.bbox{1}(:,2);
                            for ii = 2:length(obj.bbox)
                                if any(obj.bbox{ii}(:,1) < ob_min) || ...
                                        any(obj.bbox{ii}(:,2) > ob_max)
                                    error(['Outer bbox must contain all ' ...
                                        'inner bboxes: inner box #' ...
                                        num2str(ii) ' violates this'])
                                end
                            end
                        end
                        
                        % kjr 2018 June: get h0 from edge functions
                        for ee = 1:length(obj.ef)
                            if isa(obj.ef{ee},'edgefx')
                                obj.h0(ee) = obj.ef{ee}.h0;
                            end
                        end
                        
                        % kjr 2018 smooth the outer automatically
                        if length(obj.ef) > 1
                            % kjr 2020, ensure the min. sizing func is
                            % used
                            if obj.enforceMin
                                obj.ef = enforce_min_ef(obj.ef);
                            end
                            obj.ef = smooth_outer(obj.ef,obj.Fb);
                        end
                        
                        % Save the ef interpolants into the edgefx
                        for ee = 1:length(obj.ef)
                            if isa(obj.ef{ee},'edgefx')
                                obj.fh{ee} = @(p)obj.ef{ee}.F(p);
                            end
                        end
                        
                    case('plot_on')
                        obj.plot_on= inp.(fields{i});
                    case('big_mesh')
                        obj.big_mesh = inp.(fields{i});
                    case('nscreen')
                        obj.nscreen= inp.(fields{i});
                        if obj.nscreen ~=0
                            obj.nscreen = inp.(fields{i});
                            obj.plot_on = 1;
                        else
                            obj.nscreen = 5; % default
                        end
                    case('itmax')
                        obj.itmax= inp.(fields{i});
                        if obj.itmax ~=0
                            obj.itmax = inp.(fields{i});
                        else
                            obj.itmax = 100;
                            warning('No itmax specified, itmax set to 100');
                        end
                    case('qual_tol')
                        obj.qual_tol = inp.(fields{i});
                        if obj.qual_tol ~=0
                            obj.qual_tol = inp.(fields{i});
                        else
                            obj.qual_tol = 0.01;
                        end
                    case('inner')
                        if ~isa(obj.bou,'geodata')
                            obj.inner = inp.(fields{i});
                        end
                    case('outer')
                        if ~isa(obj.bou,'geodata')
                            obj.outer = inp.(fields{i});
                            if obj.inner(1)~=0
                                obj.outer = [obj.outer; obj.inner];
                            end
                        end
                    case('mainland')
                        if ~isa(obj.bou,'geodata')
                            obj.mainland = inp.(fields{i});
                        end
                    case('memory_gb')
                        if ~isa(obj.bou,'memory_gb')
                            obj.memory_gb = inp.(fields{i});
                        end
                    case('cleanup')
                        obj.cleanup = inp.(fields{i});
                        if isempty(obj.cleanup) || obj.cleanup == 0
                            obj.cleanup = 'none';
                        elseif obj.cleanup == 1
                            obj.cleanup = 'default';
                        end
                    case('dj_cutoff')
                        obj.dj_cutoff = inp.(fields{i});
                    case('direc_smooth')
                        obj.direc_smooth = inp.(fields{i});
                    case('proj')
                        obj.proj = inp.(fields{i});
                        % default CPP
                        if obj.proj == 0; obj.proj = 'equi'; end
                        if ~isempty(obj.bbox)
                            % kjr Oct 2018 use outer coarsest box for
                            % multiscale meshing
                            lon_mi = obj.bbox{1}(1,1)-obj.h0(1)/1110;
                            lon_ma = obj.bbox{1}(1,2)+obj.h0(1)/1110;
                            lat_mi = obj.bbox{1}(2,1)-obj.h0(1)/1110;
                            lat_ma = obj.bbox{1}(2,2)+obj.h0(1)/1110;
                        else
                            lon_mi = -180; lon_ma = 180;
                            lat_mi = -90; lat_ma = 90;
                        end
                        % Set up projected space
                        dmy = msh() ;
                        dmy.p(:,1) = [lon_mi; lon_ma];
                        dmy.p(:,2) = [lat_mi; lat_ma];
                        del = setProj(dmy,1,obj.proj) ;
                    case('enforceWeirs')
                        obj.enforceWeirs = inp.(fields{i});
                    case('enforceMin')
                        obj.enforceMin = inp.(fields{i});
                end
            end
            
            if isempty(varargin); return; end
            
            % error checking
            if isempty(obj.boubox) && ~isempty(obj.bbox)
                % Make the bounding box 5 x 2 matrix in clockwise order if
                % it isn't present. This case must be when the user is
                % manually specifying the PSLG.
                obj.boubox{1} = [obj.bbox(1,1) obj.bbox(2,1);
                    obj.bbox(1,1) obj.bbox(2,2); ...
                    obj.bbox(1,2) obj.bbox(2,2);
                    obj.bbox(1,2) obj.bbox(2,1); ...
                    obj.bbox(1,1) obj.bbox(2,1); NaN NaN];
            end
            if any(obj.h0==0), error('h0 was not correctly specified!'), end
            if isempty(obj.outer), error('no outer boundary specified!'), end
            if isempty(obj.bbox), error('no bounding box specified!'), end
            obj.fd = @dpoly;  % <-default distance fx accepts p and pv (outer polygon).
            % kjr build ANN object into meshgen
            obj = createANN(obj) ;
            
            global MAP_PROJECTION MAP_COORDS MAP_VAR_LIST
            obj.grd.proj    = MAP_PROJECTION ;
            obj.grd.coord   = MAP_COORDS ;
            obj.grd.mapvar  = MAP_VAR_LIST ;
            
            % if high_fidelity option was passed, form the pfix and egfix
            % for the outer geometry inside that domain.
            tpfix = []; tegfix= [];
            for box_num = 1:length(obj.h0)
                % only constrain mainland and inner geometries.
                ml = obj.mainland{box_num};
                il = obj.inner{box_num};
                if ~isempty(ml) && ~isempty(il)
                    polys = {ml; il};
                elseif ~isempty(ml) && isempty(il)
                    polys = {ml};
                else
                    polys = {il};
                end
                % user-reserves the option to choose to constrain or not.
                if obj.high_fidelity{box_num}
                    if obj.cleanup==1
                        warning('Setting clean up to 0 since high_fidelity mode is on');
                        obj.cleanup = 0;
                    end
                    disp('Redistributing vertices along 1D constraints...');
                    disp(['    forming pfix and egfix for segments in box #' num2str(box_num)]);
                    for pid = 1 : numel(polys)
                        poly = polys{pid};
                        idx = all(isnan(poly),2);
                        idr = diff(find([1;diff(idx);1]));
                        D = mat2cell(poly,idr(:),size(poly,2));
                        nd = length(D);
                        for i = 1 : 2: nd
                            poly  = D{i};
                            if length(poly) > 2
                                % Filter the constrainst to the bbox.
                                % This creates potentially disjoint linestrings
                                % that need to be induvidually meshed in 1d.
                                tedges = Get_poly_edges([poly; NaN NaN]);
                                [pts, bnde] = filter_constraints(poly,tedges,obj.boubox, box_num);
                                if length(bnde) < 1
                                    continue
                                end
                                [poly_split] = extdom_polygon(bnde,pts,1,1);
                                for ii = 1 : length(poly_split)
                                    points = unique(poly_split{ii},'rows','stable');
                                   
                                    [tmp_pfix,tmp_egfix,~]=mesh1d(points,obj.fh,obj.h0./111e3,[],obj.boubox,box_num,[]);
                                    if length(tmp_pfix) > 2
                                        [tmp_pfix,tmp_egfix]=fixgeo2(tmp_pfix,tmp_egfix);
                                    else
                                        tmp_egfix=[(1:length(tmp_pfix)-1)',(2:length(tmp_pfix))'];
                                    end
                                    tpfix = [tpfix; tmp_pfix];
                                    if isempty(tegfix)
                                        tegfix = [tegfix; tmp_egfix];
                                    else
                                        tegfix = [tegfix; tmp_egfix + max(tegfix(:))];
                                    end
                                    
                                end
                            end
                        end
                    end
                end
            end
            % remove nan 
            % Now check the full repaired geometry field
            [tpfix,tegfix]=fixgeo2(tpfix,tegfix);
            
            checkfixp = setdiff(tpfix,fixmesh(tpfix),'rows');
            if ~isempty(checkfixp)
                error('Duplicate fixed points detected, cannot proceed');
            end
            % append any user-passed fixed points.  
           if ~isempty(obj.egfix)
                obj.pfix = [obj.pfix; tpfix];
                obj.egfix = [obj.egfix;tegfix+ max(obj.egfix(:))]; % NB ADD OFFSET TO TEGFIX IF OBJ.EGFIX IS NOT EMPTY
            else
                obj.pfix = [obj.pfix; tpfix];
                obj.egfix = [obj.egfix;tegfix];
            end
            nfix = length(obj.pfix); negfix = length(obj.egfix);
            if nfix >= 0, disp(['Using ',num2str(nfix),' fixed points.']);end
            if negfix > 0
                if max(obj.egfix(:)) > length(obj.pfix)
                    error('FATAL: egfix does index correcty into pfix.');
                end
                disp(['Using ',num2str(negfix),' fixed edges.']);
                warning('Please check fixed constraints: plot(mshopts)')
            end
        end
        
        % Creates Approximate nearest neighbor objects on start-up
        function  obj = createANN(obj)
            box_vec = 1:length(obj.bbox);
            for box_num = box_vec
                if ~iscell(obj.outer)
                    dataset = obj.outer;
                    dataset(isnan(obj.outer(:,1)),:) = [];
                else
                    dataset = obj.outer{box_num};
                    dataset(isnan(obj.outer{box_num}(:,1)),:) = [];
                end
                if all(abs(obj.bbox{box_num}(1,:)) == 180)
                    % This line removes the line that can appear in the
                    % center for a global mesh
                    dataset(abs(dataset(:,1)) > 180-1e-6,:) = [];
                    dataset(abs(dataset(:,1)) < 1e-6,:) = [];
                end
                [dataset(:,1),dataset(:,2)] = m_ll2xy(dataset(:,1),dataset(:,2));
                dataset(isnan(dataset(:,1)),:) = [];
                dmy = ann(dataset');
                obj.anno{box_num} = dmy;
                obj.annData{box_num}=dataset;
            end
        end
        
        function  obj = build(obj)
            % 2-D Mesh Generator using Distance Functions.
            % Checking existence of major inputs
            %%
            warning('off','all')
            %%
            tic
            it = 1 ;
            Re = 6378.137e3;
            geps = 1e-8*min(obj.h0)/Re;
            deps = sqrt(eps);
            ttol=0.1; Fscale = 1.10; deltat = 0.1;
            delIT = 0 ;
            imp = 999; % number of iterations to do mesh improvements (delete/add)
            EXIT_QUALITY = 0.30; % minimum quality necessary to terminate if iter < itmax
            
            % unpack initial points.
            p = obj.grd.p;
            if isempty(p)
                disp('Forming initial point distribution...');
                % loop over number of boxes
                p0_boundary = []; % strip of points along all boundaries
                for box_num = 1:length(obj.h0)
                    disp(['    for box #' num2str(box_num)]);
                    % checking if cell or not and applying local values
                    h0_l = obj.h0(box_num);
                    max_r0 = 1/h0_l^2;
                    if ~iscell(obj.bbox)
                        bbox_l = obj.bbox'; % <--we must tranpose this!
                    else
                        bbox_l = obj.bbox{box_num}'; % <--tranpose!
                    end
                    if ~iscell(obj.fh)
                        fh_l = obj.fh;
                    else
                        fh_l = obj.fh{box_num};
                    end
                    % Lets estimate the num_points the distribution will be
                    num_points = ceil(2/sqrt(3)*prod(abs(diff(bbox_l)))...
                        /(h0_l/111e3)^2);
                    noblks = ceil(num_points*2*8/obj.memory_gb*1e-9);
                    len = abs(bbox_l(1,1)-bbox_l(2,1));
                    blklen = len/noblks;
                    st = bbox_l(1,1) ; ed = st + blklen; ns = 1;
                    %% 1. Create initial distribution in bounding box
                    %% (equilateral triangles)
                    for blk = 1 : noblks
                        if blk == noblks
                            ed = bbox_l(2,1);
                        end
                        ys = bbox_l(1,2);
                        ny = floor(1e3*m_lldist(repmat(0.5*(st+ed),2,1),...
                            [ys;bbox_l(2,2)])/h0_l);
                        dy = diff(bbox_l(:,2))/ny;
                        ns = 1;
                        % start at lower left and make grid going up to
                        % north latitude
                        for ii = 1:ny+1
                            if st*ed < 0
                                nx = floor(1e3*m_lldist([st;0],...
                                    [ys;ys])/(2/sqrt(3)*h0_l)) + ...
                                    floor(1e3*m_lldist([0;ed],...
                                    [ys;ys])/(2/sqrt(3)*h0_l));
                            else
                                nx = floor(1e3*m_lldist([st;ed],...
                                    [ys;ys])/(2/sqrt(3)*h0_l));
                            end
                            ne = ns+nx-1;
                            if mod(ii,2) == 0
                                % no offset
                                x(ns:ne) = linspace(st,ed,nx);
                            else
                                % offset
                                dx = (ed-st)/nx;
                                x(ns:ne) = linspace(st+0.5*dx,ed,nx);
                            end
                            % tolerance
                            if ii == (ny + 1)
                                y(ns:ne) = ys - eps;
                            else
                                y(ns:ne) = ys;
                            end
                            ns = ne+1; ys = ys + dy;
                            
                        end
                        
                        st = ed;
                        ed = st + blklen;
                        p1 = [x(:) y(:)]; clear x y
                        
                        %% 2. Remove points outside the region, apply the rejection method
                        p1 = p1(feval(obj.fd,p1,obj,box_num) < geps,:);     % Keep only d<0 points
                        r0 = 1./feval(fh_l,p1).^2;                          % Probability to keep point
                        p1 = p1(rand(size(p1,1),1) < r0/max_r0,:);          % Rejection method
                        p  = [p; p1];                                       % Adding p1 to p
                    end
                    % for now restrict to domains with a single box
                    if box_num == 1 
                      % convert to pslg 
                      [nodes, edges] = getnan2(obj.outer{box_num});
                      % find the largest polygon
                      [poly,~,~,~] = extdom_polygon(edges,nodes,1);
                      for k = 1 : length(poly)
                          [tmp_p0_boundary, layer1]=form_initial_points_near_and_along_boundary(poly{k}, ...
                              obj.boubox{box_num}, ...
                              obj.h0(box_num)/2, ...
                              obj.fh);
                          % remove all p in layer/strip along the boundary
                          [in_layer1,~]=inpoly(p,layer1);
                          remove_percent = sum(in_layer1)/length(in_layer1);
                          % if many points are to be removed consider the
                          % inverse
                          if remove_percent > 0.90
                              p = p(in_layer1,:);
                          else
                              p = p(~in_layer1,:);
                          end
                   
                          p0_boundary = [p0_boundary; tmp_p0_boundary];
                      end
                    end

                end
            else
                disp('User-supplied initial points!');
                obj.grd.b = [];
                h0_l = obj.h0(end); % finest h0 (in case of a restart of meshgen.build).
            end

            is_within = feval(obj.fd,p0_boundary,obj,[],1) <= 0; 
            p0_boundary = p0_boundary(is_within,:);

            % prune these p0_boundary point if they are too close to one
            % another 
            [~, dis_to_p0] = ourKNNsearch(p0_boundary', p0_boundary',2);

            % compute local resolution nearby these points
            local_resolution = zeros(length(p0_boundary),1);
            for box_num = 1:length(obj.h0)                             
                if ~iscell(obj.fh)                                    
                    fh_l = obj.fh;
                else
                    fh_l = obj.fh{box_num};
                end
                if box_num > 1
                    iboubox = obj.boubox{box_num}(1:end-1,:) ;
                    inside = inpoly(p0_boundary,iboubox) ;
                else
                    inside = true(size(p0_boundary,1),1);
                end
                local_resolution(inside) = feval(fh_l,p0_boundary(inside,:));      % Ideal lengths
            end
            % if these points are exceedingly close, remove them as it 
            % would otherwise lead to a poor quality element. 
            prune1 = dis_to_p0(:,2) < local_resolution./111e3./2; 

            if ~isempty(obj.pfix) 
                 % drop all p0_boundary nearby pfix 
                 % giving preference for pfix
                 [~, dis_to_p0] = ourKNNsearch(obj.pfix',p0_boundary',3);
                 prune2 = dis_to_p0(:,1) < 5*local_resolution./111e3; 
                 prune3 = dis_to_p0(:,2) < 5*local_resolution./111e3; 
                 prune4 = dis_to_p0(:,3) < 5*local_resolution./111e3;
                 prune = prune1 | prune2 | prune3 | prune4; 
            else 
                prune = prune1;
            end 
            
            p0_boundary(prune,:)=[]; 
  
            num_of_locked_boundary_nodes = length(p0_boundary);

            nfix = length(obj.pfix); negfix = length(obj.egfix);
            % first pfix associated, then p0_boundary, then regular p
            % nfix+1 : nfix+num_of_locked_boundary_nodes+1;
            if ~isempty(obj.pfix)
                p = [obj.pfix; p0_boundary; p];
            else 
                p = [p0_boundary; p];                  
            end
            N = size(p,1); % Number of points N
            disp(['Number of initial points after rejection is ',num2str(N)]);
            disp(['Number of locked boundary points is ',num2str(num_of_locked_boundary_nodes)]);

            %% Iterate
            pold = inf;                                                    % For first iteration
            if obj.plot_on >= 1
                clf,view(2),axis equal;
            end
            toc
            fprintf(1,' ------------------------------------------------------->\n') ;
            disp('Begin iterating...');
            while 1
                tic
                if ~mod(it,obj.nscreen) && delIT == 0
                    disp(['Iteration =' num2str(it)]) ;
                end
                % 3. Retriangulation by the Delaunay algorithm
                if max(sqrt(sum((p(1:size(pold,1),:)-pold).^2,2))/h0_l*111e3) > ttol         % Any large movement?
                    if it > 1
                        p = fixmesh([obj.pfix; p]); 
                    else
                        p = fixmesh(p);                                        % Ensure only unique points.
                    end
                    N = size(p,1); pold = p;                               % Save current positions
                    [t,p] = delaunay_elim(p,obj.fd,geps,0);                % Delaunay with elimination
                    if isempty(t)
                        disp('Exiting')
                        return
                    end
                    % Getting element quality and check "goodness"
                    if exist('pt','var'); clear pt; end
                    [pt(:,1),pt(:,2)] = m_ll2xy(p(:,1),p(:,2));
                    tq = gettrimeshquan( pt, t);
                    mq_m = mean(tq.qm);
                    mq_l = min(tq.qm);
                    mq_s = std(tq.qm);
                    mq_l3sig = mq_m - 3*mq_s;
                    obj.qual(it,:) = [mq_m,mq_l3sig,mq_l];
                    % If mesh quality went down "significantly" since last iteration
                    % which was a mesh improvement iteration, then rewind.
                    if ~mod(it,imp+1) && obj.qual(it,1) - obj.qual(it-1,1) < -0.10 || ...
                            ~mod(it,imp+1) && (N - length(p_before_improve))/length(p_before_improve) < -0.10
                        disp('Mesh improvement was unsuccessful...rewinding...');
                        p = p_before_improve;
                        N = size(p,1);                                     % Number of points changed
                        pold = inf;
                        it = it + 1;
                        continue
                    else
                        N = size(p,1); pold = p;                           % Assign number of points and save current positions
                    end
                    % 4. Describe each bar by a unique pair of nodes.
                    bars = [t(:,[1,2]); t(:,[1,3]); t(:,[2,3])];           % Interior bars duplicated
                    bars = unique(sort(bars,2),'rows');                    % Bars as node pairs
                    % 5. Graphical output of the current mesh
                    if obj.plot_on >= 1 && (mod(it,obj.nscreen)==0 || it == 1)
                        cla,m_triplot(p(:,1),p(:,2),t)
                        m_grid
                        title(['Iteration = ',num2str(it)]);
                        if negfix > 0
                            m_plot(reshape(obj.pfix(obj.egfix,1),[],2)',...
                                reshape(obj.pfix(obj.egfix,2),[],2)','r-')
                        end
                        if nfix > 0
                            m_plot(obj.pfix(:,1),obj.pfix(:,2),'b.')
                        end
                        m_plot(p0_boundary(:,1),p0_boundary(:,2),'r.')
                        hold on ;
                        axis manual
                        %axis([-0.0032   -0.0029    0.4818    0.4820]);
                        drawnow
                    end
                end
                % Getting element quality and check goodness
                if exist('pt','var'); clear pt; end
                [pt(:,1),pt(:,2)] = m_ll2xy(p(:,1),p(:,2));
                tq = gettrimeshquan( pt, t);
                mq_m = mean(tq.qm);
                mq_l = min(tq.qm);
                mq_s = std(tq.qm);
                mq_l3sig = mq_m - 3*mq_s;
                obj.qual(it,:) = [mq_m,mq_l3sig,mq_l];
                % Termination quality, mesh quality reached is copacetic.
                qual_diff = mq_l3sig - obj.qual(max(1,it-imp),2);
                if ~mod(it,imp)
                    if mq_l > EXIT_QUALITY
                        % Do the final elimination of small connectivity
                        disp('Quality of mesh is good enough, exit')
                        close all;
                        break;
                    end
                end
                % Saving a temp mesh
                if ~mod(it,obj.nscreen) && delIT == 0
                    disp(['Number of nodes is ' num2str(length(p))])
                    disp(['Mean mesh quality is ' num2str(mq_m)])
                    disp(['Min mesh quality is ' num2str(mq_l)])
                    disp(['3rd sigma lower mesh quality is ' num2str(mq_l3sig)])
                    tempp = p; tempt = t;
                    save('Temp_grid.mat','it','tempp','tempt');
                    clearvars tempp tempt
                end
                % 6. Move mesh points based on bar lengths L and forces F
                barvec = pt(bars(:,1),:)- pt(bars(:,2),:);                 % List of bar vectors
                if strcmp(obj.grd.proj.name,'UTM')
                    % UTM is already in meters (useful for small domains)
                    L = sqrt(sum(barvec.^2,2))*Re;
                else
                    % Get spherical earth distances
                    long   = zeros(length(bars)*2,1);
                    lat    = zeros(length(bars)*2,1);
                    long(1:2:end) = p(bars(:,1),1);
                    long(2:2:end) = p(bars(:,2),1);
                    lat(1:2:end)  = p(bars(:,1),2);
                    lat(2:2:end)  = p(bars(:,2),2);
                    L = m_lldist(long,lat); L = L(1:2:end)*1e3;            % L = Bar lengths in meters
                end
                ideal_bars = 0.5*(pt(bars(:,1),:) + pt(bars(:,2),:));      % Used to determine what bars are in bbox
                [ideal_bars(:,1),ideal_bars(:,2)] = ...                    % needs to be in non-projected
                    m_xy2ll(ideal_bars(:,1),ideal_bars(:,2));              % coordinates
                hbars = 0*ideal_bars(:,1);
                
                for box_num = 1:length(obj.h0)                             % For each bbox, find the bars that are in and calculate
                    if ~iscell(obj.fh)                                     % their ideal lengths.
                        fh_l = obj.fh;
                    else
                        fh_l = obj.fh{box_num};
                    end
                    h0_l = obj.h0(box_num);
                    if box_num > 1
                        h0_l = h0_l/111e3;                                 % create buffer to evalulate fh between nests
                        iboubox = obj.boubox{box_num}(1:end-1,:) ;
                        inside = inpoly(ideal_bars,iboubox) ;
                    else
                        inside = true(size(hbars));
                    end
                    hbars(inside) = feval(fh_l,ideal_bars(inside,:));      % Ideal lengths
                end
                
                L0 = hbars*Fscale*median(L)/median(hbars);                 % L0 = Desired lengths using ratio of medians scale factor
                LN = L./L0;                                                % LN = Normalized bar lengths
                
                % Mesh improvements (deleting and addition)
                p_before_improve = p;
                if ~mod(it,imp) % && ~HIGH_FIDELITY_MODE
                    nn = []; pst = [];
                    if qual_diff < imp*obj.qual_tol && qual_diff > 0
                        % Remove elements with small connectivity
                        nn = get_small_connectivity(p,t);
                        disp(['Deleting ' num2str(length(nn)) ' due to small connectivity'])
                        
                        % Remove points that are too close (< LN = 0.5)
                        if any(LN < 0.5)
                            % Do not delete pfix too close.
                            nn1 = setdiff(reshape(bars(LN < 0.5,:),[],1),[(1:nfix)']);
                            disp(['Deleting ' num2str(length(nn1)) ' points too close together'])
                            nn = unique([nn; nn1]);
                        end
                        
                        % Split long edges however many times to
                        % better lead to LN of 1
                        if any(LN > 2)
                            nsplit = floor(LN);
                            nsplit(nsplit < 1) = 1;
                            adding = 0;
                            % Split once
                            for jj = 2:2
                                il = find(nsplit >= jj);
                                xadd = zeros(length(il),jj-1);
                                yadd = zeros(length(il),jj-1);
                                for jjj = 1 : length(il)
                                    deltax = (p(bars(il(jjj),2),1)- p(bars(il(jjj),1),1))/jj;
                                    deltay = (p(bars(il(jjj),2),2)- p(bars(il(jjj),1),2))/jj;
                                    xadd(jjj,:) = p(bars(il(jjj),1),1) + (1:jj-1)*deltax;
                                    yadd(jjj,:) = p(bars(il(jjj),1),2) + (1:jj-1)*deltay;
                                end
                                pst = [pst; xadd(:) yadd(:)];
                                adding = numel(xadd) + adding;
                            end
                            disp(['Adding ',num2str(adding) ,' points.'])
                        end
                    end
                    if ~isempty(nn) || ~isempty(pst)
                        % Doing the actual subtracting and add
                        p(nn,:)= [];
                        p = [p; pst];
                        pold = inf;
                        it = it + 1;
                        continue;
                    end
                end
                
                F    = (1-LN.^4).*exp(-LN.^4)./LN;                         % Bessens-Heckbert edge force
                Fvec = F*[1,1].*barvec;
                
                Ftot = full(sparse(bars(:,[1,1,2,2]),ones(size(F))*[1,2,1,2],[Fvec,-Fvec],N,2));
                Ftot(1:nfix,:) = 0;                                        % Force = 0 at fixed points
                Ftot(nfix+1:nfix+1+num_of_locked_boundary_nodes,:) = 0;    % Force = 0 at boundary points.
                pt = pt + deltat*Ftot;                                     % Update node positions
                
                [p(:,1),p(:,2)] = m_xy2ll(pt(:,1),pt(:,2));
                
                %7. Bring outside points back to the boundary
                d = feval(obj.fd,p,obj,[],1); ix = d > 0;                  % Find points outside (d>0)
                ix(1:nfix) = 0;
                alpha = 1.0;
                for ib = 1 : 1  %length(obj.improve_boundary)
                    if sum(ix) > 0
                        pn = p(ix,:) + deps;
                        dgradx = (feval(obj.fd,[pn(:,1),p(ix,2)],obj,[])...%,1)...
                            -d(ix))/deps; % Numerical
                        dgrady = (feval(obj.fd,[p(ix,1),pn(:,2)],obj,[])...%,1)...
                            -d(ix))/deps; % gradient
                        dgrad2 = dgradx.^+2 + dgrady.^+2;
                        dgrad2(dgrad2 < eps) = eps;
                        p(ix,:) = p(ix,:) - alpha*[d(ix).*dgradx./dgrad2,...
                            d(ix).*dgrady./dgrad2];
                    end
                    alpha = alpha / 0.5;
                end
                
                % 8. Termination criterion: Exceed itmax
                it = it + 1 ;
                
                if ( it > obj.itmax )
                    % Do the final deletion of small connectivity
                    disp('too many iterations, exit')
                    close all;
                    break ;
                end
                toc
            end
            %%
            warning('on','all')
            %%
            disp('Finished iterating...');
            fprintf(1,' ------------------------------------------------------->\n') ;
            
            %% Doing the final cleaning and fixing to the mesh...
            % Always save the mesh!
            save('Precleaned_grid.mat','it','p','t');
            
            % Clean up the mesh if specified
            if ~strcmp(obj.cleanup,'none')
                % Put the mesh class into the grd part of meshgen and clean
                obj.grd.p = p; obj.grd.t = t;
                [obj.grd,qout] = clean(obj.grd,obj.cleanup,...
                    'nscreen',obj.nscreen,'djc',obj.dj_cutoff,...
                    'pfix',obj.pfix);
                obj.grd.pfix = obj.pfix ;
                obj.grd.egfix= obj.egfix ;
                obj.qual(end+1,:) = qout;
            else
                % Fix mesh on the projected space
                [p(:,1),p(:,2)] = m_ll2xy(p(:,1),p(:,2));
                [p,t] = fixmesh(p,t);
                [p(:,1),p(:,2)] = m_xy2ll(p(:,1),p(:,2));
                % Put the mesh class into the grd part of meshgen
                obj.grd.p = p; obj.grd.t = t;
                obj.grd.pfix = obj.pfix ;
                obj.grd.egfix= obj.egfix ;
            end
            
            % Check element order, important for the global meshes crossing
            % -180/180 boundary
            obj.grd = CheckElementOrder(obj.grd);
            
            if obj.plot_on
                figure; plot(obj.qual,'linewi',2);
                hold on
                % plot the line dividing cleanup and distmesh
                plot([it it],[0 1],'--k')
                xticks(1:5:obj.itmax);
                xlabel('Iterations'); ylabel('Geometric element quality');
                title('Geometric element quality with iterations');
                set(gca,'FontSize',14);
                legend('q_{mean}','q_{mean}-q_{3\sigma}', 'q_{min}','Location','best');
                grid minor
            end
            return;
            %%%%%%%%%%%%%%%%%%%%%%%%%%
            % Auxiliary subfunctions %
            %%%%%%%%%%%%%%%%%%%%%%%%%%
            
            function [t,p] = delaunay_elim(p,fd,geps,final)
                % Removing mean to reduce the magnitude of the points to
                % help the convex calc
                if exist('pt1','var'); clear pt1; end
                [pt1(:,1),pt1(:,2)] = m_ll2xy(p(:,1),p(:,2));
                if isempty(obj.egfix)
                    p_s  = pt1 - repmat(mean(pt1),[N,1]);
                    TR   = delaunayTriangulation(p_s);
                else
                    TR   = delaunayTriangulation(pt1(:,1),pt1(:,2),obj.egfix);
                    pt1  = TR.Points;
                end
                for kk = 1:final+1
                    if kk > 1
                        % Perform the following below upon exit from the mesh
                        % generation algorithm
                        nn = get_small_connectivity(pt1,t);
                        nn1 = [];
                        nn = unique([nn; nn1]) ;
                        TR.Points(nn,:) = [];
                        pt1(nn,:) = [];
                    end
                    t = TR.ConnectivityList;
                    pmid = squeeze(mean(reshape(pt1(t,:),[],3,2),2));      % Compute centroids
                    [pmid(:,1),pmid(:,2)] = m_xy2ll(pmid(:,1),pmid(:,2));  % Change back to lat lon
                    t    = t(feval(fd,pmid,obj,[]) < -geps,:);             % Keep interior trianglesi
                end
                if length(pt1) ~= length(p)
                    clear p
                    [p(:,1),p(:,2)] = m_xy2ll(pt1(:,1),pt1(:,2));
                end
            end
            
            function nn = get_small_connectivity(p,t)
                % Get node connectivity (look for 4)
                [~, enum] = VertToEle(t);
                % Make sure they are not boundary nodes
                bdbars = extdom_edges2(t, p);
                bdnodes = unique(bdbars(:));
                I = find(enum <= 4);
                nn = setdiff(I',[(1:nfix)';bdnodes]);                      % and don't destroy pfix or egfix!
                return;
            end
            
        end % end mesh generator
        
    end % end methods
    
end % end class

and the vectorized function to form those points

function [initial_points, layer_of_points] = form_initial_points_near_and_along_boundary(outer_polygon, ...
    bounding_polygon, ...
    minimum_resolution_in_meters, ...
    mesh_sizing_fx)
% Form a layer of points one row along the boundary rim of the domain
% so that when triangulated it forms perfectly equilateral triangles. 

minimum_resolution_in_degrees = minimum_resolution_in_meters / 111e3; 

% make sure the polygon is unique 

outer_polygon = unique(outer_polygon,'rows','stable');

%% Resample the outer polygon to match the point spacing following mesh_sizing_fx 
[outer_polygon_rs,~,~]=mesh1d(outer_polygon, ...
                                                   mesh_sizing_fx, ...
                                                   minimum_resolution_in_degrees, ...
                                                   [], ...
                                                   bounding_polygon, ...
                                                   1, ... % box_number (assumed 1). 
                                                   [] ...
                                                   );
layer_of_points = [];

%figure;

% for i = 1 : length(outer_polygon_rs)-1
%     % 1) Compute the midpoint
%     p1 = outer_polygon_rs(i,:);
%     p2 = outer_polygon_rs(i+1,:);
%     midpoint = (p1 + p2) / 2.0;
%     
%     % 2) Compute the tangent and then the unit normal vector pointing inwards
%     tangent = p2 - p1; % Tangent vector
%     tangent_length = norm(tangent); % Length of the tangent vector
%     unit_tangent = tangent / tangent_length; % Normalize to get unit tangent
%     unit_normal = [-unit_tangent(2), unit_tangent(1)]; % Rotate by 90 degrees to get the unit normal
%     
%     % 3) Calculate the distance for the equilateral triangle and compute p3
%     local_size_in_degrees = mesh_sizing_fx{1}(midpoint) / 111e3; 
%     % The distance from the midpoint to the third point to form an equilateral triangle
%     % Assuming local_size_in_degrees gives the desired edge length of the equilateral triangle
%     height_of_triangle = (sqrt(3) / 2) * local_size_in_degrees; % Height of an equilateral triangle
%     p3 = midpoint + height_of_triangle * unit_normal; % New point forming equilateral triangle
% 
%     % For debugging plot
% %     hold on; 
% %     plot([p1(1), p2(1)], [p1(2), p2(2)], 'k-'); % Plot edge
% %     plot(midpoint(1), midpoint(2), 'k.'); % Plot midpoint
% %     quiver(midpoint(1), midpoint(2), unit_normal(1), unit_normal(2), 0.05, 'r'); % Plot unit normal
% %     plot(p3(1), p3(2), 'bs'); % Plot new point
% 
%     layer_of_points = [layer_of_points; p3];
% end

% Assuming outer_polygon_rs is an Nx2 matrix of [x, y] coordinates
% Calculate midpoints for all segments
p1 = outer_polygon_rs(1:end-1, :);
p2 = outer_polygon_rs(2:end, :);
midpoints = (p1 + p2) / 2.0;

% Calculate tangent vectors and their lengths
tangents = p2 - p1;
tangent_lengths = sqrt(sum(tangents.^2, 2));

% Normalize tangents to get unit tangents and then compute unit normals
unit_tangents = tangents ./ tangent_lengths;
unit_normals = -1*[-unit_tangents(:,2), unit_tangents(:,1)]; % Rotate 90 degrees for normals

% Calculate distances for equilateral triangles and compute p3 for all segments
local_size_in_degrees = mesh_sizing_fx{1}(midpoints) ./ 111e3; % Vectorized call to the sizing function
height_of_triangle = (sqrt(3) / 2) * local_size_in_degrees; % Heights for all triangles

% Apply heights to normals to get new points
p3 = midpoints + height_of_triangle .* unit_normals;

% Optional: Plotting (assuming you want to plot all at once, adjust for large datasets)
% hold on;
% plot(outer_polygon_rs(:,1), outer_polygon_rs(:,2), 'k-'); % Plot polygon edges
% plot(midpoints(:,1), midpoints(:,2), 'k.'); % Plot all midpoints
% quiver(midpoints(:,1), midpoints(:,2), unit_normals(:,1), unit_normals(:,2), 0.05, 'r'); % Plot normals
% plot(p3(:,1), p3(:,2), 'bs'); % Plot new points (p3)

% Append new points to layer_of_points
layer_of_points = [layer_of_points; p3];

% Combine initial resampled outer_polygon with layer_of_points 
initial_points = [outer_polygon_rs; layer_of_points];

end