From 1b91d1eb198b2e28e077ee635c8ec788a30e01ff Mon Sep 17 00:00:00 2001 From: William Pringle Date: Wed, 13 Jan 2021 16:42:11 -0600 Subject: [PATCH 1/6] adding mesh-ball files from mesh2d that are sometimes used in the refine2_om routine when meshing with mesh2d (#170) --- utilities/GEOM_UTIL/mesh-ball/cdtbal2.m | 77 +++++++++++ utilities/GEOM_UTIL/mesh-ball/inv_2x2.m | 53 ++++++++ utilities/GEOM_UTIL/mesh-ball/pwrbal2.m | 163 ++++++++++++++++++++++++ utilities/GEOM_UTIL/mesh-ball/tribal2.m | 14 ++ 4 files changed, 307 insertions(+) create mode 100644 utilities/GEOM_UTIL/mesh-ball/cdtbal2.m create mode 100644 utilities/GEOM_UTIL/mesh-ball/inv_2x2.m create mode 100644 utilities/GEOM_UTIL/mesh-ball/pwrbal2.m create mode 100644 utilities/GEOM_UTIL/mesh-ball/tribal2.m diff --git a/utilities/GEOM_UTIL/mesh-ball/cdtbal2.m b/utilities/GEOM_UTIL/mesh-ball/cdtbal2.m new file mode 100644 index 00000000..a83b58e7 --- /dev/null +++ b/utilities/GEOM_UTIL/mesh-ball/cdtbal2.m @@ -0,0 +1,77 @@ +function [cc] = cdtbal2(pp,ee,tt) +%CDTBAL2 compute the modified circumballs associated with a +%constrained 2-simplex Delaunay triangulation in R^2. +% [CC] = CDTBAL2(PP,EE,TT) returns the smallest enclosing +% balls associated with the triangles in [PP,TT], such th- +% at CC = [XC,YC,RC.^2]. Such balls never lie outside the +% boundaries of the associated CDT. See TRICON2 for info- +% mation regarding the edge array EE. + +% Darren Engwirda : 2017 -- +% Email : de2363@columbia.edu +% Last updated : 01/10/2017 + +%---------------------------------------------- basic checks + if (~isnumeric(pp) || ~isnumeric(ee) || ... + ~isnumeric(tt) ) + error('cdtbal2:incorrectInputClass' , ... + 'Incorrect input class.') ; + end + +%---------------------------------------------- basic checks + if (ndims(pp) ~= +2 || ndims(ee) ~= +2 || ... + ndims(tt) ~= +2 ) + error('cdtbal2:incorrectDimensions' , ... + 'Incorrect input dimensions.'); + end + if (size(pp,2)~= +2 || size(ee,2) < +5 || ... + size(tt,2) < +6 ) + error('cdtbal2:incorrectDimensions' , ... + 'Incorrect input dimensions.'); + end + +%----------------------------------------- calc. circumballs + cc = tribal2(pp,tt); + +%------------------------ replace with face-balls if smaller + cc = minfac2(cc,pp,ee,tt,1,2,3) ; + cc = minfac2(cc,pp,ee,tt,2,3,1) ; + cc = minfac2(cc,pp,ee,tt,3,1,2) ; + +end + +function [cc] = minfac2(cc,pp,ee,tt,ni,nj,nk) +%MINFAC2 modify the set of circumballs to constrain centres +%to the boundaries of the CDT. +% [CM] = MINFAC2(CC,PP,EE,TT,NI,NJ,NK) returns the set of +% modified circmballs CM, where any ball CC lying outside +% the boundaries of the CDT [PP,EE,TT] is replaced by the +% edge-centred diametric ball. [NI,NJ] are the local inde- +% xes associated with an edge to test. NK is the local in- +% dex of the opposite vertex. + +%------------------------------------------------ outer edge + EF = ee(tt(:,ni+3),5) > +0 ; + +%------------------------------------------------ edge balls + bc = (pp(tt(EF,ni),:)+pp(tt(EF,nj),:))*.50; + +%------------------------------------------------ edge radii + br = sum((bc(:,1:2)-pp(tt(EF,ni),:)).^2,2)... + + sum((bc(:,1:2)-pp(tt(EF,nj),:)).^2,2); + br = br * +0.5 ; + +%------------------------------------------- enclosing radii + ll = sum((bc(:,1:2)-pp(tt(EF,nk),:)).^2,2); + +%------------------------------------------- replace if min. + bi = br >= ll ... + & br <= cc(EF,3) ; + ei = find(EF) ; + ti = ei (bi) ; + +%------------------------------------------- replace is min. + cc(ti,1:2) = bc(bi,:) ; + cc(ti, 3) = br(bi,:) ; + +end diff --git a/utilities/GEOM_UTIL/mesh-ball/inv_2x2.m b/utilities/GEOM_UTIL/mesh-ball/inv_2x2.m new file mode 100644 index 00000000..4b11df19 --- /dev/null +++ b/utilities/GEOM_UTIL/mesh-ball/inv_2x2.m @@ -0,0 +1,53 @@ +function [II,DA] = inv_2x2(AA) +%INV_2X2 calc. the inverses for a block of 2-by-2 matrices. +% [IA,DA] = INV_2X2(AA) returns a set of 'inverses' IA and +% an array of determinants DA for the set of 2-by-2 linear +% systems in AA. SIZE(AA), SIZE(IA) = [2,2,N], where N is +% the number of linear systems. DA is an N-by-1 array of +% determinant values. Note that each IA(:,:,K) is an 'inc- +% omplete inverse DET(A(:,:,K)) * A(:,:,K)^(-1) to improve +% numerical robustness. To solve a linear system, A*X = B, +% compute (I*B)./D, given D is non-zero. +% +% See also INV_3X3 + +% Darren Engwirda : 2018 -- +% Email : de2363@columbia.edu +% Last updated : 13/02/2020 + +%---------------------------------------------- basic checks + if ( ~isnumeric(AA)) + error('inv_2x2:incorrectInputClass' , ... + 'Incorrect input class.') ; + end + +%---------------------------------------------- basic checks + if (ndims(AA) > +3 ) + error('inv_2x2:incorrectDimensions' , ... + 'Incorrect input dimensions.'); + end + if (size(AA,1)~= +2 || ... + size(AA,2)~= +2 ) + error('inv_2x2:incorrectDimensions' , ... + 'Incorrect input dimensions.'); + end + +%---------------------------------------------- build inv(A) + II = zeros(size (AA)) ; + + DA = det_2x2(AA) ; + + II(1,1,:) = AA(2,2,:) ; + II(2,2,:) = AA(1,1,:) ; + II(1,2,:) =-AA(1,2,:) ; + II(2,1,:) =-AA(2,1,:) ; + +end + +function [DA] = det_2x2(AA) + + DA = ... + AA(1,1,:).* AA(2,2,:) - ... + AA(1,2,:).* AA(2,1,:) ; + +end diff --git a/utilities/GEOM_UTIL/mesh-ball/pwrbal2.m b/utilities/GEOM_UTIL/mesh-ball/pwrbal2.m new file mode 100644 index 00000000..a40bc640 --- /dev/null +++ b/utilities/GEOM_UTIL/mesh-ball/pwrbal2.m @@ -0,0 +1,163 @@ +function [bb] = pwrbal2(pp,pw,tt) +%PWRBAL2 compute the ortho-balls associated with a 2-simplex +%triangulation embedded in R^2 or R^3. +% [BB] = PWRBAL2(PP,PW,TT) returns the set of power balls +% associated with the triangles in [PP,TT], such that BB = +% [XC,YC,RC.^2]. PW is a vector of vertex weights. + +% Darren Engwirda : 2017 -- +% Email : de2363@columbia.edu +% Last updated : 02/05/2018 + +%---------------------------------------------- basic checks + if ( ~isnumeric(pp) || ... + ~isnumeric(pw) || ... + ~isnumeric(tt) ) + error('pwrbal2:incorrectInputClass' , ... + 'Incorrect input class.'); + end + +%---------------------------------------------- basic checks + if (ndims(pp) ~= +2 || ... + ndims(pw) ~= +2 || ... + ndims(tt) ~= +2 ) + error('pwrbal2:incorrectDimensions' , ... + 'Incorrect input dimensions.'); + end + + if (size(pp,2) < +2 || ... + size(pp,1)~= size(pw,1) || ... + size(tt,2) < +3 ) + error('pwrbal2:incorrectDimensions' , ... + 'Incorrect input dimensions.'); + end + + switch (size(pp,2)) + case +2 + %-------------------------------------------- alloc work + AA = zeros(2,2,size(tt,1)) ; + Rv = zeros(2,1,size(tt,1)) ; + bb = zeros(size(tt,1),3,1) ; + + %-------------------------------------------- lhs matrix + ab = pp(tt(:,2),:)-pp(tt(:,1),:); + ac = pp(tt(:,3),:)-pp(tt(:,1),:); + + AA(1,1,:) = ab(:,1) * +2.0 ; + AA(1,2,:) = ab(:,2) * +2.0 ; + AA(2,1,:) = ac(:,1) * +2.0 ; + AA(2,2,:) = ac(:,2) * +2.0 ; + + %-------------------------------------------- rhs vector + Rv(1,1,:) = sum(ab.*ab,2) ... + - ( pw(tt(:,2)) - pw(tt(:,1)) ) ; + + Rv(2,1,:) = sum(ac.*ac,2) ... + - ( pw(tt(:,3)) - pw(tt(:,1)) ) ; + + %-------------------------------------------- solve sys. + [II,dd] = inv_2x2(AA) ; + + bb(:,1) = ( ... + II(1,1,:).*Rv(1,1,:) + ... + II(1,2,:).*Rv(2,1,:) ) ... + ./ dd ; + + bb(:,2) = ( ... + II(2,1,:).*Rv(1,1,:) + ... + II(2,2,:).*Rv(2,1,:) ) ... + ./ dd ; + + bb(:,1:2) = ... + pp(tt(:,1),:) + bb(:,1:2) ; + + %-------------------------------------------- mean radii + r1 = sum( ... + (bb(:,1:2)-pp(tt(:,1),:)).^2,2) ; + r2 = sum( ... + (bb(:,1:2)-pp(tt(:,2),:)).^2,2) ; + r3 = sum( ... + (bb(:,1:2)-pp(tt(:,3),:)).^2,2) ; + + r1 = r1 - pw(tt(:,1)); + r2 = r2 - pw(tt(:,2)); + r3 = r3 - pw(tt(:,3)); + + bb(:,3) = ( r1+r2+r3 ) / +3.0 ; + + case +3 + %-------------------------------------------- alloc work + AA = zeros(3,3,size(tt,1)) ; + Rv = zeros(3,1,size(tt,1)) ; + bb = zeros(size(tt,1),4,1) ; + + %-------------------------------------------- lhs matrix + ab = pp(tt(:,2),:)-pp(tt(:,1),:); + ac = pp(tt(:,3),:)-pp(tt(:,1),:); + + AA(1,1,:) = ab(:,1) * +2.0 ; + AA(1,2,:) = ab(:,2) * +2.0 ; + AA(1,3,:) = ab(:,3) * +2.0 ; + AA(2,1,:) = ac(:,1) * +2.0 ; + AA(2,2,:) = ac(:,2) * +2.0 ; + AA(2,3,:) = ac(:,3) * +2.0 ; + + nv = cross(ab,ac) ; + + AA(3,1,:) = nv(:,1) * +1.0 ; + AA(3,2,:) = nv(:,2) * +1.0 ; + AA(3,3,:) = nv(:,3) * +1.0 ; + + %-------------------------------------------- rhs vector + Rv(1,1,:) = sum(ab.*ab,2) ... + - ( pw(tt(:,2)) - pw(tt(:,1)) ) ; + + Rv(2,1,:) = sum(ac.*ac,2) ... + - ( pw(tt(:,3)) - pw(tt(:,1)) ) ; + + %-------------------------------------------- solve sys. + [II,dd] = inv_3x3(AA) ; + + bb(:,1) = ( ... + II(1,1,:).*Rv(1,1,:) + ... + II(1,2,:).*Rv(2,1,:) + ... + II(1,3,:).*Rv(3,1,:) ) ... + ./ dd ; + + bb(:,2) = ( ... + II(2,1,:).*Rv(1,1,:) + ... + II(2,2,:).*Rv(2,1,:) + ... + II(2,3,:).*Rv(3,1,:) ) ... + ./ dd ; + + bb(:,3) = ( ... + II(3,1,:).*Rv(1,1,:) + ... + II(3,2,:).*Rv(2,1,:) + ... + II(3,3,:).*Rv(3,1,:) ) ... + ./ dd ; + + bb(:,1:3) = ... + pp(tt(:,1),:) + bb(:,1:3) ; + + %-------------------------------------------- mean radii + r1 = sum( ... + (bb(:,1:3)-pp(tt(:,1),:)).^2,2) ; + r2 = sum( ... + (bb(:,1:3)-pp(tt(:,2),:)).^2,2) ; + r3 = sum( ... + (bb(:,1:3)-pp(tt(:,3),:)).^2,2) ; + + r1 = r1 - pw(tt(:,1)); + r2 = r2 - pw(tt(:,2)); + r3 = r3 - pw(tt(:,3)); + + bb(:,4) = ( r1+r2+r3 ) / +3.0 ; + + otherwise + + error('pwrbal2:unsupportedDimension' , ... + 'Dimension not supported.') ; + + end + +end diff --git a/utilities/GEOM_UTIL/mesh-ball/tribal2.m b/utilities/GEOM_UTIL/mesh-ball/tribal2.m new file mode 100644 index 00000000..95ce1546 --- /dev/null +++ b/utilities/GEOM_UTIL/mesh-ball/tribal2.m @@ -0,0 +1,14 @@ +function [bb] = tribal2(pp,tt) +%TRIBAL2 compute the circumballs associated with a 2-simplex +%triangulation embedded in R^2 or R^3. +% [BB] = TRIBAL2(PP,TT) returns the circumscribing balls +% associated with the triangles in [PP,TT], such that BB = +% [XC,YC,RC.^2]. + +% Darren Engwirda : 2017 -- +% Email : de2363@columbia.edu +% Last updated : 02/05/2018 + + bb = pwrbal2(pp,zeros(size(pp,1),1),tt) ; + +end From 539a81c9cbfc7373e7d4139b7ead45e5e9900620 Mon Sep 17 00:00:00 2001 From: William Pringle Date: Thu, 21 Jan 2021 13:31:52 -0600 Subject: [PATCH 2/6] adding pfilt option into getnan2 to remove small polygons less than X edges (3 by default for triangle), made small correction in refine2_om func call to make clearer --- utilities/GEOM_UTIL/geom-util/getnan2.m | 17 ++++++++++++++--- utilities/refine2_om.m | 2 +- 2 files changed, 15 insertions(+), 4 deletions(-) diff --git a/utilities/GEOM_UTIL/geom-util/getnan2.m b/utilities/GEOM_UTIL/geom-util/getnan2.m index 9889d41f..e8afd0b4 100644 --- a/utilities/GEOM_UTIL/geom-util/getnan2.m +++ b/utilities/GEOM_UTIL/geom-util/getnan2.m @@ -1,6 +1,6 @@ function [node,edge] = getnan2(varargin) %GETNAN2 parse a NaN delimited polygon into a PSLG. -% [NODE,EDGE] = GETNAN2(NANS,FILT) converts a set of NaN +% [NODE,EDGE] = GETNAN2(NANS,FILT,PFILT) converts a set of NaN % delimited polygons to a PSLG representation. NANS is an % D-by-2 array of coordinates, with polygon vertices spec- % ified in connsecutive order, and delimited by NaN values @@ -13,6 +13,10 @@ % NODE(EDGE(JJ,1),:) and NODE(EDGE(JJ,2),:) are the coord- % inates of the endpoints of the JJ-TH edge. % +% WP: Added PFILT option that specifies minimum allowable +% length of edges for a polyon (set to 3 by +% default for triangle) + % See also FIXGEO2, BFSGEO2, REFINE2 %----------------------------------------------------------- @@ -21,10 +25,11 @@ % Last updated : 06/10/2017 %----------------------------------------------------------- - data = [] ; filt = +0. ; + data = [] ; filt = +0. ; pfilt = +3; if (nargin>=+1), data = varargin{1}; end if (nargin>=+2), filt = varargin{2}; end + if (nargin>=+3), pfilt = varargin{3}; end %---------------------------------------------- basic checks if ( ~isnumeric(data) || ... @@ -45,6 +50,12 @@ 'Incorrect input dimensions.'); end + if pfilt < +1 + warning('getnan2:inputValue' , ... + 'PFILT needs to be at least 1 (setting to 1)') + pfilt = +1; + end + %---------------------------------- parse NaN delimited data nvec = find(isnan(data(:,1))) ; @@ -72,7 +83,7 @@ pdel = pmax-pmin; - if (~isempty(pnew)) + if (size(pnew,1) > pfilt) if (any(pdel>filt)) %---------------------------------- push polygon onto output nnew = size(pnew,1); diff --git a/utilities/refine2_om.m b/utilities/refine2_om.m index bd91c08b..dea30cc7 100644 --- a/utilities/refine2_om.m +++ b/utilities/refine2_om.m @@ -194,7 +194,7 @@ end %-------------------------------- prune any non-unique topo. - [ivec,ivec,jvec] = ... + [~,ivec,jvec] = ... unique(sort(PSLG,+2),'rows') ; PSLG = PSLG(ivec,:) ; From ee895c43882e948274ba417ebaf74ad9d601f3ca Mon Sep 17 00:00:00 2001 From: Keith Roberts Date: Thu, 21 Jan 2021 18:17:57 -0300 Subject: [PATCH 3/6] Update README.md --- README.md | 3 +++ 1 file changed, 3 insertions(+) diff --git a/README.md b/README.md index 3d1c5d9c..a95d650f 100644 --- a/README.md +++ b/README.md @@ -135,6 +135,9 @@ The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.0.0/) ### Unreleased +## Added +- 'mesh2d' interface improvements to filter small polygons. + ## Changed - `msh.plot()` overhaul. All options specified via kwarg. From d3f899ec048290868b296f9db3c5e5fd3badbf88 Mon Sep 17 00:00:00 2001 From: William Pringle Date: Thu, 21 Jan 2021 22:54:03 -0600 Subject: [PATCH 4/6] - Adding '../' to path for the Examples to run correctly from Examples directory - The minimum element quality is now just > 0.5 for NZ example so updating TestSanity check --- Examples/Example_10_Multiscale_Smoother.m | 6 ++++++ Examples/Example_11_Remeshing_Patches.m | 2 ++ Examples/Example_1_NZ.m | 1 + Examples/Example_2_NY.m | 1 + Examples/Example_3_ECGC.m | 1 + Examples/Example_4_PRVI.m | 1 + Examples/Example_5_JBAY.m | 1 + Examples/Example_5b_JBAY_w_weirs.m | 1 + Examples/Example_6_GBAY.m | 1 + Examples/Example_6b_GBAY_w_floodplain.m | 1 + Examples/Example_7_Global.m | 1 + Examples/Example_8_AK.m | 1 + Examples/Example_9_TBAY.m | 5 +++++ Tests/TestSanity.m | 14 +++++++------- 14 files changed, 30 insertions(+), 7 deletions(-) diff --git a/Examples/Example_10_Multiscale_Smoother.m b/Examples/Example_10_Multiscale_Smoother.m index 6dbff6c3..217b02a5 100644 --- a/Examples/Example_10_Multiscale_Smoother.m +++ b/Examples/Example_10_Multiscale_Smoother.m @@ -1,6 +1,12 @@ % Example_10_Multiscale_Smoother: % An idealized test for multiscale nesting using boxes with a large min_el % ratio +clearvars; clc; + +addpath('..') +addpath(genpath('../utilities/')) +addpath(genpath('../datasets/')) +addpath(genpath('../m_map/')) bbox = [0, 1; 0,1]; boubox = [0,0; 1,0; 1,1; 0,1; 0,0; NaN NaN ]; diff --git a/Examples/Example_11_Remeshing_Patches.m b/Examples/Example_11_Remeshing_Patches.m index 81235cfc..21917f63 100644 --- a/Examples/Example_11_Remeshing_Patches.m +++ b/Examples/Example_11_Remeshing_Patches.m @@ -9,9 +9,11 @@ % % By Keith Roberts, 2020, USP, Brazil. clearvars; clc; +addpath('..') addpath(genpath('../utilities/')) addpath(genpath('../datasets/')) addpath(genpath('../m_map/')) + %% STEP 1: Here we mimic the steps that occurred in Example_1_NZ.m bbox = [166 176; % lon_min lon_max -48 -40]; % lat_min lat_max diff --git a/Examples/Example_1_NZ.m b/Examples/Example_1_NZ.m index 52ebedac..5f808cbd 100644 --- a/Examples/Example_1_NZ.m +++ b/Examples/Example_1_NZ.m @@ -1,6 +1,7 @@ % Example_1_NZ: Mesh the South Island of New Zealand clearvars; clc; +addpath('..') addpath(genpath('../utilities/')) addpath(genpath('../datasets/')) addpath(genpath('../m_map/')) diff --git a/Examples/Example_2_NY.m b/Examples/Example_2_NY.m index 515b775c..f8f7fc65 100644 --- a/Examples/Example_2_NY.m +++ b/Examples/Example_2_NY.m @@ -1,6 +1,7 @@ % Example_2_NY: Mesh the New York region in high resolution clearvars; clc; +addpath('..') addpath(genpath('../utilities/')) addpath(genpath('../datasets/')) addpath(genpath('../m_map/')) diff --git a/Examples/Example_3_ECGC.m b/Examples/Example_3_ECGC.m index 6cd05c3e..c410ce27 100644 --- a/Examples/Example_3_ECGC.m +++ b/Examples/Example_3_ECGC.m @@ -3,6 +3,7 @@ clearvars; clc; +addpath('..') addpath(genpath('../utilities/')); addpath(genpath('../datasets/')); addpath(genpath('../m_map/')); diff --git a/Examples/Example_4_PRVI.m b/Examples/Example_4_PRVI.m index 83c22618..dbf66b0d 100644 --- a/Examples/Example_4_PRVI.m +++ b/Examples/Example_4_PRVI.m @@ -4,6 +4,7 @@ % around 1 hr. clc; clearvars +addpath('..') addpath(genpath('../utilities/')) addpath(genpath('../datasets/')) addpath(genpath('../m_map/')) diff --git a/Examples/Example_5_JBAY.m b/Examples/Example_5_JBAY.m index a17fbc31..df57cb6d 100644 --- a/Examples/Example_5_JBAY.m +++ b/Examples/Example_5_JBAY.m @@ -3,6 +3,7 @@ clc; clearvars +addpath('..') addpath(genpath('../utilities/')); addpath(genpath('../datasets/')); addpath(genpath('../m_map/')); diff --git a/Examples/Example_5b_JBAY_w_weirs.m b/Examples/Example_5b_JBAY_w_weirs.m index e40ef094..356e644c 100644 --- a/Examples/Example_5b_JBAY_w_weirs.m +++ b/Examples/Example_5b_JBAY_w_weirs.m @@ -3,6 +3,7 @@ % high resolution with two 15-30 m wide weirs at the mouth of the estuary. clc; clearvars +addpath('..') addpath(genpath('../utilities/')); addpath(genpath('../datasets/')); addpath(genpath('../m_map/')); diff --git a/Examples/Example_6_GBAY.m b/Examples/Example_6_GBAY.m index 5f21b710..4401eff3 100644 --- a/Examples/Example_6_GBAY.m +++ b/Examples/Example_6_GBAY.m @@ -3,6 +3,7 @@ clearvars; clc; +addpath('..') addpath(genpath('../utilities/')) addpath(genpath('../datasets/')) addpath(genpath('../m_map/')) diff --git a/Examples/Example_6b_GBAY_w_floodplain.m b/Examples/Example_6b_GBAY_w_floodplain.m index 401eb70b..ffeadbfc 100644 --- a/Examples/Example_6b_GBAY_w_floodplain.m +++ b/Examples/Example_6b_GBAY_w_floodplain.m @@ -4,6 +4,7 @@ clearvars; clc; +addpath('..') addpath(genpath('../utilities/')) addpath(genpath('../datasets/')) addpath(genpath('../m_map/')) diff --git a/Examples/Example_7_Global.m b/Examples/Example_7_Global.m index f7fd958f..15b8c64a 100644 --- a/Examples/Example_7_Global.m +++ b/Examples/Example_7_Global.m @@ -2,6 +2,7 @@ clearvars; clc; +addpath('..') addpath(genpath('../utilities/')); addpath(genpath('../datasets/')); addpath(genpath('../m_map/')); diff --git a/Examples/Example_8_AK.m b/Examples/Example_8_AK.m index 3fb49a64..10f2fd70 100644 --- a/Examples/Example_8_AK.m +++ b/Examples/Example_8_AK.m @@ -9,6 +9,7 @@ clearvars; clc; +addpath('..') addpath(genpath('../utilities/')) addpath(genpath('../datasets/')) addpath(genpath('../m_map/')) diff --git a/Examples/Example_9_TBAY.m b/Examples/Example_9_TBAY.m index 006f1a70..1b28df87 100644 --- a/Examples/Example_9_TBAY.m +++ b/Examples/Example_9_TBAY.m @@ -22,6 +22,11 @@ %% Setup clearvars; clc; close all; +addpath('..') +addpath(genpath('../utilities/')) +addpath(genpath('../datasets/')) +addpath(genpath('../m_map/')) + %% Define data sources % Shoreline shoreline = 'GSHHS_f_L1'; diff --git a/Tests/TestSanity.m b/Tests/TestSanity.m index f3c3bed0..aac669fe 100644 --- a/Tests/TestSanity.m +++ b/Tests/TestSanity.m @@ -4,26 +4,26 @@ NP_TOL = 500; NT_TOL = 1500; -QUAL_TOL = 0.15; +QUAL_TOL = 0.5; % p: [5968×2 double] % t: [9530×3 double] % Element qual. metrics % MEAN LOWER 3RD MIN. % 0.9325 0.7452 0.3548 if abs(length(m.p) - 5968) > NP_TOL - error(['Incorrect number of points for Example1_NZ.m. Got ',... + error(['Incorrect number of points for Example_1_NZ.m. Got ',... num2str(length(m.p)),' expecting 5968 +- 500 points']); exit(1) end if abs(length(m.t) - 9530) > NT_TOL - error(['Incorrect number of triangles for Example1_NZ.m. Got ',... + error(['Incorrect number of triangles for Example_1_NZ.m. Got ',... num2str(length(m.t)),' expecting 9530 +- 1500 elements']); exit(1) end -if abs(mshopts.qual(end,3) - 0.3548) > QUAL_TOL - error(['Incorrect min. element quality for Example1_NZ.m. Got ',... - num2str(mshopts.qual(end,3)),' expecting 0.3548 +- 0.05']); +if abs(mshopts.qual(end,3)) < QUAL_TOL + error(['Incorrect min. element quality for Example_1_NZ.m. Got ',... + num2str(mshopts.qual(end,3)),' expecting > 0.5']); exit(1) end -disp('Passed: Example1_NZ.m'); +disp('Passed: Example_1_NZ.m'); From 3e19bccb72a257bcf344decd3ea1cc49d17d9f4a Mon Sep 17 00:00:00 2001 From: Keith Roberts Date: Sun, 24 Jan 2021 20:50:29 -0300 Subject: [PATCH 5/6] Update msh.m (#174) * Fix userval and defval order * Sometimes order in defval and userval may not be in same order leading to plotting difficulties. --- @msh/msh.m | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/@msh/msh.m b/@msh/msh.m index e4cb9ee7..6295090f 100644 --- a/@msh/msh.m +++ b/@msh/msh.m @@ -632,6 +632,7 @@ function write(obj,fname,type) disp(['Trying to plot arbitrary f13 attribute: ' type]) if ~isempty(obj.f13) ii = find(contains({obj.f13.defval.Atr(:).AttrName},type)); + jj = find(contains({obj.f13.userval.Atr(:).AttrName},type)); if isempty(ii) error(['no f13 attribute of ' type ' found']); elseif length(ii) > 1 @@ -639,7 +640,7 @@ function write(obj,fname,type) error(['more than one f13 attribute matched ' type]); end defval = obj.f13.defval.Atr(ii).Val; - userval = obj.f13.userval.Atr(ii).Val; + userval = obj.f13.userval.Atr(jj).Val; defval = reshape(defval,1,[]); values = obj.p(:,1)*0 + defval; values(userval(1,:),:) = userval(2:end,:)'; From 4adc12918d8d6c15ac75b31530006384502065e3 Mon Sep 17 00:00:00 2001 From: Qiu Jiangchao Date: Tue, 26 Jan 2021 01:44:18 +0800 Subject: [PATCH 6/6] Make f20 file with a volume flow input (#169) * Support for the creation of aperiodic boundary forcing files (fort.20) * Example 12 to demonstrate how they're created. Co-authored-by: Keith Roberts --- @msh/msh.m | 3 + @msh/private/writefort20.m | 26 +++++++ Examples/Example_12_Riverine_flow.m | 107 ++++++++++++++++++++++++++++ README.md | 3 +- Tests/test_make_f20.csv | 48 +++++++++++++ utilities/Make_f20_volume_flow.m | 71 ++++++++++++++++++ utilities/Riverflux_distribution.m | 98 +++++++++++++++++++++++++ 7 files changed, 355 insertions(+), 1 deletion(-) create mode 100644 @msh/private/writefort20.m create mode 100644 Examples/Example_12_Riverine_flow.m create mode 100644 Tests/test_make_f20.csv create mode 100644 utilities/Make_f20_volume_flow.m create mode 100644 utilities/Riverflux_distribution.m diff --git a/@msh/msh.m b/@msh/msh.m index 6295090f..342a66cc 100644 --- a/@msh/msh.m +++ b/@msh/msh.m @@ -257,6 +257,9 @@ function write(obj,fname,type) if any(contains(type,'19')) && ~isempty(obj.f19) writefort19( obj.f19, [fname '.19'] ); end + if any(contains(type,'20')) && ~isempty(obj.f20) + writefort20( obj.f20, [fname '.20'] ); + end if any(contains(type,'2001')) && ~isempty(obj.f2001) writefort19( obj.f2001, [fname '.2001'] ); end diff --git a/@msh/private/writefort20.m b/@msh/private/writefort20.m new file mode 100644 index 00000000..5348f4df --- /dev/null +++ b/@msh/private/writefort20.m @@ -0,0 +1,26 @@ +function [] = writefort20( fort20dat, finame ) +% +% + +if ( nargin == 1 ) + finputname = 'fort.20_1' ; +else + finputname = finame ; +end + +fid = fopen(finputname,'w') ; + +fprintf( fid, '%d\n', fort20dat.FTIMINC ) ; + +valpernode = size( fort20dat.Val,2 ) ; +str = []; +for ll = 1: valpernode + str = [str '%.8f '] ; +end +str = [str '\n' ] ; + +fprintf( fid, str, fort20dat.Val' ) ; + +fclose(fid) ; +%EOF +end diff --git a/Examples/Example_12_Riverine_flow.m b/Examples/Example_12_Riverine_flow.m new file mode 100644 index 00000000..079802db --- /dev/null +++ b/Examples/Example_12_Riverine_flow.m @@ -0,0 +1,107 @@ +% Example_12_Riverine_flow: Guidance for creating a riverine flow input file +% otherwise referred to as a fort.20 +% +% You should prepare a column delimited file (csv) that stores the total volume +% flow (Q) (m^3/s) time series of a cross-section where the riverine boundaries +% in the mesh are located. +% +% A command-line method or a data cursor method is available to identify the +% vstart and vend of each riverine boundary in msh.make_bc. +% +% The columns of the csv file should be organized as follows: +% year,month,day,hour,minute,second,volume_flow_1,volume_flow_2... +% The column order of the total volume flow (volume_flow_1,volume_flow_2...) +% must be specified in the order in which the riverine boundaries appear in +% the fort.14 file, or in which you make the riverine boundaries with the data +% cursor method or comannd-line method in msh.make_bc. +% +% A 'test_make_f20.csv' file has been placed in the Tests folder for reference. +% The volume flow of this file is an hourly average time series, thus DT equals +% to 3600 (s), if your volume flow is a daily average time series, DT should be +% set to 86400 (s) rather than 3600 (s). +% +% More details can be found in the two functions of Make_f20_volume_flow.m and +% Riverflux_distribution.m. +% +% Author: Jiangchao Qiu +% Created: January 7 2021 +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +addpath('..') +addpath(genpath('../utilities/')); +addpath(genpath('../datasets/')); +addpath(genpath('../m_map/')); +addpath(genpath('../Tests/')); + +%% STEP 1: set mesh extents and set parameters for mesh. +%% The Pearl River Estuary region +bbox = [112 20; 112 24; 116 24; 116 20; 112 20]; %polygon boubox +min_el = 100; % minimum resolution in meters. +max_el = 10e3; % maximum resolution in meters. +dt = 5; % Automatically set timestep based on nearshore res +grade = 0.3; % mesh grade in decimal percent. +R = 3; % Number of elements to resolve feature. + +%% STEP 2: specify geographical datasets and process the geographical data +%% to be used later with other OceanMesh classes... +dem = 'SRTM15+V2.1.nc'; +coastline = 'GSHHS_f_L1'; +gdat = geodata('shp',coastline,'dem',dem,'h0',min_el,'bbox',bbox); +%plot(gdat,'shp'); +%plot(gdat,'dem'); + +%% STEP 3: create an edge function class +fh = edgefx('geodata',gdat,'fs',R,'max_el',max_el,'dt',dt,'g',grade); + +%% STEP 4: Pass your edgefx class object along with some meshing options and... +%% build the mesh +mshopts = meshgen('ef',fh,'bou',gdat,'plot_on',1,'proj','lam'); +mshopts = mshopts.build; + +%% STEP 5: Plot and save the msh class object +m = mshopts.grd; % get out the msh object +plot(m,'type','tri'); + +m = interp(m,gdat,'mindepth',1); % interpolate bathy to the mesh with minimum depth of 1 m +m = lim_bathy_slope(m,0.1,0); +CFL = CalcCFL(m,dt); +max(CFL) +min(CFL) +m = BoundCr(m,dt); + +%% STEP 6: Make the nodestring boundary conditions +m = make_bc(m,'auto',gdat); + +%% STEP 7: Identify vstart and vend of each riverine boundary +% use the comannd-line method to identify the vstart and vend for each... +% riverine boundary based on the coordinates of river edge. +river_points1 = [112.918 22.508;112.919 22.514]; % coordinates of river edge +river_points2 = [113.283 22.365;113.287 22.367]; % coordinates of river edge +river_points3 = [113.513 22.234;113.515 22.236]; % coordinates of river edge +bc_k1 = ourKNNsearch(m.p',river_points1',1); +bc_k2 = ourKNNsearch(m.p',river_points2',1); +bc_k3 = ourKNNsearch(m.p',river_points3',1); +m = make_bc(m,'outer',0,bc_k1(1),bc_k1(2),1,22); +m = make_bc(m,'outer',1,bc_k2(1),bc_k2(2),1,22); +m = make_bc(m,'outer',1,bc_k3(1),bc_k3(2),1,22); + +% use the data cursor method to identify the vstart and vend for each... +% riverine boundary manually (GUI). +% m = make_bc(m,'outer',1); % 1,2206,2107,1(flux),22(River) +% m = make_bc(m,'outer',1); % 1,14784,14779,1(flux),22(River) +% m = make_bc(m,'outer',1); % 1,763,704,1(flux),22(River) + +plot(m,'type','bd'); % plot mesh on native projection with boundary conditions +plot(m,'type','b'); % plot bathy on native projection + +%% STEP 8: Make the riverine conditions +ts = '01-Jan-2005 00:00'; % start time of the total volume flow time series +te = '01-Jan-2005 23:00'; % end time of the total volume flow time series +DT = 3600; % time step of the total volume flow time series, + % DT=3600 for hourly data series while DT=86400 for daily data series. +m = Make_f20_volume_flow(m,'test_make_f20.csv',ts,te,DT); + +%% STEP 9: Save data +save('Riverine.mat','m'); write(m,'Riverine','20'); + + + diff --git a/README.md b/README.md index a95d650f..e8ae9788 100644 --- a/README.md +++ b/README.md @@ -136,7 +136,8 @@ The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.0.0/) ### Unreleased ## Added -- 'mesh2d' interface improvements to filter small polygons. +- `mesh2d` interface improvements to filter small polygons. +- Support for creation of `fort.20` files for forcing rivers by @Jiangchao3 ## Changed - `msh.plot()` overhaul. All options specified via kwarg. diff --git a/Tests/test_make_f20.csv b/Tests/test_make_f20.csv new file mode 100644 index 00000000..f184294f --- /dev/null +++ b/Tests/test_make_f20.csv @@ -0,0 +1,48 @@ +2005,1,19,0,0,0,1780,1960,534 +,,19,1,0,0,1960,1960,734 +,,19,2,0,0,1960,1960,794 +,,19,3,0,0,1610,1960,872 +,,19,4,0,0,549,1960,802 +,,19,5,0,0,-641,1960,714 +,,19,6,0,0,-1230,1960,401 +,,19,7,0,0,-1010,1960,-216 +,,19,8,0,0,-157,1960,-573 +,,19,9,0,0,350,1960,-531 +,,19,10,0,0,741,1960,-356 +,,19,11,0,0,912,1960,38.1 +,,19,12,0,0,934,1960,350 +,,19,13,0,0,752,1960,492 +,,19,14,0,0,199,1960,581 +,,19,15,0,0,-563,1960,571 +,,19,16,0,0,-1500,1960,401 +,,19,17,0,0,-1960,1960,-81.3 +,,19,18,0,0,-2110,1960,-667 +,,19,19,0,0,-1880,1960,-944 +,,19,20,0,0,-1280,1960,-999 +,,19,21,0,0,-279,1960,-995 +,,19,22,0,0,106,1960,-792 +,,19,23,0,0,1310,1960,-355 +2005,1,20,0,0,0,1910,1470,231 +,,20,1,0,0,1940,1470,528 +,,20,2,0,0,2100,1470,671 +,,20,3,0,0,1860,1470,706 +,,20,4,0,0,1580,1470,689 +,,20,5,0,0,1210,1470,696 +,,20,6,0,0,87.8,1470,702 +,,20,7,0,0,-638,1470,664 +,,20,8,0,0,-1200,1470,230 +,,20,9,0,0,-1370,1470,-261 +,,20,10,0,0,-1210,1470,-608 +,,20,11,0,0,-792,1470,-722 +,,20,12,0,0,-307,1470,-790 +,,20,13,0,0,96.9,1470,-645 +,,20,14,0,0,356,1470,-74.5 +,,20,15,0,0,98.3,1470,302 +,,20,16,0,0,-710,1470,523 +,,20,17,0,0,-1400,1470,347 +,,20,18,0,0,-2150,1470,-381 +,,20,19,0,0,-1960,1470,-865 +,,20,20,0,0,-1580,1470,-1040 +,,20,21,0,0,-1180,1470,-1020 +,,20,22,0,0,-632,1470,-896 +,,20,23,0,0,112,1470,-708 diff --git a/utilities/Make_f20_volume_flow.m b/utilities/Make_f20_volume_flow.m new file mode 100644 index 00000000..05da8a52 --- /dev/null +++ b/utilities/Make_f20_volume_flow.m @@ -0,0 +1,71 @@ +function obj = Make_f20_volume_flow(obj,filename,ts,te,DT) +% obj = Make_f20_volume_flow(obj,filename,ts,te,DT) +% Input a msh class object to get the values of the nodal fluxes (q) (m^2/s) for +% times between ts and te based on the input file. filename is the name of a column +% delimited file (csv) that stores the total volume flow (Q) (m^3/s) time series of +% the cross-sections where the riverine boundaries in the mesh are located. +% +% ts and te represent the start and end time of the total volume flow, respectively. +% +% The columns of the csv file should be organized as follows: +% year,month,day,hour,minute,second,volume_flow_1,volume_flow_2... +% The column order of the total volume flow (volume_flow_1,volume_flow_2...) +% must be specified in the order in which the riverine boundaries appear in the +% fort.14 file, or in which you make the riverine boundaries with the data cursor +% method in msh.make_bc. +% +% A description of transforming the volume flux (m^3/s) to areal flux (m^2/s): +% +% The function of Riverflux_distribution (has been placed in the utilities folder) +% will be called first to calculate the representative edge width and flux percentage +% for each node on the riverine boundary. Nodal flux (m^3/s) can be calculated by +% multiplying the total volume flow (Q) (m^3/s) by the percentage of flow dedicated to +% that boundary node. Then the nodal flux (q) (m^2/s) can be calculate by dividing +% nodal flux by the nodal edge width. +% +% Note: the default total volume flow (Q) (m^3/s) is an hourly average time series, +% thus the DT equals to 3600 (s). If your Q is an daily average time series, you should +% use a value of 86400(24*3600) (s) for DT. +% +% Author: Jiangchao Qiu +% Created: January 7 2021 +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +[edgesize,pertg,nvell_num] = Riverflux_distribution(obj); + +total_node_num = sum(nvell_num); +N = length(nvell_num); + +data_riverflow = csvread(filename); + +river_volumeflow = data_riverflow(:,7:7+N-1); + +if DT==3600 % the total volume flow is an hourly average time series: + time = data_riverflow(:,4); +elseif DT==86400 % the total volume flow is a daily average time series: + time = data_riverflow(:,3); +else + error('The total volume flow should be an hourly or daily average time series') +end + +T = length(time); +F = zeros(total_node_num,T); + +for i=1:T + flow_this_time = river_volumeflow(i,:); + flow_this_time1 = repmat(flow_this_time,[size(pertg,1),1]); + flow_Q = (flow_this_time1.*pertg); + flow_q = (flow_this_time1.*pertg)./edgesize; + flow_q1 = flow_q(:); + %for the condition that river boundaries have different number of nodes + flow_q2= flow_q1(~isnan(flow_q1)); + F(:,i) = flow_q2(:); +end + +%% Make into f20 struct +obj.f20.DataTitle = [datestr(ts) ' -> ' datestr(te)]; +obj.f20.FTIMINC = DT; +obj.f20.NumOfNodes = size(F,1); +obj.f20.NumOfTimes = size(F,2); +obj.f20.Val = F(:); +%EOF diff --git a/utilities/Riverflux_distribution.m b/utilities/Riverflux_distribution.m new file mode 100644 index 00000000..ca5e56f8 --- /dev/null +++ b/utilities/Riverflux_distribution.m @@ -0,0 +1,98 @@ +function [result_edgesize,result_pertg,nvell_num] = Riverflux_distribution(obj) +% [result_edgesize,result_pertgļ¼Œnvell_num] = Riverflux_distribution(obj) +% +% Input a msh class object to get the representative edge width and flux +% percentage for each node on the riverine flow boundary. +% +% This function is used to distribute a total flux of a cross-section where +% the riverine boundary located to each node of the riverine boundary. +% +% result_edgesize and result_pertgis are the result of the representative +% edge width and flux percentage for each node, respectively. nvell_num +% represents the number of nodes for each riverine boundary. +% +% Each column of result_edgesize and result_pertg represents the result of +% each riverine boundary. +% +% Nodal representative edge width equals to the sum of half the width of +% each of the two edges it is connected to. +% +% Nodeal flow percentage on the riverine boundary equals to the flow area +% of this node divided by the total flow area of the cross-section. +% +% Nodal flow area is calculated via a trapezoidal rule using its representative +% edge width and the bathymetric depths of this node and its neighboring nodes. +% +% Note that the calculation of representative edge width and flow percentage +% for the nodes at either end of the boundary is a little bit of different. +% +% User need to specify river flow boundaries (ibtype=22) before using it. +% +% Author: Jiangchao Qiu +% Created: January 7 2021 +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +bd_dat = obj.bd; +p_dat = obj.p; +b_dat = obj.b; +river_num = find(bd_dat.ibtype == 22); +N = length(river_num);% total number of river boundarys +nvell_num = bd_dat.nvell(river_num);% number of nodes for each river boundary + +% Consider the riverine boundaries may have different number of nodes, the +% number of rows for result_edgesize and result_pertg is set to max(nvell_num) +result_edgesize = NaN(max(nvell_num),N); % the final result of edgesize +result_pertg = NaN(max(nvell_num),N);% the final result of percentage + +if isempty(river_num) + error('No riverine boundary information to distribute total flow') +end + +for i=1:N + node_num = bd_dat.nbvv(:,river_num(i)); + node_num(node_num==0)=[]; + bathy = b_dat(node_num); + location = p_dat(node_num,:); + J = length(node_num);% total number of nodes for the current river boundary + flowarea = zeros(J,1); + edgesize = zeros(J,1); + %% calculate the projected distance for each edge on the boundary + proj = 'Mercator'; + m_proj(proj,'lon',[ min(obj.p(:,1))-0.25 max(obj.p(:,1))+0.25 ],... + 'lat',[ min(obj.p(:,2))-0.25 max(obj.p(:,2))+0.25]) + proj_location = zeros(J,2); + for j=1:J + [proj_location(j,1),proj_location(j,2)] = m_ll2xy(location(j,1),location(j,2)); + end + node_distance = 1000*m_xydist(proj_location(:,1),proj_location(:,2)); + %% calculate the respresentive width for each node + % note: if the current boundary has J nodes, there are J-1 edges + for j=1:J + if j==1 % the first edge (the number is J) + edgesize(j) = 0.5*node_distance(j); + elseif j==J % the last edge (the number is J-1) + edgesize(j) = 0.5*node_distance(j-1); + else % the other edges + edgesize(j) = 0.5*node_distance(j-1)+0.5*node_distance(j); + end + result_edgesize(1:J,i) = edgesize; + end + %% calculate flow area for each node on the boundary + for j=1:J + if j==1 % the first node + local_z1 = (bathy(j)+bathy(j+1))*0.5; + flowarea(j) = 0.5*(bathy(j)+local_z1)*(0.5*node_distance(j)); + elseif j==J % the last node + local_z1 = (bathy(j-1)+bathy(j))*0.5; + flowarea(j) = 0.5*(bathy(j)+local_z1)*(0.5*node_distance(j-1)); + else % the other nodes + local_z1 = (bathy(j-1)+bathy(j))*0.5; + local_z2 = (bathy(j)+bathy(j+1))*0.5; + flowarea(j) = 0.5*(bathy(j)+local_z1)*(0.5*node_distance(j-1))+... + 0.5*(bathy(j)+local_z2)*(0.5*node_distance(j)); + end + percent = flowarea/sum(flowarea); % the flux distribution percentage for each node + result_pertg(1:J,i) = percent; + end +end + +