From 341eadb5f1ea3555375dd51aeca6ea2b5bd25eb0 Mon Sep 17 00:00:00 2001
From: cblakely97 <47362791+cblakely97@users.noreply.github.com>
Date: Wed, 2 Sep 2020 14:21:38 -0400
Subject: [PATCH 1/9] Smallest Edgefunction Value

I added these two functions to handle a situation where an edgefunction placed "on top" (in the original code takes precedence) has a larger target edgelength than an edgefunction "below." This simply loops through the edgefunctions and places the smallest value "on top" so that it takes precedence during the mesh generation process.
---
 utilities/edgefx2pts.m    |  13 ++++
 utilities/edgefx_filter.m | 141 ++++++++++++++++++++++++++++++++++++++
 2 files changed, 154 insertions(+)
 create mode 100644 utilities/edgefx2pts.m
 create mode 100644 utilities/edgefx_filter.m

diff --git a/utilities/edgefx2pts.m b/utilities/edgefx2pts.m
new file mode 100644
index 00000000..673edbbc
--- /dev/null
+++ b/utilities/edgefx2pts.m
@@ -0,0 +1,13 @@
+function [pts, vals] = edgefx2pts(fh)
+% takes an edgefunction from OceanMesh2D and outputs a series of points and
+% values for use in fastscatter
+% Coleman Blakely
+% January 22, 2020
+x = fh.F.GridVectors{1};
+y = fh.F.GridVectors{2}';
+val = fh.F.Values;
+vals = reshape(val,[],1);
+[lat, lon] = meshgrid(x,y);
+c = cat(2,lat',lon');
+pts = reshape(c,[],2);
+end
\ No newline at end of file
diff --git a/utilities/edgefx_filter.m b/utilities/edgefx_filter.m
new file mode 100644
index 00000000..d5981ac5
--- /dev/null
+++ b/utilities/edgefx_filter.m
@@ -0,0 +1,141 @@
+function fh = edgefx_filter(fh,gdat)
+% This function sorts through layered edgefunctions and puts the minimum
+% target edgelength in the topmost edgefunction. This ensures that when
+% the mesh is generated, the minimum target edgelength from all boxes is
+% used.
+% Coleman Blakely
+% January 22, 2020
+if length(fh) == 1
+    disp('No need to do this silly')
+end
+% preallocate polyshape vector for finding whether boundaries are
+% overlapping
+S = cell(length(fh)-1,1);
+% build polyshapes starting with last (top) edgefunction and working down
+for kk = length(fh):-1:2
+    S{kk} = polyshape(fh{kk}.boubox);
+end
+% Loop through edgefunctions to 1) check if they overlap and 2) if they
+% overlap find the points that lie within the upper edgefunction and
+% compare them to the ones below
+for kk = length(fh):-1:2
+    fh_test = fh{kk};
+    % make the points to use with inpoly to find the overlapping points
+    [test_pts, ~] = edgefx2pts(fh_test);
+    test_vals = fh_test.F.Values;
+    x_test = fh_test.F.GridVectors{1};
+    y_test = fh_test.F.GridVectors{2};
+    [X_test, Y_test] = ndgrid(x_test,y_test);
+    % Loop through all boxes that may be below the test box
+    for ii = kk-1:-1:2
+        % check if the test box overlaps the below box
+        if ~overlaps(S{kk},S{ii})
+            continue
+        end
+        fh_lower = fh{ii};
+        % find the test points that are within the lower boundary
+        id_overlap = inpoly(test_pts,fh_lower.boubox(1:end-1,:));
+        overlap_pts = test_pts(id_overlap,:);
+        % create the grid of overlapping points to use in the gridded
+        % interpolant
+        x = unique(overlap_pts(:,1));
+        y = unique(overlap_pts(:,2));
+        [X, Y] = ndgrid(x,y);
+        % evaluate the gridded interpolant of the lower box using the
+        % points that are within the test box
+        lower_vals = fh_lower.F(X,Y);
+        % find the test points within the test box and compare with the
+        % newly found points
+        [~,idx,~] = intersect(x_test,x);
+        [~,idy,~] = intersect(y_test,y);
+        % int_vals is a dummy variable to make the indexing of the boxes
+        % the same for comparison
+        int_vals = test_vals(idx,idy);
+        % replace all test values that are greater than the lower values
+        int_vals(int_vals>lower_vals) = lower_vals(int_vals>lower_vals);
+        % place the int_vals back into the test values
+        test_vals(idx,idy) = int_vals;
+    end
+    % put the new points back in the edgefunction
+    fh{kk}.F = griddedInterpolant(X_test,Y_test,test_vals);
+    
+    
+    
+end
+% now resmooth the edgefunctions using the defined grade to eliminate
+% any discontinuities, etc.
+% most of this was pulled from stuff that William and Keith did. I do
+% not know exactly how all of it works but it appears to do the trick.
+% CPB
+for kk = length(fh):-1:1
+    obj = fh{kk};
+    feat = gdat{kk};
+    hh_m = obj.F.Values;
+    [xg, yg] = CreateStructGrid(obj);
+    if all(abs(obj.bbox(1,:)) == 180)
+        % for global mesh make it cyclical
+        hh_m = [hh_m(end,:); hh_m; hh_m(1,:)];
+    end
+    hfun = reshape(hh_m',[numel(hh_m),1]);
+    
+    dx = obj.h0*cosd(min(yg(1,:),85)); % for gradient function
+    dy = obj.h0;                       % for gradient function
+    
+    % make g a function of space
+    dmy     = xg*0 ;
+    for param = obj.g'
+        if numel(param)==1 && param~=0
+            lim   = obj.g(1);
+            dmy  = dmy + lim ;
+        else
+            lim  = param(1);
+            dp1  = param(2);
+            dp2  = param(3);
+            
+            limidx = (feat.Fb(xg,yg) < dp1 & ...
+                feat.Fb(xg,yg) > dp2) ;
+            
+            dmy( limidx ) = lim;
+        end
+    end
+    if all(abs(obj.bbox(1,:)) == 180)
+        % for global mesh make it cyclical
+        dmy = [dmy(end,:); dmy; dmy(1,:)];
+    end
+    fdfdx = reshape(dmy',[numel(dmy),1]);
+    clearvars dmy;
+    [hfun,flag] = limgradStruct(obj.ny,dx,dy,hfun,...
+        fdfdx,sqrt(length(hfun)));
+    if flag == 1
+        disp('Gradient relaxing converged!');
+    else
+        error(['FATAL: Gradient relaxing did not converge, '
+            'please check your edge functions']);
+    end
+    % reshape it back
+    hh_m = reshape(hfun,obj.ny,[])';
+    clearvars hfun fdfdx
+    if all(abs(obj.bbox(1,:)) == 180)
+        % for global mesh make it cyclical
+        hh_m = hh_m(2:end-1,:);
+        hh_m(end,:) = hh_m(1,:);
+    end
+    obj.F = griddedInterpolant(xg,yg,hh_m,'linear','nearest');
+    
+    clearvars xg yg
+    fh{kk} = obj;
+end
+end
+
+%% This is a quick little function to turn an edgefuncton into points and values that can be used with fastscatter
+function [pts, vals] = edgefx2pts(fh)
+% takes an edgefunction from OceanMesh2D and outputs a series of points and
+% minimum element sizes for each point
+x = fh.F.GridVectors{1};
+y = fh.F.GridVectors{2}';
+val = fh.F.Values;
+vals = reshape(val,[],1);
+[lat, lon] = meshgrid(x,y);
+c = cat(2,lat',lon');
+pts = reshape(c,[],2);
+end
\ No newline at end of file

From 34efeab97ab8fe478e7571f7656cc100dfb57d68 Mon Sep 17 00:00:00 2001
From: cblakely97 <47362791+cblakely97@users.noreply.github.com>
Date: Fri, 4 Sep 2020 15:15:36 -0400
Subject: [PATCH 2/9] Update edgefx_filter.m

Removed the use of polyshapes and overlaps functions in lieu of using inpoly. Created a simple example of the use of edgefx_filter.m in generating a mesh.
---
 Examples/Example_9_edgefx_filter.m | 62 ++++++++++++++++++++++++++++++
 utilities/edgefx_filter.m          | 20 ++++------
 2 files changed, 70 insertions(+), 12 deletions(-)
 create mode 100644 Examples/Example_9_edgefx_filter.m

diff --git a/Examples/Example_9_edgefx_filter.m b/Examples/Example_9_edgefx_filter.m
new file mode 100644
index 00000000..2a2eefe3
--- /dev/null
+++ b/Examples/Example_9_edgefx_filter.m
@@ -0,0 +1,62 @@
+% Example_9_edgefx_filter.m
+% 
+%   This script serves as an example of how the edgefx_filter function can
+% ensure that the smallest possible target edgelength function in a series
+% of edgefunctions is carried through to the innermost edgefunction before
+% being used to generate a mesh in meshgen. The reason that this function
+% can sometimes be necessary is because in different locations, different
+% controlling edgefunction parameters can result in larger target
+% edgelengths being put "above" or "inside" the edgefunction that meshgen
+% utilizes. A simple example, and one used here, is if the channels
+% parameter is utilized in an outer edgefunction but not an inner
+% edgefunction. As can be seen by the two meshes generated using this
+% script, the edgefx_filter fixes this problem and ensures the channels are
+% still incorporated in the final mesh.
+%
+% Coleman Blakely
+% September 4, 2020
+%% Setup
+clear;clc;hold off;close all;
+%% Add libraries
+addpath(genpath('C:\Users\Diane\Documents\GitHub\OceanMesh2D\'))
+%% Define data sources
+% Shoreline
+shoreline = 'GSHHS_f_L1';
+% Bathy data
+dem = 'SRTM15+V2.nc';
+%% Define meshing domain and parameters
+bbox{1} = [-83, -82; 27 28.5];        % outermost bbox
+min_el = 50;                       % minimum edgelength
+max_el = 1e3;                       % maximum edgelength
+grade = 0.35;                       % maximum change in size of edgelength
+R = 3;                              % elements required to resolve a feature
+load ECGC_Thalwegs.mat
+%% Build geodata and edgelength function for outermost domain
+gdat{1} = geodata('shp',shoreline,'dem',dem,'bbox',bbox{1},'h0',min_el);
+fh{1} = edgefx('geodata',gdat{1},'fs',R,'max_el',max_el,'g',grade,'Ch',1,'Channels',pts2);
+%% Build inner domain
+% To show how the edgefunction filter works, build another edgefunction
+% within the outer domain with a larger minimum edgelength
+% min_el = 150;
+bbox{2} = [-82.8, -82.4; 27.25 28.25];
+gdat{2} = geodata('shp',shoreline,'dem',dem,'bbox',bbox{2},'h0',min_el);
+fh{2} = edgefx('geodata',gdat{2},'fs',R,'max_el',max_el,'g',grade);
+% plotedgefx(fh,1);
+% bound = bbox2bound(bbox{2});
+% plot(bound(:,1),bound(:,2),'r-')
+%% Generate the mesh using the unfiltered edgefunctions
+m_unfiltered = meshgen('ef',fh,'bou',gdat,'plot_on',0);
+m_unfiltered = m_unfiltered.build;
+m_unfiltered = m_unfiltered.grd;
+%% Use the edgefunction filter
+fh = edgefx_filter(fh,gdat);
+% plotedgefx(fh,2)
+%% Generate the mesh using the filtered edgefunction
+m_filtered = meshgen('ef',fh,'bou',gdat,'plot_on',0);
+m_filtered = m_filtered.build;
+m_filtered = m_filtered.grd;
+%% Plot the two meshes
+plot(m_unfiltered,'tri',0)
+title('Unfiltered Edgefunction')
+plot(m_filtered,'tri',0)
+title('Filtered Edgefunction')
\ No newline at end of file
diff --git a/utilities/edgefx_filter.m b/utilities/edgefx_filter.m
index d5981ac5..44fdfbd0 100644
--- a/utilities/edgefx_filter.m
+++ b/utilities/edgefx_filter.m
@@ -5,16 +5,11 @@
 % used.
 % Coleman Blakely
 % January 22, 2020
+% Edited September 4, 2020 to remove the use of shapevectors and instead
+% use inpoly.m
 if length(fh) == 1
     disp('No need to do this silly')
 end
-% preallocate polyshape vector for finding whether boundaries are
-% overlapping
-S = cell(length(fh)-1,1);
-% build polyshapes starting with last (top) edgefunction and working down
-for kk = length(fh):-1:2
-    S{kk} = polyshape(fh{kk}.boubox);
-end
 % Loop through edgefunctions to 1) check if they overlap and 2) if they
 % overlap find the points that lie within the upper edgefunction and
 % compare them to the ones below
@@ -27,14 +22,15 @@
     y_test = fh_test.F.GridVectors{2};
     [X_test, Y_test] = ndgrid(x_test,y_test);
     % Loop through all boxes that may be below the test box
-    for ii = kk-1:-1:2
-        % check if the test box overlaps the below box
-        if ~overlaps(S{kk},S{ii})
-            continue
-        end
+    for ii = kk-1:-1:1
         fh_lower = fh{ii};
         % find the test points that are within the lower boundary
         id_overlap = inpoly(test_pts,fh_lower.boubox(1:end-1,:));
+        % check to see if the test box overlaps the lower box. If not, skip
+        % this box
+        if sum(id_overlap)==0
+           continue 
+        end
         overlap_pts = test_pts(id_overlap,:);
         % create the grid of overlapping points to use in the gridded
         % interpolant

From ec6ab9d5138698e1704229743970837eef7ae99e Mon Sep 17 00:00:00 2001
From: Keith Roberts <krober10@nd.edu>
Date: Sat, 5 Sep 2020 15:25:51 -0300
Subject: [PATCH 3/9] Moving CBlakely's changes to meshgen

---
 @meshgen/meshgen.m                           |   3 +
 {utilities => @meshgen/private}/edgefx2pts.m |   0
 @meshgen/private/enforce_min_ef.m            |  56 ++++++++
 Examples/Example_9_edgefx_filter.m           |  62 ---------
 utilities/edgefx_filter.m                    | 137 -------------------
 5 files changed, 59 insertions(+), 199 deletions(-)
 rename {utilities => @meshgen/private}/edgefx2pts.m (100%)
 create mode 100644 @meshgen/private/enforce_min_ef.m
 delete mode 100644 Examples/Example_9_edgefx_filter.m
 delete mode 100644 utilities/edgefx_filter.m

diff --git a/@meshgen/meshgen.m b/@meshgen/meshgen.m
index f435f6c8..81c5eb5c 100644
--- a/@meshgen/meshgen.m
+++ b/@meshgen/meshgen.m
@@ -290,6 +290,9 @@
                         
                         % kjr 2018 smooth the outer automatically
                         if length(obj.ef) > 1
+                            % kjr 2020, ensure the min. sizing func is
+                            % used
+                            obj.ef = enforce_min(obj.ef); 
                             obj.ef = smooth_outer(obj.ef,obj.Fb);
                         end
                         
diff --git a/utilities/edgefx2pts.m b/@meshgen/private/edgefx2pts.m
similarity index 100%
rename from utilities/edgefx2pts.m
rename to @meshgen/private/edgefx2pts.m
diff --git a/@meshgen/private/enforce_min_ef.m b/@meshgen/private/enforce_min_ef.m
new file mode 100644
index 00000000..4565200c
--- /dev/null
+++ b/@meshgen/private/enforce_min_ef.m
@@ -0,0 +1,56 @@
+function fh = enforce_min_ef(fh)
+% Ensures that the minimum mesh sizing value is always used in the case 
+% of many edgefxs (e.g., multiscale meshing technique).
+%
+% Loop through edgefunctions to check if they overlap. 
+% If they overlap, find the points that lie within the 
+% upper edgefunction and compare them to the ones below.
+%
+% By Coleman Blakely UND, 2020.
+%   Modifications by Keith Roberts USP, 2020. 
+
+for kk = length(fh):-1:2
+    fh_test = fh{kk};
+    % make the points to use with inpoly to find the overlapping points
+    [test_pts, ~] = edgefx2pts(fh_test);
+    test_vals = fh_test.F.Values;
+    x_test = fh_test.F.GridVectors{1};
+    y_test = fh_test.F.GridVectors{2};
+    [X_test, Y_test] = ndgrid(x_test,y_test);
+    % Loop through all boxes that may be below the test box
+    for ii = kk-1:-1:1
+        fh_lower = fh{ii};
+        % find the test points that are within the lower boundary
+        id_overlap = inpoly(test_pts,fh_lower.boubox(1:end-1,:));
+        % check to see if the test box overlaps the lower box. If not, skip
+        % this box
+        if sum(id_overlap)==0
+            continue
+        end
+        overlap_pts = test_pts(id_overlap,:);
+        % create the grid of overlapping points to use in the gridded
+        % interpolant
+        x = unique(overlap_pts(:,1));
+        y = unique(overlap_pts(:,2));
+        [X, Y] = ndgrid(x,y);
+        % evaluate the gridded interpolant of the lower box using the
+        % points that are within the test box
+        lower_vals = fh_lower.F(X,Y);
+        % find the test points within the test box and compare with the
+        % newly found points
+        [~,idx,~] = intersect(x_test,x);
+        [~,idy,~] = intersect(y_test,y);
+        % int_vals is a dummy variable to make the indexing of the boxes
+        % the same for comparison
+        int_vals = test_vals(idx,idy);
+        % replace all test values that are greater than the lower values
+        int_vals(int_vals>lower_vals) = lower_vals(int_vals>lower_vals);
+        % place the int_vals back into the test values
+        test_vals(idx,idy) = int_vals;
+    end
+    % put the new points back in the edgefunction
+    fh{kk}.F = griddedInterpolant(X_test,Y_test,test_vals);
+    
+    
+    
+end
diff --git a/Examples/Example_9_edgefx_filter.m b/Examples/Example_9_edgefx_filter.m
deleted file mode 100644
index 2a2eefe3..00000000
--- a/Examples/Example_9_edgefx_filter.m
+++ /dev/null
@@ -1,62 +0,0 @@
-% Example_9_edgefx_filter.m
-% 
-%   This script serves as an example of how the edgefx_filter function can
-% ensure that the smallest possible target edgelength function in a series
-% of edgefunctions is carried through to the innermost edgefunction before
-% being used to generate a mesh in meshgen. The reason that this function
-% can sometimes be necessary is because in different locations, different
-% controlling edgefunction parameters can result in larger target
-% edgelengths being put "above" or "inside" the edgefunction that meshgen
-% utilizes. A simple example, and one used here, is if the channels
-% parameter is utilized in an outer edgefunction but not an inner
-% edgefunction. As can be seen by the two meshes generated using this
-% script, the edgefx_filter fixes this problem and ensures the channels are
-% still incorporated in the final mesh.
-%
-% Coleman Blakely
-% September 4, 2020
-%% Setup
-clear;clc;hold off;close all;
-%% Add libraries
-addpath(genpath('C:\Users\Diane\Documents\GitHub\OceanMesh2D\'))
-%% Define data sources
-% Shoreline
-shoreline = 'GSHHS_f_L1';
-% Bathy data
-dem = 'SRTM15+V2.nc';
-%% Define meshing domain and parameters
-bbox{1} = [-83, -82; 27 28.5];        % outermost bbox
-min_el = 50;                       % minimum edgelength
-max_el = 1e3;                       % maximum edgelength
-grade = 0.35;                       % maximum change in size of edgelength
-R = 3;                              % elements required to resolve a feature
-load ECGC_Thalwegs.mat
-%% Build geodata and edgelength function for outermost domain
-gdat{1} = geodata('shp',shoreline,'dem',dem,'bbox',bbox{1},'h0',min_el);
-fh{1} = edgefx('geodata',gdat{1},'fs',R,'max_el',max_el,'g',grade,'Ch',1,'Channels',pts2);
-%% Build inner domain
-% To show how the edgefunction filter works, build another edgefunction
-% within the outer domain with a larger minimum edgelength
-% min_el = 150;
-bbox{2} = [-82.8, -82.4; 27.25 28.25];
-gdat{2} = geodata('shp',shoreline,'dem',dem,'bbox',bbox{2},'h0',min_el);
-fh{2} = edgefx('geodata',gdat{2},'fs',R,'max_el',max_el,'g',grade);
-% plotedgefx(fh,1);
-% bound = bbox2bound(bbox{2});
-% plot(bound(:,1),bound(:,2),'r-')
-%% Generate the mesh using the unfiltered edgefunctions
-m_unfiltered = meshgen('ef',fh,'bou',gdat,'plot_on',0);
-m_unfiltered = m_unfiltered.build;
-m_unfiltered = m_unfiltered.grd;
-%% Use the edgefunction filter
-fh = edgefx_filter(fh,gdat);
-% plotedgefx(fh,2)
-%% Generate the mesh using the filtered edgefunction
-m_filtered = meshgen('ef',fh,'bou',gdat,'plot_on',0);
-m_filtered = m_filtered.build;
-m_filtered = m_filtered.grd;
-%% Plot the two meshes
-plot(m_unfiltered,'tri',0)
-title('Unfiltered Edgefunction')
-plot(m_filtered,'tri',0)
-title('Filtered Edgefunction')
\ No newline at end of file
diff --git a/utilities/edgefx_filter.m b/utilities/edgefx_filter.m
deleted file mode 100644
index 44fdfbd0..00000000
--- a/utilities/edgefx_filter.m
+++ /dev/null
@@ -1,137 +0,0 @@
-function fh = edgefx_filter(fh,gdat)
-% This function sorts through layered edgefunctions and puts the minimum
-% target edgelength in the topmost edgefunction. This ensures that when
-% the mesh is generated, the minimum target edgelength from all boxes is
-% used.
-% Coleman Blakely
-% January 22, 2020
-% Edited September 4, 2020 to remove the use of shapevectors and instead
-% use inpoly.m
-if length(fh) == 1
-    disp('No need to do this silly')
-end
-% Loop through edgefunctions to 1) check if they overlap and 2) if they
-% overlap find the points that lie within the upper edgefunction and
-% compare them to the ones below
-for kk = length(fh):-1:2
-    fh_test = fh{kk};
-    % make the points to use with inpoly to find the overlapping points
-    [test_pts, ~] = edgefx2pts(fh_test);
-    test_vals = fh_test.F.Values;
-    x_test = fh_test.F.GridVectors{1};
-    y_test = fh_test.F.GridVectors{2};
-    [X_test, Y_test] = ndgrid(x_test,y_test);
-    % Loop through all boxes that may be below the test box
-    for ii = kk-1:-1:1
-        fh_lower = fh{ii};
-        % find the test points that are within the lower boundary
-        id_overlap = inpoly(test_pts,fh_lower.boubox(1:end-1,:));
-        % check to see if the test box overlaps the lower box. If not, skip
-        % this box
-        if sum(id_overlap)==0
-           continue 
-        end
-        overlap_pts = test_pts(id_overlap,:);
-        % create the grid of overlapping points to use in the gridded
-        % interpolant
-        x = unique(overlap_pts(:,1));
-        y = unique(overlap_pts(:,2));
-        [X, Y] = ndgrid(x,y);
-        % evaluate the gridded interpolant of the lower box using the
-        % points that are within the test box
-        lower_vals = fh_lower.F(X,Y);
-        % find the test points within the test box and compare with the
-        % newly found points
-        [~,idx,~] = intersect(x_test,x);
-        [~,idy,~] = intersect(y_test,y);
-        % int_vals is a dummy variable to make the indexing of the boxes
-        % the same for comparison
-        int_vals = test_vals(idx,idy);
-        % replace all test values that are greater than the lower values
-        int_vals(int_vals>lower_vals) = lower_vals(int_vals>lower_vals);
-        % place the int_vals back into the test values
-        test_vals(idx,idy) = int_vals;
-    end
-    % put the new points back in the edgefunction
-    fh{kk}.F = griddedInterpolant(X_test,Y_test,test_vals);
-    
-    
-    
-end
-% now resmooth the edgefunctions using the defined grade to eliminate
-% any discontinuities, etc.
-% most of this was pulled from stuff that William and Keith did. I do
-% not know exactly how all of it works but it appears to do the trick.
-% CPB
-for kk = length(fh):-1:1
-    obj = fh{kk};
-    feat = gdat{kk};
-    hh_m = obj.F.Values;
-    [xg, yg] = CreateStructGrid(obj);
-    if all(abs(obj.bbox(1,:)) == 180)
-        % for global mesh make it cyclical
-        hh_m = [hh_m(end,:); hh_m; hh_m(1,:)];
-    end
-    hfun = reshape(hh_m',[numel(hh_m),1]);
-    
-    dx = obj.h0*cosd(min(yg(1,:),85)); % for gradient function
-    dy = obj.h0;                       % for gradient function
-    
-    % make g a function of space
-    dmy     = xg*0 ;
-    for param = obj.g'
-        if numel(param)==1 && param~=0
-            lim   = obj.g(1);
-            dmy  = dmy + lim ;
-        else
-            lim  = param(1);
-            dp1  = param(2);
-            dp2  = param(3);
-            
-            limidx = (feat.Fb(xg,yg) < dp1 & ...
-                feat.Fb(xg,yg) > dp2) ;
-            
-            dmy( limidx ) = lim;
-        end
-    end
-    if all(abs(obj.bbox(1,:)) == 180)
-        % for global mesh make it cyclical
-        dmy = [dmy(end,:); dmy; dmy(1,:)];
-    end
-    fdfdx = reshape(dmy',[numel(dmy),1]);
-    clearvars dmy;
-    [hfun,flag] = limgradStruct(obj.ny,dx,dy,hfun,...
-        fdfdx,sqrt(length(hfun)));
-    if flag == 1
-        disp('Gradient relaxing converged!');
-    else
-        error(['FATAL: Gradient relaxing did not converge, '
-            'please check your edge functions']);
-    end
-    % reshape it back
-    hh_m = reshape(hfun,obj.ny,[])';
-    clearvars hfun fdfdx
-    if all(abs(obj.bbox(1,:)) == 180)
-        % for global mesh make it cyclical
-        hh_m = hh_m(2:end-1,:);
-        hh_m(end,:) = hh_m(1,:);
-    end
-    obj.F = griddedInterpolant(xg,yg,hh_m,'linear','nearest');
-    
-    clearvars xg yg
-    fh{kk} = obj;
-end
-end
-
-%% This is a quick little function to turn an edgefuncton into points and values that can be used with fastscatter
-function [pts, vals] = edgefx2pts(fh)
-% takes an edgefunction from OceanMesh2D and outputs a series of points and
-% minimum element sizes for each point
-x = fh.F.GridVectors{1};
-y = fh.F.GridVectors{2}';
-val = fh.F.Values;
-vals = reshape(val,[],1);
-[lat, lon] = meshgrid(x,y);
-c = cat(2,lat',lon');
-pts = reshape(c,[],2);
-end
\ No newline at end of file

From f15e97b1347fb5c34fc58f6ff584aa0861bddbbb Mon Sep 17 00:00:00 2001
From: keith roberts <krober@usp.br>
Date: Sat, 5 Sep 2020 15:29:19 -0300
Subject: [PATCH 4/9] typo

---
 @meshgen/meshgen.m | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/@meshgen/meshgen.m b/@meshgen/meshgen.m
index 81c5eb5c..1dbbcdf6 100644
--- a/@meshgen/meshgen.m
+++ b/@meshgen/meshgen.m
@@ -292,7 +292,7 @@
                         if length(obj.ef) > 1
                             % kjr 2020, ensure the min. sizing func is
                             % used
-                            obj.ef = enforce_min(obj.ef); 
+                            obj.ef = enforce_min_ef(obj.ef); 
                             obj.ef = smooth_outer(obj.ef,obj.Fb);
                         end
                         

From 572726aae513220b337910926026c8503c0bdc80 Mon Sep 17 00:00:00 2001
From: keith roberts <krober@usp.br>
Date: Sat, 19 Sep 2020 18:30:45 -0300
Subject: [PATCH 5/9] Adding Slack(#117)

* Update README.md
---
 README.md | 10 ++++++++--
 1 file changed, 8 insertions(+), 2 deletions(-)

diff --git a/README.md b/README.md
index e8d291a5..75464cd0 100644
--- a/README.md
+++ b/README.md
@@ -6,6 +6,14 @@ This is the default and recommended `PROJECTION` branch. Please use it unless yo
 
 OceanMesh2D is a set of user-friendly MATLAB functions to generate two-dimensional (2D) unstructured meshes for coastal ocean circulation problems. These meshes are based on a variety of feature driven geometric and bathymetric mesh size functions, which are generated according to user-defined parameters. Mesh generation is achieved through a force-balance algorithm combined with a number of topological improvement strategies aimed at improving the worst case triangle quality. The software embeds the mesh generation process into an object-orientated framework that contains pre- and post-processing workflows, which makes mesh generation flexible, reproducible, and script-able. 
 
+## `Slack`
+
+Besides posting 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-hcu2nag7-NUBw52cxxlYupLrc1hqvhw).
+
+Otherwise please reach out to either Dr. William Pringle (wpringle@nd.edu) or Dr. Keith Roberts (krober@usp.br) with questions or concerns or feel free to start an Issue in the issues tab above. 
+
+
+
 ## `Code framework` 
 `OceanMesh2D`  consists of four standalone classes that are called in sequence. It requires no paid toolboxes to build meshes and has been tested to work with a trial version of MATLAB.
 
@@ -39,9 +47,7 @@ Development paper[1]└── Examples/Example_6_GBAY.m %<- An example of the po
 
 All pull requests are tested with Jenkins on a local host. However, to ensure the software is fully functional on your system it's encouraged to run the tests in Tests/ yourself. 
 
-## `Questions?`
 
-Please reach out to either Dr. William Pringle (wpringle@nd.edu) or Dr. Keith Roberts (krober@usp.br) with questions or concerns or feel free to start an Issue in the issues tab above. 
 
 ## `References!`
 

From e6ed6dcc6ca081d515a92e9bc44bed652578eb00 Mon Sep 17 00:00:00 2001
From: William Pringle <wpringle@nd.edu>
Date: Sat, 19 Sep 2020 17:37:48 -0500
Subject: [PATCH 6/9] adding Coleman's example and enforceMin option

- Added modified version of Coleman's example
- In order to show effects of with and without enforce_min_ef I added the enforceMin option which is true by default.
---
 @meshgen/meshgen.m        | 12 ++++++--
 Examples/Example_9_TBAY.m | 65 +++++++++++++++++++++++++++++++++++++++
 2 files changed, 74 insertions(+), 3 deletions(-)
 create mode 100644 Examples/Example_9_TBAY.m

diff --git a/@meshgen/meshgen.m b/@meshgen/meshgen.m
index 1dbbcdf6..97527c22 100644
--- a/@meshgen/meshgen.m
+++ b/@meshgen/meshgen.m
@@ -49,6 +49,7 @@
         annData       % datat contained with KD-tree in anno
         Fb            % bathymetry data interpolant 
         enforceWeirs  % whether or not to enforce weirs in meshgen 
+        enforceMin    % whether or not to enfore minimum edgelength for all edgefxs
     end
     
     
@@ -102,7 +103,7 @@
             addOptional(p,'proj',defval);
             addOptional(p,'qual_tol',defval);
             addOptional(p,'enforceWeirs',1);
-
+            addOptional(p,'enforceMin',1);
             
             % parse the inputs
             parse(p,varargin{:});
@@ -114,7 +115,8 @@
             % 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','fh','inner','outer','mainland',...
+            inp = orderfields(inp,{'h0','bbox','enforceWeirs','enforceMin','fh',...
+                                   'inner','outer','mainland',...
                                    'bou','ef',... %<--OceanMesh classes come after
                                    'egfix','pfix','fixboxes',...
                                    'plot_on','nscreen','itmax',...
@@ -292,7 +294,9 @@
                         if length(obj.ef) > 1
                             % kjr 2020, ensure the min. sizing func is
                             % used
-                            obj.ef = enforce_min_ef(obj.ef); 
+                            if obj.enforceMin
+                                obj.ef = enforce_min_ef(obj.ef); 
+                            end
                             obj.ef = smooth_outer(obj.ef,obj.Fb);
                         end
                         
@@ -384,6 +388,8 @@
                         del = setProj(dmy,1,obj.proj) ;
                     case('enforceWeirs')
                         obj.enforceWeirs = inp.(fields{i}); 
+                    case('enforceMin')
+                        obj.enforceMin = inp.(fields{i}); 
                 end
             end
             
diff --git a/Examples/Example_9_TBAY.m b/Examples/Example_9_TBAY.m
new file mode 100644
index 00000000..c0899cec
--- /dev/null
+++ b/Examples/Example_9_TBAY.m
@@ -0,0 +1,65 @@
+% Example_9_TBay.m (Tampa Bay)
+% 
+% This script serves as an example of how the enforce_min_ef function can
+% ensure that the smallest possible target edgelength function in a series
+% of edgefunctions is carried through to the innermost edgefunction before
+% being used to generate a mesh in meshgen. The reason that this function
+% can sometimes be necessary is because in different locations, different
+% controlling edgefunction parameters can result in larger target
+% edgelengths being put "above" or "inside" the edgefunction that meshgen
+% utilizes. A simple example, and one used here, is if the channels
+% parameter is utilized in an outer edgefunction but not an inner
+% edgefunction. As can be seen by the two meshes generated using this
+% script, the enforceMin option fixes this problem and ensures the channels
+% are still incorporated in the final mesh.
+%
+% Limitation: The min_el (h0) in the inner cannot be larger than in the outer
+%
+% Coleman Blakely
+% September 4, 2020
+% Modifications by William pringle Sep 19, 2020
+
+%% Setup
+clearvars; clc; close all;
+
+%% Define data sources
+% Shoreline
+shoreline = 'GSHHS_f_L1';
+% Bathy data
+dem = 'SRTM15+V2.nc';
+%% Define meshing domain and parameters
+bbox{1} = [-83, -82; 27 28.5];      % outermost bbox
+min_el = 100;                       % minimum edgelength
+max_el = 1e3;                       % maximum edgelength
+grade = 0.35;                       % maximum change in size of edgelength
+R = 3;                              % elements required to resolve a feature
+load ECGC_Thalwegs.mat
+
+%% Build geodata and edgelength function for outermost domain
+gdat{1} = geodata('shp',shoreline,'dem',dem,'bbox',bbox{1},'h0',min_el);
+fh{1}   = edgefx('geodata',gdat{1},'fs',R,'max_el',max_el,...
+                 'g',grade,'Ch',1,'Channels',pts2);
+
+%% Build inner domain
+% To show how the enforce_min_ef function works, build another edgefunction
+% within the outer domain that does not use the channel function
+bbox{2} = [-82.8, -82.4; 27.25 28.25];
+gdat{2} = geodata('shp',shoreline,'dem',dem,'bbox',bbox{2},'h0',min_el);
+fh{2}   = edgefx('geodata',gdat{2},'fs',R,'max_el',max_el,'g',grade);
+
+%% Generate the mesh without enforcing the min of all edgefunctions
+m_nomin = meshgen('ef',fh,'bou',gdat,'enforceMin',0,'plot_on',0);
+m_nomin = m_nomin.build;
+m_nomin = m_nomin.grd;
+
+%% Generate the mesh with enforcing the min of all edgefunctions 
+%-> by default enforceMin = 1 so do not need to set usually
+m_min = meshgen('ef',fh,'bou',gdat,'enforceMin',1,'plot_on',0);
+m_min = m_min.build;
+m_min = m_min.grd;
+
+%% Plot the two meshes
+plot(m_nomin,'tri')
+title('Without Enforcing Minimum of all Edgefunctions')
+plot(m_min,'tri')
+title('Enforcing Minimum of all Edgefunctions')
\ No newline at end of file

From cb5ac5d91aed78d880929eb6f4c27edf6f121103 Mon Sep 17 00:00:00 2001
From: keith roberts <krober@usp.br>
Date: Thu, 15 Oct 2020 21:23:56 -0300
Subject: [PATCH 7/9] turn off 3d shapefile reading by default (#121)

---
 @geodata/geodata.m                | 7 +++++--
 @geodata/private/Read_shapefile.m | 5 ++++-
 2 files changed, 9 insertions(+), 3 deletions(-)

diff --git a/@geodata/geodata.m b/@geodata/geodata.m
index e49fb093..68053b8e 100644
--- a/@geodata/geodata.m
+++ b/@geodata/geodata.m
@@ -57,6 +57,7 @@
         pslg % piecewise liner straight line graph
         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
     end
     
     methods
@@ -75,7 +76,6 @@
             %addOptional(p,'weirs',defval);
             %addOptional(p,'pslg',defval);
             %addOptional(p,'boubox',defval);
-            %addOptional(p,'window',defval);
             
             % Check for m_map dir
             M_MAP_EXISTS=0 ;
@@ -119,6 +119,7 @@
             addOptional(p,'pslg',defval);
             addOptional(p,'boubox',defval);
             addOptional(p,'window',defval);
+            addOptional(p,'shapefile_3d',defval);
             
             % parse the inputs
             parse(p,varargin{:});
@@ -202,6 +203,8 @@
                             % Default value
                             obj.window = 5;
                         end
+                    case('shapefile_3d')
+                         obj.shapefile_3d = 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}))
@@ -321,7 +324,7 @@
                 end
                 
                 polygon_struct = Read_shapefile( obj.contourfile, [], ...
-                    obj.bbox, obj.gridspace, obj.boubox, 0 );
+                    obj.bbox, obj.gridspace, obj.boubox, 0, obj.shapefile_3d);
                 
                 % Unpack data from function Read_Shapefile()s
                 obj.outer     = polygon_struct.outer;
diff --git a/@geodata/private/Read_shapefile.m b/@geodata/private/Read_shapefile.m
index ff6e8c49..94a146c9 100644
--- a/@geodata/private/Read_shapefile.m
+++ b/@geodata/private/Read_shapefile.m
@@ -1,5 +1,5 @@
 function polygon_struct = Read_shapefile( finputname, polygon, bbox, ...
-                                          h0, boubox, plot_on )
+                                          h0, boubox, plot_on, shapefile_3d)
 % Read_shapefile: Reads a shapefile or a NaN-delimited vector
 % containing polygons and/or segments in the the desired region
 % of interest. Classifies the vector data as either a
@@ -174,6 +174,9 @@
         % using m_shaperead
         points = tmpC{i,1}(1:end,:) ;
         if size(points,2) == 3
+            points = points(:,1:2);
+        end
+        if shapefile_3d
             % if 3-D shapefile
             height = points(:,3); 
             points = points(:,1:2); 

From bbb070f67a85b0652bd4acdbaa3b6f334f2e667b Mon Sep 17 00:00:00 2001
From: William Pringle <wpringle@nd.edu>
Date: Sun, 18 Oct 2020 10:11:27 -0500
Subject: [PATCH 8/9] Multiscale meshing--mesh resolution transitions  (#119)

* Bug fix to dpoly.m

   * if sum(inside) == 0 then d_l will not exist, but it seems like it sum(inside) is zero then it doesn't use the d_l
     anyway so can just return then.

* use splitting at initial point generation to ensure mesh resolution transitions between largely disparate nests are sufficiently smooth for good mesh qualities.


Co-authored-by: Keith Roberts <krober10@nd.edu>
---
 @meshgen/meshgen.m                        | 90 +++++++++++------------
 @meshgen/private/dpoly.m                  |  7 +-
 @msh/msh.m                                |  8 +-
 Examples/Example_10_Multiscale_Smoother.m | 29 ++++++++
 utilities/split.m                         | 22 ++++++
 5 files changed, 105 insertions(+), 51 deletions(-)
 create mode 100644 Examples/Example_10_Multiscale_Smoother.m
 create mode 100644 utilities/split.m

diff --git a/@meshgen/meshgen.m b/@meshgen/meshgen.m
index 97527c22..22aea4b4 100644
--- a/@meshgen/meshgen.m
+++ b/@meshgen/meshgen.m
@@ -16,6 +16,30 @@
     %
     %   You should have received a copy of the GNU General Public License
     %   along with this program.  If not, see <http://www.gnu.org/licenses/>.
+    %
+    %   Available options:
+    %         ef            % edgefx class
+    %         bou           % geodata class
+    %         h0            % minimum edge length (optional if bou exists)
+    %         bbox          % bounding box [xmin,ymin; xmax,ymax] (manual specification, no bou)
+    %         proj          % structure containing the m_map projection info
+    %         plot_on       % flag to plot (def: 1) or not (0)
+    %         nscreen       % how many it to plot and write temp files (def: 5)
+    %         itmax         % maximum number of iterations.
+    %         pfix          % fixed node positions (nfix x 2 )
+    %         egfix         % edge constraints
+	%         outer         % meshing boundary (manual specification, no bou)
+    %         inner         % island boundaries (manual specification, no bou)      
+    %         mainland      % the shoreline boundary (manual specification, no bou)
+    %         fixboxes      % a flag that indicates which boxes will use fixed constraints
+    %         memory_gb     % memory in GB allowed to use for initial rejector
+    %         cleanup       % logical flag or string to trigger cleaning of topology (default is on).
+    %         direc_smooth  % logical flag to trigger direct smoothing of mesh in the cleanup
+    %         dj_cutoff     % the cutoff area fraction for disjoint portions to delete
+    %         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
+    %         big_mesh      % set to 1 to remove the bou data from memory
     properties
         fd            % handle to distance function
         fh            % handle to edge function
@@ -30,9 +54,9 @@
         bou           % geodata class
         ef            % edgefx class
         itmax         % maximum number of iterations.
-        outer         % meshing boundary
-        inner         % island boundaries      
-        mainland      % the shoreline boundary 
+        outer         % meshing boundary (manual specification, no bou)
+        inner         % island boundaries (manual specification, no bou)      
+        mainland      % the shoreline boundary (manual specification, no bou)
         boubox        % the bbox as a polygon 2-tuple
         inpoly_flip   % used to flip the inpoly test to determine the signed distance.
         memory_gb     % memory in GB allowed to use for initial rejector
@@ -40,8 +64,7 @@
         direc_smooth  % logical flag to trigger direct smoothing of mesh in the cleanup
         dj_cutoff     % the cutoff area fraction for disjoint portions to delete
         grd = msh();  % create empty mesh class to return p and t in.
-        big_mesh
-        ns_fix        % improve spacing for boundary vertices
+        big_mesh      % release bou data from memory
         qual          % mean, lower 3rd sigma, and the minimum element quality.
         qual_tol      % tolerance for the accepted negligible change in quality
         proj          % structure containing the m_map projection info
@@ -99,7 +122,6 @@
             addOptional(p,'direc_smooth',1);
             addOptional(p,'dj_cutoff',0.25);
             addOptional(p,'big_mesh',defval);
-            addOptional(p,'ns_fix',defval);
             addOptional(p,'proj',defval);
             addOptional(p,'qual_tol',defval);
             addOptional(p,'enforceWeirs',1);
@@ -122,7 +144,7 @@
                                    'plot_on','nscreen','itmax',...
                                    'memory_gb','qual_tol','cleanup',...
                                    'direc_smooth','dj_cutoff',...
-                                   'big_mesh','ns_fix','proj'});             
+                                   'big_mesh','proj'});             
             % get the fieldnames of the edge functions
             fields = fieldnames(inp);
             % loop through and determine which args were passed.
