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 <[email protected]>

Author: WPringle

@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,18 +54,17 @@
         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
         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
         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

@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

@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);

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)

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
+
+