adding high-fidelity option to meshgen for meshing hardened shorelines/coastal structures

@geodata/geodata.m

@@ -17,6 +17,7 @@
     %     addOptional(p,'pslg',defval);
     %     addOptional(p,'boubox',defval);
     %     addOptional(p,'window',defval);
+    %     addOptional(p,'high_fidelity',defval);
     %
     %   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
@@ -58,6 +59,8 @@
         spacing = 2.0 ; %Relative spacing along polygon, large effect on computational efficiency of signed distance.
         gridspace
         shapefile_3d % if the shapefile has a height attribute
+        high_fidelity % performs a 1D mesh generation step to form pfix and egfix prior to 2D meshing
+
     end
 
     methods
@@ -120,6 +123,8 @@
             addOptional(p,'boubox',defval);
             addOptional(p,'window',defval);
             addOptional(p,'shapefile_3d',defval);
+            addOptional(p,'high_fidelity',defval);
+
 
             % parse the inputs
             parse(p,varargin{:});
@@ -204,7 +209,9 @@
                             obj.window = 5;
                         end
                     case('shapefile_3d')
-                        obj.shapefile_3d = inp.(fields{i}) ;
+                         obj.shapefile_3d = inp.(fields{i}) ;
+                    case('high_fidelity')
+                        obj.high_fidelity = inp.(fields{i});
                     case('weirs')
                         if ~iscell(inp.(fields{i})) && ~isstruct(inp.(fields{i})) && inp.(fields{i})==0, continue; end
                         if ~iscell(inp.(fields{i})) && ~isstruct(inp.(fields{i}))
@@ -861,4 +868,3 @@ function plot(obj,type,projection,holdon)
     end
 
 end
-

@geodata/private/Read_shapefile.m

@@ -287,4 +287,4 @@
     end
 end
 %EOF
-end
+end

@meshgen/private/Get_line_edges.m

@@ -0,0 +1,17 @@
+function edges = Get_line_edges(nodes)
+%GET_LINE_EDGES Generate edges from a sequence of nodes in a closed loop.
+%   edges = GET_LINE_EDGES(nodes) takes a list of nodes (vertices) as input
+%   and generates a list of edges where each node is connected to the next,
+%   forming a closed loop. The output 'edges' is an Mx2 matrix, where M is
+%   the number of edges, and each row represents an edge defined by the
+%   indices of its two endpoints.
+
+    % Number of nodes
+    nNodes = length(nodes);
+
+    % Generate edges connecting each node to the next
+    edges = [(1:nNodes-1)', (2:nNodes)'];
+
+    % Add the edge connecting the last node back to the first
+    edges(end+1, :) = [nNodes, 1];
+end

@meshgen/private/filter_polygon_constraints.m