@@ -185,7 +207,7 @@
                         end
                         if obj.enforceWeirs
                             for j = 1 : length(obj.bou)
-                           	if ~isempty(obj.bou{j}.weirEgfix) && ~isempty(obj.egfix)
+                                if ~isempty(obj.bou{j}.weirEgfix) && ~isempty(obj.egfix)
                                     obj.egfix = [obj.egfix ; obj.bou{j}.weirEgfix+max(obj.egfix(:))];
                                 else
                                     obj.egfix =  obj.bou{j}.weirEgfix;
@@ -195,7 +217,6 @@
                         obj.egfix = renumberEdges(obj.egfix);
                     case('fixboxes')
                         obj.fixboxes= inp.(fields{i});
-                    
                     case('bou')
                         % got it from user arg
                         if obj.outer~=0, continue; end
@@ -311,8 +332,6 @@
                         obj.plot_on= inp.(fields{i});
                     case('big_mesh')
                         obj.big_mesh = inp.(fields{i});
-                    case('ns_fix')
-                        obj.ns_fix   = inp.(fields{i});
                     case('nscreen')
                         obj.nscreen= inp.(fields{i});
                         if obj.nscreen ~=0
@@ -455,13 +474,12 @@
             %%
             tic
             it = 1 ;
-            imp = 10;
-            qual_diff = 0;
             Re = 6378.137e3;
             geps = 1e-3*min(obj.h0)/Re; 
             deps = sqrt(eps);
             ttol=0.1; Fscale = 1.2; deltat = 0.1;
             delIT = 0 ; delImp = 2;
