update in Jan 01 2023

.gitignore

@@ -3,7 +3,7 @@ datasets/
 fort.*
 Examples/*
 !Examples/Example_*
-/*.mat
+*/*.mat
 *.15
 *.14
 *.13

@meshgen/meshgen.m

@@ -39,7 +39,10 @@
     %         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
+    % delaunay_elim_on_exit % whether or not to run delaunay_elim on exit of meshgen   
+    % improve_with_reduced_quality % whether or not to allow mesh improvements with decreases in mesh quality
     %         big_mesh      % set to 1 to remove the bou data from memory
+    %       
     properties
         fd            % handle to distance function
         fh            % handle to edge function
@@ -73,6 +76,8 @@
         Fb            % bathymetry data interpolant
         enforceWeirs  % whether or not to enforce weirs in meshgen
         enforceMin    % whether or not to enfore minimum edgelength for all edgefxs
+        delaunay_elim_on_exit % whether or not to run delaunay_elim on exit of meshgen   
+        improve_with_reduced_quality % whether or not to allow mesh improvements with decreases in mesh quality
     end
 
 
@@ -126,6 +131,8 @@
             addOptional(p,'qual_tol',defval);
             addOptional(p,'enforceWeirs',1);
             addOptional(p,'enforceMin',1);
+            addOptional(p,'delaunay_elim_on_exit',1);
+            addOptional(p,'improve_with_reduced_quality',0);
 
             % parse the inputs
             parse(p,varargin{:});
@@ -137,7 +144,9 @@
             % 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',...
