diff --git a/.gitignore b/.gitignore
index 4121b412..51c44b12 100644
--- a/.gitignore
+++ b/.gitignore
@@ -3,7 +3,7 @@ datasets/
fort.*
Examples/*
!Examples/Example_*
-/*.mat
+*/*.mat
*.15
*.14
*.13
diff --git a/@meshgen/meshgen.m b/@meshgen/meshgen.m
index 511e24c1..1b7f1120 100644
--- a/@meshgen/meshgen.m
+++ b/@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 ;
diff --git a/@msh/msh.m b/@msh/msh.m
index 52f52586..3e9d71ba 100644
--- a/@msh/msh.m
+++ b/@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
diff --git a/@msh/private/GridData.m b/@msh/private/GridData.m
index b4205326..0024e577 100644
--- a/@msh/private/GridData.m
+++ b/@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)
diff --git a/@msh/private/writefort15.m b/@msh/private/writefort15.m
index a72210d4..4b6bf0d0 100644
--- a/@msh/private/writefort15.m
+++ 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 ) ;
diff --git a/@msh/private/writefort24.m b/@msh/private/writefort24.m
index 0368bd32..4e81035c 100644
--- a/@msh/private/writefort24.m
+++ b/@msh/private/writefort24.m
@@ -1,26 +1,92 @@
-function fid = writefort24( f24dat, finame )
+function fid = writefort24( f24dat, finame, points, format )
%
%
-if ( nargin == 1 )
- finputname = 'fort.24_1' ;
+if ( nargin == 1 )
+ finputname = 'fort.24_1' ;
else
- finputname = finame ;
+ finputname = finame ;
end
-fid = fopen(finputname,'w') ;
+if nargin < 3 || isempty(format)
+ format = 'ascii';
+ warning('Using ASCII file format. Pass ''netcdf'' to write in NetCDF')
+end
+if iscell(format)
+ format = format{1};
+end
+
+if strcmp(format,'ascii') % LEGACY
+ fid = fopen(finputname,'w') ;
+
+ tipnames = f24dat.tiponame; ntip = length(tipnames) ;
+ for icon = 1: ntip
+ tipname = tipnames{icon};
+ fprintf('Writing SAL %s data \n', char(tipname)) ;
+ % The constituent details
+ fprintf(fid,'%s \n',[char(tipname) ' SAL']) ;
+ fprintf(fid,'%17.15f \n',f24dat.omega(icon)) ;
+ fprintf(fid,'%d \n',1) ;
+ fprintf(fid,'%s \n',char(tipname)) ;
+ fprintf(fid,'%d \t %12.6f %12.6f \n',f24dat.Val(icon,:,:));
+ end
+
+ fclose(fid) ;
+
+elseif strcmp(format,'netcdf') % NETCDF
+ % create nc file or overwrite it
+ file = [finputname,'.nc'];
+ if exist(file)
+ delete(file)
+ end
+ % deflate level (set zero for no deflation if necessary.. using single precision so not very different in size)
+ dl = 5;
+
+ node = length(points);
+ tipnames = char(f24dat.tiponame); ntip = size(tipnames) ;
+ fillvalue = 0.0;
+
+ nccreate(file,'x','Dimensions',{'node',node},'DataType','int32','DeflateLevel',dl) ;
+ ncwriteatt(file,'x','standard_name', 'longitude') ;
+ ncwriteatt(file,'x','units', 'degrees_east') ;
+ ncwriteatt(file,'x','positive', 'east') ;
+ ncwrite(file, 'x', points(:,1)) ;
-tipnames = f24dat.tiponame; ntip = length(tipnames) ;
-for icon = 1: ntip
- tipname = tipnames{icon};
- fprintf('Writing SAL %s data \n', char(tipname)) ;
- % The constituent details
- fprintf(fid,'%s \n',[char(tipname) ' SAL']) ;
- fprintf(fid,'%17.15f \n',f24dat.omega(icon)) ;
- fprintf(fid,'%d \n',1) ;
- fprintf(fid,'%s \n',char(tipname)) ;
- fprintf(fid,'%d \t %12.6f %12.6f \n',f24dat.Val(icon,:,:));
+ nccreate(file,'y','Dimensions',{'node',node},'DataType','int32','DeflateLevel',dl) ;
+ ncwriteatt(file,'y','standard_name', 'latitude') ;
+ ncwriteatt(file,'y','units', 'degrees_north') ;
+ ncwriteatt(file,'y','positive', 'north') ;
+ ncwrite(file, 'y', points(:,2)) ;
+
+ nccreate(file,'constituents','Dimensions',{'num_constituents',ntip(1),'char_len',ntip(2)},'DataType','char') ;
+ ncwriteatt(file,'constituents','standard_name', 'name_of_tidal_harmonic_constituents') ;
+ ncwriteatt(file,'constituents','long_name', 'name of tidal harmonic constituents') ;
+ ncwrite(file, 'constituents', tipnames) ;
+
+ nccreate(file,'frequency','dimensions',{'num_constituents',ntip(1)});
+ ncwriteatt(file,'frequency','standard_name','frequency_of_tidal_harmonic_constituents')
+ ncwriteatt(file,'frequency','long_name','frequency of tidal harmonic constituents')
+ ncwriteatt(file,'frequency','units','radians/second')
+ ncwrite(file, 'frequency', f24dat.omega) ;
+
+ nccreate(file, 'sal_amplitude','Dimensions',{'num_constituents',ntip(1),'node',node},'DataType','single','DeflateLevel',dl,'FillValue',fillvalue) ;
+ ncwriteatt(file,'sal_amplitude','standard_name','amplitude_of_self_attraction_and_loading_tide_elevation')
+ ncwriteatt(file,'sal_amplitude','long_name','amplitude of self attraction and loading tide elevation')
+ ncwriteatt(file,'sal_amplitude','units','m')
+ ncwrite(file, 'sal_amplitude', squeeze(f24dat.Val(:,2,:))) ;
+
+ nccreate(file, 'sal_phase','Dimensions',{'num_constituents',ntip(1),'node',node},'DataType','single','DeflateLevel',dl,'FillValue',fillvalue) ;
+ ncwriteatt(file,'sal_phase','long_name','phase-lag of self-attraction and loading tide elevation')
+ ncwriteatt(file,'sal_phase','standard_name','phase_lag_of_self_attraction_and_loading_tide_elevation')
+ ncwriteatt(file,'sal_phase','units','degrees with respect to GMT/UTC')
+ ncwrite(file, 'sal_phase', squeeze(f24dat.Val(:,3,:))) ;
+
+ ncwriteatt(file,'/','title','The self-attraction and loading terms for an ADCIRC simulation');
+ ncwriteatt(file,'/','creation_date',datestr(now));
+ ncwriteatt(file,'/','source',"Made by OceanMesh2D writefort24");
+ ncwriteatt(file,'/','references',"https://github.com/CHLNDDEV/OceanMesh2D/" );
+
+else
+ error(['format = ' format ' is invalid. Choose from ascii or netcdf'])
end
-
-fclose(fid) ;
%EOF
-end
\ No newline at end of file
+end
diff --git a/@msh/private/writeoffset63.m b/@msh/private/writeoffset63.m
new file mode 100644
index 00000000..40b4e1ac
--- /dev/null
+++ b/@msh/private/writeoffset63.m
@@ -0,0 +1,22 @@
+function fid = writeoffset63( f63dat, finame )
+%
+%
+if ( nargin == 1 )
+ finputname = 'offset.63' ;
+else
+ finputname = finame ;
+end
+
+fid = fopen(finputname,'w') ;
+fprintf(fid,'%s \n',f63dat.header) ;
+fprintf(fid,'%12.6f \n',f63dat.time_interval) ;
+fprintf(fid,'%12.6f \n',f63dat.default_val) ;
+
+for tt = 1: f63dat.num_times
+ fprintf(fid,'%d \t %12.6f \n',[f63dat.offset_nodes; f63dat.offset_values(tt,:)]);
+ fprintf(fid,'%s \n','##');
+end
+
+fclose(fid) ;
+%EOF
+end
diff --git a/Examples/Example_11_Remeshing_Patches.m b/Examples/Example_11_Remeshing_Patches.m
index 49b151aa..c3c4b302 100644
--- a/Examples/Example_11_Remeshing_Patches.m
+++ b/Examples/Example_11_Remeshing_Patches.m
@@ -38,7 +38,7 @@
% Get out the msh class and put on nodestrings
m = mshopts.grd;
% Interpolate topobathy data onto the vertices of the mesh.
-m = interp(m,'SRTM15+V2.1.nc');
+m = interp(m,'SRTM15+.nc');
%% Extract the region which you want to remesh
% For the purpose of this example, we have drawn a random polygon on top
% of the mesh for which we would like to remesh.
diff --git a/Examples/Example_12_Riverine_flow.m b/Examples/Example_12_Riverine_flow.m
index 38c1ba39..19218898 100644
--- a/Examples/Example_12_Riverine_flow.m
+++ b/Examples/Example_12_Riverine_flow.m
@@ -43,7 +43,7 @@
%% STEP 2: specify geographical datasets and process the geographical data
%% to be used later with other OceanMesh classes...
-dem = 'SRTM15+V2.1.nc';
+dem = 'SRTM15+.nc';
coastline = 'GSHHS_f_L1';
gdat = geodata('shp',coastline,'dem',dem,'h0',min_el,'bbox',bbox);
%plot(gdat,'shp');
diff --git a/Examples/Example_3_ECGC.m b/Examples/Example_3_ECGC.m
index c410ce27..e5ca2695 100644
--- a/Examples/Example_3_ECGC.m
+++ b/Examples/Example_3_ECGC.m
@@ -24,7 +24,7 @@
%% STEP 2: specify geographical datasets and process the geographical data
%% to be used later with other OceanMesh classes...
-dem = 'SRTM15+V2.1.nc';
+dem = 'SRTM15+.nc';
coastline = 'GSHHS_f_L1';
gdat1 = geodata('shp',coastline,'dem',dem,'h0',min_el,...
'bbox',bbox);
diff --git a/Examples/Example_4_PRVI.m b/Examples/Example_4_PRVI.m
index 2fda5a46..99130345 100644
--- a/Examples/Example_4_PRVI.m
+++ b/Examples/Example_4_PRVI.m
@@ -24,7 +24,7 @@
max_el = 10e3; % maximum resolution in meters.
coastline = 'GSHHS_f_L1';
-dem = 'SRTM15+V2.1.nc';
+dem = 'SRTM15+.nc';
gdat{1} = geodata('shp',coastline,...
'dem',dem,...
'bbox',bbox,...
diff --git a/Examples/Example_7_Global.m b/Examples/Example_7_Global.m
index b64b181f..bdbdb84e 100644
--- a/Examples/Example_7_Global.m
+++ b/Examples/Example_7_Global.m
@@ -22,7 +22,7 @@
%% STEP 2: specify geographical datasets and process the geographical data
%% to be used later with other OceanMesh classes...
-dem = 'SRTM15+V2.1.nc';
+dem = 'SRTM15+.nc';
coastline = {'GSHHS_f_L1','GSHHS_f_L6'};
gdat1 = geodata('shp',coastline,'dem',dem,...
'bbox',bbox,'h0',min_el);
diff --git a/Examples/Example_8_AK.m b/Examples/Example_8_AK.m
index 83e677fe..7a54f120 100644
--- a/Examples/Example_8_AK.m
+++ b/Examples/Example_8_AK.m
@@ -25,7 +25,7 @@
%% STEP 2: specify geographical datasets and process the geographical data
%% to be used later with other OceanMesh classes...
-dem = 'SRTM15+V2.1.nc';
+dem = 'SRTM15+.nc';
coastline = 'GSHHS_f_L1';
gdat = geodata('shp',coastline,'bbox',bbox,'h0',min_el);
% plotting the gdat to show it crosses the 180/-180 and
diff --git a/Examples/Example_9_TBAY.m b/Examples/Example_9_TBAY.m
index 61d77e03..ceeab507 100644
--- a/Examples/Example_9_TBAY.m
+++ b/Examples/Example_9_TBAY.m
@@ -31,7 +31,7 @@
% Shoreline
shoreline = 'GSHHS_f_L1';
% Bathy data
-dem = 'SRTM15+V2.1.nc';
+dem = 'SRTM15+.nc';
%% Define meshing domain and parameters
bbox{1} = [-83, -82; 27 28.5]; % outermost bbox
min_el = 100; % minimum edgelength
diff --git a/README.md b/README.md
index 7d073470..397e3a42 100644
--- a/README.md
+++ b/README.md
@@ -27,7 +27,7 @@ OceanMesh2D is a set of user-friendly MATLAB functions to generate two-dimension
Getting help
==============
-Besides posting [issues](https://github.com/CHLNDDEV/OceanMesh2D/issues) with the code on Github, you can also ask questions via our Slack channel [here](https://join.slack.com/t/oceanmesh2d/shared_invite/zt-r5ddyn78-xapbXmJmeuKwLPB5ZiSNtQ).
+Besides posting [issues](https://github.com/CHLNDDEV/OceanMesh2D/issues) with the code on Github, you can also ask questions via our Slack channel [here](https://join.slack.com/t/oceanmesh2d/shared_invite/zt-1b96mhvhw-sHUhP2emepHlGtw0~fWmAg).
Note: If the slack link invite isn't working, please send either one of us and an email and we'll fix it. By default, the invite link expires every 30 days.
@@ -142,6 +142,14 @@ GALLERY:
+
+* The following images are provided from happy users. Please send us your meshes.
+
+Jiangchao Qiu from his [publication](https://agupubs.onlinelibrary.wiley.com/doi/abs/10.1029/2021EF002638).
+
+ 
+
+
Changelog
=========
@@ -149,7 +157,21 @@ The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.0.0/)
### Unreleased (on current HEAD of the Projection branch)
## Added
+- Option `improve_with_reduced_quality` to `meshgen` for allowing for mesh improvements even when quality is decreasing or large number of nodes eliminated, which sometimes is necessary to force the advancement in mesh quality.
+- Option `delaunay_elim_on_exit` to `meshgen` to skip the last call to `delaunay_elim` to potentially avoid deleting boundary elements.
- Geoid offset nodal attribute in `Calc_f13` subroutine. https://github.com/CHLNDDEV/OceanMesh2D/pull/251
+- Support for writing Self Attraction and Loading (SAL) files in NetCDF for the ADCIRC model. https://github.com/CHLNDDEV/OceanMesh2D/pull/231
+## Changed
+- Default mesh improvement strategy is `ds` 2.
+- Retrieve boundary indices in `msh.get_boundary_of_mesh` method. https://github.com/CHLNDDEV/OceanMesh2D/pull/259
+- `msh.offset63` struct and associated write/make routines for dynamicwaterlevel offset functionality. https://github.com/CHLNDDEV/OceanMesh2D/pull/259
+- dynamicWaterLevelCorrection to fort.15 namelist, and PRBCKGRND option to met fort.15 namelist. https://github.com/CHLNDDEV/OceanMesh2D/pull/261
+## Fixed
+- Recursive cleaning issues: infinite loop and preservation of fixed points.
+- `msh.interp` method for `K` argument of length 1, and for the test to determine whether the bathymetry grid is irregular. https://github.com/CHLNDDEV/OceanMesh2D/pull/259
+- Printing of namelist character strings or numbers. https://github.com/CHLNDDEV/OceanMesh2D/pull/261
+- `Make_offset63.m` time interval computation. https://github.com/CHLNDDEV/OceanMesh2D/pull/261
+- Removed dependency on statistics toolbox when using the 'nanfill' option in `msh.interp`. https://github.com/CHLNDDEV/OceanMesh2D/pull/269
### [5.0.0] - 2021-07-29
## Added
diff --git a/Tests/RunTests.m b/Tests/RunTests.m
index e98e2a07..3ee2c400 100644
--- a/Tests/RunTests.m
+++ b/Tests/RunTests.m
@@ -18,13 +18,13 @@
TestInterp
- if exist('SRTM15+V2.1.nc')
+ if exist('SRTM15+.nc')
TestECGC
else
- warning('Need to download SRTM15+V2.1.nc to run TestECGC. Run setup.sh')
+ warning('Need to download SRTM15+.nc to run TestECGC. Run setup.sh')
end
diff --git a/Tests/TestECGC.m b/Tests/TestECGC.m
index 8752b1f2..bf0a2fc2 100644
--- a/Tests/TestECGC.m
+++ b/Tests/TestECGC.m
@@ -20,7 +20,7 @@
R = 3; % Number of elements to resolve feature.
-dem = 'SRTM15+V2.1.nc';
+dem = 'SRTM15+.nc';
coastline = 'GSHHS_f_L1';
gdat1 = geodata('shp',coastline,'dem',dem,'h0',min_el,...
'bbox',bbox);
diff --git a/Tests/TestSanity.m b/Tests/TestSanity.m
index aac669fe..101d8d99 100644
--- a/Tests/TestSanity.m
+++ b/Tests/TestSanity.m
@@ -4,7 +4,7 @@
NP_TOL = 500;
NT_TOL = 1500;
-QUAL_TOL = 0.5;
+QUAL_TOL = 0.25; %[mqa for medium cleaning option]
% p: [5968�2 double]
% t: [9530�3 double]
% Element qual. metrics
diff --git a/imgs/PRD_Mesh.png b/imgs/PRD_Mesh.png
new file mode 100644
index 00000000..2de3f146
Binary files /dev/null and b/imgs/PRD_Mesh.png differ
diff --git a/setup.sh b/setup.sh
index 3e0b5e6a..490ad878 100755
--- a/setup.sh
+++ b/setup.sh
@@ -38,11 +38,11 @@ if $gshhs; then
fi
if $srtm; then
- if [ -f "SRTM15+V2.1.nc" ]; then
- echo "SRTM15+V2.1.nc global bathymetry file already exists"
+ if [ -f "SRTM15+.nc" ]; then
+ echo "SRTM15+.nc global bathymetry file already exists"
else
- # download SRTM15+V2.1 bathymetry
- wget "ftp://topex.ucsd.edu/pub/srtm15_plus/SRTM15+V2.1.nc"
+ # download SRTM15+ bathymetry
+ wget "https://topex.ucsd.edu/pub/srtm15_plus/SRTM15_V2.4.nc"" -O SRTM15+.nc
fi
fi
@@ -51,6 +51,6 @@ if $gebco; then
echo "GEBCO_2020.nc global bathymetry file already exists"
else
# download GEBCO_2020.nc bathymetry
- wget "https://www.bodc.ac.uk/data/open_download/gebco/gebco_2020/zip/" -O GEBCO_2020.nc
+ wget "https://www.bodc.ac.uk/data/open_download/gebco/gebco_2021/zip/" -O GEBCO_2020.nc
fi
fi
diff --git a/utilities/FindLinearIdx.m b/utilities/FindLinearIdx.m
index e242c749..092dba67 100644
--- a/utilities/FindLinearIdx.m
+++ b/utilities/FindLinearIdx.m
@@ -9,6 +9,11 @@
nx = size(lon,2);
np = numel(x);
+if nx == 1 && ny == 1
+ IX = 1; IX1 = 1; IX2 = 1;
+ return
+end
+
% make sure entry points are column vectors
x = x(:);
y = y(:);
@@ -17,8 +22,8 @@
% (trying both directions so could be meshgrid or ndgrid format)
dx = diff(lon(:,1));
dy = diff(lat(1,:));
-
-if max(dx) ~= min(dx) || max(dy) ~= min(dy)
+dx_cutoff = 0.1/111e3; % approx 10 cm
+if (max(dx) - min(dx)) > dx_cutoff || (max(dy) - min(dy)) > dx_cutoff
% % IRREGULAR GRID (SLOWER)
% convert ndgrid to vector
diff --git a/utilities/Make_f15.m b/utilities/Make_f15.m
index 15568552..98fa50d8 100644
--- a/utilities/Make_f15.m
+++ b/utilities/Make_f15.m
@@ -186,16 +186,33 @@
f15dat.extraline(10).msg = '';
% control lists
+ % met
f15dat.controllist(1).type = 'met';
f15dat.controllist(1).var(1).name = 'WindDragLimit';
f15dat.controllist(1).var(1).val = 0.0025;
- f15dat.controllist(1).var(2).name = 'DragLawString';
- f15dat.controllist(1).var(2).val = 'default';
- f15dat.controllist(1).var(3).name = 'outputWindDrag';
- f15dat.controllist(1).var(3).val = 'F';
- f15dat.controllist(1).var(4).name = 'invertedBarometerOnElevationBoundary';
+ f15dat.controllist(1).var(2).name = 'PRBCKGRND';
+ f15dat.controllist(1).var(2).val = 1013;
+ f15dat.controllist(1).var(3).name = 'DragLawString';
+ f15dat.controllist(1).var(3).val = 'default';
+ f15dat.controllist(1).var(4).name = 'outputWindDrag';
f15dat.controllist(1).var(4).val = 'F';
-
+ f15dat.controllist(1).var(5).name = 'invertedBarometerOnElevationBoundary';
+ f15dat.controllist(1).var(5).val = 'F';
+ % dynamicwaterlevelcorrection control
+ f15dat.controllist(2).type = 'dynamicWaterLevelCorrection';
+ f15dat.controllist(2).var(1).name = [f15dat.controllist(2).type 'FileName'];
+ f15dat.controllist(2).var(1).val = 'offset.63';
+ f15dat.controllist(2).var(2).name = [f15dat.controllist(2).type 'Multiplier'];
+ f15dat.controllist(2).var(2).val = 1.0;
+ f15dat.controllist(2).var(3).name = [f15dat.controllist(2).type 'RampStart'];
+ f15dat.controllist(2).var(3).val = 0.0;
+ f15dat.controllist(2).var(4).name = [f15dat.controllist(2).type 'RampEnd'];
+ f15dat.controllist(2).var(4).val = 0.0;
+ f15dat.controllist(2).var(5).name = [f15dat.controllist(2).type 'RampReferenceTime'];
+ f15dat.controllist(2).var(5).val = 'coldstart';
+ f15dat.controllist(2).var(6).name = [f15dat.controllist(2).type 'SkipSnaps'];
+ f15dat.controllist(2).var(6).val = 0;
+
% Put into the msh class
obj.f15 = f15dat;
end
diff --git a/utilities/Make_f24.m b/utilities/Make_f24.m
index 2a93e294..997f8279 100644
--- a/utilities/Make_f24.m
+++ b/utilities/Make_f24.m
@@ -3,82 +3,91 @@
% Takes a msh object and interpolates the global SAL term into the f24
% struct
% Assumes that saldata is in the MATLAB path
-% The saldata required can be downloaded from:
-% saldata = 'FES2004' : Source at: ftp://ftp.legos.obs-mip.fr/pub/soa/...
-% maree/tide_model/global_solution/fes2004/
-%
+% The saldata required can be downloaded from:
+% saldata = 'FES2004' : Source at: ftp://ftp.legos.obs-mip.fr/pub/soa/...
+% maree/tide_model/global_solution/fes2004/load/
+%
% saldata = 'FES2014' : Source at: ftp://ftp.legos.obs-mip.fr/pub/...
-% FES2012-project/data/LSA/FES2014/
-% by default saldata = 'FES2014'
+% FES2012-project/data/LSA/FES2014/
+% by default saldata = 'FES2014'
%
% plot_on - 1/true: to plot and print F24 values for checking
% 0/false: no plotting by default
%
-% Created by William Pringle. July 11 2018 updated to Make_f## style
+% Created by William Pringle. July 11 2018 updated to Make_f## style
+% Modified by Keith Roberts Jun 19, 2021 to specify degree format for FES
+% file
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if isempty(obj.f15)
- error(['The msh object must have the f15 struct populated '...
- 'with tidal potential information'])
+ error(['The msh object must have the f15 struct populated '...
+ 'with tidal potential information'])
end
if nargin < 2 || isempty(saldata)
- saldata = 'FES2014';
+ saldata = 'FES2014';
end
if nargin < 3 || isempty(plot_on)
- plot_on = false;
-end
-
-ll0 = obj.f15.slam(1) ;
-if ( ll0 < 0 )
- ll0 = ll0 + 360 ;
+ plot_on = false;
end
-lon0 = ll0*pi/180 ; lat0 = obj.f15.slam(2)*pi/180 ;
-
-R = 6378206.4; % earth radius
+% parameter for cpp conversion
+R = 6378206.4; % earth radius
+lon0 = obj.f15.slam(1) ; % central longitude
+lat0 = obj.f15.slam(2) ; % central latitude
% Put in the tidal potential names
obj.f24.tiponame = {obj.f15.tipotag.name};
-ntip = length(obj.f24.tiponame) ;
+ntip = length(obj.f24.tiponame) ;
% choose tidal database file names and directories
database = strtrim(upper(saldata)) ;
direc = '';
-% % Load tide grid data
+% % Load tide grid data
if strcmp(database,'FES2004')
tide_grid = [direc 'load.k1.nc'];
tide_prefix = [direc 'load.'];
tide_suffix = '.nc';
lon = ncread(tide_grid,'lon');
lat = ncread(tide_grid,'lat');
-elseif strcmp(database,'FES2014')
+ % -180/180 degree format for 2004
+ if (lon0 > 180); lon0 = lon0 - 360 ; end
+elseif strcmp(database,'FES2014')
tide_grid = [direc 'K1_sal.nc'];
tide_prefix = direc;
tide_suffix = '_sal.nc';
lon = ncread(tide_grid,'longitude');
lat = ncread(tide_grid,'latitude');
[lon,lat] = ndgrid(lon,flipud(lat));
+ % 0-360 degree format for 2014
+ if (lon0 < 0); lon0 = lon0 + 360 ; end
+else
+ error(['database = ' database ' is invalid.'...
+ ' Choose from FES2004 or FES2014'])
end
+% Convert CPP origin parameters to radians
+lon0 = lon0*pi/180 ; lat0 = lat0*pi/180 ;
-% CPP Conversion of lat/lon
-lon = lon * pi/180; lat = lat * pi/180;
+% CPP Conversion of FES lat/lon
+lon = lon * pi/180; lat = lat * pi/180;
x = R * (lon - lon0) * cos(lat0);
y = R * lat;
-% Get mesh vertices and change to FES 0 to 360 deg style
+% Doing the CPP conversion of the mesh
VX = obj.p;
-VX(VX(:,1)<0,1) = VX(VX(:,1)<0,1) + 360;
-
-% Doing the CPP conversion
-xx = VX(:,1) * pi/180; yy = VX(:,2) * pi/180;
+if strcmp(database,'FES2004')
+ VX(VX(:,1)>180,1) = VX(VX(:,1)>180,1) - 360;
+elseif strcmp(database,'FES2014')
+ VX(VX(:,1)<0,1) = VX(VX(:,1)<0,1) + 360;
+end
+xx = VX(:,1) * pi/180; yy = VX(:,2) * pi/180;
xx = R * (xx - lon0) * cos(lat0);
yy = R * yy;
% Now interpolate onto grid and put into fort.24 struct
nnodes = length(VX) ;
-kvec = (1:nnodes)';
+kvec = (1:nnodes)';
obj.f24.Val = zeros(ntip,3,nnodes) ;
for icon = 1: ntip
% Get the frequency
@@ -106,8 +115,8 @@
% Do the gridded Interpolation
F = griddedInterpolant(x,y,z,'linear','none');
- Z = F(xx,yy);
-
+ Z = F(xx,yy);
+
% Convert back to amp and phase
amp = abs(Z);
phs = angle(Z)*180/pi;
@@ -119,20 +128,24 @@
% Plot interpolated results
if plot_on
- figure(1); fastscatter(VX(:,1),VX(:,2),amp);
- colorbar;
- constituent = obj.f24.tiponame{icon};
- title(constituent)
- print(['F24_' constituent '_check'],'-dpng')
- end
-
+ figure(1); fastscatter(VX(:,1),VX(:,2),amp);
+ colorbar;
+ constituent = obj.f24.tiponame{icon};
+ title(constituent)
+ print(['F24_' constituent '_check'],'-dpng')
+ end
+
% Put into the struct
- obj.f24.Val(icon,:,:) = [kvec'; amp'; phs'];
+ obj.f24.Val(icon,:,:) = [kvec'; amp'; phs'];
+
+ if any(isnan(amp))
+ warning('NaNs detected in amplitudes. Is degree format correct?')
+ end
end
if obj.f15.ntip ~= 2
- disp('Setting ntip = 2')
- obj.f15.ntip = 2;
+ disp('Setting ntip = 2')
+ obj.f15.ntip = 2;
end
%EOF
end
diff --git a/utilities/Make_offset63.m b/utilities/Make_offset63.m
new file mode 100644
index 00000000..bf65c337
--- /dev/null
+++ b/utilities/Make_offset63.m
@@ -0,0 +1,41 @@
+function obj = Make_offset63(obj,time_vector,offset_nodes,offset_values)
+% obj = Make_offset63(obj,time_vector,offset_nodes,offset_values)
+%
+% Inputs: obj - msh class
+% time_vector - datetime vector : size NTx1 or 1xNT
+% offset_nodes - offset nodes : size 1xNP
+% offset_values - offset values : size NTxNP
+%
+% Author: William Pringle
+% Created: Nov 22 2021
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+if size(offset_nodes,1) ~= 1
+ if size(offset_nodes,2) == 1
+ offset_nodes = offset_nodes';
+ else
+ error('offset_nodes size must be 1xNP')
+ end
+end
+if size(offset_values,1) ~= length(time_vector)
+ error('offset_values size of 1st dimension must be same as time_vector length')
+end
+if size(offset_values,2) ~= length(offset_nodes)
+ error('offset_values size of 2nd dimension must be same as offset_nodes length')
+end
+
+% enter in msh.offset63 struct
+obj.offset63.header = ['#dynamicwaterleveloffsets from ' ...
+ datestr(time_vector(1)) ' to ' datestr(time_vector(end))];
+obj.offset63.num_times = length(time_vector);
+obj.offset63.time_interval = round(...
+ seconds(time_vector(end) - time_vector(1))/ ...
+ (length(time_vector)-1)...
+ );
+obj.offset63.default_val = 0;
+obj.offset63.offset_nodes = offset_nodes;
+obj.offset63.offset_values = offset_values;
+%EOF
+end
+
+
diff --git a/utilities/refine2_om.m b/utilities/refine2_om.m
index dea30cc7..3bb05f4c 100644
--- a/utilities/refine2_om.m
+++ b/utilities/refine2_om.m
@@ -1116,7 +1116,7 @@
nold = size(vert,1);
vert = [vert; new1(:,1:2)];
vert = [vert; new2(:,1:2)];
- % make sure vertices are unique
+ % wjp: make sure vertices are unique
vert = unique(vert,'rows','stable');
nnew = size(vert,1);