Skip to content

Commit ac4c6a2

Browse files
WPringleWPringle
authored
Merge branch 'Projection' into sal_netcdf
2 parents cf1a84f + d6d17a1 commit ac4c6a2

31 files changed

+403
-181
lines changed

‎.gitignore‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -3,7 +3,7 @@ datasets/
33
fort.*
44
Examples/*
55
!Examples/Example_*
6-
/*.mat
6+
*/*.mat
77
*.15
88
*.14
99
*.13

‎@edgefx/edgefx.m‎

Lines changed: 54 additions & 44 deletions
Original file line numberDiff line numberDiff line change
@@ -5,25 +5,37 @@
55
% resolution when building a mesh.
66
% Copyright (C) 2018 Keith Roberts & William Pringle
77
%
8-
% The following properties (with default values) are available for input
9-
% to the edgefx method:
10-
% defval = 0; % placeholder value if arg is not passed.
11-
% addOptional(p,'dis',defval);
12-
% addOptional(p,'fs',defval);
13-
% addOptional(p,'wl',defval);
14-
% addOptional(p,'slp',defval);
15-
% addOptional(p,'ch',defval);
16-
% addOptional(p,'min_el_ch',100);
17-
% addOptional(p,'AngOfRe',60);
18-
% addOptional(p,'max_el',inf);
19-
% addOptional(p,'max_el_ns',inf);
20-
% addOptional(p,'g',0.20);
21-
% addOptional(p,'geodata',defval)
22-
% addOptional(p,'lmsl',defval)
23-
% addOptional(p,'dt',-1);
24-
% addOptional(p,'fl',defval);
25-
% addOptional(p,'Channels',defval);
26-
% addOptional(p,'h0',defval);
8+
% The following properties are available for input to the edgefx method:
9+
% Keyword Default Val Variable type
10+
% 'dis' defval; scalar
11+
% 'fs' defval; scalar
12+
% 'wl' defval; scalar/array
13+
% 'slp' defval; scalar/array
14+
% 'ch' defval; scalar
15+
% 'min_el_ch' 100; scalar
16+
% 'AngOfRe' 60; scalar
17+
% 'max_el' inf; scalar/array
18+
% 'max_el_ns' inf; scalar
19+
% 'g' 0.20; scalar/array
20+
% 'geodata' defval; geodata class
21+
% 'lmsl' defval; geodata class
22+
% 'dt' -1; scalar
23+
% 'fl' defval; scalar
24+
% 'Channels' defval; cell array
25+
% 'h0' defval; scalar
26+
% [defval = 0; % placeholder value if arg is not passed.]
27+
%
28+
% The 'wl','slp','max_el','g' options inputs can be scalars to apply
29+
% the parameter globally or they can be arrays to bound each parameter
30+
% value within topographic elevation bounds, using the following stucture:
31+
% [parameter_value1 z_min1 z_max1;
32+
% parameter_value2 z_min2 z_max2;
33+
% ...
34+
% parameter_valuen z_minn z_maxn];
35+
%
36+
% where, z_min is the minimum elevation bound and z_max is the maximum elevation bound.
37+
% - similar to axis bounds on matlab plots such as caxis, xlim, ylim, ...
38+
% - z is oriented same as the DEM, i.e., +ve value of z is above datum (overland)
2739
%
2840
% This program is free software: you can redistribute it and/or modify
2941
% it under the terms of the GNU General Public License as published by
@@ -393,7 +405,7 @@
393405
d = reshape(d,obj.nx,[]);
394406
end
395407

396-
%% Wavelength edgefx.
408+
%% Wavelength edge function.
397409
function obj = wlfx(obj,feat)
398410

399411
% interpolate DEM's bathy linearly onto our edgefunction grid.
@@ -413,25 +425,23 @@
413425
if numel(param)==1
414426
% no bounds specified.
415427
wlp = param(1);
416-
% set cuttof at 10 m by default
417-
dp1 = -10;
418-
dp2 = -inf;
428+
z_min = -inf;
429+
z_max = +inf;
419430
else
420431
wlp = param(1);
421-
dp1 = param(2);
422-
dp2 = param(3);
432+
z_min = param(2);
433+
z_max = param(3);
423434
end
424-
% limit to 1 m
435+
% limit min. depth to 1 m
425436
twld = period*sqrt(grav*max(abs(tmpz),1))/wlp;
426-
% Set wld with mask applied
427-
obj.wld(tmpz < dp1 & tmpz > dp2 ) = ...
428-
twld(tmpz < dp1 & tmpz > dp2);
437+
limidx = (tmpz >= z_min & tmpz < z_max);
438+
% Set wld with depth mask applied
439+
obj.wld(limidx) = twld(limidx);
429440
clearvars twld
430441
end
431442
clearvars tmpz xg yg;
432443
end
433444

434-
%% Topographic length scale/slope edge function.
435445
%% Topographic length scale/slope edge function.
436446
function obj = slpfx(obj,feat)
437447