+            imp = 10; % number of iterations to do mesh improvements (delete/add)
 
             % unpack initial points.
             p = obj.grd.p;
@@ -546,14 +564,24 @@
                         p1 = p1(rand(size(p1,1),1) < r0/max_r0,:);          % Rejection method
                         p = [p; p1];                                        % Adding p1 to p
                     end
+                    % add new points along boundary of multiscale nests
+                    if box_num < length(obj.h0)
+                        h0_rat = ceil(h0_l/obj.h0(box_num+1));
+                        nsplits = floor(log(h0_rat)/log(2));
+                        for add = 1:nsplits
+                            new_points = split(p,fh_l);
+                            p = [p; new_points];
+                        end
+                    end
                 end
             else
                 disp('User-supplied initial points!');
                 obj.grd.b = [];
-                imp = 10; % number of iterations to do mesh improvements (delete/add)
                 h0_l = obj.h0(end); % finest h0 (in case of a restart of meshgen.build).
             end
             
+     
+            
             % remove pfix/egfix outside of desired subdomain
             nfix = size(obj.pfix,1);    % Number of fixed points
             negfix = size(obj.egfix,1); % Number of edge constraints
@@ -687,8 +715,8 @@
                  end
                 
                 % Termination quality, mesh quality reached is copacetic.