@@ -0,0 +1,46 @@
+function [pfix, egfix] = filter_polygon_constraints(pfix, egfix, ibouboxes, box_number)
+%FILTER_CONSTRAINTS Removes edge constraints based on bounding box criteria.
+%   This function filters out edge constraints (egfix) where at least one
+%   endpoint (from pfix) is outside the specified bounding box (ibouboxes)
+%   for the given box_number. It also adjusts the edges based on nested
+%   bounding boxes, if applicable.
+% remove bars if one point is outside
+node1=pfix(egfix(:,1),:);
+node2=pfix(egfix(:,2),:);
+iboubox = ibouboxes{box_number};
+% to enlarge or shrink the box, you must make sure bbox is equi
+% spaced
+[ty,tx]=my_interpm(iboubox(:,2),iboubox(:,1),100/111e3);
+iboubox = [tx,ty]; 
+buffer_size = 1.0;
+iboubox(:,1) = buffer_size*iboubox(:,1)+(1-buffer_size)*mean(iboubox(1:end-1,1));
+iboubox(:,2) = buffer_size*iboubox(:,2)+(1-buffer_size)*mean(iboubox(1:end-1,2));
+inside_node1 = inpoly(node1,iboubox(1:end-1,:)) ;
+inside_node2 = inpoly(node2,iboubox(1:end-1,:)) ;
+inside = inside_node1 .* inside_node2;
+% Get all points inside inner boxes and consider these outside for
+% all the nested boxes.
+for bn = box_number+1:length(ibouboxes)
+    % enlarge nest
+    iboubox = ibouboxes{bn}(1:end-1,:);
+    [ty,tx]=my_interpm(iboubox(:,2),iboubox(:,1),100/111e3);
+    iboubox = [tx,ty];
+    iboubox(:,1) = 1.25*iboubox(:,1)+(1-1.25)*mean(iboubox(1:end-1,1));
+    iboubox(:,2) = 1.25*iboubox(:,2)+(1-1.25)*mean(iboubox(1:end-1,2));
+    inside_node1 = inpoly(node1,iboubox);
+    inside_node2 = inpoly(node2,iboubox);
+    inside2 = inside_node1 .* inside_node2;
+    inside(find(inside2)) = false;
+end
+egfix(~inside,:) = [];
+tegfix=egfix';
+uid = unique(tegfix(:));
+tuid = length(uid);
+% remove nfix outside iboubox
+if tuid > 0
+    % remove pfix if outside domain
+    pfix = pfix(uid,:);
+end
+egfix = renumberEdges(egfix);
+
+end

@meshgen/private/mesh1d.m

