update for making a riverine boundary

@msh/msh.m

@@ -257,6 +257,9 @@ function write(obj,fname,type)
                 if any(contains(type,'19')) && ~isempty(obj.f19)
                     writefort19( obj.f19, [fname '.19'] );
                 end
+                if any(contains(type,'20')) && ~isempty(obj.f20)
+                    writefort20( obj.f20, [fname '.20'] );
+                end
                 if any(contains(type,'2001')) && ~isempty(obj.f2001)
                     writefort19( obj.f2001, [fname '.2001'] );
                 end
@@ -632,14 +635,15 @@ function write(obj,fname,type)
                     disp(['Trying to plot arbitrary f13 attribute: ' type])
                     if ~isempty(obj.f13)
                         ii = find(contains({obj.f13.defval.Atr(:).AttrName},type));
+                        jj = find(contains({obj.f13.userval.Atr(:).AttrName},type));
                         if isempty(ii)
                             error(['no f13 attribute of ' type ' found']);
                         elseif length(ii) > 1
                             disp({obj.f13.defval.Atr(ii).AttrName})
                             error(['more than one f13 attribute matched ' type]);
                         end
                         defval  = obj.f13.defval.Atr(ii).Val;
-                        userval = obj.f13.userval.Atr(ii).Val;
+                        userval = obj.f13.userval.Atr(jj).Val;
                         defval = reshape(defval,1,[]);
                         values = obj.p(:,1)*0 + defval;
                         values(userval(1,:),:) = userval(2:end,:)';

@msh/private/writefort20.m