+                qual_diff = mq_l3sig - obj.qual(max(1,it-imp),2);
                 if ~mod(it,imp)
-                    qual_diff = mq_l3sig - obj.qual(max(1,it-imp),2);
                     if abs(qual_diff) < obj.qual_tol
                         % Do the final elimination of small connectivity
                         [t,p] = delaunay_elim(p,obj.fd,geps,1);
@@ -752,7 +780,7 @@
                 % Mesh improvements (deleting and addition)
                 if ~mod(it,imp)
                     nn = []; pst = [];
-                    if qual_diff < imp*0.01 && qual_diff > 0
+                    if qual_diff < imp*obj.qual_tol && qual_diff > 0
                         % Remove elements with small connectivity
                         nn = get_small_connectivity(p,t);
                         disp(['Deleting ' num2str(length(nn)) ' due to small connectivity'])
@@ -953,35 +981,7 @@
             
             
         end % end distmesh2d_plus
-        
-        function obj = nearshorefix(obj)
-            %% kjr make sure boundaries have good spacing on boundary.
-            % This is experimentary. 
-            t = obj.grd.t ; p = obj.grd.t;
-            [bnde, ~] = extdom_edges2(t,p);
-            [poly]  = extdom_polygon(bnde,p,1);
-
-            new = [];
-            for j = 1 : length(poly)
-                for i = 1 : length(poly{j})-2
-                    pt = poly{j}(i,:) ; % current point
-                    nxt= poly{j}(i+1,:) ; % next point
-                    nxt2 = poly{j}(i+2,:) ; % next next point
-
-                    dst1 = sqrt( (nxt(:,1)-pt(:,1)).^2 + (nxt(:,2)-pt(:,2)).^2 );     % dist to next point
-                    dst2 = sqrt( (nxt2(:,1)-nxt(:,1)).^2 + (nxt2(:,2)-nxt(:,2)).^2 ); % dist to next next point
-
-                    if dst2/dst1 > 2
-                        % split bar
-                        q = (nxt2 + nxt)/2;
-                        new = [new; q];
-                    end
-                end
-            end
-            p = [p; new]; % post fix new points (to avoid problems with pfix.)
-            t = delaunay_elim(p,obj.fd,geps,0);       % Delaunay with elimination
-            obj.grd.t = t ; obj.grd.p = t;
-        end
+       
         
         
     end % end methods
