Skip to content

Commit

Permalink
Merge pull request #2 from CHLNDDEV/Projection
Browse files Browse the repository at this point in the history 
merge new version
  • Loading branch information
@Jiangchao3
Jiangchao3 authored Oct 21, 2020
2 parents 18a5074 + 2f757b1 commit e4a9342
Show file tree
Hide file tree
Showing 11 changed files with 267 additions and 58 deletions.
7 changes: 5 additions & 2 deletions @geodata/geodata.m
View file
Original file line number Diff line number Diff line change
Expand Up @@ -57,6 +57,7 @@
pslg % piecewise liner straight line graph
spacing = 2.0 ; %Relative spacing along polygon, large effect on computational efficiency of signed distance.
gridspace
shapefile_3d % if the shapefile has a height attribute
end

methods
Expand All @@ -75,7 +76,6 @@
%addOptional(p,'weirs',defval);
%addOptional(p,'pslg',defval);
%addOptional(p,'boubox',defval);
%addOptional(p,'window',defval);

% Check for m_map dir
M_MAP_EXISTS=0 ;
Expand Down Expand Up @@ -119,6 +119,7 @@
addOptional(p,'pslg',defval);
addOptional(p,'boubox',defval);
addOptional(p,'window',defval);
addOptional(p,'shapefile_3d',defval);

% parse the inputs
parse(p,varargin{:});
Expand Down Expand Up @@ -202,6 +203,8 @@
% Default value
obj.window = 5;
end
case('shapefile_3d')
obj.shapefile_3d = inp.(fields{i}) ;
case('weirs')
if ~iscell(inp.(fields{i})) && ~isstruct(inp.(fields{i})) && inp.(fields{i})==0, continue; end
if ~iscell(inp.(fields{i})) && ~isstruct(inp.(fields{i}))
Expand Down Expand Up @@ -321,7 +324,7 @@
end

polygon_struct = Read_shapefile( obj.contourfile, [], ...
obj.bbox, obj.gridspace, obj.boubox, 0 );
obj.bbox, obj.gridspace, obj.boubox, 0, obj.shapefile_3d);

% Unpack data from function Read_Shapefile()s
obj.outer = polygon_struct.outer;
Expand Down
5 changes: 4 additions & 1 deletion @geodata/private/Read_shapefile.m
View file
Original file line number Diff line number Diff line change
@@ -1,5 +1,5 @@
function polygon_struct = Read_shapefile( finputname, polygon, bbox, ...
h0, boubox, plot_on )
h0, boubox, plot_on, shapefile_3d)
% Read_shapefile: Reads a shapefile or a NaN-delimited vector
% containing polygons and/or segments in the the desired region
% of interest. Classifies the vector data as either a
Expand Down Expand Up @@ -174,6 +174,9 @@
% using m_shaperead
points = tmpC{i,1}(1:end,:) ;
if size(points,2) == 3
points = points(:,1:2);
end
if shapefile_3d
% if 3-D shapefile
height = points(:,3);
points = points(:,1:2);
Expand Down
103 changes: 56 additions & 47 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,25 +54,25 @@
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
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
end


Expand Down Expand Up @@ -98,11 +122,10 @@
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);

addOptional(p,'enforceMin',1);

% parse the inputs
parse(p,varargin{:});
Expand All @@ -114,13 +137,14 @@
% kjr...order these argument so they are processed in a predictable
% manner. Process the general opts first, then the OceanMesh
% classes...then basic non-critical options.
inp = orderfields(inp,{'h0','bbox','enforceWeirs','fh','inner','outer','mainland',...
inp = orderfields(inp,{'h0','bbox','enforceWeirs','enforceMin','fh',...
'inner','outer','mainland',...
'bou','ef',... %<--OceanMesh classes come after
'egfix','pfix','fixboxes',...
'plot_on','nscreen','itmax',...
'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 @@ -183,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 @@ -193,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 @@ -290,6 +313,11 @@

% 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

Expand All @@ -304,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 @@ -381,6 +407,8 @@
del = setProj(dmy,1,obj.proj) ;
case('enforceWeirs')
obj.enforceWeirs = inp.(fields{i});
case('enforceMin')
obj.enforceMin = inp.(fields{i});
end
end

Expand Down Expand Up @@ -446,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 @@ -537,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 @@ -678,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 @@ -743,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 @@ -944,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
13 changes: 13 additions & 0 deletions @meshgen/private/edgefx2pts.m
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,13 @@
function [pts, vals] = edgefx2pts(fh)
% takes an edgefunction from OceanMesh2D and outputs a series of points and
% values for use in fastscatter
% Coleman Blakely
% January 22, 2020
x = fh.F.GridVectors{1};
y = fh.F.GridVectors{2}';
val = fh.F.Values;
vals = reshape(val,[],1);
[lat, lon] = meshgrid(x,y);
c = cat(2,lat',lon');
pts = reshape(c,[],2);
end
56 changes: 56 additions & 0 deletions @meshgen/private/enforce_min_ef.m
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,56 @@
function fh = enforce_min_ef(fh)
% Ensures that the minimum mesh sizing value is always used in the case
% of many edgefxs (e.g., multiscale meshing technique).
%
% Loop through edgefunctions to check if they overlap.
% If they overlap, find the points that lie within the
% upper edgefunction and compare them to the ones below.
%
% By Coleman Blakely UND, 2020.
% Modifications by Keith Roberts USP, 2020.

for kk = length(fh):-1:2
fh_test = fh{kk};
% make the points to use with inpoly to find the overlapping points
[test_pts, ~] = edgefx2pts(fh_test);
test_vals = fh_test.F.Values;
x_test = fh_test.F.GridVectors{1};
y_test = fh_test.F.GridVectors{2};
[X_test, Y_test] = ndgrid(x_test,y_test);
% Loop through all boxes that may be below the test box
for ii = kk-1:-1:1
fh_lower = fh{ii};
% find the test points that are within the lower boundary
id_overlap = inpoly(test_pts,fh_lower.boubox(1:end-1,:));
% check to see if the test box overlaps the lower box. If not, skip
% this box
if sum(id_overlap)==0
continue
end
overlap_pts = test_pts(id_overlap,:);
% create the grid of overlapping points to use in the gridded
% interpolant
x = unique(overlap_pts(:,1));
y = unique(overlap_pts(:,2));
[X, Y] = ndgrid(x,y);
% evaluate the gridded interpolant of the lower box using the
% points that are within the test box
lower_vals = fh_lower.F(X,Y);
% find the test points within the test box and compare with the
% newly found points
[~,idx,~] = intersect(x_test,x);
[~,idy,~] = intersect(y_test,y);
% int_vals is a dummy variable to make the indexing of the boxes
% the same for comparison
int_vals = test_vals(idx,idy);
% replace all test values that are greater than the lower values
int_vals(int_vals>lower_vals) = lower_vals(int_vals>lower_vals);
% place the int_vals back into the test values
test_vals(idx,idy) = int_vals;
end
% put the new points back in the edgefunction
fh{kk}.F = griddedInterpolant(X_test,Y_test,test_vals);



end
Loading

0 comments on commit e4a9342

Please sign in to comment.