From 18b20aac6be1827ce47b9da1d30cb9e525293bae Mon Sep 17 00:00:00 2001 From: Simon Sticko Date: Sun, 2 Jan 2022 11:27:50 +0100 Subject: [PATCH] Add a CutFEM tutorial as step-85 This step shows how the NonMatching classes: MeshClassifier, FEImmersedSurfaceValues and NonMatching::FEValues are intended to work together. This is done by solving the Poisson equation on a hypercube mesh where the domain is described by a level set function --- doc/doxygen/images/step-85-active-mesh.svg | 377 +++++++++ .../images/step-85-background-mesh.svg | 411 ++++++++++ doc/doxygen/images/step-85-ghost-faces.svg | 239 ++++++ doc/doxygen/images/step-85-level-set.png | Bin 0 -> 30128 bytes doc/doxygen/images/step-85-solution.png | Bin 0 -> 35810 bytes doc/doxygen/references.bib | 53 ++ doc/doxygen/tutorial/tutorial.h.in | 8 +- examples/step-85/CMakeLists.txt | 39 + examples/step-85/doc/builds-on | 1 + examples/step-85/doc/intro.dox | 251 ++++++ examples/step-85/doc/kind | 1 + examples/step-85/doc/results.dox | 20 + examples/step-85/doc/tooltip | 1 + examples/step-85/step-85.cc | 714 ++++++++++++++++++ 14 files changed, 2114 insertions(+), 1 deletion(-) create mode 100644 doc/doxygen/images/step-85-active-mesh.svg create mode 100644 doc/doxygen/images/step-85-background-mesh.svg create mode 100644 doc/doxygen/images/step-85-ghost-faces.svg create mode 100644 doc/doxygen/images/step-85-level-set.png create mode 100644 doc/doxygen/images/step-85-solution.png create mode 100644 examples/step-85/CMakeLists.txt create mode 100644 examples/step-85/doc/builds-on create mode 100644 examples/step-85/doc/intro.dox create mode 100644 examples/step-85/doc/kind create mode 100644 examples/step-85/doc/results.dox create mode 100644 examples/step-85/doc/tooltip create mode 100644 examples/step-85/step-85.cc diff --git a/doc/doxygen/images/step-85-active-mesh.svg b/doc/doxygen/images/step-85-active-mesh.svg new file mode 100644 index 0000000000..05f00cae1a --- /dev/null +++ b/doc/doxygen/images/step-85-active-mesh.svg @@ -0,0 +1,377 @@ + + + + + + + + + + image/svg+xml + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/doc/doxygen/images/step-85-background-mesh.svg b/doc/doxygen/images/step-85-background-mesh.svg new file mode 100644 index 0000000000..8387329a38 --- /dev/null +++ b/doc/doxygen/images/step-85-background-mesh.svg @@ -0,0 +1,411 @@ + + + + + + + + + + image/svg+xml + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/doc/doxygen/images/step-85-ghost-faces.svg b/doc/doxygen/images/step-85-ghost-faces.svg new file mode 100644 index 0000000000..396725c343 --- /dev/null +++ b/doc/doxygen/images/step-85-ghost-faces.svg @@ -0,0 +1,239 @@ + + + + + + + + + + image/svg+xml + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/doc/doxygen/images/step-85-level-set.png b/doc/doxygen/images/step-85-level-set.png new file mode 100644 index 0000000000000000000000000000000000000000..68cb21ea404a6fa30d3e939fe3b0c9492cf4ee12 GIT binary patch literal 30128 zcmeFYRcxG3^ftD3Y{$$Dv14XtX2;BoYv#nvkeJz7GshTXOfmDCIc8>NX4vje?I*QW z)vI3oFS^pJl}4kP=R9+0&dfxqsmP%r6C(ov05th8(i#8&OyYkZA}q88V=X@m0H6TK zOG{{Zo1Av|WYfrb9J~nq^eSsGJXV&(USr_}ic6Z~P+Ts4s@@^;vvy1bz`&w@+9a~P ze9VC9B8rW3;gs`Y9cQ5Ec82S?DGaGn|Iv7vr@#MeQ{auAO!< z6=_S$b-xPa&{n33$p<&FH6JP}steFl6%HNM4mOvRI&Lv~dK@Etbb%wDmR@whECOC? z0eyo+STVhVlNyjAT7%?&WOM!}u zGrD`*+TPxtZpX2Nap3gjgz<76$Gg;eLJ>GY2wfF zXKUIB^sAUNfd1}peLJq~JuH**!~>RgN!lT%SF~P#jrPOX{I}+GO+`%3{%T?Xk*fXH zPD+qYxbmLHz+iB##h~PHQ|DckY~gZd&&}!P7OWby)bL<~k)qKql;c*|_3tXj4%& zyuUHJ(Di;NC-yuAkS)Y&#L9xe*^lgRY`(ol^K})bpHz3f*Nw9Jt?qnkQpuU==f#?z zp4ND+LiD>o74SS=?0YuM81G73ac(=y+IZc6I>_5_ZX$objDQU1y&`cHzB|6-lrm%m zg*@Cat-r>@>Zf+>kX<~EZ@iA%;i<5#Nb8MdQ_C0EOxeDQZ=B-;Z!Y$r&=f(tJLHhZ zk&zLP1!J!@kA;rsxoX+BrB^yNOk4;W~c+ALq%IJDq>FPMG8x7py3)}&tu&vr^7Jtfqu6}g->bKIcVB#=_ zk7m`lI{Lh}@wP_34KC-P9j#LUc5Q6b<@qV?v3z~+Am#rij&8N3p5QwxLtyMaqs*7T z|9_P$bTa7u2D;KOC?=Obq8yrF4bKt1-%d>Pd^+zNs-3jfYHcQS3*{09XvSX8z6QNy zyG>*?4$tRbr}90Stl#8582Ozw?4>4J=~v82jxX)fshRRVAzVOo+lR^E|I^IQ$9^^w zJ-zg1PepdmwQnb!qLn=cSyn_`(0TLnva$?KOwcPjpT@ruTb#GwOxQCS=H6{6dXMt0 zd+)@P37xkbhZ2RUBrB5He*}?<1^NdDz7LRm&0e2Q4fbv`V${L4TB+XyokthY&-}moDSHNH z%H`4#&mC0%QEWy;VPnPt@NLy#5Klb3(} zpPD`1u%dszwFOaF3}OD{2p=M(1YrX-_8DS23<2ZYu{ z7yn>LHHL4Jr%!wMD`dt*9g zwmN}&1=@dI50~T6*urIlAeSck_BD5^OEst15lKsL1kO2d53(-i&u^zWU9xq(HTinY zHK#wKBOX?)CGgyzmD=TRF-@zmYEk3J zchP|Xpuw%q4EWcg3VD>XnH2FFDT%%gF2qXJ*v<<);OX6?nc*z4k__G3uB1|gJXYFZ zq;gMasE{x+7n?mPnhM_UWetgcDv?ha{N{6nqFYg)JQMDTgg0nrcN8<}H(_SCv-^D7 zav-|*fjoD~YeE(eZ%-~3SrI0<=Vx(NE@bt63sl7@P??_LHW6@b8X%rioNcb9hYVP~ zzD^_UdKx8}WBdho7H}YX^!xxFTr+CF2A8yJZ$a$z2_Hx0U7r+vE42;vCc+#<@Nl1k>_ z==*gyR=+_%DCr;KkWa?;cE#4N?=}LT1o3!dYqq6+cRXY(ZUlA+D*A}419S3G(-KyF z&tf9vVkS1jeapqAQ=nbEHfX*EBtW}tssVGNVp9`VpUcxw(-a9FX^-s#V&cdJ&pHl_ zu2Z{g&(?+V8w*yyYU$ZfZg`(D={UbB6#sh-#=XoVAJl~>7R5ln8`{sHg{+N}_BDd>@UO$C+8nnTd@Hmlla$<}AuRuy&%1K)D$ z)@5*f?d8+mOVFPQjljR9H4_xFgYkhk4Qz^O3E`kR7uuEMkIEJ>6y2N#e)l5BXcvxW z8a51KoKRSp8@|02I$vMPIZoCR?m|qMl6!6ZnGinMo-Mx&EkIX}oCr#HSi{%%nWi8} z6F{pqfW-5IopH}q@~aMir9s=nr`Nj;6aU#_oC<~!j8Np0STpz6!zRcI6BL8EhITl6 z8K2xQk7*5@GH&b$16RvhD?m&_=jo>6yW7pFh+6W?cVvf`q)(>r_ZJ=J?*96@@p4){ zoHOQC9zQBQ(DM{7=T;jYu&Up1JW^eQ)PPQI(>SA$>u?9T!9r=F;VvA~UEMgkqCmm{ zznF1jKl1m?rZOmMV!V0YgLpp3#IjhBaw-j|;7UANGL%PP0Q>Is+tVZRc>RO$sn>b| zlrVlRy_JlKCg+YFmM$)rOYrmklf)V~meAgi{rU1u%|@W&P32mWF4KrRv}5XR-?i8K zVG|UD%jb+?U3g;*7b;~ZZ+VUZl@kR{P`#Jfqt61X@HS`|57!ErvGPpb*Ze7-@9n^> zg5Ml@I>U^dA)WKDJ>?0h$O)-_Dp5CBCe7$_RoA|>*zpq-G_{D&DHOM-=k}rW(Z9cK z*?cziu{5W@$;3}c4W2FhY;+yday?WR#>3|$^=!!C zcUZ>ac2Z*S52fLrTz1=H+9;jS;0-}`F#MS-e81N8rsC^L6#N;ST$vmBXm+F~lA zvwiZ8*&#Vfda#>ibJ22%7xZtj#lh7{=Bl)&5#wmsV#DvnVH#%;kQn8vogYpP(Vqq- z#4(DhK?nJ; zG*T-BIfV@A#J3#4CC_Od%`3OT(^a}!C>Kguhye2yzRZnI^(TcVs>Id+*v&jo1tCNb zyZ6rtnR*feEdCE}rjM^e>mpwFe@tI@AiPZ{r*Y+C(2?y`8(m-joCNZT_@L5z+@enH zI`~WDeC@g7J68((O_M;q!ciTKx-L&bD-XCA9!CX>=vM z1eem54`~$?v0vb+AO4ODv!<^vLwcJhqH0}NtBk4{A7q|itUYNDf;hMl+lyST`s@&h zL2QEVp;P>1skdqbbf>cUA8VRkMv?F?KSQ0kwqo7fPQ>~Pf=$KMtqyjHIr<1Pk^ zRX$P=ZFb(4L+NxM$ysZn88+O33jPu2annZ&jORT9yV5(r)may0cbwcCl2gDaFDg6H zy1zr-aTz^BY9bM)Emt>(;{_>h(xxeATp37*vyyfa&TTi?&?MQb=f%XQrMq#$mFG~X zQi|-Vg#QNxy=v+^_1Yt+v%o0E9MS2&qbv5>vm$mkEHM$)g_OTtiIV5rfAjV6P?jYJ zSxb+H%234R;pLseP$XVqm0E?t>sG;7vYPs+lc>Oh`}RYhKESr4kiIlx01p>jdQ7b6 z@Ea|>X6)u8xz9UWM`>C9{M*hs+lKqh-G$(7F$v7=3@K`3HfEgAi4XbPVVG*bF~WEb zww9g;C3K!i(=HA}5tGg&AJ-tv35DX6Q=yQ#Tl@kF+k&Baw7PkL@NK=aY3bPuCNu|e zfk)Wr%~8(E20*6Jgy_XahniFsIaC#2j;aGrrJ#~g84|!^)Mc}%y3hCg7QXW#e;Sfi z&ZX4SYxry^1f2^4Z1ionEak^_-4iX~D6yw0xO5`uS-Lf|ziB#h)XTDgCw7?QKg>$&dnH%z(a2zMfCN?QF0X$Z zz4c&JZ#E* z1uFcW;s{LrlcJvGb2iEJ^vRu9CNZ(;#0@xC4);5FHegs&e?R#Wps& zP7Ez0#ZgUERrfdwBeAuGb(YkI-LE8S-oGvpB=hi{v1bwjqf)%rLM9DJsuNa*>#xuY z;r|&4zhno!jgesloh~kO#kX}nF0DW|y4r4T@?ic=Uc=Fc2Xg}04WsH2THxzby30vz z$eLdnCIcW5^J|l3vSAG#!}BVuFI#*z4s_CJaNX&b9&smS-EH!r+EIbjN200aY`ORj z;QZRDz9RYuJ|-nkd_cTT++6JxCbK{eENCm!`_Nh+eLp9m+5KTB$b4(;gnHfUgaz{2 z#JBnsHVGSSWqodg@0jtK%80-7b{NB@99LTFZvU)E;cY` z0~Di@0WPTxV`FAUZ=D#Q+Fuz?POtucV+_5^;Oo5GoXPdyvg;UbGr7fh8=ZpO^$4rf zlE4n&aQ_#X^*_f}> z_hAA!Ah<2qf|Y0?FNos7d651hByuh7t%3|=1-&bDQOw2h3G%XRZ2EXh#uu>89+H8& zNK843Pa!o>2aomr6=E>-SPQdmF_ns*r4N~^Wpup z%g|6z3CZ~JmMJ&zj0xSORy*9okpXJwFV5bEjh=OuiQDme*GSB!t3{cYB;wCrf4to! z@Oin$0>I%Xq_Q}!dVhkMcB-~r{-hE)R#1I`-#rLu;UW0sd-Z3~sIxwN zzx5%NZ8c!VCJplB@MbXzNN-8okLqirBt_}|u}1fA(sg0Y%|U+zhte8#Iy=KyX#~Jk zKiY5+w4Ida_V%(6!16mz^D{G`;$i(E4YHrMfW9FYUMfK#dKOJ*;xsCAEghQj^8!A*N#sWxH*5f<=1{)% zQ7Waz>Bqay%M=~jbcO|YaQT2{Y&6qC#23eAXLN@jDWO*xmT<6!qn_=`rMUWew?M4E7U2N{AwYHgcw{!Z(KIT5|+H|~B51P#CY|SE}bRP#iyxp-u zBIbG@qH2^R*)1bk{YNxegN{bY+qH||gQvR=?y2&gznoViwE9HTdEU|zm(QXx8xFu* z;#9$5D!1Oy!s+9{7T?|IpODIFM(~6hexqHBl#CWPi4s08P67V9NHGYkNQL#%Vgw@Ky+e z-<)9czml5tP35p8mX>DkgX830a5Vb%ltlhgSwViCB-Wu$vJzV9c73=J^cSjVVZoBB zEO4?rAV!Z5Xq%Mb8q)YjW&R8CZs~rUB1a~V_602^YNDypXrtirPhTyT=1V#;7FcAT7%ocpRBLYmQ?u(Gn(dyRo>H*WI`7JFJlk2T?YHSSN^U-0{RPKH2NI@2vCt|e^$)!faQV!)bp$(FvFv@CQ_0`ISfUaZU3&?ezQ$nM--T1@f1|-i&61YF zf7YKE+h$%Ude9XcUy{?tvtcEtLmjELH1t_`RPyLhL6FuaRvbg06qBlQhF_+N|jLx74KUM_>N`~t6D2aLRP0hMF zxE$qA6vUD(wcu%;B%P|>P^sasd{5@&!QSTK{22Xv+D2ri#0cN|JruIie8Z-+`@s@( ztUzm)#a062_DzkKLiK&o>`l?8d3l?}y>nesn!VjrDJsA}l6X+*&*(jl^<_BDVD0$!qUG(8CGUJA zKnPW$HQ?uOKw1bT?sruW`VS?9j4FTAPP7K;V-4-5DG zb+G_I%G7Pko^X>z=3hjO>k8NXx?F|l%N2ccleh|HeCA;F%r=Y?5X#NubMMdFiD>+h zebC`B8pBpW7)8lCMT*njjU%pgJ)pIhxjikq9zZA=l#JH|^yL`NeJ&h>H^#H!xEFiB zu2_FcTfENx1ccEO@_!7`eOK_Z8%AQ{u#7~k9}hCz1GaWq==>3Km14`1bF{fZCOaa?5m{I zW3@+Rrg6D9pYD2RpQDWA5WCvX1D}8X$SCeUQP((RlEOyJsi0a>@4%79O;D#Brag-C zeq>8+JPVE#pEC*pMqfTyf&N>(3280r*|+^Q2?NJ|REj=xjV3>7Zlr`Mltt0T(;{S4(3Z{>AA=-WPU+ zEiK^Q-dt`{_fv68sDN|pK>_^lq8|Dbt6_Q{a=hhZG+9j3$A1jXvdf2Pa;)$Kc1%AJ zGo4gpcZxJ`xvLs#R!O3m%66H{99#1C#lJWT0Igc!ce;1-ZY#ZJDcZqDc=*8$y%A!! zr#9~Ktlu$1;R}|mH4sM$QM_odBrFjq+_WIRadW&D!(>cn8d{1Mxs()NK_!$BPZj{)5tjFz^xz!&xs4txzY)l`htvO zVg3N`b{9?*Nd1EFzwQF4#g^JrBdAB>4p2rq8;zbCTHMMq@K0kZR{v|DttF{%GuLy0 z+Nym2W~H!=waH@Jh0n$N<#wl=+DDc88bBor+54tj?vjXWdaV10&u=|TvXonRh?LBj z<4swBA2mY&c1v4GQr`O})BD{yx!>v7QXAj$)`!t$vZku{n{kJYHw1CC7FIGpM&Y~T zpA~3VjFP8pHhj8al=_ch(Yeu98?4$q9|1~zq{{Ng6{);YNNje`SKyfic|z?R;3q(V zUsjcjKJ&K@Fw*$4FjIqm)sCd#iQ5f8fvqDuWb#^ABPBs@G@wI0Gt4O_fsWhRhkc|y zkJ&#F&!Iwy6qZttt}UG7a}$l|_NH=eOqY8)=2J`v>b<7E^tJ1Xf#(EB_b_Kpm8-cu zwcT`demufsZlzkRX3y>H`wSa1z$@TrAmMHKhMptNBXsXM78+}lo!}`)zCU^5NOQH_ zp9bz)Jc!&!*^v1^f}cL}F5e`e9+U|=$AnPO&L}C}2~w1cx(+j&5_$|O=rgvYGs4R> zYqWCO0)mbQ1~QVp(hOAeaBj|b$m9nXVRlS=1qU^R%gt89`}CB%(FU<7kjM z6QB{O#M4BhafM<|3P1}vrwF4Tlav1g%VkL&+H4Q1f3~Bo;B1gcc{}t6+ zL3zzc`;;8dq1Ier*#r)kLFqSWQhQclT5J0yJY~Qsp)K=^&uWf)UM$u8F4;uWJ;IG= z4sUImlp<*K4>^sgqoh`II8mwqDl)xL%#zm}z|zO$6>dpSwFcJ&uMj64TrBPYBfXV^ z>v}cCDOIPNqEY}4hynv!gcw8jvMEY-rnXJj{@e>8C3V zcu{9OjWihtZMvvR1}rjy4GaZe^AC)d1X`ux6i~5vP&2|6yFUVnEX>3$_YLH0@M3Gn z5=AX4h#a-u!(ub+Xenm0P19HW)?a;I2qwf9(@kyv+9ch_+KqEBJI&Z;+$&9I$}e+0 zeKV5_@#p@#G*m)mX(vj&?4F3>`&PI?&Px6%mp^lqLY%Y8di(wdC?uOdGBYY;v#gXj zc$+|VPn%%jLcy6bnm0l9iQRrQCb;D=bmRc@fwaFUQijOsB=?ET(bNA!ut$z|&OD<% z@9KfDDB9+u6Cc3!#P^$o%~iNGoI?hXN&-N?i`Owh-;S)DlKVr-hTVSLU_tr8fQaAf zkr>2~pg_Hbz8~F4EDdPwUk3dKi^#6^(0_*YD^APNzWunROdwFfrKG;QoMFte9(06L z?Z4uu#?%=7zdMaf9$k-kQd&wVVM^wpzQ5Qw&mq#gVRP<{ux7C|0&{kqYqw$2(&kK* zma$=;nv9A9@TP<_&b*cauGyx%@`{qwQh-QQNisk=kKIi;!|#1=bJg^p?qy3SAIVlrOEJIIV$}ww0A&>7P|hl|u{~<9rkqBTHUh z1q1S4jEczgYwC#7sW3|IXWpy3DD}*23OUS0oME+2d+)y<43Z0oLr<|W6&Y`lY4p#J z`~0rRnxbp??=u^WA+WT96?6{*F$q8F{Y{&>l+fgb#>O0-O8xf?A?6S(59I2aVrU}c zbaW|a@=g0-NCOiem1D}k0t?Okn-8nlNDt=ch?zM43Uj8!;YylD@cI?Uv%OispV?gpybDnN8ZAmfrH9GuR@1A8FSfX8vw_DJb-{LtMp&p1T zh=n?=F1X^Xtp!gAxOR}|sF2Ae;Dxb`jza)R=k9PJPNxUTsmn79W0^puXrei2}qaU2>mirdpWw2_PkdMk?sGPhRHoPY}VRqaG>4(W}68OujY1Zcd zjK&${W3OK^ek>>yXW-x#kim0>0D>3Ylae(KopWFccA9X)hl&Lk6)w-|?_laRx~*b{_m-$` zxwS{L1W!<;h&vY4x1KK-xas5`m@RzSh|qmE4bG|0y(;lK>$;69{cNcI9}A?NTm6$q zq{hhLANiQ7QBk12FKig0MuDFOR*Gwpj<5}$B4Es-@&b_xYIPWayb&~j!?Mf8B%UcT$7K#YWl7|4JUA89V%GA_|y z?|Mw;Djq?gFAmsc+{$$;F|d0&d>BMgltWEPosX*-*o6g{eQV?CMdf7mT6E+dCbWpn z4~21Xb8u_)Q9r$u;uyI4F4L-T%GT0}bvv?N*GRp{H1#ri z=%tvF3Q$wXT@gRYsBQh_dAOP-E+W!@_RA}PF`jnI0^#2ghSWI!kgJ2@Qbl>aY?<^g z0JysXtpCQoQ)-%GO++~9JJW_fThB$G9i8{|obPkF!xr?Q$cLhCaX9e zRXj&saKjnKDBD_^;@J90QV8Zrs%M6*7=_EJOsZ~vWkkmQ@*Ot23Kiz?cg6{!#xAa= zw@u3ReiuD7pMNPPBSA_Ma2kkoudP#K&qUg zbe7X@*?1x~wV5jN*`yh$cVLl*${4*M+=qb4_%j0mW4*B{eCE{^e(AdNnFY!3D;>VD6C`YKiSNnFUBed(f#kw^6#IYLY zyfU*m#}vY6_YrY@(uR0E#Rk6@J{Dqym%jXyhb^8k(^D8o^9x(+K{V5|{uWv8=uCt< z+5u#|((0uM|)hn^DEe3~-Q36@8wrSV&Fawt6h( z_(lLzMqMc_=qnS186`i|{qI#t(N(Xus>oDe!?{PzT740Dg6)RONnI9Qwg$vOF$TRv ztB*|pWO^U~Emws@DhxHWOP3rUPyMt^jvk&~tQC5+t|+&BGY_Oe1_rG{%pb@5ePPvp z;$%$MXyCGil+|suz9TQ%b%~0tGocQCiWt68z~sW(iS4$; z*V_1C^P02U9FW2U=$*<<*VS-cpVubT?o)S14`OW8gsY}#jj~mNB|`era>e^)>xVQ^ z7@pq<|CK1*j|$lWoNb;EX=atMf5xER74-pj!#A+T8!|a1;toSLw`tRu$qD8iN!W?} z!%-%5C^GXhuyl807j+8KmogfY_3ufD*?2f9?D?`>CAGBRIV770Yl;{N+smXTH4jfxn2_hX22vDv?x<-SrUG=D)qHZX${>=x&9!4BdY`e zOdn~Q>N|wk?%K?cdJqCQh~oy7{GfM}C&nkZ+?p1M(4TWAGSl1`f{9~gr{lAPF&b_E zHEo<=m;;AhE2~t32v7Hqe_`tU=0oIpi;SimI=+v{UwsRNTU5R=YY}q)Kqjr1=#p=y zcb>3X`oKL&Q-cjW>)O52ecQ?_e4mxl!?fu*^~u?zNg^biQ9>o+TCi6$cSQ=uUcgJX zBl%HNgp$EEIr&%O^E;A9ZnOMg$h3j zXTE5uyu=eJLt)$PvGoctrGUbVUdV3{FB^*NLT%-iwI_ggfhYbs{@WD?@xfkGT>Ke3 zCHLetwY;^>zV&ide;AoV&YVQo#?c?VuTfREAB*8FY=z%ly*Lw)m~1Pk23*2YsfVI} zl&J&(024K*saroPor4gq`SkLPg|;i#zd`G8rqC)#ddc-m87h05KypU4|4ts6x8qqa z)9jc3G8xu@)4j?T~>hbDWk^6%X=!Ph{=sd4IZg&UcMY zF1p)J5vOI%!!L$E^H z>AcR&8Ai`AmLbl*KiV?mCZRbp9UmXxW`!#TC8o!yP-|N1PkfB{&Pi}h$NRsa5#BOS zq<%RmWLk5C6>CuUM`bedoEA`BJ;GGG?v+IVP{R3p>O8pH%ev8j6EtZ0G$7X7_)F!# zRHpT|cEdXmYCN$V+m+VS_S{4KXi+so=`9(wGR`9s6wG6d7QnDi<|Cdio+UMSm@f{2 zjX8!w&P>q6$Js#E4kx&H6e(C~vuK@FBypu%0YU^Eg2sJMZ zPRgtLFg1J^*B}1=N_UHy-SB4)$#!>-*ZCA#JaUlC#K^V(rW1`R-OHgP@L|1T!#`jR zS{2GfLJ>bWoJOl|wmw@_hi+?<2~n0gRx%lx`Os0YDOtuEm>OX}H|S+M5R6dZ`ch`* z*+t}Eh9Sm6H#{ZUv15$W+JTSSLanBBr3$GTM3$>H%o~pNJS({wPLpF82p#Hz#ON=l zgN1&;r+16P2Y7wLhZhYTG;O3xNrGt}*08H5m1(DV1u)3n<9GL7XKuu1+t;nc3d_oC zjmTY{ma^&rBiSj8puMS%85bT`{BN%E^}^hSMv?A3Yco^y#zuPCSvhMh76&6E!zx}* z|1Myspy8cbEi~hyuaBxqD^8s5Ai{cIcfOa(9iI1rq7LfR|qf#KCA(Vq&v}aCmfuHY1>)G_arHh-na9?ubXuU5po; z!^{lO!qTxL^KREoUuoy8_v?3oqb+DIjEJsZwL>#6!Y)bikguWsnncw7jUTnblbLy9 zauNsVX{NR zZ#C&qae7SIXo|c zrh4pEX$A@9l9Q6t8t6&8h-|6R7z{JDZPd1U=ccZ0eFk22#FfVj=v`fIS<=^C`lut} z#EI?WX^rr++_L9x2s98+W-~d#ii*&&Wy?;Cw)SnHlz_pB*rJaZcm7$nV7E031T;tg z6p4#(63h^)Er5%fp&hWnit9aGUc(H|*K{>FMd~H}0W2-1jjrW3dIzm|6U~9SM$g+|Ept zDB6$uDORR(p;n{QFaS*Dy1B*DKS-5ANlLZl1-ZCY9_$^s!-cN=fFVc^U53>GU(H<%b&g10pBy)vn=7tEN6CRQ&wV33DK>r znGBQV9M-->F6syF+Ao-n}3-x@-nL<&EL(q1u_VSDaw{ zQZP|Vt<>tTEa&t1WKH7SC(4+}}sZWPxSMoW>#@V8h@dtHLIOng~6udUB0W=$O)NmmJ0)hhq9UYZF9`l<{GT)$Y zarh`M$`V_V6%0;p0nRkm1@->ASU5RzWU+N)XuyvI0v<%0Bb^mPsd9D$W0P-hUZ4C= z@1cYov;0R06<_XmS!^}LmzI>zy80`sT78F zw>dMgV~A?kBxlS^xyQj-g{7o`VdJYBPL2=dIbfAcaE z((nVq95RzXEmj0p`Ms9PBMu~TxCMd`_{gQqEoJ{nPaE#eoDDFKHBkHEQ8N25Po!X9 zW~6K>TT490_UMEIic)3dZ``;5KdJyp^(opE21a`6S%$dUXb~T1gBleB(ya@zc4dA& zI!ou>b6z)s)9^>uZZI*GIcNV`aoG`z9SmMw%U`KCOi%-Bz3Gq(Q5BT?D;T+M;?nqI z& zbUe%CSjdsYNr&6Yi=u2-_U(eKDJd;L)G+d#4oZ3h-oSALsOGzM899pY)?KXY#8OHi zn3=T30nE~e3{W^}qZaLyG5@v*W4J(@%G;R0xwNrZEj+u$L}zxVlhUx@#35-=XNJE0V>gH>%jEq ziTmC>EmVQ2v&)ut(rzp>;5U zc08N|&SHO=g(G>@P%0dn*7lgon^Q>4JNgNOVW+lVNdfvQRd|NqBV}x+;Cse=NfI-% zvKVAPVd{@^)w{Xf#H^V#i5+rHjg7@>g{gY@Ni?(no{EidP7!nvedpYY2L2TwG|>Yx zLDo{y@NRbpIKz-gX!WAzkFc-*_RaSpDh>^F;B@dFvFi-*d8cW!q^{6IpM$os^!1@q zD9_YHC`^kuvUbJnh|qB8YJK-1ly$_VxX&m z)4!ZJ9t!Ipa`$2F8fJU|_{8X+Y%{YHzV#r|JJNs!y{~ruCh~#}fMrbfW z1H?^?ggV|}D4qO)k~tB{B{c!01foW%KU097+wkx)g^#$RH^L%@zsg7*u^Z^8Hn_d} z7RVC|d7Z-^K90DQ^|IG5!6) zCfEI%(2uKP~-x#40?ZFz)z z4rtqUl~?TWxy~?%{qo@hCAY@Zyg}0j~?v4;QP&_n9FtF@v=Q)<&XUN z)&MqQ-RWEPDw3wXh5hS#iwQN(FFeBVWa+PqR)PC(E-w1?7#X4PWAG0MAZbB*9?7;&pqhT z!PXm}oJcRKg}yn+ah8K+icE+k@27xqOvQ9&A&vqVq(OYg=T z0UR#8gwd)sp*BD~g%w!z>l=)VaC&q?)`;N3_go_{Fq;=rc;eo?3J39sMYC&`+>CWk z0(u-OA_BX@n$y@;oEDe<*jrC(IVVZp?DEwFj&!Q~afBJ66+K2?gzXR{t*;$o=`v=k@8?HY@P+ylrI zZxbH-y;s9I^C|RnFRR*gv4g-7=#;vWyE;m2<|(GAZ{4t9V{MHKF2`I#md<(qm-fq? z>4@%J*x_L~C8>8x%^%khzTll6o@i-A0}KkFWhbKT{ROnn7|IMC)>9zl1jhhCRNPIkYFGoN>fxAF!ANXSd^)QxuJ_Kj^*02y?pUMk;Tvm7P0i<)K7a2JIdb)lew$} zbni_96;2&d5cD6t`Vc6r>)5vnYb?XQ}lRIxOho^Mcl&;$)e$+&+ej#ey5zsc~9urqB`vXV z{HN$}gs{7boX`p>9H&DPUNF)Bx(o1@FZ4)h2m>ug33;B223~txvxM+dcfIyR53%8C z>w|E{ME1a%b-X|}xnNgo0j{+tB~2soTCjp(WE^Jd+NYe@((A>4LL{P$_Dnp^0!=3Z zTDjJZyF8xPu5l6buX1PYyh}!376aVtsPq{45a4JHPVTLGx$B+1Bn$^2Sg>H-_T{5> zw~c^pW`wGHgp^(?euw0x?5eKn^^|&oBr)2U--W{DRxPiIs#kpJMBO3eH_a$bUkgO6 z+dPmg&O6VT@_Y|Sga*@=2o;nnmTZTsm($@w;WKi3ypXsYfMeUI%&1{8A+CuT)OJm4 zYIk~iN`KkKMt`lk@{=YNwb)E{77iv@`G^Yfp>%U@TM&N>=)PAX8FTSc)uLg`B^uXZ zGSDvii;4Dw@BNeG(Yd;NpDsc3Y*(}@YQgbFRZXLR^_?|3-j74Gwc5g+Fu2Bf>%cKV zthnOZXe3KlL(B8tmCb?#1lYEGROhiYu`6S*$&N;eYI1{b^?hGdF=b#YR*oL4nwM;Y ztEzhe^Yg9wUk^2O$HuLv6^IFh4<#bb1^(8}!xNqkSjsmuEj6h$MVK2{wwNFK#{A); z7m>IZOBh0Jl;XQoplz+|Rd02B=uB%rT3hiTl?priFSWFg6MfH2`}4_^H%bB&S+i@( z2?}OILiyGI)WVMrLN3hqVjWY^VOmsp|HL{rX*D0e5RS~F*-o7QbSCbj67PqGDmd&yYWXq z*m<2ny(d@07RZ3wFlf8X9ABLA7A`dUTui=64!K@>iU$cHSfcBrIJQ!MEXvPD!JJPp zm%S93aiIV)P%ccRN9&(aG)%qhIY&akc#+Gxz=8jP-v&vXY>)o@ejk8YPs0LLc39sL z1^6B4WXbd?e)Vk<lALHf*3>A=xp}N5GjhkC6}xV{ZU(wOlTn#?-B#NVDD4Qx4&cu3$B^CxBP>( zrFq}0F?O?%@B!5m2aiv`cEcZR_$(L~x>L_YJoBSCbbR9I79jjqZ`DFk!AFSsVLo0% z%(BM-(3{`&lyPR3{QYdg*w5s8&_=c zGZH}?^@^okZ@4n=-?N%h{=ePJbcfb;0=0!xBwKX#h`;`C)m`IzUER|@vC-ISY&1z@ z8;x_~6B|!#+qP|^F&f)WPSU7xz=jOo~3DcluoxpA(Tem zsxV&Ejqe?nG&@cJeZu3>aV<%Vw=%S@itI9gzmN5XyWCiZ(m;J|0!aK9CEsqc7FlrC zaN2;2)wp{4@5NtIpyhkbhOF6zCYKypSHP#ob(hOGB4w;}xOKiMKC*niTozaCIRBvantH2E+BMTqnI;$Ag^#ihS@R z*5DFt)E=(bqNdH~iBBehEny4UOsiYSQ?E(T_jh`fw_j)8{7D)UE-2jkDxD@nfASkP zj8?K$5&JAZ%rv5$EsFT~yvzzC;XR2N)MsICoS^u-)h{=)6@>`wmW$hQ7xXCKb8@d- z+m2klMC;J^XjuY!L%y&cAJ11<%pbkWvZzKWl6FT8fGcNYQN>YvP2Z@I`uz^*Fs$*1 zZsTlL7B-Q?0jSVl@U5%l0FBch6`ieS(x->U{{}ekR7w{_n$w|5)kNopH za9=sfsc@>Tzv*r`E;NxpN&QU@nCq;Td=mKP_ViPIO&h)9zY&#Ut_Ahk)QlRPfX%Ul zsKqF9Q1eASj|LeYtt<T|hI^`o^3l@F$~lGto1%@h z(V%X_*kM2O?JU*~HphZZ$B6-69~@GgaZm}J)-U}3-VWsg-SgbiNr-K)p1KcJI0FA9 z14#5P;}?;T`i!%qs>4fAQ3J`HZw|;H�|~&nyn)lh1bVQ?QMn}mpUsFV4bsO2!r}B|&qGd~QZXsvLI~cDrVxug!%|G12rNv1!1Qkp6*OJqk zv}wRa7tON9x6p?eN*ZLlaMW!EFP4)ugP&Y zyWrPO6f9E)gB7A#hcpX=0;l%NoX~`X@xm7r2=?_@Y%q|!ED&_SrBXZ#VBTyC@ zPCIDt_A1YTSAM-)A+$cUbA#PBFQ1rrzy5{N42;y%yh_s!>%fimb1uu{dfV%2x0e^} z&6F!Z^J6Z#kfQzFq~1ee8l;Q}GQ7MWK+XB3G}?GV6X3cNf61Nj=a?^qs&?RGOz^6D zt>49cY(pm-y=D4+-xB5N6xMa!&_N zQj<_SwY(HOKb({=*+?CnL~U6-O}E#CrGNy25?P-^0n)!2`)KJ&^)_R&4~?= zn1LCi!s@Ta5px_-l#iYv75zyI6&jQnYtKT6Vu^;Dd*yekxed;=<>N9fV64(Ct&+>< zQD>ywm}bFW>!VL;?0g4!yc3`%U{D%clUSGWoK$#RsP{fZesF)de0vDNo& zcz37!=pUi^#;~gx*c~GXe_3t+y68h*5Y3A%NYFbos7{8 z5LM=!ZN3=Rp8>?6k!?TUn*S+zziSuK9g z78JD+>$}>StPwxAxa}+_S9a*y2`hIPrskWv1^oK$ou~lAUf;dytx})7Rn%#0{a6fB zu4IBj&)t06v(%t^a>}ym(AYRr(T??cCRxbT{;2cPwOjw&?*^Y#wM0ZZO-|dRQUXXp07wBJFbh4C7Qvxyz^a#*pl=(yH-w!)bw47n8e0k?HU=Z7*p_zu; zqf0`XTw>QBN{g*zuQ+OR7yWB;m9`zRShSP3y`E2pGg&PUV|zH^gV%fx6|>NK#Sd50 zpE%{zZkUn`3FR%Ozmj6ee041Q!_5`knW25j#K=kEC% zO=(WjkRLLBU-kCCbae3sM*GFJC^yOm+(~*=Ez@+e~8cAYY)UJx@6qwm|W)8cqmr?+nrMfL>`#lg5_D1UjV+ z@oT2|UtBZe7_3fO(R4n`=IN%)-@Ot3WWPyPSp*x|c>nLt53{DAv8`fCdDI*-oQ*4!=b z7ydp|YCj@0KE9f`a7@rj+h>0-7e@R+*Y(@e8dv5Z-%of1Ytqs`kqEMZrS#G5xt>$; zOk)mZzfCpW?KM}^>i#scC@lNUubNME(J-cCOw1ZPnq;IXk!UwZUHpqQCB;KZ0w5-) z_9WMaPzMes4(w%T@GC;SPV#{N+w=2;g-KXG!Ro9&me%o*S-EiXn@J@Oeseh z@*`?;LY?{JY-!G1OP6(Zx2`%}*=#t+cZ0vqkH{6kY{k^sTbp?u(e(bUfn2h2_NNLC zou6$Xe=Zr>$K|R{kdCBR9JbJ z*4Hcfc>y_E0~Z(i27bLSe;b_*#hv4##KM9b()|#VTqKY_;rUCl!;l}YH!;NT+E zY?=q9U%Nqkt$Iqalg87ip20#;=q!WGBJYUh@baQ01Y{&$X9Y$(H@SA+H~1nxb6vX+ zCH8=U-Rqp*v{z;JsYu6eCOq@8n zw?;?K4n&p6thn;rEgRX1{5?1J=Q)DShtH1?9 z`RG=8H@0}C#0judq*!ubThrQP)l-u5`^5^4ILWgbm5DrQ0k~>cR--U=FpD+RPi`EiFj!^63+~?)V63m z7kj0d2sVFk|HpmB5?2c|Y!d(li&>`f~Gp;&oV*WyZ2Dwj3*a(k<;mm z+y@u-z~l1rtj^9VYbhqJ>*7FHJ@MAoG7*<2(^;TqiVG0M`sMHN(^apvgEv-FCIwrP z;;607o|ht7&w8}^6o&a|*)PTy8a}6dd2WuwZ8)wpC^xh?EUMb-3)o7HroMV*O1W-m zG&U0zd;7Zh2tGy;dgp2ls-)KMWzdnbN+B!%)MLz=ZaGKgIBM)hu3+J5_~A7^-Yg%g zv|=O4RNMjl;v(5gy+T(lf)utxQQG}#G;!+dwy@FI1d z*=$YSE$~1;)R3U3r_x#LlvXd+g#Z>_CA1v)C}GLhDfHYO?E$L`He!4rQmG}G=FgX5 zi3mdm~96vS`yE-f!j58c_4~@_)LaY;A02R7tHbSI9A0HO}gm2VT1K3D|zt6IvXcnNiGXL2dti^V2#v zpT9U+1GjM6K|H6ns;Y`pBW<~88wUzDDRQI`=CIVC7-D(w{!Xd{(YHEl`_Q8uTtMy<)(aFTVq$AIHn_Ip5#SPG3xZ2dsb)yC zL57-oA|%7kR%gKtT3KF%rVwX>e;&aKAtK5FX|9?FV%G8yOqi(^(1`S5*|W3UDm~Kw zG^U4Uz%FV3MGz_YQ$-pmSBcJKn%SieUJ1ieUtPh+T=V*8Kz20Ti9T z<#UriAIv`0=jY^vj}m?P)@V}xA3v`gdp`lsbhX?vW@3)U_pc9g#x6NS z*n(eokSWDNuwt!ym)#Yx$$JFhv(WGPdGTA4Lq>}i)U1+Hl7&jB4ghPy)4@4;KstYJD8GjhrXsYdZ~l(>UbdNMTD2YKQIt;m`3$pAfw@><=nyEp-%8= zr%E4U$GWWTWeHvb%@A?M6~Jw#;Qokj`SEco8gqDtjWUua0JRYy&)u=v3R$UtX;iYa z^;j6ch)RJGC`8fFlz659sjw2=j26_%a;8<|j86pk@Bqd96GBJa$yMmb4SOh*3v8yz z!kKhUSA4+e)3yJprHH4?(~`R)qQ}afn?|`N|AKzqdr#UF>kFGf%?oO~jm>WWV(ru= zRCsHF$+(RWp}&p}C0r}?`qHSpFoK=CMR{TYooIculsM+kk~Z^=MH-QSM?Ng7V^!5< zefM81*9J>^AEp8<6bQB03;R(wjzYfIA3+8u1^u=4teqFR<_!sPbGXrNZf-Y{SZCIG zWGQQ@@=)~u`X~~z4qL+^)i|REY!usuig983uBQmMP^%rS7G40` z8|nB*53)=6!`qlo9Xe5|x-bu@RZN|1o;W)n_w_W19*2ba z`oEvuJzeZ%6cjK$2b;Ea2>%0a5%%8?A&oK=tNc}%Kc$30i&Gjk=19Rr z6}9@^n0-XiF(ac$sTc^b(3&vtI+YjTS(i~#FO6+$+~UM4cD4GG_4!2p`OOeYQkjbA5H{s*`J*$=h0S}Sv@I?>R zi9=+rIJIC-q2`EViUFPMPLh@@Pl%I?Ysv4iGs7srWiDLFiC>iu)T^l!o8WiZ*ZqWc zb$j=`+MZtDlFrlCzP!4IqO%*Uog(|KOzjH@sr)z{Pp-hEi5@i#NSeHa4@AR(=9j5) zlqJlK!jiZNa|PJg*kl=lEx&?IQbe$l7@9I97&U%1uOE6fMb}3Yy5IfDx7=O7SspAh z{hi9kf&5jkwGgEylUmC}3s$p8hzfd4Y36()=b0kCE58{Us?Yny&&RH1ciU4b6^S_i_pF zNHeR=%$(HwI=c9D&vUh8_Ka$9iHbM)&}5{_x3$WG2z_W2G?ZBC7#~b;{R+724i~yc zaR3&GVVF`(A*G(htIIo=4VZwiKHbU|nEf-4NP0`q5K#5gld)IuF|g%N{~rl6xyd#&uxBQV z89BHM6K5tvpm-DacL;x{iowcwd)p>?rBP)(c2Rff)UN`48*33S)@9JRL*=QZ&o zLDaL(CLrH>?3`*oC~?MeQSh;SYIfP`>OT5uynwG%g$FNBdyf^@YG4x#P83!${t6k#wClR*3 z%w!o+1TAw%$%f3OXv0mMM5F*V1-h@F&r)WAdCYvEJVl733Hg5t{pxlv)}XQo9Nb=A z*=f7j?poyu0y3^l16sxc@qe&W?5(VD%jz~M!w_O#aPxLF`qLYWG6C$Ob@YX`5aGFK zlG~0dIm!J9WH+!WCMn=zlopLA&ST>QUG~Dr6VPP%%r%-6q&(H})Wc*@K$FtIz!JpF z$rkZ8QxlwakS=xja%GV=l3L|71A#Uim6+k&s8T_{bxz5OTq?}!#9F?=dji3*vKRn&u1g!x0 z!WBwSlWaqZXvHu4l4y`%`S4#Hb)fWQ7ze&Zd3+B)Bvj#>el^eK-(_tDgba-xy47T} zOAaxW$nVoxVGtElekuNv1XakP4yZ|C7QwZ{6#T&sK~cnGsojm3EpSOL4Q6$#IO^s` zui<5S|3=jWYS6juD^kv5nPu6WPP)9u+VV<^K1+%!I5U{!Lp|jnDhsxxs&nN@fk@8G ziY8^a>qXkS)@)=F$O2cp37?}5aE*#q{AvrzB1oZw&McTz{;-()_i>WY!=E~5t~4x{ zMD)s3l9c@0%oBXE?sKHwsO!2|@mPE#IarK7%CRAWo0JC6Sgm=S`U1u-T%Dcu`FuN{ z1@@9U4|sQfOCtQ7^m{_ZxS&#{0z6892SLiqt7ZIMTeF6{&G)&!Ma!G3tdt`*)L#O} z3SYu$kC9tVr2QmHI`KDwEOM-r@~F$pt1yT&qsV%UGYq{Tyu>+kIwRuQo_}vi+E1{F zi*%WNYH(a*)v^UOku1l^xb>pA7;wWJXiz2QHKfK)npL;B(zuGk1jjhvPRC8N;2-UO zo|QvFg+unr78OEmhbj#i{B+R&x!eEt`hg#C$pHf`z9>J0vT5lwmF$uYZt3wm(%XWCQ&M+ws`H>VSHobF{7S%Sw z>V9wE+db~t9Y1_!!yxgp{|?JgcW_RP8-g$T0~32wJ)?<;0ZYtW&zvP6jZ}&5ZlpzD z6H&48Tz6Hwv#=$O5i>GV5URDk&EUPFK=(cA$GEX%MTZ1B(6Qmzvf|`PN&q z+ME6#W~D;LLp2!R@{M>DXDMt^$l=-n*>+pOu)lN}eH4jwk&whj|C5fqL9es3z37Xw zV-88^o5vjRVt#+l?1B+SH zE6Y5owA7jC({0Y{Yo1f6B0ek<`PC1CeA$JSccxQM<^5;@41++3d9ycQbJVfU}!rjB)H%T zz*y;`Xrjl*8F#_UJ6PvEDM7&W0Ug3$@=~TT9l;-^J(_5@`#Om~x{*^9eoyvrLrsC` zNA21w36?;HXs^sVueVsR(Z%IhRg*c&Do5IkEqWy7s)Vcbw>Isl-mMCYqjDYR0Z#0~ zm6YaUEPEQ`hD|+~Uj#7noYsX{R5s>L z-|>Yr(sUQ#b#QzvnWxJd+qLIKj%Xs;8DZ|Iu{9lys{r#Kf>A<*GNrU$7HE9D8{4Db z;}(Fhlue1{k?Q#bQHpb!wgbkY8t6_$vAP|8GgSt`H%lJsobQ#jFcvvSJYTSpM)UO8 zglvB~cDS;HEeS9^jqMw#D6B1(gWuN?B8=olvfhWrW7zuNHi?PU7IBf-Q?03ra&Z9W zsjVl_1ByjIeul;gN_t;H!HasUM#wf&OA?^*uvO2oR1R=-SKaCivc#^iqDcNxr=j+x zZ?0sxgcRszWaV-_f6fbw!n%k;iS8avZ_2`HdNE12)d1Xvkx>w2!1r)30ItCK z;n|LPW5!&a@WMYwOhYx*y;K0hO5w7t2ly{xzoP(+zcZEJT~YQ7=DWSe^8JFX?~Z54 zyAM-pyuJ7CuO^_*ch)3n@BWl*$-ui$wSk3u%pt{5U-bK7u?L-GisbR@pf!><+G8~d zfUDiiRd|JQLo&t;t?DC>BoW_z#p>q6rbzHKH%tH8It?swVTd9YC=m4!sij|X#&O1q zk^6tr(!M-d8iA+=msb~_2P81cm##K=uA503m2uI~&>aE}lYU_nvSV|J`Fc)`hIegd z82PyTE6c#`1dK4VOC{p8>jHqlA6b=lRg8VlD|>K~mNKLr)vj%O)!A0Klwkw%40uyX*5wG?FDC;WfYTzQ{&dq$t6Q8M@JjCc4%X@qCB?%65 z%8HeA4V9dNAz~o-04H_u2=N;Fh2QA_@OXt6og^z89#H9~P$c|}x=P%*kFBDZ*_x!2 zgI^;uu9QQzzgx^zZ+a%KKLUb6X4TSeWg`Zb&M)k#4YF3rA!!MoGgNv@uUy;h{tZI$ z4%iU6zsn0bv!s&M;6OY~%Y+;8_Ia(}@h&S_siH-wRoBRx9{#%3-)w1xISpG>Kscg8 zE(-55=aZ4QlV{ioIivdyDU`hfcMjC`}>g(IJ?J&$z?)Y_@3K65KUeV%r(ij;z2?>>K4P0`< zFLeinIX1XVNcsV^57-NK%UP|5EeV+GGkU~;ziw+b7!Aw1E^lXU(AX8C=J<-H`Jb6d zH3os($4fDnVml<#|H(qX0JBBl3~DXt3=eb{ zlNHvkeQgA8u^H*C3V$gNf7zb4bMSV*o@x;8wgqy=jOY4aFvgrTaxuo#fEof>&V6Q$ z{FO~53{3rm0#7$KvNk&38+vZG2O^_a*X|F6zG-O+w=&tMGNybN z?)Ki$*>RhYwCB{xQCz0Q1{~joJ^LC97Zl9l8+sk=-G%rtJ=1kL?NpOWyiUXMKGh4K z2<=QrDsd*U8$lhSjQAEMvf*r+AUuV$;uAhB37ic{0tC}G#007W_5Zr2 zAv3?HzD*V;@45PnJ*@d!^eVvMCqSSpfQNq){P$;#j+PUa_V<9LoQBxifxW4D0v0r^ z#o*9%K5w*pIEdFIF}`fxLBsl<&9%2#nLrs{g3o|HT+F$=lA%KZrc6Rm$2{YiwZ;er z$pVUPCx7>2-e4BU!V1Ac9)g87*h!C;YyS&k@VEsYbvylE9Z{8AtvAkA#)557{V%9f z01~<$-KGoGiU(d&4uV_vv9jwkV`e~Uvo^fiBRwUyYz_m!CY}Z2c;fRf59ba&$H+1@ zQgjHlBh`wuug~X>8Nv!CC|!&!{X?Ft?|8&7nKxn^Sg81C*q7Q{ifR9q*_?2=)*R=@ zW^Uy76%y5XqNGMkTIkL7QfYZ;XD`XI_7=K`vLhH5YqE99q!9E{K%x+7L)P5%3Om2DJ06t47UznY!6zmvf=!&N)xy*gTXEb{Mb@XN zI%J_2uWIf;$^r2}VP?_PT(kfz_{NYK`QX!Fxt3n~+IS7lV^^rsVOERR>Vx)Q$)BZZ zQ)o;Q)-`-;y#+;WEUNhF_A+9Z>$H#wy`&KJ^y73?s~mCAdB6>7Pbhv8r^|?xy4L<2 zg$_TfWSP(}G^>VBU0=#=**r#Vnx>%6UFSV%eU98bVlA7=nx-+iD_U2Qe{wKT4&9O2 zS!iY|5nZhBEfDCaX!2`|3dj5A4Cq0%_aaj=h%*6|Cq?7Lrl~_Q6(J+mBy2AvIfN`x zWx%Xj7ela^!Vqcl0ZBeHbN|Lj*5E{hX`$bsW3+%POcZ8XW?E)MlD7y!^!;X~Ib=v- zf~U35V5QBL&nL-G)u>AZL{@9AB}2vMln;hm8hx;J?7VW8sab+>mx?rZwAk&tzt~xd ztLp2}g}&);!%P;G&DF8ALE-Qys_4aV<(wo4trzZjLFpy(mzhXftExF#-mvWJfTY^} zsc(PN&sn#A?y}hUzCif3>k>r!kyU{)E}gj$W^1@gPHxXod0kj@byOBA*6GOVW+@(6 ze8k?=7#btV?-rE0a14~#Q+jzJCjdgw7sTRWn(|rK*R^)=?rH4GDGfsHl8R~8g%F1D zM!(PcS5LA@dEI%g#!$*oArAGO4Q%Pm1!!RuL^KMmGu_7p`?!|H0~;6UUh!Xltzsh zZk8I(mcVE(A!8jLO1Ok;KcVeIoPz8|^mY;udtG@ppIs-nx z96|{4xYH*Ly;lql4s_;B=e{EQ;sH`vbhTncL0>vOvZphUtGAW0V$;ePu+Kd4k^yM7 zg$D2AfIJGM1&i`S~?AZeUMyxit? zK6wX|%$W%#nA+5fsw5z242?5~vO6##*CJNX-(%sqxcz|798reafsq(I^m0F)1fW22 zlhmDd>QK;`Y{!oKfq4#3ZP3ywD?%-ol{GDW(;0ZBWuQ7Q4?P>MnpY7CcxGS(kvG< z6vJ&~PB7{$Tp`TaiK}*bm6oRw!;xSSuw})Pq>nrlFL~= z)^4W~1(-P%h&~HMq;n44=TW3=sADFFUO3R{I;$hMI?ZopVvMfb=VYgo`7y@%e_%o{ zabY;v@Ft+jCa@s&wg>Y_FTo$-kU^~YKz#epA4!h>uMq<7ah~~>X|1&j62Y~@GTt## zti5vjaaxU+ntRpvlt|2I3=*t()6+CTF|5G5GW!eyk9Q&qu`2|IQl;*5G5wW*>3nVM<1{RE=i9h$0kUAh~7bF_qJl%b`X+OHn{ZZZ(^ zqHGMAOmn5-MaH4`JZhnd8DJrX*E`vyX&}a}khZ*1GfP2uI08i6V1{+TSG1xBhO~2Vn`mmZ#KbFe2mo##Z9`Tuf99tVlsNXok*?(wg^uACb?Gm{qBU^TSriLdsQEL6D&neyPgA zAR$;if%`qRe0D^%!>@W)5Tc4+|knHN{ zUh6RM$Q8@=ugY=4$155Fy2W9=&TlBl&!z@%eY=+!`j=Km_jpH;HqQTR)9-iD1rtGA z=$Af?kl|j@A>#Qka(yl|`K+uwMgajRr5pxe#O?EgGH(D88=@%dC``T6&0H{gjhgU#T#)sgEWM4H3HRx2Dq#HEReiBlYi l@1>+H6({{ac`bY*R4h1|pnkdsoDtQQ9b{~y3mXKw%i literal 0 HcmV?d00001 diff --git a/doc/doxygen/images/step-85-solution.png b/doc/doxygen/images/step-85-solution.png new file mode 100644 index 0000000000000000000000000000000000000000..3b8932a25b1119abfb8cd160902c8b0e6fdcfa29 GIT binary patch literal 35810 zcmdSARZv`A5H32n```q3hrylT1b2rZgS$h30TO}}2<}dR!7T(24uk8UNstg+0t7$v zpIhhcR^9t@s?K@XwP&XGT0QIQ?yq~T?j5J2rGke|g$)1z@YGZl^#A~*Pybz*$cSGE zHcRpV0D6F$qP#(n)p1X7K8uR~?+3}$oWO<(bDlUQC!`h%Hz|H@cYYf;{=&PU#-gW- zD{PT>F^*DpjktFE0%JBUt3jy)w73-Lgta_Dx$?Fc@Ku{nsA+g>-^@38+7|07Iu7ek zPAPvsc{4IHWtkgpt{nypYgZq_#>U3FgPD4Rr{l)H4#6<7v9WW}tFHx7(a_M`$HAFleO~RA+t|yFpFIOU-Rwm)Jymj^Uw2}eG z{tcm3(BYKTRow!ZEHUVKW=txp9HWuY&GAvSwFO~N&Un!w)9_metM97 z{7En~GZXQPq=?;d-NttvcP+91{(hxz>yKYDo$RSzz}eQT(F??rD`JPcWhU~wAD(Mqv?TGt@l`#^(k%JLx z&vz60g~&EyH1%1535WZ6o-!|wS9G^uM?2Gh%iMMbZg`JQOvD$R{LNh-teC=t4t^wf z=+JJp8%^FBO)=ib;HOP_xo>#6zYRaTWp%)m81zHP+Vq*^NeTa3c65EeL-+D8u;<7U z*1aWf{zmo+B-;{rW+l59y8W>A`1r1C^M*vbtW!D>fr;%tdpvfWH%)^slk8Q_471wh zd{(cYUQ97aLv+(mAE;FH{I`mnoE$Oow)eK@EbQsXy6?YPb%9GR!X3fKt4*TaKk?>E zemY7f?ceZoTJa6~9sH*(zTSVf6(UR#(X_Z>C$84(O?DqCaQMP65*oTMC;mgEA$ZD_ z_uE0~HU4J2nTbizT!Pi#9ogp{&!9o-<;kw?D zeXw>aNE2a6Vj+LO6n%$<-W-eHZuoLt{b#fdAxpIq1zLMs#UO3EfBhJ|z5jl_{uLy3 zBV9~CCw;qd6-D0tLrptm9$^lt%xdC;)gm1qjq`|1X5Q_C8anr>WcNE`NICZ`EP-o*w_*+yS?}SIThXAn96DIgGa~F{ZT{TPO@a-)3KwrObFNE zwB*^l=c}SV-`~^PVTbzO`3Wyy^Ev~3*WLaz+a|h)rUJ>p-t!N9-S>yJ1hi8B?q(vM zD-bd!t!|oec)tSrKO!t1@ys{aS@>;N#Od$w)5oisA`7coO2*rzvyINH-x1HhtuNEE z-XU}$>31+&CHwTJ=m}mFwve~ww-x-m&Ul1uzBlZl5aE`xzr|-d|IWY3%E}_Ivv26i zaTV>Ds??pQHe-MUB)v_SYs$5zPE4(;vNtiH#+EDIks@%1TO*e-GlavhI%R^V+Al(Pm!4ez&vF?c1f7ZsHA?mikp@{qYNll zg*#kkU&&4XVo`Ta!{b%LMmo<-Nz$NAvn_FB@YsYBHs_;KvM_O#A7git8oh@S=Dlv@ zLS+=8?6*Si@-PMvZt!N}a+-@$!Weir|GSmOWoKjMW@BkUHEEZ@3l-~bHh6L(17LEK zJUa&xK;vZrU;n8v*BTA-_-E~p<*s;H0~N0Nv8Hn{DAKCtR6F5vCE@bEo@vV@eH>I% z=fzfc$1H;S|%9oKdI>Okp{? zI^USn7hI5c_{l9&iyt85NaibSnAoV0~rlVKd^WNJ`woin0!hZyZwQ;kj{^V*KD7 z@d;nW*)u&^xdjC;Dt;$JnOL;+i%Y5lbIkyr9xlIM7(R!fZ@cv*yDHbi3FnWY7{ z^Lg-1_&IsR&nDR?{m(d(#G!;PLIz-rcJIXE-0$F7^ZmB2ily@U*^4>I z{f<7%aTNR<0CyctjeSE$lsRewl-2yNZX#@`8U^>JAbW*uo4i>X7gK;uxUH=SvpfLj|KZ3viNq3rH;;Os`nCf9ese_=g26bl@1SO{pNnh8 zcR53&!V86kY&;A)!}nSRzrQ!ld2=wQZpA5CSjlOW629DGeV!@nQFMigTJQ(EUt-Yz z_IhRWyyTYde%jwwxe2@sgOxjn5|b47VC%4EF;5SXV0@@g;3-++r`xR&X_%m8tTrqt zg>X`s!4>`8q;i9InF?PrWrG(%uce-GjRDJ)wM$GN`<|Cyyl(^R7tBy>(6iA@0Pl?e z78)+Zo#5Ki{yg+Ii_8e5)4k)_Ob^PL(!J z=>dkLHh6l)uCKqoB|mgSg2U;@suz&O`**%?;Vs5l(C$LfhhmT}V$Nd3oV}MFR_!6X zR?9YEBr*UpxQv~_XFJ#fVi;%9vP`l_UCS)9#_7yLm9a6v8RY;@FKXs|6&q}Z1R#Q> zf3a>L!01tt-Z7QE8R=-7#;;t?sRTelnzbqf*4O5I=HfmSP`pN)L~g7m6i8dlGX08$ z#k&(QQ)OoHRFz?lB>TtKn#*(&!*T#Lx#f$I8+tDCv@|mfU6TIm-y6O!JE5{<+5Ryc zfTP5D>6iNg&&LZEu8#yD6eyAy(zAr#2H{dx`k$2v7P&e(K; zs&CadQ?>a6N<;nMorLTQgNK+3A4WQ)e>u-^lamTnAe67r`fm3;>t&Z9v+vfDe{ziD zZ6?6goOtwRNf!R*{RV}$XtBI6Nbh3z7Ig))$xc(dQSUocuy6u1p&O{T$ny}BuC$kP zu!~AVDVr2`H(CE^D}jo$ll`SQzWYrDb=*^ zvk?IC1kG%3F_^U^Nm~WakKb*0Q0_99^9oK#X0# zo7{Z0KQs7?<&kJ4gcKESz(w6=1A1W1wtlUhJW+-nFhoW|hRzfZBuXu93E!oh3E#}i zb1Xy>ag28-xyV3o@+<&#*S-VZ^)wvi8@Kl9!ZSMFwo| zU%E~+Yc~q!_Qs%^3pz%3xfSVzN0b74CA22LT=F#74H6Mm$EsHEC;M8RJJMmoa;?y8 zG;1?cXrJeco(BqE0-l)G>-_k=zC?xHp+m%?;i$M<#+EV0pm^3^29p}*a79>09cT1j zRe9vHcl~F94p9}kJk4WnI>WI?MP#pc{KPhqU;J;xXCi)Yqk}c85%B0u{s&K6pSXot z?t)#TR=RDwQ3mB`KdigQiV4zlaL*jn_Uqc(tzM%jXR#>g+kASUbVTj$U zcLu4~c-pwi^LEc;ezPYbj zx+ZpRj{ed)r(|?S5GR{R*%!Zq>T#76UpdoCmu!tM+?&L8lK&w>6!70aPJL}R+Bp;Y z5Q@@LsfA>sU%6+eB6HW(2e%e$l2*2aVN4~P0)`gjNvD>GKUUE=sh_(Hr`A&eC8|BX z&gT}Qd>v7ZEoDXR)(e&3h-D+h4sMyWKC7gXDIY7IK=1T|$};N#IS7RV^Pjpqo;2-3 zZI*>N0Is>MKF4Eh5oe@tZnsv{N?WLeFOJH&2&X3(sQTN@lGsP9v4fleUem>%@9uJk zn!$f64%Aw@J(>R&hFs779Iw!1G?1ZfK=|t@nVz^oCv5bd35krOH6Z_LLpDU)W*LZ( z;qsr1?Y?KjXuFT%OOm<5O7>JjE$`y+&DW2faZ5dh#UWefV!7rYs_kmN7>+j!Hri`A zD;g9gzIAtKH~drg8eN+a$enWh&9o`Nh{%w&rPQv6A=aQ=wZwZ_d?sx2#m#Lyve;P% zVWdubtErO_F9=C0_7;Y=1@fcZCi+KE@T*?R3`!$Y=R$2``OLVm4^~)H;~;w^v5NZF z!3j;5PmX12lnyIB+~T93Hn)&#n9rKF`jK$Ta~%t=@pS<%fVuvzlGqO8>d4{024@H>BIGPq?dXfIY z@)XS=P0dE+D29UFj1o~KV)m!9i1X8t_I~J8JpgfK?5OZY4Ay=eP2q2`m8mt=S4`H= z0gd6xIen1w;QTR6Mce51`U7XJ3MiTNC$w0&QAxdm^{JDE8rMo_ zke~k#h*CU4K+l!t88&Qf7ZE5I%`C<4|GYdO4t&?u5(asCF1skQ!;P9d%S>;l3b>v2 ztSJ8W9vWq=8^x6%Pc8?b>)KjaD3R=z9h=4o9}5=Ql1y+z8Ws{O7b-m>$g&rT@^&LOJ7%hxQL1O`~e{ z$jz&n=d%c+>^5gaGNRV@&FgPV*5gWsok})kVr+2N38fz_>Rbr~Inms~1?yYv6@Pl; z#}!4Sq2PVsosf&+Y24nhN$($?Af&FgqV_I`6| zN5#~YKU$nY@YD)UX}_F!7h;RwZtdWZzKLE^8Up2fKxM=4i}K^B$AjW6F``fz=sMRu z59y!dsm+byPt6FjA12QW*8Iq^WV=UUREh^MDqWvQ7G8e~Lc9rY1RFr@QvErMGGB=vSvTw<+0pWGG}1Ir9myn5 zsg^Cp`ffR)Ycxd=)N~iNe=-%qqi1sh<-8+@$hmLflK2N*y4NsvQxTE&x#0n>rBbl@ zenH*M-MgQGf@uRoW~w!+Gt{=kGmLsnplc^0JwUT@aoo4ErW8&w*96S|xUJba_59KM z_Lic%F3*M;ca_5)cNs-tWD~P*Rz5^2v<#${FyR#VaTi_j)ey;Q zidk?p-9nbZ!0!9w*wjMa1t;C0&+y*kS58r{c^#;llO1D;A(Gy{b#NJBjKk*K*H}({ zC|qjpl^h4jA)^v@+uf^~cPHydE$f7>SO_(8{QVQy>%Et}86xz4`8zk-+FU&s;|^eA zW!K_0M0iGYN>XBetVJx4v^!UAnvOgCjJAX-AgcCKm*DGT(6ZvTmyhGG*$44qTXHZiL?{ts52(UZ|-TPn{aofJ4!i58a zc%%FC#J3Cg{4k(3Pm|e(mcfeWZIBnlvy*^4#|;ug!`;7}+#|f^fo69aOXHpLSQ))P z72{3KI=QO2MY2*SDJ{~ku#x;!sPGMcfrukb^W_0--xS_(5prNT`yaLNJ#KciM-^N* zaw3+v_@63!o^D;ANTq4TaJU$AJ(tr|?;17%;%+%6p`E_7=7G$d1loL2mGQRMVtFC^ z2Coz%MtH8ycI)N3B3^~7v34Bh%d4;}6X~(*@V!wC6r0vTbG6`z!U12$$kfF;iL(f; z$p&Ann-5Hf`~yO+8=akhuCi_2->-d9DPt%`(ppT}!JndrUDap#+Z8s`lW}xWT(Bh>0;mtG^rb-3$68*sOHDyW~V2=6LqX zdEB>josAR@kaC`#%dk3vwnSVaQ2r1;_)!hXgdk7mj?Ggl;Gr-lmJ+v)i7Z^KCo6B%--^Q&aAcv>K$|1c?tEtN?nsm}sVX2YR)Lg;yR{$UR zKtU5OhaUKSk1VoA^{JQrL=z?MbSxnM7=Nf6Ikr@lW`?QuY5C@r?44u;MfP>$7l|*F ze{GJY&HEy}Zl`z>KByM%7BeTicq^7Cf^5n4?7@s#tmWO7MxC$;+9)>$($vy?MhErJ zEPAm%zpp0}LpW~rD%nff87r2pK3OV;SC3(t+Pl5>;IJqxC2EOJkMB+?VbmiUD-_UU zE(onm(Vo|g(@`gjXUR*n{KNqUreAvgY1dj1iz7?s)Jw^B57g$hxcx(*aSx8BH7HmE zfew)7=D#=by*p?KfcY;5VbaIg@(UI7FMP@S*u89v|1SWW113s6QXYSte@Z{D*jQFH zMln5QhMVP7&Q5NJOxhl)Z1-yy1PIw_%@|^*lwq8)aZI=9^dDX>Ev>m!sb^ZHhC4~fvUDk0nS7nv82j=2zB8|y#D#BE5P9&ONmVIe&?NOy$v3Y+g z)Kz@2qYGOyj9jVsYR*ml)IYTTG@NW-lrkw!QhjaT1$lq*fQ|Jmw=H;vAzC%gLvIVyF6mR6>#qgJXx(C@}^@Zt+!O zF)$kEi{N9F*}+Y0I@T*rG(9$>2b7Uc+4Ehxh~Xomzd9U4jkFT4R{a+Wx;$lSGlurY z2`E}lVw{DR`Q4=bdRC#-=B>1%nS}RaS;VFCX&r1y zL*rFrCCSyh+Ct1>DY{Q@r|7^zd=vI_?4JLkS_MBA#+MX+=t$*~E1xW*#I%?0_OFTm z`mxG5XM(q2VfE?sue}IT1)~NSU2rG)@>7&`x!&eaOu)K~{725urgPubQ`nn1M6%l! z-903alp5|M%m)PCsUjFo$LY_|#pKqVYzSMjr>p3ppe3##m?bsWa|vY7*+-m+Q5pdma23;{jE#*C;*`kl z?CE6~*|IuYdfg=n`U-@c+6+x7WRyz=5{orSK#7tLDwH(<5dNfxN1l6u$;tX4R{8$=6w{kcrofGWf63Cnt^O*aN{?Za zEl1a*gnB$Lc5M|BB5)!WXfVl_v2!etUTGTepM=5$Q2Vb=Vpetbbv65UvJO>>Vl10r15fwQ~`)T=qIDr`qD#-PctSKjc4(SBXqR zO<$3UHC+1r?W?uTOaGy;SzwXr-ue6y4$BHPM6$BT)m`zVh##!NxKg&l?sgc z#u!eJ-YIBW0SKmupDf%h?izp|F3jI7VY7FtdIV-8 z5xBk)mcXTgkFc%&EnmKOQc<)t|6w>kKmH-b>RzR1c!#L5Zoxp&1`ll8KEM%nV}bDr z-LpVG$b! z#9uEc`9t)`#!bCAU7 zDEc+ZJd%m)Nc~hmaV?Q(9~#>J^0w(oyZcwGBFX&Hgt$E1Wa;mBe_~i zW!zHPI0WfE%~3SwUpjqoYv=yFFr%Q3)VfOjuho_jLY8T)Pg>wRM-p+F6@ep95T$iN zNzp2>*%`R9`4C8=ld!`aS7Q2Wm!3uS%P>*wh7Kw~2V%5iK<$>A>%IgRkK<)C_f-1w zc2q*hMWFwtO-#tu8((^UOKO%-Y1X2gv^jHH+_9N{yN9Bo3o#iw*{dC1oQ`yh$I=Mq2e=+arm!@@sL`$P4E}je=rwOYEylp&+4n?n6br_ zzwn{M2H%ma47TuX+NQeZOWK?Pdghn*Q#lO1LE%#4JCbT- z55>Ldfeot`I9+lJN>5HB^Jh&_FG1zZV}@upN!vmIcfbR>NhKMPNVi~&_VU^cTtrlhr)swe=JqUc6o;v6|A2Qw5vU>9 z43nJNM^zO6s$-8=zX}FJ5_mQkDuiMtClnyF?t=x8!SnY@bGVoMi-jV)kY%6Uadl8~qfS3cfrp$=7gYa8@b@{+_d2*V|_Z@n%oa zpa1r#Ij!{daXXkdUeHeyM#8Z2wgU;F#pUdq6e;!k0r^n9ll4H7b9;NWMY^OC#&Hi% zTR0`NeiW4byNS2~9EXW88#spNO;4i2WLv6L)w6MeBIdJla=YpOz%(7f9Yss;{uhg~ zwLX(v``zi8EINQ5y{$Z&z$W_uJv#f#f{S_^SP$S3N$i<(B6RX?vvz_D%y?yi#c7>8 zan_vIqx(&Wf*~GO>h!(zu}nQdV#$Ghut0s1e~C3EzPR*oL9?6K;f{AiE8VjFodhG{ zwMz8{duHfqbLaOb!4r#1jYX$LY&Kg)qY1G%pQ{p^*Jt8ug4Ld3do7z86+Z4*64Hp| z&*t0VWeN{p;U|APtKctDXk?-8(UusW7el)e72cfWbq7~iQrS9>S}rWWa^)^v;SQ_^qa{HDd4 zD4wuyho|SBw7y-SG)Q3jbl7VCy8d)xzVF>$i>pz7vA&KM79GW_)rj@Tr;Y*KYUM6n zj94t@(*jClRDO%SfW9mnMotlK;RTPJuVR=mjYQR)o*PkCn7m?+EQL}q0q;cA!LPlq zV?bT$IF%h>zg~pUdQBPhmt^$xU9Ob!+npZKbhJ7o+)*pwMpo{*(6!)^8dA{N>a!Pa z#E7CDi9^viToG*G0cOle%eyuUbW(BSh$5d!gMmON6-$M|+(ZJ|(vnBL3K*x-L8DA` zzC=-{9f;^o5w>&4*LBStx5JF&`j(8Mardw_7i^yU(Ix>?M9%@BA_B$dnTll@+8IMT zlcvpELXCbU27Jr-Jsqvw>JSf{hX*WQ#!l#U87I7&b4p=k-|PLDu%?|ruJ>D!{WZ_Hx%QcDuCIzaOG z*lHY-ypcX5z&p$Ls6ZlCl50~E4&+853VGkQ>J0BOf_P&G)E_q8g)uves6Sp=EdMZI zp_}5JcwQ@i7sL5^mNC)N*jbrt)sl1Fpau@TMOR?Gj6P~+%-IL8u^ET*C-(_PzR?G6 zh930yU=T=M3Xn4%JaZ%te(N(cGkff}e!06GCs9)z`bJiJr$r6LvRNKEv(uETf&;R% zKQoh=5lLz%ArW4=b>+;B`tkWF)9B~I-4=Z^(FF9De_YBG-14k^Mun<44;NR#!7MF4 zxwr&nNr1_*elKZSclw@rf#Y>6-J69ca^mscWYcb9+s0900=q;` zh@wg-UVd`3mvy+Ds-QFLdET>q?G@@K0h-OGbDMhZ_c2%nTMeo-l6bK zlZUct$bcvk?D$~EN8mR-W&*Q!f^aO1{5Ru;3%Ii9>l#X|{LMu5f-_Iujl)%_s0BI1 zKLQBxzhCR~JS*wKx6hD`|0VYJ0`&@(_7xojBPjr|v)Y*;0gm8W)V-R!0g-NBrJ0Ns zw($PK6^UE2jMyTYy>Zj(3WGd>xo+Yq^FmT!KthjvvMI1mASscLAXC{uovMPvTpgcM zjU~hB`ybZPZXqWpYd5g$+J3y=K8VYa4hhh#i;1ssDc`Lr-mB}0sMxBI9$(qCFIZXp z8JCoBQx+HFBmQ{WdvvN0D1G03MdW9OD9VO~-A!8CdWsp=!UYYQqbYD(+nB)Y!w3BEP_2WE~BWp2f+okFJWN zD?svT9ALNxY}K?>KV{(etX=3u^JcL14D>J_433y5D6u=N+^RY>D>m^K-C-^4m41B6}fv>ADDBz55%SB0pQr1vv8a`?VYQU+<>}D85v=7VAOXz0efZ z(HRNfeA@jjB->FEe4{%l`i2)+^n@#fDy_(Cb&OH7>t#&Lc-!k&@{{4VSKt$#>16`4 zoXH{HZ%Lllkxrk5QY8Z-Xi))>l6_`geK#LnmfdS~1#eI-{a%LMsZ^(b7F>o6ws)6D zW;N}bo@%5t81KvFL)J3?r!~vy;j&`WYcLR{QA-;EW7(w4ebY}b0AZvd5ad^9^jUWS zcAPz|hvss-0+DGGP9}tw#Tb1M17iB21=y}GoB7#UZ>$Ii*|NzV7-Rizjr7>zWug5x8GBGtt!}x3YY|;Bexh${H0%OBqAT*%PkGyJCLA&f60!d zV)w!H{aL3+j}CUq=If*3TEH}9{xkU`?}$?Nlg(LZ*-I4Lr|k~6H#Nhdwtei~a?ouN zAW?W!Hy!7uE7Zy9*V9y1znitCh;ruD8S~s;;Mq1seiKe45mI)}9~#LfcE|DH;_`km zzw_A=g^!44Rw>kic!UdScHJJsdqP(p-7t#T^#E*Y#B6c=@vxm&w*P=p+GM1aS~!Ma zNgk7fMo<>F->TL+{Y9B_4^j~SM9+3rY!bnsA}=+24V(ACr_@0>e)ZT%G3MA^ToXg= z&oT9g4QF@QH0O-Px5!AMvh>fF1YKnHf_=3sbd0hZNMqM1Azgy?4v8oh1PZ*MNO=dX zPN{3KFF(b;&Ijo3mY3KvlyNw~ge-ckVEQs%pt;xg2@ipC4FIaYK;HYe4W45vSIdr^ z9T`$KfQ;4L-562+Eb(slyNTBb8oclH?P;j{YhJJ-zfAcV3brCPCJ60*=MTtyFgo_i z;S*bWM}j5#wi(A<1M`BPZ~MGcKAD_8Pic|{_s-%8!Q@?zQ@OueOh``I<@?L}sg4#H zH~z>N@o(B-txkni+!jF@n+P%6&;|-@9YkY#0#_}hl3nCuW4bZ`pH|{OkVjF${}W*O zu7K@!*W1%Oczay{J$>~O0Bp-Y4l?~RB-gbm*cXI-`3WhzMtcYbhg|CkO)zM<%H7jc z!0A1^cZtAZN~|PlGdZ-qNJ`Itc8j7k71e0qznWoOPYc-K3~z-!h6(_`s=qgK(B(*a zuq`9+Jm)%16>(TrzLf}LOmNw3pEz&I*W*Te4qt2t<=ONiUCudWG`abt$*V@0Gixxt`pZvDWhzpl{OgK;!m@xKRt0&XGj82^@ zNtu0|Z&@i_eMzy1xMc*sODNt8)xFyt9k>LivN6U6@3D-_-mnq_Ft|UyQK3^ZD9$!_ zSm$kRN<$p1u3CM+6$BEQ@lb#>)xzDW;cCe7p_(p>u(i?8#Y~#(GA$5}~7=d@mU< zIP%mwjg2DZw-E=Q-@HjE2zX;lJ@0drRJkgLt2&&1Q8I5n3V?AB3;vKxrF>fq(?Q|L z)tJdyVHlGG2jn=m%8Et{WpGz&C0ie_8!5VBXaHI2(mO+jAP!$bBUa@XpYt+0=U=t1 z*cnBtE4ex~$5dUnMi_HYtwCabqiCo%bx%Nfw8ckp_%-!bkIq8&xvH*lHp zcT$kuah+p2I6F)W(^kzR4QhN$q3JpLA^lsHmBRFh3*Z{I*?iP)9(aVj}|$2(Ar_~7deUmxy>6)JY6e~u{L6RbGV zG&r0iqA(^A-@%8-nT7mQAgP-Cfx*ywIJb!aR}W*$<4;4TyQDa*=S?cXq=?rkC04$y z|A-l%V42^y7&7!F2~*JpOw@g23@HuO>O^Rgv12k?zEcQ+B)h)0P5|gF1w9((R$lhL zZDT4gqa*6cz4ed3WLwRW%<9k^oyX+#-n8Nw^9wjsX1>jg#2go(Y#*sc=a*dYma|DH zTpG{`+BHinE;o83Xvaq0Q7YJ1UKYO^x(>JhB#S1ZhLQkbocqN*(GE97m%(ICm77zK z*d6^5R$3<S1Sx8e>S*c^a0zaPX>TH@e=RU=!MaVtNtBRP2%SA46Q|@H2qlf*R6P zKE!UvINA&@o8J02$r-MekQS3}M+j)&NpbVrrVCO3#PwI?=^EGRR6p446i=d#^^&TB zB>cOprYflrRn>QH(by^7>e#CB3+|s^E!oC+8W4>@7;+t=03yMo>R8S@?zRwjewvaY zj`1H0*j_Nkmqg>Ry}X9UrL2JOgdc_dCt~i67BZ~=jfU{V?_{$G4g&~>)qrdxMV}9| zmz1iZ@g=)IrL-0P&abN5$OV%4V+@|!d{GmiPssF@ z7f5HldQ;jKk`$K>%sUm~eB>U+rV>up5^o@vW$2J5W1|r=koyU0HJ>9UK&URdHAF372l=Qyar0dauRP~3- zciUaB2AJe7--zFt+sibIW!UIAf|r*FO_j?oV8>7c5SY6GY8AdGYp^s>tTFK2>`lk< zR+sC{7938A4wXj*Tz|tgrROKn-I^98a6wKPwe2HBVeXO#1cD`4r3mC`SD0^5^^XTTk{U})w2!n!OnCc*ZjKYOxq&01XDt` z%1$f)YL8y!Ul#?KQ}&@5jY0@I2X&>5;lUb`ENTwAFI{T0k%pAqaZers4XM!)h$EpGG5>ytW zvXNfxo$!pIB7fVOJ9^USA1CN=h@fzWnnDbDeClhB@+2fQJ64^c{U+$DXs?vQNJpHI zL5J+Y4F<%7=5n{He9q=&VKWR0|zqAErt0=!-{LGiu;4f`URWHdmdR6Ikkm*t8~tHW^6eIR744Q!>@Vw z^wrQ~)nhGb%?jYAMMmGdZ(uacL;$aA(R#j>LHa72x0;JQUALRxLeF3@+@Uf0rmCUc zNDve|VA;uSYZxg}cZ9xq87r?>!OjKA1{9ddU>%&U0H%ky(O!-S5M~=cpaAN?1*-*& zEmWKu)uZe?y-sF!xU8KcF>l%kV1=L~zd0wqvY8y^BoBzQxyylt$NG(){ zG$6w~jUI}4JjJt$S%mCOc`ZG1N`Wg z`dBxuG0OPCCoB()7S#iF8kQbQeo`Aulh&=7`yFY7pyhEQD4%m0ePLfai}d?QRErfeT&rKIYC~5d z%VQPQRR}Aaxd#O5YHf8xR4bNqa?H)lz;A|X`+^5+!%rLu4JrUom!99zS$y2wt!A1K z`(7Qaw=CW?!&tjXcqyjDLi79NzJi}jwH#J`&M&9;?P8@l;hAc?P6N}^uam8L1>$#> zo>gP~?tFf*uPe@P*+&qMit`4`v){g+r}AN^bV;5BnQM-96D5t>>2WS;==O4XY1Y=s z$)Tr8bC7U!oU$NOBrDTjk{G@hT5>ShBg!nnzBy_UN&i4Dd4-Kjl`_hy-u%!%gy2=) zS7Jfi2$M&RQkQr|@|3tT)pfJk9O_B4;&RQZF$JnwqK;u2DaVk@1$|wO6oDi;G$twE zRS~*%vS*kh?T%Y$I8xIgv6Hto2~d$VFpxx=y~4d}s&Si(Z9S<pw7sV58^O`44#KlEdGJzD8`tls%AQg(S%<{SoeHwjW14% z`e0bXC3BmimaphtUYZSG6}6paR+}=@rf0=tb#-i$1X1CRDHysq&#mR=IJ+K+-c1*f zIF69V$&M|D*xF&_^h~e*?0y=KaYgv<{Ktb=8RG>VUPx=c7x5~()t_FcPk^4Xo$_6- zXstQ2GbFz-)-m>7l?nTLmwvc(w_xjbNM7oFQU&n0+EqWI?3wFeo4IDColHp}wfJo} zPAtPoai{}7p76d`CX1k{{Cn~8biI8y^@DB7Ru8nRh)gz9xqgUGB)oKvFf(` zHgEYoe`b_}e*P^6&Pn$_bZt{I5)E~d0rbw5bGQ;~=`elT#g3925oJ2R z8%qXC*c~`aIfnRWX%FEb9|oUTxLU19sxYB%hp5w97TDtGT9a+sKYy=*m9Ns`>@D1d zgrkB2-HK_-+kj}iM#iU^qO!w4EL$*rVl~yloW!d4zui(4dMXnzrtIxa<=b{!IuVGT zZp}5-mx*Pb8VhWiI*OC&yPMf!CtDFYKtPk$zF{v#s0teo*pon(?Vt%#x+lMZyp-R7VEbq#ufwK@Xamp+IVr4;EJh z0T7-GWU;YTb_54_d&fo2IZ+I^+1O{0LNXQ{>+e-e<9?xzEHFWkYL@bBoH7AI-6kY7 z?ZDN36x>fX!-JUu4XbdK=T6E1K6norDv$){G5(h-gkq1HPr?5|fqZ6S9t$AU^p`sG zaUlOu%&fvbOOG=zZ>0C~%L?_GRpJNDtw}|2>-zwO<(l(d%disaxf9B&>iFhuk0lnf zft`fCRH;wPuQSH?@rgu&wK3pJhcF&KK2SfdEmptd3hPlp&K=*&!Yy*+J2V@IZ_X>v z<~KAts7>XG3651Q@vd!)Xm!w3v}Oh#bbuzSS*ncotnmNp1<28+Z5P6%N1HxoO4BfK zdTZ-TI`7MxF?Db?=|ArpKyv1wTjk*iD2BV4M9h&QOH%mo;EW}ZqkMiE$Y0;nC^t|c zRiJsyCqWSVa}Y0`mHk<}WH@o&F}To`9xXa!zzKH-nRDM*9jxI1F-Wh3lk5IS??`5@ z#{0-b^76|+txvAi4BpIzxwAh#|#v8y6micm+0}-D^@JrkVLdCS2Ac< z4t&Q-_J+KJH!mVrQlaL4CAL$(yZY=964zn!`V z&$TcB1gE1Vlv$8!3ST1$RU=j?R5dsE#N{2n`2^vRE3VH7P3I>jEqt;C+> zOoLW?VpHRRC9TXVLX3CAi$KrAFXN6DO`)98y;Ixasom7CI zDrJggz`}%8k1f+k_Fv#zGK{hx#K`>@A1eCp`$?OvYDIJ6;;?MciR6HzZ5Kh9a;D_7 ze+h`&9Km{!{^7R^FkGxg^V-pS5H9Go6%E&paBK(WCz#8jgTnN6O9XSYrMDwSfB3C! zrqdO2{gF=Yx1i!|?%U_2p>8%r$M0*Ol)<^|9kHyacm?cGoV@a@JN+qxe0u(@$SWnF z)+jyLopr>=+IrwYSjdmV6#O6@hv~6_;z&sov&%kx<{k((Hp<*!xa?TUF}Jd6Em!E7 z+v>SJwG2Ib9vH&C`98o=IX4-gh7ZtDt2e9LzyJPj1j!`YdOpT)-1?>Zzp<}I0$_pg-3 zT#BN!f?@E?mf5Mj12Ab+BEGDRVK3oL)g&EP(eV0_Wl+xl0KGs$zfI*$HUmMa5>hs+ zT!Y7SCZ!^oIkv#5tmLdf%PiTnhLu~_na^l5+lycTOi`S~Nqp7UZV13a02YAn-?IOv zT?b?qfTY$iRu#a~7#PfEfR*s0z%0!~HYCktHl+x`mQ=wMX8J4T7&C!fQhj0Z#2!J) zHW~nc;z+Sg2avP27?@yf^?HM5v%#RX+xD9Lg5#1Ie&wlP4cXl5MA zRiM?GTUlBfkA^C$1E@9Hx9+?9*1dN>dCN0v{n~TJ>+&oA*{w%$4N|w{p=Z{!N+wlP zM$j;$)t+S&>>y@~p#}XESEH28cy6t+#*g=KZoB#|aEOV0|BS=C*qkF2^kG8TQpx zoWYD%XVx_#IgzWHWpX7+jcv>Eah7#ad3dtL0&rNd+3Y11KCj=BWl6OKXg)SyxfFup z?sP;^*53SvSGCAL(Bsdjn=}6<-hsZ)!5qF>Mr5s zSAEN^C$)xVD3C1BfZo^OVm)Q_S8!yx# z#$}U4vd!Gk31H7RjAP|Rpa_5hY{6!w(yg<_QmWZajOPxlJgTtmWiI({U{nzkz%Z?T zU&j|uIM(dq;w~Y3AzLAXzsvtEHO=meaZi+7W&j#)kbjnk*}T=D?N?kl);Ww?o3B+k zWfJWTZnfuDmKWg=&*aE891WwW#u`zSBuVYIJ5S!S_nxObe(U~%;g$y0jk^xqyyu{- zM^UUIGqsyV4D8E+=FvwCHqk5)q1~QgAY~9Zhb<;eGV0wtt4+MAid-yh+mVZ1i^$$R z$@t}ls7hfZr%*0Ilmf+f6b;C7jhY?Zn9`K}$vpiGPGpRg3YFo)`HstivWXl#xu#wj zy)t<~T)G4eXrH-xoE6Chuy>ak#PDxp?htn>SEj6(Tc_f2S{(Q;()K##Z(OU7hC^pq znP3@boOOR4t(ehj&o3`8x-x5xR)2N%|9s{90a#m6$4@N3>jR&2o3(n2Ef@%31`}GX zSq6zcW0)m-iNI`n2>~Y3%z&}Bm_gD3PdwA^NWGM&X_;J3Y<_=Gsy_pext6WWpa=H* z{g`ZvFj?I$*LZv@$+Q#iq)`;5fcect)8^TgR*)@|kuO$F_~pT}9aP?1$-gT6W%4^T z{hU5ZybL)RvUvGL+T6X|C3yxHK*XA@wdFE^^X6VI25-M*S=Tnv8CF6oW-wE$J>TuE za^|I4mStuaqBvoqUHcB-a@)O6dBPT*LwV~0>$1Q9+OdUYc`8RytfK5z2)R@g%Y0ZM zi9Ub`0xdG`Av6dCH(`Td%l;1nrIs%SZEKP%dyFG9${ZRrGPB9qN>P+p@4=L3oz1K$ zxYM~z`E0IfY0u6oyusrgEhT%TuNFe!N|a@wYmS5*fNK4ZkOo|;4HbYU={)p2H_(C5mfiuqXq_B329t zc-TJEZO|RCm5InS5@HiLSIn#l1Q0hfs4;9+&ykC=S}Xn_dnOgbwjT>Gz~h%3&3-zZ zbcs-~ZOLVeqDY#Gg7XC72jB|_+70%NFH)0!Vh7An{xwzRG0mI2l!|ErC!Zo$KBugd z^6|v^BefU3LUwlZElG-hOKlqn~fcP=+`_MlF!3gN|9P_KX?OV6kP&B7nhw4OIZv z@Tn#j2tqD0KSP-LXbw-7Uk)^YFIdJ%Lhd2w7ix(yj4}ym$b8L}B8ED{#*g2Q(5U&AJQ970{ z5U>}+oTHVTG)=QKM=06WK)uoK^;QkG03N4uU|ees`dzk$f1^2n^hEcCJ0Fj&r!6;d z^#A+ZzW})N&{5g+ok6LZ%o1)-vU(;J#oyF2E-@Q`(Y0z+#u5mp&8wi?TSouS3k-yuppK+{N?z~ScHe);G+a$T)nmJ|IF`~rTB6Wa( z>^^vdNM$kCueE>Wa#GTo<64aYUAWca1{bQ3Hua4eTBoM|S>3$*;Jrr{zkbbLY%y&u z3+sj(_kHo|?>yt*J>j~4-3vgesNS3%4px~>v z9iNp&`Bv>E16V8E?x})7l^TswKl9opr={Hm*a>PqUr!gnmnRgt6WY(qZrrZ3$lQuy z4analQd%&Y8F2)FVPfh4#7v9h`fxa4*SR1T8M`Zh23UKv0+`Wk&8@C1%GYa4r|WOJ zYyX|c{^CtPyM4Zmt)VUC62A4%|L(}~rR9|#GoJH|ClWu zv^I6o`D>Y1$dqcb9wM5pS*@*wkF2!JQ}DdSIIfE?<9m5&=JJp&!nmTGd1rR<@|bLT zZGR99I2UQsk}}3Ai%SKVg&&+nVC+uKRZFE;SbnLSk#q4F@NZYjD)nEw)F#>}<(o3C zzn|J0t~E)&ZgzI_g=QI;k(wGJHwp}-sg^vFK|3!Ga1N}Ym91d4dK<3#xte1u9&2XJ zPkQV*i_1MlIDTUFHE;Ul(c@d?rRc2!tn0tG?}i)qN>xybUi|zYQHr1v(QdJ&Hn`(do9FgML5936{{Pf2`)pS9Mv!lb$SF2sHFpNKAd2;oU$ z-7HG|co2mdL}sMwLh4gmRg58u3+%nxNf40XRu3`5hVSMhz?#CsOG~ z#a8NFI(zQ{IdAg8QQUslufjEchU#y^B^=x{^ocvZlS(%TAMAm8n7|*CdY6o%|H0rnph(3d1C}xzja{E#vMINXZ~zie-jqg$3Fe- z!M(%cVwvOq_WljK28T{wHV>lCiORu=6WmywliH%{jbDd!O)@GI3~P%M%AS`>)p5_` z9_JqSJWis@%3NWaD;7_Xa=1dj*wVbx$2E*E;eIMp1)e@DwRmOE{2V-UWkCqDrCR5m zvU)@|F2)V13lAT!%0frF0A|-=byF*li}OIQ94R}?2Or)VOQZoX#|6=v4NpG3YyX=sdU3Y74x}0( z03eZSfTEAsE-N%ZKq}K{ZqXp~Jdbm)?0IF+gR=Q`P1WS`I4^T4&m+$x&*R*vvP7ta zm)S@}nU3Rv{CJjFJYGq$iZ)LS>RR-|iIqwo50(}vnH;f{Js5Pu^p+|#tZ-dD4bX}c z;9M@jHY{jT5z->xaU=CbR3nt{>uE><cc#R4r2~u`-etv`TTv8Hiem>84j%7dxmbj@Og69}wov^t%|!%d09Qk#+M><@e5%Dijf%rXo+<&A zU?Pfu7jt4EJedhX+0_m7x9Be-Ky!f9w`?SkXG>smqJLw%T+}sk^U_p3QQ7_U3E#*f z@mD11STuKbTjOv3{D`RGL&Jw(eO6CHeTq5G!|R@YXx;XRW-x;^!+~|fr@pg$@IWk{ zf-tWD0FL8s+BUFu^V406=X5NX1=V8RL|Tk2;u<#`Fsl*3L!7E;Y>8wB@G${hlqMHf z8v!R47oUU$=7-xO;DrFZjMnQ(f2V?CN|AO`$Y!ZVlS);3YLSZ$eovXwzb<`6%FnGl zABEtEWBF(xknrz}69|}H0i^o7L&8Xp0DM{r z=n#-d)d^IY)zM$kbxy(s+Q#9k5CN$Gxu6;hKtP#$9`{u3kA|1y#E2pXF8B>zP%P5_ zGFuv{rptL%ngl!Y#A%RupeYubdhvOrPNuRSQ(Uj^P!IgN3naRpGZV17$+)>n4?V>Ri2o5>qXsd zxw@bI{!!r%isiBG!^! zYOv_O10!Gg*6(MvHedOUbBRa}zamw2h@ea@k8|QY7EhV+EF*EC5`$zsf|ERAfLw-o z#S!6jaAGpv1$ch555r3%Ebwa|q*!Gr+F`CXd?FiprJWMEOp$;yBBzQ27+yGsu|Aa- zX>dfZO@6{5KuWI3Bg%4!hw-HExju8ic=v~&Se0Y+A%@D)I0GDUjm&WgZunRlRqCTg zl}bk9JdKp(?*fQp@ri{cCl)q7x_9f7;z?Lp!#KxRzrWUD&Yd@1T*x&%v~K(AHBVuN z%Dw5Po2E)rllsv0U*GolfE#CtWE~U3#TjFaCDV<7c%GQ;*?Vy5!r47bj@|souEXPu zbD5VEYZzmkI0FcX#gZgu5fOU>UYE?3NM-a031-x9wNW~BB5_>2p~uLsRzxxqCLtWz zl0_Meg@0p2M!(sTOeqHl#;U^Zhk}ll*n&d`#ZV)5ZxsNUL9QY;5bAKKE5f`h(d)95;Q^?hS`Ia}7J$mo8|2`vwSu@bW z+C4aQ#Wi0W8!yG;t`mzhLjYp2q#h8+7#kli{qO&+DV4qR&+LA3$AQtYA{ZkX7-wJr zfibx8Bmiw;nFCD$XlqNpD2%{Fons(n#xRv$RAa`(=mf^JfdGFOAPdxo5I@;5H*VSp zmD~szB#wBTST3#zeo-L3TnKQSaCSpJIKaS|BP=z>hzL{530Xt>s;Qwp@fnxP4BHRnMx(k}d~c%mKvfyk&9L~QdD!?~7* z8?QJw+nfgae#w+0W(Vc=fn~+jHQ2aRrpB< zswyon0@G12u>_+wdJ-6a9FX4*{k8{n`QgPkD-gYi7UHRe!_hr@N3`&1ES@x}{1euT zxTugK3js5qVG5eA#+~8j9Z(mm=a%7|3Gn$Hg`{&a5?3Vd4+n9G$zX($6*`>QBI}Pj z9H~Pq6U>iZH8mt9t?vH)z5}DB(l{r8hzN;vgIE(obBJi|=D~Ga_5f(jH{AW1e`u&n zt=qVx;y%^3W*V?ouHLqC^|t1wx=UVn#zZkts~yJ`ks9J8Dm)QDeDbf?P&j*W-Nt(B%>;gi)0E_+#~41rQ~{ZW9mbfz6ijl z@A<>Np)u(LO2f^3vm(Yv_6u}<|EHV8mezd3hu?eF{i_}wJ~)OMq#5NB?j9Wa@b%wl zYO4G2`!2ly!6)|YJLGvqFvb`dSS*5>&6BO$Nof0Z|OH9seBhSEy0F%nD zOVDm3b7G7+B#sY82)J+}i(Vn}F2waBAjngKeBGwLoIgYsJe&ND8GLJ=m`DWl;7p2pi8AWhAB*MU+k&cD=wtO9cR5HF`Uh94TyS{VD%z4IeMgi-Z8~@+_p~L@j=ex&?{M$cV z#TesWQ3tSMiIfOrNoFy~RwMh}*w}b!RsR!n=j6K<&tBen+^RKCj*gWWLoA+DbGfTEtZsYs@mfQW*ZWGI9oAw9vFDRbahpkGg3tQ1!nBWWcl z12AC_kq93KIi)`u=f;R=5>x>Sa9lqxK`o7`%P;A@;Pm+;W2HOpUA2FBf^+h`qOi)< zMDN5oh^SPQ-G#AuN*Gyt_K&{eC3D5$Q_1*Dq7`Nmu=>|MJ}`LT6CZu!F|)FFeEAnU zcMmhhxHlnw;Nz)!N4|=fF$Tt1p{4%X_ndS5oc!9&J5-ZkW8|>b1)E@e;te=YgqOY&JteES?}$DH;66rV3>cdl?9DTM3L}h2tZ>6dV@&ns8Mddec66DfCIwP^^;*y&)gKg=j1yq~7PDj`v5_%{l#xYX zB-o-+1%q59KXRhxhQzy-VfDY{&m;feyql-~Rqz7M|FaO4U5GW(csZy8i29 z3TzNDbs5W7g^-z2Y?lKS2mr(| zI#NfhL=^#x=C)k+`tIXrH*I=s-Dd7b&}~u$8i`#A>Z(k_ny<+d928I``v$BwRXq&$N}jx$pDe` z85u(=)7Y_~?d`8$(OziS_{i@6`o-GyTX!*qKBn4YGMy1QRT;V{0O~ams%6FFGJs`R zd~%on!U(-ZZO||8lcrXX7~&-wT7H493&3?^0gxHFOo}SgYJSGiq2ZA~-oN3BcfY#1sqTUPCmwkCappKo zyoSLrQHZFJYxvS9FDc|2bl$)8Kg}2;fN|mu1BU~ULm&naz(ClrZTBbd`t9X!I=N#( z+tT@GZP+^S-@jV7VcRYM+IkQbQM`PV3v6*f)HDa4u$4#x!f1zjO*U!IKWCX(^h7z1 z>*?%bKMx?tQwWoVtHdDT56Wz276Nn+f;cD(0YZIrSTLv98BLR;^ZUoElwp<~Scm_J zsFp57b4e2_RJ?yU@^_gVwH!AliJ-qS;lU6JHpC%4C@}X1FW5gZiaGzNo=ZvxaL!+Pi@^3p9wqE$k zr5y{~I~KNY*f#Lrzg)Lr+b$v`QyETR48)~>!W^a>i6oZ;K$43az8;eZBkf~pA14A{ z0y;V+0K*OvtxV$tfs4a^L9_D47gZKmEhwn;$y>z;WZ8Pl#ig zv=RUkrIAT1xdrh=!}#d%_{f1os#y+s06nyB*Td@vdOK$nvuYk#Z@%PovGTyfkFR>@ zapp+liinis=Ja+Rd;JH_Gie3je}1!tK;VpX&Kw5L;7C(TFe{mA6qU{Mh%vTd>)@J= z+dCJ|y7*OHor`95ENI`bZQw`0TEFGVA?OAynt`hL%wclXlbCQ52~l?ys$>pNzqcM$ z;3c<#@j#o`)RP}Y$*KtnLm{~}o#hLlP!XCDGML?B)5GB~Bn;mZ4yf0-kk#A>M0CX(#_13=8A6#>L4 z7SBwK9WL_;$4Q7SB!Xi4vZQq<*6H8=bFBr|3shhGw?Bvm3qa6{2<`3dH(h0TAzE=L82~W~pCP004PJoYVTPyVh^rWnydNBZI%b|B=<}cY*=({q-`V zW4el}9xJm+Zl|f>#T=K22=F8N2ESy+Ai}g7h5g6S z&LEa10bn&%0}vS}r^8sxG>$3~GsopzzVkh}+gufl2+&RMmjKc0GR6cbPWflWfq;1h zD0uZODxChJIT5Jh1SU}hz;P2|m;w-{l`y6VAP#cW7ESJ_~@`oFTwTD+MN%r z+tt@K3o}JE46L8s_uzTwtPty~9(sJ$>c<_I5hr30v}(;aeC^IlTk`?KO6Yq(Sx+3| zbvLPW^Z3|MC9^n@BSus+4F%$y04A}uXx8h_>FHR|wsd|gfYlpz{%Q3St2gWf08Wl@ z3>m=Xis6dgC~?561^jbkfN6pmF)W=#thi#jG+{T1m%^F1S8+-e10v0-{K;MI zJ&W>USA3vFBI+cOf5w_t;t=2bLH|qw*3pTQ{?+x) z4}9Z+ho4{$7;_w#F~-_jo9_MGTO(<;X4Ak&Z~v7Ts0DyI7#}?#J}`%gKzb_E0M#4l z$e=X5@&o|KfvD+@!{Ui_$AZ?gPg}5}%O6j?dc)3N|7qLCM+QY2$mEE{61wGTk4%HW z*u~6Uj^moJ2=(QK1{j2}5{o4xVNn^zf+$(_d+It7i5kY_ycClt5Xn-u;gpK5ja5)3 z-n7h|VCP^|0iIXZ>j4ij$Z4NAdW92Im#E%bE|*k>U@ozZ6M*jdEhlxgX-)y?-@NBH zE4Oc|Z~$STic11fxjdm!B~l7n36K@XM;tdHK5`Yx;+zOHF*2M;H;Z0E=^c9fb2p#6 zynBYc&VD8UYt`z+_FvghUoYZ&1P1n4pHQyMS{^l?I?#~ZB@no&x~6h*0zmZGg5b6=3L{yi6gjLRBwe+RY-EE`AtTK=gC^XZ zNhY!;OHW*lY%o)p$ph;_cMmFhrByB!RTdw_Dw^K>vc6f}3v)Vg|K`1aUbCxz^In3s zg@UL8Kt8I3VO1^`i3F$QWocRo`;L2M;S72nX{mj?|5x% zekQ!mekK6x9q<3w2OfHyxs1VSZ)^I*br-&<_k>B(pL^*y0Wg0!6&oKL(w+ALTd7Q= z_+Vm-dqp|&QvoBH$^g0?kps&+k2(9a1sw}oMVi$ccKm7elMij!;kXHCn=1-WmEA`v>|r?|~wk zjLmaza#V>kD3^+wR-|1~MJo^$$4A^)LJVurv?3YNVm`1q5rK1WdDrYO-ui0PzG^sz zCA4z&qbnbLlraX*+S(ew@#VL-6`~C*q5D@qMWTb+5j}|jL``LoA#+Nk8%h0uast#) zapLNw%$zYtzRUtDaIahTgX?xYxNe7b6_KXzq&6WBQiP>RUPsz|g_* zEsyW*-@IGEiy65n2ro$tqaV;S)4z5;>&T& zRzI2iUxjHZSS(dJ3a?0s$SnW>bCwMgaRG{n4~i^8j^PXnd#630KDMoX@d=Idk8NtI zPw4FZTlW3#fgM|(81mBxQ3U`5#vB!PibfSgBun~p%ya{xa}l|*wC9yU6|^F9U1aK# znN%#Dj3p9oD(R+@@k~1QEf{VW>%s|QhHa&HX*mj5 zsk$aVjl;^yH#DWQI^joiCgY3_t{r`;3R$fmQdd}}J54oiDkMn36^sL^0(O>$v|kbt z95>cfpO}Aa(27(P_o`qa?sjeFK?*=KMH00_jrvgVZDF4^L{c=A9d zs-VuZTq>HhB1tiv(sxYP;zC^}o=L@1NjH^@B@?ksI+jVtGpTqy9`T6Mj8Sa^>pdU( z#uHEN{rX*(wYB;!a}r&3^UojNu*+exiEC_4j1CE|Fb4Tr9R&pwTjf%bInskjq|(sd zgd>tvAs2%iOZXpwmZsF=dHLQYvnGd^=>KvIZZH;jVlD)fAjLREsn*U?keDrq`~ZY<8dva(~M5*ynkuuMW3dx`wZ6EUf#`ox0cnmXra=N}g^ zqX6vMJNneFkqwXR+xpbd@S!3ZFao4ME231KmqVz6Hn_MKfhzHz{-40TqDiY4_s%}^ z)VZhh)x4GTY*fR*>Re5+L>)B27(;VoYVo{6?~?XPcx`=R-^NGx?$~{3aR122 zSP1~DFl|B{k)PCeny!zyGv(4mr0wLmv5?&skx1PT_eg-{#N^Z{K*-36Q04x}90FJr zg||wE5%BIdR9KOXl`=|Qp1Of+nzeO>9OJM1BZ9+85JA#$OdIA0g6-5 zR4R+nsy?UWUk;Lc-oD^`{4(1ONvjjIK?s;`oq2S3;Oq%8aVSmdA?#7+5@!Hfa?>E@Qi}umEG&lwyYF zp@LVw5?%mC#!9>QjqDse^yv1Xor8x5_l=H>l|)N%LYQHG8%|{!gFB@}$NVHl({%dv zrtcMsE=A8PVLGZp7E|I>xmQLISdr-o!3i4(nRKkRB|R&jnK!4gJ)ddKrZcI)wFpS9 z-*#ZzQwN6+7iBp{;l{`Y-&j+~QE80R)2ISKfEeFEnphR60>DAp#^jZ6?dx7#fR$1W z1M5ij^4Hx1X+=q_XlzZ44>JbFSVEXrOtPjDTf(Ih6C&$mh{P6ra=Ed%j*aWRCcIeS ze=vrY=FH*~a!cnI=Cn1m*usJEI2-UOj~1HzHSnPJ+*82 z!0`B?k>a7z;sj4l6n*=`_8;T_)%rVi^rPL$bLck!wmG+g?{_NoO0A*`|oMA8R~4_VmG>dqxNMkB*F$0WdJZ4%dx`2&Gqr z-{*N{(`W)f2B&m#alIi_iN9i)zWy(s#3oLOK;or`&WhR+mYpld)97$)sYf*<`jU z*<7D!ZcL^U(Y<5)h9~w9jqlz+`oxYyJNJz4IWTs3)VK8lT;c@45kOQs0*^zhjq?5S z27t_>#0n8{?u9s|_i`D@q7bTx(7m{DWyW`3mehirAl4F`$6B9!W zES_v|9G_E&!WJ(VV`_SM^%ut9Fu5T7L6bN^c=-W%?ID8cCr*rw2J9R+9Jd0X5!PCN z-0a4d##Cc{vavqdRG&;%X~{wpMepE9@z6-|&`5FkaA|P=*ywn9e4;#2@`@#1EP16e zFO_-OBhTZWM`e$DoEU@SOaoW|0LCD?AYF$s#$AWGE{nNNJSLZ9+({(dRKm@qV!7sI zCKXF1+;l3IPP*x&4DKJ>H$1U#NRAYV9NA!o=S4$s z+!$jl5>jSHPT+F9yGbGzNH5=o6M-tRg}4XuT}UG2Thdp)wbur%CehUOUb7lcVj;{2++ejci{N3}yqcI(a9Udzkdg{Og3;y?fY`i==Q68J%W8>wq32%JD8yznX>>b@V zJU(==xc}h9_(Z^V3;hri1@HueEKt#?Y*X@_mmL4y`yMlMnF-7fqZW29EVJYGF zg*XMZiNFYcBrD&Xy6svUw3A(bIX804ntrnQL#a$glGtiIgZQyDY1M`LGi0{_&XpG zW;Q~{VWLTm7rpWQV}tw05QM8xGU3*x&`=j|tV^~urCYO^OgffIx~YVlh`Win6OXy^ zm=lY!nCl43$eo7U7%7jFw9z~&dtBIV#S)(=c@srGK2aVWFYgIRWH2##3QIQ;X!F7|9@wEyt0Z(%6c{ z)AEN}d3=CKq7^W>h?57BLpnJLCzMD!hL!LFk>C?)%L~G|Qze3!w!a84Tt#se8YbK$ z=%;goGCj)iqBmZ|?tP;W4v}xtO(w=zWrW_uC5#h;>VNcGJ1oB;vQ+c95AuPJ%IVx;si0DsN*m-VvxkNe31qh zQ&Ar6jQH>(1uqT}OXH+i$ECZcEiVN5#d&DbVV$1C^f4h@Fe>D!XY2g9n42hILdGCd4D$rU3F4_?{%FF#JKEIFZh z@!ZBuj}FOXMk)iiWTMI(KI-=yikm+vs)$&wIeGhaFSoRsre-*>9_-)#z{A^RdscN? zho}Gu4>2CEBTj&P-$N=HYi>+8HKv*xQ%#NO=Eii#qPCu;vjME#G|<0s*N$DoyZ0U* z+&2OZqbP%mh&<1Lf`~v>z{^pP5Vnn|r5|>=fW>Vz< zVp7duqt{$R&|DJYQvfTFz9Kgmi&Yz$6$x1}%}dB$^aq#w7q9*g@rwrt6$-a6ou3uv zmC!pbUUc)_{riU}%rgmQ0XNRZL={OAk@%pKnVeFi87QAk-*KI_t&Tu5yqECutN!D` z{_Wz4LB?deVT?7@CtLINt%Zi!?TsfKn{94NHVHyCq?#JjscG~{3GLo z*{0-0ubBJ8Uq7~gIQn@*A?-~%)#vc5jZ-5m8x<|f8^lW`3<59K&~bInU7jp9hs)v zGHX<`+8guPI)E)-Hm9$9&g@x@t@(zwLVaB(@I>n~R$HOItx(^+r2V`z7Ktr~MkaO) z9DHiW@NfU}#F|aJC84y<1wg)LB4lb0ED{M?bixh^4qyeErP#(Y>78n4j zBGl#s)mm>ncYg1Zd~-wcj(htjL7Q)981`4C!WNlfC;HnQw+gcipa`SbjfX8JAs!5k z1JY6sN_uBsYWkf{%j;-jFHixP25yY z*x)P3l&^pCbg^>#zdpS2(F0=LK_HE($tG-IJ`@0eCQBT}ETi!Qs4_@ZGQ#L6G#o>@ z5*$$#u)VZH_tAItEuEi5g*H5L;MRNk%`9es-!sNx`KKzt#03orda%2JDskLft~sq| zi5+is#G2t2)?ZgWOWJ#uw6EDT@Vy_e@8388e^>(pfQn%afP8(iX-1J6az!fUrGL_m zkUW`$AkW|Bc~t3vLS=&jU*A%{@x0@8OK<+lgZqca0eFVB5}Jey#x7IL+~WhJ>hB^C z?jj-tx45zJB*lUB>KK9eTht>yfoc3yDIE*4k+cHPF~4QWi7o55?iU{e=fH_Em-DhR zsT8Nw`RHSXP$hm>dX}_WT0JvaU_IxO*xESo-Jh&qvvB~TPCE?n=Z%wZN8~)RTTt^$ zz{{Y4K$!G*7}JK0z?saeO^_%@5F*PfUhL1Zlx<4hc=?J?-1UH9mY*)zx9IUw8q&~DSI*h$^3E7uN=jYq_o zL^(tdCShPPt|^d(b7pv|oCk7&Ibz;m`O#R2lkS3z7IETaX5renA`<(F2n031r=PWB ze)jLqIWeM4rPf?})1yO!d&l&xHIKzT@y=4oDWU+aoXp%+D^J-?Dz24B}MwEk-(el~Tv5?$-f~xq4PoeQtnw^>$c+AyL{> z3&C1AVn$JB42y`m7Unt^Mn|tqbuP?yE@)Y|WxvR!w@c&6;1rZUD7`;UU}}Id5rC`S z`65fJXDJJ;znXk({b=2qjRTOQ6QYK>2l?fX8YehR_VWhQ2l^1j^uUor833$UI$8Ah z@=H2D_Jx&V&VxW)h%r>5Ni|WG#F{udn5EOXZY(rhvr>{6!O;ytrw{<3;iAsuqO<2= zdb;@Rd22WC3tYZ9D3La)U~)MQzGdTTqENdz~F9wzyzsr|9-0AfkLr%vA1hcuFH5o^W^GUl<-D9d)uG)9vC<6EU#QdHBQ0*m~T#e{^nN}veSL} zfYl_jLH8GmGU@nPFP+!5sI`CNz{qGxwhIx>_3bQ3(?UF74c{LWi4hgEcVbDBJJq6e z=|zi=X$#fnH>Lmp3wcRIK~#s^E0s>g<{s1V+m$;)J&bJfhkXD5q(kq)KG0}!;+`kY zO}@vr%A#N6z=`zcV8q?07ZHK{C#65~t`&3VO!30qkTU7m+&PWEzJI&kLYx?LGz5KA zsQ|qATWCpt{u8e#Xy0%M%a>@o%^(PWhJf=Ra|#rPWcC1=b5r zg>2^Y|8(xMj&>q3$1XSFW%A2phiah_{324=;xLz!y4t2Xmw>K#OHZyTrMbyHpEDlk zp2s5*Gx-up9|kwHMV!Bwa5%Iu=7#V@M4rbrp*SaDw23!Gzy93$Q(b6Pz2W@%@|5)P z#YcJnIgx~y_+7c?@)ucJ9bH*qy%1H%X70HD>^p8aJKs{L?IR5@BJOd{392#sZM6!& z{Jmy^^G-j0suKw)*POcbs#9~#DZjl23CPQ2E|Y+25&)TtCeGxF!c(9SA|e8D>Rz1h zUYwumLaVB0NxpkYj&o&Rm4jzyz>A1_J6d}?kL)MGYD!sPy>QjrG3$#rpWEBn?h`@5 z3k0(r2YtY3ewiCN=ba0)&%iUvwWL0M^-FTisX&QUfMY_0_z?i;Bk%0_`1|@=8WUy- z(Tm5q$D{Bh6}U4*uRm|WRF^n~E_w9=pS;E)#BJ9%ypT;_``(kMD$MGLWr6jARv}mS znHyjAnH$f}XEUmK5OI)ft47P5gkMoM&iSX$o$BDx%C)54e`!y^nusbtv~b}wPJ=$3 z^P87<<~tX*T=ABUYNnR~!hk{-8)A3&qTI7FuY`J*7J8QEiO7>4VCW2pIOSSWU;5Np z)){;*WPw!!>g}BMrB9ycGfOlTWyMt5Um}gU&N&Z7x?kMeKGk8Xs(W$1Yf+AK9loOJ zhiqNaeBk~&(^#WdgE&ragU4pXdI77 z(ek?YonmSA9LWN!HdM&fefd-8_jb&ZZNjvnxZgz=%`ch7BhM52PFePBV=EtcTW_v8 zt?=igL@c4$mcIm5NuSPz*{((Y<1+8Rv@_e3@_PoWns5O;5d};w(HqdcB;T{-*_u~E z%Q_0nItu22WQPOiAnNOA?K{HW46EtO0;^_J$kl!JrgJ`f(>aA)hAT6xvVB(Lml+p7 z^VG>kw@*j;mh@+?e`&5M7D*6A7|-KD%y4jyNdWZva|4gd|s)mr6aA^v$?Yn&j)v zrIJgRF}4`Rny>4QBAannln@bCV#3p*SNQYfMXJHSx@BJZJC!h+&b{Ks|xL|6dodm z&Yu&HNQMWz_U(s0yO{8VsMk zfsq1c^ewzr&iPP#+hAXhqJe2q#vbDZK4(#mT(13L^?urUMe55eNMH>ANmuUUdrZQZ zTWScnxFv zFbiJ*zWW2}$>24n#{IwxEh63!z!$h|sq5mDJE|J?pwyk0{oyfQdmiz_LO@Z=ea?a? z0sCDBF^{I1Hi$iBWj0Ea>n!V6oFZ?6e@Cy5RQq8G9+B2R^YxDh8KH`74rjv_`=5I~ zmdk0vFY`aT;WpnhKrt0xx?#E-rf&N6P~}@4+?8N^@GjaQ`Imaxyo7^qhr|t4 zW$ik|z>h*g7G*R7dKZl8BBQCfTDm&Sn~|L9>9TCB?<=*${Us4C^YbDEn1T4DWtE`v z1H{Al_#J1h(pbT~>1{;7!fRQ8yX}b*QtS=%xK!uYVM{iiY>3acF7X)sM1qOzl%2aA zuMPX?fYSZ9>bR8)eQ0J`PHo#IwcGFb4eGb|H8-BDU-wg+X7ylxLUy`(W-npwRWcVd z_U=G9t#^kq*PzaEl~h@gp%|+H6c<+9l1kc_pBCA-PF|I(B!#c1J|^mu& zy#ijy6izlGjtGZ#gp#1KZznE0M5 z$ZO<8g=4>*o{KyhlylYXtcD`TVZS`={Ls7If=aWe8$H-GnCTtKo6+M0e0p*5=n=W% z?ybCN*AyO84P4Svvu974j^%peAWj87rpGKZV(|kxS}*jPRM6Q#FN6_bzoQQzQlmv3 z3rGhMtDQnZ%3^!`oCueU9~+NL-#^c5`YcU_P1xx|ndwUM_T4!)KKi~G$QSZ2 z&Ac5!*mzens_WRCvVHvW?y!Erdl%JYTiv91Gnl@5d?In5Msd=~{-bq=( zWuR(6dbJU~L&=(kA?fvzNKd3SP+oSjnP>|-8w#8opW4bTmAxM1F-@$hQ*@hpuFpC( zv3zC{Iuncwd8-@MB;KRS@YaW&v^e%0+jj*5ptYpN!(QM9V*Hjgo1tw2aYNBlI=Q5CZJ3Zp42y!H8McdO<-Y;m# zJ&TZqXj4U@h+u?+A4F~ypFAmbJ{l5dq*%B70&p7f{CzNS^7fWuBSlt=YL;6e!1#8s z+F5vgF-0U*iU7ly@RjmX>snNtz-vm&UpByN;oh^+J~p11aA7GOJ=j9&{yH2hrc^H} zsKEC5>yvf?*}9Ar$GC-XhQ_u}Pm@hn{gN-GHK4(@T11BiE!NUKt?oQs8|^DBIdABh z=n6WLNni&#wCJxa_abYf)ak_t@Q0C*+-uvMHIgnvkXzA`SJe(gS_&L$VR!MCSXs_@ zFuBSdKYMHGmCk{#S#}d|o=H}|@!gLP4`48 Solving the fourth-order biharmonic equation using a lifting operator approach. * * + * + * step-85 + * Solving the Poisson equation using the cut finite element method. + * + * * * * @@ -1004,7 +1009,8 @@ * * * step-60, - * step-70 + * step-70, + * step-85 * * * diff --git a/examples/step-85/CMakeLists.txt b/examples/step-85/CMakeLists.txt new file mode 100644 index 0000000000..dbc0f7dfb4 --- /dev/null +++ b/examples/step-85/CMakeLists.txt @@ -0,0 +1,39 @@ +## +# CMake script for the step-85 tutorial program: +## + +# Set the name of the project and target: +SET(TARGET "step-85") + +# Declare all source files the target consists of. Here, this is only +# the one step-X.cc file, but as you expand your project you may wish +# to add other source files as well. If your project becomes much larger, +# you may want to either replace the following statement by something like +# FILE(GLOB_RECURSE TARGET_SRC "source/*.cc") +# FILE(GLOB_RECURSE TARGET_INC "include/*.h") +# SET(TARGET_SRC ${TARGET_SRC} ${TARGET_INC}) +# or switch altogether to the large project CMakeLists.txt file discussed +# in the "CMake in user projects" page accessible from the "User info" +# page of the documentation. +SET(TARGET_SRC + ${TARGET}.cc + ) + +# Usually, you will not need to modify anything beyond this point... + +CMAKE_MINIMUM_REQUIRED(VERSION 3.1.0) + +FIND_PACKAGE(deal.II 10.0.0 + HINTS ${deal.II_DIR} ${DEAL_II_DIR} ../ ../../ $ENV{DEAL_II_DIR} + ) +IF(NOT ${deal.II_FOUND}) + MESSAGE(FATAL_ERROR "\n" + "*** Could not locate a (sufficiently recent) version of deal.II. ***\n\n" + "You may want to either pass a flag -DDEAL_II_DIR=/path/to/deal.II to cmake\n" + "or set an environment variable \"DEAL_II_DIR\" that contains this path." + ) +ENDIF() + +DEAL_II_INITIALIZE_CACHED_VARIABLES() +PROJECT(${TARGET}) +DEAL_II_INVOKE_AUTOPILOT() diff --git a/examples/step-85/doc/builds-on b/examples/step-85/doc/builds-on new file mode 100644 index 0000000000..a22675975f --- /dev/null +++ b/examples/step-85/doc/builds-on @@ -0,0 +1 @@ +step-12 step-46 diff --git a/examples/step-85/doc/intro.dox b/examples/step-85/doc/intro.dox new file mode 100644 index 0000000000..dbd5feca4c --- /dev/null +++ b/examples/step-85/doc/intro.dox @@ -0,0 +1,251 @@ + +This program was contributed by Simon Sticko. + +The material is based upon work partially supported by +eSSENCE of e-Science and the Swedish Research Council +under grants 2014-6088 (Kreiss) and 2017-05038 (Massing). + + + +

Introduction

+ +

The Cut Finite Element Method

+ +In this example, we show how to use the cut finite element method (CutFEM) in deal.II. +For illustration, we want to solve the simplest possible problem, +so we again consider Poisson's equation: +@f{align*} + -\Delta u &= f \qquad && \text{in }\, \Omega, + \\ + u &= u_D \qquad && \text{on }\, \Gamma = \partial \Omega, +@f} +where we choose $f(x) = 4$ and $u_D(x) = 1$. +CutFEM is an immersed method. +In this context, +"immersed" means that the mesh is unfitted to the geometry of the domain, $\Omega$. +Instead, $\Omega$ floats freely on top of a uniform background mesh, $\mathcal{T}^h$. +@image html step-85-background-mesh.svg +Since we no longer use the mesh to describe the geometry of the domain, +we need some other way to represent it. +This can be done in several ways but here we assume that $\Omega$ is described by a level set function, +$\psi : \mathbb{R}^{\text{dim}} \to \mathbb{R}$ such that +@f{align*} + \Omega &= \{x \in \mathbb{R}^{\text{dim}} : \psi(x) < 0 \}, \\ + \Gamma &= \{x \in \mathbb{R}^{\text{dim}} : \psi(x) = 0 \}. +@f} +For simplicity, we choose $\Omega$ to be a unit disk, so that +@f{equation*} + \psi(x) = \| x \| - 1. +@f} +As can be seen from the figure below, +the level set function is negative for points in $\Omega$, +zero on the boundary, and positive everywhere else. +@image html step-85-level-set.png +To solve this problem, +we want to distribute degrees of freedom over the smallest submesh, $\mathcal{T}_\Omega^h$, +that completely covers the domain: +@f{equation*} + \mathcal{T}_\Omega^h = \{ T \in \mathcal{T}^{h} : T \cap \Omega \neq \emptyset \}. +@f} +This is usually referred to as the "active mesh". +@image html step-85-active-mesh.svg +The finite element space where we want to find our numerical solution, $u_h$, is now +@f{equation*} + V_\Omega^h = \{ v \in C(\mathcal{N}_\Omega^h) : v \in Q_p(T), \, T \in \mathcal{T}_\Omega^h \}, +@f} +where +@f{equation*} + \mathcal{N}_\Omega^h = \bigcup_{T \in \mathcal{T}_\Omega^h} \overline{T}, +@f} +and $\overline{T}$ denotes the closure of $T$. +The set $\mathcal{N}_\Omega^h$ is sometimes referred to as the "fictitious domain". +Since $\Omega \subset \mathcal{N}_\Omega^h$, +we see that the numerical solution is defined over a slightly larger region than the analytical solution. + +In this type of immersed finite element method, +the standard way to apply boundary conditions is using Nitsche's method. +Multiplying the PDE with a test function, $v_h \in V_\Omega^h$, +and integrating by parts over $\Omega$, as usual, gives us +@f{equation*} + (\nabla u_h, \nabla v_h)_\Omega - (\partial_n u_h, v_h)_\Gamma = (f,v)_\Omega. +@f} +Let $\gamma_D > 0$ be a scalar penalty parameter and let $h$ be some measure of the local cell size. +We now note that the following terms are consistent with the Dirichlet boundary condition: +@f{align*} + -(u_h, \partial_n v_h)_\Gamma &= -(u_D, \partial_n v_h)_\Gamma, \\ + \left (\frac{\gamma_D}{h} u_h, v_h \right )_\Gamma &= \left (\frac{\gamma_D}{h}u_D, v_h \right )_\Gamma. +@f} +Thus, we can add these to the weak formulation to enforce the boundary condition. +This leads to the following weak formulation: +Find $u_h \in V_\Omega^h$ such that +@f{equation*} + a_h(u_h, v_h) = L_h(v_h), \quad \forall v_h \in V_\Omega^h, +@f} +where +@f{align*} + a_h(u_h, v_h) &= (\nabla u_h, \nabla v_h)_\Omega + - (\partial_n u_h, v_h)_\Gamma + - (u_h, \partial_n v_h)_\Gamma + + \left (\frac{\gamma_D}{h} u_h, v_h \right )_\Gamma, + \\ + L_h(v_h) &= (f,v)_\Omega + + \left (u_D, \frac{\gamma_D}{h} v_h -\partial_n v_h \right )_\Gamma. +@f} +In this formulation, there is one big difference, +compared to a standard boundary-fitted finite element method. +On each cell, +we need to integrate over the part of the domain and the part of the boundary that falls within the cell. +Thus, on each cell intersected by $\Gamma$, +we need special quadrature rules that only integrate over these parts of the cell, +that is, over $T \cap \Omega$ and $T \cap \Gamma$. +@image html immersed_quadratures.svg +Since $\Omega \cap T$ is the part of the cell that lies inside the domain, +we shall refer to the following regions +@f{align*} + \{x \in T : \psi(x) < 0 \}, \\ + \{x \in T : \psi(x) > 0 \}, \\ + \{x \in T : \psi(x) = 0 \}, +@f} +as the "inside", "outside" and the "surface region" of the cell $T$. + +The above finite element method that uses the bilinear form $a_h(\cdot, \cdot)$ +is sometimes referred to as the "naive weak formulation" +because it suffers from the so-called "small cut problem". +Depending on how $\Omega$ is located relative to $\mathcal{T}_h$, +a cut between a cell, $T \in \mathcal{T}_h$, and $\Omega$ can become arbitrarily small: +$|\Omega \cap T | \rightarrow 0$. +For Neumann boundary conditions, +the consequence is that the stiffness matrix can become arbitrarily ill-conditioned +as the cut-size approaches zero. +For a Dirichlet condition, the situation is even worse. +For any finite choice of Nitsche constant, $\gamma_D$, +the bilinear form $a_h(\cdot,\cdot)$ loses coercivity as the size of a cell cut approaches zero. +This makes the above weak formulation essentially useless +because as we refine we typically can not control how the cells intersect $\Gamma$. +One way to avoid this problem is to add a so-called ghost penalty term, $g_h$, +to the weak formulation (see e.g. @cite burman_hansbo_2012 and @cite cutfem_2015). +This leads to the stabilized cut finite element method, +which reads: Find $u_h \in V_\Omega^h$ such that +@f{equation*} + A_h(u_h, v_h) = L_h(v_h), \quad \forall v_h \in V_\Omega^h, +@f} +where +@f{equation*} + A_h(u_h,v_h) = a_h(u_h,v_h) + g_h(u_h, v_h). +@f} +The point of this ghost penalty is that it makes the numerical method essentially independent +of how $\Omega$ relates to the background mesh. +In particular, $A_h$ can be shown to be continuous and coercive, +with constants that do not depend on how $\Omega$ intersects $\mathcal{T}^h$. +To define the ghost penalty, let $\mathcal{T}_\Gamma^h$ be the set of intersected cells: +@f{equation*} + \mathcal{T}_{\Gamma}^h = \{ T \in \mathcal{T}_{\Omega}^{h} : T \cap \Gamma \neq \emptyset \}, +@f} +and let $\mathcal{F}_h$ denote the interior faces of the intersected cells in the active mesh: +@f{equation*} + \mathcal{F}_h = \{ F = \overline{T}_+ \cap \overline{T}_- : \, + T_+ \in \mathcal{T}_{\Gamma}^h, \, + T_- \in \mathcal{T}_{\Omega}^h + \}. +@f} +@image html step-85-ghost-faces.svg +The ghost penalty acts on these faces and reads +@f{equation*} + g_h(u_h,v_h) = \gamma_A \sum_{F \in \mathcal{F}_h} g_F(u_h, v_h), +@f} +where $g_F$ is the face-wise ghost penalty: +@f{equation*} + g_F(u_h, v_h) = \gamma_A \sum_{k=0}^p \left(\frac{h_F^{2k-1}}{k!^2}[\partial_n^k u_h], [\partial_n^k v_h] \right)_F. +@f} +Here, $\gamma_A$ is a penalty parameter and $h_F$ is some measure of the face size. +We see that $g_F$ penalizes the jumps in the face-normal derivatives, $\partial_n^k$, +over $F = \overline{T}_+ \cap \overline{T}_-$. +Since we include all normal derivatives up to the polynomial degree, +we weakly force the piecewise polynomial to behave as a single polynomial over $\overline{T}_+ \cup \overline{T}_-$. +Hand-wavingly speaking, +this is the reason why we obtain a cut-independent method when we enforce $g_F(u_h, v_h) = 0$ over the faces in $\mathcal{F}_h$. +Here, we shall use a continuous space of $Q_1$-elements, +so the ghost penalty is reduced to +@f{equation*} + g_h(u_h,v_h) = \gamma_A \sum_{F \in \mathcal{F}_h} (h_F [\partial_n u_h], [\partial_n v_h])_F. +@f} + +

