Skip to content

Commit

Permalink
Multiscale meshing--mesh resolution transitions (#119)
Browse files Browse the repository at this point in the history 
* 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]>
  • Loading branch information
@WPringle
WPringle and Keith Roberts authored Oct 18, 2020
1 parent cb5ac5d commit bbb070f
Show file tree
Hide file tree
Showing 5 changed files with 105 additions and 51 deletions.
90 changes: 45 additions & 45 deletions @meshgen/meshgen.m
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand All @@ -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
Expand Down Expand Up @@ -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);
Expand All @@ -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.
Expand Down Expand Up @@ -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;
Expand All @@ -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
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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;
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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);
Expand Down Expand Up @@ -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'])
Expand Down Expand Up @@ -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
Expand Down
7 changes: 2 additions & 5 deletions @meshgen/private/dpoly.m
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
8 changes: 7 additions & 1 deletion @msh/msh.m
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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';
Expand All @@ -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')
Expand Down Expand Up @@ -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);
Expand Down
29 changes: 29 additions & 0 deletions Examples/Example_10_Multiscale_Smoother.m
View file
Original file line number Diff line number Diff line change
@@ -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)
22 changes: 22 additions & 0 deletions utilities/split.m
View file
Original file line number Diff line number Diff line change
@@ -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


0 comments on commit bbb070f

Please sign in to comment.