@@ -0,0 +1,26 @@
+function [] = writefort20( fort20dat, finame )
+%
+%
+
+if ( nargin == 1 ) 
+  finputname = 'fort.20_1' ;
+else
+  finputname = finame ;   
+end
+
+fid = fopen(finputname,'w') ;
+
+fprintf( fid, '%d\n', fort20dat.FTIMINC ) ; 
+
+valpernode = size( fort20dat.Val,2 ) ;
+str = [];
+for ll = 1: valpernode
+    str = [str '%.8f '] ;
+end
+str = [str '\n' ] ;
+
+fprintf( fid, str, fort20dat.Val' ) ;
+
+fclose(fid) ; 
+%EOF
+end

Examples/Example_10_Multiscale_Smoother.m

@@ -1,6 +1,12 @@
 % Example_10_Multiscale_Smoother: 
 % An idealized test for multiscale nesting using boxes with a large min_el
 % ratio
+clearvars; clc;
+
+addpath('..')
+addpath(genpath('../utilities/'))
+addpath(genpath('../datasets/'))
+addpath(genpath('../m_map/'))
 
 bbox = [0, 1; 0,1]; 
 boubox = [0,0; 1,0; 1,1; 0,1; 0,0; NaN NaN ];

Examples/Example_11_Remeshing_Patches.m

@@ -9,9 +9,11 @@
 %
 % By Keith Roberts, 2020, USP, Brazil.
 clearvars; clc;
+addpath('..')
 addpath(genpath('../utilities/'))
 addpath(genpath('../datasets/'))
 addpath(genpath('../m_map/'))
+
 %% STEP 1: Here we mimic the steps that occurred in Example_1_NZ.m
 bbox = [166 176;		% lon_min lon_max
     -48 -40]; 		% lat_min lat_max

Examples/Example_12_Riverine_flow.m

@@ -0,0 +1,107 @@
+% Example_12_Riverine_flow: Guidance for creating a riverine flow input file 
+% otherwise referred to as a fort.20
+%
+% You should prepare a column delimited file (csv) that stores the total volume 
+% flow (Q) (m^3/s) time series of a cross-section where the riverine boundaries 
+% in the mesh are located. 
+%
+% A command-line method or a data cursor method is available to identify the 
+% vstart and vend of each riverine boundary in msh.make_bc.
+%
+% The columns of the csv file should be organized as follows:  
+% year,month,day,hour,minute,second,volume_flow_1,volume_flow_2... 
+% The column order of the total volume flow (volume_flow_1,volume_flow_2...) 
+% must be specified in the order in which the riverine boundaries appear in 
+% the fort.14 file, or in which you make the riverine boundaries with the data 
+% cursor method or comannd-line method in msh.make_bc.
+%
+% A 'test_make_f20.csv' file has been placed in the Tests folder for reference.
+% The volume flow of this file is an hourly average time series, thus DT equals 
+% to 3600 (s), if your volume flow is a daily average time series, DT should be  
+% set to 86400 (s) rather than 3600 (s).
+%
+% More details can be found in the two functions of Make_f20_volume_flow.m and 
+% Riverflux_distribution.m.
+%
+% Author:      Jiangchao Qiu                                
+% Created:     January 7 2021                                      
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+addpath('..')
+addpath(genpath('../utilities/'));
+addpath(genpath('../datasets/'));
+addpath(genpath('../m_map/'));
+addpath(genpath('../Tests/'));
+
+%% STEP 1: set mesh extents and set parameters for mesh. 
+%% The Pearl River Estuary region
+bbox = [112 20; 112 24; 116 24; 116 20; 112 20]; %polygon boubox
+min_el    = 100;  		       % minimum resolution in meters.
+max_el    = 10e3; 		       % maximum resolution in meters. 
+dt        = 5;                 % Automatically set timestep based on nearshore res
+grade     = 0.3;               % mesh grade in decimal percent. 
+R         = 3; 			       % Number of elements to resolve feature.
+
+%% STEP 2: specify geographical datasets and process the geographical data
+%% to be used later with other OceanMesh classes...
+dem       = 'SRTM15+V2.1.nc';
+coastline = 'GSHHS_f_L1';
+gdat = geodata('shp',coastline,'dem',dem,'h0',min_el,'bbox',bbox);
+%plot(gdat,'shp'); 
+%plot(gdat,'dem');  
+
+%% STEP 3: create an edge function class
+fh = edgefx('geodata',gdat,'fs',R,'max_el',max_el,'dt',dt,'g',grade);
+
+%% STEP 4: Pass your edgefx class object along with some meshing options and... 
+%% build the mesh         
+mshopts = meshgen('ef',fh,'bou',gdat,'plot_on',1,'proj','lam');
+mshopts = mshopts.build; 
+
+%% STEP 5: Plot and save the msh class object
+m = mshopts.grd; % get out the msh object
+plot(m,'type','tri');
+
+m = interp(m,gdat,'mindepth',1); % interpolate bathy to the mesh with minimum depth of 1 m
+m = lim_bathy_slope(m,0.1,0);
+CFL = CalcCFL(m,dt); 
+max(CFL)
+min(CFL)
+m = BoundCr(m,dt);
+
+%% STEP 6: Make the nodestring boundary conditions
+m = make_bc(m,'auto',gdat); 
+
+%% STEP 7: Identify vstart and vend of each riverine boundary
+% use the comannd-line method to identify the vstart and vend for each...
+% riverine boundary based on the coordinates of river edge.
+river_points1 = [112.918 22.508;112.919 22.514]; % coordinates of river edge
+river_points2 = [113.283 22.365;113.287 22.367]; % coordinates of river edge
+river_points3 = [113.513 22.234;113.515 22.236]; % coordinates of river edge            
+bc_k1 = ourKNNsearch(m.p',river_points1',1);
+bc_k2 = ourKNNsearch(m.p',river_points2',1);
+bc_k3 = ourKNNsearch(m.p',river_points3',1);
+m = make_bc(m,'outer',0,bc_k1(1),bc_k1(2),1,22);
+m = make_bc(m,'outer',1,bc_k2(1),bc_k2(2),1,22);
+m = make_bc(m,'outer',1,bc_k3(1),bc_k3(2),1,22);
+
+% use the data cursor method to identify the vstart and vend for each...
+% riverine boundary manually (GUI).
+% m = make_bc(m,'outer',1); % 1,2206,2107,1(flux),22(River)
+% m = make_bc(m,'outer',1); % 1,14784,14779,1(flux),22(River)
+% m = make_bc(m,'outer',1); % 1,763,704,1(flux),22(River)
+
+plot(m,'type','bd');  % plot mesh on native projection with boundary conditions
+plot(m,'type','b');   % plot bathy on native projection
+
+%% STEP 8: Make the riverine conditions
+ts = '01-Jan-2005 00:00'; % start time of the total volume flow time series
+te = '01-Jan-2005 23:00'; % end time of the total volume flow time series
+DT = 3600; % time step of the total volume flow time series,
+           % DT=3600 for hourly data series while DT=86400 for daily data series.
+m = Make_f20_volume_flow(m,'test_make_f20.csv',ts,te,DT);
+
+%% STEP 9: Save data
+save('Riverine.mat','m'); write(m,'Riverine','20');
+
+
+

Examples/Example_1_NZ.m

@@ -1,6 +1,7 @@
 % Example_1_NZ: Mesh the South Island of New Zealand
 clearvars; clc;
 
+addpath('..')
 addpath(genpath('../utilities/'))
 addpath(genpath('../datasets/'))
 addpath(genpath('../m_map/'))

Examples/Example_2_NY.m

@@ -1,6 +1,7 @@
 % Example_2_NY: Mesh the New York region in high resolution
 clearvars; clc;
 
+addpath('..')
 addpath(genpath('../utilities/'))
 addpath(genpath('../datasets/'))
 addpath(genpath('../m_map/'))

Examples/Example_3_ECGC.m

@@ -3,6 +3,7 @@
 
 clearvars; clc;
 
+addpath('..')
 addpath(genpath('../utilities/'));
 addpath(genpath('../datasets/'));
 addpath(genpath('../m_map/'));

Examples/Example_4_PRVI.m

@@ -4,6 +4,7 @@
 % around 1 hr. 
 clc; clearvars
 
+addpath('..')
 addpath(genpath('../utilities/'))
 addpath(genpath('../datasets/'))
 addpath(genpath('../m_map/'))

Examples/Example_5_JBAY.m

@@ -3,6 +3,7 @@
 
 clc; clearvars
 
+addpath('..')
 addpath(genpath('../utilities/'));
 addpath(genpath('../datasets/'));
 addpath(genpath('../m_map/')); 

Examples/Example_5b_JBAY_w_weirs.m

@@ -3,6 +3,7 @@
 % high resolution with two 15-30 m wide weirs at the mouth of the estuary.
 clc; clearvars
 
+addpath('..')
 addpath(genpath('../utilities/'));
 addpath(genpath('../datasets/'));
 addpath(genpath('../m_map/'));

Examples/Example_6_GBAY.m

@@ -3,6 +3,7 @@
 
 clearvars; clc;
 
+addpath('..')
 addpath(genpath('../utilities/'))
 addpath(genpath('../datasets/'))
 addpath(genpath('../m_map/'))

Examples/Example_6b_GBAY_w_floodplain.m

@@ -4,6 +4,7 @@
 
 clearvars; clc;
 
+addpath('..')
 addpath(genpath('../utilities/'))
 addpath(genpath('../datasets/'))
 addpath(genpath('../m_map/'))

Examples/Example_7_Global.m

@@ -2,6 +2,7 @@
 
 clearvars; clc;
 
+addpath('..')
 addpath(genpath('../utilities/'));
 addpath(genpath('../datasets/'));
 addpath(genpath('../m_map/'));

Examples/Example_8_AK.m

@@ -9,6 +9,7 @@
 
 clearvars; clc;
 
+addpath('..')
 addpath(genpath('../utilities/'))
 addpath(genpath('../datasets/'))
 addpath(genpath('../m_map/'))

Examples/Example_9_TBAY.m

@@ -22,6 +22,11 @@
 %% Setup
 clearvars; clc; close all;
 
+addpath('..')
+addpath(genpath('../utilities/'))
+addpath(genpath('../datasets/'))
+addpath(genpath('../m_map/'))
+
 %% Define data sources
 % Shoreline
 shoreline = 'GSHHS_f_L1';

README.md

@@ -135,6 +135,10 @@ The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.0.0/)
 
 ### Unreleased
 