Discrete Level Set Function

+A typical use case of a level set method is a problem where the domain is advected in a velocity field, +such that the domain deforms with time. +For such a problem, +one would typically solve for an approximation of the level set function, $\psi_h \in V^h$, +in a separate finite element space over the whole background mesh: +@f{equation*} + V^h = \{ v \in C(\mathcal{N}^h) : v \in Q_p(T), \, T \in \mathcal{T}^h \}, +@f} +where $\mathcal{N}^h = \bigcup_{T \in \mathcal{T}^h} \overline{T}$. +Even if we solve a much simpler problem with a stationary domain in this tutorial, +we shall, just to illustrate, still use a discrete level set function for the Poisson problem. +Technically, +this is a so-called "variational crime" because we are actually not using the bilinear form $a_h$ but instead +@f{equation*} + a_h^\star(u_h, v_h) = (\nabla u_h, \nabla v_h)_{\Omega_h} + - (\partial_n u_h, v_h)_{\Gamma_h} + \ldots +@f} +This is an approximation of $a_h$ since we integrate over the approximations of the geometry that we get via the discrete level set function: +@f{align*} + \Omega_h &= \{x \in \mathbb{R}^{\text{dim}} : \psi_h(x) < 0 \}, \\ + \Gamma_h &= \{x \in \mathbb{R}^{\text{dim}} : \psi_h(x) = 0 \}. +@f} +Using $\Omega_h$ instead of $\Omega$ in the method will give rise to a larger error in the numerical solution. +This is often referred to as the "geometrical error". +However, when the same element order, $p$, is used in $V^h$ and $V_\Omega^h$, +one can often show that the method gives the same order of convergence +as if the exact domain would have been used. +Furthermore, deal.II allows us to independently choose a more accurate geometry representation +with a higher-order level set function, compared to the function space for solving the Poisson equation. + +