@@ -0,0 +1,207 @@
+
+function [pout,t,converged]=mesh1d(poly,fh0,h,fix,boubox,box_num0,varargin)
+% Mesh Generator using Distance Functions.
+%   [P,T]=mesh1d(FDIST,FH,H,BOX,FIX,FDISTPARAMS)
+%
+%   POUT:           Resampled Node positions (N X 2)
+%      T:           Edge indices (NTx2))
+%   POLY:           The polygon (NX X 2
+%      FH:          Edge length function
+%      H:           Smallest edge length
+%      FIX:         Indices of Poly Fixed node positions (NFIX X 1)
+
+% Preprocessing steps - move into parametric space
+%poly = add_flairs(poly, 50); 
+
+X = poly(:,1);
+Y = poly(:,2);
+
+% find sharp corners & add flairs 
+% [sharp_corners] = find_sharp_corners(poly, 40);
+% if ~isempty(sharp_corners)
+%     fix = [fix; sharp_corners]; 
+% end
+
+
+% add ending point to emulate periodic mesh generation
+%X(end+1) = X(1,:) + eps;
+%Y(end+1) = Y(1,:) + eps;
+
+xd = diff(X);
+yd = diff(Y);
+dist = sqrt(xd.^2 + yd.^2);
+u = cumsum(dist);
+u = [0; u];
+totaldist = u(end);
+
+np = totaldist / min(h);
+% kjr check if total np will be too little 
+% i.e., less than 5 
+if np < 5
+    pout = [];
+    t = [];
+    converged=1; 
+    return 
+end 
+t = linspace(0,max(u),ceil(np));
+
+xn = interp1(u,X,t);
+yn = interp1(u,Y,t);
+
+u_fix = [];
+if ~isempty(fix) 
+    % Convert fixed points to parametric space
+    % determine the distance these are points are from the zero point. 
+    u_fix = [];
+    for i = 1:length(fix)
+        xd_tmp = diff(X(1:fix(i)));
+        yd_tmp = diff(Y(1:fix(i)));
+        dist_tmp = sqrt(xd_tmp.^2 + yd_tmp.^2);
+        if isempty(dist_tmp)
+            continue
+        elseif cumsum(dist_tmp) == totaldist
+            continue
+        else
+            csum = cumsum(dist_tmp);
+        end
+        u_fix = [u_fix; csum(end)];
+    end
+end
+
+A = 0;
+B = totaldist;
+
+box  = [A,B]';
+
+fdist = @(x) my_1d_sdf(x, A, B);
+
+tmph = nestedHFx([xn',yn'],boubox,fh0,h);
+
+    function [h_inner] = nestedHFx(x,boubox,fh0,h)
+        h_inner= x(:,1)*0;
+        for box_num = 1:length(h)  % For each bbox, find the points that are in and calculate
+            fh_l = fh0{box_num};
+            if box_num > 1
+                iboubox = boubox{box_num}(1:end-1,:) ;
+                inside = inpoly(x,iboubox) ;
+            else
+                inside = true(size(h_inner));
+            end
+            h_inner(inside) = feval(fh_l,x(inside,:))./111e3; 
+        end
+    end
+
+F = griddedInterpolant(t,tmph,'linear','linear');
+
+fh = @(x) F(x);
+
+dim=size(box,2);
+ptol=.01;
+L0mult=1+.4/2^(dim-1);
+deltat=.10; geps=1e-3*h(box_num0);
+deps=sqrt(eps)*h(box_num0);
+
+% 1. Create initial distribution in bounding box
+if dim==1
+    p=(box(1):h(box_num0):box(2))';
+else
+    disp('Only 1D is supported')
+    %error
+end
+
+% 2. Remove points outside the region, apply the rejection method
+p=p(feval(fdist,p)<geps,:);
+r0=feval(fh,p);
+% Add the fixed points.
+if ~isempty(u_fix)
+    p=[u_fix; p(rand(size(p,1),1)<min(r0)^dim./r0.^dim,:)];
+else 
+    p = p(rand(size(p,1),1)<min(r0)^dim./r0.^dim,:);
+end
+% find their indicies
+p = unique(p,'stable'); 
+[p,I] = sort(p);
+[~,fixed_point_ind] = sort(I);
+fixed_point_ind = fixed_point_ind(2:length(u_fix)-1);
+% find the position of the fixed points  
+N=size(p,1);
+
+count=0;
+while 1
+    % 3. Retriangulation
+    t=[(1:length(p)-1)',(2:length(p))'];
+    pmid=zeros(size(t,1),dim);
+    for ii=1:dim+1
+        pmid=pmid+p(t(:,ii),:)/(dim+1);
+    end
+    t=t(feval(fdist,pmid)<-geps,:);
+    % 4. Describe each edge by a unique pair of nodes
+    pair=zeros(0,2);
+    localpairs=nchoosek(1:dim+1,2);
+    for ii=1:size(localpairs,1)
+        pair=[pair;t(:,localpairs(ii,:))];
+    end
+    pair=unique(sort(pair,2),'rows');
+    % 5. Graphical output of the current mesh
+    %fprintf('Retriangulation #%d\n',count)
+    count = count + 1;
+    % 6. Move mesh points based on edge lengths L and forces F
+    bars=p(pair(:,1),:)-p(pair(:,2),:);
+    L=sqrt(sum(bars.^2,2));
+    L0=feval(fh,(p(pair(:,1),:)+p(pair(:,2),:))/2);
+    L0=L0*L0mult*(sum(L.^dim)/sum(L0.^dim))^(1/dim);
+    F=max(L0-L,0);
+    Fbar=[bars,-bars].*repmat(F./L,1,2*dim);
+    dp=full(sparse(pair(:,[ones(1,dim),2*ones(1,dim)]), ...
+        ones(size(pair,1),1)*[1:dim,1:dim], ...
+        Fbar,N,dim));
+    %dp(1:size(fix,1),:)=0;
+    dp(fixed_point_ind,:)=0;
+    % lock the first and last point 
+    dp(1) = 0; dp(end) = 0;
+
+    p=p+deltat*dp;
+
+    % 7. Bring outside points back to the boundary
+    d=feval(fdist,p); 
+    ix=d>0;
+    %ix(1:length(fix),:) = 0;
+    ix(fixed_point_ind,:) = 0;
+    gradd=zeros(sum(ix),dim);
+    for ii=1:dim
+        a=zeros(1,dim);
+        a(ii)=deps;
+        d1x=feval(fdist,p(ix,:)+ones(sum(ix),1)*a);
+        gradd(:,ii)=(d1x-d(ix))/deps;
+    end
+    p(ix,:)=p(ix,:)-d(ix)*ones(1,dim).*gradd;
+
+    % 8. Termination criterion
+    maxdp=max(deltat*sqrt(sum(dp(d<-geps,:).^2,2)));
+
+    %disp(maxdp)
+    if maxdp<ptol*h
+        %disp('Converged...')
+        converged=1;
+        % put back in the real space
+        pout(:,1) = interp1(u,X,p);
+        pout(:,2) = interp1(u,Y,p);
+        %figure; plot(pout(:,1),pout(:,2),'rs')
+        if any(isnan(pout(:,1)))
+            warning('Nans detected in 1d mesh')
+        end
+        break
+    end
+    if count > 5000
+        warning('Mesh1d: some line segments did not converge...')
+        converged=0;
+        pout(:,1) = interp1(u,X,p);
+        pout(:,2) = interp1(u,Y,p);
+        if any(isnan(pout(:,1)))
+            warning('Nans detected in 1d mesh')
+        end
+        break
+    end
+
+end
+end

@meshgen/private/my_1d_sdf.m

@@ -0,0 +1,10 @@
+function [dist] = my_1d_sdf(p,A,B)
+distA = abs(p - A);
+distB = abs(p - B);
+dist = -min(distA,distB);
+if p < A
+    dist = -1*dist;
+elseif p > B
+    dist = -1*dist;
+end
+end

@meshgen/private/nandelim_to_cell.m

@@ -0,0 +1,28 @@
+function polygons = nandelim_to_cell(data)
+%NANDELIM_TO_CELL Convert NaN-delimited data into a cell array of polygons.
+%   polygons = NANDELIM_TO_CELL(data) takes an Nx2 matrix 'data', where
+%   each polygon is separated by a row with NaN in the first column, and
+%   converts it into a cell array. Each cell contains the vertices of a
+%   polygon defined by consecutive rows up to the next NaN row.
+
+    % Identify the indices of rows that have NaN in the first column
+    nanRows = find(isnan(data(:, 1)));
+
+    % Preallocate the cell array for efficiency
+    polygons = cell(length(nanRows) + 1, 1);
+
+    % Set the starting index for the first polygon
+    startIndex = 1;
+
+    % Iterate over each NaN row to extract polygons
+    for i = 1:length(nanRows)
+        % Extract the polygon data up to the row before the NaN
+        polygons{i} = data(startIndex:nanRows(i) - 1, :);
+
+        % Update the start index for the next polygon
+        startIndex = nanRows(i) + 1;
+    end
+
+    % Handle the last polygon after the final NaN row
+    polygons{end} = data(startIndex:end, :);
+end

@meshgen/private/renumberEdges.m

@@ -0,0 +1,22 @@
+function edges = renumberEdges(edges)
+%RENUMBEREDGES Renumber edges to maintain consistency after point modifications.
+%   edges = RENUMBEREDGES(edges) updates the indices of edge endpoints in
+%   the 'edges' matrix to reflect a new, sequential numbering based on
+%   their occurrence in 'edges'. This is useful after points have been
+%   filtered or reordered, ensuring that edge references are consistent
+%   with the updated point list.
+
+% Find unique indices of edge endpoints and create a direct mapping
+idx = unique(edges(:));
+for i = 1 : length(idx)
+    map(idx(i),1) = i;
+end
+
+% renumber bnde to correspond with pfix
+for i = 1 : size(edges,1)
+    edges(i,1) = map(edges(i,1));
+    edges(i,2) = map(edges(i,2));
+end
+
+
+end