adding high-fidelity option to meshgen for meshing hardened shoreline…

…s/coastal structures (#264)

- geodata.m: Fixed an issue where errors occurred if obj.contourfile contained multiple items.
- Read_shapefile.m: Addressed a bug where errors were thrown if finputname contained more than one item.
- Make_f15.m: General update.
- README.md: Updated documentation.
- meshgen: Introduced a high-fidelity option for improved meshing of hardened shorelines and coastal structures. This includes adding 1D meshing routines and tests, fixing issues related to merging with the master branch, and implementing high fidelity meshing routines. Also addressed bugs in subsetting with ocean boundary and handling cases where LN goes to zero, causing F to go to infinity (now zeros out).
- Example scripts: Added examples demonstrating the high fidelity meshing option, included a suffix for images to distinguish between different meshing strategies, and resolved issues from a bad merge with meshgen. Updated the setup in the runall example script and made format changes to align with recent updates, including the new Example13.
- meshgen.m: Restored mesh improvement strategy when not in high_fidelity mode.
- Test updates: Modified the JBAY test to reflect updated expected mesh outcomes and integrated changes into the documentation and example scripts to showcase the new features and fixes.

---------

Co-authored-by: omouraenko <44512173+omouraenko@users.noreply.github.com>
Co-authored-by: Keith Roberts <kroberts@baird.com>

@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