The MeshClassifier Class

+Even if we have used $\mathcal{T}_\Omega^h$ to define the finite element space, +we will not create this submesh in practice. As in step-46, we shall +instead use the hp-framework. To create $V_\Omega^h$, we first add an FE_Q and an +FE_Nothing element to an hp::FECollection. We then iterate over each cell, +$T$, and depending on whether $T$ belongs to $\mathcal{T}_\Omega^h$ or not, +we set the active_fe_index to either 0 or 1. To do so, we need to +determine if a given cell is in $\mathcal{T}_\Omega^h$ or not. +For this purpose, we will use the class NonMatching::MeshClassifier. +The MeshClassifier takes the discrete level set function, +described as a (DoFHandler, Vector)-pair, as arguments to its constructor: +@code + MeshClassifier(const DoFHandler &level_set_dof_handler, + const VectorType & level_set); +@endcode +When we call the reclassify() function on an object of this class, +each active cell and face is associated with one of the values +{inside, outside, intersected} of the enum NonMatching::LocationToLevelSet. +Here, "inside" means that the level set function is negative over the whole cell +so that it lies completely inside the domain. +Analogously, "outside" means that $\psi$ is positive over the whole cell, +and "intersected" means that $\psi$ varies in sign over $T$ +so that the zero-contour of $\psi$ goes through $T$. + +| LocationToLevelSet | $\psi(x)$ for $x \in T$ | Relation to $\Omega$ | +|:------------------:|:------------------------:|:------------------------------:| +| inside | $\psi(x) < 0$ | $T \cap \Omega = T$ | +| outside | $0 < \psi(x)$ | $T \cap \Omega = \emptyset$ | +| intersected | $\psi(x)$ varies in sign | $T \cap \Gamma \neq \emptyset$ | + +Each active face is classified in the same way, according to how the sign of $\psi$ varies over the face. +MeshClassifier lets you query this information for a given cell/face via its MeshClassifier::location_to_level_set() methods: +@code + NonMatching::MeshClassifier mesh_classifier(dof_handler, level_set); + mesh_classifier.reclassify(); + + for (const auto &cell : triangulation.active_cell_iterators()) + { + NonMatching::LocationToLevelSet cell_location = + mesh_classifier.location_to_level_set(cell); + + for (const unsigned int f : cell->face_indices()) + { + NonMatching::LocationToLevelSet face_location = + mesh_classifier.location_to_level_set(cell, f); + } + } +@endcode diff --git a/examples/step-85/doc/kind b/examples/step-85/doc/kind new file mode 100644 index 0000000000..c1d9154931 --- /dev/null +++ b/examples/step-85/doc/kind @@ -0,0 +1 @@ +techniques diff --git a/examples/step-85/doc/results.dox b/examples/step-85/doc/results.dox new file mode 100644 index 0000000000..b8cfe758ff --- /dev/null +++ b/examples/step-85/doc/results.dox @@ -0,0 +1,20 @@ +