+## Added
+- `mesh2d` interface improvements to filter small polygons. 
+- Support for creation of `fort.20` files for forcing rivers by @Jiangchao3
+
 ## Changed
 - `msh.plot()` overhaul. All options specified via kwarg.
 

Tests/TestSanity.m

@@ -4,26 +4,26 @@
 
 NP_TOL = 500;
 NT_TOL = 1500;
-QUAL_TOL = 0.15;
+QUAL_TOL = 0.5;
 % p: [5968�2 double]
 % t: [9530�3 double]
 % Element qual. metrics
 % MEAN      LOWER 3RD  MIN.
 % 0.9325    0.7452    0.3548
 if abs(length(m.p) - 5968) > NP_TOL
-    error(['Incorrect number of points for Example1_NZ.m. Got ',...
+    error(['Incorrect number of points for Example_1_NZ.m. Got ',...
         num2str(length(m.p)),' expecting 5968 +- 500 points']);
     exit(1)
 end
 if abs(length(m.t) - 9530) > NT_TOL
-    error(['Incorrect number of triangles for Example1_NZ.m. Got ',...
+    error(['Incorrect number of triangles for Example_1_NZ.m. Got ',...
         num2str(length(m.t)),' expecting 9530 +- 1500 elements']);
     exit(1)
 end
-if abs(mshopts.qual(end,3) - 0.3548) > QUAL_TOL
-    error(['Incorrect min. element quality for Example1_NZ.m. Got ',...
-        num2str(mshopts.qual(end,3)),' expecting 0.3548 +- 0.05']);
+if abs(mshopts.qual(end,3)) < QUAL_TOL
+    error(['Incorrect min. element quality for Example_1_NZ.m. Got ',...
+        num2str(mshopts.qual(end,3)),' expecting > 0.5']);
     exit(1)
 end
-disp('Passed: Example1_NZ.m');
+disp('Passed: Example_1_NZ.m');
 

Tests/test_make_f20.csv

@@ -0,0 +1,48 @@
+2005,1,19,0,0,0,1780,1960,534
+,,19,1,0,0,1960,1960,734
+,,19,2,0,0,1960,1960,794
+,,19,3,0,0,1610,1960,872
+,,19,4,0,0,549,1960,802
+,,19,5,0,0,-641,1960,714
+,,19,6,0,0,-1230,1960,401
+,,19,7,0,0,-1010,1960,-216
+,,19,8,0,0,-157,1960,-573
+,,19,9,0,0,350,1960,-531
+,,19,10,0,0,741,1960,-356
+,,19,11,0,0,912,1960,38.1
+,,19,12,0,0,934,1960,350
+,,19,13,0,0,752,1960,492
+,,19,14,0,0,199,1960,581
+,,19,15,0,0,-563,1960,571
+,,19,16,0,0,-1500,1960,401
+,,19,17,0,0,-1960,1960,-81.3
+,,19,18,0,0,-2110,1960,-667
+,,19,19,0,0,-1880,1960,-944
+,,19,20,0,0,-1280,1960,-999
+,,19,21,0,0,-279,1960,-995
+,,19,22,0,0,106,1960,-792
+,,19,23,0,0,1310,1960,-355
+2005,1,20,0,0,0,1910,1470,231
+,,20,1,0,0,1940,1470,528
+,,20,2,0,0,2100,1470,671
+,,20,3,0,0,1860,1470,706
+,,20,4,0,0,1580,1470,689
+,,20,5,0,0,1210,1470,696
+,,20,6,0,0,87.8,1470,702
+,,20,7,0,0,-638,1470,664
+,,20,8,0,0,-1200,1470,230
+,,20,9,0,0,-1370,1470,-261
+,,20,10,0,0,-1210,1470,-608
+,,20,11,0,0,-792,1470,-722
+,,20,12,0,0,-307,1470,-790
+,,20,13,0,0,96.9,1470,-645
+,,20,14,0,0,356,1470,-74.5
+,,20,15,0,0,98.3,1470,302
+,,20,16,0,0,-710,1470,523
+,,20,17,0,0,-1400,1470,347
+,,20,18,0,0,-2150,1470,-381
+,,20,19,0,0,-1960,1470,-865
+,,20,20,0,0,-1580,1470,-1040
+,,20,21,0,0,-1180,1470,-1020
+,,20,22,0,0,-632,1470,-896
+,,20,23,0,0,112,1470,-708