diff --git a/@geodata/geodata.m b/@geodata/geodata.m index c6ebd5ca..c20bb9d3 100644 --- a/@geodata/geodata.m +++ b/@geodata/geodata.m @@ -17,6 +17,7 @@ % addOptional(p,'pslg',defval); % addOptional(p,'boubox',defval); % addOptional(p,'window',defval); + % addOptional(p,'high_fidelity',defval); % % This program is free software: you can redistribute it and/or modify % it under the terms of the GNU General Public License as published by @@ -58,6 +59,8 @@ 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 + high_fidelity % performs a 1D mesh generation step to form pfix and egfix prior to 2D meshing + end methods @@ -120,6 +123,8 @@ addOptional(p,'boubox',defval); addOptional(p,'window',defval); addOptional(p,'shapefile_3d',defval); + addOptional(p,'high_fidelity',defval); + % parse the inputs parse(p,varargin{:}); @@ -204,7 +209,9 @@ obj.window = 5; end case('shapefile_3d') - obj.shapefile_3d = inp.(fields{i}) ; + obj.shapefile_3d = inp.(fields{i}) ; + case('high_fidelity') + obj.high_fidelity = 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})) @@ -861,4 +868,3 @@ function plot(obj,type,projection,holdon) end end - diff --git a/@geodata/private/Read_shapefile.m b/@geodata/private/Read_shapefile.m index ca807b25..ac748107 100644 --- a/@geodata/private/Read_shapefile.m +++ b/@geodata/private/Read_shapefile.m @@ -287,4 +287,4 @@ end end %EOF -end +end \ No newline at end of file diff --git a/@meshgen/meshgen.m b/@meshgen/meshgen.m index 5d4efee7..3c3a35dc 100644 --- a/@meshgen/meshgen.m +++ b/@meshgen/meshgen.m @@ -1,1027 +1,1173 @@ -classdef meshgen - % MESHGEN: Mesh generation class - % Handles input parameters to create a meshgen class object that can be - % used to build a msh class. - % Copyright (C) 2018 Keith Roberts & William Pringle - % - % This program is free software: you can redistribute it and/or modify - % it under the terms of the GNU General Public License as published by - % the Free Software Foundation, either version 3 of the License, or - % (at your option) any later version. - % - % This program is distributed in the hope that it will be useful, - % but WITHOUT ANY WARRANTY; without even the implied warranty of - % MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the - % GNU General Public License for more details. - % - % You should have received a copy of the GNU General Public License - % along with this program. If not, see . - % - % 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 - % delaunay_elim_on_exit % whether or not to run delaunay_elim on exit of meshgen - % improve_with_reduced_quality % whether or not to allow mesh improvements with decreases in mesh quality - % big_mesh % set to 1 to remove the bou data from memory - % - properties - fd % handle to distance function - fh % handle to edge function - h0 % minimum edge length - edgefx % edgefx class - bbox % bounding box [xmin,ymin; xmax,ymax] - pfix % fixed node positions (nfix x 2 ) - egfix % edge constraints - fixboxes % a flag that indicates which boxes will use fixed constraints - plot_on % flag to plot (def: 1) or not (0) - nscreen % how many it to plot and write temp files (def: 5) - bou % geodata class - ef % edgefx class - itmax % maximum number of iterations. - 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 % 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 - anno % Approx. Nearest Neighbor search object. - 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 - delaunay_elim_on_exit % whether or not to run delaunay_elim on exit of meshgen - improve_with_reduced_quality % whether or not to allow mesh improvements with decreases in mesh quality - end - - - methods - - % class constructor/default grd generation options - function obj = meshgen(varargin) - % Check for m_map dir - M_MAP_EXISTS=0; - if exist('m_proj','file')==2 - M_MAP_EXISTS=1 ; - end - if M_MAP_EXISTS~=1 - error('Where''s m_map? Chief, you need to read the user guide') - end - - % Check for utilties dir - UTIL_DIR_EXISTS=0 ; - if exist('inpoly.m','file') - UTIL_DIR_EXISTS=1; - end - if UTIL_DIR_EXISTS~=1 - error('Where''s the utilities directory? Chief, you need to read the user guide') - end - - p = inputParser; - % unpack options and set default ones, catch errors. - - defval = 0; % placeholder value if arg is not passed. - % add name/value pairs - addOptional(p,'h0',defval); - addOptional(p,'bbox',defval); - addOptional(p,'fh',defval); - addOptional(p,'pfix',defval); - addOptional(p,'egfix',defval); - addOptional(p,'fixboxes',defval); - addOptional(p,'inner',defval); - addOptional(p,'outer',defval); - addOptional(p,'mainland',defval); - addOptional(p,'bou',defval); - addOptional(p,'ef',defval); - addOptional(p,'plot_on',defval); - addOptional(p,'nscreen',defval); - addOptional(p,'itmax',defval); - addOptional(p,'memory_gb',1); - addOptional(p,'cleanup',1); - addOptional(p,'direc_smooth',1); - addOptional(p,'dj_cutoff',0.25); - addOptional(p,'big_mesh',defval); - addOptional(p,'proj',defval); - addOptional(p,'qual_tol',defval); - addOptional(p,'enforceWeirs',1); - addOptional(p,'enforceMin',1); - addOptional(p,'delaunay_elim_on_exit',1); - addOptional(p,'improve_with_reduced_quality',0); - - % parse the inputs - parse(p,varargin{:}); - - %if isempty(varargin); return; end - % store the inputs as a struct - inp = p.Results; - - % 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','enforceMin',... - 'delaunay_elim_on_exit','improve_with_reduced_quality',... - 'fh',... - 'inner','outer','mainland',... - 'bou','ef',... %<--OceanMesh classes come after - 'egfix','pfix','fixboxes',... - 'plot_on','nscreen','itmax',... - 'memory_gb','qual_tol','cleanup',... - 'direc_smooth','dj_cutoff',... - 'big_mesh','proj'}); - % get the fieldnames of the edge functions - fields = fieldnames(inp); - % loop through and determine which args were passed. - % also, assign reasonable default values if some options were - % not assigned. - for i = 1 : numel(fields) - type = fields{i}; - switch type - % parse aux options first - case('h0') - % Provide in meters - obj.h0 = inp.(fields{i}); - case('fh') - if isa(inp.(fields{i}),'function_handle') - obj.fh = inp.(fields{i}); - end - % can't check for errors here yet. - case('bbox') - obj.bbox= inp.(fields{i}); - if iscell(obj.bbox) - % checking bbox extents - ob_min = obj.bbox{1}(:,1); - ob_max = obj.bbox{1}(:,2); - for ii = 2:length(obj.bbox) - if any(obj.bbox{ii}(:,1) < ob_min) || ... - any(obj.bbox{ii}(:,2) > ob_max) - error(['Outer bbox must contain all ' ... - 'inner bboxes: inner box #' ... - num2str(ii) ' violates this']) - end - end - end - - % if user didn't pass anything explicitly for - % bounding box make it empty so it can be populated - % from ef as a cell-array - if obj.bbox(1)==0 - obj.bbox = []; - end - case('pfix') - obj.pfix= inp.(fields{i}); - if obj.pfix(1)~=0 - obj.pfix(:,:) = inp.(fields{i}); - else - obj.pfix = []; - end - if obj.enforceWeirs - for j = 1 : length(obj.bou) - if ~isempty(obj.bou{j}.weirPfix) - obj.pfix = [obj.pfix ; obj.bou{j}.weirPfix]; - end - end - end - case('egfix') - obj.egfix= inp.(fields{i}); - if ~isempty(obj.egfix) && obj.egfix(1)~=0 - obj.egfix = inp.(fields{i}); - else - obj.egfix = []; - end - if obj.enforceWeirs - for j = 1 : length(obj.bou) - 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; - end - end - end - 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 - - obj.outer = {} ; - obj.inner = {} ; - obj.mainland = {} ; - - obj.bou = inp.(fields{i}); - - % handle when not a cell - if ~iscell(obj.bou) - boutemp = obj.bou; - obj.bou = cell(1); - obj.bou{1} = boutemp; - end - - % then the geodata class was provide, unpack - for ee = 1:length(obj.bou) - try - arg = obj.bou{ee} ; - catch - arg = obj.bou; - end - if isa(arg,'geodata') - obj.outer{ee} = obj.bou{ee}.outer; - obj.inner{ee} = obj.bou{ee}.inner; - - % save bathy interpolant to meshgen - if ~isempty(obj.bou{ee}.Fb) - obj.Fb{ee} = obj.bou{ee}.Fb ; - end - - if ~isempty(obj.inner{ee}) && ... - obj.inner{ee}(1)~= 0 - obj.outer{ee} = [obj.outer{ee}; - obj.inner{ee}]; - end - obj.mainland{ee} = obj.bou{ee}.mainland; - obj.boubox{ee} = obj.bou{ee}.boubox; - obj.inpoly_flip{ee} = obj.bou{ee}.inpoly_flip; - if obj.big_mesh - % release gdat's - obj.bou{ee}.mainland= []; - obj.bou{ee}.outer= []; - if ~isempty(obj.bou{ee}.inner) - obj.bou{ee}.inner= []; - end - end - end - end - - case('ef') - tmp = inp.(fields{i}); - if isa(tmp, 'function_handle') - error('Please specify your edge function handle through the name/value pair fh'); - end - obj.ef = tmp; - - % handle when not a cell - if ~iscell(obj.ef) - eftemp = obj.ef; - obj.ef = cell(1); - obj.ef{1} = eftemp; - end - - % Gather boxes from ef class. - for ee = 1 : length(obj.ef) - if isa(obj.ef{ee},'edgefx') - obj.bbox{ee} = obj.ef{ee}.bbox; - end - end - - % checking bbox extents - if iscell(obj.bbox) - ob_min = obj.bbox{1}(:,1); - ob_max = obj.bbox{1}(:,2); - for ii = 2:length(obj.bbox) - if any(obj.bbox{ii}(:,1) < ob_min) || ... - any(obj.bbox{ii}(:,2) > ob_max) - error(['Outer bbox must contain all ' ... - 'inner bboxes: inner box #' ... - num2str(ii) ' violates this']) - end - end - end - - % kjr 2018 June: get h0 from edge functions - for ee = 1:length(obj.ef) - if isa(obj.ef{ee},'edgefx') - obj.h0(ee) = obj.ef{ee}.h0; - end - end - - % kjr 2018 smooth the outer automatically - if length(obj.ef) > 1 - % kjr 2020, ensure the min. sizing func is - % used - if obj.enforceMin - obj.ef = enforce_min_ef(obj.ef); - end - obj.ef = smooth_outer(obj.ef,obj.Fb); - end - - % Save the ef interpolants into the edgefx - for ee = 1:length(obj.ef) - if isa(obj.ef{ee},'edgefx') - obj.fh{ee} = @(p)obj.ef{ee}.F(p); - end - end - - case('plot_on') - obj.plot_on= inp.(fields{i}); - case('big_mesh') - obj.big_mesh = inp.(fields{i}); - case('nscreen') - obj.nscreen= inp.(fields{i}); - if obj.nscreen ~=0 - obj.nscreen = inp.(fields{i}); - obj.plot_on = 1; - else - obj.nscreen = 5; % default - end - case('itmax') - obj.itmax= inp.(fields{i}); - if obj.itmax ~=0 - obj.itmax = inp.(fields{i}); - else - obj.itmax = 100; - warning('No itmax specified, itmax set to 100'); - end - case('qual_tol') - obj.qual_tol = inp.(fields{i}); - if obj.qual_tol ~=0 - obj.qual_tol = inp.(fields{i}); - else - obj.qual_tol = 0.01; - end - case('inner') - if ~isa(obj.bou,'geodata') - obj.inner = inp.(fields{i}); - end - case('outer') - if ~isa(obj.bou,'geodata') - obj.outer = inp.(fields{i}); - if obj.inner(1)~=0 - obj.outer = [obj.outer; obj.inner]; - end - end - case('mainland') - if ~isa(obj.bou,'geodata') - obj.mainland = inp.(fields{i}); - end - case('memory_gb') - if ~isa(obj.bou,'memory_gb') - obj.memory_gb = inp.(fields{i}); - end - case('cleanup') - obj.cleanup = inp.(fields{i}); - if isempty(obj.cleanup) || obj.cleanup == 0 - obj.cleanup = 'none'; - elseif obj.cleanup == 1 - obj.cleanup = 'default'; - end - case('dj_cutoff') - obj.dj_cutoff = inp.(fields{i}); - case('direc_smooth') - obj.direc_smooth = inp.(fields{i}); - case('proj') - obj.proj = inp.(fields{i}); - % default CPP - if obj.proj == 0; obj.proj = 'equi'; end - if ~isempty(obj.bbox) - % kjr Oct 2018 use outer coarsest box for - % multiscale meshing - lon_mi = obj.bbox{1}(1,1)-obj.h0(1)/1110; - lon_ma = obj.bbox{1}(1,2)+obj.h0(1)/1110; - lat_mi = obj.bbox{1}(2,1)-obj.h0(1)/1110; - lat_ma = obj.bbox{1}(2,2)+obj.h0(1)/1110; - else - lon_mi = -180; lon_ma = 180; - lat_mi = -90; lat_ma = 90; - end - % Set up projected space - dmy = msh() ; - dmy.p(:,1) = [lon_mi; lon_ma]; - dmy.p(:,2) = [lat_mi; lat_ma]; - del = setProj(dmy,1,obj.proj) ; - case('enforceWeirs') - obj.enforceWeirs = inp.(fields{i}); - case('enforceMin') - obj.enforceMin = inp.(fields{i}); - case('delaunay_elim_on_exit') - obj.delaunay_elim_on_exit = inp.(fields{i}); - case('improve_with_reduced_quality') - obj.improve_with_reduced_quality = inp.(fields{i}); - end - end - - if isempty(varargin); return; end - - % error checking - if isempty(obj.boubox) && ~isempty(obj.bbox) - % Make the bounding box 5 x 2 matrix in clockwise order if - % it isn't present. This case must be when the user is - % manually specifying the PSLG. - obj.boubox{1} = [obj.bbox(1,1) obj.bbox(2,1); - obj.bbox(1,1) obj.bbox(2,2); ... - obj.bbox(1,2) obj.bbox(2,2); - 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 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). - % kjr build ANN object into meshgen - obj = createANN(obj) ; - - global MAP_PROJECTION MAP_COORDS MAP_VAR_LIST - obj.grd.proj = MAP_PROJECTION ; - obj.grd.coord = MAP_COORDS ; - obj.grd.mapvar = MAP_VAR_LIST ; - - end - - % Creates Approximate nearest neighbor objects on start-up - function obj = createANN(obj) - - box_vec = 1:length(obj.bbox); - - for box_num = box_vec - if ~iscell(obj.outer) - dataset = obj.outer; - dataset(isnan(obj.outer(:,1)),:) = []; - else - dataset = obj.outer{box_num}; - dataset(isnan(obj.outer{box_num}(:,1)),:) = []; - end - if all(abs(obj.bbox{box_num}(1,:)) == 180) - % This line removes the line that can appear in the - % center for a global mesh - dataset(abs(dataset(:,1)) > 180-1e-6,:) = []; - dataset(abs(dataset(:,1)) < 1e-6,:) = []; - end - [dataset(:,1),dataset(:,2)] = m_ll2xy(dataset(:,1),dataset(:,2)); - dataset(isnan(dataset(:,1)),:) = []; - dmy = ann(dataset'); - obj.anno{box_num} = dmy; - obj.annData{box_num}=dataset; - end - end - - function obj = build(obj) - %DISTMESH2D 2-D Mesh Generator using Distance Functions. - % Checking existence of major inputs - %% - warning('off','all') - %% - tic - it = 1 ; - 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; - if isempty(p) - disp('Forming initial point distribution...'); - % loop over number of boxes - for box_num = 1:length(obj.h0) - disp([' for box #' num2str(box_num)]); - % checking if cell or not and applying local values - h0_l = obj.h0(box_num); - max_r0 = 1/h0_l^2; - if ~iscell(obj.bbox) - bbox_l = obj.bbox'; % <--we must tranpose this! - else - bbox_l = obj.bbox{box_num}'; % <--tranpose! - end - if ~iscell(obj.fh) - fh_l = obj.fh; - else - fh_l = obj.fh{box_num}; - end - % Lets estimate the num_points the distribution will be - num_points = ceil(2/sqrt(3)*prod(abs(diff(bbox_l)))... - /(h0_l/111e3)^2); - noblks = ceil(num_points*2*8/obj.memory_gb*1e-9); - len = abs(bbox_l(1,1)-bbox_l(2,1)); - blklen = len/noblks; - st = bbox_l(1,1) ; ed = st + blklen; ns = 1; - %% 1. Create initial distribution in bounding box - %% (equilateral triangles) - for blk = 1 : noblks - if blk == noblks - ed = bbox_l(2,1); - end - ys = bbox_l(1,2); - ny = floor(1e3*m_lldist(repmat(0.5*(st+ed),2,1),... - [ys;bbox_l(2,2)])/h0_l); - dy = diff(bbox_l(:,2))/ny; - ns = 1; - % start at lower left and make grid going up to - % north latitude - for ii = 1:ny - if st*ed < 0 - nx = floor(1e3*m_lldist([st;0],... - [ys;ys])/(2/sqrt(3)*h0_l)) + ... - floor(1e3*m_lldist([0;ed],... - [ys;ys])/(2/sqrt(3)*h0_l)); - else - nx = floor(1e3*m_lldist([st;ed],... - [ys;ys])/(2/sqrt(3)*h0_l)); - end - ne = ns+nx-1; - if mod(ii,2) == 0 - % no offset - x(ns:ne) = linspace(st,ed,nx); - else - % offset - dx = (ed-st)/nx; - x(ns:ne) = linspace(st+0.5*dx,ed,nx); - end - y(ns:ne) = ys; - ns = ne+1; ys = ys + dy; - end - st = ed; - ed = st + blklen; - p1 = [x(:) y(:)]; clear x y - - %% 2. Remove points outside the region, apply the rejection method - p1 = p1(feval(obj.fd,p1,obj,box_num) < geps,:); % Keep only d<0 points - r0 = 1./feval(fh_l,p1).^2; % Probability to keep point - p1 = p1(rand(size(p1,1),1) < r0/max_r0,:); % Rejection method - p = [p; p1]; % Adding p1 to p - end - if box_num == 1 - % add points along the outermost polygon to fill - % outer extent more quickly. - outer_temp = obj.outer{1}; - Inan = find(isnan(outer_temp(:,1)),1,'first'); - p1 = outer_temp(1:Inan-1,:); - p1 = p1(feval(obj.fd,p1,obj,box_num) < geps,:); % Keep only d<0 points - r0 = 1./feval(fh_l, p1).^2; % Probability to keep point - p1 = p1(rand(size(p1,1),1) < r0/max_r0,:); % Rejection method - p = [p; p1]; % Adding p1 to p - end - end - else - disp('User-supplied initial points!'); - obj.grd.b = []; - 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 - if negfix > 0 - if length(obj.fixboxes)==1 && obj.fixboxes(1)==0 - obj.fixboxes(1)=1 ; - end - pfixkeep = setdiff([1:nfix]',unique(obj.egfix(:))); - % remove bars if midpoint is outside domain - egfix_mid = (obj.pfix(obj.egfix(:,1),:) + ... - obj.pfix(obj.egfix(:,2),:))/2; - for jj = 1 : length(obj.fixboxes) - if obj.fixboxes(jj) - iboubox = obj.boubox{jj}; - inbar(:,jj) = inpoly(egfix_mid,iboubox(1:end-1,:)); - end - end - inbar = sum(inbar,2) ; - obj.egfix(~inbar,:) = []; - tmppfix = obj.pfix([unique(obj.egfix(:)); pfixkeep],:); - obj.pfix = tmppfix; - obj.egfix = renumberEdges(obj.egfix); - negfix = size(obj.egfix,1); % Number of edge constraints. - end - if nfix > 0 - if length(obj.fixboxes)==1 && obj.fixboxes(1)==0 - obj.fixboxes(1)=1 ; - end - % remove pfix if outside domain - for jj = 1 : length(obj.fixboxes) - if obj.fixboxes(jj) - inbox(:,jj) = inpoly(obj.pfix,obj.boubox{jj}(1:end-1,:)); - end - end - inbox = sum(inbox,2) ; - inbox(unique(obj.egfix(:))) = 1; - obj.pfix(~inbox,:) = []; - nfix = size(obj.pfix,1); % Number of fixed points - end - if nfix >= 0, disp(['Using ',num2str(nfix),' fixed points.']);end - if negfix > 0 - if max(obj.egfix(:)) > length(obj.pfix) - error('FATAL: egfix does index correcty into pfix.'); - end - disp(['Using ',num2str(negfix),' fixed edges.']); - end - - if ~isempty(obj.pfix); p = [obj.pfix; p]; end - N = size(p,1); % Number of points N - disp(['Number of initial points after rejection is ',num2str(N)]); - %% Iterate - pold = inf; % For first iteration - if obj.plot_on >= 1 - clf,view(2),axis equal; - end - toc - fprintf(1,' ------------------------------------------------------->\n') ; - disp('Begin iterating...'); - while 1 - tic - if ~mod(it,obj.nscreen) && delIT == 0 - disp(['Iteration = ' num2str(it)]) ; - end - - % 3. Retriangulation by the Delaunay algorithm - if max(sqrt(sum((p(1:size(pold,1),:)-pold).^2,2))/h0_l*111e3) > ttol % Any large movement? - p = fixmesh(p); % Ensure only unique points. - N = size(p,1); pold = p; % Save current positions - [t,p] = delaunay_elim(p,obj.fd,geps,0); % Delaunay with elimination - - if isempty(t) - disp('Exiting') - return - end - - % Getting element quality and check "goodness" - if exist('pt','var'); clear pt; end - [pt(:,1),pt(:,2)] = m_ll2xy(p(:,1),p(:,2)); - tq = gettrimeshquan( pt, t); - mq_m = mean(tq.qm); - mq_l = min(tq.qm); - mq_s = std(tq.qm); - mq_l3sig = mq_m - 3*mq_s; - obj.qual(it,:) = [mq_m,mq_l3sig,mq_l]; - - - % If mesh quality went down "significantly" since last iteration - % ..or.. - % If not allowing improvements with reduction in quality - % then if the number of points significantly decreased - % due to a mesh improvement iteration, then rewind. - if ~mod(it,imp+1) && ((obj.qual(it,1) - obj.qual(it-1,1) < -0.10) || ... - (~obj.improve_with_reduced_quality && ... - (N - length(p_before_improve))/length(p_before_improve) < -0.10)) - disp('Mesh improvement was unsuccessful...rewinding...'); - p = p_before_improve; - N = size(p,1); % Number of points changed - pold = inf; - it = it + 1; - continue - else - N = size(p,1); pold = p; % Assign number of points and save current positions - end - % 4. Describe each bar by a unique pair of nodes. - bars = [t(:,[1,2]); t(:,[1,3]); t(:,[2,3])]; % Interior bars duplicated - bars = unique(sort(bars,2),'rows'); % Bars as node pairs - - % 5. Graphical output of the current mesh - if obj.plot_on >= 1 && (mod(it,obj.nscreen)==0 || it == 1) - cla,m_triplot(p(:,1),p(:,2),t) - m_grid - title(['Iteration = ',num2str(it)]); - if negfix > 0 - m_plot(reshape(obj.pfix(obj.egfix,1),[],2)',... - reshape(obj.pfix(obj.egfix,2),[],2)','r-') - end - if nfix > 0 - m_plot(obj.pfix(:,1),obj.pfix(:,2),'b.') - end - plt = cell2mat(obj.boubox'); - % reduce point spacing for asecthics - [plt2(:,2),plt2(:,1)] = my_interpm(plt(:,2),plt(:,1),0.1) ; - hold on ; axis manual - m_plot(plt2(:,1),plt2(:,2),'g','linewi',2) - drawnow - end - end - - % Getting element quality and check goodness - if exist('pt','var'); clear pt; end - [pt(:,1),pt(:,2)] = m_ll2xy(p(:,1),p(:,2)); - tq = gettrimeshquan( pt, t); - mq_m = mean(tq.qm); - mq_l = min(tq.qm); - mq_s = std(tq.qm); - mq_l3sig = mq_m - 3*mq_s; - obj.qual(it,:) = [mq_m,mq_l3sig,mq_l]; - - % Improve the quality of triangles next to fixed edges by - % deleting the point part of thin triangles without the fixed - % point in it. Thin triangles have poor geometric quality < - % 10%. - if ~isempty(obj.egfix) && ~mod(it,delImp) - del = heal_fixed_edges(p,t,obj.egfix) ; - if ~isempty(del) - delIT = delIT + 1 ; - if delIT < 5 - p(del,:)= []; - pold = inf; - disp(['Deleting ',num2str(length(del)),... - ' points close to fixed edges']); - continue; - else - % Abandon strategy..if it will not terminate - disp('Moving to next iteration'); - end - end - delIT = 0 ; - end - - % Termination quality, mesh quality reached is copacetic. - qual_diff = mq_l3sig - obj.qual(max(1,it-imp),2); - if ~mod(it,imp) - if abs(qual_diff) < obj.qual_tol - % Do the final elimination of small connectivity - if obj.delaunay_elim_on_exit - [t,p] = delaunay_elim(p,obj.fd,geps,1); - end - disp('Quality of mesh is good enough, exit') - close all; - break; - end - end - - % Saving a temp mesh - if ~mod(it,obj.nscreen) && delIT == 0 - disp(['Number of nodes is ' num2str(length(p))]) - disp(['Mean mesh quality is ' num2str(mq_m)]) - disp(['Min mesh quality is ' num2str(mq_l)]) - disp(['3rd sigma lower mesh quality is ' num2str(mq_l3sig)]) - tempp = p; tempt = t; - save('Temp_grid.mat','it','tempp','tempt'); - clearvars tempp tempt - end - - % 6. Move mesh points based on bar lengths L and forces F - barvec = pt(bars(:,1),:)- pt(bars(:,2),:); % List of bar vectors - if strcmp(obj.grd.proj.name,'UTM') - % UTM is already in meters (useful for small domains) - L = sqrt(sum(barvec.^2,2))*Re; - else - % Get spherical earth distances - long = zeros(length(bars)*2,1); - lat = zeros(length(bars)*2,1); - long(1:2:end) = p(bars(:,1),1); - long(2:2:end) = p(bars(:,2),1); - lat(1:2:end) = p(bars(:,1),2); - lat(2:2:end) = p(bars(:,2),2); - L = m_lldist(long,lat); L = L(1:2:end)*1e3; % L = Bar lengths in meters - end - ideal_bars = 0.5*(pt(bars(:,1),:) + pt(bars(:,2),:)); % Used to determine what bars are in bbox - [ideal_bars(:,1),ideal_bars(:,2)] = ... % needs to be in non-projected - m_xy2ll(ideal_bars(:,1),ideal_bars(:,2)); % coordinates - hbars = 0*ideal_bars(:,1); - - for box_num = 1:length(obj.h0) % For each bbox, find the bars that are in and calculate - if ~iscell(obj.fh) % their ideal lengths. - fh_l = obj.fh; - else - fh_l = obj.fh{box_num}; - end - h0_l = obj.h0(box_num); - if box_num > 1 - h0_l = h0_l/111e3; % create buffer to evalulate fh between nests - iboubox = obj.boubox{box_num}(1:end-1,:) ; - inside = inpoly(ideal_bars,iboubox) ; - else - inside = true(size(hbars)); - end - hbars(inside) = feval(fh_l,ideal_bars(inside,:)); % Ideal lengths - end - - L0 = hbars*Fscale*median(L)/median(hbars); % L0 = Desired lengths using ratio of medians scale factor - LN = L./L0; % LN = Normalized bar lengths - - % Mesh improvements (deleting and addition) - if ~mod(it,imp) - p_before_improve = p; - nn = []; pst = []; - if abs(qual_diff) < imp*obj.qual_tol && ... - (obj.improve_with_reduced_quality || qual_diff > 0) - % Remove elements with small connectivity - nn = get_small_connectivity(p,t); - disp(['Deleting ' num2str(length(nn)) ' due to small connectivity']) - - % Remove points that are too close (< LN = 0.5) - if any(LN < 0.5) - % Do not delete pfix too close. - nn1 = setdiff(reshape(bars(LN < 0.5,:),[],1),[(1:nfix)']); - disp(['Deleting ' num2str(length(nn1)) ' points too close together']) - nn = unique([nn; nn1]); - end - - % Split long edges however many times to - % better lead to LN of 1 - if any(LN > 2) - nsplit = floor(LN); - nsplit(nsplit < 1) = 1; - adding = 0; - % Probably we can just split once? - for jj = 2:2 - il = find(nsplit >= jj); - xadd = zeros(length(il),jj-1); - yadd = zeros(length(il),jj-1); - for jjj = 1 : length(il) - deltax = (p(bars(il(jjj),2),1)- p(bars(il(jjj),1),1))/jj; - deltay = (p(bars(il(jjj),2),2)- p(bars(il(jjj),1),2))/jj; - xadd(jjj,:) = p(bars(il(jjj),1),1) + (1:jj-1)*deltax; - yadd(jjj,:) = p(bars(il(jjj),1),2) + (1:jj-1)*deltay; - end - pst = [pst; xadd(:) yadd(:)]; - adding = numel(xadd) + adding; - end - disp(['Adding ',num2str(adding) ,' points.']) - end - end - if ~isempty(nn) || ~isempty(pst) - % Doing the actual subtracting and add - p(nn,:)= []; - p = [p; pst]; - pold = inf; - it = it + 1; - continue; - end - end - - F = (1-LN.^4).*exp(-LN.^4)./LN; % Bessens-Heckbert edge force - Fvec = F*[1,1].*barvec; - - Ftot = full(sparse(bars(:,[1,1,2,2]),ones(size(F))*[1,2,1,2],[Fvec,-Fvec],N,2)); - Ftot(1:nfix,:) = 0; % Force = 0 at fixed points - pt = pt + deltat*Ftot; % Update node positions - - [p(:,1),p(:,2)] = m_xy2ll(pt(:,1),pt(:,2)); - - %7. Bring outside points back to the boundary - d = feval(obj.fd,p,obj,[],1); ix = d > 0; % Find points outside (d>0) - ix(1:nfix) = 0; - if sum(ix) > 0 - pn = p(ix,:) + deps; - dgradx = (feval(obj.fd,[pn(:,1),p(ix,2)],obj,[])...%,1)... - -d(ix))/deps; % Numerical - dgrady = (feval(obj.fd,[p(ix,1),pn(:,2)],obj,[])...%,1)... - -d(ix))/deps; % gradient - dgrad2 = dgradx.^+2 + dgrady.^+2; - p(ix,:) = p(ix,:) - [d(ix).*dgradx./dgrad2,... - d(ix).*dgrady./dgrad2]; - end - - % 8. Termination criterion: Exceed itmax - it = it + 1 ; - - if ( it > obj.itmax ) - % Do the final deletion of small connectivity - if obj.delaunay_elim_on_exit - [t,p] = delaunay_elim(p,obj.fd,geps,1); - end - disp('too many iterations, exit') - close all; - break ; - end - toc - end - %% - warning('on','all') - %% - disp('Finished iterating...'); - fprintf(1,' ------------------------------------------------------->\n') ; - - %% Doing the final cleaning and fixing to the mesh... - % Clean up the mesh if specified - if ~strcmp(obj.cleanup,'none') - % Put the mesh class into the grd part of meshgen and clean - obj.grd.p = p; obj.grd.t = t; - [obj.grd,qout] = clean(obj.grd,obj.cleanup,... - 'nscreen',obj.nscreen,'djc',obj.dj_cutoff,... - 'pfix',obj.pfix); - obj.grd.pfix = obj.pfix ; - obj.grd.egfix= obj.egfix ; - obj.qual(end+1,:) = qout; - else - % Fix mesh on the projected space - [p(:,1),p(:,2)] = m_ll2xy(p(:,1),p(:,2)); - [p,t] = fixmesh(p,t); - [p(:,1),p(:,2)] = m_xy2ll(p(:,1),p(:,2)); - % Put the mesh class into the grd part of meshgen - obj.grd.p = p; obj.grd.t = t; - obj.grd.pfix = obj.pfix ; - obj.grd.egfix= obj.egfix ; - end - - % Check element order, important for the global meshes crossing - % -180/180 boundary - obj.grd = CheckElementOrder(obj.grd); - - if obj.plot_on - figure; plot(obj.qual,'linewi',2); - hold on - % plot the line dividing cleanup and distmesh - plot([it it],[0 1],'--k') - xticks(1:5:obj.itmax); - xlabel('Iterations'); ylabel('Geometric element quality'); - title('Geometric element quality with iterations'); - set(gca,'FontSize',14); - legend('q_{mean}','q_{mean}-q_{3\sigma}', 'q_{min}','Location','best'); - grid minor - end - return; - %%%%%%%%%%%%%%%%%%%%%%%%%% - % Auxiliary subfunctions % - %%%%%%%%%%%%%%%%%%%%%%%%%% - - function [t,p] = delaunay_elim(p,fd,geps,final) - % Removing mean to reduce the magnitude of the points to - % help the convex calc - if exist('pt1','var'); clear pt1; end - [pt1(:,1),pt1(:,2)] = m_ll2xy(p(:,1),p(:,2)); - if isempty(obj.egfix) - p_s = pt1 - repmat(mean(pt1),[N,1]); - TR = delaunayTriangulation(p_s); - else - TR = delaunayTriangulation(pt1(:,1),pt1(:,2),obj.egfix); - pt1 = TR.Points; - end - for kk = 1:final+1 - if kk > 1 - % Perform the following below upon exit from the mesh - % generation algorithm - nn = get_small_connectivity(pt1,t); - nn1 = heal_fixed_edges(pt1,t,obj.egfix) ; - nn = unique([nn; nn1]) ; - TR.Points(nn,:) = []; - pt1(nn,:) = []; - end - t = TR.ConnectivityList; - pmid = squeeze(mean(reshape(pt1(t,:),[],3,2),2)); % Compute centroids - [pmid(:,1),pmid(:,2)] = m_xy2ll(pmid(:,1),pmid(:,2)); % Change back to lat lon - t = t(feval(fd,pmid,obj,[]) < -geps,:); % Keep interior triangles - if kk == 1 - % Deleting very straight triangles - tq_n = gettrimeshquan( pt1, t); - bad_ele = any(tq_n.vang < 1*pi/180 | ... - tq_n.vang > 179*pi/180,2); - t(bad_ele,:) = []; - end - end - if length(pt1) ~= length(p) - clear p - [p(:,1),p(:,2)] = m_xy2ll(pt1(:,1),pt1(:,2)); - end - end - - function nn = get_small_connectivity(p,t) - % Get node connectivity (look for 4) - [~, enum] = VertToEle(t); - % Make sure they are not boundary nodes - bdbars = extdom_edges2(t, p); - bdnodes = unique(bdbars(:)); - I = find(enum <= 4); - nn = setdiff(I',[(1:nfix)';bdnodes]); % and don't destroy pfix or egfix! - return; - end - - - function del = heal_fixed_edges(p,t,egfix) - % kjr april2019 - % if there's a triangle with a low geometric quality that - % contains a fixed edge, remove the non-fixed vertex - % perform this on every other iteration to allow non-fixed - % points to create equilateral triangles nearby the locked - % edge. - % returns points IDs that should be deleted. - TR = triangulation(t,p) ; - elock = edgeAttachments(TR,egfix) ; - tq = gettrimeshquan(p,t); - elock = unique(cell2mat(elock')); - dmy = elock(tq.qm(elock) < 0.25); - badtria = t(dmy,:); - del = badtria(badtria > nfix) ; - end - - - end % end distmesh2d_plus - - - - end % end methods - -end % end class +classdef meshgen + % MESHGEN: Mesh generation class + % Handles input parameters to create a meshgen class object that can be + % used to build a msh class. + % Copyright (C) 2018 Keith Roberts & William Pringle + % + % This program is free software: you can redistribute it and/or modify + % it under the terms of the GNU General Public License as published by + % the Free Software Foundation, either version 3 of the License, or + % (at your option) any later version. + % + % This program is distributed in the hope that it will be useful, + % but WITHOUT ANY WARRANTY; without even the implied warranty of + % MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + % GNU General Public License for more details. + % + % You should have received a copy of the GNU General Public License + % along with this program. If not, see . + % + % 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 + % delaunay_elim_on_exit % whether or not to run delaunay_elim on exit of meshgen + % improve_with_reduced_quality % whether or not to allow mesh improvements with decreases in mesh quality + % improve_boundary % run a gradient desc. to improve the boundary conformity + % high_fidelity % flag to form pfix and egfix for this domain + % + properties + fd % handle to distance function + fh % handle to edge function + h0 % minimum edge length + edgefx % edgefx class + bbox % bounding box [xmin,ymin; xmax,ymax] + pfix % fixed node positions (nfix x 2 ) + egfix % edge constraints + fixboxes % a flag that indicates which boxes will use fixed constraints + plot_on % flag to plot (def: 1) or not (0) + nscreen % how many it to plot and write temp files (def: 5) + bou % geodata class + ef % edgefx class + itmax % maximum number of iterations. + 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. + 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 + anno % Approx. Nearest Neighbor search object. + 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 + improve_boundary % improve the boundary representation + high_fidelity % flag to form pfix and egfix for this domain + delaunay_elim_on_exit % whether or not to run delaunay_elim on exit of meshgen + improve_with_reduced_quality % whether or not to allow mesh improvements with decreases in mesh quality + end + + + + + methods + + + function obj = plot(obj) + if ~isempty(obj.pfix) && ~isempty(obj.egfix) + figure; hold on; % Initialize figure and hold for multiple plots + + % Plot the bounding boxes for visual aid + for box_number = 1:length(obj.boubox) + iboubox = obj.boubox{box_number}; + plot(iboubox(:,1),iboubox(:,2), 'g-', 'LineWidth', 2, 'DisplayName', 'Bounding Box'); + + % Plot outer boundaries + touter = obj.outer(box_number); + if ~isempty(touter) + tedges = Get_poly_edges(touter{1}); % Assuming Get_poly_edges is defined elsewhere + [touter, ~] = filter_polygon_constraints(touter{1}, tedges, obj.boubox, box_number); + plot(touter(:,1), touter(:,2), 'bx', 'DisplayName', 'Outer Boundary'); + end + + % Plot inner boundaries + tinner = obj.inner(box_number); + if ~isempty(tinner) && ~isempty(tinner{1}) + tedges = Get_poly_edges(tinner{1}); % Assuming Get_poly_edges is defined elsewhere + [tinner, ~] = filter_polygon_constraints(tinner{1}, tedges, obj.boubox, box_number); + plot(tinner(:,1), tinner(:,2), 'rx', 'DisplayName', 'Inner Boundary'); + end + end + + % Plot the fixed constraints as squares and edges in black + if exist('drawedge2', 'file') == 2 % Check if drawedge2 exists + drawedge2(obj.pfix, obj.egfix, [0,0,0]); % Assuming drawedge2 is a custom function + else + % Alternative plotting if drawedge2 is not available + plot(obj.pfix(:,1), obj.pfix(:,2), 'ks', 'MarkerSize', 8, 'MarkerFaceColor', 'k', 'DisplayName', 'Fixed Points'); + for i = 1:size(obj.egfix, 1) + plot(obj.pfix(obj.egfix(i,:), 1), obj.pfix(obj.egfix(i,:), 2), 'k-', 'LineWidth', 2, 'DisplayName', 'Fixed Edges'); + end + end + + axis equal; + title('Mesh Generation Constraints are black lines'); + xlabel('Longitude'); + ylabel('Latitude'); + legend('show', 'Location', 'bestoutside'); + else + disp('No constraints to plot!'); + end + end + + + + + % class constructor/default grd generation options + function obj = meshgen(varargin) + % Check for m_map dir + M_MAP_EXISTS=0; + if exist('m_proj','file')==2 + M_MAP_EXISTS=1 ; + end + if M_MAP_EXISTS~=1 + error('Where''s m_map? Chief, you need to read the user guide') + end + + + % Check for utilties dir + UTIL_DIR_EXISTS=0 ; + if exist('inpoly.m','file') + UTIL_DIR_EXISTS=1; + end + if UTIL_DIR_EXISTS~=1 + error('Where''s the utilities directory? Chief, you need to read the user guide') + end + + + p = inputParser; + % unpack options and set default ones, catch errors. + + + defval = 0; % placeholder value if arg is not passed. + % add name/value pairs + addOptional(p,'h0',defval); + addOptional(p,'bbox',defval); + addOptional(p,'fh',defval); + addOptional(p,'pfix',defval); + addOptional(p,'egfix',defval); + addOptional(p,'fixboxes',defval); + addOptional(p,'inner',defval); + addOptional(p,'outer',defval); + addOptional(p,'mainland',defval); + addOptional(p,'bou',defval); + addOptional(p,'ef',defval); + addOptional(p,'plot_on',defval); + addOptional(p,'nscreen',defval); + addOptional(p,'itmax',defval); + addOptional(p,'memory_gb',1); + addOptional(p,'cleanup',1); + addOptional(p,'direc_smooth',1); + addOptional(p,'dj_cutoff',0.25); + addOptional(p,'big_mesh',defval); + addOptional(p,'proj',defval); + addOptional(p,'qual_tol',defval); + addOptional(p,'enforceWeirs',1); + addOptional(p,'enforceMin',1); + addOptional(p,'delaunay_elim_on_exit',1); + addOptional(p,'improve_with_reduced_quality',0); + + + % parse the inputs + parse(p,varargin{:}); + + + %if isempty(varargin); return; end + % store the inputs as a struct + inp = p.Results; + + + % 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','enforceMin',... + 'delaunay_elim_on_exit','improve_with_reduced_quality',... + 'fh',... + 'inner','outer','mainland',... + 'bou','ef',... %<--OceanMesh classes come after + 'egfix','pfix','fixboxes',... + 'plot_on','nscreen','itmax',... + 'memory_gb','qual_tol','cleanup',... + 'direc_smooth','dj_cutoff',... + 'big_mesh','proj'}); + % get the fieldnames of the edge functions + fields = fieldnames(inp); + % loop through and determine which args were passed. + % also, assign reasonable default values if some options were + % not assigned. + for i = 1 : numel(fields) + type = fields{i}; + switch type + % parse aux options first + case('h0') + % Provide in meters + obj.h0 = inp.(fields{i}); + case('fh') + if isa(inp.(fields{i}),'function_handle') + obj.fh = inp.(fields{i}); + end + % can't check for errors here yet. + case('bbox') + obj.bbox= inp.(fields{i}); + if iscell(obj.bbox) + % checking bbox extents + ob_min = obj.bbox{1}(:,1); + ob_max = obj.bbox{1}(:,2); + for ii = 2:length(obj.bbox) + if any(obj.bbox{ii}(:,1) < ob_min) || ... + any(obj.bbox{ii}(:,2) > ob_max) + error(['Outer bbox must contain all ' ... + 'inner bboxes: inner box #' ... + num2str(ii) ' violates this']) + end + end + end + + + % if user didn't pass anything explicitly for + % bounding box make it empty so it can be populated + % from ef as a cell-array + if obj.bbox(1)==0 + obj.bbox = []; + end + case('pfix') + obj.pfix= inp.(fields{i}); + if obj.pfix(1)~=0 + obj.pfix(:,:) = inp.(fields{i}); + else + obj.pfix = []; + end + if obj.enforceWeirs + for j = 1 : length(obj.bou) + if ~isempty(obj.bou{j}.weirPfix) + obj.pfix = [obj.pfix ; obj.bou{j}.weirPfix]; + end + end + end + case('egfix') + obj.egfix= inp.(fields{i}); + if ~isempty(obj.egfix) && obj.egfix(1)~=0 + obj.egfix = inp.(fields{i}); + else + obj.egfix = []; + end + if obj.enforceWeirs + for j = 1 : length(obj.bou) + 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; + end + end + end + 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 + + obj.outer = {} ; + obj.inner = {} ; + obj.mainland = {} ; + + obj.bou = inp.(fields{i}); + + % handle when not a cell + if ~iscell(obj.bou) + boutemp = obj.bou; + obj.bou = cell(1); + obj.bou{1} = boutemp; + end + + % then the geodata class was provide, unpack + for ee = 1:length(obj.bou) + try + arg = obj.bou{ee} ; + catch + arg = obj.bou; + end + if isa(arg,'geodata') + + obj.high_fidelity{ee} = obj.bou{ee}.high_fidelity; + + obj.outer{ee} = obj.bou{ee}.outer; + obj.inner{ee} = obj.bou{ee}.inner; + + % save bathy interpolant to meshgen + if ~isempty(obj.bou{ee}.Fb) + obj.Fb{ee} = obj.bou{ee}.Fb ; + end + + if ~isempty(obj.inner{ee}) && ... + obj.inner{ee}(1)~= 0 + obj.outer{ee} = [obj.outer{ee}; + obj.inner{ee}]; + end + obj.mainland{ee} = obj.bou{ee}.mainland; + obj.boubox{ee} = obj.bou{ee}.boubox; + obj.inpoly_flip{ee} = obj.bou{ee}.inpoly_flip; + + end + end + + + case('ef') + tmp = inp.(fields{i}); + if isa(tmp, 'function_handle') + error('Please specify your edge function handle through the name/value pair fh'); + end + obj.ef = tmp; + + + % handle when not a cell + if ~iscell(obj.ef) + eftemp = obj.ef; + obj.ef = cell(1); + obj.ef{1} = eftemp; + end + + + % Gather boxes from ef class. + for ee = 1 : length(obj.ef) + if isa(obj.ef{ee},'edgefx') + obj.bbox{ee} = obj.ef{ee}.bbox; + end + end + + + % checking bbox extents + if iscell(obj.bbox) + ob_min = obj.bbox{1}(:,1); + ob_max = obj.bbox{1}(:,2); + for ii = 2:length(obj.bbox) + if any(obj.bbox{ii}(:,1) < ob_min) || ... + any(obj.bbox{ii}(:,2) > ob_max) + error(['Outer bbox must contain all ' ... + 'inner bboxes: inner box #' ... + num2str(ii) ' violates this']) + end + end + end + + + % kjr 2018 June: get h0 from edge functions + for ee = 1:length(obj.ef) + if isa(obj.ef{ee},'edgefx') + obj.h0(ee) = obj.ef{ee}.h0; + end + end + + + % kjr 2018 smooth the outer automatically + if length(obj.ef) > 1 + % kjr 2020, ensure the min. sizing func is + % used + if obj.enforceMin + obj.ef = enforce_min_ef(obj.ef); + end + obj.ef = smooth_outer(obj.ef,obj.Fb); + end + + + % Save the ef interpolants into the edgefx + for ee = 1:length(obj.ef) + if isa(obj.ef{ee},'edgefx') + obj.fh{ee} = @(p)obj.ef{ee}.F(p); + end + end + + + case('plot_on') + obj.plot_on= inp.(fields{i}); + case('nscreen') + obj.nscreen= inp.(fields{i}); + if obj.nscreen ~=0 + obj.nscreen = inp.(fields{i}); + obj.plot_on = 1; + else + obj.nscreen = 5; % default + end + case('itmax') + obj.itmax= inp.(fields{i}); + if obj.itmax ~=0 + obj.itmax = inp.(fields{i}); + else + obj.itmax = 100; + warning('No itmax specified, itmax set to 100'); + end + case('qual_tol') + obj.qual_tol = inp.(fields{i}); + if obj.qual_tol ~=0 + obj.qual_tol = inp.(fields{i}); + else + obj.qual_tol = 0.01; + end + case('inner') + if ~isa(obj.bou,'geodata') + obj.inner = inp.(fields{i}); + end + case('outer') + if ~isa(obj.bou,'geodata') + obj.outer = inp.(fields{i}); + if obj.inner(1)~=0 + obj.outer = [obj.outer; obj.inner]; + end + end + case('mainland') + if ~isa(obj.bou,'geodata') + obj.mainland = inp.(fields{i}); + end + case('memory_gb') + if ~isa(obj.bou,'memory_gb') + obj.memory_gb = inp.(fields{i}); + end + case('cleanup') + obj.cleanup = inp.(fields{i}); + if isempty(obj.cleanup) || obj.cleanup == 0 + obj.cleanup = 'none'; + elseif obj.cleanup == 1 + obj.cleanup = 'default'; + end + case('dj_cutoff') + obj.dj_cutoff = inp.(fields{i}); + case('direc_smooth') + obj.direc_smooth = inp.(fields{i}); + case('proj') + obj.proj = inp.(fields{i}); + % default CPP + if obj.proj == 0; obj.proj = 'equi'; end + if ~isempty(obj.bbox) + % kjr Oct 2018 use outer coarsest box for + % multiscale meshing + lon_mi = obj.bbox{1}(1,1)-obj.h0(1)/1110; + lon_ma = obj.bbox{1}(1,2)+obj.h0(1)/1110; + lat_mi = obj.bbox{1}(2,1)-obj.h0(1)/1110; + lat_ma = obj.bbox{1}(2,2)+obj.h0(1)/1110; + else + lon_mi = -180; lon_ma = 180; + lat_mi = -90; lat_ma = 90; + end + % Set up projected space + dmy = msh() ; + dmy.p(:,1) = [lon_mi; lon_ma]; + dmy.p(:,2) = [lat_mi; lat_ma]; + del = setProj(dmy,1,obj.proj) ; + case('enforceWeirs') + obj.enforceWeirs = inp.(fields{i}); + case('enforceMin') + obj.enforceMin = inp.(fields{i}); + case('delaunay_elim_on_exit') + obj.delaunay_elim_on_exit = inp.(fields{i}); + case('improve_with_reduced_quality') + obj.improve_with_reduced_quality = inp.(fields{i}); + end + end + + + if isempty(varargin); return; end + + + % error checking + if isempty(obj.boubox) && ~isempty(obj.bbox) + % Make the bounding box 5 x 2 matrix in clockwise order if + % it isn't present. This case must be when the user is + % manually specifying the PSLG. + obj.boubox{1} = [obj.bbox(1,1) obj.bbox(2,1); + obj.bbox(1,1) obj.bbox(2,2); ... + obj.bbox(1,2) obj.bbox(2,2); + 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 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). + % kjr build ANN object into meshgen + obj = createANN(obj) ; + + + global MAP_PROJECTION MAP_COORDS MAP_VAR_LIST + obj.grd.proj = MAP_PROJECTION ; + obj.grd.coord = MAP_COORDS ; + obj.grd.mapvar = MAP_VAR_LIST ; + + % Check and update geometry for high fidelity option + tpfix = []; + tegfix = []; + for box_num = 1:length(obj.h0) + % Define polygons based on available geometry data + ml = obj.mainland{box_num}; + il = obj.inner{box_num}; + polys = {}; + if ~isempty(ml), polys{end+1} = ml; end + + % High fidelity - formation of point & edge constraints + if obj.high_fidelity{box_num} + if obj.cleanup == 1 + warning('Setting cleanup to 0 since high_fidelity mode is on'); + obj.cleanup = 0; + end + disp(['Redistributing vertices for box #' num2str(box_num)]); + + % Convert polygon data, removing NaNs and duplicates + poly = cell2mat(polys'); + D = nandelim_to_cell(poly); % Custom function to handle NaN-delimited cells + D = D(~cellfun(@(p) all(isnan(p(:))), D)); % Remove cells with all NaNs + areas = cellfun(@(d) polyarea(d(:, 1), d(:, 2)), D); + [~, uniqueIdx] = uniquetol(areas, 1e-12); % Remove duplicate polygons + D = D(uniqueIdx); + + % Process each polygon + for i = 1:length(D) + my_poly = D{i}; + if size(my_poly, 1) > 2 + tedges = Get_line_edges(my_poly); + [pts, bnde] = filter_polygon_constraints(my_poly, tedges, obj.boubox, box_num); + if isempty(bnde), continue; end + poly_split = extdom_polygon(bnde, pts, 0, 1); + + for ii = 1:length(poly_split) + points = unique(poly_split{ii}, 'rows', 'stable'); + % Option 2 simply constrains the vector as + % it is in the file. + if obj.high_fidelity{box_num} == 2 + tmp_pfix = points; + tmp_egfix = Get_poly_edges([points; NaN,NaN]); + % Option 1: resamples based on the local + % mesh resolution + elseif obj.high_fidelity{box_num} == 1 + [tmp_pfix, tmp_egfix] = mesh1d(points, obj.fh, obj.h0./111e3,... + [], obj.boubox, box_num, []); + end + + if size(tmp_pfix, 1) > 2 + [tmp_pfix, tmp_egfix] = fixgeo2(tmp_pfix, tmp_egfix); + if max(tmp_egfix(:)) ~= size(tmp_pfix, 1), continue; end + tpfix = [tpfix; tmp_pfix]; + if isempty(tegfix) + tegfix = tmp_egfix; + else + tegfix = [tegfix; tmp_egfix + max(tegfix(:))]; + end + end + end + end + end + end + end + + % Final geometry check and update + [tpfix, tegfix] = fixgeo2(tpfix, tegfix); + checkfixp = setdiff(tpfix, fixmesh(tpfix), 'rows'); + if ~isempty(checkfixp), error('Duplicate fixed points detected, cannot proceed'); end + + % Update object with user-passed fixed points + obj.pfix = [obj.pfix; tpfix]; + + if isempty(obj.egfix) + obj.egfix = tegfix; + else + obj.egfix = [obj.egfix; tegfix + max(obj.egfix(:))]; + end + + if ~isempty(obj.pfix) + disp(['Using ', num2str(length(obj.pfix)), ' fixed points.']); + end + + if ~isempty(obj.egfix) + if max(obj.egfix(:)) > length(obj.pfix) + error('FATAL: egfix does not index correctly into pfix.'); + end + disp(['Using ', num2str(length(obj.egfix)), ' fixed edges.']); + warning('Please check fixed constraints: plot(mshopts)'); + end + end + + + % Creates Approximate nearest neighbor objects on start-up + function obj = createANN(obj) + box_vec = 1:length(obj.bbox); + for box_num = box_vec + if ~iscell(obj.outer) + dataset = obj.outer; + dataset(isnan(obj.outer(:,1)),:) = []; + else + dataset = obj.outer{box_num}; + dataset(isnan(obj.outer{box_num}(:,1)),:) = []; + end + if all(abs(obj.bbox{box_num}(1,:)) == 180) + % This line removes the line that can appear in the + % center for a global mesh + dataset(abs(dataset(:,1)) > 180-1e-6,:) = []; + dataset(abs(dataset(:,1)) < 1e-6,:) = []; + end + [dataset(:,1),dataset(:,2)] = m_ll2xy(dataset(:,1),dataset(:,2)); + dataset(isnan(dataset(:,1)),:) = []; + dmy = ann(dataset'); + obj.anno{box_num} = dmy; + obj.annData{box_num}=dataset; + end + end + + + function obj = build(obj) + % 2-D Mesh Generator using Distance Functions. + % Checking existence of major inputs + %% + warning('off','all') + %% + tic + it = 1 ; + Re = 6378.137e3; + geps = 1e-12*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) + EXIT_QUALITY = 0.30; % minimum quality necessary to terminate if iter < itmax + + % unpack initial points. + p = obj.grd.p; + if isempty(p) + disp('Forming initial point distribution...'); + % loop over number of boxes + for box_num = 1:length(obj.h0) + disp([' for box #' num2str(box_num)]); + % checking if cell or not and applying local values + h0_l = obj.h0(box_num); + max_r0 = 1/h0_l^2; + if ~iscell(obj.bbox) + bbox_l = obj.bbox'; % <--we must tranpose this! + else + bbox_l = obj.bbox{box_num}'; % <--tranpose! + end + if ~iscell(obj.fh) + fh_l = obj.fh; + else + fh_l = obj.fh{box_num}; + end + % Lets estimate the num_points the distribution will be + num_points = ceil(2/sqrt(3)*prod(abs(diff(bbox_l)))... + /(h0_l/111e3)^2); + noblks = ceil(num_points*2*8/obj.memory_gb*1e-9); + len = abs(bbox_l(1,1)-bbox_l(2,1)); + blklen = len/noblks; + st = bbox_l(1,1) ; ed = st + blklen; ns = 1; + %% 1. Create initial distribution in bounding box + %% (equilateral triangles) + for blk = 1 : noblks + if blk == noblks + ed = bbox_l(2,1); + end + ys = bbox_l(1,2); + ny = floor(1e3*m_lldist(repmat(0.5*(st+ed),2,1),... + [ys;bbox_l(2,2)])/h0_l); + dy = diff(bbox_l(:,2))/ny; + ns = 1; + % start at lower left and make grid going up to + % north latitude + for ii = 1:ny+1 + if st*ed < 0 + nx = floor(1e3*m_lldist([st;0],... + [ys;ys])/(2/sqrt(3)*h0_l)) + ... + floor(1e3*m_lldist([0;ed],... + [ys;ys])/(2/sqrt(3)*h0_l)); + + else + nx = floor(1e3*m_lldist([st;ed],... + [ys;ys])/(2/sqrt(3)*h0_l)); + end + ne = ns+nx-1; + if mod(ii,2) == 0 + % no offset + x(ns:ne) = linspace(st,ed,nx); + else + % offset + dx = (ed-st)/nx; + x(ns:ne) = linspace(st+0.5*dx,ed,nx); + end + % tolerance + if ii == (ny + 1) + y(ns:ne) = ys - eps; + else + y(ns:ne) = ys; + end + ns = ne+1; ys = ys + dy; + + end + + st = ed; + ed = st + blklen; + p1 = [x(:) y(:)]; clear x y + + + %% 2. Remove points outside the region, apply the rejection method + p1 = p1(feval(obj.fd,p1,obj,box_num) < geps,:); % Keep only d<0 points + r0 = 1./feval(fh_l,p1).^2; % Probability to keep point + p1 = p1(rand(size(p1,1),1) < r0/max_r0,:); % Rejection method + p = [p; p1]; % Adding p1 to p + end + if box_num == 1 + % add points along the outermost polygon to fill + % outer extent more quickly. + outer_temp = obj.outer{1}; + Inan = find(isnan(outer_temp(:,1)),1,'first'); + p1 = outer_temp(1:Inan-1,:); + p1 = p1(feval(obj.fd,p1,obj,1) < geps,:); % Keep only d<0 points + r0 = 1./feval(fh_l, p1).^2; % Probability to keep point + p1 = p1(rand(size(p1,1),1) < r0/max_r0,:); % Rejection method + p = [p; p1]; % Adding p1 to p + end + end + else + disp('User-supplied initial points!'); + obj.grd.b = []; + h0_l = obj.h0(end); % finest h0 (in case of a restart of meshgen.build). + end + + nfix = length(obj.pfix); negfix = length(obj.egfix); + if ~isempty(obj.pfix); p = [obj.pfix; p]; end + % kjr July 2023, set to these values for better convg. + if nfix > 0 + Fscale=1.1; + deltat=0.10; + end + + % Check if any boxes are set to high-fidelity + % If so turn off heal_fixed_edges + HIGH_FIDELITY_MODE = 0; + for i = 1 : length(obj.h0) + if obj.high_fidelity{i} + HIGH_FIDELITY_MODE = 1; + end + end + + N = size(p,1); % Number of points N + disp(['Number of initial points after rejection is ',num2str(N)]); + %% Iterate + pold = inf; % For first iteration + if obj.plot_on >= 1 + clf,view(2),axis equal; + end + toc + fprintf(1,' ------------------------------------------------------->\n') ; + disp('Begin iterating...'); + while 1 + tic + if ~mod(it,obj.nscreen) && delIT == 0 + disp(['Iteration = ' num2str(it)]) ; + end + % 3. Retriangulation by the Delaunay algorithm + if max(sqrt(sum((p(1:size(pold,1),:)-pold).^2,2))/h0_l*111e3) > ttol % Any large movement? + if it > 1 + p = fixmesh([obj.pfix; p]); + else + p = fixmesh(p); % Ensure only unique points. + end + N = size(p,1); pold = p; % Save current positions + [t,p] = delaunay_elim(p,obj.fd,geps,0); % Delaunay with elimination + if isempty(t) + disp('Exiting') + return + end + % Getting element quality and check "goodness" + if exist('pt','var'); clear pt; end + [pt(:,1),pt(:,2)] = m_ll2xy(p(:,1),p(:,2)); + tq = gettrimeshquan( pt, t); + mq_m = mean(tq.qm); + mq_l = min(tq.qm); + mq_s = std(tq.qm); + mq_l3sig = mq_m - 3*mq_s; + obj.qual(it,:) = [mq_m,mq_l3sig,mq_l]; + + % If not allowing improvements with reduction in quality + % ..or.. + % If not allowing improvements with reduction in quality + % then if the number of points significantly decreased + % due to a mesh improvement iteration, then rewind. + if ~mod(it,imp+1) && ((obj.qual(it,1) - obj.qual(it-1,1) < -0.10) || ... + (~obj.improve_with_reduced_quality && ... + (N - length(p_before_improve))/length(p_before_improve) < -0.10)) + + disp('Mesh improvement was unsuccessful...rewinding...'); + p = p_before_improve; + N = size(p,1); % Number of points changed + pold = inf; + it = it + 1; + continue + else + N = size(p,1); pold = p; % Assign number of points and save current positions + end + % 4. Describe each bar by a unique pair of nodes. + bars = [t(:,[1,2]); t(:,[1,3]); t(:,[2,3])]; % Interior bars duplicated + bars = unique(sort(bars,2),'rows'); % Bars as node pairs + % 5. Graphical output of the current mesh + if obj.plot_on >= 1 && (mod(it,obj.nscreen)==0 || it == 1) + cla,m_triplot(p(:,1),p(:,2),t) + m_grid + title(['Iteration = ',num2str(it)]); + if negfix > 0 + m_plot(reshape(obj.pfix(obj.egfix,1),[],2)',... + reshape(obj.pfix(obj.egfix,2),[],2)','r-') + end + if nfix > 0 + m_plot(obj.pfix(:,1),obj.pfix(:,2),'b.') + end + hold on ; + axis manual + drawnow + end + end + % Getting element quality and check goodness + if exist('pt','var'); clear pt; end + [pt(:,1),pt(:,2)] = m_ll2xy(p(:,1),p(:,2)); + tq = gettrimeshquan( pt, t); + mq_m = mean(tq.qm); + mq_l = min(tq.qm); + mq_s = std(tq.qm); + mq_l3sig = mq_m - 3*mq_s; + obj.qual(it,:) = [mq_m,mq_l3sig,mq_l]; + + % Improve the quality of triangles next to fixed edges by + % deleting the point part of thin triangles without the fixed + % point in it. Thin triangles have poor geometric quality < + % 10%. + if ~isempty(obj.egfix) && ~mod(it,delImp) && ~HIGH_FIDELITY_MODE + del = heal_fixed_edges(p,t,obj.egfix) ; + if ~isempty(del) + delIT = delIT + 1 ; + if delIT < 5 + p(del,:)= []; + pold = inf; + disp(['Deleting ',num2str(length(del)),... + ' points close to fixed edges']); + continue; + else + % Abandon strategy..if it will not terminate + disp('Moving to next iteration'); + end + end + delIT = 0 ; + end + + % Termination quality, mesh quality reached is copacetic. + qual_diff = mq_l3sig - obj.qual(max(1,it-imp),2); + if ~mod(it,imp) + if mq_l > EXIT_QUALITY + % Do the final elimination of small connectivity + if obj.delaunay_elim_on_exit + end + disp('Quality of mesh is good enough, exit') + close all; + break; + end + end + % Saving a temp mesh + if ~mod(it,obj.nscreen) && delIT == 0 + disp(['Number of nodes is ' num2str(length(p))]) + disp(['Mean mesh quality is ' num2str(mq_m)]) + disp(['Min mesh quality is ' num2str(mq_l)]) + disp(['3rd sigma lower mesh quality is ' num2str(mq_l3sig)]) + tempp = p; tempt = t; + save('Temp_grid.mat','it','tempp','tempt'); + clearvars tempp tempt + end + % 6. Move mesh points based on bar lengths L and forces F + barvec = pt(bars(:,1),:)- pt(bars(:,2),:); % List of bar vectors + if strcmp(obj.grd.proj.name,'UTM') + % UTM is already in meters (useful for small domains) + L = sqrt(sum(barvec.^2,2))*Re; + else + % Get spherical earth distances + long = zeros(length(bars)*2,1); + lat = zeros(length(bars)*2,1); + long(1:2:end) = p(bars(:,1),1); + long(2:2:end) = p(bars(:,2),1); + lat(1:2:end) = p(bars(:,1),2); + lat(2:2:end) = p(bars(:,2),2); + L = m_lldist(long,lat); L = L(1:2:end)*1e3; % L = Bar lengths in meters + end + ideal_bars = 0.5*(pt(bars(:,1),:) + pt(bars(:,2),:)); % Used to determine what bars are in bbox + [ideal_bars(:,1),ideal_bars(:,2)] = ... % needs to be in non-projected + m_xy2ll(ideal_bars(:,1),ideal_bars(:,2)); % coordinates + hbars = 0*ideal_bars(:,1); + + + for box_num = 1:length(obj.h0) % For each bbox, find the bars that are in and calculate + if ~iscell(obj.fh) % their ideal lengths. + fh_l = obj.fh; + else + fh_l = obj.fh{box_num}; + end + h0_l = obj.h0(box_num); + if box_num > 1 + h0_l = h0_l/111e3; % create buffer to evalulate fh between nests + iboubox = obj.boubox{box_num}(1:end-1,:) ; + inside = inpoly(ideal_bars,iboubox) ; + else + inside = true(size(hbars)); + end + hbars(inside) = feval(fh_l,ideal_bars(inside,:)); % Ideal lengths + end + + + L0 = hbars*Fscale*median(L)/median(hbars); % L0 = Desired lengths using ratio of medians scale factor + LN = L./L0; % LN = Normalized bar lengths + + + % Mesh improvements (deleting and addition) + p_before_improve = p; + if ~mod(it,imp) % + nn = []; pst = []; + if abs(qual_diff) < imp*obj.qual_tol && ... + (obj.improve_with_reduced_quality || qual_diff > 0) + + % Remove elements with small connectivity + nn = get_small_connectivity(p,t); + disp(['Deleting ' num2str(length(nn)) ' due to small connectivity']) + + + % Remove points that are too close (< LN = 0.5) + if any(LN < 0.5) + % Do not delete pfix too close. + nn1 = setdiff(reshape(bars(LN < 0.5,:),[],1),[(1:nfix)']); + disp(['Deleting ' num2str(length(nn1)) ' points too close together']) + nn = unique([nn; nn1]); + end + + + % Split long edges however many times to + % better lead to LN of 1 + if any(LN > 2) + nsplit = floor(LN); + nsplit(nsplit < 1) = 1; + adding = 0; + % Split once + for jj = 2:2 + il = find(nsplit >= jj); + xadd = zeros(length(il),jj-1); + yadd = zeros(length(il),jj-1); + for jjj = 1 : length(il) + deltax = (p(bars(il(jjj),2),1)- p(bars(il(jjj),1),1))/jj; + deltay = (p(bars(il(jjj),2),2)- p(bars(il(jjj),1),2))/jj; + xadd(jjj,:) = p(bars(il(jjj),1),1) + (1:jj-1)*deltax; + yadd(jjj,:) = p(bars(il(jjj),1),2) + (1:jj-1)*deltay; + end + pst = [pst; xadd(:) yadd(:)]; + adding = numel(xadd) + adding; + end + disp(['Adding ',num2str(adding) ,' points.']) + end + end + if ~isempty(nn) || ~isempty(pst) + % Doing the actual subtracting and add + p(nn,:)= []; + p = [p; pst]; + pold = inf; + it = it + 1; + continue; + end + end + + + F = (1-LN.^4).*exp(-LN.^4)./LN; % Bessens-Heckbert edge force + F(isinf(F)) = 0; + Fvec = F*[1,1].*barvec; + + + Ftot = full(sparse(bars(:,[1,1,2,2]),ones(size(F))*[1,2,1,2],[Fvec,-Fvec],N,2)); + Ftot(1:nfix,:) = 0; % Force = 0 at fixed points + + pt = pt + deltat*Ftot; % Update node positions + + + [p(:,1),p(:,2)] = m_xy2ll(pt(:,1),pt(:,2)); + + + %7. Bring outside points back to the boundary + d = feval(obj.fd,p,obj,[],1); ix = d > 0; % Find points outside (d>0) + ix(1:nfix) = 0; + alpha = 1.0; + for ib = 1 : 1 % length(obj.improve_boundary) + if sum(ix) > 0 + pn = p(ix,:) + deps; + dgradx = (feval(obj.fd,[pn(:,1),p(ix,2)],obj,[])...%,1)... + -d(ix))/deps; % Numerical + dgrady = (feval(obj.fd,[p(ix,1),pn(:,2)],obj,[])...%,1)... + -d(ix))/deps; % gradient + dgrad2 = dgradx.^+2 + dgrady.^+2; + dgrad2(dgrad2 < eps) = eps; + p(ix,:) = p(ix,:) - alpha*[d(ix).*dgradx./dgrad2,... + d(ix).*dgrady./dgrad2]; + end + alpha = alpha / 0.5; + end + + + % 8. Termination criterion: Exceed itmax + it = it + 1 ; + + + if ( it > obj.itmax ) + % Do the final deletion of small connectivity + if obj.delaunay_elim_on_exit + end + disp('too many iterations, exit') + close all; + break ; + end + toc + end + %% + warning('on','all') + %% + disp('Finished iterating...'); + fprintf(1,' ------------------------------------------------------->\n') ; + + + %% Doing the final cleaning and fixing to the mesh... + % Always save the mesh! + save('Precleaned_grid.mat','it','p','t'); + + % Clean up the mesh if specified + if ~strcmp(obj.cleanup,'none') + % Put the mesh class into the grd part of meshgen and clean + obj.grd.p = p; obj.grd.t = t; + [obj.grd,qout] = clean(obj.grd,obj.cleanup,... + 'nscreen',obj.nscreen,'djc',obj.dj_cutoff,... + 'pfix',obj.pfix); + obj.grd.pfix = obj.pfix ; + obj.grd.egfix= obj.egfix ; + obj.grd.egfix= obj.egfix ; + obj.qual(end+1,:) = qout; + else + % Fix mesh on the projected space + [p(:,1),p(:,2)] = m_ll2xy(p(:,1),p(:,2)); + [p,t] = fixmesh(p,t); + [p(:,1),p(:,2)] = m_xy2ll(p(:,1),p(:,2)); + % Put the mesh class into the grd part of meshgen + obj.grd.p = p; obj.grd.t = t; + obj.grd.pfix = obj.pfix ; + obj.grd.egfix= obj.egfix ; + end + + + % Check element order, important for the global meshes crossing + % -180/180 boundary + obj.grd = CheckElementOrder(obj.grd); + + + if obj.plot_on + figure; plot(obj.qual,'linewi',2); + hold on + % plot the line dividing cleanup and distmesh + plot([it it],[0 1],'--k') + xticks(1:5:obj.itmax); + xlabel('Iterations'); ylabel('Geometric element quality'); + title('Geometric element quality with iterations'); + set(gca,'FontSize',14); + legend('q_{mean}','q_{mean}-q_{3\sigma}', 'q_{min}','Location','best'); + grid minor + end + return; + %%%%%%%%%%%%%%%%%%%%%%%%%% + % Auxiliary subfunctions % + %%%%%%%%%%%%%%%%%%%%%%%%%% + + + function [t,p] = delaunay_elim(p,fd,geps,final) + % Removing mean to reduce the magnitude of the points to + % help the convex calc + if exist('pt1','var'); clear pt1; end + [pt1(:,1),pt1(:,2)] = m_ll2xy(p(:,1),p(:,2)); + if isempty(obj.egfix) + p_s = pt1 - repmat(mean(pt1),[N,1]); + TR = delaunayTriangulation(p_s); + else + TR = delaunayTriangulation(pt1(:,1),pt1(:,2),obj.egfix); + pt1 = TR.Points; + end + for kk = 1:final+1 + if kk > 1 + % Perform the following below upon exit from the mesh + % generation algorithm + nn = get_small_connectivity(pt1,t); + nn1 = []; + nn = unique([nn; nn1]) ; + TR.Points(nn,:) = []; + pt1(nn,:) = []; + end + t = TR.ConnectivityList; + pmid = squeeze(mean(reshape(pt1(t,:),[],3,2),2)); % Compute centroids + [pmid(:,1),pmid(:,2)] = m_xy2ll(pmid(:,1),pmid(:,2)); % Change back to lat lon + t = t(feval(fd,pmid,obj,[]) < -geps,:); % Keep interior trianglesi + end + if length(pt1) ~= length(p) + clear p + [p(:,1),p(:,2)] = m_xy2ll(pt1(:,1),pt1(:,2)); + end + end + + + function nn = get_small_connectivity(p,t) + % Get node connectivity (look for 4) + [~, enum] = VertToEle(t); + % Make sure they are not boundary nodes + bdbars = extdom_edges2(t, p); + bdnodes = unique(bdbars(:)); + I = find(enum <= 4); + nn = setdiff(I',[(1:nfix)';bdnodes]); % and don't destroy pfix or egfix! + return; + end + + function del = heal_fixed_edges(p,t,egfix) + % kjr april2019 + % if there's a triangle with a low geometric quality that + % contains a fixed edge, remove the non-fixed vertex + % perform this on every other iteration to allow non-fixed + % points to create equilateral triangles nearby the locked + % edge. + % returns points IDs that should be deleted. + TR = triangulation(t,p) ; + elock = edgeAttachments(TR,egfix) ; + tq = gettrimeshquan(p,t); + elock = unique(cell2mat(elock')); + dmy = elock(tq.qm(elock) < 0.25); + badtria = t(dmy,:); + del = badtria(badtria > nfix) ; + end + + end % end mesh generator + + + end % end methods + + +end % end class + diff --git a/@meshgen/private/Get_line_edges.m b/@meshgen/private/Get_line_edges.m new file mode 100644 index 00000000..bb666b2f --- /dev/null +++ b/@meshgen/private/Get_line_edges.m @@ -0,0 +1,17 @@ +function edges = Get_line_edges(nodes) +%GET_LINE_EDGES Generate edges from a sequence of nodes in a closed loop. +% edges = GET_LINE_EDGES(nodes) takes a list of nodes (vertices) as input +% and generates a list of edges where each node is connected to the next, +% forming a closed loop. The output 'edges' is an Mx2 matrix, where M is +% the number of edges, and each row represents an edge defined by the +% indices of its two endpoints. + + % Number of nodes + nNodes = length(nodes); + + % Generate edges connecting each node to the next + edges = [(1:nNodes-1)', (2:nNodes)']; + + % Add the edge connecting the last node back to the first + edges(end+1, :) = [nNodes, 1]; +end diff --git a/@meshgen/private/filter_polygon_constraints.m b/@meshgen/private/filter_polygon_constraints.m new file mode 100644 index 00000000..b1c50d86 --- /dev/null +++ b/@meshgen/private/filter_polygon_constraints.m @@ -0,0 +1,46 @@ +function [pfix, egfix] = filter_polygon_constraints(pfix, egfix, ibouboxes, box_number) +%FILTER_CONSTRAINTS Removes edge constraints based on bounding box criteria. +% This function filters out edge constraints (egfix) where at least one +% endpoint (from pfix) is outside the specified bounding box (ibouboxes) +% for the given box_number. It also adjusts the edges based on nested +% bounding boxes, if applicable. +% remove bars if one point is outside +node1=pfix(egfix(:,1),:); +node2=pfix(egfix(:,2),:); +iboubox = ibouboxes{box_number}; +% to enlarge or shrink the box, you must make sure bbox is equi +% spaced +[ty,tx]=my_interpm(iboubox(:,2),iboubox(:,1),100/111e3); +iboubox = [tx,ty]; +buffer_size = 1.0; +iboubox(:,1) = buffer_size*iboubox(:,1)+(1-buffer_size)*mean(iboubox(1:end-1,1)); +iboubox(:,2) = buffer_size*iboubox(:,2)+(1-buffer_size)*mean(iboubox(1:end-1,2)); +inside_node1 = inpoly(node1,iboubox(1:end-1,:)) ; +inside_node2 = inpoly(node2,iboubox(1:end-1,:)) ; +inside = inside_node1 .* inside_node2; +% Get all points inside inner boxes and consider these outside for +% all the nested boxes. +for bn = box_number+1:length(ibouboxes) + % enlarge nest + iboubox = ibouboxes{bn}(1:end-1,:); + [ty,tx]=my_interpm(iboubox(:,2),iboubox(:,1),100/111e3); + iboubox = [tx,ty]; + iboubox(:,1) = 1.25*iboubox(:,1)+(1-1.25)*mean(iboubox(1:end-1,1)); + iboubox(:,2) = 1.25*iboubox(:,2)+(1-1.25)*mean(iboubox(1:end-1,2)); + inside_node1 = inpoly(node1,iboubox); + inside_node2 = inpoly(node2,iboubox); + inside2 = inside_node1 .* inside_node2; + inside(find(inside2)) = false; +end +egfix(~inside,:) = []; +tegfix=egfix'; +uid = unique(tegfix(:)); +tuid = length(uid); +% remove nfix outside iboubox +if tuid > 0 + % remove pfix if outside domain + pfix = pfix(uid,:); +end +egfix = renumberEdges(egfix); + +end diff --git a/@meshgen/private/mesh1d.m b/@meshgen/private/mesh1d.m new file mode 100644 index 00000000..ca9f9a46 --- /dev/null +++ b/@meshgen/private/mesh1d.m @@ -0,0 +1,207 @@ + +function [pout,t,converged]=mesh1d(poly,fh0,h,fix,boubox,box_num0,varargin) +% Mesh Generator using Distance Functions. +% [P,T]=mesh1d(FDIST,FH,H,BOX,FIX,FDISTPARAMS) +% +% POUT: Resampled Node positions (N X 2) +% T: Edge indices (NTx2)) +% POLY: The polygon (NX X 2 +% FH: Edge length function +% H: Smallest edge length +% FIX: Indices of Poly Fixed node positions (NFIX X 1) + +% Preprocessing steps - move into parametric space +%poly = add_flairs(poly, 50); + +X = poly(:,1); +Y = poly(:,2); + +% find sharp corners & add flairs +% [sharp_corners] = find_sharp_corners(poly, 40); +% if ~isempty(sharp_corners) +% fix = [fix; sharp_corners]; +% end + + +% add ending point to emulate periodic mesh generation +%X(end+1) = X(1,:) + eps; +%Y(end+1) = Y(1,:) + eps; + +xd = diff(X); +yd = diff(Y); +dist = sqrt(xd.^2 + yd.^2); +u = cumsum(dist); +u = [0; u]; +totaldist = u(end); + +np = totaldist / min(h); +% kjr check if total np will be too little +% i.e., less than 5 +if np < 5 + pout = []; + t = []; + converged=1; + return +end +t = linspace(0,max(u),ceil(np)); + +xn = interp1(u,X,t); +yn = interp1(u,Y,t); + +u_fix = []; +if ~isempty(fix) + % Convert fixed points to parametric space + % determine the distance these are points are from the zero point. + u_fix = []; + for i = 1:length(fix) + xd_tmp = diff(X(1:fix(i))); + yd_tmp = diff(Y(1:fix(i))); + dist_tmp = sqrt(xd_tmp.^2 + yd_tmp.^2); + if isempty(dist_tmp) + continue + elseif cumsum(dist_tmp) == totaldist + continue + else + csum = cumsum(dist_tmp); + end + u_fix = [u_fix; csum(end)]; + end +end + +A = 0; +B = totaldist; + +box = [A,B]'; + +fdist = @(x) my_1d_sdf(x, A, B); + +tmph = nestedHFx([xn',yn'],boubox,fh0,h); + + function [h_inner] = nestedHFx(x,boubox,fh0,h) + h_inner= x(:,1)*0; + for box_num = 1:length(h) % For each bbox, find the points that are in and calculate + fh_l = fh0{box_num}; + if box_num > 1 + iboubox = boubox{box_num}(1:end-1,:) ; + inside = inpoly(x,iboubox) ; + else + inside = true(size(h_inner)); + end + h_inner(inside) = feval(fh_l,x(inside,:))./111e3; + end + end + +F = griddedInterpolant(t,tmph,'linear','linear'); + +fh = @(x) F(x); + +dim=size(box,2); +ptol=.01; +L0mult=1+.4/2^(dim-1); +deltat=.10; geps=1e-3*h(box_num0); +deps=sqrt(eps)*h(box_num0); + +% 1. Create initial distribution in bounding box +if dim==1 + p=(box(1):h(box_num0):box(2))'; +else + disp('Only 1D is supported') + %error +end + +% 2. Remove points outside the region, apply the rejection method +p=p(feval(fdist,p)0; + %ix(1:length(fix),:) = 0; + ix(fixed_point_ind,:) = 0; + gradd=zeros(sum(ix),dim); + for ii=1:dim + a=zeros(1,dim); + a(ii)=deps; + d1x=feval(fdist,p(ix,:)+ones(sum(ix),1)*a); + gradd(:,ii)=(d1x-d(ix))/deps; + end + p(ix,:)=p(ix,:)-d(ix)*ones(1,dim).*gradd; + + % 8. Termination criterion + maxdp=max(deltat*sqrt(sum(dp(d<-geps,:).^2,2))); + + %disp(maxdp) + if maxdp 5000 + warning('Mesh1d: some line segments did not converge...') + converged=0; + pout(:,1) = interp1(u,X,p); + pout(:,2) = interp1(u,Y,p); + if any(isnan(pout(:,1))) + warning('Nans detected in 1d mesh') + end + break + end + +end +end diff --git a/@meshgen/private/my_1d_sdf.m b/@meshgen/private/my_1d_sdf.m new file mode 100644 index 00000000..1ede3d29 --- /dev/null +++ b/@meshgen/private/my_1d_sdf.m @@ -0,0 +1,10 @@ +function [dist] = my_1d_sdf(p,A,B) +distA = abs(p - A); +distB = abs(p - B); +dist = -min(distA,distB); +if p < A + dist = -1*dist; +elseif p > B + dist = -1*dist; +end +end \ No newline at end of file diff --git a/@meshgen/private/nandelim_to_cell.m b/@meshgen/private/nandelim_to_cell.m new file mode 100644 index 00000000..0c12fd0b --- /dev/null +++ b/@meshgen/private/nandelim_to_cell.m @@ -0,0 +1,28 @@ +function polygons = nandelim_to_cell(data) +%NANDELIM_TO_CELL Convert NaN-delimited data into a cell array of polygons. +% polygons = NANDELIM_TO_CELL(data) takes an Nx2 matrix 'data', where +% each polygon is separated by a row with NaN in the first column, and +% converts it into a cell array. Each cell contains the vertices of a +% polygon defined by consecutive rows up to the next NaN row. + + % Identify the indices of rows that have NaN in the first column + nanRows = find(isnan(data(:, 1))); + + % Preallocate the cell array for efficiency + polygons = cell(length(nanRows) + 1, 1); + + % Set the starting index for the first polygon + startIndex = 1; + + % Iterate over each NaN row to extract polygons + for i = 1:length(nanRows) + % Extract the polygon data up to the row before the NaN + polygons{i} = data(startIndex:nanRows(i) - 1, :); + + % Update the start index for the next polygon + startIndex = nanRows(i) + 1; + end + + % Handle the last polygon after the final NaN row + polygons{end} = data(startIndex:end, :); +end diff --git a/@meshgen/private/renumberEdges.m b/@meshgen/private/renumberEdges.m new file mode 100644 index 00000000..a227ab3e --- /dev/null +++ b/@meshgen/private/renumberEdges.m @@ -0,0 +1,22 @@ +function edges = renumberEdges(edges) +%RENUMBEREDGES Renumber edges to maintain consistency after point modifications. +% edges = RENUMBEREDGES(edges) updates the indices of edge endpoints in +% the 'edges' matrix to reflect a new, sequential numbering based on +% their occurrence in 'edges'. This is useful after points have been +% filtered or reordered, ensuring that edge references are consistent +% with the updated point list. + +% Find unique indices of edge endpoints and create a direct mapping +idx = unique(edges(:)); +for i = 1 : length(idx) + map(idx(i),1) = i; +end + +% renumber bnde to correspond with pfix +for i = 1 : size(edges,1) + edges(i,1) = map(edges(i,1)); + edges(i,2) = map(edges(i,2)); +end + + +end diff --git a/Examples/Example_13_New_Zealand_High_Fidelity.m b/Examples/Example_13_New_Zealand_High_Fidelity.m new file mode 100644 index 00000000..e687f0e4 --- /dev/null +++ b/Examples/Example_13_New_Zealand_High_Fidelity.m @@ -0,0 +1,69 @@ +% Example_13_NZ_high_fidelity: High-Fidelity Mesh Generation for the South Island of New Zealand +% This script demonstrates the application of multiscale meshing techniques with a focus on high-fidelity options. +% High-fidelity meshing is particularly crucial for areas within a larger regional model that contain engineered structures +% or coastlines. Such areas require precise representation in the mesh connectivity to ensure accurate simulation outcomes. +% This approach enables detailed analysis and modeling of geophysical processes near these critical features, enhancing the +% reliability of simulations in coastal engineering, environmental assessment, and resource management applications. +clearvars; clc; close all + +PREFIX = '13_NZ_highfidelity'; + +%% STEP 1: Define mesh boundaries and resolution parameters. +% Bounding box coordinates (longitude and latitude) define the area of interest. +bbox1 = [166, 176; -48, -40]; % lon_min, lon_max; lat_min, lat_max +min_el1 = 1e3; % Minimum element size (m). +max_el1 = 100e3; % Maximum element size (m). +max_el_ns1 = 5e3; % Maximum nearshore element size (m). +grade1 = 0.15; % Mesh grading factor. +R1 = 3; % Feature resolution factor. + +%% STEP 2: Initialize geographic data for mesh generation. +% Specify coastline dataset for mesh boundary definition. +coastline = 'GSHHS_f_L1'; +gdat{1} = geodata('shp', coastline, 'bbox', bbox1, 'h0', min_el1); + +%% STEP 3: Generate edge functions for mesh refinement. +% Edge function defines transition of mesh sizes across the domain. +fh{1} = edgefx('geodata', gdat{1}, 'fs', R1, 'max_el_ns', max_el_ns1, 'max_el', max_el1, 'g', grade1); + +%% STEP 4: Define a refined mesh area. +% Additional bounding box for targeted refinement (e.g., a specific bay or coastal region). +% In this case Christchurch Bay with a more detailed shoreline vector +% file that contains hardened shoreline features. +coastline2 = 'Coastline_epsg4326'; + + +bbox2= [172.61348324,172.85932924; + -43.67102284,-43.51383230]; + +min_el2 = 50.0; % Refined minimum element size (m). +max_el_ns2 = 250.0; % Refined maximum nearshore element size (m). +max_el2 = 500; % Refined maximum element size (m). +grade2 = 0.05; % Mesh grading factor for refined area. +R2 = 5; % Feature resolution factor for refined area. + +% High-fidelity option for enhanced mesh detail in the refined area. +gdat{2} = geodata('shp', coastline2, 'bbox', bbox2, 'h0', min_el2, 'high_fidelity', 1); +fh{2} = edgefx('geodata', gdat{2}, 'fs', R2, 'max_el_ns', max_el_ns2, 'max_el', max_el2, 'g', grade2); + +%% STEP 5: Mesh generation with defined edge functions. +% Configure mesh generation options, including visualization during the process. +mshopts = meshgen('ef', fh, 'bou', gdat, 'plot_on', 1, 'nscreen', 5); + +% You can extract the developed constraints +pfix = mshopts.pfix; +egfix = mshopts.egfix; + +mshopts = mshopts.build(); % Execute mesh building with specified options. + +%% STEP 6: Visualize and export the generated mesh. +% Extract and visualize the generated mesh, highlighting refined areas. +m = mshopts.grd; +m.plot('type', 'tri', 'proj', 'lamb'); % Plot entire mesh. +hold on; m_plot(m.pfix(:,1),m.pfix(:,2),'r.') + +m.plot('type', 'tri', 'proj', 'none', 'subdomain', bbox2); % Highlight refined area. +hold on; drawedge2(pfix, egfix, 'r'); % Overlay edge constraints. + +% Optional: Export mesh to a .2dm file for editing connectivity in GIS software like QGIS. +write(m,sprintf('%s_mesh.2dm',PREFIX),'2dm'); diff --git a/Examples/runall.m b/Examples/runall.m index 31c00823..83016063 100644 --- a/Examples/runall.m +++ b/Examples/runall.m @@ -1,4 +1,5 @@ % run all examples +run('../setup_oceanmesh2d.m') Example_1_NZ Example_2_NY @@ -13,4 +14,5 @@ Example_9_TBAY Example_10_Multiscale_Smoother Example_11_Remeshing_Patches -Example_12_Riverine_flow \ No newline at end of file +Example_12_Riverine_flow +Example_13_New_Zealand_High_Fidelity diff --git a/README.md b/README.md index 8eaff1f5..3b8b4cb8 100644 --- a/README.md +++ b/README.md @@ -1,3 +1,12 @@ +# `OceanMesh2D:` +## `Precise distance-based two-dimensional automated mesh generation toolbox intended for coast ocean/shallow water flow models` + +## IMPORTANT NOTE: +This is the legacy `MASTER` branch, and is not recommended unless for example you want to reproduce the results from the GMD paper referenced at the bottom here. Otherwise, it is recommended to use the default `PROJECTION` branch. + +## `DISCLAIMER: ` +The boundary of the meshing domain must be a polygon (first point equals the last and non-self intersecting) but it does not need to be simplified. Read the user guide for more information about the inputs. +

OceanMesh2D

Precise distance-based two-dimensional automated mesh generation toolbox intended for coastal ocean/shallow water flow models.

@@ -101,7 +110,7 @@ If you make use of `OceanMesh2D` please include a reference to [1], and to any o [1] - Roberts, K. J., Pringle, W. J., and Westerink, J. J., 2019. OceanMesh2D 1.0: MATLAB-based software for two-dimensional unstructured mesh generation in coastal ocean modeling, - Geoscientific Model Development, 12, 1847-1868. https://doi.org/10.5194/gmd-12-1847-2019. + Geosci. Model Dev. (GMD), https://doi.org/10.5194/gmd-12-1847-2019. [2] - Roberts, K. J., Pringle, W. J, 2018. OceanMesh2D: User guide - Precise distance-based two-dimensional automated mesh generation toolbox intended for coastal ocean/shallow water. https://doi.org/10.13140/RG.2.2.21840.61446/2. @@ -140,6 +149,7 @@ GALLERY:                         +        

@@ -157,6 +167,7 @@ The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.0.0/) ### Unreleased (on current HEAD of the Projection branch) ## Added +- Added ability to generate constrainsts for shorelines automatically using the new high-fidelity method. https://github.com/CHLNDDEV/OceanMesh2D/pull/264 - auxiliary file reader function that can be used standalone e.g., m = m.read({'fort.13','fort.15'})), as well as called from within msh() function. https://github.com/CHLNDDEV/OceanMesh2D/pull/287 - swanoutput namelist for SWAN model outputs. https://github.com/CHLNDDEV/OceanMesh2D/pull/287 - Added new function in `msh` called `remesh_patch` to remesh within specified polygonal domains and insert back into parent mesh. https://github.com/CHLNDDEV/OceanMesh2D/pull/301 @@ -168,6 +179,7 @@ The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.0.0/) - Option `delaunay_elim_on_exit` to `meshgen` to skip the last call to `delaunay_elim` to potentially avoid deleting boundary elements. - Geoid offset nodal attribute in `Calc_f13` subroutine. https://github.com/CHLNDDEV/OceanMesh2D/pull/251 - Support for writing Self Attraction and Loading (SAL) files in NetCDF for the ADCIRC model. https://github.com/CHLNDDEV/OceanMesh2D/pull/231 +- Added 'high-fidelity' option for automatically forming and constraining edges into the mesh. https://github.com/CHLNDDEV/OceanMesh2D/pull/264 ## Changed - treating logicals in namelists as a logicals type within MATLAB but ouput as a string. https://github.com/CHLNDDEV/OceanMesh2D/pull/287 - `Calc_Sponge.m` to allow for sponge in polygon region and coefficents based on spatially varying depth. https://github.com/CHLNDDEV/OceanMesh2D/pull/287 diff --git a/Tests/TestJBAY.m b/Tests/TestJBAY.m index 39403bcb..3cbc8822 100644 --- a/Tests/TestJBAY.m +++ b/Tests/TestJBAY.m @@ -8,18 +8,18 @@ ERR2_TOL = 0.01; QUAL_TOL = 0.25; Re2 = 111^2; -% p: [23957×2 double] -% t: [42853×3 double] -% Area = 211.19 km^2 -% Volume = 2.0782 km^3 -TARGET = 23957; +% p: [25190x2 double] +% t: [45270x3 double] +% Area = 211 km^2 +% Volume = 2.075 km^3 +TARGET = 25190; VALUE = length(m.p); if abs(VALUE - TARGET)/TARGET > ERR_TOL error(['Incorrect number of points for %s. ',... 'Got %i, expecting %i +- %4.2f%%.'],... PREFIX,VALUE,TARGET,ERR_TOL*100); end -TARGET = 42853; +TARGET = 45270; VALUE = length(m.t); if abs(VALUE - TARGET)/TARGET > ERR_TOL error(['Incorrect number of triangles for %s. ',... @@ -27,7 +27,7 @@ PREFIX,VALUE,TARGET,ERR_TOL*100); end -TARGET = 211.19; +TARGET = 211; X = m.p(:,1); Y = m.p(:,2); Area = sum(polyarea(X(m.t)',Y(m.t)').*cosd(mean(Y(m.t)')))*Re2; if abs(Area - TARGET)/TARGET > ERR2_TOL @@ -36,7 +36,7 @@ PREFIX,Area,TARGET,ERR2_TOL*100); end -TARGET = 2.0782; +TARGET = 2.075; bc = (m.b(m.t(:,1))+m.b(m.t(:,2))+m.b(m.t(:,3)))/3; bc = bc/1e3; % convert to km Volume = sum(polyarea(X(m.t)',Y(m.t)').*cosd(mean(Y(m.t)')).*bc')*Re2; diff --git a/Tests/test_1d_original.m b/Tests/test_1d_original.m new file mode 100644 index 00000000..a3babd49 --- /dev/null +++ b/Tests/test_1d_original.m @@ -0,0 +1,31 @@ +clearvars ; close all; clc ; + +minh = 0.01; +maxh = 0.20; +grade = 0.15; + +poly = [0.2483 0.4942 + 0.1838 0.6080 + 0.2851 0.7861 + 0.4741 0.8445 + 0.5616 0.8883 + 0.6953 0.7277 + 0.6745 0.3891 + 0.8474 0.3161 + 0.9211 0.6168 + 0.9764 0.6752 + 0.9741 0.1876 + 0.7967 0.0241 + 0.6838 0.0358 + 0.3658 0.1905 + 0.3520 0.345 + 0.2483 0.4942]; + +fh = @(x) min(grade*abs(x(:,1)) + minh,maxh); + +ID_pfix = []; +[p,t] = mesh1d(poly,fh,minh,[]); + +figure; drawedge2(p,t); +hold on; plot(p(:,1),p(:,2),'s','MarkerFaceColor','r'); +axis equal; axis off; drawnow; diff --git a/datasets/New_Zealand_Marine_Environment_Classification_Feature_Layers-shp/Coastline_epsg4326.cpg b/datasets/New_Zealand_Marine_Environment_Classification_Feature_Layers-shp/Coastline_epsg4326.cpg new file mode 100644 index 00000000..3ad133c0 --- /dev/null +++ b/datasets/New_Zealand_Marine_Environment_Classification_Feature_Layers-shp/Coastline_epsg4326.cpg @@ -0,0 +1 @@ +UTF-8 \ No newline at end of file diff --git a/datasets/New_Zealand_Marine_Environment_Classification_Feature_Layers-shp/Coastline_epsg4326.dbf b/datasets/New_Zealand_Marine_Environment_Classification_Feature_Layers-shp/Coastline_epsg4326.dbf new file mode 100644 index 00000000..665ec149 Binary files /dev/null and b/datasets/New_Zealand_Marine_Environment_Classification_Feature_Layers-shp/Coastline_epsg4326.dbf differ diff --git a/datasets/New_Zealand_Marine_Environment_Classification_Feature_Layers-shp/Coastline_epsg4326.prj b/datasets/New_Zealand_Marine_Environment_Classification_Feature_Layers-shp/Coastline_epsg4326.prj new file mode 100644 index 00000000..f45cbadf --- /dev/null +++ b/datasets/New_Zealand_Marine_Environment_Classification_Feature_Layers-shp/Coastline_epsg4326.prj @@ -0,0 +1 @@ +GEOGCS["GCS_WGS_1984",DATUM["D_WGS_1984",SPHEROID["WGS_1984",6378137.0,298.257223563]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]] \ No newline at end of file diff --git a/datasets/New_Zealand_Marine_Environment_Classification_Feature_Layers-shp/Coastline_epsg4326.shp b/datasets/New_Zealand_Marine_Environment_Classification_Feature_Layers-shp/Coastline_epsg4326.shp new file mode 100644 index 00000000..2eeef954 Binary files /dev/null and b/datasets/New_Zealand_Marine_Environment_Classification_Feature_Layers-shp/Coastline_epsg4326.shp differ diff --git a/datasets/New_Zealand_Marine_Environment_Classification_Feature_Layers-shp/Coastline_epsg4326.shx b/datasets/New_Zealand_Marine_Environment_Classification_Feature_Layers-shp/Coastline_epsg4326.shx new file mode 100644 index 00000000..1668156e Binary files /dev/null and b/datasets/New_Zealand_Marine_Environment_Classification_Feature_Layers-shp/Coastline_epsg4326.shx differ diff --git a/imgs/multiscale_nz_ex_13.png b/imgs/multiscale_nz_ex_13.png new file mode 100644 index 00000000..3daecf25 Binary files /dev/null and b/imgs/multiscale_nz_ex_13.png differ diff --git a/utilities/drawedge2.m b/utilities/drawedge2.m index 834d35c8..ff72d721 100644 --- a/utilities/drawedge2.m +++ b/utilities/drawedge2.m @@ -1,9 +1,9 @@ -function drawedge2(pp,e2) +function drawedge2(pp,e2,color) %DRAWEDGE2 draw EDGE2 elements defined by [PP,E2]. Here, PP %is an array of positions & E2 is an array of edge-indexing. - -ec = [0 0 0 ]; - +if nargin < 3 + color = 'k'; +end %-- draw lines in R^2 xx = [pp(e2(:,1),1), ... pp(e2(:,2),1), ... @@ -14,7 +14,7 @@ function drawedge2(pp,e2) line('xdata',xx(:), ... 'ydata',yy(:), ... - 'color',ec, ... + 'color',color, ... 'linewidth',2.0);