@@ -444,7 +454,7 @@
444454
clear tmpz2
445455
end
446456
tmpz(tmpz > 50) = 50; % ensure no larger than 50 m above land
447-
% use a harvestine assumption
457+
% use a Harvestine assumption
448458
dx = obj.h0*cosd(min(yg(1,:),85)); % for gradient function
449459
dy = obj.h0; % for gradient function
450460
% lets filter the bathy to get only relevant features
@@ -613,19 +623,19 @@
613623
if numel(param)==1
614624
% no bounds specified. valid in this range.
615625
slpp = param(1);
616-
% default cutoff is 10 m
617-
dp1 = -10;
618-
dp2 = -inf;
626+
z_min = -inf;
627+
z_max = +inf;
619628
else
620629
slpp = param(1);
621-
dp1 = param(2);
622-
dp2 = param(3);
630+
z_min = param(2);
631+
z_max = param(3);
623632
end
624633
% Calculating the slope function
625634
dp = max(1,-tmpz);
626635
tslpd = (2*pi/slpp)*dp./(bs+eps);
627-
obj.slpd(tmpz < dp1 & tmpz > dp2 ) = ...
628-
tslpd(tmpz < dp1 & tmpz > dp2);
636+
% apply slope with mask
637+
limidx = (tmpz >= z_min & tmpz < z_max);
638+
obj.slpd(limidx) = tslpd(limidx);
629639
clearvars tslpd
630640
end
631641
clearvars tmpz xg yg
@@ -850,10 +860,10 @@ function plot(obj,type)
850860
hh_m( limidx ) = mx;
851861
else
852862
mx = param(1);
853-
dp1 = param(2);
854-
dp2 = param(3);
863+
z_min = param(2);
864+
z_max = param(3);
855865

856-
limidx = (tmpz < dp1 & tmpz > dp2) & ...
866+
limidx = (tmpz >= z_min & tmpz < z_max) & ...
857867
(hh_m > mx | isnan(hh_m));
858868

859869
hh_m( limidx ) = mx;
@@ -893,10 +903,10 @@ function plot(obj,type)
893903
dmy = dmy + lim ;
894904
else
895905
lim = param(1);
896-
dp1 = param(2);
897-
dp2 = param(3);
906+
z_min = param(2);
907+
z_max = param(3);
898908

899-
limidx = (tmpz < dp1 & tmpz > dp2) ;
909+
limidx = (tmpz >= z_min & tmpz < z_max) ;
900910

901911
dmy( limidx ) = lim;
902912
end

‎@meshgen/meshgen.m‎

Lines changed: 27 additions & 11 deletions
Original file line numberDiff line numberDiff line change
@@ -28,7 +28,7 @@
2828
% itmax % maximum number of iterations.
2929
% pfix % fixed node positions (nfix x 2 )
3030
% egfix % edge constraints
31-
% outer % meshing boundary (manual specification, no bou)
31+
% outer % meshing boundary (manual specification, no bou)
3232
% inner % island boundaries (manual specification, no bou)
3333
% mainland % the shoreline boundary (manual specification, no bou)
3434
% fixboxes % a flag that indicates which boxes will use fixed constraints
@@ -564,15 +564,6 @@
564564
p1 = p1(rand(size(p1,1),1) < r0/max_r0,:); % Rejection method
565565
p = [p; p1]; % Adding p1 to p
566566
end
567-
% add new points along boundary of multiscale nests
568-
if box_num < length(obj.h0)
569-
h0_rat = ceil(h0_l/obj.h0(box_num+1));
570-
nsplits = floor(log(h0_rat)/log(2));
571-
for add = 1:nsplits
572-
new_points = split_bars(p,fh_l);
573-
p = [p; new_points];
574-
end
575-
end
576567
end
577568
else
578569
disp('User-supplied initial points!');
@@ -656,7 +647,31 @@
656647
disp('Exiting')
657648
return
658649
end
659-
N = size(p,1);
650+
651+
% Getting element quality and check "goodness"
652+
if exist('pt','var'); clear pt; end
653+
[pt(:,1),pt(:,2)] = m_ll2xy(p(:,1),p(:,2));
654+
tq = gettrimeshquan( pt, t);
655+
mq_m = mean(tq.qm);
656+
mq_l = min(tq.qm);
657+
mq_s = std(tq.qm);
658+
mq_l3sig = mq_m - 3*mq_s;
659+
obj.qual(it,:) = [mq_m,mq_l3sig,mq_l];
660+
661+
662+
% If mesh quality went down "significantly" since last iteration
663+
% which was a mesh improvement iteration, then rewind.
664+
if ~mod(it,imp+1) && obj.qual(it,1) - obj.qual(it-1,1) < -0.10 || ...
665+
~mod(it,imp+1) && (N - length(p_before_improve))/length(p_before_improve) < -0.10
666+
disp('Mesh improvement was unsuccessful...rewinding...');
667+
p = p_before_improve;
668+
N = size(p,1); % Number of points changed
669+
pold = inf;
670+
it = it + 1;
671+
continue
672+
else
673+
N = size(p,1); pold = p; % Assign number of points and save current positions
674+
end
660675
% 4. Describe each bar by a unique pair of nodes.
661676
bars = [t(:,[1,2]); t(:,[1,3]); t(:,[2,3])]; % Interior bars duplicated
662677
bars = unique(sort(bars,2),'rows'); % Bars as node pairs
@@ -779,6 +794,7 @@
779794

780795
% Mesh improvements (deleting and addition)
781796
if ~mod(it,imp)
797+
p_before_improve = p;
782798
nn = []; pst = [];
783799
if qual_diff < imp*obj.qual_tol && qual_diff > 0
784800
% Remove elements with small connectivity

0 commit comments

Comments
 (0)