diff --git a/@meshgen/private/dpoly.m b/@meshgen/private/dpoly.m
index bb8cb138..f376eb3d 100644
--- a/@meshgen/private/dpoly.m
+++ b/@meshgen/private/dpoly.m
@@ -96,15 +96,12 @@
             in_outer = ~in_outer;
         end
     else
-        in_boubox   = d_l*0;
-        in_outer    = d_l*0;
+        return
     end
     
     % d is signed negative if inside and vice versa.
     d_l = (-1).^( in_outer & in_boubox).*d_l;
-    
-    if sum(inside)==0; return; end
-    
+   
     d(inside) = d_l;
 end
 % EOF
diff --git a/@msh/msh.m b/@msh/msh.m
index f91b5230..618f25a5 100644
--- a/@msh/msh.m
+++ b/@msh/msh.m
@@ -279,7 +279,9 @@ function write(obj,fname,type)
             %    i) 'numticks': number of colorbar tickmarks and (optional) range:
             %           [numticks] or [numticks caxis_lower caxis_upper]
             %    ii) 'fontsize': figure fontsize
-            %   iii) 'holdon'  : =1 to plot on existing figure (otherwise
+            %   iii) 'backcolor': figure background RGB color (where mesh
+            %                  doesn't exist), default is [1 1 1] => white
+            %   iv)  'holdon'  : =1 to plot on existing figure (otherwise
             %                    will use new figure)
             if nargin < 2
                 type = 'tri';
