Some fixes for the examples

- Now all examples include the makens function and the interp function at the end (except for the NZ example which doesn't include bathy)
- Other various fixes to ensure everything works and plots smoothly and is clean

Author: WPringle

.gitignore

@@ -1,4 +1,4 @@
 m_map/
 datasets/
 fort.*
-./*.mat
+/*.mat

@meshgen/meshgen.m

@@ -307,7 +307,7 @@
                     obj.bbox(1,2) obj.bbox(2,1); ...
                     obj.bbox(1,1) obj.bbox(2,1); NaN NaN];
             end
-            if any(obj.h0==0), error('h0 was not correctly specified!'), end;
+            if any(obj.h0==0), error('h0 was not correctly specified!'), end
             if isempty(obj.outer), error('no outer boundary specified!'), end
             if isempty(obj.bbox), error('no bounding box specified!'), end
             obj.fd = @dpoly;  % <-default distance fx accepts p and pv (outer polygon).
@@ -681,6 +681,13 @@
                 % Perform the direct smoothing
                 [obj.grd.p,obj.grd.t] = direct_smoother_lur(obj.grd.p,...
                     obj.grd.t,obj.pfix,obj.nscreen);
+                tq = gettrimeshquan( obj.grd.p, obj.grd.t);
+                if min(tq.qm) < 0
+                    % Need to clean it again
+                    disp('Overlapping elements due to smoother, cleaning again')
+                    obj.grd = Make_Mesh_Boundaries_Traversable(...
+                                       obj.grd,obj.dj_cutoff,obj.nscreen);
+                end
             end
 
             % Checking and displaying element quality

@msh/msh.m

@@ -203,7 +203,7 @@ function write(obj,fname,type)
                     end
                 case('b')
                     figure;
-                    cmap = cmocean('deep');
+                    cmap = cmocean('deep',numticks-1);
                     if logaxis
                         q = log10(max(1,abs(obj.b))); % plot on log scale
                     else
@@ -215,9 +215,9 @@ function write(obj,fname,type)
                         trisurf(obj.t,obj.p(:,1),obj.p(:,2),q);
                         view(2); shading interp;
                     end
-                    cmocean('deep',numticks-1); cb = colorbar;
+                    cb = colorbar;
                     if logaxis
-                        desiredTicks = round(10.^(linspace(min(q),max(q),numticks)));
+                        desiredTicks = round(10.^(linspace(min(q),max(q),numticks)),1);
                         caxis([log10(min(desiredTicks)) log10(max(desiredTicks))]);
                         cb.Ticks     = log10(desiredTicks);
                         for i = 1 : length(desiredTicks)
@@ -508,6 +508,10 @@ function write(obj,fname,type)
             if nargin < 2
                error('Needs type: one of auto, islands or outer') 
             end
+            trim = 0;
+            if strcmp(type(max(1,end-3):end),'trim')
+                type = type(1:end-4); trim = 1;
+            end
             switch type
                 case('auto')
                     if ~isa(dir,'geodata')
@@ -520,11 +524,7 @@ function write(obj,fname,type)
                     [~,poly_idx] = extdom_polygon(etbv,obj.p,1);
 
                     % Get geodata outer and mainland polygons
-                    outer = [gdat.outer(~isnan(gdat.outer(:,1)),:);
-                             gdat.inner(~isnan(gdat.inner(:,1)),:)];
-                    main = [gdat.mainland(~isnan(gdat.mainland(:,1)),:);
-                            gdat.inner(~isnan(gdat.inner(:,1)),:)];
-                    
+                    outer = gdat.outer(~isnan(gdat.outer(:,1)),:);
                     nope = 0; neta = 0; nbou  = 0; nvel  = 0;
                     % Find the polygon that is a combination of ocean
                     % and mainland boundaries.
@@ -537,13 +537,23 @@ function write(obj,fname,type)
                             break;
                         end
                     end
+                    % Get geodata outer and mainland polygons
+                    outer = [gdat.outer(~isnan(gdat.outer(:,1)),:);
+                             gdat.inner(~isnan(gdat.inner(:,1)),:)];
+                    main = [gdat.mainland(~isnan(gdat.mainland(:,1)),:);
+                            gdat.inner(~isnan(gdat.inner(:,1)),:)];
+                    [~, d_out] = ourKNNsearch(outer',obj.p(idv,:)',1);    
                     [~, d_main] = ourKNNsearch(main',obj.p(idv,:)',1);
 
                     % Mainland are nodes where shortest distance to
                     % mainland and outer are the same and where absolute
                     % distance to mainland is relatively small
-                    mainland = abs(d_out - d_main) < gdat.h0/111e3 & ...
-                               d_main < 5*gdat.h0/111e3;
+                    if trim
+                        mainland = abs(d_out - d_main) < gdat.h0/111e3 & ...
+                                   d_main < 5*gdat.h0/111e3;
+                    else
+                        mainland = abs(d_out - d_main) < gdat.h0/111e3; 
+                    end
 
                     % indices of switch
                     Cuts  = find(diff(mainland ~= 0));

Example_1_NZ.m

@@ -1,4 +1,5 @@
 % Example_1_NZ: Mesh the South Island of New Zealand
+clearvars; clc;
 addpath(genpath('utilities/'))
 addpath(genpath('datasets/'))
 addpath(genpath('m_map/'))
@@ -20,8 +21,12 @@
             'max_el',max_el,'g',grade);
 %% STEP 4: Pass your edgefx class object along with some meshing options and
 % build the mesh...
-mshopt = meshgen('ef',{fh},'bou',{gdat},'nscreen',5,'plot_on',1,'itmax',75);
-mshopt = mshopt.build; 
+mshopts = meshgen('ef',fh,'bou',gdat,'plot_on',1);
+mshopts = mshopts.build; 
+
 %% STEP 5: Plot it and write a triangulation fort.14 compliant file to disk.
-plot(mshopt.grd,'tri');
-write(mshopt.grd,'South_Island_NZ');
+% Get out the msh class and put on nodestrings
+m = mshopts.grd;
+m = makens(m,'auto',gdat); % make the nodestring boundary conditions
+plot(m,'bd');
+write(m,'South_Island_NZ');

Example_2_NY.m

@@ -1,10 +1,9 @@
 % Example_2_NY: Mesh the New York region in high resolution
+clearvars; clc;
 addpath(genpath('utilities/'))
 addpath(genpath('datasets/'))
 addpath(genpath('m_map/'))
 %% STEP 1: Set mesh extents and set parameters for mesh.
-bbox = [-74.5 -73.8 	% lon_min lon_max
-        40.5 40.9]; 	% lat_min lat_max
 min_el    = 30;  	% Minimum resolution in meters.
 max_el    = 1e3; 	% Maximum resolution in meters. 
 max_el_ns = 240;        % Maximum resolution nearshore in meters.
@@ -17,16 +16,19 @@
 dem       = 'PostSandyNCEI.nc';
 gdat = geodata('shp',coastline,...
                'dem',dem,...
-               'bbox',bbox,...
                'h0',min_el); 
 %% STEP 3: Create an edge function class.
 fh = edgefx('geodata',gdat,...
             'fs',R,'max_el_ns',max_el_ns,...
             '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},'nscreen',1,'plot_on',1,'itmax',150);
+mshopts = meshgen('ef',fh,'bou',gdat,'plot_on',1);
 mshopts = mshopts.build; 
 %% STEP 5: Plot it and write a triangulation fort.14 compliant file to disk.
-plot(mshopts.grd,'tri');
-write(mshopts.grd,'NY_HR');
+% Get out the msh class and put on bathy and nodestrings
+m = mshopts.grd;
+m = interp(m,gdat); m.b = max(m.b,1); % interpolate bathy to the mesh
+m = makens(m,'auto',gdat); % make the nodestring boundary conditions
+plot(m,'bd'); plot(m,'blog'); % plot triangulation and bathy
+write(m,'NY_HR');

Example_3_ECGC.m

@@ -7,7 +7,7 @@
 
 %% STEP 1: set mesh extents and set parameters for mesh. 
 %% The greater US East Coast and Gulf of Mexico region
-bbox      = [-100 -50; 10  60]; % lon min lon max; lat min lat max
+bbox      = [-100 -50; 5  55];  % lon min lon max; lat min lat max
 min_el    = 1e3;  		        % minimum resolution in meters.
 max_el    = inf; 		        % maximum resolution in meters. 
 wl        = 60;                 % 60 elements resolve M2 wavelength.
@@ -44,12 +44,13 @@
 
 %% STEP 4: Pass your edgefx class object along with some meshing options 
 %% and build the mesh...
-mshopts = meshgen('ef',{fh1 fh2},'bou',{gdat1 gdat2},...
-                 'nscreen',1,'plot_on',1,'itmax',50);
+mshopts = meshgen('ef',{fh1 fh2},'bou',{gdat1 gdat2},'plot_on',1);
 mshopts = mshopts.build; 
 
 %% Plot and save the msh class object/write to fort.14
 m = mshopts.grd; % get out the msh object
+m = interp(m,{gdat1 gdat2}); m.b = max(m.b,1); % interpolate bathy to the mesh
 m = makens(m,'auto',gdat1); % make the nodestring boundary conditions
 plot(m,'bd',1,'Mollweide'); % plot on Mollweide projection with nodestrings
+plot(m,'b',1,'Mollweide'); % plot bathy on Mollweide projection
 save('ECGC_w_NYHR.mat','m'); write(m,'ECGC_w_NYHR');

Example_4_PRVI.m

@@ -62,29 +62,28 @@
 %% Pass your edgefx class objects along with some meshing options 
 %% and build the mesh... 
 % (note that the nested edgefxs will be smoothed together with this call)
-mshopts = meshgen('ef',fh,'bou',gdat,'nscreen',5,'plot_on',1,'itmax',50);  
+mshopts = meshgen('ef',fh,'bou',gdat,'plot_on',1,'itmax',50);  
 
 % now build the mesh with your options and the edge function.
 mshopts = mshopts.build; 
 
-%% Get out the msh class from meshgen
+% Get out the msh class from meshgen
 m = mshopts.grd;
 
-%plot(m,'tri'); 
-
-%plot(m,'reso');
-
 %% Interpolate on the bathy and gradients (automatically loops over all data)
-% m = interp(m,gdat); 
+m = interp(m,gdat); 
 % % ensure max depth in domain is 1 m (we find this step useful for coastal
 % % meshes to help with connectivity through narrow channels)
-% m.b = max(m.b,1); 
-% 
-% plot(m,'b');
+m.b = max(m.b,1); 
 
+%% Make the nodestrings
+m = makens(m,'auto',gdat{1}); % make the nodestring boundary conditions
+
+%% Plot and save the msh class object/write to fort.14
+plot(m,'bd',1,'Sinusoidal'); % plot on Sinusoidal projection with nodestrings
+plot(m,'b',1,'Sinusoidal'); % plot the bathy on Sinusoidal projection
 % Save as a msh class
 save('PRVI_msh.mat','m');
-
 % Write an ADCIRC fort.14 compliant file to disk.
-%write(m,'PRVI_mesh')
+write(m,'PRVI_mesh')
 

Example_5_JBAY.m

@@ -1,4 +1,8 @@
-%Example_2_NY: Mesh the New York Jamaica bay (JBAY region in high resolution.
+% Example_5_JBAY: Mesh the New York Jamaica bay (JBAY) region in 
+% high resolution.
+
+clc; clearvars
+
 addpath(genpath('utilities/'));
 addpath(genpath('datasets/'));
 addpath(genpath('m_map/')); 
@@ -26,8 +30,13 @@
             '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);
+mshopts = meshgen('ef',fh,'bou',gdat,'plot_on',1);
 % now build the mesh with your options and the edge function.
 mshopts = mshopts.build; 
-%% STEP 5: Plot it 
-plot(mshopts.grd,'tri');
+%% STEP 5: Plot it and save the msh file
+% Get out the msh class and put on bathy and nodestrings
+m = mshopts.grd;
+m = interp(m,gdat); m.b = max(m.b,1); % interpolate bathy to the mesh
+m = makens(m,'auto',gdat); % make the nodestring boundary conditions
+plot(m,'bd'); plot(m,'blog'); % plot triangulation and bathy
+save('JBAY_HR.mat','m')

Example_6_GBAY.m

@@ -1,6 +1,12 @@
+% Example_6_GBAY: Mesh the Galveston bay (GBAY) region in 
+% high resolution. 
+
+clearvars; clc;
+
 addpath(genpath('utilities/'))
 addpath(genpath('datasets/'))
 addpath(genpath('m_map/')) 
+
 %% STEP 1: set mesh extents and set parameters for mesh.
 bbox = [-95.40 -94.4;  
          29.14  30.09];
@@ -14,7 +20,7 @@
                'dem',demfile,...
                'bbox',bbox,...
                'h0',min_el);
-load ECGC_Thalwegs.mat
+load ECGC_Thalwegs.mat % Load the Channel thalweg data
 %% STEP 3: create an edge function class
 fh = edgefx('geodata',gdat,...
             'fs',3,...
@@ -28,5 +34,8 @@
 % now build the mesh with your options and the edge function.
 mshopts = mshopts.build; 
 %% STEP 5: Plot it and write a triangulation fort.14 compliant file to disk.
-plot(mshopts.grd,'tri',0);
-write(mshopts.grd,'HoustonShipChannel');
+m = mshopts.grd;
+m = interp(m,gdat); m.b = max(m.b,1); % interpolate bathy to the mesh
+m = makens(m,'auto',gdat); % make the nodestring boundary conditions
+plot(m,'bd'); plot(m,'blog'); % plot triangulation and bathy
+write(m,'HoustonShipChannel');

utilities/m_trisurf.m

@@ -6,14 +6,13 @@
 
 % Have to have initialized a map first
 
-if isempty(MAP_PROJECTION),
+if isempty(MAP_PROJECTION)
     disp('No Map Projection initialized - call M_PROJ first!');
     return;
-end;
-m_grid('linest','-');
+end
 
 [X,Y]=m_ll2xy(long,lat,'clip','on');
 colormap(cmap); 
-hold on; h=trisurf(tri,X,Y,z,'facecolor', 'flat', 'edgecolor', 'none');
-shading interp;
+hold on; h=trisurf(tri,X,Y,z,'facecolor', 'interp', 'edgecolor', 'none');
+m_grid('linest','-');
 end