From e078ab465249f5e99c57e45b42ac0cfd1d9dbdcc Mon Sep 17 00:00:00 2001 From: BillSenior <89970704+BillSenior@users.noreply.github.com> Date: Thu, 30 May 2024 10:12:48 +0200 Subject: [PATCH] Casulli derefinement (#329|gridedit 778) --- .../CasulliRefinement/casulli_deref_21_21.nc | Bin 0 -> 8824 bytes .../CasulliRefinement/casulli_deref_21_22.nc | Bin 0 -> 8392 bytes .../CasulliRefinement/casulli_deref_22_21.nc | Bin 0 -> 9820 bytes .../CasulliRefinement/casulli_deref_22_22.nc | Bin 0 -> 9436 bytes .../casulli_derefine_polygon.nc | Bin 0 -> 25324 bytes .../casulli_dref_two_polygon.nc | Bin 0 -> 23644 bytes .../casulli_polygon_then_all.nc | Bin 0 -> 12316 bytes libs/MeshKernel/CMakeLists.txt | 4 + .../MeshKernel/CasulliDeRefinement.hpp | 191 ++++ .../include/MeshKernel/Entities.hpp | 20 + libs/MeshKernel/include/MeshKernel/Mesh.hpp | 27 +- .../include/MeshKernel/Operations.hpp | 4 +- .../include/MeshKernel/Polygons.hpp | 5 + .../UndoActions/FullUnstructuredGridUndo.hpp | 68 ++ .../MeshKernel/Utilities/NumericFunctions.hpp | 33 + libs/MeshKernel/src/CasulliDeRefinement.cpp | 1013 +++++++++++++++++ libs/MeshKernel/src/Mesh.cpp | 22 + libs/MeshKernel/src/MeshRefinement.cpp | 3 +- libs/MeshKernel/src/Operations.cpp | 25 +- libs/MeshKernel/src/Polygons.cpp | 27 + .../UndoActions/FullUnstructuredGridUndo.cpp | 20 + .../tests/src/MeshRefinementTests.cpp | 368 ++++++ .../include/MeshKernelApi/MeshKernel.hpp | 30 + libs/MeshKernelApi/src/MeshKernel.cpp | 152 +++ libs/MeshKernelApi/tests/src/ApiTest.cpp | 281 +++++ tools/test_utils/CMakeLists.txt | 2 + .../include/TestUtils/MakeMeshes.hpp | 25 +- .../include/TestUtils/MeshReaders.hpp | 38 + tools/test_utils/src/MakeMeshes.cpp | 32 +- tools/test_utils/src/MeshReaders.cpp | 51 + 30 files changed, 2423 insertions(+), 18 deletions(-) create mode 100644 data/test/data/CasulliRefinement/casulli_deref_21_21.nc create mode 100644 data/test/data/CasulliRefinement/casulli_deref_21_22.nc create mode 100644 data/test/data/CasulliRefinement/casulli_deref_22_21.nc create mode 100644 data/test/data/CasulliRefinement/casulli_deref_22_22.nc create mode 100644 data/test/data/CasulliRefinement/casulli_derefine_polygon.nc create mode 100644 data/test/data/CasulliRefinement/casulli_dref_two_polygon.nc create mode 100644 data/test/data/CasulliRefinement/casulli_polygon_then_all.nc create mode 100644 libs/MeshKernel/include/MeshKernel/CasulliDeRefinement.hpp create mode 100644 libs/MeshKernel/include/MeshKernel/UndoActions/FullUnstructuredGridUndo.hpp create mode 100644 libs/MeshKernel/src/CasulliDeRefinement.cpp create mode 100644 libs/MeshKernel/src/UndoActions/FullUnstructuredGridUndo.cpp create mode 100644 tools/test_utils/include/TestUtils/MeshReaders.hpp create mode 100644 tools/test_utils/src/MeshReaders.cpp diff --git a/data/test/data/CasulliRefinement/casulli_deref_21_21.nc b/data/test/data/CasulliRefinement/casulli_deref_21_21.nc new file mode 100644 index 0000000000000000000000000000000000000000..7ce38031ba3687e758064fa8bf0535032bf1f213 GIT binary patch literal 8824 zcmeI0TX)k|8itkBIS~Q`Xn+Q&dWI4Qp9o2MKxv%N6hZ?eoQtB^vLj+kMwa7*c3PK~ zmX@++u4ZoMV%GeQT+HSCgt_3K%)IZ`r;r3Pht3Lo*YaKCWjM)$xLb}ot_@h$7)4Ak=9d*RC+*9o*hl4 zMw6KVy^L!O^;*ZGQw=c^x`*%fu6B;)&teWO7t>Kb1%&5*vnA^V*$9 z4>K2blzATBt~1W}sP>qIQ*i9OQmA$MXMMWIJ$Ko5xvcB=plj~@m8onzF*p)BVzrt> z9q-SHGlk@EI-SYIQ@QlncqWq>iI1d)lkwa_+8VNkbGei?9Cp*XksY9Y+O969vhrzo zPuqDeS1^6Aq;3M;J!P~r1FzzhyrNn;P21;UP0gI=`FY31a&VGtdq7IdT+yRdGH zKJS!mm$wcr!}c3%vCqyIZBxXqYv%%I*$Kh{>F$uSWZyli?KZx!mbLxSl+DfUSmWFi zgzrHYYq=|$@d^&`iuMg_{T0wwbDeNUgmZ7fX5g!*8dz@L^7E!^mBZ3iYjLRJd-r%1 zo@=gbSS$AFDtT^kOZ(j`@#hPQHnvNO!-^sbH{yTG9ed*kK0?e&nuG7wM z5Si80YQV=u_)^z^w=r2F-p}4ScggkYuD;$sY#XyH)gs%1=sJ6gzLPi0R;A*+ux{IV z`nkquntO z#Wv^ua`5OM_M**x2V2*#+Ge3-1%d5i#;UJ!8UCM+{XQr^#wT49Ggr@J$<%`w$6S+k zem0s=?c0{EQf)&$^{Zpusp|IqnpM&RPiJdAGDt1={*CUamLDk9`gQU?aPMzd%bgXk z>R{h-Sk9-lR>8m9*ErV9!gb;~VCL;gfO5syFRAVL0-muYd)Y4OR-n5!z;##k6X;vl zZG5rdVTUX6U98+mTJYWSd^fB=1Gjs=#Q2OtAGL-Xo@<>YFK5Y@6VJssEK>GX?=!sL z_*`Kd)&w)!_uK974qQe_eQmz{nvd|LzF-G++jjM!?(w?SWkZfvH=Y}QU^;HzR!^L$ z-Z?(rcj}tf_bJ=qUhYF}IMW5Pb zR7M_;j${^}v5g*WTYRyUHhhdd6*xJoR9ZE$aJi*TeS8mFyja6_^PhJK=NsFnX~!2E zoA)|S{qMev!EkOnmT!wMBRSSkZ;YMFkHwpnydeIa$|><)lDx9UH!Vh(+ON*b4f-@= zR)uC(EUz^}-%S;oac3<&-l@FZ2z_@|XvX}r@Oa142z^!~%-sCW=E7R&%d5}~ z`&oFrV`+rGLL*E^g=Q@NUA=yF7S}d;Lp$ne>%QB%PqxyLD%>)@_!7osD~M?t7OyPTpt#O z@qb!BYE$s~!~NCrFZNf9@gMH52EGH=emks&X=ZaBSP0EJXx71hqYmKYZ_Yl@D#_o2 zhd@lS4S5*wHzSk35s!dQ$u8t;fPaT(^6!xSfWKpzJ>V$VBgwy|90U0Mrm`1#9PF1o zfP5VslzbJb3vVDf*Ea#@_!i)t_&+kv=^em1ybCz?dw^q~037oq;8^bij&TamIUvN5 zbn}-Q-KPQF{L2m9NkDfB(47W!X8_%20Nq1??z4bS7BAs92Hm58?hgRn9|F2R0(A3# zlysj1bQ^%~F+lexfbMZXXI7X%(me_2J`d>T-y7(j0(5@{=)M39a1qcw4d~`CGP*wp zbj}Euk#t`HbYBH@^A{G~bAaw^fbQ#n?i+yan}F_HfbQFX&O5>vNV@L=x=ld$JfPbG zbT0t9bAawVpxXv?7XaNwK8JW{S~15H-PT10o~sKIvc{bNV*>ay1xT-e-G$>1nB+& z(ETHz`!S&VCqVblfbJ)N&R>M5$QIBlyn-Y#p$$peg$^X?6uOXPm(YzQyM-Pk*(2;l zl6^uilI#}_Ajv`DRpjeH2Yuj>a2QFB2(KYYzi<>ujtR$+gf8+1I0Xj4o8T?sZ6tX| zco#|D6HXw>N#T7Yp^L=9SuhL+!D%6ZBuOEKBxxaoBxi&nB%zCpfN_uoqu>MKLnQe~ z_!vpf2?mmk37;SdU1S1$4rah4I4^vPBvZm?NOD2Ah$Pd(B_yGXTn4wn9dHF)6=soS zPPm37*M%EMa#Og4By^E4KoKm0yTBCYk;D=fkR&JMk;E1XNJ1BJKotbw9#|4eNKzJD zB=Lj_lH3=3B%zCNFMk6XU>Vef6(r#vCJ%%!k>o3ZdzgGJ&_%ul55afB_eklNplUFHPEa9*1$gi!#JO= literal 0 HcmV?d00001 diff --git a/data/test/data/CasulliRefinement/casulli_deref_21_22.nc b/data/test/data/CasulliRefinement/casulli_deref_21_22.nc new file mode 100644 index 0000000000000000000000000000000000000000..eece9b64a021b45763ade9f559de25fe470cc95a GIT binary patch literal 8392 zcmeI0+j8618HFL4*1@uD%V$S+f+S92$0AKqvQ9Y|iYPgXEl0K`pJNCjHYrFXKm(vC zIVZ(&W-`-@_Nu+@MW)Z;i(d4i50Pi+3$$x({FLR!YBN(_G&91?$Au02d+h`My-9Gg zFx_S($3933botj+d_8a#&Yv~Q*F1j}=l=jbywCDn9HD5m16^mTqN;|`_7lqcc@6Y< zeiVCgJ@x`0#U!?%DzOtPR^At?k_wey)TQM9@#Q$KjbFUjXf!gegl7B-s~a6rP!Egj zSrR+b^cpY=>FG+)FsEnD+bWFEo|((!Gr7@hCO4EH9vqwUea*D|yP*7J8okQ=Oj4z+Y27r*-QYbT3@FIa5eyGh>NI zs$TbK$LDiy%pJPq=0=^glN)x@${ot4oe^g&J?xAOIeC>+rQu6SHKPaG2|8v~w45)< z)#RCW@LsWGhe1U@1bX*Y5ogCiEvN)#J#&U}$j6-8#UKb>&&PDIlZ-tgsb#O5%`;u& zfR1;)s`B~h5Shep&&A`)Eh}4$eP0!0Z_SI72I<`?Y00_=bnFg(lFNsiUAAA2OQAm3 z1mQX8W-fO(d%TVUj^?^yuI~YR>%Nz)h@|grSPVRRquB9XCvy3|2Ub zz_-`8%@yl(R|3EMZ}I!r(=R3zF?bd!6QOK>t8X&@t?b*zuj$*yEBn?vfWGz1IF-B2 z$c{FmnA=5i=rQ0lB~!%b*|*@Y`a#1tx9TUcDOHQgj0Mec_LW1=wX068=KW&ccJk@x z9Aj@nq6Q_q79#oVs)CQ_=mcc0H~rc2!AMm-w6PL|wzKZ(b*Ez{?#}5!IbHZ7-h~$2d|B1 zb^P8nr&8aRPyM=VdXcFf*PV(P2PUJ)Tc$pJD$k(#1 zlT6Py?>E;7<2V<*wihtISI1YC!wlCo62BkWpe^IJa|OMnN9Idc^UbTOip+4ab>s@T z#u&4?Z1KfX*>)LgD)x#_rLtjT;)<&;)d>SE@p2u@EqqZW>F=Ju8O9l7^X0XUTmS!u zl}Y;8a-G=J##XEKs`$2AZ$4*#+pI*57HOMxNBbHr*N5WUYF+#p(_k1@tF^dDTf-_Y zQd+0#=ftno>S3LFiwHTJNl5%M?(em>|o8w-wbb0vizk7P}hpA<2OcxDndQ+cRb5^1Y zqW48@(IwO1VYz=!v}`7t__po(TAYvFdH9XygHC3A$hv6zFYBVs`gQBV*hEe||NW@8 z-__lDU?S}1!EPRG&jXzE-RuNt7v+1?Nstod`_n1FcPlF2rCtMFqTSGb!2fHh{J-ur z;QJo67n}hHMERe?>)?=RAM_1ySo8?g1V=@WLEi+&McLQ40Q>khV4uzc_URqKKD-Oq z_IrSBp95_3JYZYz0~1^T^d128OaprIO@^Kq0X?&To;g6zJfPoq@}cKApyvmGo*x2wegx>rUvTui4CrYAdQJd(egfz@3Fuh>^qd0poCfr~ z0_gcEpyv#r=V!nIR{=e*0ea2?dh(YdJwFHZoCEZ{0qA)X(DN3cCx1!Oa{10Q7ta==lYp=OaK*8_;tR(9;3*Tmtkg0(!cDo(j;j1n5}?^jrq? z^Z-3qz#>=$^sE4SRslVIK+gcsvj*t-7|=5W^o#&KV?fV3pywK(X9Li49nfp`=6Tgpw|y8%p*G zJy5b==!KF4!a*oGB=kYaVc`gr92Jg1$#HN3oD@z$$!kJCl$;jMK*{UE8&F~jZ$jS! zZwqIkjL3<;N@WLOx1j)E~^97;YA zK7^8wgpZ-*vS2~UgzyP;5)_0fD47H;N{prkA;Ly0G>K*_4WHAt!g z*I)o@!eb~2g$PPwp$;W$0@ohl+9Mm_3HVa@3QE2fzJZc&g-s~=PIw9>-wV&6?Yr#G L!EO%xwsYWLH+y%Y literal 0 HcmV?d00001 diff --git a/data/test/data/CasulliRefinement/casulli_deref_22_21.nc b/data/test/data/CasulliRefinement/casulli_deref_22_21.nc new file mode 100644 index 0000000000000000000000000000000000000000..d6d4f0dbbd916dfcea39a67c83b771495bc816a0 GIT binary patch literal 9820 zcmeI1X>$`-7KUZVE5>XlU}8v`tPm_C$-A@Aju+y9i2<8EwA)fkZd>XWt(F(EDI_E$ zBq1wg51~@2`2q6-{9!)MkH|lndEeWojWKv?s^Sk!ja2pU_Hyof&%J%_mD~&r4p&5i z&jv^psB?~4{xLUY;rLE?e$;kmBaw2u$_W$ z`$gY&9h?TSgI3lzJ&TPucvjl-EGMbXg8MtBeLvrK@ZkLXd|OIF+ng+$M`{Xg(MvM3 zAa*45YSjj#!&!G;8y?k8SzZDCY4Ns1Tf8UM7Vk)OwYQJ9YW>BG7HiiMu|#{T7Vqv$ zB>LhVty&Ia4eTSe({{mky~WVKZNRfk-%4q)L$i)(tTP(x8IQ;N;yryGT?b-`SZqbt zDzA+xv@nuEL&>vwt;V?DsN#`?owlu%vQYE%pV~AGxXzsAa9Wq2LH)$x6C;DsSX*!4 zQB^G3wBz&H*A<&ec67$m(e~c1_Gou!M{jg0k?M(>UEK+5s=cdcDxD6RMH-;Bpk~x6 zOeY5AYVb^JcrTeYJU6Qz0=*k^h%wNyVq^oaOM*Au@_8Rr7YQYlZsu#FR10nwlHPm z-nc5X=b9iq2KB7v`Y_`KoZu6#8`kZXmnhvGAG7JO_|6(DOp^k>hT3s%lXA2Y6Jm#hAZ_$VHB z=QRsoEZFQY63@ZR?WYUATkfugNm**g z3G1Ex5qJfD^&fe3;bFsuTMJzY(Zj&J2s=r4fd;eCt zcdg^r!~4VcnG3N$&nzB;uO_@$IY zbnDOHeBD$J{L-aBr;UH{qfbxD`Crt`Qs84NiRd$WaQ={MkS*$d_5QGWJ+EGWtJllw z`o6CoWj@0Cx$k+){qXs$n$NC(Ucb_Ob!b}e4Kmh;OX@%$lXerL=R^(BDJ{Z}p1K`q zMhi52PNHT1^8Bzc^m6?K>tmhQmonx2KIHpM3BM_pZuF>h%UYfX7DBlW%60JnTn89C z-!<$2m7;tPuoF~?RzsTs-xW~#ZeSOv6Rn3n2KfIERK9<10epW!Z3KLexJi`n(Dnel z`&7CW`Z(Ao+5~+9Y!`hLstJ3c%=JmY9G?QriQhcT=^4Nro(1gtIl#X60rt5cu&?I< z`#1pT*(yY#^yXja^gama&G$(3jstop0KMA*y*mKCI|04B0KK~bJ-NsTf8o-*577Gs zK<^g;y0KLto z3XecZl~4^OH9{?v)Cu)avO#Eol8r(mlxz|eM zt57l`yapvlgriV0Djb6ndXeMcJh%W(fRn;FluQVxpyaf021?Ef=b(gM&1P@B+MT+2EPy%ycURZ#V zMd31(yeYf|C2tGwKncCbd*C|w6kGxC3m-tqRpCP@`AE11B_9i)KncCbXW(;iL-+zp zz7)QKlCOntpyXTOCY0O~zJrqQg&&~gN8u+Z`C0e{N`4h?L&d+M9SZ0O literal 0 HcmV?d00001 diff --git a/data/test/data/CasulliRefinement/casulli_deref_22_22.nc b/data/test/data/CasulliRefinement/casulli_deref_22_22.nc new file mode 100644 index 0000000000000000000000000000000000000000..8860af0d5df0bb80f70cb7ec2f07e8e9ea716a73 GIT binary patch literal 9436 zcmeI0%X1S)9>>Ry@B?EWhVY0XnIt43!7`Hkz&H=av25aii2<8;7>}ivG`2J&MkC86 zFNK7Jge2q%38x&k_Oyq+=Efn%J?tN_haC66*!_NIelqq3Z*A48ww5ZX>f?`kru*AH zJ>4Tqqsj3mEqrZ+BB0GVsr-{}TH*b@TKRXy3V34{`C)_Ww(h?L7V<9L4#45sD=28V`@MEm=!R5TOoPerZ6R$sJ#D4vO>6a51!GZZ${ zI^b5&GO0?rL{g51_q2uQQW?{8^TADEc4qx7S8d4V+Irjr+!KWFK|9xSdu_x^DBxAwH(cwFfemHH4tGR2_9kowzB)_3<)key zZ8}yVT)M$r>?(ThJg>ra%%xRx#WwAE*UA1%{Y^{J2NQ}qxEGPXqAL4e8Jo<1J;%2E zpNwt!p=0a(8^+cx*QwI0g3Qu#$>-xDd5xn=y+X7N1#A$g!l zYsn7wT~KG_{PZ6U)>x@yJHM}U^w%AG?`}bzE*o3dD=9ORw|rkY7;(@x$TGY??7QO^ zZsJJ$*wmS0n6ltOj56n>g~x_tLH!MjR=&JyJ#`0V-7e|sYT3%`zN@pgUQ3X*-1R4B z*Vgj2AZq10c^|m<*IUbNMYm*Q-?3Ycr@36j-}Tox=FGx%;yGZZRngC}?!(r0bPmth zyjoOw-STxu0bF-UKY+FsUEzZT4?A3m=V0buWC}i89_)sdXW)9zmk6IxXcNrgs^?l; z-c4EZ;ly(hb_>hi>iRR@uYRtu4%Yd?3+9{WYkP$CI2Sy&0kD2&P`{uouJE8m zn0J#MnsVJ{JQ+OFOZu9liP}+Bz{*JD(GgGL80+Y@ZHo_<{Hnv)Q@)+D^7&;G6PH?8 z)5deL#j|B>H}64{aJ&&;)AU?*?Ow;J|NocK8;-H=9;pVwM!kJpa_jBW_c`7sBeYg0 zHyP)GT&?cDEV=de>~EL`O*884)GT8)BQ;Ca&+4}@uinm?|KQ8W&4&5^wIy%f2=g8h z2ZO~QV}JP0=oh~xPKqbQC&j0#!5d>*%4el~LCTk;d_~H$M)1b4b)kBB5O^sZG z^U;`x#(G$@FM$kBE8tdRuUI!REzpL*8&0>C!-VGvR zexKe0_+6Xkck3rWo46hB2K@b==I^(A0l)v#oq)gFZx-`6s3*Zz@izD=&?Vju>tKiY zarkMlQ_QhG131QK0mrl-a7@nuj^TO0zFz?B`v72{{ON{$y##b{5HPz3Ff$66$#1sI zJOr4@U%Z$Z2h2qF!MTK z=0|{;9|LCI5Pt$Q`&02vn0nAyM;Y4c|z!h$zEX}Or8{;f{89X4L<{(752mAIpKMjydWHa$&12E zFgYmnz){dE9D+$qh{Gfy^ueTG7=X#3a2Or}M}%RRyezx|lUIe;U~*J2U@{`S4v&JQ zFb0!x;TTNb5GG*qrf?i4Cxl7(BzQ}hg2^f2G)&G2(=eG4&cfuJa2~z@E((`ma#?sA zCRYR#CbNPClQ|&;r-2eOFv$uznApNROcsPZObUVnyPzmsg^4GWVB!m9m@Envm@Em) z@HOy`@GeZ=6W)i(2SODl9}3rD@{#Z{d;@$Ud(z;cG*jYfJQ|Dp{0Xa6@nBH8L~k@L)x5{K$?^^rGTq$R1{p`78P;g z#J#HGUbwg7=I(v3ci-=O@=L{g;ZeCI zlrR||T6c8%vHJLnvPVYa)&7y={L|b-CF}X~a`JQX2Ic1D^~)bnP%y5Ke^^zSKexc2 zpPOIM$DcQ_Fh9RAuU{X3IbBl^?`u3aQkje=7N;LSXL!O5CEc)}hvRPjvUB@q=jKnz z%P;I-P?($FJ2yW!cZ;QUH(&EG_Xw7x?#OM;*=bvM=6fx?;~64S8gavXfsgOHPqX3i z*g`j^ZLNC_vZjujG`c7|H)n9_F&b1=Mf7lbte&|A`6Z#;{sq|s!=nfO;Np@2`GZrt%5%o?jnQbFTRAtsh|f|VdQKamaY<<~5s!L@>#FuL;rT$kMr- zm7aJJH~2`uZnQ7;(Xc8ONxdRcvpae*@X?|&8H$BNiEuC$Dz{D6yd5hN@p)ta3K@D>jg ztXy1~^ez`1k-A=FkqvB%v_Ea8#^%T3)iHnF{dFYAxD}OUN~C+XDN97c!SYZ=MdUyB zZCj7iMSIM(i4!$m8mvfg@&(IX9w!p=PC%XO>Br|&U#m*DJi;3D;)!5rQN(-QdB+T! zF=cH)N7Nk?iT&p>LpQYVx~dh&w2ko6O;ozU(r74|bYo2H-Pe1pKw?d|mC5p@d?jnd zgvq1WGF^!Ep*?A&uNCHcw;L`DMXR>>UKj7SKT_#;7gU9!{$$**qy6beuyeV`GVWb> zzAW(`s;;NbCB62ybuKrlh*w5<-Dz3PQ?IImKev63vu7cCavcbU-HPN~9lN!2JG+=` zY}8%oM*X3rKjwnomHwXGx7zRW#=^ypCW#n3*GrGWb?Tc}L)|rSTi2Hcx<+vyZx6S) zt~H6qOG5VMr0+G*vZTGXdaU5{Ew3v&M(2c5d!^U?wqD=F@$) zw%@Ii9PE#mYRUAB%{YH}5hX&gGVgJ8-*jJ{6U{hZ%YxoX=Iyt>pMH(#IPHbL)(|?r zg?D_p8`2r>-6QpQOSNFodI$50yi0ne|65BuYMYdXfCe%bQI65Zm< zyrzUMX5^-MQWmUY~65resF5_UV6}(`;Y1H>caaY;SJwx<}jo?ChLn`?aHTw(Z{!&Fx+PSlhQv9q#dPHMev<&K&o>@DL9>~84>Qv0)~Wob>S_nwwz zwW;2FSbC-Tqwe>}eyOk2{hVH8x6?ltrhnc`|NN8wxhDPopMD=tzptg=pVIFS>FaFz z`jWmbq+e(0^Dg~7L~Q+e0q>K*{F+ooYg6g{wIB1x1F72;>D-jMy`UzQ3BP5nm+3z; zy;D1pdA@BQezfQLYHz^6o^pQgs-80gseKMw7W-4T!>dzSX161DJ8HLMc6(RNy|dPx z8t~~Ou)FNjZnH1mw(k7*u^-F(=XRe@yJ+A3{pX&{K6+0?;J^L(@yD`{-}Q6scfBv{ zu0N;lu5)L8uSmWAGQUry_F*Ue9$CYGgRWgZ$Xhd$nSIcSkU0-B=RxK?$h;0R*MZD+ zAafnaTn94Ofy{Lva~;TB2Qt@ztzQRtZ8-3w9`=L&UHblLfCJDF-JyTGejxO3)ek}w z^gtGRLjONpFX;b=IvD!*^u5sneb5rwXoVcKMqji+F4`gw?U0Z5D8OFO-)we3e{{qE zbizP%hRW@Y!RUfQbj2a)hC{Ir4#U0}0zZaAodZzkVNmD8q0Ymh&P7n?5m4umQ0Gxl z=Odub`Y)L3d?eI)4Agln)Oj4#c|6qlD5&!UsPjaq^CYPAWT^8LsPj~)^E9aQbg1(T zsI&g&RGnu+oo7LvXG5KjfjaA71_Tg(^Saj0_z)Oi8aIRSO9ggPgo&Q(z7g-~auI#)xT7eSpDL!FmE zosWY$9}jgt0dsI7)cGW+a}Ct_WT^8gQ0G&j&Zj}0Plr060d+nT>bw-{d=}LCY^d`& zQ0HY(=X0UX=Ruv9L!DPZozI6lUjTK!5bAgl)bV1d<4UOGDyZWnP{&K5j+a3lFNZo_ zff`&1b-W7dcs0~Ut;C^)9IE-B8!{P}h5)uJ>X!?t{AC z4|RP2TJu3@&4-{hABNU^1X}Y^XwAo3TFFPVwuvboN%Qhr? z=Cmc*Gp8NNo;mGF_RQIfWY3%qBzxv`B-t~k6L>(UGnn1k8_e!>0kb<@!R$^qFuSu4 znBCbI%IvJ;CfwFKC8?p&5EZGxULG z$cARffoA9n&5#StkO$3>56w^j&Cm~;p+7Xk0BDAR&{TW|#%dFdLfT z7-)uLp&5eE40E6vLeLDw&MlQhQlySw?D4=R$iLf%bGBw5Ri-J&jtHlj;n(*6;P)IP^SdcsS@gxggRA0ofbl!s-aGcpiYaSPD`Lp$3dNrhdP}AbvhC1bQ095 z2I_P&)aevx=2M}WPlIMY9h&(JXy!AanU_K{p9Rf)HZ=1&(9FxAna_o0J`b9CIW+SM zXy)^wnJ<85z7U%EB53A|p_x}gGp~YXz66^2Qs}w3%<^(l&&d^*SCV>euClzE)N{1j z@)}Z|uZ22a2X$Tpb-o_zycX(w1JwCOsPj!w=bNF!Hr~K%MV}I^PF%z8~uR0Mz+GsPjWm=ZB%rk3gLtg*rb5b*_avKMr+% z0_yxE)OiEc`6;OL(@^JUpw7=iou7j`KM!?&0qXoB)cGZ-^UF}@SD?-tq0X;DonM1G zzYcYN1M2)H)cGx_^V?A8cc9MiLY?1(I=>Hf{s8LyA=LRJsPo5A=S@)OPoU19LY+T@ zI)4sz{sQX!CDi#VsPkqYf1sYJPf7#Rkd#KIF)2+<7AZ|lGg6wH7NoQ^tw?EY+K|%L zv?Hax*^878rXwkxOlMN|HeE>RYPylKkJ*>>V?VP$DF>MDq#S4tBBh7vNlGtsFe$xF zA5yYS4k>+2E-85?pOgaAkCgsq04W2_AW{aKLQ)Pfhmvxb8A1+4zzieha5J2gA~S-N zk!BPrN0`y19BIaoGS-YEWxP3xlnG`cDU-}(Ql^-xq)ao@Ntt1eCS|6XMapb*40$Yq zW)3MKQ%p*U36tWQQc}vyTv8%t9x3xpl$3H4BPDJsNLgSKq*R(DDOF}6Db;2XDT~b# zQjWv%IKiAq%1Ne%l#|UVq?~F_Bjt2+1}SHnrKFr?&L-s?vy7B;&3U9OH!Dav-&{b- zg(lteB5upYW+f@B%q8TdxXfHm$`$5HQm!&rld{@eL&~-0I`U@RVs0g6ow<#a+i?f( zGihJ?36g?lbq3@_^C0vc_Cb%35;+DL0y%NYSkFpm~Uths`6TJc`FqYaS=% z3G*Z=8_ZLrJZ+vKU%^K6Dk-m-*GYK;Z{jWUHYx9zcS(8AyidvpM(fD4<~dTHH!qO# zqIrpwmyKqW56wrUd~7z6@(Di0XXbNKzA#^s@|D?4%1_Ks$=~3&=69t0-u!`-Kbk+0 z@-_a9znH(0@;CE$QvPB7Nx~t-`~6#M|IGZHlwX)%lJYC_Yf>~Tyx%2xAn$j{dSrdm zfRu)&5h;z)1X-pjDa}lCQd*dnq_i@vNoix+lAX}m>`h7+)0LENW*<`ag&+Hw{Yg2% zbSLFNa}X&#OixmJ8OOnXxHG95_iXf(4NY!{-g{r z14$WV29r`~4k3r(a5J2gA~S-Nk!BPrM_@FLG-F5^YsQf>-W)~B1T&G8Nk((Xq2@4B zhM1wG1dL{u$z}>EQ_VC|rkfe09F3WnWoDCdj5(H+pqWET$P|-OV#1`jrj(R2Gnbr? zs3|8UX5yq&m<6OHOeHBvRAHg1CS{RXOv(~-94W_}6G%DHoJ2~EIhmAGjOLbznMaDc z$f@QuQcgE#kaDJ3O3GR0Y*Nm_GMsD9BW1Z+LCX2&0#YtC7m;$YSxL$&a|tPzn#)MJ z++0D*mF6l^t~RSlxyD>e%5`Q9Dc75|q}*U`B;_V^GkFVcHS0*Z&D>7P9p+9_?lO0i zvfkW7%Dv`3Qtmeokn*5;h?IxTBcwcP9wVjJJWk3J=1EdEn5Rg2+B`$bv*tPSdAwj= zB;_UZGAXZ^jikJ4UL)mo^9Ctznzu-K+q^@{yXHMo-ZvkR@}c>Nl#k6OQa&-ClJc4P zoRlxjm!y1UHv7K$owi=)^Sa;fGkcLaFEZyq<{ZeJ1DSIma}H$Afy_COIR`T5K;|6C KoC80CbKqauX__Gb literal 0 HcmV?d00001 diff --git a/data/test/data/CasulliRefinement/casulli_dref_two_polygon.nc b/data/test/data/CasulliRefinement/casulli_dref_two_polygon.nc new file mode 100644 index 0000000000000000000000000000000000000000..fae5af3f4b1d2ae8e300b66d1717c12419bc7d53 GIT binary patch literal 23644 zcmeI(2Y3{98prV+2qlCLQUr8WP!JMQ=qfA{niT0>mrZsO7B;(Kvmt;jcC26nyQql0 zH|$Zl-h02}dcX5}FW=wHJBXZlq8tb3!{&KD{I<-@|9$`SpVl3e+j2LVjb&Wz+#%yuxfwTA={@t`@3ktMtsOUfctb-&QQS5yN+q>>sBK-kK2u3$ z{Ak?i!Lf#g%L}I@(+%M%v%-tqOdaEeONvU1O2!lym5eAIRaQ1@Sa@Q6b-1`JTv}XO zHY{8+dR%GgxRMdW!Zq~jA3oG>RiZAN&TI&dUo#g?0p&b-&j+nyp4RS7rF3;eXSPr;;gYMqv+gWbis?u_M)jy|g^2~`UC3N?QB_U-_G8)%^c5P%c1T;$XRJ@)i`~=sH z)_QBSsZS;RE5fhc(Z#?|r@Cw`6^~`&(NwI)U%K90^sddMSF4M3D!P83xzaxQ$#klE zzwJA&FMNAK*@kn`In1USe{XFz{|B?S4ezD44ewf8?_t!|&-!%Z`y41*x1lcU-7eU| zUt9TPo=uU?XSao^HK}w%D!g}pZOJxnZC$kzL7&~KGl_V#CRSUU_>XyO?l}E4$K0D( zsOhR`ZHC1ct#LU{BIYeXt@YseJPdWIb88Zeu{xcJ#?~ji>(1L|;*teBBidr`Hi^`K zZZlzb&;HxGuub;_mu{xcjaDUN*{qwQVsBsXSdq+*0b8>*n|LLE@|^ionKIpo^`kjy ztJlUAd(YdfizVy#d9R=MJe;TtyKC!X$#6Cu*4E+R2)337&1CPr>t&;NsJ(q^E$P~C zZY{U3P1hy3?li3GX;oj#pUtmv<}5~^+y|m@w>G;<+rHb{F09}ln{?N?$#5(iPPw4( zy6|B3Z3w%3uyC`ZlT3=48=|A|+4A<)u=gI=-2Ej__bB%9=5U|;TKi;U(zjY!G>CIHRyd&ecI>U=DSOljl=0G4H?u} z!201L%EVID-f^^V(681+0qeCo>Mb&FzV&`^jc7Z~gtV*V}1FC^UOl<2$c)-a2xAA_M(ikGLFReea)R?e&&&|NRd4^YK>q zp~{Jk-hh$!I$pWe*HTWg;|4i3?2psZeuw(~>DEt6Im7ySzuD4pgZF1!9WCXYh__*6 z&K5sg+J67{=QjHLwUqOk{Ot}jTRML5{sM29P~?A|iyHlY4m6j(A>=U$MuE_ zMV20-T;AyScc2-K`1`MDdVhJKWF*8gd!ibbR+)!9BEnA==Xi3Wp$I^_Ysy})BLgb_oNa2D|^5H zmfQ3FznAU*eI@uj4?dTJ&#&NfDEK@G?z6%DCAcpH*IBUcg7c8D@mEK@2az?6ekPmz zOnDKCq$B?G+93D%&(}8knF(7qc^UkX>)W4+Nbq}4@O?h`o*R6>3%M|-gF+i`}(2&90}Ipl9dOzZVvJM?EiUxsOM|{ z@!$7#c9`cYINyhR-hw(>`kvO(efxv@J?w*8ckX+$f1k|#zTnTpA^JV3k^hF+w0Vp- zW-fE{poNfI2f1~STL-!8Aa@_g-3M~_f!uu{cOS^z2Xgm;BlkYQxpxpoD;x*?yV>KB zhZE2S1EGIgI|%x>v?roHPC`BgL;t_d5a|EcI2ro)zC+Oo!_XOp=z=2XFZG9`8;a2# zCFp@t^h6nsg8s{9FN{QQj6xrbMqlV6J{n`u5981ur(gh1#W6Sy$6`Ffm;m)0fqI?@ z^?W+ivr;`zf_he}=W?iLrFxzW^{iCSQ=p!e>Uk>E^E9Yu{mr|2J_G7`2GsLRsOMQw z&$FSP&xCrO1NA%?>Ukd2^L(i11yIimp`I5(JuilOUIO*3zd?*(8PxM~sKymgjb}kM z>c0T1#wb+dN~p#dRAU8HVW4NjR~m6)liLVpc<1 zDX7LYbWUn5*OEFX8R%TpLFXU~&3(P)I#RtgK)tPpdfNc?wh=3F4%FMZP;ci!y`2yB zb^+8|Bh=f4P%jrjy<80SatYMSrBE-ILA_iK^|A@-+xE~th0ceZ| zp<@Qee~8DLqlclnc?6o1N1?fR44Q+-q28Z>dVdnz@f1|k(@;&%Ks7xJ)$|-x)ALYG zyP%q0fNFXXs_7-D_m?eSAywn6&={{lV>Cfyybg`=1~kT-&=|X+F+K#1@nL9;k3eI5 z6dL1W&=?E8vhQ}_z$SYe?m3>3)JWks?pH{YIItW)abM(snN+J zsnKaeGB-|J%XZ{jP~Ftvv?r;-$tS77=|ED0Q$SLK(~+bGrxQsHPG^!DoG##@(-k~$ zx`7&;?w|&z2dKg632JbT0yQ|jKn+fBc4p3;J|uJI^d*@y=V+2SbNZ3YnbV(S&YS_D zR_7Q{t8*-<)d_=Io#Q~Q&hemD=LArzGZ3nE5LD}lP^~9HwGM`A9Rk&QGF0nOsMcYi zR;Q3Um@}scs$n=(LorlC2~IOCxjCRjZY%Za?kTso&)P9m98ryR_wGa1aOGX>14GZmWCX_nJT^>hZ*(+sGm znNUx&pq^$!J)H^lGzaQwF4WUJsHgc*PYa-)7D7EOf_hpE^|S=)X({N*Sq6G?mP0+Q zfOgjB#<|tJ2N~q=-RC5JXb0t)B9IDxcYOaE6u7+w}1=XB@YF-W1yauW{3DsN! z)trKAPD3@Uk5?^A%9fS3*5s1@(M2)blk^&(}gd zZ-#o_0`+_y)bsUF&o@9l-w5@56V&rosON1^&o@Io-vae~E7bFCP|vqRJ>LQKd?(cN zT~N>4p`P!CdcFth`Ch2!9Z=8rK|Sw;dcGg(`2ncs2ce!Hf_i=!>iH3<=SQKQAA@>+ z9P0TAsOKl4o}Yqxej4ie8K~!Hp`M?EdVU`2c^B043sBE5LOs6(_53o_^D9u#uR=Y) z2KC$o_53>2^BYjlZ$drq4)F(CnbxG_nKq=fHSI`gZ}Lg$UoYLtLa8c zchiHEp5`c0dYRs&^f7%&Iok9irN0?K$}#3xGK}NQ@uZw!29h$!oJh(^W-uv3%*mt- zHN!|LG)1HgH^rosm{L;8%m`9Ono*>THe*N`YsQgsiaC{()695s0wQK2DW_u!rkZJ_ zOgCqcGQ-RyWtN#u%9&;kDRa#{Qs$clq%1UxNLg%_kV~=5EGI>E$|O@x%4DM+WQ944 zl(S8gl$9n%N` zq?~KcBjtQ^0V$2xG5&3#FUa!W=4=Q5~DEMj3H&L8Ar+~=2TKnGvi5_U?QYUG^dj?$&`~*G0jXT zGCQgcLsz|BEDkRKmQr4IxDK#cVO4`(tveslssWVwp>diV*8q9i9Hkgg1 zoMX-<DQeH8ylJc5qBIR}S1}Sfv y-J!RCM{AXPy!ZFj++5_=MQ#n`)*1%y|1OEmS#=uel literal 0 HcmV?d00001 diff --git a/data/test/data/CasulliRefinement/casulli_polygon_then_all.nc b/data/test/data/CasulliRefinement/casulli_polygon_then_all.nc new file mode 100644 index 0000000000000000000000000000000000000000..23f38cd1fd0eff2c52d2b0227ac74ef1a3c4b128 GIT binary patch literal 12316 zcmeI1`F9glAI1kPg;L6H0RbI%kftQvXzRitrHCxDlpPtS$uu3B%!J9L3vLlbQB+VA za7RVl_YL=|h#M-d{K}v5KF`cIr4;%e&w1a|_q>PXoDZK&<}S~@x%Y;oa~m25DZyhn zQ~-wAOH6Nx6EpGrMDBW#Ww$BHpgp)Q=6zl-_0qUh2ItDmOPEQNIg0BMd<_&^cG|PN zjAuDEE`zcSX2LUEla+_NW{c^XcGSNL?w7WDUaEfP%#MzZ@|YAYw-cS)rkh58YNWiPJW?Ajk5pCER8}r7Q|DyjYPeFZ2v<~=sgYUr z6&3Z7sxmc+u?Ed6LtCx1=eS+D_RHtGrs0_}70sbqRVZ8?3eRc^*VR{4)mPWd2v>x| zy}I^ujfkO!9t|o=&WZh1#{Ck%JVID4mKpP#@aO5@)hV9q*zKmxY3+FridN2FzMvr# zF0TvPD9B_iHplm8T60TPolzYzL$!6am7(fzRZGaIsf>i`YK&-YozXli60HuZDaG(m zP`KDkw^lUBtHFC($a~Qi-E|WFo4~e9k|?KpPRdC*aew9%Wwebs)uWE%#w;7t!H_IF zPL`J5spg*fHAed7V^-3%`PQLiP`>wCEHz_sQH2dqA>}%oISR+tJA2I)>l7s%JKm@Kh|bWd2}K#a7lnsWl>PT)lllK|X6yPV z*}Be}ZPXEDn=IorW4#*E(_Lwg$3^h)$AG5=GDUnp$E>v5Y^TFkd+G;e3(QnH&axmk z&M|S`pkgwH;{qjTFjYOu`dYbHCt5#Yyw`7cj>NzS)tGO09U!e3> zxM#g=^_%J$C*KF|{r=YSu#}Uwuo84Q^*S7m z+aER>$D%N%VRdd&;w-Mk3LrW8$K1 zr`2&CZ1H#o+s!>yCCFR2NKvLA&z`>5aqHZ9Xyw6J2i%uQzg6NnS^tSP(7sOkZxAn( z90TnI(*F|i@@_sEqy<(6%*(RuSJ9e$Dcb68|H*)RRW|5%tuIB}&>dV4wAW{YemDA3 zv`r_2>jAeG^lM}TM>qX%+}s`X8}p@T=E>lCz^w)SwqygxyU$7cgT_Bv-}AKZxz_i5 z>ATPS?zc16y%x-GZvCCH4&;8{_rvn5DcWYa-w2omdw)s|p`l>_FhS8rYKSgk0Nlo!TAT z-*P;#+kKkzIr8(t1)=5FN3bvQ`z4r%e|UcQg}HvtvGY5N|J!!&;#XKj$vZy}`E`)r z2l;a#e-8XVoCCJSj^u)F9rqTA@C)Dzmd}XP4rSQOk4!>6Q6%~rTKTa%K?A4 zrAxpSV5FFTJGc^z7LS3afwAIo@N_U2kDXMl;~NpKmsKs*@^fhpptFn<}m zP<#L=bX1f+J z6Tita%WS~DHDM0S4A%o@mBJ60%lkUm|+oMhQ)vx zmH=kB2{6M_z`mCW%VB0%0hpl)FvCi~466V$tOm@m1~9`~zzpjEGpq;9umLc`&43v; z0%p(wGi(CPU;y^rEJR^uhyi9W0W-7!W{3l3Xa&q*0cO|?n4t|YLjo{E5-@`en85+e zkOG^)7QhTHVBcxMgP9=%n4ujoLkD1nPQVOZfEl&|X1E0~!>xcBZUfA4J79)805fC( zGu#Q7VH;qEy8!#XTi6aW!##i*?gh-S12DsVfEn%w%&-$M!vlaB9t6zr5MYK~fEjj! zEO;0&!y|wh_5fyh6tM5dguO5`>;ufOA27q?fEk_u%&hNl2CJPnxP8NdwB0%kY> znBh6V49^2*cmXiOLBI?z0``4KcnM~PmjN@p0+``dzznYeW_TSi!yAAZ-UQ6>7GQ?A z!2$3NV1{=AGrR|w;eEgi9{^_f5U}r$gu^g1d@TM1W|mLIpTW%Zx%dm1*}fEa!_4@V z_-mM1zY!mSnfY7sQJC4k6MqkL3_pl}ggKU<#K&Nc>1XjTFvs?*_&Cfl{wDq%=2-s_ z{|R%3zxv;QsFX~TrOMzlPiU3Fqtk~ z1(O*<8B9U~`y^p78)zT`DuhayR0-8EsS#$uq*kbdNxg72Os)~Gg~@dS`y_L~Qm_nM z59SIDFqtPb!eqX1156eOH^OA0um~oLg(Wb#NnoF3IWRyoSOJ=Zl`vT)tcJ-NVJ%G7 z3F~39LAV(v8wDLEn*{brqQC`d5Cf*r0+YDV3KL7%43jn?0h6R)!^9C%FxeuoPvU{Q zz}+AN+Jz37bP8QC*(%%ulUs$`U~;=~2TZcUoiN!ZuurlbJPIBI_keqa9Wc31xF04v zg$H2rpzsh(b_u&-^04p-O!f%ulk5dAf#q z5k7^. +// +// contact: delft3d.support@deltares.nl +// Stichting Deltares +// P.O. Box 177 +// 2600 MH Delft, The Netherlands +// +// All indications and logos of, and references to, "Delft3D" and "Deltares" +// are registered trademarks of Stichting Deltares, and remain the property of +// Stichting Deltares. All rights reserved. +// +//------------------------------------------------------------------------------ + +#pragma once + +#include +#include + +#include "MeshKernel/Mesh2D.hpp" +#include "MeshKernel/Polygons.hpp" +#include "MeshKernel/UndoActions/UndoAction.hpp" + +namespace meshkernel +{ + /// @brief Compute the Casulli de-refinement for a mesh. + class CasulliDeRefinement + { + public: + /// @brief Compute the Casulli de-refinement for the entire mesh. + /// + /// @param [in, out] mesh Mesh to be de-refined + [[nodiscard]] std::unique_ptr Compute(Mesh2D& mesh) const; + + /// @brief Compute the Casulli de-refinement for the part of the mesh inside the polygon + /// + /// @param [in, out] mesh Mesh to be de-refined + /// @param [in] polygon Area within which the mesh will be de-refined + [[nodiscard]] std::unique_ptr Compute(Mesh2D& mesh, const Polygons& polygon) const; + + /// @brief Compute centres of elements to be deleted. + /// + /// Requires that the element mass-centres are pre-computed. + std::vector ElementsToDelete(const Mesh2D& mesh, const Polygons& polygon) const; + + private: + /// @brief Maximum number of iterations allowed when initialising the element mask + static const UInt maxIterationCount = 1000; + + /// @brief Maximum reserve size for arrays used in the de-refinement + static const UInt maximumSize = 1000; + + /// @brief Label for the element in the de-refined grid. + /// + /// The prefix, e.g. WasNode, indicates the original mesh entity around/from which the de-refined element was constructed. + /// First and after suffix, indicates the order in which the elements are to be processed. + /// Enumeration values and comments are from the original Fortran code. + enum class ElementType + { + WasNodeFirst = 1, //< front, 'A' cell (used to be node, delete it): 1 + WasEdgeFirst = 2, //< front, 'B' cell (used to be link, keep it): 2 + WasCell = 3, //< 'C' cell (used to be cell, keep it): 3 + WasNodeAfter = -1, //< not in front, 'A' cell: -1 + WasEdgeAfter = -2, //< not in front, 'B' cell: -2 + Unassigned = 0 //< not assigned a value 0 + }; + + /// @brief Indicate if the element can be a seed element or not. + bool ElementIsSeed(const Mesh2D& mesh, + const std::vector& nodeTypes, + const UInt element) const; + + /// @brief Find the seed element id to start the mesh de-refinement. + UInt FindElementSeedIndex(const Mesh2D& mesh, + const std::vector& nodeTypes) const; + + /// @brief Find all elements that are connected along edges to elementId. + void FindDirectlyConnectedCells(const Mesh2D& mesh, + const UInt elementId, + std::vector& connected) const; + + /// @brief Find all elements that are connected by nodes to elementId. + void FindIndirectlyConnectedCells(const Mesh2D& mesh, + const UInt elementId, + const std::vector& directlyConnected, + std::vector& indirectlyConnected) const; + + /// @brief Find element id's + void FindAdjacentCells(const Mesh2D& mesh, + const std::vector& directlyConnected, + const std::vector& indirectlyConnected, + std::vector>& kne) const; + + /// @brief Find the elements that are connected to the elementId. + void FindSurroundingCells(const Mesh2D& mesh, + const UInt elementId, + std::vector& directlyConnected, + std::vector& indirectlyConnected, + std::vector>& kne) const; + + /// @brief Initialise the element mask. + std::vector InitialiseElementType(const Mesh2D& mesh, + const std::vector& nodeTypes) const; + + /// \brief Determine if the element can be deleted from the mesh or not. + bool ElementCannotBeDeleted(const Mesh2D& mesh, + const std::vector& nodeTypes, + const Polygons& polygon, + const UInt elementId) const; + + /// @brief Compute coordinates of the new node. + Point ComputeNewNodeCoordinates(const Mesh2D& mesh, + const std::vector& nodeTypes, + const UInt nodeId) const; + + /// @brief Update the mesh members for the mesh description and connectivity. + [[nodiscard]] bool UpdateDirectlyConnectedElements(Mesh2D& mesh, + const UInt elementId, + const std::vector& directlyConnected, + const std::vector>& kne) const; + + /// @brief Update the mesh members for the mesh description and connectivity for triangle elements + bool UpdateDirectlyConnectedTriangleElements(Mesh2D& mesh, + const UInt index, + const UInt connectedElementId, + const std::vector>& kne) const; + + /// @brief Update the mesh members for the mesh description and connectivity for non-triangle elements + /// + /// That is, element with 4 or more edges. + void UpdateDirectlyConnectedNonTriangleElements(Mesh2D& mesh, + const UInt index, + const UInt elementId, + const UInt connectedElementId) const; + + /// @brief Get the most significant node type for all nodes of the element. + int GetNodeCode(const Mesh2D& mesh, + const std::vector& nodeTypes, + const UInt elementId) const; + + /// @brief Add element id to the list of id's + /// + /// only added is it is not already on the list and the element is a quadrilateral + void AddElementToList(const Mesh& mesh, const std::vector& frontList, std::vector& frontListCopy, const UInt elementId) const; + + /// @brief Redirect nodes of connected cells, deactivate polygons of degree smaller than three + void RedirectNodesOfConnectedElements(Mesh2D& mesh, const UInt elementId, const UInt nodeId, const std::vector& indirectlyConnected) const; + + /// @brief Removes nodes from the boundary that will not be part of the de-refined mesh. + [[nodiscard]] bool RemoveUnwantedBoundaryNodes(Mesh2D& mesh, + const std::vector& nodeTypes, + const Polygons& polygon, + const std::vector& indirectlyConnected) const; + + /// @brief Delete an element + [[nodiscard]] bool DeleteElement(Mesh2D& mesh, + std::vector& nodeTypes, + const Polygons& polygon, + const UInt elementId, + const std::vector& directlyConnected, + const std::vector& indirectlyConnected, + const std::vector>& kne) const; + + /// @brief Clean up the edge + /// + /// @returns Indicates if the cleanp-up was successful or not + [[nodiscard]] bool CleanUpEdge(Mesh2D& mesh, const UInt edgeId) const; + + /// @brief Do the Casullu de-refinement + [[nodiscard]] bool DoDeRefinement(Mesh2D& mesh, const Polygons& polygon) const; + + /// @brief Compute the mesh node types. + /// + /// Uses the m_nodeTypes has been generated in the mesh. + std::vector ComputeNodeTypes(const Mesh2D& mesh, const Polygons& polygon) const; + }; + +} // namespace meshkernel diff --git a/libs/MeshKernel/include/MeshKernel/Entities.hpp b/libs/MeshKernel/include/MeshKernel/Entities.hpp index 3b4a3a957..3dda58893 100644 --- a/libs/MeshKernel/include/MeshKernel/Entities.hpp +++ b/libs/MeshKernel/include/MeshKernel/Entities.hpp @@ -31,6 +31,7 @@ #include "MeshKernel/Constants.hpp" #include "MeshKernel/Definitions.hpp" +#include "MeshKernel/Exceptions.hpp" #include "MeshKernel/Point.hpp" namespace meshkernel @@ -57,6 +58,25 @@ namespace meshkernel return node == edge.first ? edge.second : edge.first; } + /// @brief Get the node at the position + /// + /// Enables easier accessing of nodes of edge in loop. + static UInt EdgeNodeIndex(const Edge& edge, UInt position) + { + if (position == 0) + { + return edge.first; + } + else if (position == 1) + { + return edge.second; + } + else + { + throw ConstraintError("Position out of bounds: {} not in [0 .. 1]", position); + } + } + /// @brief A struct describing the three coordinates in a cartesian projection. struct Cartesian3DPoint { diff --git a/libs/MeshKernel/include/MeshKernel/Mesh.hpp b/libs/MeshKernel/include/MeshKernel/Mesh.hpp index a28a224ce..4a9c122c0 100644 --- a/libs/MeshKernel/include/MeshKernel/Mesh.hpp +++ b/libs/MeshKernel/include/MeshKernel/Mesh.hpp @@ -37,6 +37,7 @@ #include "MeshKernel/UndoActions/CompoundUndoAction.hpp" #include "MeshKernel/UndoActions/DeleteEdgeAction.hpp" #include "MeshKernel/UndoActions/DeleteNodeAction.hpp" +#include "MeshKernel/UndoActions/FullUnstructuredGridUndo.hpp" #include "MeshKernel/UndoActions/MeshConversionAction.hpp" #include "MeshKernel/UndoActions/NodeTranslationAction.hpp" #include "MeshKernel/UndoActions/ResetEdgeAction.hpp" @@ -189,12 +190,18 @@ namespace meshkernel /// @brief Set all nodes to a new set of values. void SetNodes(const std::vector& newValues); + /// @brief Set a node to a new value, bypassing the undo action. + void SetNode(const UInt index, const Point& newValue); + /// @brief Set the node to a new value, this value may be the in-valid value. [[nodiscard]] std::unique_ptr ResetNode(const UInt index, const Point& newValue); - /// @brief Get the edge + /// @brief Get constant reference to an edge const Edge& GetEdge(const UInt index) const; + /// @brief Get a non-constant reference to an edge + Edge& GetEdge(const UInt index); + /// @brief Get all edges // TODO get rid of this function const std::vector& Edges() const; @@ -379,6 +386,9 @@ namespace meshkernel /// @brief Apply the delete edge action void CommitAction(const DeleteEdgeAction& undoAction); + /// @brief Set the node and edge values. + void CommitAction(FullUnstructuredGridUndo& undoAction); + /// @brief Undo the reset node action /// /// Restore mesh to state before node was reset @@ -419,6 +429,11 @@ namespace meshkernel /// Restore mesh to state before edge was deleted void RestoreAction(const DeleteEdgeAction& undoAction); + /// @brief Undo entire node and edge values + /// + /// Restore mesh to previous state. + void RestoreAction(FullUnstructuredGridUndo& undoAction); + /// @brief Get a reference to the RTree for a specific location RTreeBase& GetRTree(Location location) const { return *m_RTrees.at(location); } @@ -527,6 +542,16 @@ inline const meshkernel::Edge& meshkernel::Mesh::GetEdge(const UInt index) const return m_edges[index]; } +inline meshkernel::Edge& meshkernel::Mesh::GetEdge(const UInt index) +{ + if (index >= GetNumEdges()) + { + throw ConstraintError("The edge index, {}, is not in range.", index); + } + + return m_edges[index]; +} + inline const std::vector& meshkernel::Mesh::Edges() const { return m_edges; diff --git a/libs/MeshKernel/include/MeshKernel/Operations.hpp b/libs/MeshKernel/include/MeshKernel/Operations.hpp index 7538932ac..374fbc68e 100644 --- a/libs/MeshKernel/include/MeshKernel/Operations.hpp +++ b/libs/MeshKernel/include/MeshKernel/Operations.hpp @@ -241,13 +241,13 @@ namespace meshkernel return f1 < f2 ? x1 : x2; } - /// @brief Get the next forward index. + /// @brief Get the next forward index, for a range in 0 .. size-1. /// @param[in] currentIndex The current index. /// @param[in] size The size of the vector. /// @returns The next forward index. [[nodiscard]] UInt NextCircularForwardIndex(UInt currentIndex, UInt size); - /// @brief Get the next backward index. + /// @brief Get the next backward index, for a range in 0 .. size-1. /// @param[in] currentIndex The current index. /// @param[in] size The size of the vector. /// @returns The next backward index. diff --git a/libs/MeshKernel/include/MeshKernel/Polygons.hpp b/libs/MeshKernel/include/MeshKernel/Polygons.hpp index 58e8568b1..4d2645dc8 100644 --- a/libs/MeshKernel/include/MeshKernel/Polygons.hpp +++ b/libs/MeshKernel/include/MeshKernel/Polygons.hpp @@ -102,6 +102,11 @@ namespace meshkernel /// @return True if it is included, false otherwise [[nodiscard]] bool IsPointInPolygon(Point const& point, UInt polygonIndex) const; + /// @brief Checks if a point is included in any of the polygonal enclosures contained. + /// @param[in] point The point to check + /// @return True if it is included, false otherwise + [[nodiscard]] bool IsPointInAnyPolygon(const Point& point) const; + // TODO can reduce the result of this function to only a UInt (valid value => found, invalid => not found) /// @brief Checks if a point is included in any of the polygons (dbpinpol_optinside_perpol) /// @param[in] point The point to check diff --git a/libs/MeshKernel/include/MeshKernel/UndoActions/FullUnstructuredGridUndo.hpp b/libs/MeshKernel/include/MeshKernel/UndoActions/FullUnstructuredGridUndo.hpp new file mode 100644 index 000000000..0439e2840 --- /dev/null +++ b/libs/MeshKernel/include/MeshKernel/UndoActions/FullUnstructuredGridUndo.hpp @@ -0,0 +1,68 @@ +//---- GPL --------------------------------------------------------------------- +// +// Copyright (C) Stichting Deltares, 2011-2024. +// +// 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 version 3. +// +// 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 . +// +// contact: delft3d.support@deltares.nl +// Stichting Deltares +// P.O. Box 177 +// 2600 MH Delft, The Netherlands +// +// All indications and logos of, and references to, "Delft3D" and "Deltares" +// are registered trademarks of Stichting Deltares, and remain the property of +// Stichting Deltares. All rights reserved. +// +//------------------------------------------------------------------------------ + +#pragma once + +#include +#include + +#include "MeshKernel/Entities.hpp" +#include "MeshKernel/Point.hpp" +#include "MeshKernel/UndoActions/BaseMeshUndoAction.hpp" + +namespace meshkernel +{ + /// @brief Forward declaration of the unstructured mesh + class Mesh; + + /// @brief Action to save all the node and edge values from a mesh + /// + /// The undo will simply swap these values + class FullUnstructuredGridUndo : public BaseMeshUndoAction + { + public: + /// @brief Allocate a DeleteNodeAction and return a unique_ptr to the newly created object. + static std::unique_ptr Create(Mesh& mesh); + + /// @brief Constructor + explicit FullUnstructuredGridUndo(Mesh& mesh); + + /// @brief Swap the nodes and edges with the saved values. + void Swap(std::vector& nodes, std::vector& edges); + + /// \brief Compute the approximate amount of memory being used, in bytes. + std::uint64_t MemorySize() const override; + + private: + /// @brief The saved node values. + std::vector m_savedNodes; + + /// @brief The saved edge values. + std::vector m_savedEdges; + }; + +} // namespace meshkernel diff --git a/libs/MeshKernel/include/MeshKernel/Utilities/NumericFunctions.hpp b/libs/MeshKernel/include/MeshKernel/Utilities/NumericFunctions.hpp index 5026b4476..1988f09b1 100644 --- a/libs/MeshKernel/include/MeshKernel/Utilities/NumericFunctions.hpp +++ b/libs/MeshKernel/include/MeshKernel/Utilities/NumericFunctions.hpp @@ -27,6 +27,7 @@ #pragma once +#include "MeshKernel/Exceptions.hpp" #include #include #include @@ -65,4 +66,36 @@ namespace meshkernel return lowerBound <= value && value <= upperBound; } + /// \brief Get the next index, wrapping around if index is at + /// + /// Range is in [0, size - 1] + static UInt RotateIndex(const UInt index, const UInt size, const bool forward) + { + if (size == 0) + { + throw ConstraintError("Invalid range for index rotation"); + } + + if (index >= size) + { + throw ConstraintError("Index is out of range of array: {} not in [0 .. {}]", index, size - 1); + } + + if (forward) + { + return index == size - 1 ? 0 : index + 1; + } + else + { + return index == 0 ? size - 1 : index - 1; + } + } + + /// \brief Get the next index, wrapping around if index is at + template + UInt RotateIndex(const UInt index, const std::vector& vec, const bool forward) + { + return RotateIndex(index, static_cast(vec.size()), forward); + } + } // namespace meshkernel diff --git a/libs/MeshKernel/src/CasulliDeRefinement.cpp b/libs/MeshKernel/src/CasulliDeRefinement.cpp new file mode 100644 index 000000000..2ffdcf92d --- /dev/null +++ b/libs/MeshKernel/src/CasulliDeRefinement.cpp @@ -0,0 +1,1013 @@ +//---- GPL --------------------------------------------------------------------- +// +// Copyright (C) Stichting Deltares, 2011-2024. +// +// 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 version 3. +// +// 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 . +// +// contact: delft3d.support@deltares.nl +// Stichting Deltares +// P.O. Box 177 +// 2600 MH Delft, The Netherlands +// +// All indications and logos of, and references to, "Delft3D" and "Deltares" +// are registered trademarks of Stichting Deltares, and remain the property of +// Stichting Deltares. All rights reserved. +// +//------------------------------------------------------------------------------ + +#include + +#include "MeshKernel/CasulliDeRefinement.hpp" +#include "MeshKernel/Exceptions.hpp" +#include "MeshKernel/Operations.hpp" +#include "MeshKernel/UndoActions/AddNodeAction.hpp" +#include "MeshKernel/UndoActions/CompoundUndoAction.hpp" +#include "MeshKernel/UndoActions/FullUnstructuredGridUndo.hpp" + +std::unique_ptr meshkernel::CasulliDeRefinement::Compute(Mesh2D& mesh) const +{ + Polygons emptyPolygon; + return Compute(mesh, emptyPolygon); +} + +std::unique_ptr meshkernel::CasulliDeRefinement::Compute(Mesh2D& mesh, const Polygons& polygon) const +{ + std::unique_ptr refinementAction = FullUnstructuredGridUndo::Create(mesh); + const std::vector meshNodes(mesh.Nodes()); + const std::vector meshEdges(mesh.Edges()); + + if (DoDeRefinement(mesh, polygon)) + { + mesh.DeleteInvalidNodesAndEdges(); + mesh.Administrate(); + } + else + { + // The de-refinement failed, restore the original mesh + refinementAction->Restore(); + // Reset the undo action to a null action + refinementAction.reset(nullptr); + throw AlgorithmError("Unable to compute the Casulli de-refinement"); + } + + return refinementAction; +} + +void meshkernel::CasulliDeRefinement::FindDirectlyConnectedCells(const Mesh2D& mesh, + const UInt elementId, + std::vector& connected) const +{ + connected.clear(); + + for (UInt i = 0; i < mesh.m_numFacesNodes[elementId]; ++i) + { + UInt edgeId = mesh.m_facesEdges[elementId][i]; + + if (mesh.m_edgesNumFaces[edgeId] < 2) + { + continue; + } + + UInt neighbouringElementId = elementId == mesh.m_edgesFaces[edgeId][0] ? mesh.m_edgesFaces[edgeId][1] : mesh.m_edgesFaces[edgeId][0]; + connected.push_back(neighbouringElementId); + } +} + +void meshkernel::CasulliDeRefinement::FindIndirectlyConnectedCells(const Mesh2D& mesh, + const UInt elementId, + const std::vector& directlyConnected, + std::vector& indirectlyConnected) const +{ + indirectlyConnected.clear(); + + for (UInt i = 0; i < mesh.m_numFacesNodes[elementId]; ++i) + { + UInt nodeId = mesh.m_facesNodes[elementId][i]; + + for (UInt j = 0; j < mesh.m_nodesNumEdges[nodeId]; ++j) + { + UInt edgeId = mesh.m_nodesEdges[nodeId][j]; + + for (UInt k = 0; k < mesh.m_edgesNumFaces[edgeId]; ++k) + { + UInt otherElementId = mesh.m_edgesFaces[edgeId][k]; + + if (elementId == otherElementId || + std::ranges::find(directlyConnected, otherElementId) != directlyConnected.end() || + std::ranges::find(indirectlyConnected, otherElementId) != indirectlyConnected.end()) + { + continue; + } + + // Add new cell + indirectlyConnected.push_back(otherElementId); + } + } + } +} + +void meshkernel::CasulliDeRefinement::FindAdjacentCells(const Mesh2D& mesh, + const std::vector& directlyConnected, + const std::vector& indirectlyConnected, + std::vector>& kne) const +{ + std::ranges::fill(kne, std::array{constants::missing::intValue, constants::missing::intValue}); + + for (UInt i = 0; i < directlyConnected.size(); ++i) + { + UInt elementId = directlyConnected[i]; + + for (UInt j = 0; j < mesh.m_numFacesNodes[elementId]; ++j) + { + UInt edgeId = mesh.m_facesEdges[elementId][j]; + + if (mesh.m_edgesNumFaces[edgeId] < 2) + { + continue; + } + + UInt neighbouringElementId = elementId == mesh.m_edgesFaces[edgeId][0] ? mesh.m_edgesFaces[edgeId][1] : mesh.m_edgesFaces[edgeId][0]; + + for (UInt k = 0; k < directlyConnected.size(); ++k) + { + if (directlyConnected[k] == neighbouringElementId) + { + if (kne[i][0] == constants::missing::intValue) + { + kne[i][0] = -static_cast(neighbouringElementId); + } + else + { + kne[i][1] = -static_cast(neighbouringElementId); + } + + neighbouringElementId = constants::missing::uintValue; + } + } + + if (neighbouringElementId == constants::missing::uintValue) + { + continue; + } + + for (UInt k = 0; k < indirectlyConnected.size(); ++k) + { + if (indirectlyConnected[k] == neighbouringElementId) + { + if (kne[i][0] == constants::missing::intValue) + { + kne[i][0] = static_cast(neighbouringElementId); + } + else + { + kne[i][1] = static_cast(neighbouringElementId); + } + } + } + } + } +} + +void meshkernel::CasulliDeRefinement::FindSurroundingCells(const Mesh2D& mesh, + const UInt elementId, + std::vector& directlyConnected, + std::vector& indirectlyConnected, + std::vector>& kne) const +{ + FindDirectlyConnectedCells(mesh, elementId, directlyConnected); + FindIndirectlyConnectedCells(mesh, elementId, directlyConnected, indirectlyConnected); + FindAdjacentCells(mesh, directlyConnected, indirectlyConnected, kne); +} + +bool meshkernel::CasulliDeRefinement::ElementIsSeed(const Mesh2D& mesh, + const std::vector& nodeTypes, + const UInt element) const +{ + bool isSeed = true; + + for (UInt i = 0; i < mesh.m_numFacesNodes[element]; ++i) + { + if (nodeTypes[mesh.m_facesNodes[element][i]] == 0) + { + isSeed = false; + break; + } + } + + return isSeed; +} + +meshkernel::UInt meshkernel::CasulliDeRefinement::FindElementSeedIndex(const Mesh2D& mesh, + const std::vector& nodeTypes) const +{ + UInt seedIndex = constants::missing::uintValue; + + for (UInt e = 0; e < mesh.Edges().size(); ++e) + { + if (mesh.m_edgesNumFaces[e] != 1) + { + continue; + } + + if (const Edge& edge = mesh.GetEdge(e); nodeTypes[edge.first] != 2 || nodeTypes[edge.second] != 2) + { + continue; + } + + UInt elementId = mesh.m_edgesFaces[e][0]; + + if (mesh.m_numFacesNodes[elementId] != constants::geometric::numNodesInQuadrilateral) + { + continue; + } + + if (ElementIsSeed(mesh, nodeTypes, elementId)) + { + seedIndex = elementId; + break; + } + } + + // No seed index found, select the first quadrilateral inside the selecting polygon + if (seedIndex == constants::missing::uintValue) + { + for (UInt face = 0; face < mesh.GetNumFaces(); ++face) + { + + if (mesh.m_numFacesNodes[face] != constants::geometric::numNodesInQuadrilateral) + { + continue; + } + + if (!ElementIsSeed(mesh, nodeTypes, face)) + { + continue; + } + + seedIndex = face; + break; + } + } + + // still no element found, so take the first + if (seedIndex == constants::missing::uintValue) + { + seedIndex = 0; + } + + return seedIndex; +} + +void meshkernel::CasulliDeRefinement::AddElementToList(const Mesh& mesh, const std::vector& frontList, std::vector& frontListCopy, const UInt elementId) const +{ + if (elementId != constants::missing::uintValue) + { + + if (mesh.m_numFacesNodes[elementId] != constants::geometric::numNodesInQuadrilateral) + { + return; + } + + // If the id is not already on the frontList then add it to the copy. + if (std::ranges::find(frontList, elementId) == frontList.end()) + { + frontListCopy.push_back(elementId); + } + } +} + +std::vector meshkernel::CasulliDeRefinement::InitialiseElementType(const Mesh2D& mesh, + const std::vector& nodeTypes) const +{ + UInt seedElement = FindElementSeedIndex(mesh, nodeTypes); + UInt iterationCount = 0; + + std::vector directlyConnected; + std::vector indirectlyConnected; + std::vector> kne(maximumSize); + std::vector frontIndex; + std::vector frontIndexCopy; + + directlyConnected.reserve(maximumSize); + indirectlyConnected.reserve(maximumSize); + frontIndex.reserve(maximumSize); + frontIndexCopy.reserve(maximumSize); + + std::vector cellMask(mesh.GetNumFaces(), ElementType::Unassigned); + + cellMask[seedElement] = ElementType::WasNodeFirst; + frontIndex.push_back(seedElement); + + while (frontIndex.size() > 0 && iterationCount < maxIterationCount) + { + ++iterationCount; + frontIndexCopy.clear(); + + for (UInt i = 0; i < frontIndex.size(); ++i) + { + UInt elementId = frontIndex[i]; + + // get the connected cells + FindSurroundingCells(mesh, elementId, directlyConnected, indirectlyConnected, kne); + + if (cellMask[elementId] == ElementType::WasNodeFirst) + { + + for (UInt j = 0; j < directlyConnected.size(); ++j) + { + UInt connectedElementId = directlyConnected[j]; + + if (mesh.m_numFacesNodes[connectedElementId] != constants::geometric::numNodesInQuadrilateral) + { + continue; + } + + if ((cellMask[connectedElementId] != ElementType::WasNodeFirst && cellMask[connectedElementId] != ElementType::WasNodeAfter) && + (cellMask[connectedElementId] != ElementType::WasEdgeFirst && cellMask[connectedElementId] != ElementType::WasEdgeAfter)) + { + cellMask[connectedElementId] = ElementType::WasEdgeFirst; + AddElementToList(mesh, frontIndex, frontIndexCopy, connectedElementId); + } + } + + for (UInt j = 0; j < indirectlyConnected.size(); ++j) + { + UInt connectedElementId = indirectlyConnected[j]; + + if (mesh.m_numFacesNodes[connectedElementId] != constants::geometric::numNodesInQuadrilateral) + { + continue; + } + + if (cellMask[connectedElementId] != ElementType::WasCell) + { + cellMask[connectedElementId] = ElementType::WasCell; + } + } + + cellMask[elementId] = ElementType::WasNodeAfter; + } + else if (cellMask[elementId] == ElementType::WasEdgeFirst) + { + + for (UInt j = 0; j < directlyConnected.size(); ++j) + { + UInt connectedElementId = directlyConnected[j]; + + if (mesh.m_numFacesNodes[connectedElementId] != constants::geometric::numNodesInQuadrilateral) + { + continue; + } + + if ((cellMask[connectedElementId] != ElementType::WasCell) && + (cellMask[connectedElementId] != ElementType::WasNodeFirst && cellMask[connectedElementId] != ElementType::WasNodeAfter) && + (cellMask[connectedElementId] != ElementType::WasEdgeFirst && cellMask[connectedElementId] != ElementType::WasEdgeAfter)) + { + cellMask[connectedElementId] = ElementType::WasNodeFirst; + AddElementToList(mesh, frontIndex, frontIndexCopy, connectedElementId); + } + } + + for (UInt j = 0; j < indirectlyConnected.size(); ++j) + { + UInt connectedElementId = indirectlyConnected[j]; + + if (mesh.m_numFacesNodes[connectedElementId] != constants::geometric::numNodesInQuadrilateral) + { + continue; + } + + if ((cellMask[connectedElementId] != ElementType::WasEdgeFirst && cellMask[connectedElementId] != ElementType::WasEdgeAfter) && + (cellMask[connectedElementId] != ElementType::WasNodeFirst && cellMask[connectedElementId] != ElementType::WasNodeAfter) && + cellMask[connectedElementId] != ElementType::WasCell) + { + cellMask[connectedElementId] = ElementType::WasEdgeFirst; + AddElementToList(mesh, frontIndex, frontIndexCopy, connectedElementId); + } + } + + cellMask[elementId] = ElementType::WasEdgeAfter; + } + } + + frontIndex = frontIndexCopy; + } + + return cellMask; +} + +bool meshkernel::CasulliDeRefinement::DoDeRefinement(Mesh2D& mesh, const Polygons& polygon) const +{ + std::vector directlyConnected; + std::vector indirectlyConnected; + std::vector> kne(maximumSize); + std::vector nodeTypes(ComputeNodeTypes(mesh, polygon)); + directlyConnected.reserve(maximumSize); + indirectlyConnected.reserve(maximumSize); + + std::vector cellMask(InitialiseElementType(mesh, nodeTypes)); + mesh.ComputeCircumcentersMassCentersAndFaceAreas(true); + + for (UInt k = 0; k < cellMask.size(); ++k) + { + if (cellMask[k] == ElementType::WasNodeAfter && mesh.m_numFacesNodes[k] > 0) + { + FindSurroundingCells(mesh, k, directlyConnected, indirectlyConnected, kne); + + if (!DeleteElement(mesh, nodeTypes, polygon, k, directlyConnected, indirectlyConnected, kne)) + { + return false; + } + } + } + + return true; +} + +bool meshkernel::CasulliDeRefinement::ElementCannotBeDeleted(const Mesh2D& mesh, + const std::vector& nodeTypes, + const Polygons& polygon, + const UInt elementId) const +{ + bool noGo = false; + + // Firstly, check and see if + // the cell to be deleted is not a corner cell, but has a link that + // is internal and whose both nodes are marked as non-internal nodes + // return if so + + for (UInt i = 0; i < mesh.m_numFacesNodes[elementId]; ++i) + { + UInt edgeId = mesh.m_facesEdges[elementId][i]; + + if (mesh.GetEdge(edgeId).first == constants::missing::uintValue || + mesh.GetEdge(edgeId).second == constants::missing::uintValue) + { + noGo = true; + break; + } + + if (mesh.m_edgesNumFaces[edgeId] == 2 && + nodeTypes[mesh.GetEdge(edgeId).first] != 1 && + nodeTypes[mesh.GetEdge(edgeId).second] != 1) + { + noGo = true; + break; + } + } + + //-------------------------------- + + for (UInt i = 0; i < mesh.m_numFacesNodes[elementId]; ++i) + { + UInt nodeId = mesh.m_facesNodes[elementId][i]; + + if (nodeTypes[nodeId] == 3 && mesh.m_nodesNumEdges[nodeId] <= 2) + { + noGo = false; + break; + } + } + + //-------------------------------- + + // check if all nodes are in the selecting polygon + for (UInt i = 0; i < mesh.m_numFacesNodes[elementId]; ++i) + { + UInt nodeId = mesh.m_facesNodes[elementId][i]; + + if (!polygon.IsPointInAnyPolygon(mesh.Node(nodeId))) + { + noGo = true; + break; + } + } + + return noGo; +} + +meshkernel::Point meshkernel::CasulliDeRefinement::ComputeNewNodeCoordinates(const Mesh2D& mesh, + const std::vector& nodeTypes, + const UInt elementId) const +{ + + Point newNode(0.0, 0.0); + double factor = 0.0; + + for (UInt i = 0; i < mesh.m_numFacesNodes[elementId]; ++i) + { + double fac = 1.0; + UInt nodeId = mesh.m_facesNodes[elementId][i]; + + if (nodeTypes[nodeId] == 2 || nodeTypes[nodeId] == 4) + { + fac = 1.0e45; + } + else if (nodeTypes[nodeId] == 3) + { + factor = 1.0; + newNode = mesh.Node(nodeId); + break; + } + + newNode += fac * mesh.Node(nodeId); + factor += fac; + } + + newNode /= factor; + + return newNode; +} + +bool meshkernel::CasulliDeRefinement::UpdateDirectlyConnectedTriangleElements(Mesh2D& mesh, + const UInt index, + const UInt connectedElementId, + const std::vector>& kne) const +{ + + for (UInt j = 0; j < mesh.m_numFacesNodes[connectedElementId]; ++j) + { + UInt edgeId = mesh.m_facesEdges[connectedElementId][j]; + + if (mesh.m_edgesNumFaces[edgeId] < 2 && !CleanUpEdge(mesh, edgeId)) + { + return false; + } + } + + // Find adjacent direct neighbours + UInt otherEdgeId = constants::missing::uintValue; + + for (UInt i = 0; i < 2; ++i) + { + UInt leftElementId = kne[index][i] == constants::missing::intValue || kne[index][i] < 0 ? constants::missing::uintValue : static_cast(kne[index][i]); + + if (leftElementId == constants::missing::uintValue) + { + continue; + } + + UInt oppositeSide = 1 - i; + UInt rightElementId = kne[index][oppositeSide] == constants::missing::intValue || kne[index][oppositeSide] < 0 ? constants::missing::uintValue : static_cast(kne[index][oppositeSide]); + + if (leftElementId == constants::missing::uintValue || rightElementId == constants::missing::uintValue) + { + continue; + } + + UInt j; + UInt edgeId = constants::missing::uintValue; + + // find the common link + for (j = 0; j < mesh.m_numFacesNodes[leftElementId]; ++j) + { + edgeId = mesh.m_facesEdges[leftElementId][j]; + + if (mesh.m_edgesNumFaces[edgeId] < 2) + { + continue; + } + + if (mesh.m_edgesFaces[edgeId][0] == connectedElementId && mesh.m_edgesFaces[edgeId][1] == leftElementId) + { + if (rightElementId != constants::missing::uintValue) + { + mesh.m_edgesFaces[edgeId][0] = rightElementId; + } + else + { + mesh.m_edgesFaces[edgeId][0] = mesh.m_edgesFaces[edgeId][1]; + mesh.m_edgesFaces[edgeId][1] = constants::missing::uintValue; + mesh.m_edgesNumFaces[edgeId] = 1; + } + + break; + } + else if (mesh.m_edgesFaces[edgeId][1] == connectedElementId && mesh.m_edgesFaces[edgeId][0] == leftElementId) + { + if (rightElementId != constants::missing::uintValue) + { + mesh.m_edgesFaces[edgeId][1] = rightElementId; + } + else + { + mesh.m_edgesFaces[edgeId][1] = constants::missing::uintValue; + mesh.m_edgesNumFaces[edgeId] = 1; + } + + break; + } + } + + if (otherEdgeId != constants::missing::uintValue) + { + mesh.m_facesEdges[leftElementId][j] = otherEdgeId; + + if (!CleanUpEdge(mesh, edgeId)) + { + return false; + } + } + + otherEdgeId = edgeId; + } + + // deactivate cell + mesh.m_numFacesNodes[connectedElementId] = 0; + return true; +} + +void meshkernel::CasulliDeRefinement::UpdateDirectlyConnectedNonTriangleElements(Mesh2D& mesh, + const UInt index, + const UInt elementId, + const UInt connectedElementId) const +{ + // polygons of degree higher than three: remove node and link + + for (UInt j = 0; j < mesh.m_numFacesNodes[connectedElementId]; ++j) + { + UInt edgeId = mesh.m_facesEdges[connectedElementId][j]; + + if (mesh.m_edgesNumFaces[edgeId] < 2) + { + continue; + } + + if (mesh.m_edgesFaces[edgeId][0] == elementId || mesh.m_edgesFaces[edgeId][1] == elementId) + { + UInt ndum = mesh.m_numFacesNodes[connectedElementId] - 1; + + std::shift_left(mesh.m_facesEdges[connectedElementId].begin() + j, mesh.m_facesEdges[connectedElementId].begin() + ndum, 1); + + // remove one node per removed link + // take the first node that has not been removed before, but not the node that is kept, + // which is the first of the center cell + + UInt i = 0; + + while (i < mesh.m_numFacesNodes[connectedElementId] && + mesh.m_facesNodes[connectedElementId][i] != mesh.GetEdge(edgeId).first && + mesh.m_facesNodes[connectedElementId][i] != mesh.GetEdge(edgeId).second && + mesh.m_facesNodes[connectedElementId][i] != mesh.m_facesNodes[elementId][0]) + { + ++i; + } + + if (index < mesh.m_numFacesNodes[connectedElementId]) + { + std::shift_left(mesh.m_facesNodes[connectedElementId].begin() + i, mesh.m_facesNodes[connectedElementId].begin() + ndum, 1); + } + else + { + throw AlgorithmError("No node found in Casulli de-refinement"); + } + + mesh.m_numFacesNodes[connectedElementId] = ndum; + } + } +} + +bool meshkernel::CasulliDeRefinement::UpdateDirectlyConnectedElements(Mesh2D& mesh, + const UInt elementId, + const std::vector& directlyConnected, + const std::vector>& kne) const +{ + for (UInt k = 0; k < directlyConnected.size(); ++k) + { + UInt connectedElementId = directlyConnected[k]; + + if (mesh.m_numFacesNodes[connectedElementId] < constants::geometric::numNodesInQuadrilateral) + { + if (!UpdateDirectlyConnectedTriangleElements(mesh, k, connectedElementId, kne)) + { + return false; + } + } + else + { + UpdateDirectlyConnectedNonTriangleElements(mesh, k, elementId, connectedElementId); + } + } + + return true; +} + +int meshkernel::CasulliDeRefinement::GetNodeCode(const Mesh2D& mesh, + const std::vector& nodeTypes, + const UInt elementId) const +{ + int nodeCode = std::numeric_limits::lowest(); + + for (UInt i = 0; i < mesh.m_numFacesNodes[elementId]; ++i) + { + nodeCode = std::max(nodeCode, nodeTypes[mesh.m_facesNodes[elementId][i]]); + } + + return nodeCode; +} + +void meshkernel::CasulliDeRefinement::RedirectNodesOfConnectedElements(Mesh2D& mesh, const UInt elementId, const UInt nodeId, const std::vector& indirectlyConnected) const +{ + for (UInt i = 0; i < indirectlyConnected.size(); ++i) + { + UInt connectedElementId = indirectlyConnected[i]; + + if (mesh.m_numFacesNodes[connectedElementId] < 3) + { + mesh.m_numFacesNodes[connectedElementId] = 0; + } + + for (UInt j = 0; j < mesh.m_numFacesNodes[connectedElementId]; ++j) + { + // Another node id of the element + UInt faceNodeId = mesh.m_facesNodes[connectedElementId][j]; + + for (UInt k = 1; k < mesh.m_numFacesNodes[elementId]; ++k) + { + if (faceNodeId == mesh.m_facesNodes[elementId][k]) + { + mesh.m_facesNodes[elementId][j] = nodeId; + } + } + } + } +} + +bool meshkernel::CasulliDeRefinement::RemoveUnwantedBoundaryNodes(Mesh2D& mesh, + const std::vector& nodeTypes, + const Polygons& polygon, + const std::vector& indirectlyConnected) const +{ + for (UInt kk = 0; kk < indirectlyConnected.size(); ++kk) + { + UInt connectedElementId = indirectlyConnected[kk]; + bool continueOuterLoop = false; + + if (mesh.m_numFacesNodes[connectedElementId] < 3) + { + mesh.m_numFacesNodes[connectedElementId] = 0; + } + + for (UInt i = 0; i < mesh.m_numFacesNodes[connectedElementId]; ++i) + { + UInt edgeId = mesh.m_facesEdges[connectedElementId][i]; + + if (mesh.m_edgesNumFaces[edgeId] == 1) + { + UInt im1 = NextCircularBackwardIndex(i, static_cast(mesh.m_facesNodes[connectedElementId].size())); + UInt previousEdgeId = mesh.m_facesEdges[connectedElementId][im1]; + + if (mesh.m_edgesNumFaces[previousEdgeId] == 1) + { + UInt nodeId = mesh.FindCommonNode(edgeId, previousEdgeId); + + // [sic] weird + if (nodeId == constants::missing::uintValue) + { + continue; + } + + // this node may be outside polygon: ignore + if (nodeTypes[nodeId] != 3 && polygon.IsPointInAnyPolygon(mesh.Node(nodeId))) + { + if (mesh.m_numFacesNodes[connectedElementId] > 3) + { + std::shift_left(mesh.m_facesNodes[connectedElementId].begin() + i, mesh.m_facesNodes[connectedElementId].end(), 1); + std::shift_left(mesh.m_facesEdges[connectedElementId].begin() + i, mesh.m_facesEdges[connectedElementId].end(), 1); + + --mesh.m_numFacesNodes[connectedElementId]; + + // redirect node of the link that is kept + if (mesh.GetEdge(previousEdgeId).first == nodeId) + { + mesh.GetEdge(previousEdgeId).first = mesh.GetEdge(edgeId).first + mesh.GetEdge(edgeId).second - nodeId; + } + else + { + mesh.GetEdge(previousEdgeId).second = mesh.GetEdge(edgeId).first + mesh.GetEdge(edgeId).second - nodeId; + } + + // delete other link + + if (CleanUpEdge(mesh, edgeId)) + { + return false; + } + + mesh.SetNode(nodeId, {constants::missing::doubleValue, constants::missing::doubleValue}); + continueOuterLoop = true; + break; + } + else + { + mesh.m_numFacesNodes[connectedElementId] = 0; + + if (!CleanUpEdge(mesh, edgeId) || !CleanUpEdge(mesh, previousEdgeId)) + { + return false; + } + + // previous-previous edgeId (not a real word) + UInt antePreviousEdgeId = std::accumulate(mesh.m_facesEdges[connectedElementId].begin(), mesh.m_facesEdges[connectedElementId].begin() + 3, 0) - edgeId - previousEdgeId; + + if (mesh.m_edgesNumFaces[antePreviousEdgeId] > 1) + { + if (mesh.m_edgesFaces[antePreviousEdgeId][0] == connectedElementId) + { + mesh.m_edgesNumFaces[antePreviousEdgeId] = 1; + mesh.m_edgesFaces[antePreviousEdgeId][0] = mesh.m_edgesFaces[antePreviousEdgeId][1]; + } + else if (mesh.m_edgesFaces[antePreviousEdgeId][1] == connectedElementId) + { + mesh.m_edgesNumFaces[antePreviousEdgeId] = 1; + } + } + + continueOuterLoop = true; + break; + } + } + } + } + } + + if (continueOuterLoop) + { + continue; + } + } + + return true; +} + +bool meshkernel::CasulliDeRefinement::DeleteElement(Mesh2D& mesh, + std::vector& nodeTypes, + const Polygons& polygon, + const UInt elementId, + const std::vector& directlyConnected, + const std::vector& indirectlyConnected, + const std::vector>& kne) const +{ + if (directlyConnected.empty() || indirectlyConnected.empty()) + { + return true; + } + + if (ElementCannotBeDeleted(mesh, nodeTypes, polygon, elementId)) + { + return true; + } + + Point newNode(ComputeNewNodeCoordinates(mesh, nodeTypes, elementId)); + + for (UInt i = 0; i < mesh.m_numFacesNodes[elementId]; ++i) + { + mesh.SetNode(mesh.m_facesNodes[elementId][i], newNode); + } + + UInt N = mesh.m_numFacesNodes[elementId]; + + if (!UpdateDirectlyConnectedElements(mesh, elementId, directlyConnected, kne)) + { + return false; + } + + //-------------------------------- + // Set the node code + + UInt nodeId = mesh.m_facesNodes[elementId][0]; + nodeTypes[nodeId] = GetNodeCode(mesh, nodeTypes, elementId); + + // merge nodes + mesh.SetNode(nodeId, newNode); + + for (UInt kk = 1; kk < N; ++kk) + { + [[maybe_unused]] auto undoAction = mesh.MergeTwoNodes(mesh.m_facesNodes[elementId][kk], nodeId); + } + + // redirect nodes of indirectly connected cells, deactivate polygons of degree smaller than three + RedirectNodesOfConnectedElements(mesh, elementId, nodeId, indirectlyConnected); + + // remove unwanted boundary node: a non-corner node that is shared by two boundary links + if (!RemoveUnwantedBoundaryNodes(mesh, nodeTypes, polygon, indirectlyConnected)) + { + return false; + } + // redirect nodes of directly connected cells and deactivate polygons of degree smaller than three + RedirectNodesOfConnectedElements(mesh, elementId, nodeId, directlyConnected); + + // deactivate links + for (UInt i = 0; i < mesh.m_numFacesNodes[elementId]; ++i) + { + UInt L = mesh.m_facesEdges[elementId][i]; + + if (!CleanUpEdge(mesh, L)) + { + return false; + } + } + + // deactivate cell + mesh.m_numFacesNodes[elementId] = 0; + return true; +} + +bool meshkernel::CasulliDeRefinement::CleanUpEdge(Mesh2D& mesh, const UInt edgeId) const +{ + + for (UInt i = 0; i < 2; ++i) + { + UInt nodeId = EdgeNodeIndex(mesh.GetEdge(edgeId), i); + + if (nodeId == constants::missing::uintValue) + { + // already cleaned up + continue; + } + + UInt startOffset = constants::missing::uintValue; + + for (UInt j = 0; j < mesh.m_nodesNumEdges[nodeId]; ++j) + { + if (mesh.m_nodesEdges[nodeId][j] == edgeId) + { + startOffset = j; + break; + } + } + + if (startOffset != constants::missing::uintValue) + { + std::shift_left(mesh.m_nodesEdges[nodeId].begin() + startOffset, mesh.m_nodesEdges[nodeId].end(), 1); + } + else + { + // Error: no link found + return false; + } + + --mesh.m_nodesNumEdges[nodeId]; + } + + mesh.GetEdge(edgeId).first = constants::missing::uintValue; + mesh.GetEdge(edgeId).second = constants::missing::uintValue; + return true; +} + +std::vector meshkernel::CasulliDeRefinement::ComputeNodeTypes(const Mesh2D& mesh, const Polygons& polygon) const +{ + std::vector nodeTypes(mesh.GetNumNodes(), 0); + + for (UInt i = 0; i < mesh.GetNumNodes(); ++i) + { + if (polygon.IsPointInAnyPolygon(mesh.Node(i))) + { + nodeTypes[i] = mesh.m_nodesTypes[i]; + } + } + + return nodeTypes; +} + +std::vector meshkernel::CasulliDeRefinement::ElementsToDelete(const Mesh2D& mesh, const Polygons& polygon) const +{ + std::vector nodeTypes(ComputeNodeTypes(mesh, polygon)); + std::vector cellMask(InitialiseElementType(mesh, nodeTypes)); + std::vector elementCentres; + elementCentres.reserve(cellMask.size()); + + for (UInt k = 0; k < cellMask.size(); ++k) + { + if (cellMask[k] == ElementType::WasNodeAfter && mesh.m_numFacesNodes[k] > 0) + { + bool toDelete = false; + + for (UInt j = 0; j < mesh.m_numFacesNodes[k]; ++j) + { + if (nodeTypes[mesh.m_facesNodes[k][j]] > 0) + { + toDelete = true; + break; + } + } + + if (toDelete) + { + elementCentres.push_back(mesh.m_facesMassCenters[k]); + } + } + } + + return elementCentres; +} diff --git a/libs/MeshKernel/src/Mesh.cpp b/libs/MeshKernel/src/Mesh.cpp index c92d2de83..d69864ad0 100644 --- a/libs/MeshKernel/src/Mesh.cpp +++ b/libs/MeshKernel/src/Mesh.cpp @@ -507,6 +507,16 @@ std::unique_ptr Mesh::ResetNode(const UInt nodeId, return undoAction; } +void Mesh::SetNode(const UInt nodeId, const Point& newValue) +{ + if (nodeId >= GetNumNodes()) + { + throw ConstraintError("The node index, {}, is not in range.", nodeId); + } + + m_nodes[nodeId] = newValue; +} + std::unique_ptr Mesh::DeleteEdge(UInt edge) { if (edge == constants::missing::uintValue) [[unlikely]] @@ -1155,6 +1165,12 @@ void Mesh::CommitAction(MeshConversionAction& undoAction) undoAction.Swap(m_nodes, m_projection); } +void Mesh::CommitAction(FullUnstructuredGridUndo& undoAction) +{ + undoAction.Swap(m_nodes, m_edges); + Administrate(); +} + void Mesh::RestoreAction(const AddNodeAction& undoAction) { m_nodes[undoAction.NodeId()] = Point(constants::missing::doubleValue, constants::missing::doubleValue); @@ -1202,3 +1218,9 @@ void Mesh::RestoreAction(MeshConversionAction& undoAction) { undoAction.Swap(m_nodes, m_projection); } + +void Mesh::RestoreAction(FullUnstructuredGridUndo& undoAction) +{ + undoAction.Swap(m_nodes, m_edges); + Administrate(); +} diff --git a/libs/MeshKernel/src/MeshRefinement.cpp b/libs/MeshKernel/src/MeshRefinement.cpp index 712bb609b..5a1f08d38 100644 --- a/libs/MeshKernel/src/MeshRefinement.cpp +++ b/libs/MeshKernel/src/MeshRefinement.cpp @@ -824,7 +824,8 @@ void MeshRefinement::FindHangingNodes(UInt face) std::fill(m_isHangingNodeCache.begin(), m_isHangingNodeCache.end(), false); std::fill(m_isHangingEdgeCache.begin(), m_isHangingEdgeCache.end(), false); - auto kknod = numFaceNodes; + auto kknod = numFaceNodes - 1; + for (UInt n = 0; n < numFaceNodes; n++) { const auto edgeIndex = m_mesh.m_facesEdges[face][n]; diff --git a/libs/MeshKernel/src/Operations.cpp b/libs/MeshKernel/src/Operations.cpp index 05dabd740..3b22cc257 100644 --- a/libs/MeshKernel/src/Operations.cpp +++ b/libs/MeshKernel/src/Operations.cpp @@ -115,21 +115,32 @@ namespace meshkernel UInt NextCircularForwardIndex(UInt currentIndex, UInt size) { - UInt index = currentIndex + 1; - if (index >= size) + if (size == 0) { - index = index - size; + throw ConstraintError("Invalid rotation range "); } - return index; + + if (currentIndex >= size) + { + throw ConstraintError("Index is out of range: {} not in [0 .. {}]", currentIndex, size - 1); + } + + return currentIndex == size - 1 ? 0 : currentIndex + 1; } UInt NextCircularBackwardIndex(UInt currentIndex, UInt size) { - if (currentIndex == 0) + if (size == 0) { - return currentIndex + size - 1; + throw ConstraintError("Invalid rotation range "); } - return currentIndex - 1; + + if (currentIndex >= size) + { + throw ConstraintError("Index is out of range: {} not in [0 .. {}]", currentIndex, size - 1); + } + + return currentIndex == 0 ? size - 1 : currentIndex - 1; } bool IsPointOnPole(const Point& point) diff --git a/libs/MeshKernel/src/Polygons.cpp b/libs/MeshKernel/src/Polygons.cpp index 3cd9de4b1..2021c372e 100644 --- a/libs/MeshKernel/src/Polygons.cpp +++ b/libs/MeshKernel/src/Polygons.cpp @@ -325,6 +325,33 @@ std::tuple Polygons::IsPointInPolygons(const Point& poin return {false, constants::missing::uintValue}; } +bool Polygons::IsPointInAnyPolygon(const Point& point) const +{ + // empty polygon means everything is included + if (IsEmpty()) + { + return true; + } + + for (UInt i = 0; i < m_enclosures.size(); ++i) + { + PolygonalEnclosure::Region containingRegion = m_enclosures[i].ContainsRegion(point); + + if (containingRegion == PolygonalEnclosure::Region::Exterior) + { + // Point was found in + return true; + } + else if (containingRegion == PolygonalEnclosure::Region::Interior) + { + // Point can be found in an hole in the polygon + break; + } + } + + return false; +} + std::vector Polygons::PointsInPolygons(const std::vector& points) const { std::vector result(points.size(), false); diff --git a/libs/MeshKernel/src/UndoActions/FullUnstructuredGridUndo.cpp b/libs/MeshKernel/src/UndoActions/FullUnstructuredGridUndo.cpp new file mode 100644 index 000000000..1b67c6dad --- /dev/null +++ b/libs/MeshKernel/src/UndoActions/FullUnstructuredGridUndo.cpp @@ -0,0 +1,20 @@ +#include "MeshKernel/UndoActions/FullUnstructuredGridUndo.hpp" +#include "MeshKernel/Mesh.hpp" + +std::unique_ptr meshkernel::FullUnstructuredGridUndo::Create(Mesh& mesh) +{ + return std::make_unique(mesh); +} + +meshkernel::FullUnstructuredGridUndo::FullUnstructuredGridUndo(Mesh& mesh) : BaseMeshUndoAction(mesh), m_savedNodes(mesh.Nodes()), m_savedEdges(mesh.Edges()) {} + +void meshkernel::FullUnstructuredGridUndo::Swap(std::vector& nodes, std::vector& edges) +{ + std::swap(nodes, m_savedNodes); + std::swap(edges, m_savedEdges); +} + +std::uint64_t meshkernel::FullUnstructuredGridUndo::MemorySize() const +{ + return sizeof(FullUnstructuredGridUndo) + sizeof(Point) * m_savedNodes.capacity() + sizeof(Edge) * m_savedEdges.capacity(); +} diff --git a/libs/MeshKernel/tests/src/MeshRefinementTests.cpp b/libs/MeshKernel/tests/src/MeshRefinementTests.cpp index c4620668a..8dc4f17c9 100644 --- a/libs/MeshKernel/tests/src/MeshRefinementTests.cpp +++ b/libs/MeshKernel/tests/src/MeshRefinementTests.cpp @@ -1,4 +1,32 @@ +//---- GPL --------------------------------------------------------------------- +// +// Copyright (C) Stichting Deltares, 2011-2024. +// +// 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 version 3. +// +// 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 . +// +// contact: delft3d.support@deltares.nl +// Stichting Deltares +// P.O. Box 177 +// 2600 MH Delft, The Netherlands +// +// All indications and logos of, and references to, "Delft3D" and "Deltares" +// are registered trademarks of Stichting Deltares, and remain the property of +// Stichting Deltares. All rights reserved. +// +//------------------------------------------------------------------------------ + #include "MeshKernel/BilinearInterpolationOnGriddedSamples.hpp" +#include "MeshKernel/CasulliDeRefinement.hpp" #include "MeshKernel/CasulliRefinement.hpp" #include "MeshKernel/SamplesHessianCalculator.hpp" @@ -12,6 +40,7 @@ #include "MeshKernel/Polygons.hpp" #include "TestUtils/Definitions.hpp" #include "TestUtils/MakeMeshes.hpp" +#include "TestUtils/MeshReaders.hpp" #include "TestUtils/SampleFileReader.hpp" #include "TestUtils/SampleGenerator.hpp" @@ -21,6 +50,8 @@ using namespace meshkernel; +namespace mk = meshkernel; + TEST(MeshRefinement, MeshRefinementRefinementLevels_OnFourByFourWithFourSamples_ShouldRefinemesh) { auto mesh = MakeRectangularMeshForTesting(5, 5, 10.0, Projection::cartesian); @@ -1847,3 +1878,340 @@ TEST(MeshRefinement, CasulliPatchRefinement) } } } + +void TestDerefinedMesh(const mk::UInt nx, const mk::UInt ny, const std::string& interactorFileName) +{ + constexpr double tolerance = 1.0e-12; + + auto interactorMesh = ReadLegacyMesh2DFromFile(interactorFileName); + + auto curviMesh = MakeRectangularMeshForTesting(nx, ny, 10.0, Projection::cartesian, {0.0, 0.0}, + true /*ewIndexIncreasing*/, + true /*nsIndexIncreasing*/); + Mesh2D mesh(curviMesh->Edges(), curviMesh->Nodes(), Projection::cartesian); + mesh.Administrate(); + + const std::vector originalNodes(mesh.Nodes()); + const std::vector originalEdges(mesh.Edges()); + + meshkernel::CasulliDeRefinement casulliDerefinement; + + auto undoAction = casulliDerefinement.Compute(mesh); + + const std::vector refinedNodes(mesh.Nodes()); + const std::vector refinedEdges(mesh.Edges()); + + //-------------------------------- + // Now compare de-refined mesh with one produced by interactor. + + ASSERT_EQ(mesh.GetNumNodes(), interactorMesh->GetNumNodes()); + ASSERT_EQ(mesh.GetNumEdges(), interactorMesh->GetNumEdges()); + + for (UInt i = 0u; i < mesh.GetNumNodes(); ++i) + { + EXPECT_NEAR(interactorMesh->Node(i).x, mesh.Node(i).x, tolerance); + EXPECT_NEAR(interactorMesh->Node(i).y, mesh.Node(i).y, tolerance); + } + + for (UInt i = 0u; i < mesh.GetNumEdges(); ++i) + { + EXPECT_EQ(interactorMesh->GetEdge(i).first, mesh.GetEdge(i).first); + EXPECT_EQ(interactorMesh->GetEdge(i).second, mesh.GetEdge(i).second); + } + + //-------------------------------- + // Now test undo + undoAction->Restore(); + + ASSERT_EQ(originalNodes.size(), mesh.Nodes().size()); + ASSERT_EQ(originalEdges.size(), mesh.Edges().size()); + + for (mk::UInt i = 0; i < originalNodes.size(); ++i) + { + EXPECT_NEAR(originalNodes[i].x, mesh.Node(i).x, tolerance); + EXPECT_NEAR(originalNodes[i].y, mesh.Node(i).y, tolerance); + } + + for (mk::UInt i = 0; i < originalEdges.size(); ++i) + { + EXPECT_EQ(originalEdges[i].first, mesh.GetEdge(i).first); + EXPECT_EQ(originalEdges[i].second, mesh.GetEdge(i).second); + } + + //-------------------------------- + // Now test redo + undoAction->Commit(); + + ASSERT_EQ(refinedNodes.size(), mesh.Nodes().size()); + ASSERT_EQ(refinedEdges.size(), mesh.Edges().size()); + + for (mk::UInt i = 0; i < refinedNodes.size(); ++i) + { + EXPECT_NEAR(refinedNodes[i].x, mesh.Node(i).x, tolerance); + EXPECT_NEAR(refinedNodes[i].y, mesh.Node(i).y, tolerance); + } + + for (mk::UInt i = 0; i < refinedEdges.size(); ++i) + { + EXPECT_EQ(refinedEdges[i].first, mesh.GetEdge(i).first); + EXPECT_EQ(refinedEdges[i].second, mesh.GetEdge(i).second); + } +} + +TEST(MeshRefinement, CasulliDeRefinement) +{ + // de-refine the entire mesh (of different sizes) then compare results with precomputed data. + + const std::string prefix(TEST_FOLDER + "/data/CasulliRefinement/"); + TestDerefinedMesh(21, 21, prefix + "casulli_deref_21_21.nc"); + TestDerefinedMesh(21, 22, prefix + "casulli_deref_21_22.nc"); + TestDerefinedMesh(22, 21, prefix + "casulli_deref_22_21.nc"); + TestDerefinedMesh(22, 22, prefix + "casulli_deref_22_22.nc"); +} + +TEST(MeshRefinement, CasulliDeRefinementPolygon) +{ + // de-refine the mesh inside a polygon then compare results with precomputed data. + + auto interactorMesh = ReadLegacyMesh2DFromFile(TEST_FOLDER + "/data/CasulliRefinement/casulli_derefine_polygon.nc"); + + std::vector points{{55.0, 55.0}, {155.0, 105.0}, {175.0, 225.0}, {25.0, 274.0}, {55.0, 55.0}}; + meshkernel::Polygons polygon(points, Projection::cartesian); + + auto curviMesh = MakeRectangularMeshForTesting(20, 31, 10.0, Projection::cartesian, {0.0, 0.0}, + true /*ewIndexIncreasing*/, + true /*nsIndexIncreasing*/); + Mesh2D mesh(curviMesh->Edges(), curviMesh->Nodes(), Projection::cartesian); + mesh.Administrate(); + + meshkernel::CasulliDeRefinement casulliDerefinement; + + auto undoAction = casulliDerefinement.Compute(mesh, polygon); + + //-------------------------------- + // Now compare de-refined mesh with one produced by interactor. + + ASSERT_EQ(mesh.GetNumNodes(), interactorMesh->GetNumNodes()); + ASSERT_EQ(mesh.GetNumEdges(), interactorMesh->GetNumEdges()); + + constexpr double tolerance = 1.0e-10; + + for (UInt i = 0u; i < mesh.GetNumNodes(); ++i) + { + EXPECT_NEAR(interactorMesh->Node(i).x, mesh.Node(i).x, tolerance); + EXPECT_NEAR(interactorMesh->Node(i).y, mesh.Node(i).y, tolerance); + } + + for (UInt i = 0u; i < mesh.GetNumEdges(); ++i) + { + EXPECT_EQ(interactorMesh->GetEdge(i).first, mesh.GetEdge(i).first); + EXPECT_EQ(interactorMesh->GetEdge(i).second, mesh.GetEdge(i).second); + } +} + +TEST(MeshRefinement, CasulliDeRefinementPolygonThenAll) +{ + // Step 1 de-refine the mesh inside a polygon + // Step 2 de-refine the entire mesh + // compare results with precomputed data. + + auto interactorMesh = ReadLegacyMesh2DFromFile(TEST_FOLDER + "/data/CasulliRefinement/casulli_polygon_then_all.nc"); + + std::vector points{{55.0, 55.0}, {155.0, 105.0}, {175.0, 225.0}, {25.0, 274.0}, {55.0, 55.0}}; + meshkernel::Polygons polygon(points, Projection::cartesian); + + auto curviMesh = MakeRectangularMeshForTesting(20, 31, 10.0, Projection::cartesian, {0.0, 0.0}, + true /*ewIndexIncreasing*/, + true /*nsIndexIncreasing*/); + + Mesh2D mesh(curviMesh->Edges(), curviMesh->Nodes(), Projection::cartesian); + mesh.Administrate(); + + meshkernel::CasulliDeRefinement casulliDerefinement; + + // Derefine on polygon + auto undoAction = casulliDerefinement.Compute(mesh, polygon); + + // Derefine on all + undoAction = casulliDerefinement.Compute(mesh); + + //-------------------------------- + // Now compare de-refined mesh with one produced by interactor. + + ASSERT_EQ(mesh.GetNumNodes(), interactorMesh->GetNumNodes()); + ASSERT_EQ(mesh.GetNumEdges(), interactorMesh->GetNumEdges()); + + constexpr double tolerance = 1.0e-10; + + for (UInt i = 0u; i < mesh.GetNumNodes(); ++i) + { + EXPECT_NEAR(interactorMesh->Node(i).x, mesh.Node(i).x, tolerance); + EXPECT_NEAR(interactorMesh->Node(i).y, mesh.Node(i).y, tolerance); + } + + for (UInt i = 0u; i < mesh.GetNumEdges(); ++i) + { + EXPECT_EQ(interactorMesh->GetEdge(i).first, mesh.GetEdge(i).first); + EXPECT_EQ(interactorMesh->GetEdge(i).second, mesh.GetEdge(i).second); + } +} + +TEST(MeshRefinement, CasulliTwoPolygonDeRefinement) +{ + // Step 1 de-refine the mesh inside a polygon + // Step 2 de-refine the mesh inside an overlapping polygon + // compare results with precomputed data. + + constexpr double tolerance = 1.0e-10; + auto interactorMesh = ReadLegacyMesh2DFromFile(TEST_FOLDER + "/data/CasulliRefinement/casulli_dref_two_polygon.nc"); + + // Centre of element to be deleted + std::vector elementCentreX{25.0, 25.0, 25.0, 25.0, 25.0, 25.0, 45.0, + 45.0, 45.0, 46.66666666666666, 65.0, 65.0, 66.6666666666666, + 85.0, 103.3333333333333, 105.0, 125.0, 125.0}; + std::vector elementCentreY{15.0, 35.0, 55.0, 75.0, 95.0, 115.0, 35.0, + 55.0, 75.0, 96.6666666666666, 35.0, 55.0, 76.6666666666666, + 55.0, 76.6666666666666, 55.0, 98.3333333333333, 75.0}; + + std::vector centrePoints{{55.0, 55.0}, {155.0, 105.0}, {175.0, 225.0}, {25.0, 274.0}, {55.0, 55.0}}; + meshkernel::Polygons centrePolygon(centrePoints, Projection::cartesian); + + std::vector lowerPoints{{14.5, 125.0}, {15.5, 14.5}, {95.5, 53.5}, {165.5, 115.5}, {14.5, 125.0}}; + meshkernel::Polygons lowerPolygon(lowerPoints, Projection::cartesian); + + // Ensure that the edges are numbered in the correct order. + auto curviMesh = MakeRectangularMeshForTesting(20, 31, 10.0, Projection::cartesian, {0.0, 0.0}, + true /*ewIndexIncreasing*/, + true /*nsIndexIncreasing*/); + + Mesh2D mesh(curviMesh->Edges(), curviMesh->Nodes(), Projection::cartesian); + mesh.Administrate(); + + const std::vector originalNodes(mesh.Nodes()); + const std::vector originalEdges(mesh.Edges()); + + meshkernel::CasulliDeRefinement casulliDerefinement; + + // Derefine on centre polygon + auto undoCentrePolygon = casulliDerefinement.Compute(mesh, centrePolygon); + + const std::vector centreRefinedNodes(mesh.Nodes()); + const std::vector centreRefinedEdges(mesh.Edges()); + + mesh.ComputeCircumcentersMassCentersAndFaceAreas(true); + + // Get the element centres of the elements to be deleted. + std::vector toDelete(casulliDerefinement.ElementsToDelete(mesh, lowerPolygon)); + + // Compare the elements to be deleted by the lowerPolygon with expected data. + ASSERT_EQ(toDelete.size(), elementCentreX.size()); + + for (size_t i = 0u; i < elementCentreX.size(); ++i) + { + EXPECT_NEAR(elementCentreX[i], toDelete[i].x, tolerance); + EXPECT_NEAR(elementCentreY[i], toDelete[i].y, tolerance); + } + + // Derefine on lower polygon + auto undoLowerPolygon = casulliDerefinement.Compute(mesh, lowerPolygon); + + const std::vector lowerRefinedNodes(mesh.Nodes()); + const std::vector lowerRefinedEdges(mesh.Edges()); + + //-------------------------------- + // Now compare de-refined mesh with one produced earlier + // Originally compared with interactor results. + + ASSERT_EQ(mesh.GetNumNodes(), interactorMesh->GetNumNodes()); + ASSERT_EQ(mesh.GetNumEdges(), interactorMesh->GetNumEdges()); + + for (UInt i = 0u; i < mesh.GetNumNodes(); ++i) + { + EXPECT_NEAR(interactorMesh->Node(i).x, mesh.Node(i).x, tolerance); + EXPECT_NEAR(interactorMesh->Node(i).y, mesh.Node(i).y, tolerance); + } + + for (UInt i = 0u; i < mesh.GetNumEdges(); ++i) + { + EXPECT_EQ(interactorMesh->GetEdge(i).first, mesh.GetEdge(i).first); + EXPECT_EQ(interactorMesh->GetEdge(i).second, mesh.GetEdge(i).second); + } + + //-------------------------------- + // Now test undo + // First undo the de-refinement inside the lower polygon + undoLowerPolygon->Restore(); + + ASSERT_EQ(centreRefinedNodes.size(), mesh.Nodes().size()); + ASSERT_EQ(centreRefinedEdges.size(), mesh.Edges().size()); + + for (mk::UInt i = 0; i < centreRefinedNodes.size(); ++i) + { + EXPECT_NEAR(centreRefinedNodes[i].x, mesh.Node(i).x, tolerance); + EXPECT_NEAR(centreRefinedNodes[i].y, mesh.Node(i).y, tolerance); + } + + for (mk::UInt i = 0; i < centreRefinedEdges.size(); ++i) + { + EXPECT_EQ(centreRefinedEdges[i].first, mesh.GetEdge(i).first); + EXPECT_EQ(centreRefinedEdges[i].second, mesh.GetEdge(i).second); + } + + // Next undo the de-refinement inside the centre polygon + undoCentrePolygon->Restore(); + + ASSERT_EQ(originalNodes.size(), mesh.Nodes().size()); + ASSERT_EQ(originalEdges.size(), mesh.Edges().size()); + + for (mk::UInt i = 0; i < originalNodes.size(); ++i) + { + EXPECT_NEAR(originalNodes[i].x, mesh.Node(i).x, tolerance); + EXPECT_NEAR(originalNodes[i].y, mesh.Node(i).y, tolerance); + } + + for (mk::UInt i = 0; i < originalEdges.size(); ++i) + { + EXPECT_EQ(originalEdges[i].first, mesh.GetEdge(i).first); + EXPECT_EQ(originalEdges[i].second, mesh.GetEdge(i).second); + } + + //-------------------------------- + // Now test redo + + // First redo the de-refinement inside the centre polygon + undoCentrePolygon->Commit(); + + ASSERT_EQ(centreRefinedNodes.size(), mesh.Nodes().size()); + ASSERT_EQ(centreRefinedEdges.size(), mesh.Edges().size()); + + for (mk::UInt i = 0; i < centreRefinedNodes.size(); ++i) + { + EXPECT_NEAR(centreRefinedNodes[i].x, mesh.Node(i).x, tolerance); + EXPECT_NEAR(centreRefinedNodes[i].y, mesh.Node(i).y, tolerance); + } + + for (mk::UInt i = 0; i < centreRefinedEdges.size(); ++i) + { + EXPECT_EQ(centreRefinedEdges[i].first, mesh.GetEdge(i).first); + EXPECT_EQ(centreRefinedEdges[i].second, mesh.GetEdge(i).second); + } + + // Next redo the de-refinement inside the lower polygon + + undoLowerPolygon->Commit(); + + ASSERT_EQ(lowerRefinedNodes.size(), mesh.Nodes().size()); + ASSERT_EQ(lowerRefinedEdges.size(), mesh.Edges().size()); + + for (mk::UInt i = 0; i < lowerRefinedNodes.size(); ++i) + { + EXPECT_NEAR(lowerRefinedNodes[i].x, mesh.Node(i).x, tolerance); + EXPECT_NEAR(lowerRefinedNodes[i].y, mesh.Node(i).y, tolerance); + } + + for (mk::UInt i = 0; i < lowerRefinedEdges.size(); ++i) + { + EXPECT_EQ(lowerRefinedEdges[i].first, mesh.GetEdge(i).first); + EXPECT_EQ(lowerRefinedEdges[i].second, mesh.GetEdge(i).second); + } +} diff --git a/libs/MeshKernelApi/include/MeshKernelApi/MeshKernel.hpp b/libs/MeshKernelApi/include/MeshKernelApi/MeshKernel.hpp index 337d8dc7a..589112497 100644 --- a/libs/MeshKernelApi/include/MeshKernelApi/MeshKernel.hpp +++ b/libs/MeshKernelApi/include/MeshKernelApi/MeshKernel.hpp @@ -867,11 +867,41 @@ namespace meshkernelapi /// @returns Error code MKERNEL_API int mkernel_mesh2d_translate(int meshKernelId, double translationX, double translationY); + /// @brief Get list of elements that will be removed after the Casulli de-refinement algorithm + /// @param[in] meshKernelId The id of the mesh state + /// @param[in,out] elements List of elements to be removed. + /// @returns Error code + MKERNEL_API int mkernel_mesh2d_casulli_derefinement_elements(int meshKernelId, GeometryList& elements); + + /// @brief Get list of elements that will be removed after the Casulli de-refinement algorithm + /// @param[in] meshKernelId The id of the mesh state + /// @param[in] polygonGeometry The input polygon geometry + /// @param[in,out] elements List of elements to be removed. + /// @returns Error code + MKERNEL_API int mkernel_mesh2d_casulli_derefinement_elements_on_polygon(int meshKernelId, const GeometryList& polygonGeometry, GeometryList& elements); + + /// @brief De-refine mesh using the Casulli de-refinement algorithm + /// @param[in] meshKernelId The id of the mesh state + /// @returns Error code + MKERNEL_API int mkernel_mesh2d_casulli_derefinement(int meshKernelId); + + /// @brief De-refine mesh using the Casulli de-refinement algorithm + /// @param[in] meshKernelId The id of the mesh state + /// @param[in] polygons The input polygons + /// @returns Error code + MKERNEL_API int mkernel_mesh2d_casulli_derefinement_on_polygon(int meshKernelId, const GeometryList& polygons); + /// @brief Refine mesh using the Casulli refinement algorithm /// @param[in] meshKernelId The id of the mesh state /// @returns Error code MKERNEL_API int mkernel_mesh2d_casulli_refinement(int meshKernelId); + /// @brief Refine mesh using the Casulli refinement algorithm + /// @param[in] meshKernelId The id of the mesh state + /// @param[in] polygons The input polygons + /// @returns Error code + MKERNEL_API int mkernel_mesh2d_casulli_refinement_on_polygon(int meshKernelId, const GeometryList& polygons); + /// The function modifies the mesh for achieving orthogonality between the edges and the segments connecting the face circumcenters. /// The amount of orthogonality is traded against the mesh smoothing (in this case the equality of face areas). /// The parameter to regulate the amount of orthogonalization is contained in \ref meshkernel::OrthogonalizationParameters::orthogonalization_to_smoothing_factor diff --git a/libs/MeshKernelApi/src/MeshKernel.cpp b/libs/MeshKernelApi/src/MeshKernel.cpp index d7f02b53b..2c15afebd 100644 --- a/libs/MeshKernelApi/src/MeshKernel.cpp +++ b/libs/MeshKernelApi/src/MeshKernel.cpp @@ -35,6 +35,7 @@ #include #include +#include #include #include #include @@ -2394,6 +2395,157 @@ namespace meshkernelapi return lastExitCode; } + MKERNEL_API int mkernel_mesh2d_casulli_refinement_on_polygon(int meshKernelId, const GeometryList& polygons) + { + lastExitCode = meshkernel::ExitCode::Success; + + try + { + if (!meshKernelState.contains(meshKernelId)) + { + throw meshkernel::MeshKernelError("The selected mesh kernel id does not exist."); + } + + // Convert polygon date from GeometryList to Polygons + auto polygonPoints = ConvertGeometryListToPointVector(polygons); + const meshkernel::Polygons meshKernelPolygons(polygonPoints, + meshKernelState[meshKernelId].m_mesh2d->m_projection); + + meshKernelState[meshKernelId].m_undoStack.Add(meshkernel::CasulliRefinement::Compute(*meshKernelState[meshKernelId].m_mesh2d, + meshKernelPolygons)); + } + catch (...) + { + lastExitCode = HandleException(); + } + return lastExitCode; + } + + MKERNEL_API int mkernel_mesh2d_casulli_derefinement_elements(int meshKernelId, GeometryList& elements) + { + lastExitCode = meshkernel::ExitCode::Success; + + try + { + if (!meshKernelState.contains(meshKernelId)) + { + throw meshkernel::MeshKernelError("The selected mesh kernel id does not exist."); + } + + if (elements.coordinates_x == nullptr || elements.coordinates_y == nullptr) + { + throw meshkernel::MeshKernelError("element coordinate list is null."); + } + + meshkernel::CasulliDeRefinement casulliDerefinement; + meshkernel::Polygons polygons; // empty polygonal region, i.e. the entire mesh + + std::vector elementCentres(casulliDerefinement.ElementsToDelete(*meshKernelState[meshKernelId].m_mesh2d, polygons)); + + elements.num_coordinates = static_cast(elementCentres.size()); + + for (size_t i = 0; i < elementCentres.size(); ++i) + { + elements.coordinates_x[i] = elementCentres[i].x; + elements.coordinates_y[i] = elementCentres[i].y; + } + } + catch (...) + { + lastExitCode = HandleException(); + } + return lastExitCode; + } + + MKERNEL_API int mkernel_mesh2d_casulli_derefinement_elements_on_polygon(int meshKernelId, const GeometryList& polygonGeometry, GeometryList& elements) + { + lastExitCode = meshkernel::ExitCode::Success; + + try + { + if (!meshKernelState.contains(meshKernelId)) + { + throw meshkernel::MeshKernelError("The selected mesh kernel id does not exist."); + } + + if (elements.coordinates_x == nullptr || elements.coordinates_y == nullptr) + { + throw meshkernel::MeshKernelError("element coordinate list is null."); + } + + // Convert polygon date from GeometryList to Polygons + auto polygonPoints = ConvertGeometryListToPointVector(polygonGeometry); + const meshkernel::Polygons polygons(polygonPoints, + meshKernelState[meshKernelId].m_mesh2d->m_projection); + + meshkernel::CasulliDeRefinement casulliDerefinement; + std::vector elementCentres(casulliDerefinement.ElementsToDelete(*meshKernelState[meshKernelId].m_mesh2d, polygons)); + + elements.num_coordinates = static_cast(elementCentres.size()); + + for (size_t i = 0; i < elementCentres.size(); ++i) + { + elements.coordinates_x[i] = elementCentres[i].x; + elements.coordinates_y[i] = elementCentres[i].y; + } + } + catch (...) + { + lastExitCode = HandleException(); + } + return lastExitCode; + } + + MKERNEL_API int mkernel_mesh2d_casulli_derefinement(int meshKernelId) + { + lastExitCode = meshkernel::ExitCode::Success; + + try + { + if (!meshKernelState.contains(meshKernelId)) + { + throw meshkernel::MeshKernelError("The selected mesh kernel id does not exist."); + } + + meshkernel::CasulliDeRefinement casulliDerefinement; + + meshKernelState[meshKernelId].m_undoStack.Add(casulliDerefinement.Compute(*meshKernelState[meshKernelId].m_mesh2d)); + } + catch (...) + { + lastExitCode = HandleException(); + } + return lastExitCode; + } + + MKERNEL_API int mkernel_mesh2d_casulli_derefinement_on_polygon(int meshKernelId, const GeometryList& polygons) + { + lastExitCode = meshkernel::ExitCode::Success; + + try + { + if (!meshKernelState.contains(meshKernelId)) + { + throw meshkernel::MeshKernelError("The selected mesh kernel id does not exist."); + } + + // Convert polygon date from GeometryList to Polygons + auto polygonPoints = ConvertGeometryListToPointVector(polygons); + const meshkernel::Polygons meshKernelPolygons(polygonPoints, + meshKernelState[meshKernelId].m_mesh2d->m_projection); + + meshkernel::CasulliDeRefinement casulliDerefinement; + + meshKernelState[meshKernelId].m_undoStack.Add(casulliDerefinement.Compute(*meshKernelState[meshKernelId].m_mesh2d, + meshKernelPolygons)); + } + catch (...) + { + lastExitCode = HandleException(); + } + return lastExitCode; + } + MKERNEL_API int mkernel_splines_snap_to_landboundary(int meshKernelId, const GeometryList& land, GeometryList& splines, diff --git a/libs/MeshKernelApi/tests/src/ApiTest.cpp b/libs/MeshKernelApi/tests/src/ApiTest.cpp index 91703527a..c608ae6d5 100644 --- a/libs/MeshKernelApi/tests/src/ApiTest.cpp +++ b/libs/MeshKernelApi/tests/src/ApiTest.cpp @@ -3401,3 +3401,284 @@ TEST(Mesh2d, GetFacePolygons_OnAValidMesh_ShouldGetFacePolygons) ASSERT_THAT(expectedFacesXCoordinates, ::testing::ContainerEq(facesXCoordinates)); ASSERT_THAT(expectedFacesYCoordinates, ::testing::ContainerEq(facesYCoordinates)); } + +TEST(Mesh2D, CasulliRefinementErrorCases) +{ + // Prepare + int meshKernelId; + const int isGeographic = 0; + auto errorCode = meshkernelapi::mkernel_allocate_state(isGeographic, meshKernelId); + ASSERT_EQ(meshkernel::ExitCode::Success, errorCode); + + errorCode = meshkernelapi::mkernel_mesh2d_casulli_refinement(meshKernelId + 1); + ASSERT_EQ(meshkernel::ExitCode::MeshKernelErrorCode, errorCode); + + errorCode = meshkernelapi::mkernel_mesh2d_casulli_derefinement(meshKernelId + 1); + ASSERT_EQ(meshkernel::ExitCode::MeshKernelErrorCode, errorCode); +} + +TEST(Mesh2D, CasulliRefinementWholeMesh) +{ + // Prepare + int meshKernelId; + const int isGeographic = 0; + meshkernelapi::mkernel_allocate_state(isGeographic, meshKernelId); + + // Set-up new mesh + auto [num_nodes, num_edges, node_x, node_y, edge_nodes] = MakeRectangularMeshForApiTesting(10, 10, 1.0); + + std::vector originalNodesX(node_x); + std::vector originalNodesY(node_y); + std::vector originalEdges(edge_nodes); + + meshkernelapi::Mesh2D mesh2d{}; + mesh2d.num_edges = static_cast(num_edges); + mesh2d.num_nodes = static_cast(num_nodes); + mesh2d.node_x = node_x.data(); + mesh2d.node_y = node_y.data(); + mesh2d.edge_nodes = edge_nodes.data(); + + auto errorCode = mkernel_mesh2d_set(meshKernelId, mesh2d); + ASSERT_EQ(meshkernel::ExitCode::Success, errorCode); + + errorCode = meshkernelapi::mkernel_mesh2d_casulli_refinement(meshKernelId); + ASSERT_EQ(meshkernel::ExitCode::Success, errorCode); + + meshkernelapi::Mesh2D refinedMesh2d{}; + + // Just do a rudimentary check that the number of noeds and edges is greater in the refined mesh. + + errorCode = mkernel_mesh2d_get_dimensions(meshKernelId, refinedMesh2d); + ASSERT_EQ(meshkernel::ExitCode::Success, errorCode); + + EXPECT_GT(refinedMesh2d.num_nodes, mesh2d.num_nodes); + EXPECT_GT(refinedMesh2d.num_edges, mesh2d.num_edges); +} + +TEST(Mesh2D, CasulliDeRefinementWholeMesh) +{ + // Prepare + int meshKernelId; + const int isGeographic = 0; + meshkernelapi::mkernel_allocate_state(isGeographic, meshKernelId); + + // Set-up new mesh + auto [num_nodes, num_edges, node_x, node_y, edge_nodes] = MakeRectangularMeshForApiTesting(20, 20, 1.0); + + std::vector originalNodesX(node_x); + std::vector originalNodesY(node_y); + std::vector originalEdges(edge_nodes); + + std::vector edgeCentresX(num_edges); + std::vector edgeCentresY(num_edges); + std::vector edgeFaces(2 * num_edges); + std::vector faceCentresX(20 * 20); + std::vector faceCentresY(20 * 20); + std::vector faceNodes(20 * 20 * 4); + std::vector faceEdges(20 * 20 * 4); + std::vector nodesPerFace(20 * 20); + + meshkernelapi::Mesh2D mesh2d{}; + mesh2d.num_edges = static_cast(num_edges); + mesh2d.num_nodes = static_cast(num_nodes); + mesh2d.node_x = node_x.data(); + mesh2d.node_y = node_y.data(); + mesh2d.edge_nodes = edge_nodes.data(); + + mesh2d.edge_x = edgeCentresX.data(); + mesh2d.edge_y = edgeCentresY.data(); + mesh2d.edge_faces = edgeFaces.data(); + mesh2d.face_x = faceCentresX.data(); + mesh2d.face_y = faceCentresY.data(); + mesh2d.face_nodes = faceNodes.data(); + mesh2d.face_edges = faceEdges.data(); + mesh2d.nodes_per_face = nodesPerFace.data(); + + auto errorCode = mkernel_mesh2d_set(meshKernelId, mesh2d); + ASSERT_EQ(meshkernel::ExitCode::Success, errorCode); + + errorCode = meshkernelapi::mkernel_mesh2d_casulli_derefinement(meshKernelId); + ASSERT_EQ(meshkernel::ExitCode::Success, errorCode); + + meshkernelapi::Mesh2D derefinedMesh2d{}; + + // Just do a rudimentary check that there are fewer nodes and edges in the de-refined mesh. + + errorCode = mkernel_mesh2d_get_dimensions(meshKernelId, derefinedMesh2d); + ASSERT_EQ(meshkernel::ExitCode::Success, errorCode); + + EXPECT_LT(derefinedMesh2d.num_nodes, mesh2d.num_nodes); + EXPECT_LT(derefinedMesh2d.num_edges, mesh2d.num_edges); + + // Now check undo + + bool didUndo = false; + errorCode = meshkernelapi::mkernel_undo_state(meshKernelId, didUndo); + ASSERT_EQ(meshkernel::ExitCode::Success, errorCode); + EXPECT_TRUE(didUndo); + + errorCode = mkernel_mesh2d_get_dimensions(meshKernelId, mesh2d); + ASSERT_EQ(meshkernel::ExitCode::Success, errorCode); + + errorCode = mkernel_mesh2d_get_data(meshKernelId, mesh2d); + ASSERT_EQ(meshkernel::ExitCode::Success, errorCode); + + ASSERT_EQ(static_cast(originalNodesX.size()), mesh2d.num_nodes); + ASSERT_EQ(static_cast(originalEdges.size()), 2 * mesh2d.num_edges); + + constexpr double tolerance = 1.0e-12; + + for (size_t i = 0; i < originalNodesX.size(); ++i) + { + EXPECT_NEAR(originalNodesX[i], mesh2d.node_x[i], tolerance); + EXPECT_NEAR(originalNodesY[i], mesh2d.node_y[i], tolerance); + } + + for (size_t i = 0; i < static_cast(2 * mesh2d.num_edges); ++i) + { + EXPECT_EQ(originalEdges[i], mesh2d.edge_nodes[i]); + } +} + +TEST(Mesh2D, CasulliDeRefinementElementsWholeMesh) +{ + // Prepare + int meshKernelId; + const int isGeographic = 0; + meshkernelapi::mkernel_allocate_state(isGeographic, meshKernelId); + + const int numElementsX = 10; + const int numElementsY = 10; + const int numElements = numElementsX * numElementsY; + + // Set-up new mesh + auto [num_nodes, num_edges, node_x, node_y, edge_nodes] = MakeRectangularMeshForApiTesting(numElementsX, numElementsY, 1.0); + + std::vector originalNodesX(node_x); + std::vector originalNodesY(node_y); + std::vector originalEdges(edge_nodes); + + std::vector edgeCentresX(num_edges); + std::vector edgeCentresY(num_edges); + std::vector edgeFaces(2 * num_edges); + std::vector faceCentresX(numElements); + std::vector faceCentresY(numElements); + std::vector faceNodes(numElements * 4); + std::vector faceEdges(numElements * 4); + std::vector nodesPerFace(numElements); + + std::vector removedElementCentresX(num_nodes); + std::vector removedElementCentresY(num_nodes); + meshkernelapi::GeometryList elementsToRemove; + + elementsToRemove.coordinates_x = removedElementCentresX.data(); + elementsToRemove.coordinates_y = removedElementCentresY.data(); + + meshkernelapi::Mesh2D mesh2d{}; + mesh2d.num_edges = static_cast(num_edges); + mesh2d.num_nodes = static_cast(num_nodes); + mesh2d.node_x = node_x.data(); + mesh2d.node_y = node_y.data(); + mesh2d.edge_nodes = edge_nodes.data(); + + mesh2d.edge_x = edgeCentresX.data(); + mesh2d.edge_y = edgeCentresY.data(); + mesh2d.edge_faces = edgeFaces.data(); + mesh2d.face_x = faceCentresX.data(); + mesh2d.face_y = faceCentresY.data(); + mesh2d.face_nodes = faceNodes.data(); + mesh2d.face_edges = faceEdges.data(); + mesh2d.nodes_per_face = nodesPerFace.data(); + + auto errorCode = mkernel_mesh2d_set(meshKernelId, mesh2d); + ASSERT_EQ(meshkernel::ExitCode::Success, errorCode); + + errorCode = meshkernelapi::mkernel_mesh2d_casulli_derefinement_elements(meshKernelId, elementsToRemove); + ASSERT_EQ(meshkernel::ExitCode::Success, errorCode); + + const std::vector removedElementCentres{{1.5, 0.5}, {1.5, 2.5}, {1.5, 4.5}, {1.5, 6.5}, {1.5, 8.5}, {3.5, 0.5}, {3.5, 2.5}, {3.5, 4.5}, {3.5, 6.5}, {3.5, 8.5}, {5.5, 0.5}, {5.5, 2.5}, {5.5, 4.5}, {5.5, 6.5}, {5.5, 8.5}, {7.5, 0.5}, {7.5, 2.5}, {7.5, 4.5}, {7.5, 6.5}, {7.5, 8.5}, {9.5, 0.5}, {9.5, 2.5}, {9.5, 4.5}, {9.5, 6.5}, {9.5, 8.5}}; + + ASSERT_EQ(elementsToRemove.num_coordinates, static_cast(removedElementCentres.size())); + constexpr double tolerance = 1.0e-12; + + for (int i = 0; i < elementsToRemove.num_coordinates; ++i) + { + EXPECT_NEAR(removedElementCentres[i].x, elementsToRemove.coordinates_x[i], tolerance); + EXPECT_NEAR(removedElementCentres[i].y, elementsToRemove.coordinates_y[i], tolerance); + } +} + +TEST(Mesh2D, CasulliDeRefinementElementsMeshRegion) +{ + // Prepare + int meshKernelId; + const int isGeographic = 0; + meshkernelapi::mkernel_allocate_state(isGeographic, meshKernelId); + + const int numElementsX = 10; + const int numElementsY = 10; + const int numElements = numElementsX * numElementsY; + + // Set-up new mesh + auto [num_nodes, num_edges, node_x, node_y, edge_nodes] = MakeRectangularMeshForApiTesting(numElementsX, numElementsY, 1.0); + + std::vector originalNodesX(node_x); + std::vector originalNodesY(node_y); + std::vector originalEdges(edge_nodes); + + std::vector edgeCentresX(num_edges); + std::vector edgeCentresY(num_edges); + std::vector edgeFaces(2 * num_edges); + std::vector faceCentresX(numElements); + std::vector faceCentresY(numElements); + std::vector faceNodes(numElements * 4); + std::vector faceEdges(numElements * 4); + std::vector nodesPerFace(numElements); + + std::vector polygonPointsX({2.5, 7.5, 5.5, 2.5}); + std::vector polygonPointsY({2.5, 4.5, 8.5, 2.5}); + meshkernelapi::GeometryList polygon; + polygon.num_coordinates = 4; + polygon.coordinates_x = polygonPointsX.data(); + polygon.coordinates_y = polygonPointsY.data(); + + std::vector removedElementCentresX(num_nodes); + std::vector removedElementCentresY(num_nodes); + meshkernelapi::GeometryList elementsToRemove; + + elementsToRemove.coordinates_x = removedElementCentresX.data(); + elementsToRemove.coordinates_y = removedElementCentresY.data(); + + meshkernelapi::Mesh2D mesh2d{}; + mesh2d.num_edges = static_cast(num_edges); + mesh2d.num_nodes = static_cast(num_nodes); + mesh2d.node_x = node_x.data(); + mesh2d.node_y = node_y.data(); + mesh2d.edge_nodes = edge_nodes.data(); + + mesh2d.edge_x = edgeCentresX.data(); + mesh2d.edge_y = edgeCentresY.data(); + mesh2d.edge_faces = edgeFaces.data(); + mesh2d.face_x = faceCentresX.data(); + mesh2d.face_y = faceCentresY.data(); + mesh2d.face_nodes = faceNodes.data(); + mesh2d.face_edges = faceEdges.data(); + mesh2d.nodes_per_face = nodesPerFace.data(); + + auto errorCode = mkernel_mesh2d_set(meshKernelId, mesh2d); + ASSERT_EQ(meshkernel::ExitCode::Success, errorCode); + + errorCode = meshkernelapi::mkernel_mesh2d_casulli_derefinement_elements_on_polygon(meshKernelId, polygon, elementsToRemove); + ASSERT_EQ(meshkernel::ExitCode::Success, errorCode); + + const std::vector removedElementCentres{{2.5, 2.5}, {4.5, 4.5}, {4.5, 6.5}, {6.5, 4.5}, {6.5, 6.5}}; + + ASSERT_EQ(elementsToRemove.num_coordinates, static_cast(removedElementCentres.size())); + constexpr double tolerance = 1.0e-12; + + for (int i = 0; i < elementsToRemove.num_coordinates; ++i) + { + EXPECT_NEAR(removedElementCentres[i].x, elementsToRemove.coordinates_x[i], tolerance); + EXPECT_NEAR(removedElementCentres[i].y, elementsToRemove.coordinates_y[i], tolerance); + } +} diff --git a/tools/test_utils/CMakeLists.txt b/tools/test_utils/CMakeLists.txt index 07f115e54..b481eae50 100644 --- a/tools/test_utils/CMakeLists.txt +++ b/tools/test_utils/CMakeLists.txt @@ -30,6 +30,7 @@ set( SRC_LIST ${SRC_DIR}/MakeCurvilinearGrids.cpp ${SRC_DIR}/MakeMeshes.cpp + ${SRC_DIR}/MeshReaders.cpp ${SRC_DIR}/SampleFileReader.cpp ${SRC_DIR}/SampleFileWriter.cpp ${SRC_DIR}/SampleGenerator.cpp @@ -41,6 +42,7 @@ set( ${DOMAIN_INC_DIR}/Definitions.hpp ${DOMAIN_INC_DIR}/MakeCurvilinearGrids.hpp ${DOMAIN_INC_DIR}/MakeMeshes.hpp + ${DOMAIN_INC_DIR}/MeshReaders.hpp ${DOMAIN_INC_DIR}/MockUndoAction.hpp ${DOMAIN_INC_DIR}/SampleFileReader.hpp ${DOMAIN_INC_DIR}/SampleFileWriter.hpp diff --git a/tools/test_utils/include/TestUtils/MakeMeshes.hpp b/tools/test_utils/include/TestUtils/MakeMeshes.hpp index 28198fad8..4671353fd 100644 --- a/tools/test_utils/include/TestUtils/MakeMeshes.hpp +++ b/tools/test_utils/include/TestUtils/MakeMeshes.hpp @@ -50,20 +50,41 @@ ComputeEdgesAndNodes( std::filesystem::path const& file_path, meshkernel::Mesh::Type meshType); +/// @brief Generate a regular grid using the unstructured grid class. +/// +/// @param n number of points in x-direction +/// @param m number of points in y-direction +/// @param dim_x grid extent in x-direction +/// @param dim_y grid extent in y-direction +/// @param origin grid origin +/// @param ewIndexIncreasing Increasing (or decreasing) node indices for edges in east-west direction +/// @param nsIndexIncreasing Increasing (or decreasing) node indices for edges in north-south direction std::unique_ptr MakeRectangularMeshForTesting( meshkernel::UInt n, meshkernel::UInt m, double dim_x, double dim_y, meshkernel::Projection projection, - meshkernel::Point const& origin = {0.0, 0.0}); + meshkernel::Point const& origin = {0.0, 0.0}, + const bool ewIndexIncreasing = true, + const bool nsIndexIncreasing = false); +/// @brief Generate a regular grid using the unstructured grid class. +/// +/// @param n number of points in x-direction +/// @param m number of points in y-direction +/// @param delta grid delta (in both directions) +/// @param origin grid origin +/// @param ewIndexIncreasing Increasing (or decreasing) node indices for edges in east-west direction +/// @param nsIndexIncreasing Increasing (or decreasing) node indices for edges in north-south direction std::unique_ptr MakeRectangularMeshForTesting( meshkernel::UInt n, meshkernel::UInt m, double delta, meshkernel::Projection projection, - meshkernel::Point const& origin = {0.0, 0.0}); + meshkernel::Point const& origin = {0.0, 0.0}, + const bool ewIndexIncreasing = true, + const bool nsIndexIncreasing = false); std::unique_ptr MakeRectangularMeshForTestingRand( meshkernel::UInt n, diff --git a/tools/test_utils/include/TestUtils/MeshReaders.hpp b/tools/test_utils/include/TestUtils/MeshReaders.hpp new file mode 100644 index 000000000..0e306ddda --- /dev/null +++ b/tools/test_utils/include/TestUtils/MeshReaders.hpp @@ -0,0 +1,38 @@ +//---- GPL --------------------------------------------------------------------- +// +// Copyright (C) Stichting Deltares, 2011-2024. +// +// 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 version 3. +// +// 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 . +// +// contact: delft3d.support@deltares.nl +// Stichting Deltares +// P.O. Box 177 +// 2600 MH Delft, The Netherlands +// +// All indications and logos of, and references to, "Delft3D" and "Deltares" +// are registered trademarks of Stichting Deltares, and remain the property of +// Stichting Deltares. All rights reserved. +// +//------------------------------------------------------------------------------ + +#pragma once + +#include "MeshKernel/Entities.hpp" +#include "MeshKernel/Point.hpp" + +#include +#include + +std::vector ReadPointFile(std::string const& filePath); + +std::vector ReadEdgeFile(std::string const& filePath); diff --git a/tools/test_utils/src/MakeMeshes.cpp b/tools/test_utils/src/MakeMeshes.cpp index f50d3b62b..d7d424840 100644 --- a/tools/test_utils/src/MakeMeshes.cpp +++ b/tools/test_utils/src/MakeMeshes.cpp @@ -176,7 +176,9 @@ std::unique_ptr MakeRectangularMeshForTesting( double dim_x, double dim_y, meshkernel::Projection projection, - meshkernel::Point const& origin) + meshkernel::Point const& origin, + const bool ewIndexIncreasing, + const bool nsIndexIncreasing) { std::vector> node_indices(n, std::vector(m)); std::vector nodes(n * m); @@ -205,7 +207,15 @@ std::unique_ptr MakeRectangularMeshForTesting( { for (meshkernel::UInt j = 0; j < m; ++j) { - edges[index] = {node_indices[i][j], node_indices[i + 1][j]}; + if (ewIndexIncreasing) + { + edges[index] = {node_indices[i][j], node_indices[i + 1][j]}; + } + else + { + edges[index] = {node_indices[i + 1][j], node_indices[i][j]}; + } + index++; } } @@ -214,7 +224,15 @@ std::unique_ptr MakeRectangularMeshForTesting( { for (meshkernel::UInt j = 0; j < m - 1; ++j) { - edges[index] = {node_indices[i][j + 1], node_indices[i][j]}; + if (nsIndexIncreasing) + { + edges[index] = {node_indices[i][j], node_indices[i][j + 1]}; + } + else + { + edges[index] = {node_indices[i][j + 1], node_indices[i][j]}; + } + index++; } } @@ -228,7 +246,9 @@ std::unique_ptr MakeRectangularMeshForTesting( meshkernel::UInt m, double delta, meshkernel::Projection projection, - meshkernel::Point const& origin) + meshkernel::Point const& origin, + const bool ewIndexIncreasing, + const bool nsIndexIncreasing) { double const dim_x = delta * static_cast(n - 1); double const dim_y = delta * static_cast(m - 1); @@ -238,7 +258,9 @@ std::unique_ptr MakeRectangularMeshForTesting( dim_x, dim_y, projection, - origin); + origin, + ewIndexIncreasing, + nsIndexIncreasing); } std::unique_ptr MakeRectangularMeshForTestingRand( diff --git a/tools/test_utils/src/MeshReaders.cpp b/tools/test_utils/src/MeshReaders.cpp new file mode 100644 index 000000000..86e5691d3 --- /dev/null +++ b/tools/test_utils/src/MeshReaders.cpp @@ -0,0 +1,51 @@ +#include "MeshReaders.hpp" + +#include +#include +#include + +std::vector ReadPointFile(std::string const& filePath) +{ + std::vector points; + + // read point file + std::string line; + std::ifstream infile(filePath.c_str()); + + while (std::getline(infile, line)) + { + std::istringstream iss(line); + double pointX; + double pointY; + + iss >> pointX; + iss >> pointY; + + points.push_back({pointX, pointY}); + } + + return points; +} + +std::vector ReadEdgeFile(std::string const& filePath) +{ + std::vector edges; + + // read edge file + std::string line; + std::ifstream infile(filePath.c_str()); + + while (std::getline(infile, line)) + { + std::istringstream iss(line); + meshkernel::UInt edgeStart; + meshkernel::UInt edgeEnd; + + iss >> edgeStart; + iss >> edgeEnd; + + edges.push_back({edgeStart, edgeEnd}); + } + + return edges; +}