@@ -292,11 +294,14 @@ function write(obj,fname,type)
             end
             np_g = length(obj.p) ;
             fsz = 12; % default font size
+            bgc = [1 1 1]; % default background color
             numticks = 10; % default num of ticks
             holdon = 0;
             for kk = 1:2:length(varargin)
                 if strcmp(varargin{kk},'fontsize')
                     fsz = varargin{kk+1};
+                elseif strcmp(varargin{kk},'backcolor')
+                    bgc = varargin{kk+1};
                 elseif strcmp(varargin{kk},'numticks')
                     numticks = varargin{kk+1};
                 elseif strcmp(varargin{kk},'holdon')
@@ -790,6 +795,7 @@ function write(obj,fname,type)
             end
             ax = gca;
             ax.FontSize = fsz;
+            ax.Color = bgc;
             if proj == 1
                 % now add the box
                 m_grid('FontSize',fsz);
diff --git a/Examples/Example_10_Multiscale_Smoother.m b/Examples/Example_10_Multiscale_Smoother.m
new file mode 100644
index 00000000..5895ee7f
--- /dev/null
+++ b/Examples/Example_10_Multiscale_Smoother.m
@@ -0,0 +1,29 @@
+% Example_10_Multiscale_Smoother: 
+% An idealized test for multiscale nesting using boxes with a large min_el
+% ratio
+
+bbox = [0, 1; 0,1]; 
+boubox = [0,0; 1,0; 1,1; 0,1; 0,0; NaN NaN ];
+min_el = 1e3; 
+
+gdat1 = geodata('pslg',boubox,'bbox',bbox,'h0',min_el);
+fh1 = edgefx('geodata',gdat1,'g',0.2);
+
+bbox2 = [-1, 2; -1,2]; 
+boubox2 = [-1,-1; 2,-1; 2,2; -1,2; -1,-1; NaN NaN ];
+min_el2 = min_el*10; 
+
+gdat2 = geodata('pslg',boubox2,'bbox',bbox2,'h0',min_el2);
+fh2 = edgefx('geodata',gdat2,'g',0.2);
+
+mshopts = meshgen('ef',{fh2, fh1},'bou',{gdat2,gdat1},...
+                  'plot_on',1,'qual_tol',0.0025,'cleanup',0);
+mshopts = mshopts.build;
+
+m = mshopts.grd; 
+
+plot(m)
+
+m = m.clean;
+
+plot(m)
\ No newline at end of file
diff --git a/utilities/split.m b/utilities/split.m
new file mode 100644
index 00000000..3e0db6fb
--- /dev/null
+++ b/utilities/split.m
@@ -0,0 +1,22 @@
+% Function for adding points to an untriangulated set of points based on
+% nearest distance versus the desired edgefx 
+function [new_points]=split(points,fh)
+    % find the point closest to each point (that isn't the *same* point)
+    [idx, ~] = ourKNNsearch(points',points',2);
+    % the ideal spacing between points
+    ideal_dst = fh(points);
+    % where the dst is 2*ideal_dist, add a point in between
+    long   = zeros(length(points)*2,1);
+    lat    = zeros(length(points)*2,1);
+    long(1:2:end) = points(idx(:,1),1);
+    long(2:2:end) = points(idx(:,2),1);
+    lat(1:2:end)  = points(idx(:,1),2);
+    lat(2:2:end)  = points(idx(:,2),2);
+    L = m_lldist(long,lat);
+    L = L(1:2:end)*1e3;            % L = lengths in meters
+    splits = L > 1.5*ideal_dst;
+    disp(['Number of splits near multiscale nest: ' num2str(sum(splits))])
+    new_points = (points(idx(splits,1),:) + points(idx(splits,2),:))./2.0;
+end
+
+

From 2f757b1a6c532d9aea42e5f5f50821a06dffc2dc Mon Sep 17 00:00:00 2001
From: keith roberts <krober@usp.br>
Date: Tue, 20 Oct 2020 14:52:12 -0300
Subject: [PATCH 9/9] Update README.md (#124)

---
 README.md | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/README.md b/README.md
index 75464cd0..820f0e83 100644
--- a/README.md
+++ b/README.md
@@ -8,7 +8,7 @@ OceanMesh2D is a set of user-friendly MATLAB functions to generate two-dimension
 
 ## `Slack`
 
-Besides posting 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-hcu2nag7-NUBw52cxxlYupLrc1hqvhw).
+Besides posting 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-igvr8rvo-OctG9t0QEYNyyvVr55ynUQ).
 
 Otherwise please reach out to either Dr. William Pringle (wpringle@nd.edu) or Dr. Keith Roberts (krober@usp.br) with questions or concerns or feel free to start an Issue in the issues tab above.