Results

+ +The numerical solution for one of the refinements is shown in the below figure. +The zero-contour of the level set function is shown as a white line. +On the intersected cells, +we see that the numerical solution has a value also outside $\overline{\Omega}$. +As mentioned earlier, this extension of the solution is artificial. + +The results of the convergence study is shown in the table below. +We see that the $L^2$ error decreases as we refine and that the estimated +order of convergence, EOC, is close to 2. + +@image html step-85-solution.png + +| Cycle | Mesh size | $L^2$-Error | EOC | +|:-----:|:---------:|:-----------:|:----:| +| 0 | 0.3025 | 8.0657e-02 | - | +| 1 | 0.1513 | 1.8711e-02 | 2.11 | +| 2 | 0.0756 | 4.1624e-03 | 2.17 | +| 3 | 0.0378 | 9.3979e-04 | 2.15 | diff --git a/examples/step-85/doc/tooltip b/examples/step-85/doc/tooltip new file mode 100644 index 0000000000..927b4dc20e --- /dev/null +++ b/examples/step-85/doc/tooltip @@ -0,0 +1 @@ +Solving Poisson's equation using the cut finite element method. diff --git a/examples/step-85/step-85.cc b/examples/step-85/step-85.cc new file mode 100644 index 0000000000..a3a8e8be15 --- /dev/null +++ b/examples/step-85/step-85.cc @@ -0,0 +1,714 @@ +/* --------------------------------------------------------------------- + * + * Copyright (C) 2021 - 2021 by the deal.II authors + * + * This file is part of the deal.II library. + * + * The deal.II library is free software; you can use it, redistribute + * it, and/or modify it under the terms of the GNU Lesser General + * Public License as published by the Free Software Foundation; either + * version 2.1 of the License, or (at your option) any later version. + * The full text of the license can be found in the file LICENSE.md at + * the top level directory of deal.II. + * + * --------------------------------------------------------------------- + */ + +// @sect3{Include files} + +// The first include files have all been treated in previous examples. + +#include + +#include +#include +#include +#include +#include + +#include + +#include +#include +#include +#include +#include + +#include +#include +#include + +#include +#include + +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include + +#include +#include + +// The first new header contains some common level set functions. +// For example, the spherical geometry that we use here. +#include + +// We also need 3 new headers from the NonMatching namespace. +#include +#include +#include + +// @sect3{The LaplaceSolver class Template} +// We then define the main class that solves the Laplace problem. + +namespace Step85 +{ + using namespace dealii; + + template + class LaplaceSolver + { + public: + LaplaceSolver(); + + void run(); + + private: + void make_grid(); + + void setup_discrete_level_set(); + + void distribute_dofs(); + + void initialize_matrices(); + + void assemble_system(); + + void solve(); + + void output_results() const; + + double compute_L2_error() const; + + bool face_has_ghost_penalty( + const typename Triangulation::active_cell_iterator &cell, + const unsigned int face_index) const; + + const unsigned int fe_degree; + + const Functions::ConstantFunction rhs_function; + const Functions::ConstantFunction boundary_condition; + + Triangulation triangulation; + + // We need two separate DoFHandlers. The first manages the DoFs for the + // discrete level set function that describes the geometry of the domain. + const FE_Q fe_level_set; + DoFHandler level_set_dof_handler; + Vector level_set; + + // The second DoFHandler manages the DoFs for the solution of the Poisson + // equation. + hp::FECollection fe_collection; + DoFHandler dof_handler; + Vector solution; + + NonMatching::MeshClassifier mesh_classifier; + + SparsityPattern sparsity_pattern; + SparseMatrix stiffness_matrix; + Vector rhs; + }; + + + + template + LaplaceSolver::LaplaceSolver() + : fe_degree(1) + , rhs_function(4.0) + , boundary_condition(1.0) + , fe_level_set(fe_degree) + , level_set_dof_handler(triangulation) + , dof_handler(triangulation) + , mesh_classifier(level_set_dof_handler, level_set) + {} + + + + // @sect3{Setting up the Background Mesh} + // We generate a background mesh with perfectly Cartesian cells. Our domain is + // a unit disc centered at the origin, so we need to make the background mesh + // a bit larger than $[-1, 1]^{\text{dim}}$ to completely cover $\Omega$. + template + void LaplaceSolver::make_grid() + { + std::cout << "Creating background mesh" << std::endl; + + GridGenerator::hyper_cube(triangulation, -1.21, 1.21); + triangulation.refine_global(2); + } + + + + // @sect3{Setting up the Discrete Level Set Function} + // The discrete level set function is defined on the whole background mesh. + // Thus, to set up the DoFHandler for the level set function, we distribute + // DoFs over all elements in $\mathcal{T}_h$. We then set up the discrete + // level set function by interpolating onto this finite element space. + template + void LaplaceSolver::setup_discrete_level_set() + { + std::cout << "Setting up discrete level set function" << std::endl; + + level_set_dof_handler.distribute_dofs(fe_level_set); + level_set.reinit(level_set_dof_handler.n_dofs()); + + const Functions::SignedDistance::Sphere signed_distance_sphere; + VectorTools::interpolate(level_set_dof_handler, + signed_distance_sphere, + level_set); + } + + + + // @sect3{Setting up the Finite Element Space} + // To set up the finite element space $V_\Omega^h$, we will use 2 different + // elements: FE_Q and FE_Nothing. For better readability we define an enum for + // the indices in the order we store them in the hp::FECollection. + enum ActiveFEIndex + { + lagrange = 0, + nothing = 1 + }; + + // We then use the MeshClassifier to check LocationToLevelSet for each cell in + // the mesh and tell the DoFHandler to use FE_Q on elements that are inside or + // intersected, and FE_Nothing on the elements that are outside. + template + void LaplaceSolver::distribute_dofs() + { + std::cout << "Distributing degrees of freedom" << std::endl; + + fe_collection.push_back(FE_Q(fe_degree)); + fe_collection.push_back(FE_Nothing()); + + for (const auto &cell : dof_handler.active_cell_iterators()) + { + const NonMatching::LocationToLevelSet cell_location = + mesh_classifier.location_to_level_set(cell); + + if (cell_location == NonMatching::LocationToLevelSet::outside) + cell->set_active_fe_index(ActiveFEIndex::nothing); + else + cell->set_active_fe_index(ActiveFEIndex::lagrange); + } + + dof_handler.distribute_dofs(fe_collection); + } + + + + // @sect3{Sparsity Pattern} + // The added ghost penalty results in a sparsity pattern similar to a DG + // method with a symmetric-interior-penalty term. Thus, we can use the + // make_flux_sparsity_pattern() function to create it. However, since the + // ghost-penalty terms only act on the faces in $\mathcal{F}_h$, we can pass + // in a lambda function that tells make_flux_sparsity_pattern() over which + // faces the flux-terms appear. This gives us a sparsity pattern with minimal + // number of entries. When passing a lambda function, + // make_flux_sparsity_pattern requires us to also pass cell and face coupling + // tables to it. If the problem was vector-valued, these tables would allow us + // to couple only some of the vector components. This is discussed in step-46. + template + void LaplaceSolver::initialize_matrices() + { + std::cout << "Initializing matrices" << std::endl; + + const auto face_has_flux_coupling = [&](const auto & cell, + const unsigned int face_index) { + return this->face_has_ghost_penalty(cell, face_index); + }; + + DynamicSparsityPattern dsp(dof_handler.n_dofs(), dof_handler.n_dofs()); + + const unsigned int n_components = fe_collection.n_components(); + Table<2, DoFTools::Coupling> cell_coupling(n_components, n_components); + Table<2, DoFTools::Coupling> face_coupling(n_components, n_components); + cell_coupling[0][0] = DoFTools::always; + face_coupling[0][0] = DoFTools::always; + + const AffineConstraints constraints; + const bool keep_constrained_dofs = true; + + DoFTools::make_flux_sparsity_pattern(dof_handler, + dsp, + constraints, + keep_constrained_dofs, + cell_coupling, + face_coupling, + numbers::invalid_subdomain_id, + face_has_flux_coupling); + sparsity_pattern.copy_from(dsp); + + stiffness_matrix.reinit(sparsity_pattern); + solution.reinit(dof_handler.n_dofs()); + rhs.reinit(dof_handler.n_dofs()); + } + + + + // The following function describes which faces are part of the set + // $\mathcal{F}_h$. That is, it returns true if the face of the incoming cell + // belongs to the set $\mathcal{F}_h$. + template + bool LaplaceSolver::face_has_ghost_penalty( + const typename Triangulation::active_cell_iterator &cell, + const unsigned int face_index) const + { + if (cell->at_boundary(face_index)) + return false; + + const NonMatching::LocationToLevelSet cell_location = + mesh_classifier.location_to_level_set(cell); + + const NonMatching::LocationToLevelSet neighbor_location = + mesh_classifier.location_to_level_set(cell->neighbor(face_index)); + + if (cell_location == NonMatching::LocationToLevelSet::intersected && + neighbor_location != NonMatching::LocationToLevelSet::outside) + return true; + + if (neighbor_location == NonMatching::LocationToLevelSet::intersected && + cell_location != NonMatching::LocationToLevelSet::outside) + return true; + + return false; + } + + + + // @sect3{Assembling the System} + template + void LaplaceSolver::assemble_system() + { + std::cout << "Assembling" << std::endl; + + const unsigned int n_dofs_per_cell = fe_collection[0].dofs_per_cell; + FullMatrix local_stiffness(n_dofs_per_cell, n_dofs_per_cell); + Vector local_rhs(n_dofs_per_cell); + std::vector local_dof_indices(n_dofs_per_cell); + + const double ghost_parameter = 0.5; + const double nitsche_parameter = 5 * (fe_degree + 1) * fe_degree; + + // Since the ghost penalty is similar to a DG flux term, the simplest way to + // assemble it is to use an FEInterfaceValues object. + const QGauss face_quadrature(fe_degree + 1); + FEInterfaceValues fe_interface_values(fe_collection[0], + face_quadrature, + update_gradients | + update_JxW_values | + update_normal_vectors); + + + // As we iterate over the cells in the mesh, we would in principle have to + // do the following on each cell, $T$, + // + // 1. Construct one quadrature rule to integrate over the intersection with + // the domain, $T \cap \Omega$, and one quadrature rule to integrate over + // the intersection with the boundary, $T \cap \Gamma$. + // 2. Create FEValues-like objects with the new quadratures. + // 3. Assemble the local matrix using the created FEValues-objects. + // + // To make the assembly easier, we use the class NonMatching::FEValues, + // which does the above steps 1 and 2 for us. The algorithm @cite saye_2015 + // that is used to generate the quadrature rules on the intersected cells + // uses a 1-dimensional quadrature rule as base. Thus, we pass a 1D + // Gauss--Legendre quadrature to the constructor of NonMatching::FEValues. + // On the non-intersected cells, a tensor product of this 1D-quadrature will + // be used. + // + // As stated in the introduction, each cell has 3 different regions: inside, + // surface, and outside, where the level set function in each region is + // negative, zero, and positive. We need an UpdateFlags variable for each + // such region. These are stored on an object of type + // NonMatching::RegionUpdateFlags, which we pass to NonMatching::FEValues. + const QGauss<1> quadrature_1D(fe_degree + 1); + + NonMatching::RegionUpdateFlags region_update_flags; + region_update_flags.inside = update_values | update_gradients | + update_JxW_values | update_quadrature_points; + region_update_flags.surface = update_values | update_gradients | + update_JxW_values | update_quadrature_points | + update_normal_vectors; + + NonMatching::FEValues non_matching_fe_values(fe_collection, + quadrature_1D, + region_update_flags, + mesh_classifier, + level_set_dof_handler, + level_set); + + // As we iterate over the cells, we don't need to do anything on the cells + // that have FENothing elements. To disregard them we use an iterator + // filter. + for (const auto &cell : + dof_handler.active_cell_iterators() | + IteratorFilters::ActiveFEIndexEqualTo(ActiveFEIndex::lagrange)) + { + local_stiffness = 0; + local_rhs = 0; + + const double cell_side_length = cell->minimum_vertex_distance(); + + // First, we call the reinit function of our NonMatching::FEValues + // object. In the background, NonMatching::FEValues uses the + // MeshClassifier passed to its constructor to check if the incoming + // cell is intersected. If that is the case, NonMatching::FEValues calls + // the NonMatching::QuadratureGenerator in the background to create the + // immersed quadrature rules. + non_matching_fe_values.reinit(cell); + + // After calling reinit, we can retrieve a dealii::FEValues object with + // quadrature points that corresponds to integrating over the inside + // region of the cell. This is the object we use to do the local + // assembly. This is similar to how hp::FEValues builds dealii::FEValues + // objects. However, one difference here is that the dealii::FEValues + // object is returned as an optional. This is a type that wraps an + // object that may or may not be present. This requires us to add an + // if-statement to check if the returned optional contains a value, + // before we use it. This might seem odd at first. Why does the function + // not just return a reference to a const FEValues? The reason is + // that in an immersed method, we have essentially no control of how the + // cuts occur. Even if the cell is formally intersected: $T \cap \Omega + // \neq \emptyset$, it might be that the cut is only of floating point + // size $|T \cap \Omega| \sim \epsilon$. When this is the case, we can + // not expect that the algorithm that generates the quadrature rule + // produces anything useful. It can happen that the algorithm produces 0 + // quadrature points. When this happens, the returned optional will not + // contain a value, even if the cell is formally intersected. + const std_cxx17::optional> &inside_fe_values = + non_matching_fe_values.get_inside_fe_values(); + + if (inside_fe_values) + for (const unsigned int q : + inside_fe_values->quadrature_point_indices()) + { + const Point &point = inside_fe_values->quadrature_point(q); + for (const unsigned int i : inside_fe_values->dof_indices()) + { + for (const unsigned int j : inside_fe_values->dof_indices()) + { + local_stiffness(i, j) += + inside_fe_values->shape_grad(i, q) * + inside_fe_values->shape_grad(j, q) * + inside_fe_values->JxW(q); + } + local_rhs(i) += rhs_function.value(point) * + inside_fe_values->shape_value(i, q) * + inside_fe_values->JxW(q); + } + } + + // In the same way, we can use NonMatching::FEValues to retrieve an + // FEFaceValues-like object to integrate over $T \cap \Gamma$. The only + // thing that is new here is the type of the object. The transformation + // from quadrature weights to JxW-values is different for surfaces, so + // we need a new class: NonMatching::FEImmersedSurfaceValues. In + // addition to the ordinary functions shape_value(..), shape_grad(..), + // etc., one can use its normal_vector(..)-function to get an outward + // normal to the immersed surface, $\Gamma$. In terms of the level set + // function, this normal reads + // @f{equation*} + // n = \frac{\nabla \psi}{\| \nabla \psi \|}. + // @f} + // An additional benefit of std::optional is that we do not need any + // other check for whether we are on intersected cells: In case we are + // on an inside cell, we get an empty object here. + const std_cxx17::optional> + &surface_fe_values = non_matching_fe_values.get_surface_fe_values(); + + if (surface_fe_values) + { + for (const unsigned int q : + surface_fe_values->quadrature_point_indices()) + { + const Point &point = + surface_fe_values->quadrature_point(q); + const Tensor<1, dim> &normal = + surface_fe_values->normal_vector(q); + for (const unsigned int i : surface_fe_values->dof_indices()) + { + for (const unsigned int j : + surface_fe_values->dof_indices()) + { + local_stiffness(i, j) += + (-normal * surface_fe_values->shape_grad(i, q) * + surface_fe_values->shape_value(j, q) + + -normal * surface_fe_values->shape_grad(j, q) * + surface_fe_values->shape_value(i, q) + + nitsche_parameter / cell_side_length * + surface_fe_values->shape_value(i, q) * + surface_fe_values->shape_value(j, q)) * + surface_fe_values->JxW(q); + } + local_rhs(i) += + boundary_condition.value(point) * + (nitsche_parameter / cell_side_length * + surface_fe_values->shape_value(i, q) - + normal * surface_fe_values->shape_grad(i, q)) * + surface_fe_values->JxW(q); + } + } + } + + cell->get_dof_indices(local_dof_indices); + + stiffness_matrix.add(local_dof_indices, local_stiffness); + rhs.add(local_dof_indices, local_rhs); + + // The assembly of the ghost penalty term is straight forward. As we + // iterate over the local faces, we first check if the current face + // belongs to the set $\mathcal{F}_h$. The actual assembly is simple + // using FEInterfaceValues. Assembling in this we will traverse each + // internal face in the mesh twice, so in order to get the penalty + // constant we expect, we multiply the penalty term with a factor 1/2. + for (unsigned int f : cell->face_indices()) + if (face_has_ghost_penalty(cell, f)) + { + const unsigned int invalid_subface = + numbers::invalid_unsigned_int; + + fe_interface_values.reinit(cell, + f, + invalid_subface, + cell->neighbor(f), + cell->neighbor_of_neighbor(f), + invalid_subface); + + const unsigned int n_interface_dofs = + fe_interface_values.n_current_interface_dofs(); + FullMatrix local_stabilization(n_interface_dofs, + n_interface_dofs); + for (unsigned int q = 0; + q < fe_interface_values.n_quadrature_points; + ++q) + { + const Tensor<1, dim> normal = fe_interface_values.normal(q); + for (unsigned int i = 0; i < n_interface_dofs; ++i) + for (unsigned int j = 0; j < n_interface_dofs; ++j) + { + local_stabilization(i, j) += + .5 * ghost_parameter * cell_side_length * normal * + fe_interface_values.jump_in_shape_gradients(i, q) * + normal * + fe_interface_values.jump_in_shape_gradients(j, q) * + fe_interface_values.JxW(q); + } + } + + const std::vector + local_interface_dof_indices = + fe_interface_values.get_interface_dof_indices(); + + stiffness_matrix.add(local_interface_dof_indices, + local_stabilization); + } + } + } + + + // @sect3{Solving the System} + template + void LaplaceSolver::solve() + { + std::cout << "Solving system" << std::endl; + + const unsigned int max_iterations = solution.size(); + SolverControl solver_control(max_iterations); + SolverCG<> solver(solver_control); + solver.solve(stiffness_matrix, solution, rhs, PreconditionIdentity()); + } + + + + // @sect3{Data Output} + // Since both DoFHandler instances use the same triangulation, we can add both + // the level set function and the solution to the same vtu-file. Further, we + // do not want to output the cells that have LocationToLevelSet value outside. + // To disregard them, we write a small lambda function and use the + // set_cell_selection function of the DataOut class. + template + void LaplaceSolver::output_results() const + { + std::cout << "Writing vtu file" << std::endl; + + DataOut data_out; + data_out.add_data_vector(dof_handler, solution, "solution"); + data_out.add_data_vector(level_set_dof_handler, level_set, "level_set"); + + data_out.set_cell_selection( + [this](const typename Triangulation::cell_iterator &cell) { + return cell->is_active() && + mesh_classifier.location_to_level_set(cell) != + NonMatching::LocationToLevelSet::outside; + }); + + data_out.build_patches(); + std::ofstream output("step-85.vtu"); + data_out.write_vtu(output); + } + + + + // @sect3{$L^2$-Error} + // To test that the implementation works as expected, we want to compute the + // error in the solution in the $L^2$-norm. The analytical solution to the + // Poisson problem stated in the introduction reads + // @f{align*} + // u(x) = 1 - \frac{2}{\text{dim}}(\| x \|^2 - 1) , \qquad x \in + // \overline{\Omega}. + // @f} + // We first create a function corresponding to the analytical solution: + template + class AnalyticalSolution : public Function + { + public: + double value(const Point & point, + const unsigned int component = 0) const override; + }; + + + + template + double AnalyticalSolution::value(const Point & point, + const unsigned int component) const + { + AssertIndexRange(component, this->n_components); + (void)component; + + return 1. - 2. / dim * (point.norm_square() - 1.); + } + + + + // Of course, the analytical solution, and thus also the error, is only + // defined in $\overline{\Omega}$. Thus, to compute the $L^2$-error we must + // proceed in the same way as when we assembled the linear system. We first + // create an NonMatching::FEValues object. + template + double LaplaceSolver::compute_L2_error() const + { + std::cout << "Computing L2 error" << std::endl; + + const QGauss<1> quadrature_1D(fe_degree + 1); + + NonMatching::RegionUpdateFlags region_update_flags; + region_update_flags.inside = + update_values | update_JxW_values | update_quadrature_points; + + NonMatching::FEValues non_matching_fe_values(fe_collection, + quadrature_1D, + region_update_flags, + mesh_classifier, + level_set_dof_handler, + level_set); + + // We then iterate iterate over the cells that have LocationToLevelSetValue + // value inside or intersected again. For each quadrature point, we compute + // the pointwise error and use this to compute the integral. + const AnalyticalSolution analytical_solution; + double error_L2_squared = 0; + + for (const auto &cell : + dof_handler.active_cell_iterators() | + IteratorFilters::ActiveFEIndexEqualTo(ActiveFEIndex::lagrange)) + { + non_matching_fe_values.reinit(cell); + + const std_cxx17::optional> &fe_values = + non_matching_fe_values.get_inside_fe_values(); + + if (fe_values) + { + std::vector solution_values(fe_values->n_quadrature_points); + fe_values->get_function_values(solution, solution_values); + + for (const unsigned int q : fe_values->quadrature_point_indices()) + { + const Point &point = fe_values->quadrature_point(q); + const double error_at_point = + solution_values.at(q) - analytical_solution.value(point); + error_L2_squared += + std::pow(error_at_point, 2) * fe_values->JxW(q); + } + } + } + + return std::sqrt(error_L2_squared); + } + + + + // @sect3{A Convergence Study} + // Finally, we do a convergence study to check that the $L^2$-error decreases + // with the expected rate. We refine the background mesh a few times. In each + // refinement cycle, we solve the problem, compute the error, and add the + // $L^2$-error and the mesh size to a ConvergenceTable. + template + void LaplaceSolver::run() + { + ConvergenceTable convergence_table; + const unsigned int n_refinements = 3; + + make_grid(); + for (unsigned int cycle = 0; cycle <= n_refinements; cycle++) + { + std::cout << "Refinement cycle " << cycle << std::endl; + triangulation.refine_global(1); + setup_discrete_level_set(); + std::cout << "Classifying cells" << std::endl; + mesh_classifier.reclassify(); + distribute_dofs(); + initialize_matrices(); + assemble_system(); + solve(); + if (cycle == 1) + output_results(); + const double error_L2 = compute_L2_error(); + const double cell_side_length = + triangulation.begin_active()->minimum_vertex_distance(); + + convergence_table.add_value("Cycle", cycle); + convergence_table.add_value("Mesh size", cell_side_length); + convergence_table.add_value("L2-Error", error_L2); + + convergence_table.evaluate_convergence_rates( + "L2-Error", ConvergenceTable::reduction_rate_log2); + convergence_table.set_scientific("L2-Error", true); + + std::cout << std::endl; + convergence_table.write_text(std::cout); + std::cout << std::endl; + } + } + +} // namespace Step85 + + + +// @sect3{The main() function} +int main() +{ + const int dim = 2; + + Step85::LaplaceSolver laplace_solver; + laplace_solver.run(); +} -- 2.39.5