+            inp = orderfields(inp,{'h0','bbox','enforceWeirs','enforceMin',...
+                                   'delaunay_elim_on_exit','improve_with_reduced_quality',...
+                                   'fh',...
                                    'inner','outer','mainland',...
                                    'bou','ef',... %<--OceanMesh classes come after
                                    'egfix','pfix','fixboxes',...
@@ -409,6 +418,10 @@
                         obj.enforceWeirs = inp.(fields{i});
                     case('enforceMin')
                         obj.enforceMin = inp.(fields{i});
+                    case('delaunay_elim_on_exit')
+                        obj.delaunay_elim_on_exit = inp.(fields{i});
+                    case('improve_with_reduced_quality')
+                        obj.improve_with_reduced_quality = inp.(fields{i});
                 end
             end
 
@@ -660,9 +673,13 @@
 
 
                     % 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
+                    % ..or.. 
+                    % If not allowing improvements with reduction in quality 
+                    % then if the number of points significantly decreased
+                    % due to a mesh improvement iteration, then rewind.
+                    if ~mod(it,imp+1) && ((obj.qual(it,1) - obj.qual(it-1,1) < -0.10)  || ...
+                        (~obj.improve_with_reduced_quality && ...
+                        (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
@@ -734,7 +751,9 @@
                 if ~mod(it,imp)
                     if abs(qual_diff) < obj.qual_tol
                         % Do the final elimination of small connectivity
-                        [t,p] = delaunay_elim(p,obj.fd,geps,1);
+                        if obj.delaunay_elim_on_exit
+                            [t,p] = delaunay_elim(p,obj.fd,geps,1);
+                        end
                         disp('Quality of mesh is good enough, exit')
                         close all;
                         break;
@@ -796,7 +815,8 @@
                 if ~mod(it,imp)
                     p_before_improve = p;
                     nn = []; pst = [];
-                    if qual_diff < imp*obj.qual_tol && qual_diff > 0
+                    if abs(qual_diff) < imp*obj.qual_tol && ...
+                        (obj.improve_with_reduced_quality || qual_diff > 0)
                         % Remove elements with small connectivity
                         nn = get_small_connectivity(p,t);
                         disp(['Deleting ' num2str(length(nn)) ' due to small connectivity'])
@@ -870,7 +890,9 @@
 
                 if ( it > obj.itmax )
                     % Do the final deletion of small connectivity
-                    [t,p] = delaunay_elim(p,obj.fd,geps,1);
+                    if obj.delaunay_elim_on_exit
+                        [t,p] = delaunay_elim(p,obj.fd,geps,1);
+                    end
                     disp('too many iterations, exit')
                     close all;
                     break ;

@msh/msh.m

@@ -59,6 +59,7 @@
         f24 % A struct of the fort24 SAL values
         f2001 % A struct for the fort2001 non-periodic flux/ele sponge bc
         f5354 % A struct for the fort53001/54001 tidal ele/flux sponge bc
+        offset63 % A struct for the offset.63 dynamicwaterlevelcorrection file
         proj   % Description of projected space (m_mapv1.4)
         coord  % Description of projected space (m_mapv1.4)
         mapvar % Description of projected space (m_mapv1.4)
@@ -175,9 +176,9 @@
         end
 
         % write mesh to disk
-        function write(obj,fname,type)
+        function write(obj,fname,type,varargin)
             % Usage:
-            % write(obj,fname,type)
+            % write(obj,fname,type,varargin)
             %
             % Examples:
             % write(obj);       % writes all available data to fort_1.xx (ADCIRC) files
@@ -186,6 +187,7 @@ function write(obj,fname,type)
             % write(obj,fname,'gr3'); % writes mesh data to fname.gr3 (SCHISM) file
             % write(obj,fname,'ww3'); % writes mesh data to fname.ww3 (WaveWatchIII) file
             % write(obj,fname,{'13','14'}); % writes mesh data and f13 attribute data to fname.14 and fname.13 (ADCIRC) files
+            % write(obj,fname,'24','netcdf'); % writes fort.24 SAL data to fname.24.nc netcdf file
             if nargin == 1
                 fname = 'fort_1';
             end
@@ -211,11 +213,14 @@ function write(obj,fname,type)
                     writefort15( obj.f15, [fname '.15'], obj.bd );
                 end
                 if ~isempty(obj.f24)
-                    writefort24( obj.f24, [fname '.24'] );
+                    writefort24( obj.f24, [fname '.24'], obj.p, varargin);
                 end
                 if ~isempty(obj.f5354)
                     writefort5354( obj.f5354, fname );
                 end
+                if ~isempty(obj.offset63)
+                    writeoffset63( obj.offset63, [fname '.offset.63'] );
+                end
             else
                 if any(contains(type,'14')) || any(contains(type,'ww3')) || ...
                         any(contains(type,'gr3'))
@@ -264,11 +269,14 @@ function write(obj,fname,type)
                     writefort19( obj.f2001, [fname '.2001'] );
                 end
                 if any(contains(type,'24')) && ~isempty(obj.f24)
-                    writefort24( obj.f24, [fname '.24'] );
+                    writefort24( obj.f24, [fname '.24'], obj.p, varargin);
                 end
                 if any(contains(type,'5354')) && ~isempty(obj.f5354)
                     writefort5354( obj.f5354, fname );
                 end
+                if any(contains(type,'offset')) && ~isempty(obj.offset63)
+                    writeoffset63( obj.offset63, [fname '.offset.63'] );
+                end
             end
         end
 
@@ -1057,13 +1065,13 @@ function plotter(cmap,round_dec,yylabel,apply_pivot)
                 varargin(strcmp(varargin,'passive')) = [];
             elseif any(strcmp(varargin,'aggressive'))
                 disp('Employing aggressive option')
-                opt.db = 0.5; opt.ds = 1; opt.con = 9; opt.djc = 0.25;
+                opt.db = 0.5; opt.ds = 2; opt.con = 9; opt.djc = 0.25;
                 opt.sc_maxit = inf; opt.mqa = 0.5; opt.renum = 1;
                 varargin(strcmp(varargin,'aggressive')) = [];
             else
                 disp('Employing default (medium) option or user-specified opts')
-                opt.db = 0.25; opt.ds = 1; opt.con = 9; opt.djc = 0.1;
-                opt.sc_maxit = 1; opt.mqa = 0.25; opt.renum = 1;
+                opt.db = 0.25; opt.ds = 2; opt.con = 9; opt.djc = 0.1;
+                opt.sc_maxit = 0; opt.mqa = 0.25; opt.renum = 1;
                 varargin(strcmp(varargin,'default')) = [];
                 varargin(strcmp(varargin,'medium')) = [];
             end
@@ -1189,18 +1197,22 @@ function plotter(cmap,round_dec,yylabel,apply_pivot)
             qual = [mq_m,mq_l3sig,mq_l];
             LTn = size(obj.t,1);
 
-            if mq_l < opt.mqa && (opt.ds || LT ~= LTn)
+            if mq_l < opt.mqa && LT ~= LTn
                 % Need to clean it again
                 disp('Poor or overlapping elements, cleaning again')
                 disp(['(Min Qual = ' num2str(mq_l) ')'])
                 % repeat without projecting (already projected)
-                ii = find(strcmp(varargino,'proj'));
+                ii = find(strcmp(varargino,'proj'), 1);
                 if ~isempty(ii)
                     varargino{ii+1} = 0;
                 else
                     varargino{end+1} = 'proj';
                     varargino{end+1} = 0;
                 end
+                ii = find(strcmp(varargino,'pfix'), 1);
+                if ~isempty(ii)
+                    varargino{ii+1} = pfixV;
+                end
                 obj = clean(obj,varargino(:));
             elseif opt.nscreen
                 disp(['number of nodes is ' num2str(length(obj.p))])
@@ -3854,10 +3866,11 @@ function plotter(cmap,round_dec,yylabel,apply_pivot)
             egfix = renumberEdges(egfix) ;
         end
 
-        function boundary = get_boundary_of_mesh(obj,ascell)
-            % boundary = get_boundary_of_mesh(obj,ascell)
+        function [boundary, bou_index] = get_boundary_of_mesh(obj,ascell)
+            % [boundary, bou_index] = get_boundary_of_mesh(obj,ascell)
             %
-            % Returns the boundary of the mesh
+            % Returns the boundary of the mesh and/or the mesh indices of
+            % the boundary
             %
             % INPUTS: msh_obj
             % OUTPUTS: msh boundary in one of two forms:
@@ -3872,13 +3885,14 @@ function plotter(cmap,round_dec,yylabel,apply_pivot)
             end
             bnde = extdom_edges2(obj.t,obj.p) ;
             try
-                boundary = extdom_polygon(bnde,obj.p,-1) ;
+                [boundary,bou_index] = extdom_polygon(bnde,obj.p,-1) ;
             catch
                 warning('ALERT: Boundary of mesh is not walkable. Returning polylines.');
-                boundary = extdom_polygon(bnde,obj.p,-1,1) ;
+                [boundary,bou_index] = extdom_polygon(bnde,obj.p,-1,1) ;
             end
             if ascell; return; end
             boundary = cell2mat(boundary');
+            bou_index = cell2mat(bou_index');
         end
 
         function obj = map_mesh_properties(obj,varargin)
@@ -3946,8 +3960,13 @@ function plotter(cmap,round_dec,yylabel,apply_pivot)
                     % Remove open boundary info...
                     obj.op = [];
                 else
-                    % Remove unnessary part from the nbdv
+                    % Remove zero length boundary and unnessary part from the nbdv
                     obj.op.nbdv = obj.op.nbdv(1:max(obj.op.nvdll),:);
+                    zero_bound = obj.op.nvdll == 0;
+                    obj.op.nope = sum(~zero_bound);
+                    obj.op.ibtype(zero_bound) = [];
+                    obj.op.nvdll(zero_bound) = [];
+                    obj.op.nbdv(:,zero_bound) = [];
                 end
             end
             % land boundary info

@msh/private/GridData.m

@@ -135,7 +135,7 @@
     disp(['Bounding maximum depth to ',num2str(maxdepth), ' meters.']) ;  
 end
 
-if strcmp(type,'slope')
+if strcmp(type,'slope') 
     disp(['Computing cell-averaged slope using the ' slope_calc ' method'])
     if isempty(obj.b) || all(obj.b == 0) 
         error(['You must have the bathymetry on the mesh before ' ...
@@ -145,6 +145,11 @@
                  'affect the computation of the slopes'])
     end
 end
+if length(K) == 1
+    if strcmp(type,'slope') || strcmp(type,'all')
+        error('Cannot compute slopes on bathymetry with K of length 1')      
+    end
+end
 
 if strcmp(NaNs,'fill') || strcmp(NaNs,'fillinside')
     disp('Fill in NaNs using nearest neighbour interpolation.')
@@ -266,19 +271,22 @@
 
 % delta is max and min bounds of the connecting element centers (or if only
 % one connecting element then do difference with the vertex itself
-
-% kjr 10/17/2018 enlarge the cell-averaging stencil by a factor N 
-pcon_max = max(squeeze(max(pcon,[],1)),2*obj.p(K,:)-squeeze(min(pcon,[],1)));
-pcon_min = min(squeeze(min(pcon,[],1)),2*obj.p(K,:)-squeeze(max(pcon,[],1)));
+if length(K) == 1
+    pcon_max = max(max(squeeze(pcon)),2*obj.p(K,:)-max(squeeze(pcon)));
+    pcon_min = min(min(squeeze(pcon)),2*obj.p(K,:)-min(squeeze(pcon)));
+else
+    pcon_max = max(squeeze(max(pcon,[],1)),2*obj.p(K,:)-squeeze(min(pcon,[],1)));
+    pcon_min = min(squeeze(min(pcon,[],1)),2*obj.p(K,:)-squeeze(max(pcon,[],1)));
+end
 delta = pcon_max - pcon_min;
 
+% kjr 10/17/2018 enlarge the cell-averaging stencil by a factor N 
 pcon_max = N*pcon_max+(1-N)*obj.p(K,:) ;
 pcon_min = N*pcon_min+(1-N)*obj.p(K,:) ;
 
-
 %% Now read the bathy (only what we need)
 if ~isa(geodata,'geodata')
-    BufferL = max(delta); 
+    BufferL = max(delta,[],1); 
     lon_min = min(obj.p(K,1)) - BufferL(1);
     lon_max = max(obj.p(K,1)) + BufferL(1);
     lat_min = min(obj.p(K,2)) - BufferL(2);
@@ -339,7 +347,6 @@
 end
 
 %% Calculate N - number of DEM grids to average - for each node
-tic
 if strcmp(interp,'CA')
     [~,IDXX,IDXY] = FindLinearIdx(obj.p(K,1),obj.p(K,2),DEM_X,DEM_Y);
     [~,IDXR,IDXT] = FindLinearIdx(pcon_max(:,1),pcon_max(:,2),DEM_X,DEM_Y);
@@ -382,7 +389,8 @@
         if nanfill
             if sum(isnan(b)) > 0
                 localcoord = obj.p(K,:);
-                KI = knnsearch(localcoord(~isnan(b),:),localcoord(isnan(b),:));
+                [KI, ~] = ourKNNsearch(localcoord(~isnan(b),:),',localcoord(isnan(b)',1); 
+                %KI = knnsearch(localcoord(~isnan(b),:),localcoord(isnan(b),:));
                 bb = b(~isnan(b));
                 b(isnan(b)) = bb(KI); clear bb localcoord
             end
@@ -411,13 +419,15 @@
             % Try and fill in the NaNs
             if ~isempty(find(isnan(bx),1))
                 localcoord = obj.p(K,:);
-                KI = knnsearch(localcoord(~isnan(bx),:),localcoord(isnan(bx),:));
+                [KI, ~] = ourKNNsearch(localcoord(~isnan(bx),:)',localcoord(isnan(bx),:)',1); 
+                %KI = knnsearch(localcoord(~isnan(bx),:),localcoord(isnan(bx),:));
                 bb = bx(~isnan(bx));
                 bx(isnan(bx)) = bb(KI); clear bb localcoord
             end
             if ~isempty(find(isnan(by),1))
                 localcoord = obj.p(K,:);
-                KI = knnsearch(localcoord(~isnan(by),:),localcoord(isnan(by),:));
+                [KI, ~] = ourKNNsearch(localcoord(~isnan(by),:)',localcoord(isnan(by),:)',1); 
+                %KI = knnsearch(localcoord(~isnan(by),:),localcoord(isnan(by),:));
                 bb = by(~isnan(by));
                 by(isnan(by)) = bb(KI); clear bb localcoord
             end
@@ -437,7 +447,6 @@
         b = F(obj.p(K,1),obj.p(K,2));
     end
 end
-toc
 %% Put in the global msh class
 if strcmp(type,'depth') || strcmp(type,'all')
     if isempty(obj.b)

@msh/private/writefort15.m

@@ -175,7 +175,7 @@
 
     nm = 0 ; 
     for ib = 1: boudat.nbou
-        ibty = boudat.ibtype(ib)
+        ibty = boudat.ibtype(ib);
 
         switch ibty
             case {2,12,22,32,52}
@@ -197,7 +197,7 @@
         % val = fscanf(fid, '%f ' ) ; % Must be revisit
         icnt = 0 ;
         for ib = 1: boudat.nbou
-            ibty = boudat.ibtype(ib) 
+            ibty = boudat.ibtype(ib); 
 
             switch ibty
                 case {2,12,22,52}
@@ -316,8 +316,9 @@
     fprintf( fid, '! -- Begin %s Control Namelist -- \n', f15dat.controllist(k).type ) ;
     fprintf( fid, '&%sControl\n', f15dat.controllist(k).type ) ;
     for m = 1:length(f15dat.controllist(k).var)
-        fprintf( fid, '%s = %s,\n',f15dat.controllist(k).var(m).name,...
-                 f15dat.controllist(k).var(m).val) ;
+        val = f15dat.controllist(k).var(m).val;
+        if ~ischar(val); val = num2str(val); end
+        fprintf( fid, '%s = %s,\n',f15dat.controllist(k).var(m).name,val) ;
     end
     fprintf( fid, '/\n') ;
     fprintf( fid, '! -- End %s Control Namelist -- \n', f15dat.controllist(k).type ) ;