From 61be2b17fcb220a5606d15d30dec6ed6298b3cc1 Mon Sep 17 00:00:00 2001 From: Gilles Daviet Date: Mon, 22 Jul 2024 00:15:09 -0700 Subject: [PATCH] warp.fem: New field types, examples and visualization options --- CHANGELOG.md | 4 +- README.md | 4 +- docs/img/examples/fem_magnetostatics.png | Bin 0 -> 40510 bytes docs/modules/fem.rst | 142 ++++- docs/modules/sparse.rst | 1 + warp/examples/fem/example_burgers.py | 11 +- .../fem/example_convection_diffusion.py | 4 +- .../fem/example_convection_diffusion_dg.py | 36 +- .../examples/fem/example_deformed_geometry.py | 6 +- warp/examples/fem/example_diffusion.py | 5 +- warp/examples/fem/example_diffusion_3d.py | 7 +- warp/examples/fem/example_diffusion_mgpu.py | 2 +- warp/examples/fem/example_magnetostatics.py | 190 ++++++ warp/examples/fem/example_mixed_elasticity.py | 13 +- warp/examples/fem/example_navier_stokes.py | 44 +- .../fem/example_nonconforming_contact.py | 290 +++++++++ warp/examples/fem/example_stokes.py | 24 +- warp/examples/fem/example_stokes_transfer.py | 19 +- warp/examples/fem/example_streamlines.py | 68 ++- warp/examples/fem/utils.py | 573 +++++++++++------- warp/fem/__init__.py | 12 +- warp/fem/cache.py | 100 ++- warp/fem/domain.py | 18 +- warp/fem/field/__init__.py | 17 +- warp/fem/field/field.py | 465 +++++++++++++- warp/fem/geometry/deformed_geometry.py | 65 +- warp/fem/integrate.py | 223 +++---- warp/fem/space/basis_space.py | 14 + warp/fem/space/collocated_function_space.py | 2 + warp/fem/space/function_space.py | 9 +- warp/fem/space/shape/cube_shape_function.py | 88 +++ warp/fem/space/shape/shape_function.py | 7 + warp/fem/space/shape/square_shape_function.py | 32 + warp/fem/space/shape/tet_shape_function.py | 9 + .../space/shape/triangle_shape_function.py | 9 + warp/fem/types.py | 7 + warp/tests/test_examples.py | 12 + warp/tests/test_fem.py | 88 +++ 38 files changed, 2090 insertions(+), 530 deletions(-) create mode 100644 docs/img/examples/fem_magnetostatics.png create mode 100644 warp/examples/fem/example_magnetostatics.py create mode 100644 warp/examples/fem/example_nonconforming_contact.py diff --git a/CHANGELOG.md b/CHANGELOG.md index 03e1f4d75..a188be8ab 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -22,7 +22,9 @@ - Support for variable number of nodes per element - Global `wp.fem.lookup()` operator now supports `wp.fem.Tetmesh` and `wp.fem.Trimesh2D` geometries - Simplified defining custom subdomains (`wp.fem.Subdomain`), free-slip boundary conditions - - New `streamlines` example, updated `mixed_elasticity` to use a nonlinear model + - New field types: `wp.fem.UniformField`, `wp.fem.ImplicitField` and `wp.fem.NonconformingField` + - New `streamlines`, `magnetostatics` and `nonconforming_contact` examples, updated `mixed_elasticity` to use a nonlinear model + - Function spaces can now export VTK-compatible cells for visualization - Fixed edge cases with Nanovdb function spaces - Fixed differentiability of `wp.fem.PicQuadrature` w.r.t. positions and measures - Improve error messages for unsupported constructs diff --git a/README.md b/README.md index 3c43f93b7..c05a1f2f1 100644 --- a/README.md +++ b/README.md @@ -199,13 +199,13 @@ Built-in unit tests can be run from the command-line as follows: - + convection diffusion navier stokes burgers - deformed geometry + magnetostatics diff --git a/docs/img/examples/fem_magnetostatics.png b/docs/img/examples/fem_magnetostatics.png new file mode 100644 index 0000000000000000000000000000000000000000..6175b7a47587e3462b7e7638169d36e6b799cf6f GIT binary patch literal 40510 zcmXt=1yEc|w}o*A7+`|C2M=z+-66QUySuw)cC&VLAXBrpH;0mRu_j~vOeX7YN5;m&#zDry$;-~h%fd~jq(CMn zp{kbm!ImHRN&NqQQq|GI&C}S$9Lmzc(cb)%tEr2*xr3{fquXEDK0zobGAJ2w5jC&e zGkteAwXyYp3)x0z&qzFSqz@qop7up$(x}Lv$|O!M%6{UUDA6j>8Wi%tp};AKZftHg z{i^x9#kWlgExX$FRrvhjDtBmiKVc)Q=Vqo`5Q{u)hn(a$xz@nH3dBnY`jz9T7)-|( ze*U2MwyPxJ3(UUWcfMu8XZDi`0#c00cd~DhY*7IaQ8LWQxjFiX2JNtj80uzjH(~m8 zN{4sUMJ(?Ai7lAqJ+YrO!S=*QO*NSH?9a!YB-C#*F|v6+#dZSV;{sS zq)dLMvSjLHTHVgFfT4nhyuC8t39l~U=0YPTEZkaUc8-iuD1$%5!TC9PZ}1(ASRXVv zw?@>aVDoGU;dbbhMYQXO6gRtU#L0-gO=Ix?d$R-QW1x||&)ySioD1L9h=2$;-yQQ- z#PA{_;jX!hXj{NtLPBxye+6R#LSHZV6E}TPwKcCVPnelExEV%&f_c|RtZdJOb;S^e z&rPJo*?>1E9xVA!3P$FR&Nq$qeapK8!po#yc z5Z6ZSm(arccX0Gar{*iyN-O!B8NJGkxhVbDfyD`|dxzE?!W$zCYZz5IRi_Rjw`Pyx zk}^!!F6NeY@@puJ@ZnnMQn6?sF+E`%|XO}w1k;s**k91=q71j2a3u$YK zqAlm=eldlY<$aD&M3GtyQ(EQGdF!K2kE9vzgQDTelhQRs#0g*7RD*;xA*PVn#FFLC=4?rKWr6OejyoG=1~b`Ib4~I)6EtL%NoXZJK2H$E*nMQ# z`8zt-8F?&xP&xi0)(WLi2N?cG*|(|caCwTz%5Wt4PjD=4T%l#P_@xy1i}VEG=^#i5 zu)v+#^=^rG!r#PRAfLJLmJ`wi0{TL`KJgHt?PI`Cx#hwkXJPaOctN>D!;qe77d0`7 zv9CTKIO}bYCf3h(!$w&A)LTdB5QmN+XJKWJhz^IMwt&gh$?hK=Xw0#LgBqolJmMb5 zeMR^u%A{UM)r5R_uQ#-ddiddUs}QvSHfme1qCo!NK%{r(FOZfPN!=2NwX?T$Uv6Yf zw@#*R76nTcdI-YvAon7s#y!rC)91D@o0HkiuTXGsVL&*on-^P)s7f&4>CC0$r$>X5 zUs>)t<}SfbOTGjf5$$3jUr3&D{$3cM@!=B2Unk&v00qVR)|F6IW7ON}flV6vel_LiEg%ZVYw9PQ7AQTHme3W-5$IJM)hx_$+Za7hp2(W|L z>JBL1)OfRP*IAtWSVD@QioZl1!hVpHs5Q}hqx|5Eh`ED%5!zAz#Z>?DLA18`YG#qs z>)`OC(}U|A?e_JBAaKvoV~Q{5-qJt>zbllb=rk#wCCTfE6h$Q*P`xDe@Q~tu&QT0r zX2D)9b|sz=iAh!CagN-5WBev2WATv|Ej*(yKxI#@gAlm`gYJ6=Pm2IrIS#eg*#QzL z>XOyV^pSZdFmeBHu#N&H@Oznv<;B%S1pGoaS>1PO;lv;NrzZA;ujkS|{$zHz@sG}u zCM>M)UH%FwXfV6egG0Yr`A5Q(2y2Td*Sv(HjWjH61ZUg?lBD))<^pj}R3o1@?Cpu| z;NP&=`B3G*I&lH}MMN1Z)#am6KT+-3l(G3?_BJONeoI-Kb2 z@r?FNujXp?I{Fgz?xi3}`HcoC3<6W87rwu0N3`zup;b`>Pf&nJPA~@+h0xWBC-6Ly zy@g}Il%b-Zdx|h;dn@)>Zq($usic_+FQ#~B$071hw#KrHC$!>df^g>nk*-f6TJMGrWh5BOqhbAJD^}xDMOnzKe^j$J-pz^``N&a&JdQne`c{FXR-w1Q48RsV!96a-~Ict8~6t< z*<-@#pER4v!`RT~EIzZ}6`UWUSMcPZnCAY{PFGsj4qrbeY%y%1#@sjh8ZmE!Pr&v7 zDpHV}GqoNG?r6E>j!Vap)j6CPyO9+_0*MC@H1Yj*qxR>|O}&UvXZrS6C$q6vz`8Ly zwTr)GJrmUAnpviuu9`Qpmu;jMskR1Z^Y6^Cr>VpnZv8?8Hfi_$4Mr`DFR)Fozsf0R ztB0Qy`!p}3h&a>YC5v?K`K-W24HzQ_@1AKKNH zJ^p$)(~T-uafpF1zpE(PyEdI=F( ze0r1f=%q@At|vypEueDkt>a-vT(7Q_-l@)@#p?>#{YY!Ul z=Gx^evzW%Y`}518)71mqw^WO9?4?516-3W{pW|QlRTy8KaGucb zU%#|JzSAIy1Ai|oPt*C_i+*&LFeds7+RCa-Fq2gOa8p*4KqknQK{d{g3WI#;#rV`L z{Ogpg*JFpUQTeaBHXxl21^+eqxwRM^PC6qg9F{m>f0H76`TV5%xZhvJ7hA_tOT;H6 z6tWc9jf|ENn*9~Q@AUCk{=8M$c+LhXP z88SE@FQ!O!Xy@b{Gyye&2&_czn11!CAA3OA{n-x|anM5K)k9V^xFgxie44 z`fjv(%5Hz-a~l`P*<;pyc&v1Ydl5yo06$o8fFv-MT&`A)@2ki;sp4zs_E=C6{*S}gV6^CP3o(#zJP?k)EGBm? zAd8Xkj;_O7F!rZSr|#z=meR7`yt98#aM}c~tiL+$@vLhcK5{779xG$srA?bXXwpdA zYnvvnZK6aNODkBzYIDA!xY)TeD$_0btS_c>ebpTGm&A+cxH%j8yq$dZC@Kmfcd5)h zGYtGXUKV4Jcq$TVj7M>i^tGgK_-ol5X$o<5V$z?^RtMLg2x%%33634o>G6+YrTQ}q zkNNL7FEolt>75=gt$*-C!H%hq=DGZz*v&QxE zxJDwLaGoe-@U!O>h|txK6v9i>v*$uuWT3aGk9#H(GI#?~TcAsR``C6S#`c{^*Hw8%`$eY+wb94Iq8JyflKS%%$Hk`*81tP2ehz*39Col=LBHU zu$W%qp>~`UjP|Ka%Yj4`WGE+YY?N}fqjL@$S^j(%$-Mgw2a?Ga>jJDd*q;c`wR@%jvp=R3ki*&bVw3z!Kxy#qO44aSKVW1w)mJj3F`lYKuLi8M%uM_thM(B^&`cnb z--urVE%Pp@+l=iUgDbZ4qmgbRb{Fp_dAbgH|Maj)ICcST>*b&F!-lJAf2%tCp}_}B z0i}P(0ax|>B0z1`U|x$Wz4-+N&wq14nc3xXorf*~5&^oRu4NK2ZADa^ZdT$CC82H? zf?*J83AR$b39-Q}aO%M_8QBJ5D0M_3klJg`rhqA7AZv4l)p~%X{iy*ZA3W`bzrkAy z#A2jnkM0Odmi73@ZsyY&NMl;w+*n|laVrG9wGYZu$Wj9)J?6m6d+KmCznK6_!uW4z zN<$<2*zMuy?NJ+Lm@u_6*=H`qD^`5k*=h;lFm8l+^a#yD21fhofi82b{i}Ae(n;b- z4R?GAiE|T5kipA6^ghQWa!gMs8?_wNWYbKU`sE@bIPc(a@rN5(Nt_9`nF@cf&WdWslG zf|5hSPztGfRsEK;b|FASBNXw2py!~`#Yw+_MD(>W*RsWds8rvT7<|r-BdlP*932fO zxF8eGSRemQPP)^!$hqUh^JuqgF%{s<9}@XUz>gCyA3f}<&C_{+(G8==Vke6?C<%9~ zqbd#(5i>`buivLUcM&A0(t$vBLxJcPJC28H`T+Ce+3OP(@&bLs$+IPg)^Q-rJ1dsMe+81q|6 z-z7>9zDy6fmzC6l(YyHqo*Yi9ZIpahkDA2SnTSFPKG^YPM~Fc}xuLPGKg?T<)x?6z z1tZqk?)2l?`vurSkQ;fZBC=_8;{l6b5Ep1vLV&?-nekUh8y9m#MPGU3&31SoeDH>a zYOih!3MEPO@YF}si?c>xn)H)LxUb=h_O6#7IQkX(_7yW|qp@9%&m7VsLuQj9%$Qo- zBnIOg2`>Kw15Cw(KZelKxhx==9@YWP9nH50<>}N)6KPhqSOXI?FnakHy=8eWXIj26 z)9&Ywmn*nNSKOE9SP}!nQ+UdIA*@bzfmr5hK>Fgz{pAXKW0)7RL*Dyvem1_VgvUV+ z*9_Dg^YtrWDFK7qtCZp!@S1_5Zv1*p45TbC@!6>M8?kR#^k>i12w9qluOdA=SA}DvkNo1*8Rau9sR7O;G zKmiv7&cch^8R%_rZSDN~J?odQSXt?u+5)mKE<@fLCGh<@f;@2Sa+uvc&P5_FtDy)I z${^6gD520N@zFccGyOp1eCRqalM$-g8nZ7?B)sNOi?W={%+$5O|>s^+nMP{f6Ec!Hp zjKqrF#S20Z;YH+U=7dmh^sSq#;@ZqM6PSDZsF6YL@!A^C+RZUfAd0~H?`*A9YSN_YNbFtDnGtOhL{@C&CW!2X`t~3}>(e}xg+>D{YtTcW zK5{IHm+S`m+Xon6h`YRNJja6&2W82PzI+sU3xi;ls_%@KYkgY+Q%{q|K+nAe#%8zm zXooF4VN5obbKfNr;wV6(Y7(be>W7lZC;_3Bt^1ZWy2g*w5KmDvPOHqHAb_F$T1He`pCObk4&nIutUJxO7|5bVv(fVk>y8ucl%EgYiE2#H5EGZWuIl0SLgo743{gm;W#uDDkx)lDl{$6yMC-pqw?yty6eA=B%E{# z<}dS1;o$3j$``?Za@$V0fUxG<%Zrqyl~7h);I%UhOtWTs48#8wg1jb94gVV|zmSc0 z42nfWv`EqIfI26arjgkK;?TA}b|SqRcc!&w!jXnnZYgVnPt;Q+sE);s*XI*Ah&-@?hVqh}H!b7W3BNWg^3Y}k-(+y)cPZdk*YVeu1 z?pfN{ue!GZsy5fPs*3LY&dQ?w%O^HqV!>Qye1os| zdjR_(T797fWOG_H^UJ!ypWv@de~iXXOtr4B?7&3p@4=v*{d6fiIWfd*l?i3_<*@Q> z59whSX>pftYBrbfk~@S~M)sCywz}+a=qI6b3kX(yEm=_X_%A zMC5YtS`*T>-k!{s1H?t?1oc~ROQ<7_M!N^zzl|{rv9x;&68iQ<5A?IXvReLw3vCV! zgX~wN?nbUWnEcbcap(0uIXzs9R?FRTh9|wyyYFY)*pBbZ6$Vo7vS1%PW&u%xa-yVl zZG}qJcgI(2`Ai{jogkXevwF1^(Kb$G7<-o1)4#okydQ}Wm`bgL+GYTTCdKLo`#jxSTP{VhTn(Gnf z?OX|CV^~o{AJbiwl9rh+sy;Mm@{4z5PbVlqN*exdpP`(p9E~gJ4bFdCi_i5eA*&@4 z>iLcyE9ixlY#QRjNg83%a4;0QbKS0WJ9mc2)cjsnMVxJM`TW4@^x!^vdHg}=^n)nm zK?-BSrOLP@RH|>GK6|v+(NsxWE3u9H6EF_F;wJ#Bf{9q>zf5;|=`OjL8Q}GjuD&TP z6_!`lLbriR_GcXk5sOXZeIXrFVrYT~($uGxYE^GPCKxdRx8BI4y5TjfZ+y{h@3JP1 zjZOWkj`BkXls~lEOlWX=`7o19o{#-4jl$5=>vIxoNnfEUp-IWbUVZo#0hs{aQGP7@ zqizx;Wo}yiSY(I!V-jJ-yi%@mG<=Sc!fZ|ljM zJkB>uCj}Dnr_&>M!=3m?#?1%84bs0ax8PiXvo|(lO;Vtmgo)BIf-isf5V~|=EMb45 z9j_B{4iKQ1%^|~-ElQ$?hQfk!Wx&k5?6|!5j0pUO2s6ZQ{MXFdDMDIoRPg>vY zG|}-~O4Xf#8YowqLFjRp11JU41Bs3oBkMW&p_)OqiQD#oxrN)88QUV)+7D>&$57HdM7ddFgZgnft#Cu z-8Ns9I(cufh)R$DYlt%YAghf9oZUOTX8iejLD076R{#3n>Hf>)YzJ!3Pza22B>0Lu z#%Wg+jWPb;<$Y_UVIhN`AMvi2T-nr5dCjnbY)N3xO(35rr>e&n=Mp&4A`3##;C;4- zf)N!U)#7h&=)m0W>BWUDe!UuZExGgg2rCS9$Y}OFiT@39e1}_`OJ3*V^d^CcV_oa% z9w^hjq%R+7qQt3&M%x}tFRFg=T6VM<1q=Ri0NXLCCiFx%M+SE)3)B^$LbI zW^j0{gvH8ENCP)Lzkh82X0k6g_O9AvXk#7@^AFnGvoAQ04GLBL)9U4NKzD~T`Dxzu z6N{e@k3+8F*e`*V%W+uQ_}l(3kO=bC_i4KEGFiqKqjU{QVsbb@9^`o(()PQ5M$WWB z3dC-ScHR{REDxHZ4#g-(#cnju?ffo*Nh;)*tI3ZW-ASPp200iB~6B# zp&t%#0MkXxoME}Q8TLj~VK6D0+$Rtg(3I^9ND#y@8DGzwFPpuDl6sU<0> z+nz7ra6U)ov|oxACp=WYj7I9K2LlrIIpq*IZX{NwW=q9nAzvgm7$Epcb~H^zUS%*)pi+jVOgksRcG%drwQn?QG#jfnp;~G zxp#>ASUckGKNy;NX?%Ny@l%+PeMUPzd>zDj3X{HyF?RE0)0oad8LR|Tjr-oZ~ zo}OJE{Au0nr@9m=hf>laI-ip#!L#(^EZ}IN42kpveK8fL!#i5l zV!d*z$W{YtX2th(yXPkk>nhtm=2+$-7OJ+Ie828H@COfB?bBiGw6Ebpq-3y$&-p|$ zKaQ2(Wy`XO1hZskhAmRQHkLn6AbUS#>JsveN1(rx)eGFi(*3ACb~ofuNP!i~|7h6P z9r`Ub;`{`E1$lfz#t_5UmgaP|j!2H0l>e$^nG@)vh0`SiYND1uBnsdHzjO&yAmQu2o8qK|r`_Q>Sn@yGbL)_e$_&*`@lnKz+=EI51 z#&CC!JuQM9OZCMlY%upCxgU13+AIgV3}+j)hl%AR&WdupUb|?3;W1ONr>+WnyVo5G zFx@7F^5DM^cTthUai)EG^nS%J9nCNq&86k(IPoy$hJz7A7tJCXF;~#yV8ehR?wj3I zs#;X4Jm%};1H`2kb}8{x^oRsdkV?^L(RN))WBZ3!1h0EurqSg&j1QDB4X*y+VPZBo z@a5U};JK+mK(QGe^+Yf3Z2nNUh1Ss8_VZ~SAOB`|a~8GF&kH0*9d1)9!q&zXp4<PvYsd7%Diq^}V;sCC?xB0cs*BrItCIic^u$s(z@#%R*BUAOKKuxBq{GJXYr0cuG3>6S$z|$S$q(9}PZ2SEIVAo*^GdK? zug~KCnD-V+X&O_F)~oX?%pc;zTo?~42>6)ive}4F9=6B{`oT>&n_OYIDvMLQA+q)Q z8`P1nB(C<5YXKxi#pKkNr&sSrlL6|l0Zq5P97gBU6xVabpUl>?_j=|rmee_tf6Ord z|6|3XX65-ANZ~Vo6KQNS92PPO$-@Bisqy~br~Q*FtmR%a40!LuIY;S}viu+JZO+?p zE;eK4jRF5W^2V5kSj^VP0+S4d_=+D8WqLEkw!HQU=jl@M_z^0#6lcz2X(p-tWMZ1C zQEi4;z#Yy^u{D1CzABmOWv+U_3!Rvkyxy88HhcbZ7~P;?ce&~lq91K-_89R=ym;|^ z&_I7%GC$BC_p=9e9w;1xB7;H3M+4jY^flCM)00)4%0v^+j4C5}1;^oIUrj7|G098s z10eRv_WIfi{d1p2T9&9ramfnd$9^NpBE9DN+;USTLI_q^*zchUAOm*r_w>UpMry|? zEvaGRvm3>(T&b@}fw3qM;Z7 ztf2bziI!6VGH3FE>T)u^zX7G?!x-UhC_7lJZp7>>HhV8+Sx1Y{sbxiK3brfv*PAhw zceZC(^02ecEXiVSJ94(P`#kMBrQ|eQ^yKzo(fAHoD&$7TKz{Fsw4Tz+=4@JAlZWH5?vv$-{ZTi&m=yo6fyM=>3DoGlO21HHS9LsprjtA6fUsuZO zDJjyYWy;i2+ z+|Y}bI6VBlF^0;+2ZS~z3t4}1gvJp{(Eg~4=tx6Ch(}pwKw$w79AF;vv>#(^y*NNW z!C8AInBpO1%$gbRWn+51O;_RxQNA#vZkGuM7enbgM$Ixsrj;Q{N?QB$W@i)ZMt=_d zz`2NPE5sWz>SNsYi~!v~u?GIq%Y!ht_iTE{@T{!4RHZxy$;<|3_EYb67dGOn!$$B$ zXfdtV6W@@~6iOl^NqoUJA5@xXS(`rKGdv-5?M>yM68!yr)dBDy(LB9(7+7p&lDfbk z;$8n`XgB^4Fo?l9JLjzeRsOgDvXbOSPMFxkAyOfngOvk22oaCFsAypp1$|{cC|0E% zs43e6<9Nlx`FX*&_m3=~4TL@s?9&JC+2JVujUQVV;zvsCT@&L}r{{cIbw;f~-x!wK z779-q&mR%=)K4koY8QlZ^2X{n_bh^f;#AnP5=FC6z$;X5nL2Q}9N$yUWU+_R=nKSi zZ#l1VFQ1AYoPrm-T+&5UHdy`1K=IwRaj=THp#QCasgSI;QOce!X0Q}~$jHtl43n>J zt0PB(5Ey_+3EGs8>8Pv;J=*`!2Up-AESAyp`s$UH>V0I^E=p*l*{*9zFXjW2!)Cz4BqghY%h= zZU4afR7Sx&Wvdt6u18+buqCpY;3qTTL@ePrAM+WRrl~8cM(hQHBIVpOw{=UZCU;B9 z$SFD`S*;HO;C>x$f;rvYRt$8XEHp#sImGJnX33ca?lI{mj(f;i?qmhE1>K1&>pA2< zyKn=B$+&K!zZ5<43AifIM3+dr%@6Uk)7JKkY*I2nR_El`((|$A^^a%Q-oJC49zDDw zR229Q`ICc3zm$1BpZ8II0m_hslA2-OtG-}53?et-Eoo`j)+VGyuHmLS#p}dxu9;tPa`-wwYFcAgY8xlq8t+C^b&TZnDT1f2nCmsqf z#BDvLA5HoY$GqBdWT_^~UC*jry2#NKYa4vUV$WL&-N2vk; z(Dt_KMe?O4w{*c)W089F%0A=!PD>>MO;X7cM{9rwmZ=4*_~lNok2af31z4!(9kE>C zGZuMh3nn=^oQx*2lhLGk_xmf%>Suz0=c)$%ka(fo(4kr`KT47(1Rm3WgsKgPMigo` zI)OA;WDR6z50Q|N+-wB6PTRJ*)wL|sg$NWY`{T2dgkp#QFJ#X2$#@8trYOnw`U&s` z`*#GrTz+Z2T+S(TETN$0-)#EFP>QAH`DyUU@ApJZ@d0h!fwRP@}m9W~pc;UV${5cCIzY20W& z(=yxj2iIy|o}7&Nh!fdx#uKmC8Qrdu-F7kP(@H2GogTtTSu#+S9c{PJLfH)V#ok{0N#g7P{4EdoY%Z0{OFncs#nKtr()LYZb=?fZqIv zVTEOaH7KHp)=38~KG?_-fC+LGVJaCMMSuUAa$uH?pA2ou`}`EqFebsR*18ib(YJGZ z&j3(wa~)*Wv4lwk8Ie+kEH5O2-=!mjxC{q#&G*2K$lGz#66VUV`TqXMkGxD?%Jhq? ztxO0_`S@r;u#IF}@wYkVR+aW+HDid@smfU4GJB!Uf?Q*T`f@9}B~R@`3mR>n{mogs zfH5s)V{)8l=_iOj+sBIpXnF00|`Mv#Ms5&i3@Ua>`Ae$#5H;Fv}`cB+v8QA3WE@ByTaZB zz#VU$Mnik*Kz&220S@PX+;o~=T0Y3O3Cik-qY)oe}IilV`Q=zIDR88`( zduQ>*_O~d3)NgI^;qmyRGx%|s>%TcB08H0{2|rY|Pv02(Wx7`WKVq#i!{shP=G12V z3sZPnD7>CSSqphvF2<1L5t@1}opq2rL<|-yKNeexz1*CSJ69b? zxW^0U|8}QwS^XxU+PzK_DQD&Aq|}{Zsb@*6n(~30GxK-W^kROl;dr+v^1Rk0S0Fut zCZD)2qapWcO$S=xDoo5nOI=tU#_>EsD!32{%m)D*V(OIlaB#NWHH-Qn~^i7ob&?@AN5)*%f zwTqRvy=9Ie6fssr+zlgq?fR@7tY96A>jSPrk+y{<=Cc|19= z!K2W z{M20;90k83Jnft~E0`R*e8yw@jJ1+b)|fB(_CXM+;YxuSpaBF8m=7({%1J5_A}Me3 zE*iPo@PMnGr)ea^cOV@jx3J?^4_YMy0PiT#zM~2ETu93SZg4Kz zo@q9u+`6eN1Ut(Etn#6$ScY-q1OCksxJp>z_=)J^u40JDuAR1NMEW&X;~_<>O~)7v zD5A<(PvO2ItJETC(2df7pWNyfP_w>~UhfUVTpvKxA!AM^p!Dp*g9Z4t(d@o;VrBDR zaxrnhLur&hYS!v^7bpNbC`z$Ntz>FZG$?9Lq&?tkUeHBcrvo=f7Cy(*ve^qf2UKJNc1$IncDMJG3c<6Refn&5GzyLn+Vs7%OR>NI%wj+1&)M(Fko+Ap88q9 z6MQ3=m-B(HOZ`?ti)yWfb>)1%mS$68YRy8wO`7=#5;BQG;vP++_&7Xwx{Lv+?#E5U z&e47K<`YIH!S38hiIbD5U83TL7%}5Nzvp;5_F@d&StBP0qNFIC52PDM`2e1$*iv!_ z7`x+9VW0rG62P@YRO8DXSMuKe`v2Ub64m#{grT^=NTNe#==DbC;*GOh8bO7il<}g5 z{^uVXQH)M4%bWP?hVz}w1#c?S8Ozbk*I^l*nLhoM%5a;2`10Vh7c4PX=<((wc#hZfZO|ssTn4lAr1f3Kk zfOI)HqkxVSA>lM}m_1A9d4(+LdV_$+d)2)waYn0O>WC#UkCgeSFgPdvDZgtw@6&~?~_>cjD zr1N~42n{ij2|xS>5kuSa>QyLL3?sHDWO5HkSVDk~w#d+iBp?4Dk_kY#dgb1ei|cG+e;VCXM{z|Lb! zD4Tb0L}=l2DPoYo{LMlHsC<9o-5tB|5b!JCzo=tl+L)I z-a3lbY;z@Y+WNIL&jnB!$-f~Qv_Sn-nggeaqJBwt%W4Bih2b$XG9RDzavRP6KohLE zEaD}5XWUYa{^1S0-hs}AS6E?X4Q#7yNg)9okCvFxk&;ThXR9!Tx#T*K+_`HZelA~j z%&}`R$67gcI+=|$sg|{}wF;IIIg4bZ{PP$!R2(?!b7d{=ph7G~X-zvc4cf+{QAN2B zE9eUMhNEV^SH$5Py>*!UV1Ouq%Y*T^NpTa)LyX;g*LBvjOo7x5T=>@cD}Zo?Ng$HF z)s%N*LN4MHG*}6=O!xst^}p=V*sYWImy;FOXz<0j!M$w zr?(_O9O^~v)M%K^uC?r2#77OgVlW-hRyUXi>Ln6m)<=LNxPf)+@U z->R)7jlw+J4p$T_Ovi_nLSfRFNacox#}H3@2-;hV@d41`b+Cs6)KcCJPjgRYbtn0W z*(T^a>zBO%G*OU)0U(vg6|@0SsZoGGO0uDWAVbl1C=akre zpLs=V)k-+}y5fv%EWqy&PqVGX#oG-M&s7#fwG6JJ7Lm51lww-(1mxlPChg6Dp>-x& zU~c22c$;E`1L2><4x{|4Ne&q*58GXf=Ex(=09Z1EHygqKcEewb%2DoyX3j`mQmOrl z20}r^FJOtMH2!5zr=U7?qMLQ>Gu=Q8wrEjHEfZA*2CCuOM4qV5`-qJ zV!(9^uh|#a?4KrR1Lo9al^n zhJGB|Jd6>~1zYl(?2?x!h1jVW(a*w6Wnf|`?m1MO2+LTbpzM^mwzM#dG)nE4YUbJY z1L!b}Q^bH>Jz4b~nG`#!=8-~U<=&Uo$^_()KWp*{g>VTSU-L$L@`ch*Np%Ctx}yB< z4J{Y|%Y%n@|D)_$@;3)%K?9D4L}Jc=3SO)<)A=lUjv&QIZ{WZ!6sw$te0y%D3M*7$3`tAxg6hyi3lSAaN~xn`2a{k& z0Xz_a&->+Wz_<)kq>7v!kAyNXynrP1lgbd- zKf!f`8iVMy?bs47TAi*MrHnh@oPxhOzb5A?fjiVv^oe*{h34!3G+`8F{dY1@{-;{) z@o|5LK;@xOAdG{5-sQ~l$ICero)o#X|H$Nmo%~A45n>f!+8W?)envapu?WprnF%ke z{hTXsF6?mC4!9Zlv!h0^is+)Cx?K6OZ&J60J!xL7xLa>CPR7m#b~Uz%f_o=RTKI81cLRak{iDAdk`X0iGT=nWU;Ww1KKCessRWL=BMZb+8EKtonBgI z`}x&#eSSOW7|~6Em)@4z&vwM0lnlNJNDeZu5(Bj%JgYO@9STp-8>$ggj+1dUlm@PZ zqh3R-c5=^^)k1^nE3*wRx(ZhK0H&DpUTVM!rkt*8`x8%b?njc59}6Ua!cY?f)60qu z9HOC}Ud~-7#xLzC$JQf3oMT5EoJjCrarkBxk(Sr7+{NEqL{*Ip5c32N3c%@9;B4nm z1!_M{SNZ_ZlMCZcswHUzMt-z0+h$g)4wr9IMuGJz@3>`COI13bzup(7pBttS$EZms)_>r0xTf!65p`G@D}{s z4SS`es()`mEP$y*XJo!B8u?6lXn|Cn{b&#dzA;6uJhxoy5;ounrds--Suz*?{afdSE(^=Yjoch()?XIqaL2!hLKxJ_mQyj$OYN zydMb`FDQ(>S>xX!*S>7a^1>b$078hc$3+w>Dq|HFza|qH77r{SqTh;4w207$l-7$|@Rhdq%CIGfl?Oi)T#9fm zB_^YR4cswLUk9MLB>TMdg2>zPb3#Lf#!eJgp`b>7EnOlx8J{;5b)q(B&z+=vOrZxj z)}f3MQKk&KBDE@G>$2@=6{l0UxsOK#r$0OyRI6_PoO3oof?$6Bpm2M7GrwVP5@P#A;y^9?YBS+1F9dL9`*xtJ-qS|-=C}6x>xs)8y*<0^1(wkgJ-Swe@ z7CE3OpLSyL@4~_E-?-KYG->f{*9VjvQQ(x(B4aD9lqcUf@4*q42Pre+kqbr*V~tcF z?KH=ATTv37vhEtFfMkNS7F86$?UxJ708?bO>sLc(=_v{t{Y=$pr{=vV#n-0J4;!&l zVHjY%iJ1s&3-zI~S((%0Ut1eU9cu8pILoT!n!yklRoIu|IB)*EbT_EN=!;G5IyJGv zk1khSfE)shdCG<*m3H5oxnLF8t-5Ip(>uzLYV4+!RqHlibBbwiV=qNXP-5Mf8EfQ? zZJ31j1O$U4M=t)y0vO#hEdxiymBxNsknk>|YhPvSa23Y7Pb91q5RYz5U;(z`S`LT$ zyvNM4(zseu!*DieDsC71C}aY#eDM&WsEAdf!c};DtULn%=E&I22CPf`rWLPVOORkc z-Et0KP3Ar9rHu6mMME#XY3}P=BLqX#ks)J;CmVM`brO1++ToB-_)0NNwmYO%d6#m? zC%b!phffx1Gw4Ms7a^HPv!gSIE`=6&0w(PawD<&nG@i_baV_OIzEj2TInh4CcUX@T z0#xGwtbbCPTDj{7N%3b&%3ksMxg^Ko(03^G{q2+*;LjC*Ng}LHR0q^&+KDQY=>$B; zZhTA#uLkPUCU7`m_xC|Ry)Ez$1cDu>#B5=^mNrJlluL*aq1Fia$XM8yf6Q~Nw1l#N z<)>SkFuM8#ql)f`7hK>71;8jcw4V+`v_ERY4I$6=J--$4;b3}2iaKE^1>vZ>4FGftHCzXE4GN>F zkt04%9~hq)9Z0#*^gTO0k9Ki@JijAMx?f?1!J9EAu?<1TD%WfXfCY}w z#M{vj6~>}AdNe>qR#T?7wnJa^D(cj(g+bK*tPE_7cXP#j&;#odFuy9Srw?ht|#Slr#+ z-Q5!0ZE+{K1`QV6-QDxe^HzOTyMJb@rgr9b_dR|3bYH}9bEC4MRzwI!`Rt|_JGbaP z_SSiS=!I?*do4>|Y-fh}@7#6SpkB_c?`R4(oRKRf)yRFT(}-z9afD_U!e~Id4@LLH zxg>?wcF19im@^Y(XjmEpAl_BeD;2~-hYovU`|ATL8!{-tM z3WQ5vl&b!~7?yqmY4^A2OgsoeZrO^5z?5$4r*hYBcNw)!g3t70gt5u=+U$^2pRXmOfyG{&k{*fWG&mbJSPRmp&R>$Ll@= zmjKO9icE?Y11UinAV$1-gOHNF@nXl{;VxD@m1U(hUI`{_ zu^_2sBB==6>cod(=q}YaUSME$c$Bbp;2w>F-t<`z~GaNPhP$R{DTeLGl2$g(y>r>masG;Gc_9I4Ex}Z9`T* zL=Zn9E7s=<#z6o(x_4rxkgwTnsLn>u_Z84;o0t*7?bCcE(N>5c;Ri*(4yqcgVXINj zf?TM!EI?ehVX3H&N0u#=ntSal+0`>Ap%QJv4iLlmc9M_GY3yRCls5!u{USYbRbp;3 zc-gk@<&xTzK}HQ2l!A8p$w}h_IjJD&uS%BU_A`ltTalp({Uo0pL|ccV0@=C(snyTC z$;(l#<#Rk&L>=b(7mxFkjbuiMK&MSIgDaVtpwiFjqu8GZX|u{7VPIgJ98wi1zm7UC ztT##|+mTG-YOe;wQ|YaC6lhEoU@=)Z)KhZP%`4he+DnKNuDcuhFDWWs{CQ*zJy=2G zV`y>JOM!J_*oG~>Fm}|Yj(Z=tMDN|bj>(q&VYgl=)SZszvrIWeDx0N1xc9k8+hi_< zosk7mP8zzFgiiTu@ggJW9P6McTH$kHO%x^hu{qb0NBHvg^%(7PNA8vWiW|28j;osW zch<1*!8?4eui{|YyO$;-af)=oy^}^2+!{iK^^ci843UO|( zAaw$efx;q(n@JFaizw8n^=N&D+k4Qy`o9X=i_u{z40vGEF`)nS@MjOwL%kIzL@@ zDtZ?d?Pa5`B*8f!DRTWJbn340uhA(O8J8Uj;?!c>-P3dbVBG*s(!e~m5*{rlW~mqF zsuA4e@%5JQ9dzzUZjdj&1agsJOcSLt7I$lLA7GyGfrL=*lu6C4Wur!crlchO9EFWQ zC9FoseP{h`8QdUA86+$pV$bB)8pm_G-PGyTEQ2yYF8p_%s2M64TU5DIuvmzOg^ps- z2@h(3JyPLwK@G=Uqjw7Kxf1U8;))|}8U%6eV~2h$mRjGW=%vh zm}=5Vw;%X|R1x0Q&3(bk)J5c&yhK|@=dO{O=zJte$dRc@fmHpJJBPlY7#yAW5axb0 zN4O%mbzHn1hMsIi))3t)$=1oZW&g`rqq%VQAF^gxmiLUVdit^*^BC5jPgQlOH;`6C zn)(pCtl_itk(!lZI~NH0CBS#Bf3_cAC;&e>?{@JzjE&Xf6lG-7{1;|t*NU7i53%L)N|7Y%xr@-znaiHZY8Zkcz6Eag%&sUJf^$(fttRC{ zoex=*MsDUUhG#8bjf~UktWX+l$%frT|F)lE}N~6 zX#V1hYgnY!x1Wz3$to|wVt&@%)7;&;bhn3$AhXrfF3uHMGr+3LaeDjICyPG{tP2Iq zIxZc`FVPUbj`c)Y8r4V6AYz`fquLET_Z%W`NJyuZb@mw3nJ@N+;3zp!3?BjVg;asY zk1d3t$kPSlm0|4 z1$=tSvk{Krqw;Hi^3ev$6uF9$%?^6hvv=G?4rx4N_v6BOM^~9rF)jfdeQYH&u0N9{ z`B#XDkvxHFn&gSsK%BbbZOQ5wO43Fqft+sKjA=7~2C4`s@i4*2 za(E?Z0NK%wjo$r0B=*|FmTj{FDmp^~*5I4VvylHNRb=2jUs@9(GkNAQo{1PDQly4z z$WvEPBprhk!?)46_^>=j*_o5X%0nC0ba0~FAnf3DJvrv0yF7KQQ9~jmG+@jkG#LfF zm7vyeR`L`bbeoJK51B;*;P3lT5_VPK%#Z`dWSg2;lNmTjjX$8N0?|hvrd|`6qppoxUGl%yQ43w6TBFI^al%qwxaFmbmrdWvflSHs38sYFVD&eUOPT)|~ zDzQP|%#TEQnttR=#GBi>w!%AX)IwfsaZ2SyWq81k=!?O{cM+O(rPhh%_$1sWq-#w_ z(^=^$k(*z9?3Qhy#9U(29W~1LRIO5j|GHLaNQzok1)+iq?OH)#r^CNG@(Vp(4158}h`{?{^aKJh)GY~1 zoN|0FB1D`9o~d`OR&~60wn+q9)^u)TQgmg(vbk7)2x^}d@5G)LE@t<03`3fe_lWj0C=E*ar zmTo0F&g{9qwHBy{hW|8akEG-mCl3_`8b2ttMkyJWj!tc7w{Br!2n3W#7-fdSj3{j< zv}DTj_&vIiP@74DF);}TH|$ls#3p_ynR@fANMAGp($s!>nDh&QbGtCIWNN$M4$F9yw`d0nbI6(YH4=x3T4TWFQ) zOj|hQkdZv69LD0bg^Ic;r!tkO5Mgvx+a0F6q^pL6YHhswqFh*zYNMt{j-k`gxqRAa z39InoAV){K05A&q*@Rn*+klnFqd(o2MSZ}8POgCVMh@}5Hd-h-Wjt4Xd6DF}VXu4| zREQKG6s6P*?*TF}PK*u;GM!MLoEY-d!Y72c%V8qzv-WByCW50X&?ZT}ig&UE&Wx-45pPC;Z-yf-wU}!Tc-Q+8?3TLY-nw|)MFrNKJfZz+L#sY!_nqAo^ z!m4KZKCf*EaxK~IGR11xKmz25-R9|l1TlU@&zpTBN?m&95iv?*Q_N62j1+0um$J54GR*6wR9sMPnQ60-0oNyvu z%ebg#Bdl5hJFFu*689w-lIA06^-e#eD_tuwBSOc0CrKWgC!Fs}g7>4)8WrWE6kwUR zTtwCcNoU<_2;hatA-4_;xD+qwmK;UGk&kMnrP_!06^r>sM=`1zoUikJRi!v88^xJZ zm_>Q9uUTZ=zJH+_r-c>`bWHqA>4k#@fqy+m+^vq-$=e{7&swIHI8$5LhS`j_>Ilr* zd{-XJ(q+ejV~it~q6edkN#7KXdz4tVV zQTMGvGY}F6{$$it9r#$Y28Bb6z>?TPA~M_;=MO5dOFiveY0a% ze~Ic?>A+-BXC_3m4Gc5LC$8w6kh8nt%tpoMXTfVP;sozK5C}ouwEL-+nS!4&6- z9#PJa{~8^{;$K`am+<=qdXR2JQ)(u`!~_ow5|zhs%6Q=PnLG8LUYo48j+E<))~E&U zA~b(5f82;8v~Gt(Wgaou>A@loG}6)TSL|D0#BucQpTtM(zqF`$zay4OOv!_)`5x9a zdn~lUXj4XID%QvxVIT9UbUDxCWNOh~Qih(KgoWoP96-dhC`TOa?KwIwuc^|T%_Hb| zWJrLM*J$T9Ot*R1xu#eP@Sf>i~3e5 zn91$#4Jf)*G$1)>B*U&chG}KKwk!MOtBx0_j?GfAXCWb&8>U&qa_g4Fs7pA^T_dfu z__^UjfKMZK`8Mg&*^*=vGwr26)g1pK+euc`a1K7whL%?R9XwSYq|f@vmFJVXQejPr zKsVB2?dW(9zc@zCT%k1$K+Fww2opGEHs9JDPCt? z?wQqH3VITjq|$(lQizh)%Nu?Qbp?`JmS7!=a8R%!CieI!scZ?BC8{;!X^SF7$R^%a zmJUkbZ_t=5QEo4jgNkB8fim(*0+1Ms-73{4ij`@X0@v(|TjUGKTQ%_K^G!ZFL1W}`4NNtg+X1E|d=`o6hq;5;$k9L+1<; zB*7iqHdCB`KA3&XFKt4E=b4qesf;0lgkIIm!cvcg zP93PqXJ*0_k=qGjP@h<1a(ZPNu>=^jLtLJ0dCuPx?EvkB$)Am9|61?jte2RNT<3Za z&g;^k%nCHAo14MdME|y|#o;N~!HTmidtw(MOAyZ)z3#bkVm4LZhO3rMGMXu%iG{1z z(U1I$*ea+eKq5;pc*xwTJr_F9?9{WH4u_0FYbV*Hv~O{mn8>>JNQTyIHYq!?|Az-t z9N@|##jqFs`Z4Fz{`@*@=-h5<$Dk7>SeRXgb#cZ=HRd|qzI~DKg>orGB*E6}f!~IS z1mPE>7c?G6k^DLN2BfQ9D^*Q>*MKO3=hd4$O)fDLHzZuX<2i$ynZPYIDC3K)VodnL zN#O(F0Jh@TD|)t=QuDsd$YGIiQB63E6i}#h3m_S;Fr@}~JOx5bm01WJJi^3J=CbPR zLz=+y4Vcd&QqT+N#C6s0Ec6&})-hU(nv>=cE?y$%m()U6@cU96weVqHHEN*F&6-NN z1}Ew6ju=E&hb@o;RPh1zrLS;`^U+fOuPN{(XLBKzS$eGaEC>fSya%E(TMy!mt2MDsQeD^mf(MGly%k-cmzk^z%b97I=7)@F9o0F$~>De1Lz&d{&m~oP$ z^<48;qYMOKTRB$#x|R+0l?%1)F*3vx6!RS=>enQ3s-%MMZ@EO36|=$1lHORg6}8B; zjX$-5jgP?q;J3sKM*GDUxDve>iRE;yyN!D1Jd3^dqh>OF`ky7_m5LQ757`O{ul!qye zgzB5a43V42As=iYjyKFsdcR*?b*i7X2by8R<0pID{cn%)!}~_Z{?ZHia#VYnC`Se_ z?UrdmFnx12-KQydioFUPrJzz*)0Mr<{#v;w<{XLuGr0dChKVqsiv{UwO+!o>r%F6h zP%=x&=ha$UiE5QqFC=DHaT&}nEbls(S+t5;H^MfqwH6YCN0`Jx%}^LDO;}ChUfv|n z$f+#xjYz%W$mE)kNiRv5I;jJbd#Q=)OUBfq;a9-Hyb2CVwtgJOkBR<0KG_&5pM^#t?+G|c27x5;W>XZxAXpdm^=Yp)T=HOzc6DnJ^jAVO@i zrpcv*+HrafN%4R!)c$5RIa+*$qm|K0?$Go^e|6BIJV=w%zX-aIU7jG~J!;@Ke{L5* z8CM8boT}ri<QgACDK zF&I-fgs}OY(&!+R9tMiJd8d`F_%nmF^4RdKfz-SYQEpyrOlp?j7xrNo1$5?lZw^F7 zS`q<<=N}(ea8Mw*GO8fun!o{;TzXgjZm|`jAli}`n#s|Nr4Ed%Mn0!}VhLDUr~9qt zpICeGipZ<<`|h5gW3b7L*$Wx*$o5&E`ku<33hjf}TRE$}B%E1O-s*+4f*pJ^LlPfA zk+OD1$TXlq_}7u4;^}!EV*ARM=ms4<(Q_@qxuBv97>zC~rvAjn4IGTs1~%yvjA2sK z{jz+ZfnLiQ&;&Y5mlV|&hFJ{DQKJk`5k{x@D0n_*eqn~C0fk~g3@c9d5d?Acj}e&D z&m;B4w#(DV=$aE}(|0of#P1g-G#1L~HE0!b@=7s$Noj<#uT0xFPcU?o9mL-y5o6OD zT?HWEWY04sPi;@b2QCiyFZ{&#f*T{4*u5G3a@`H*Q+Y$EVCti`U1;nhB9>}S++JrK z>v4b;#N-B^j5H%dK)0)V3@S=PJkE&wTQdctve~uMcc7$fR* zLd?j@i#|z~?wi#E{+XL2)Cp5)b1a2!#PZ$|hDZ@S;SVG${wR9EuLFWS#y3XTq}#6) zy9ayQ%QXL_P~-?WfE2Twaa(W|D+)$Jx^#Q=C)PM`NKkO@8#b|^sxl6E5c^Z*__GAZ zD^+fwV`!M?k~aoYJ9bF2@fMz?YE_MW$VPr?gUWYit6=52gQVjhB5HBf6j6mj1I=kQ zK;9rbh9KhEY~({hJ4r;VGsoii|5|{PX^6h)HS|1#ksJljdN3(Rl(aP^cux|iA)-uu zH7q7U%lPWZ_4~n3jNqN!H2lvzDJ#~ie@(@l&GiwA=*cO`fv2@_g&TimUTn&RQ8Wi>IA9*UU9fPUEbbba|( zE{I$vdI@KDsLQ{Sd9I_&VwMhPCzHl9Lq{Bs;c|s+^a1KT$`?K4*}X)ZHuUn!nrg4i z!_!j|e~)#4O0cFVa9HcIDxWc+J^zxf>(3-;=W9#bm|6vf3&!=r#^%rCXIiT)=!}lJeu`sL7xA|RUlP))rqJy&!&d6 zV0SsALmzmy0IhoDT;iJ!XLxk%potEej0c{{N$?x^QPZ~)pyWxCVQ8TrhuVXkftZDc zxbIBd*B<;i_U9BHxFwvGc+EL`1-jn}=9^U9=x>#MlG|4(%bLslsU4IjXnOwdQQIl- zw;L7Z|FN*3b8w9L=qnWiCmO(bpU3s@6J-sk&&=9Nnc)1=pcdwagApC_5}uri-8TUT z{5i(h2D)-xY6&9za+D*K7Ki9^{)^jYn7=;jFs2yw|JEn|YoN%T2A~OA#p$6WzmhH! z=^|XJVe>$oE5tyBuWi+6mcme*14$%dnN3EXqZUL}EBUHpNh=O#ZOz?8bsN8s zNFy>jAbIs|gBYene;E8#{fp#M7=;G`6+9Hvf+qOQQxaCP+ zKnuVc>DVGqjTM{4k6reFdmArjh6VHt62pr8++~^nT_Crm9dFTZhN>XeJ6+hKzx;UI zm-_s?{PFJeGh~7tWq2At7@5un#^s9%t=abPdif1Mxniuf+fi{#Ni1Vx3P;aw_Hwqo zb73qX8Ckn|&6@xGNcAeu0rD>QlHv>gL8M~xaHAG{NF8dlnckTB*-+o-Nyy*#VH`J? zr96_J|LDxUsJ07MVgwE8`QNsk1&8mP(RcpQsk(&xCr>l>U6pH4f~{c0o7~^|&&frr z)1ad#=9BvI+xtkRsPLc3iP$(8;UIQ^=2HHO+pz30JUuZo_C>0Y=efzO$jSE5i}cHv z0d{;5UUZtWN{JQ!AS0Er;^${RFn`mV4kd`)X@i}3RjRD@I7ZIpR#k5{`aqy3Tme9~ z(q)onh`(r2Am$p33l(HCYP|)JzCWS1p(!pVz z>)4%IG7#?~!^Gitji6QowisZr)$v0ho+vN~+^JwwK4TT%cv7@{l}uYxMvuSXOVC-0 z-SzXIwfmC@(C?=gc?qfh!0;SrzMEYgMT6nTB>wn;F#|ky;GmEbzuYgX(Mf+IRP$JD zqd9>9#>h7R^-vZ8bT#VArP-{OCFz~W!v{$$Gb77D8LJX64dg>obd{cJRdUGTyKl2G zUBdTS6Vki=*W=(usyk^^+R4cQl~d&_U;u53@$CI>{q8$qDmw4pdSWzN{t$kqkqGABb6>gHv+w6I><5N42dN32UcXLooL6#C zj94oWQY6WH7L5iONwGsFDi_pHhHmZr)kP`V(+dO;9iSi3cXpNKVj@PF)ddt2J8vL= zu1a|MI4G*JgRQ69Ta6U6;Tw)+U}_CKfIH*bQjZeRB*>Tda5Avs23DUBPXTtC-R4*} zS4dmgj*66%pXf);dO1egUhDpS+&wh4jT}a*6o4BjgV$lVlxl zC6E&r9m;@mkM>Ot!}9ey#}hp#-3uk_-%x-UT5b&^dseT+X*+0PczQ^A*AU~Bv|@Ml zM2bfSX2fF{>0n@}q4r}TcG;Sa^#|p@7&D!%lg8qwY`ltyBg9{@Nh~H(&vsPyo21PU z&W~c#Pai>!I4^%ps{JV%B`kO7X@ifdZz0mI z`Qf#~W;bN+%GX%Am?v6tA|t1_vfIoqOX~aY@s~mgU5I~!pR*$?M`@Yy!|mn` z2Z7&OzdBhS!dfT1HNg&mjvcr9j|o3KXi#uj$>O>NSrHy0XSQw|XTYD!C4CbO2L^Y5 z(3CS}Jhih;F~(UplHj+BRKNFv7~HG{sG80<$ws1A%feY$tSZk%O}-HgPnTZq_|2*> zMMd`uC>jb7q71;~Y>8%gshA@kD8!JVdC|j10mKI9aA#8K;dGYdf!U*}kOy(&Nks zw?HeNi38#5-*#p#cN{osgOmxdV}M|-{m%??b1;Qn+vybB_Tz<~+bJ`N1ioEA{00j| z8S`P;;3~oPEtH_wYu~QA7EN$Gnvlaw+7iogIw1@Y7l2gb9R<&V><)Y3_+m3hRl4p0 zVqQ-Z&uGK#Q+~mflw`=fuyX84X3nEK!18$et{XHKD$ua&#YSFPuUv0Hjb9x^w%ga_!{413y1 zc-I8!svjT%{fJaR%WYy5w9dWLJ7;~L-=5zEuMX1b#8*fVD761wlQMU>G2Xw?Jz`S` zPsVDD`=$S^3>T*vAs^ShdWGWq_)1aAOu*7%CB;l$+D1QEmg7}&Vw^Z;mg#oxdG@^B z=BaFYY^j0EozDQ3o-CV~iW#F3{bhv=$e~t(eSu2GwgH~!$guL1w6isezqU(e8N<9a z>_wLkt?;WSnb6Ox9$ASL*er7aHBoc~%n$zGgs3MSHJCQc@l;R73y{C{_=x5FWD&9n zI~?Cs`D@$V8n*h&neD}A$E=xu*+pAA01}G`e<3`svO_dvAC>e6d)x1B5L)XTRC_A` z)}dmQj_9NMlYgWo{%3I}>{e!+21Y2pv^(EO^z)n^ zNumX*yDjjCyRRx}SI+e66)>%vyznv)X@#kuBOLBBQ3f;JmT}AjRQbeg)Z6Q5vqe-SGFf!On>R54V>c2KZpWMP-9m98W1PtM4!cuD41`cKXTZyu{ z$VH;@Iv6-;8mGdgk-e}R(SyJbiHAmZ{I(G?d%}a*EvE?49`A3@JpFV zVR6DjJ4=f(*Ky&)Rkm*lJX*BX+zLo{4`)+CeEj*sG5~QkRh^>!kz>zaXLv~+{8u_j_*^oUCwld{( zM#q%~(W}sJn=Gt#QD{{r`#d&@aZH0m^ff{l20l}|%)pWwFoG~@=6eq?@c>swJ7Wuz zrx^*A)9GGf=HColU!C4XM`z1dexM26I_gm;UiKpO_|n5yZQ{Qzi6cL~XtIB45urp5 z_Xg>Ed!fwS=_gS5q84@#t2mN!xh}T9JSSy!^(6R9%9qjbwq?<6#~4In{(S%)h0m6i z_^&%IvcnJQ?mUSc6)H{t4hz8X$R`EVyu$E`sBI zGsoKaXN0H!POsd0-VJyc5wpASAk$x!e@WX6zIVO9!RhRhiwTO0n@QbGadnxG>A z6coh7<oe6Ks5*>I?hRGor=7>*+e z#iUZd-n{R7;fSYNZn~gGJ7?)zhpOxm>G#qYNWoHFo1@J{=lAbI!j1XFyTOhLoZ__0 z=X<|^>B)!2yTOx2nB?5)>-icMm&s0r$H`Da=+{?IEbL>(lng&)e693pAUhUVRXR3% zNg%O8!8?K=rbnQ7+V1|5>9Q_)GGA#Th=fVvAo+mMGkm~@=CedHva@O-DP;cA07Fe9 z2nldU*tn_J@u+!1_&_#p@sM%F=ydp91QQqpZ&`^83d)dk8jwJ)w^M?wTk94v#wrBg zpu2z%ChaH)Wf;@8N?#8`P#pou(mM2b9YFC!aYh(jEEHY8{Hu{MF|d&Ag9ry*t*(}W z#&$5f|NXH{^V~@qicH{^+zWuwzRsg!VMJn$!PCTBM<);#!ch@!#|FgY#4=!~T`TOL zbtB!6lc-tevmy35$d_Gh1Zo@t`ToGktG^kwN@{@bZfT{0&{k0j?E1Wy z(XnoD!8qMRbXDyCq22g3st?GsKGkg}7Ht zI91k;5fVi&l<=>s{+LDU#QZb*m$Sn*J@}LS0)lJAm+cWPEBd5gxBg zvZIr@@Q_VJPvRo)kX<2WQu;o6M;!`HJP)ww?Vn#=Rz#Q^I?JRb^h^f-eQL@zUy|7u;HV~0J^94uKn<* z&xLb|0PSq33;e{^?NbuBP^uiIfH2a&h4nVLBk{}UIFIU@@!Y8@CknJnA?AT3m9Vcp^qdP8CE**E^Ti-N3Tbkv zgRn@bausrj%O^>tzh7YIV}!0D_b z8(m>1kNT{O0Zl~e?bs6+xnp|pBucox*E`Z{xBd5#g^GXds%?i}0xtWz**36e6;m+u zYbF1yYbOb@AChU<&1kGkU*z4Y%xSC-3og#uQlPW=Mrq*jMqZe*MnnktoW7baJ&yEW zc-og`bh-Y+GhJZHK_oJgh<75E0U$`S4_>!4njIVU8U#1wh*NRnfbxc?6|&de2b-Lg z(=N^ut%KpR3MuO{#;gUAUd>)uYiqL4M^pCN>ruvDLP8C#p|62xu5!=t1#-ICK#xq` zk_t}^D|x`nWJGaA3svFk_J1bem@MRduUFh$aN5Ati~&>=N4q({=j&oIPb4M$e2%a7 z^NZPX@OGYtZ7?}6nCN(5!^KE98}m`|m*Vm5($l_QeSjSuu$-}{r7KLx`&`g;A{z?ph)yyI*oU0=RQcQg zU8sy@e}k`E7B4?;O7!7<&!~PnKyE0jP6lJAmHF`A72NAB{W`ykH@%bu zf7mT+`{N<9P7vdb@)I&bHao-L>!%09Rd!=*v^zgaTA*yfLWe2%O8})?Ehg1pNM5(9{CZoFucsv)B@pt4}d-WW~Sy^cRxP)_yNm5jn~Wr~t4V zJM;Xfc@tQ%9VO1jZkdUH?6Ts>6k^dP&tPURHqYlX*bEK0eg zTXM@E-CK_a6pW2H!IYtY+v+~#F*>^=UVB~^KSJ5jagGN#e#^enZ}8xkI6KdmkW~(P zMK@b?Us92c_gh{do5U9iR@Rn5=)1T~C~W*2okOPg+h-bJTORY?kQD6mg4SjUFbzC< z+pyG!orrZF`;>21{PNU2=95#t?S1?Wpac3bJgQfB$p16e+!Z})Z&aS%M>Ap9Q z_n~G3mv44)27C23|D@dgJyN~9R(Gf`KkP`n-b&GBpPyfttg-abg#1O?A9lG-2V*W1 zKQ-U~`(x9YOWJVoZm?HyxJd=&%m>c!iHb824IfLKs{uOCR6Yp24)VPxaJ z41BTMeRV%kei}!(_V10m?Rw*{=WCsa(RpX5q*Zh+3|~^MXo+I@b%=VR>5PdjW4Ay$ zMsRXC51&wvA+a-UXs@}G@%sss&`nN2$UegNw=)lX9-f9md7j>2v-ap{UhI8&1b_>OYFs+jg^T`e|7L5dG_;7JMkP zutRQ=udf?7O->qSDc<`4yYw+7<4x^{@$ocr!KZy>(*U)VWI z#L~vPD8CVhnbM!n#RB)i$C4aVOU*{w)5-iK>blFJI4p)0;*QF;xjsVqYz4c&J88C!Ao!u?9PrNHCRkcF zsrZ5{cyL9yn|@7rc*+aCxfo$b+V41C&yeXut9rD7U z41m&xFcYkE7B-o^Ky2X+Zh4wBwdARkL<__gxYNlrfS~8)4uZ$m1zVs7>sPx6&jAZr z9DNGw&S8_Az}3{(&0thV&!;q^q->G|r^X0rY+<`<`+-*#zZk&}=jQThN8=^fGXoE( zgdR5{x52&^HuS=x*82`FXsZxJb_9r@JQFOWi4XMAo zL0H0`CqzPxCqFTq!AE@1JL<-&zFk;Dbc~o7MiCbHa9+Ed^{219e9d*(z+nXvQ?l3U z&S}y?L>D;`xk(kC2d!kF(M(qql9ZsNuW*l88UGcuayhZ7=*$X?v@nIAUZsx;xiXtg zVSrjZBnGovLOf^ac^X8pVDrf7Nd0Z|GPn-oI*Yc;>GiYaw|)sBo&GR0kBo_8Z6N`? zKhB(vSGF{Zwjxm;J%WKz1PqZ2Hb(5_m?c_NTQ^1Rv%#gbOdU1KO&T8Lg1nF}xCy;i z9LK6Qsi3e*HYiA300x?z|7tNgXqeT0Sheczi*xJ;P8$T29^V756so%%NTrqcJs=+b zlQ>lwuE$1fwx#L98R{*D*5r}(-WrNPgvrA8ra?*gS8yow9lsB0zu(Xb{d$1>zODdJ zcr`LTlJf3C_THXOvY}S@RXMYNjfaL@+&rHq^*FSRwljWI^ubU2{wMmgytjVa4m+i+ zef-qRR zX3U=+GISzoQ3hSIP_yjbg-sg*$8;Fizua|YlJf6kfOY^Z zcTx8wnqDiBeJi@IbnuZS306vMWDGAk4L6+iWG~$rmiP?7pF3>4BG#bdu07AI17azn z8LtzzmYqo%O_OH>P7a$PzCL-~go&7dFD-CDS;uI$=1tfoVlu3i^aZQ2E~mGM^u75? z_7mf02nZ68oTQjW6#@W2wQo;wU7NM9XNX#FY)1jwI)TL-HpV_Dyu@Sb?EDk3^#Iu` zUq=zfzxPM`w|x%ALw(UAv@&a zs}&d<0w*Btris+{0~C|V0k;Wk02RP?_4MEXo7#WT*^5>WSpPqVOFv!&;x1v!@dbZdrn+y5^iLYi<=ynTkddCta1Cp9EdKpwQsuVICFM+j)1* zUJ_mp+0|n$HGqNKb7KB?=cYHMn}+cr0j06KcCMYw&nSoE?Xa-yKal zKO8MmFr?v;mx3E{G96igx>+LHpm48p>7Bk)xT-=#pvaI{GXM8OBe!P1d#gYAbHnHx z{QtfW_bgqP05f(|H=>y`{Zaq(sjI)xpe6xWofO8t2X2fW+)M>67cZuiy3(IKPuma-T!@WEe&p z4h^0Z{cv>r!|?Vr1PO$`jl{nz1~oigHxViz<aN#S*O^1m{L##)=I}I zQiVgTMlCq$O)+Z4&ZPhdS1%{-8iuw5b2cNO86$s?;j|5IXPHQM4`eMX??%A2`N?%l zi+x+}chgV`A|tkEg>AxoTLWc{G_LLbC~NE=mEpQYq8C>hzrS<9U1#V}Y(2tW?Vf}z z-wq(IUEev!=B6_dqJZ^&cpjtY8GoRru%GG*PssY1;C=YN{&IS}W#TZle#(l_BK%bu z2+Ud8ar}jMEX2Z;}?LSvBoV?Yt{ejl?e3daukQNSL zM9LEjjQIQ~Q=A@OH{=GH;_ZV7yl@gvl`rTNjw@e|Y)_??bdXuUTZ$T(5iuK}@frAo znc>BQ-v-owKp~FN7Ml%H`p8mBCy)&iaS|hdpRRu7u{(f-!=u@Sh2L%jJXuKJUeol- zN7#U3K0F}8Fw8Ru+QZ(KuCCc685yoi=qDx%Ya65={J)B>IgDx?5>RO1IKAdUObgG)Ti7J^H)efA9Uh@7?{Kcy`|BJm-ZRBKJLF z$}IfH^h0PWxp5;2z-U%tqyq|Z>C~d4v@VeT1c^zR$)(M%@Nqv9Ac#FN?tr|ySq`3> zo)P)cNDY>1i_m4JN7utP5aWr0V%S`$^l4v{)=ufD+WKSDLeYFuk zOJWs7Z2ZeLQ9S^Z-QO$?05jS}zAgy~$^m_)K(GAKo_~3QxEL#jMc#Y`5F{{^^1DW< zXDvcwG8qgaB2P}9XgSG?TZ3|WxHD%$1KV3DxAd@4&4E=r{t~C-f7j|@Xykq)-ts_4 zI~(RMvj8(N3_kVbosyb!@W^YQ32XQOb;hdZeNgm{*dm5llixrt6QEoOoZnG(f>wH5 zK-Bs-!;L^?oUp3hi?~q0MsAy8x&lC53x^Y#-xH8UOd+OojfyoE4i9roe5Jpe+4fc8 z-sv`E0jdN!O?yJefK`_eLscBay^r_6Fk)X4pziQz9=gnK(%$?9X1KJ`{Cu~Xts=<=XI@~Jkl^by+Db$}HfF8Md?K8Hy|`bfye zDYoMG`5YDyUkdt^b=sQ}Z@f1<2X}ATI)96AiNs~io2W9RU6a^e_g05;Wk8+VG1wucM zw56CFN&!tmXcT;Q{esspyy6BeTErE4pZMB*iltHa0*14Nk^0|yOb#c z^jKT`jA6$7Y4&*VP~hbK8#)Pq>UiW9mRDIzU;^k>1}hdJDA>oSS&~%}B_$}yw>AN8 z<70T|cbY!-GSi4zgypQsK%oPN(^L?0`oJY(S5_2HOKKMd;p{5=`}S_EguLFi-&iJ<8(hA zG2o}_d11*cjC9v8A_yS`TCzDr-$(mrj(6L-nVoAF@cZr#l2&TEw@``CX7}C(Le<$I z0Lf?3d&$@IhE8hbMyWTtJ+A@#KERvH+`?TiV1c5LM7?P}=r!)=o~_qwYxp;(hnf!_ zMLVU39i4>9r)1r-{qxg1c&CDaI+mkf$O`$f-A15Pa+(=6Ay98-W9Q79;w-H%C7&9B zt>20t&GwiutrD1S-k!Yi^a02!jwn7*Vdl{$SEjSQH8aX@>v?bS)V9%Ss6lye1?2{x z2a-dJVs;IE{`$b!i_AElM8D@9pZg+{CBTgk63;0!0s_Un=dPi#d0sxwE(=X0@~$ zXbo{*x^MK_dZ<04Pj@c?i3oYKx`s(0S97NBvS(?;e-HJ-YL7^a)fdcYwh(y1 z?|}TMFwPLtt+2e8`ClX%PJ7~>GS4ouxsWZJbze*B$8+p_4m0cOuGwWAs4r~=<*Zgp zX{q-9?#@P3DqNg-%?K}=SYdb5d@{F;>qpzs_L*victj#OPk9;kPhjTxRpL0G9bv`N z4Mv-Q;b*+1Zv)7MLzq%g5!xMg#^kESq3krLlGw==U(NYG8E*7iri53x0e|Gf=Y#i; z3+aYd-OY*fA-P^!V1@Ch-r@1vusP$8cJB5tFL8diC*;5O@4Fda1V_${J9YKzAbOle z8q#L2eVtVZ3mG6un9Fk5xu0*^vpK1dVv}K3nfDZM-i>-nY10YX=KIZ3ccOMF$rtTo z5yp01xq^9fH2ga;5&0qelxvgU5ABN4azNnD-mg_+7$Q*A1dm^a~Wc-!_zI@ch$)oxS;Z~HNJbt}OWDZ?F8o36RRMDZ3_w;+NxHv0E#S5Rztr`2sceEQ- z`LW#2`{vNQaZx0fD2C}wDXBMc=|8{o#=n=-nQlAJ=v|E_sUb-dHm3)7SERMHB$KQX zpODdt6KuT~-y%Bjf|PWW8EAF;SRTf=31h>@DjD0S;^Es1LcXPKZ~N-v88AzP&ka)9 zJz3*4jVsoZH#47(V8ypfR?FIWc!Xw<<@v7RP!c+sR~|T9a?s)X2R+9p<>=kzQ z*OA?4>cO{{G)t=;&n#1GuA0wI85xESiZ-5WWh~qN9>IMV4}ws9_;-islx8aXa$K#& z%;nY~Hm4#*BJNPUXWnAFWM*%dk^VJi@a=0>ET|+6#R%4XF@Z?KU&2UV4GTM`o*LvD zA)G8=wER@@JZBA^1GpL2W-Vp?jyaHAL6F`6x3sRD?44OvZ8!)5-Ylx-tX%tps1J$7 z{qbuaNv04S=_LJ2{S9Qo_#2N zB>TtGPmJ;2gIq;t@y)M@2`HJWO{B05(^*>VhYq9W=t2|q+K=j)_x#dB-6C>zM`Dpy(cDvyfu;R^8t$$t zA`kl3ckz~j7m+R@A^KGVsi4S8yNjd%?$m9a6>FED(>6pvxMd1}3S6c?!`(A#e6TnD zgEds2?ErK0z>7ELL)VV6NV6ytK1-0k478v)Ur!!LGtJicrKLaL8dtDb$}lBvTulhF z)byRM*4-pP7P_#IE}AP{l#{QA9&~$Q%P5ZFNJL5cvqIgrsiV*tkk8Q zuR8lE3Ohy{#?N!35!3resDpnyfkorX;M8!nF}4lZOL@^EitVq7Ze>;5t*xno$|m57 ziprY0x^5tsJ+T%S^SWTpO=ZjuY6GT@JJlh#zWlU3J-G3czHWRW;CT18f02AT*1q_K4(sYx&uWjf1;ub^~J~hh>>n^D}63M*MvXjYDjl;!4iKNB;ty%8f1(h&}X~ zA{>DjF%Tmzs)lBET&A8N1ovrPRgoxO6biWdJi5)#tj>8~{=Fi=5qr!rj$5C>LYY>xHFh8YsfM5MAp$ifU*!vn++8YXmIovq-L zoE(_}uGx2hY2CyR{KDL+Z%KnAcYvs<{q#t|1xl?X8N^iQ>o0LF%xHBgsTq?&InO~? zGQBf@G#{|uz`(Xw6%dogAY05dE_}l0Zh-BnJNOp}l;s4oHLI7A2Uh|)5oybW>7=vk z9j8BVX4^X`wbL^OR?@@E>s*{-xF%lgw~XFX=6!5U&s|tp4^ecPgd;&k5I?=Ew0Iv zBg6H&$0KF(2bfoJpK`K+9MzE}iN6|_#PokQrtn9LqO%_s;3vD%RK zB76RHFyCsLcIoD7WZ|@sE?_A|ra#h*h~wZCaLI`1mv<4NjD=9xU>(bVrJm)^IsbXkY#SLmH`u>b>eVcO-T&z0>& z(I(k0qw~050T@EiYy~<{XEvL|lapofQ5FI8b?pUO(!5zxJ*u7h5c4_myc$0~$~0m( zR(_W8zis~fckZV9GhxMz+`_2kA%Px!>=&mn|_ zZLjoaO1@hb{UH}-`DxDC0fGPs$k4eeO>r~*m+Pz&~ z{f;}A6e%L2l6<|s&iiAb@rn+3dnVb$>`vPKL9fC!seEG7-?#~1=L-cXU`3_;-aAIJ z`7?Qmh6a3={xfECd$`J)t;+TLFQv8nrqS&d*?MnoB`Yf@-|Of>rk3A>28zsn($a{0 z15YOR;I-{2=(h||_HlMxJA;dft?JO5rEEoNYK{7T0Y}?cL5S1b36WGfcO55)r?Dre@2hm)O19OsF71QYN%XTrL6cNFL~)Om|+jrZ%WC`P?-wO>&P9EXfrEFF)w zb?idc8$e0G_HU}kZNF!8L>(A5#)lZOdv6X8WN<&6up_dh8)tLT;4CH+aPf5Os;_@y zoc#p=X{&1+QHIRLEw4v2o^|vuBfc_VVJ-eGsK=h;0<=XPXZ|)56#~prPS$^EYQ7t=3B(2(R|Bcs1!C72X?61d)0(p{ zufs;jx4rKQ!qKRFdl%_}i{jdoX_1OeofJx)JWU%K@jwnQo``pnrE7@4GwOIV{wkm9b|)A$&}w zJ1DPNSdp7_c)v8@RqlmN%HKeu5#S8(%k9m6y#UElmlzgcoSoH?Q^PHQ-dFkE|J5_jI3_kL@DcPzQcK(}X13#~ z%z^M(aIzz_eoCm*kAWor%ZxIoyp!~cOb4eK2M4)0*A`<8X@Z!(r`5pNod;U=c+7;; zk>ea&_Q#@AahA#1N#T|4Se>SfmYpM4H7X5R-9x_(DaOlv)(pd73l%b9=SG zoHj2#+{tMc^k0;lwh4L^U@gm9TOx+1d2`~sqI`AT?#l)7npAspg??(;VyrftKz-kG z=bod{9;_u!)-7j$qgqo5bX@H!Dl7q0=c9ehsZ#bIZ&`kRNz@P-Wj0>2p{kEpUkAN&sZf5c+`<&Pp?M;ltQF)U$(NAQE{LZBbdI0HjvqgX#N z_gut(hr7sov(oYSeLQG(CaZB;;!v}K3y6&f>tXQ5)}b<>WiQDRgeF4w78RU$l$L85 zmzeJl!IeZ4`h#>oknwi}*1J8>B%SYM`@AMt@Ejx?6{IwkL*W@vm-_oTy|*<%BD2L9 z!Fy_v23`g!$6sQzYleH?SA+yszsbE;NJFU#mX$O;@rJG26p@#%mI5?YIb*OZu6^YXl@+9_SpST9 zQ>V$@U4f4PcICSo$iG(iJ^^fOlrrbR;3am5Gllgb`ICAEQXsI`aRq9(!EqU$8 zs38K*7$5RwBXK~ht+a9y00#VfYScs?wj4&UZ+F1??6-;+{iD;HVS7^PgfDbPOd?Id zjOJ8?= 1.0: # boundary side + return phi(s) * (-psi(s) * vel_n + 0.5 * psi(s) * wp.abs(vel_n)) + + # interior side return fem.jump(phi, s) * (-fem.average(psi, s) * vel_n + 0.5 * fem.jump(psi, s) * wp.abs(vel_n)) @@ -63,7 +66,7 @@ def sip_diffusion_form( class Example: - def __init__(self, quiet=False, degree=2, resolution=50, mesh="grid", viscosity=0.001, ang_vel=1.0): + def __init__(self, quiet=False, degree=2, resolution=50, mesh="grid", viscosity=0.0001, ang_vel=1.0): self._quiet = quiet res = resolution @@ -108,37 +111,30 @@ def __init__(self, quiet=False, degree=2, resolution=50, mesh="grid", viscosity= side_test = fem.make_test(space=scalar_space, domain=sides) side_trial = fem.make_trial(space=scalar_space, domain=sides) - bsr_axpy( - fem.integrate( - upwind_transport_form, - fields={"phi": side_trial, "psi": side_test}, - values={"ang_vel": ang_vel}, - ), - y=matrix_transport, + matrix_transport += fem.integrate( + upwind_transport_form, + fields={"phi": side_trial, "psi": side_test}, + values={"ang_vel": ang_vel}, ) matrix_diffusion = fem.integrate( diffusion_form, fields={"u": trial, "v": self._test}, ) - bsr_axpy( - fem.integrate( - sip_diffusion_form, - fields={"phi": side_trial, "psi": side_test}, - ), - y=matrix_diffusion, + + matrix_diffusion += fem.integrate( + sip_diffusion_form, + fields={"phi": side_trial, "psi": side_test}, ) - self._matrix = matrix_inertia - bsr_axpy(x=matrix_transport, y=self._matrix) - bsr_axpy(x=matrix_diffusion, y=self._matrix, alpha=viscosity) + self._matrix = matrix_inertia + matrix_transport + viscosity * matrix_diffusion # Initial condition self._phi_field = scalar_space.make_field() fem.interpolate(initial_condition, dest=self._phi_field) self.renderer = fem_example_utils.Plot() - self.renderer.add_surface("phi", self._phi_field) + self.renderer.add_field("phi", self._phi_field) def step(self): self.current_frame += 1 @@ -156,7 +152,7 @@ def step(self): def render(self): self.renderer.begin_frame(time=self.current_frame * self.sim_dt) - self.renderer.add_surface("phi", self._phi_field) + self.renderer.add_field("phi", self._phi_field) self.renderer.end_frame() diff --git a/warp/examples/fem/example_deformed_geometry.py b/warp/examples/fem/example_deformed_geometry.py index 4cfd80662..8c4e99c45 100644 --- a/warp/examples/fem/example_deformed_geometry.py +++ b/warp/examples/fem/example_deformed_geometry.py @@ -35,7 +35,7 @@ def deformation_field_expr( r = x[1] + 0.5 t = 0.5 * 3.1416 * x[0] - return r * wp.vec2(wp.sin(t), wp.cos(t)) - x + return r * wp.vec2(wp.sin(t), wp.cos(t)) @fem.integrand @@ -82,7 +82,7 @@ def __init__( deformation_field = deformation_space.make_field() fem.interpolate(deformation_field_expr, dest=deformation_field) - self._geo = deformation_field.make_deformed_geometry() + self._geo = deformation_field.make_deformed_geometry(relative=False) # Scalar function space on deformed geometry element_basis = fem.ElementBasis.SERENDIPITY if serendipity else None @@ -123,7 +123,7 @@ def step(self): self._scalar_field.dof_values = x def render(self): - self.renderer.add_surface("solution", self._scalar_field) + self.renderer.add_field("solution", self._scalar_field) if __name__ == "__main__": diff --git a/warp/examples/fem/example_diffusion.py b/warp/examples/fem/example_diffusion.py index adcedc296..aaaa08f9b 100644 --- a/warp/examples/fem/example_diffusion.py +++ b/warp/examples/fem/example_diffusion.py @@ -20,7 +20,6 @@ import warp.examples.fem.utils as fem_example_utils import warp.fem as fem from warp.fem.utils import array_axpy -from warp.sparse import bsr_axpy @fem.integrand @@ -130,7 +129,7 @@ def step(self): else: # Weak BC: add together diffusion and boundary condition matrices boundary_strength = 1.0 / self._boundary_compliance - bsr_axpy(x=bd_matrix, y=matrix, alpha=boundary_strength, beta=1) + matrix += bd_matrix * boundary_strength array_axpy(x=bd_rhs, y=rhs, alpha=boundary_strength, beta=1) # Solve linear system using Conjugate Gradient @@ -141,7 +140,7 @@ def step(self): self._scalar_field.dof_values = x def render(self): - self.renderer.add_surface("solution", self._scalar_field) + self.renderer.add_field("solution", self._scalar_field) if __name__ == "__main__": diff --git a/warp/examples/fem/example_diffusion_3d.py b/warp/examples/fem/example_diffusion_3d.py index 110f65521..4174bd0ed 100644 --- a/warp/examples/fem/example_diffusion_3d.py +++ b/warp/examples/fem/example_diffusion_3d.py @@ -31,7 +31,8 @@ def vert_boundary_projector_form( v: fem.Field, ): # Non-zero mass on vertical sides only - w = 1.0 - wp.abs(fem.normal(domain, s)[1]) + nor = fem.normal(domain, s) + w = wp.abs(nor[0]) return w * u(s) * v(s) @@ -51,7 +52,7 @@ def __init__( self._viscosity = viscosity self._boundary_compliance = boundary_compliance - res = wp.vec3i(resolution, resolution // 2, resolution * 2) + res = wp.vec3i(resolution, max(1, resolution // 2), resolution * 2) if mesh == "tet": pos, tet_vtx_indices = fem_example_utils.gen_tetmesh( @@ -122,7 +123,7 @@ def step(self): self._scalar_field.dof_values = x def render(self): - self.renderer.add_volume("solution", self._scalar_field) + self.renderer.add_field("solution", self._scalar_field) if __name__ == "__main__": diff --git a/warp/examples/fem/example_diffusion_mgpu.py b/warp/examples/fem/example_diffusion_mgpu.py index 2e4b35fe7..281f75eff 100644 --- a/warp/examples/fem/example_diffusion_mgpu.py +++ b/warp/examples/fem/example_diffusion_mgpu.py @@ -159,7 +159,7 @@ def step(self): array_cast(in_array=global_res, out_array=self._scalar_field.dof_values) def render(self): - self.renderer.add_surface("solution", self._scalar_field) + self.renderer.add_field("solution", self._scalar_field) def _assemble_local_system(self, geo_partition: fem.GeometryPartition): scalar_space = self._scalar_space diff --git a/warp/examples/fem/example_magnetostatics.py b/warp/examples/fem/example_magnetostatics.py new file mode 100644 index 000000000..f197ba669 --- /dev/null +++ b/warp/examples/fem/example_magnetostatics.py @@ -0,0 +1,190 @@ +# Copyright (c) 2024 NVIDIA CORPORATION. All rights reserved. +# NVIDIA CORPORATION and its licensors retain all intellectual property +# and proprietary rights in and to this software, related documentation +# and any modifications thereto. Any use, reproduction, disclosure or +# distribution of this software and related documentation without an express +# license agreement from NVIDIA CORPORATION is strictly prohibited. + +########################################################################### +# Example Magnetostatics +# +# This example demonstrates solving an in-plane magnetostatics problem +# using a curl-curl formulation +# +# 1/mu Curl B + j = 0 +# Div. B = 0 +# +# solved over field A such that B = Curl A, A = (0, 0, a_z), +# and a_z = 0 on the domain boundary +# +# This example also illustrates using an ImplictField to warp a square mesh +# to a circular domain +########################################################################### + +import numpy as np + +import warp as wp +import warp.examples.fem.utils as fem_example_utils +import warp.fem as fem + +# Vacuum and copper magnetic permeabilities +MU_0 = wp.constant(np.pi * 4.0e-7) +MU_C = wp.constant(1.25e-6) + + +@wp.func +def square_to_disk(x: wp.vec2): + # mapping from unit square to unit disk + return wp.normalize(x) * wp.max(wp.abs(x)) + + +@wp.func +def square_to_disk_grad(x: wp.vec2): + # gradient of mapping from unit square to unit disk + if x == wp.vec2(0.0): + return wp.mat22(0.0) + + d = wp.normalize(x) + d_grad = (wp.identity(n=2, dtype=float) - wp.outer(d, d)) / wp.length(x) + + ax = wp.abs(x) + xinf = wp.max(ax) + + xinf_grad = wp.select(ax[0] > ax[1], wp.vec2(0.0, wp.sign(x[1])), wp.vec(wp.sign(x[0]), 0.0)) + + return d_grad * xinf + wp.outer(d, xinf_grad) + + +@wp.func +def permeability_field(pos: wp.vec2, coil_height: float, coil_internal_radius: float, coil_external_radius: float): + # space-varying permeability + + x = wp.abs(pos[0]) + y = wp.abs(pos[1]) + return wp.select( + y < coil_height and x > coil_internal_radius and x < coil_external_radius, + MU_0, + MU_C, + ) + + +@wp.func +def current_field( + pos: wp.vec2, coil_height: float, coil_internal_radius: float, coil_external_radius: float, current: float +): + # space-varying current direction along z axis (0, +1 or -1) + x = wp.abs(pos[0]) + y = wp.abs(pos[1]) + return ( + wp.select( + y < coil_height and x > coil_internal_radius and x < coil_external_radius, + 0.0, + wp.sign(pos[0]), + ) + * current + ) + + +@fem.integrand +def curl_z(u: fem.Field, s: fem.Sample): + # projection of curl((0, 0, u)) over z axis + du = fem.grad(u, s) + return wp.vec2(du[1], -du[0]) + + +@fem.integrand +def curl_curl_form(s: fem.Sample, domain: fem.Domain, u: fem.Field, v: fem.Field, mu: fem.Field): + return wp.dot(curl_z(u, s), curl_z(v, s)) / mu(s) + + +@fem.integrand +def mass_form(s: fem.Sample, domain: fem.Domain, v: fem.Field, u: fem.Field): + return u(s) * v(s) + + +class Example: + def __init__(self, quiet=False, degree=2, resolution=32, domain_radius=2.0, current=1.0e6): + # We mesh the unit disk by first meshing the unit square, then building a deformed geometry + # from an implicit mapping field + square_geo = fem.Grid2D( + bounds_lo=wp.vec2(-domain_radius, -domain_radius), + bounds_hi=wp.vec2(domain_radius, domain_radius), + res=wp.vec2i(resolution, resolution), + ) + + def_field = fem.ImplicitField(domain=fem.Cells(square_geo), func=square_to_disk, grad_func=square_to_disk_grad) + disk_geo = def_field.make_deformed_geometry(relative=False) + + coil_config = {"coil_height": 1.0, "coil_internal_radius": 0.1, "coil_external_radius": 0.3} + + domain = fem.Cells(disk_geo) + self._permeability_field = fem.ImplicitField(domain, func=permeability_field, values=coil_config) + self._current_field = fem.ImplicitField(domain, func=current_field, values=dict(current=current, **coil_config)) + + z_space = fem.make_polynomial_space(disk_geo, degree=degree, element_basis=fem.ElementBasis.LAGRANGE) + xy_space = fem.make_polynomial_space( + disk_geo, degree=degree, element_basis=fem.ElementBasis.LAGRANGE, dtype=wp.vec2 + ) + + self.A_field = z_space.make_field() + self.B_field = xy_space.make_field() + + self.renderer = fem_example_utils.Plot() + + def step(self): + z_space = self.A_field.space + disk_geo = z_space.geometry + + u = fem.make_trial(space=z_space) + v = fem.make_test(space=z_space) + lhs = fem.integrate(curl_curl_form, fields={"u": u, "v": v, "mu": self._permeability_field}) + rhs = fem.integrate(mass_form, fields={"v": v, "u": self._current_field}) + + # Dirichlet BC + boundary = fem.BoundarySides(disk_geo) + u_bd = fem.make_trial(space=z_space, domain=boundary) + v_bd = fem.make_test(space=z_space, domain=boundary) + dirichlet_bd_proj = fem.integrate(mass_form, fields={"u": u_bd, "v": v_bd}, nodal=True) + fem.project_linear_system(lhs, rhs, dirichlet_bd_proj) + + x = wp.zeros_like(rhs) + fem_example_utils.bsr_cg(lhs, b=rhs, x=x, tol=1.0e-8, quiet=False) + + # make sure result is exactly zero outisde of circle + wp.sparse.bsr_mv(dirichlet_bd_proj, x=x, y=x, alpha=-1.0, beta=1.0) + wp.utils.array_cast(in_array=x, out_array=self.A_field.dof_values) + + # compute B as curl(A) + fem.interpolate(curl_z, dest=self.B_field, fields={"u": self.A_field}) + + def render(self): + self.renderer.add_field("A", self.A_field) + self.renderer.add_field("B", self.B_field) + + +if __name__ == "__main__": + import argparse + + wp.set_module_options({"enable_backward": False}) + + parser = argparse.ArgumentParser(formatter_class=argparse.ArgumentDefaultsHelpFormatter) + parser.add_argument("--device", type=str, default=None, help="Override the default Warp device.") + parser.add_argument("--resolution", type=int, default=32, help="Grid resolution.") + parser.add_argument("--degree", type=int, default=2, help="Polynomial degree of shape functions.") + parser.add_argument("--radius", type=float, default=2.0, help="Radius of simulation domain.") + parser.add_argument( + "--headless", + action="store_true", + help="Run in headless mode, suppressing the opening of any graphical windows.", + ) + parser.add_argument("--quiet", action="store_true", help="Suppresses the printing out of iteration residuals.") + + args = parser.parse_known_args()[0] + + with wp.ScopedDevice(args.device): + example = Example(quiet=args.quiet, degree=args.degree, resolution=args.resolution, domain_radius=args.radius) + example.step() + example.render() + + if not args.headless: + example.renderer.plot({"A": {"contours": {"levels": 30}}, "B": {"streamlines": {"density": 1.0}}}) diff --git a/warp/examples/fem/example_mixed_elasticity.py b/warp/examples/fem/example_mixed_elasticity.py index 65ec1af72..de361bd90 100644 --- a/warp/examples/fem/example_mixed_elasticity.py +++ b/warp/examples/fem/example_mixed_elasticity.py @@ -25,7 +25,6 @@ import warp as wp import warp.examples.fem.utils as fem_example_utils import warp.fem as fem -from warp.sparse import bsr_mm, bsr_mv, bsr_transposed @fem.integrand @@ -193,9 +192,7 @@ def step(self): tau_test = fem.make_test(space=self._tau_space, domain=domain) tau_trial = fem.make_trial(space=self._tau_space, domain=domain) - gradient_matrix = bsr_transposed( - fem.integrate(displacement_gradient_form, fields={"u": u_trial, "tau": tau_test}) - ) + gradient_matrix = fem.integrate(displacement_gradient_form, fields={"u": u_trial, "tau": tau_test}).transpose() # Compute inverse of the (block-diagonal) tau mass matrix tau_inv_mass_matrix = fem.integrate(tensor_mass_form, fields={"sig": tau_trial, "tau": tau_test}, nodal=True) @@ -217,10 +214,10 @@ def step(self): ) # Assemble system matrix - u_matrix = bsr_mm(gradient_matrix, bsr_mm(tau_inv_mass_matrix, stress_matrix)) + u_matrix = gradient_matrix @ tau_inv_mass_matrix @ stress_matrix # Enforce boundary conditions (apply displacement only at first iteration) - u_rhs = bsr_mv(gradient_matrix, bsr_mv(tau_inv_mass_matrix, stress_rhs, alpha=-1.0)) + u_rhs = -gradient_matrix @ (tau_inv_mass_matrix @ stress_rhs) fem.project_linear_system(u_matrix, u_rhs, u_bd_matrix, u_bd_rhs if newton_iteration == 0 else None) x = wp.zeros_like(u_rhs) @@ -232,7 +229,7 @@ def step(self): fem.utils.array_axpy(x=delta_u, y=self._u_field.dof_values) def render(self): - self.renderer.add_surface_vector("solution", self._u_field) + self.renderer.add_field("solution", self._u_field) if __name__ == "__main__": @@ -273,4 +270,4 @@ def render(self): example.render() if not args.headless: - example.renderer.plot(displacement="solution") + example.renderer.plot(options={"solution": {"displacement": {}}}) diff --git a/warp/examples/fem/example_navier_stokes.py b/warp/examples/fem/example_navier_stokes.py index b4759a54b..730583542 100644 --- a/warp/examples/fem/example_navier_stokes.py +++ b/warp/examples/fem/example_navier_stokes.py @@ -21,16 +21,15 @@ import warp.examples.fem.utils as fem_example_utils import warp.fem as fem from warp.fem.utils import array_axpy -from warp.sparse import bsr_copy, bsr_mm, bsr_mv -@fem.integrand -def u_boundary_value(s: fem.Sample, domain: fem.Domain, v: fem.Field, top_vel: float): +@wp.func +def u_boundary_value(x: wp.vec2, box_height: float, top_velocity: float): # Horizontal velocity on top of domain, zero elsewhere - if domain(s)[1] == 1.0: - return wp.dot(wp.vec2f(top_vel, 0.0), v(s)) + if x[1] >= box_height: + return wp.vec2(top_velocity, 0.0) - return wp.dot(wp.vec2f(0.0, 0.0), v(s)) + return wp.vec2(0.0, 0.0) @fem.integrand @@ -124,10 +123,14 @@ def __init__(self, quiet=False, degree=2, resolution=25, Re=1000.0, top_velocity u_bd_test = fem.make_test(space=u_space, domain=boundary) u_bd_trial = fem.make_trial(space=u_space, domain=boundary) u_bd_projector = fem.integrate(mass_form, fields={"u": u_bd_trial, "v": u_bd_test}, nodal=True) + + # Define an implicit field for our boundary condition value and integrate + u_bd_field = fem.ImplicitField( + domain=boundary, func=u_boundary_value, values={"top_velocity": top_velocity, "box_height": 1.0} + ) u_bd_value = fem.integrate( - u_boundary_value, - fields={"v": u_bd_test}, - values={"top_vel": top_velocity}, + mass_form, + fields={"u": u_bd_field, "v": u_bd_test}, nodal=True, output_dtype=wp.vec2d, ) @@ -137,12 +140,8 @@ def __init__(self, quiet=False, degree=2, resolution=25, Re=1000.0, top_velocity u_bd_rhs = wp.zeros_like(u_bd_value) fem.project_linear_system(u_matrix, u_bd_rhs, u_bd_projector, u_bd_value, normalize_projector=False) - # div_bd_rhs = div_matrix * u_bd_rhs - div_bd_rhs = wp.zeros(shape=(div_matrix.nrow,), dtype=div_matrix.scalar_type) - bsr_mv(div_matrix, u_bd_value, y=div_bd_rhs, alpha=-1.0) - - # div_matrix = div_matrix - div_matrix * bd_projector - bsr_mm(x=bsr_copy(div_matrix), y=u_bd_projector, z=div_matrix, alpha=-1.0, beta=1.0) + div_bd_rhs = -div_matrix @ u_bd_value + div_matrix -= div_matrix @ u_bd_projector # Assemble saddle system self._saddle_system = fem_example_utils.SaddleSystem(u_matrix, div_matrix) @@ -158,7 +157,7 @@ def __init__(self, quiet=False, degree=2, resolution=25, Re=1000.0, top_velocity self._p_field = p_space.make_field() self.renderer = fem_example_utils.Plot() - self.renderer.add_surface_vector("velocity", self._u_field) + self.renderer.add_field("velocity", self._u_field) def step(self): self.current_frame += 1 @@ -172,7 +171,7 @@ def step(self): # Apply boundary conditions # u_rhs = (I - P) * u_rhs + u_bd_rhs - bsr_mv(self._u_bd_projector, x=u_rhs, y=u_rhs, alpha=-1.0, beta=1.0) + wp.sparse.bsr_mv(self._u_bd_projector, x=u_rhs, y=u_rhs, alpha=-1.0, beta=1.0) array_axpy(x=self._u_bd_rhs, y=u_rhs, alpha=1.0, beta=1.0) p_rhs = self._div_bd_rhs @@ -197,7 +196,7 @@ def step(self): def render(self): self.renderer.begin_frame(time=self.current_frame * self.sim_dt) - self.renderer.add_surface_vector("velocity", self._u_field) + self.renderer.add_field("velocity", self._u_field) self.renderer.end_frame() @@ -243,7 +242,12 @@ def render(self): example.step() example.render() - example.renderer.add_surface_vector("velocity_final", example._u_field) + example.renderer.add_field("velocity_final", example._u_field) if not args.headless: - example.renderer.plot(streamlines={"velocity_final"}) + example.renderer.plot( + options={ + "velocity": {"arrows": {"glyph_scale": 0.25}}, + "velocity_final": {"streamlines": {"density": 2.0}}, + } + ) diff --git a/warp/examples/fem/example_nonconforming_contact.py b/warp/examples/fem/example_nonconforming_contact.py new file mode 100644 index 000000000..c30f482c1 --- /dev/null +++ b/warp/examples/fem/example_nonconforming_contact.py @@ -0,0 +1,290 @@ +# Copyright (c) 2024 NVIDIA CORPORATION. All rights reserved. +# NVIDIA CORPORATION and its licensors retain all intellectual property +# and proprietary rights in and to this software, related documentation +# and any modifications thereto. Any use, reproduction, disclosure or +# distribution of this software and related documentation without an express +# license agreement from NVIDIA CORPORATION is strictly prohibited. + +########################################################################### +# Example Nonconforming Contact +# +# This example demonstrates using nonconforming fields (warp.fem.NonconformingField) +# to solve a no-slip contact problem between two elastic bodies discretized separately. +# +# Div[ E: D(u) ] = g over each body +# u_top = u_bottom along the contact surface (top body bottom boundary) +# u_bottom = 0 along the bottom boundary of the bottom body +# +# with E the rank-4 elasticity tensor +# +# Below we use a simple staggered scheme for solving bodies iteratively, +# but more robust methods could be considered (e.g. Augmented Lagrangian) +########################################################################### + +import numpy as np + +import warp as wp +import warp.examples.fem.utils as fem_example_utils +import warp.fem as fem + + +@wp.func +def compute_stress(tau: wp.mat22, E: wp.mat33): + """Strain to stress computation (using Voigt notation to drop tensor order)""" + tau_sym = wp.vec3(tau[0, 0], tau[1, 1], tau[0, 1] + tau[1, 0]) + sig_sym = E * tau_sym + return wp.mat22(sig_sym[0], 0.5 * sig_sym[2], 0.5 * sig_sym[2], sig_sym[1]) + + +@fem.integrand +def stress_form(s: fem.Sample, u: fem.Field, tau: fem.Field, E: wp.mat33): + """Stress inside body: (E : D(u)) : tau""" + return wp.ddot(tau(s), compute_stress(fem.D(u, s), E)) + + +@fem.integrand +def boundary_stress_form( + s: fem.Sample, + domain: fem.Domain, + u: fem.Field, + tau: fem.Field, +): + """Stress on boundary: u' tau n""" + return wp.dot(tau(s) * fem.normal(domain, s), u(s)) + + +@fem.integrand +def symmetric_grad_form( + s: fem.Sample, + u: fem.Field, + tau: fem.Field, +): + """Symmetric part of gradient of displacement: D(u) : tau""" + return wp.ddot(tau(s), fem.D(u, s)) + + +@fem.integrand +def gravity_form( + s: fem.Sample, + v: fem.Field, + gravity: float, +): + return -gravity * v(s)[1] + + +@fem.integrand +def bottom_boundary_projector_form( + s: fem.Sample, + domain: fem.Domain, + u: fem.Field, + v: fem.Field, +): + # non zero on bottom boundary only + nor = fem.normal(domain, s) + return wp.dot(u(s), v(s)) * wp.max(0.0, -nor[1]) + + +@fem.integrand +def tensor_mass_form( + s: fem.Sample, + sig: fem.Field, + tau: fem.Field, +): + return wp.ddot(tau(s), sig(s)) + + +class Example: + def __init__( + self, + degree=2, + resolution=16, + young_modulus=1.0, + poisson_ratio=0.5, + nonconforming_stresses=False, + ): + self._geo1 = fem.Grid2D(bounds_hi=wp.vec2(1.0, 0.5), res=wp.vec2i(resolution)) + self._geo2 = fem.Grid2D(bounds_lo=(0.33, 0.5), bounds_hi=(0.67, 0.5 + 0.33), res=wp.vec2i(resolution)) + + # Strain-stress matrix + young = young_modulus + poisson = poisson_ratio + self._elasticity_mat = wp.mat33( + young + / (1.0 - poisson * poisson) + * np.array( + [ + [1.0, poisson, 0.0], + [poisson, 1.0, 0.0], + [0.0, 0.0, (2.0 * (1.0 + poisson)) * (1.0 - poisson * poisson)], + ] + ) + ) + + # Displacement spaces and fields -- S_k + self._u1_space = fem.make_polynomial_space( + self._geo1, degree=degree, dtype=wp.vec2, element_basis=fem.ElementBasis.SERENDIPITY + ) + self._u2_space = fem.make_polynomial_space( + self._geo2, degree=degree, dtype=wp.vec2, element_basis=fem.ElementBasis.SERENDIPITY + ) + self._u1_field = self._u1_space.make_field() + self._u2_field = self._u2_space.make_field() + + # Stress spaces and fields -- Q_{k-1}d + # Store stress degrees of freedom as symmetric tensors (3 dof) rather than full 2x2 matrices + self._tau1_space = fem.make_polynomial_space( + self._geo1, + degree=degree - 1, + discontinuous=True, + element_basis=fem.ElementBasis.LAGRANGE, + dof_mapper=fem.SymmetricTensorMapper(wp.mat22), + ) + self._tau2_space = fem.make_polynomial_space( + self._geo2, + degree=degree - 1, + discontinuous=True, + element_basis=fem.ElementBasis.LAGRANGE, + dof_mapper=fem.SymmetricTensorMapper(wp.mat22), + ) + + self._sig1_field = self._tau1_space.make_field() + self._sig2_field = self._tau2_space.make_field() + self._sig2_field_new = self._tau2_space.make_field() + + self.renderer = fem_example_utils.Plot() + + def step(self): + # Solve for the two bodies separately + # Body (top) is 25x more dense and 5x stiffer than top body + # Body 1 affects body 2 through bottom displacement dirichlet BC + # Body 2 affects body 1 through applied strain on top + self.solve_solid(self._u1_field, self._sig1_field, self._u2_field, self._sig2_field, gravity=1.0, stiffness=1.0) + self.solve_solid( + self._u2_field, self._sig2_field_new, self._u1_field, self._sig1_field, gravity=25.0, stiffness=5.0 + ) + + # Damped update of coupling stress (for stability) + alpha = 0.1 + fem.utils.array_axpy( + x=self._sig2_field_new.dof_values, y=self._sig2_field.dof_values, alpha=alpha, beta=1.0 - alpha + ) + + def solve_solid( + self, + u_field, + stress_field, + other_u_field, + other_stress_field, + gravity: float, + stiffness: float, + ): + u_space = u_field.space + stress_space = stress_field.space + geo = u_field.space.geometry + + domain = fem.Cells(geometry=geo) + boundary = fem.BoundarySides(geometry=geo) + + u_test = fem.make_test(space=u_space, domain=domain) + u_trial = fem.make_trial(space=u_space, domain=domain) + tau_test = fem.make_test(space=stress_space, domain=domain) + tau_trial = fem.make_trial(space=stress_space, domain=domain) + + u_bd_test = fem.make_test(space=u_space, domain=boundary) + u_bd_trial = fem.make_trial(space=u_space, domain=boundary) + + # Assemble stiffness matrix + # (Note: this is constant per body, this could be precomputed) + sym_grad_matrix = fem.integrate(symmetric_grad_form, fields={"u": u_trial, "tau": tau_test}) + + tau_inv_mass_matrix = fem.integrate(tensor_mass_form, fields={"sig": tau_trial, "tau": tau_test}, nodal=True) + fem_example_utils.invert_diagonal_bsr_matrix(tau_inv_mass_matrix) + + stress_matrix = tau_inv_mass_matrix @ fem.integrate( + stress_form, fields={"u": u_trial, "tau": tau_test}, values={"E": self._elasticity_mat * stiffness} + ) + stiffness_matrix = sym_grad_matrix.transpose() @ stress_matrix + + # Right-hand-side + u_rhs = fem.integrate(gravity_form, fields={"v": u_test}, values={"gravity": gravity}, output_dtype=wp.vec2d) + + # Add boundary stress from other solid field + other_stress_field = fem.field.field.NonconformingField(boundary, other_stress_field) + fem.utils.array_axpy( + y=u_rhs, + x=fem.integrate( + boundary_stress_form, fields={"u": u_bd_test, "tau": other_stress_field}, output_dtype=wp.vec2d + ), + ) + + # Enforce boundary conditions + u_bd_matrix = fem.integrate( + bottom_boundary_projector_form, fields={"u": u_bd_trial, "v": u_bd_test}, nodal=True + ) + + # read displacement from other body set create bottom boundary Dirichlet BC + other_u_field = fem.field.field.NonconformingField(boundary, other_u_field) + u_bd_rhs = fem.integrate( + bottom_boundary_projector_form, fields={"u": other_u_field, "v": u_bd_test}, nodal=True + ) + + fem.project_linear_system(stiffness_matrix, u_rhs, u_bd_matrix, u_bd_rhs) + + # solve + x = wp.zeros_like(u_rhs) + wp.utils.array_cast(in_array=u_field.dof_values, out_array=x) + fem_example_utils.bsr_cg(stiffness_matrix, b=u_rhs, x=x, tol=1.0e-6, quiet=True) + + # Extract result + stress = stress_matrix @ x + wp.utils.array_cast(in_array=x, out_array=u_field.dof_values) + wp.utils.array_cast(in_array=stress, out_array=stress_field.dof_values) + + def render(self): + self.renderer.add_field("u1", self._u1_field) + self.renderer.add_field("u2", self._u2_field) + self.renderer.add_field("sig1", self._sig1_field) + self.renderer.add_field("sig2", self._sig2_field) + + +if __name__ == "__main__": + import argparse + + wp.set_module_options({"enable_backward": False}) + + parser = argparse.ArgumentParser(formatter_class=argparse.ArgumentDefaultsHelpFormatter) + parser.add_argument("--device", type=str, default=None, help="Override the default Warp device.") + parser.add_argument("--resolution", type=int, default=32, help="Grid resolution.") + parser.add_argument("--degree", type=int, default=2, help="Polynomial degree of shape functions.") + parser.add_argument("--young_modulus", type=float, default=10.0) + parser.add_argument("--poisson_ratio", type=float, default=0.9) + parser.add_argument("--num_steps", type=int, default=50) + parser.add_argument( + "--headless", + action="store_true", + help="Run in headless mode, suppressing the opening of any graphical windows.", + ) + + args = parser.parse_known_args()[0] + + with wp.ScopedDevice(args.device): + example = Example( + degree=args.degree, + resolution=args.resolution, + young_modulus=args.young_modulus, + poisson_ratio=args.poisson_ratio, + ) + + for i in range(args.num_steps): + print("Step", i) + example.step() + example.render() + + if not args.headless: + example.renderer.plot( + { + "rows": 2, + "u1": {"displacement": {}}, + "u2": {"displacement": {}}, + }, + ) diff --git a/warp/examples/fem/example_stokes.py b/warp/examples/fem/example_stokes.py index 23495fc86..31971d5c4 100644 --- a/warp/examples/fem/example_stokes.py +++ b/warp/examples/fem/example_stokes.py @@ -19,15 +19,9 @@ import warp as wp import warp.examples.fem.utils as fem_example_utils import warp.fem as fem -import warp.sparse as sparse from warp.fem.utils import array_axpy -@fem.integrand -def constant_form(val: wp.vec2): - return val - - @fem.integrand def viscosity_form(s: fem.Sample, u: fem.Field, v: fem.Field, nu: float): return nu * wp.ddot(fem.D(u, s), fem.D(v, s)) @@ -103,11 +97,7 @@ def __init__( self._u_field = u_space.make_field() self._p_field = p_space.make_field() - # Interpolate initial condition on boundary (for example purposes) - self._bd_field = u_space.make_field() - f_boundary = fem.make_restriction(self._bd_field, domain=fem.BoundarySides(geo)) - top_velocity = wp.vec2(top_velocity, 0.0) - fem.interpolate(constant_form, dest=f_boundary, values={"val": top_velocity}) + self._bd_field = fem.UniformField(domain=fem.BoundarySides(geo), value=wp.vec2(top_velocity, 0.0)) self.renderer = fem_example_utils.Plot() @@ -132,9 +122,7 @@ def step(self): # Weak velocity boundary conditions u_bd_test = fem.make_test(space=u_space, domain=boundary) u_bd_trial = fem.make_trial(space=u_space, domain=boundary) - u_rhs = fem.integrate( - top_mass_form, fields={"u": self._bd_field.trace(), "v": u_bd_test}, output_dtype=wp.vec2d - ) + u_rhs = fem.integrate(top_mass_form, fields={"u": self._bd_field, "v": u_bd_test}, output_dtype=wp.vec2d) u_bd_matrix = fem.integrate(mass_form, fields={"u": u_bd_trial, "v": u_bd_test}) # Pressure-velocity coupling @@ -143,7 +131,7 @@ def step(self): # Define and solve the saddle-point system u_matrix = u_visc_matrix - sparse.bsr_axpy(x=u_bd_matrix, y=u_matrix, alpha=self.boundary_strength, beta=1.0) + u_matrix += self.boundary_strength * u_bd_matrix array_axpy(x=u_rhs, y=u_rhs, alpha=0.0, beta=self.boundary_strength) p_rhs = wp.zeros(p_space.node_count(), dtype=wp.float64) @@ -163,8 +151,8 @@ def step(self): wp.utils.array_cast(in_array=x_p, out_array=self._p_field.dof_values) def render(self): - self.renderer.add_surface("pressure", self._p_field) - self.renderer.add_surface_vector("velocity", self._u_field) + self.renderer.add_field("pressure", self._p_field) + self.renderer.add_field("velocity", self._u_field) if __name__ == "__main__": @@ -212,4 +200,4 @@ def render(self): example.render() if not args.headless: - example.renderer.plot() + example.renderer.plot(options={"velocity": {"streamlines": {}}, "pressure": {"contours": {}}}) diff --git a/warp/examples/fem/example_stokes_transfer.py b/warp/examples/fem/example_stokes_transfer.py index 34aadc3da..8569b0cc6 100644 --- a/warp/examples/fem/example_stokes_transfer.py +++ b/warp/examples/fem/example_stokes_transfer.py @@ -22,7 +22,6 @@ import warp.examples.fem.utils as fem_example_utils import warp.fem as fem from warp.fem.utils import array_axpy -from warp.sparse import bsr_axpy, bsr_mm, bsr_mv, bsr_transposed from warp.utils import array_cast @@ -174,12 +173,12 @@ def step(self): # Assemble linear system u_matrix = u_visc_matrix - bsr_axpy(u_bd_matrix, u_matrix, alpha=self.bd_strength) + u_matrix += u_bd_matrix * self.bd_strength - div_matrix_t = bsr_transposed(div_matrix) - gradient_matrix = bsr_mm(div_matrix_t, inv_p_mass_matrix) - bsr_mm(gradient_matrix, div_matrix, u_matrix, alpha=1.0 / self.compliance, beta=1.0) + gradient_matrix = div_matrix.transpose() @ inv_p_mass_matrix + u_matrix += gradient_matrix @ div_matrix / self.compliance + # scale u_rhs array_axpy(u_rhs, u_rhs, alpha=0.0, beta=self.bd_strength) # Solve for displacement @@ -187,8 +186,8 @@ def step(self): fem_example_utils.bsr_cg(u_matrix, x=u_res, b=u_rhs, quiet=self._quiet) # Compute pressure from displacement - div_u = bsr_mv(A=div_matrix, x=u_res) - p_res = bsr_mv(A=inv_p_mass_matrix, x=div_u, alpha=-1) + div_u = div_matrix @ u_res + p_res = -inv_p_mass_matrix @ div_u # Copy to fields u_nodes = wp.indexedarray(self._u_field.dof_values, indices=self._active_space_partition.space_node_indices()) @@ -198,8 +197,8 @@ def step(self): array_cast(in_array=p_res, out_array=p_nodes) def render(self): - self.renderer.add_surface("pressure", self._p_field) - self.renderer.add_surface_vector("velocity", self._u_field) + self.renderer.add_field("pressure", self._p_field) + self.renderer.add_field("velocity", self._u_field) def _gen_particles(self, circle_radius, c1_center, c2_center): """Generate some particles along two circles defining velocity boundary conditions""" @@ -252,4 +251,4 @@ def _gen_particles(self, circle_radius, c1_center, c2_center): example.render() if not args.headless: - example.renderer.plot(streamlines=["velocity"]) + example.renderer.plot(options={"velocity": {"streamlines": {}}}) diff --git a/warp/examples/fem/example_streamlines.py b/warp/examples/fem/example_streamlines.py index dacf5bd6a..d24f8015f 100644 --- a/warp/examples/fem/example_streamlines.py +++ b/warp/examples/fem/example_streamlines.py @@ -1,9 +1,24 @@ +# Copyright (c) 2024 NVIDIA CORPORATION. All rights reserved. +# NVIDIA CORPORATION and its licensors retain all intellectual property +# and proprietary rights in and to this software, related documentation +# and any modifications thereto. Any use, reproduction, disclosure or +# distribution of this software and related documentation without an express +# license agreement from NVIDIA CORPORATION is strictly prohibited. + +########################################################################### +# Example Streamlines +# +# Shows how to generate 3D streamlines by tracing through a velocity field +# using the `warp.fem.lookup` operator. +# Also illustrates using `warp.fem.Subdomain` to define subsets of elements. +# +########################################################################### + import numpy as np import warp as wp import warp.examples.fem.utils as fem_example_utils import warp.fem as fem -import warp.sparse as sp from warp.examples.fem.example_apic_fluid import divergence_form, solve_incompressibility @@ -70,14 +85,6 @@ def mass_form( return u(s) * v(s) -@fem.integrand -def velocity_norm( - s: fem.Sample, - u: fem.Field, -): - return u(s)[0] # wp.length(u(s)) - - @fem.integrand def spawn_streamlines(s: fem.Sample, domain: fem.Domain, jitter: float): rng = wp.rand_init(s.qp_index) @@ -95,6 +102,7 @@ def gen_streamlines( s: fem.Sample, domain: fem.Domain, u: fem.Field, + spawn_points: wp.array(dtype=wp.vec3), point_count: int, dx: float, pos: wp.array2d(dtype=wp.vec3), @@ -102,7 +110,8 @@ def gen_streamlines( ): idx = s.qp_index - p = domain(s) + p = spawn_points[idx] + s = fem.lookup(domain, p) for k in range(point_count): v = u(s) pos[idx, k] = p @@ -175,7 +184,8 @@ def __init__(self, quiet=False, degree=2, resolution=16, mesh="grid", headless: camera_front=(-0.66, 0.0, -1.0), draw_axis=False, ) - except Exception: + except Exception as err: + wp.utils.warn(f"Could not initialize OpenGL renderer: {err}") pass def step(self): @@ -198,7 +208,6 @@ def step(self): # here we use a fixed number of points per streamline, otherwise we would need to # do a first pass to count points, then array_scan the offsets, then a second pass # to populate the per-point data - spawn_qp = fem.PicQuadrature(fem.Cells(self._geo), positions=spawn_points) point_count = self._streamline_point_count points = wp.array(dtype=wp.vec3, shape=(n_streamlines, point_count)) @@ -206,9 +215,11 @@ def step(self): fem.interpolate( gen_streamlines, - quadrature=spawn_qp, + domain=fem.Cells(self._geo), + dim=n_streamlines, fields={"u": self.velocity_field}, values={ + "spawn_points": spawn_points, "point_count": self._streamline_point_count, "dx": self._streamline_dx, "pos": points, @@ -220,9 +231,9 @@ def step(self): self._speed = speed def render(self): - # self.renderer.add_volume("solution", self.pressure_field) - self.plot.add_volume("pressure", self.pressure_field) - self.plot.add_volume("vel_norm", self.velocity_norm_field) + # self.renderer.add_field("solution", self.pressure_field) + self.plot.add_field("pressure", self.pressure_field) + self.plot.add_field("velocity", self.velocity_field) if self.renderer is not None: streamline_count = self._points.shape[0] @@ -252,7 +263,6 @@ def _generate_incompressible_flow(self): self.pressure_field = p_space.make_field() self.velocity_field = u_space.make_field() - self.velocity_norm_field = s_space.make_field() # Boundary condition projector and matrices inflow_test = fem.make_test(u_space, domain=self._inflow) @@ -263,14 +273,11 @@ def _generate_incompressible_flow(self): freeslip_test = fem.make_test(u_space, domain=self._freeslip) freeslip_trial = fem.make_trial(u_space, domain=self._freeslip) - sp.bsr_axpy( - y=dirichlet_projector, - x=fem.integrate( - freeslip_projector_form, - fields={"u": freeslip_test, "v": freeslip_trial}, - nodal=True, - output_dtype=float, - ), + dirichlet_projector += fem.integrate( + freeslip_projector_form, + fields={"u": freeslip_test, "v": freeslip_trial}, + nodal=True, + output_dtype=float, ) fem.normalize_dirichlet_projector(dirichlet_projector) @@ -294,7 +301,7 @@ def _generate_incompressible_flow(self): output_dtype=float, ) - # Solve unilateral incompressibility + # Solve incompressibility solve_incompressibility( divergence_matrix, dirichlet_projector, @@ -304,8 +311,6 @@ def _generate_incompressible_flow(self): quiet=self._quiet, ) - fem.interpolate(velocity_norm, dest=self.velocity_norm_field, fields={"u": self.velocity_field}) - if __name__ == "__main__": import argparse @@ -335,4 +340,9 @@ def _generate_incompressible_flow(self): example.render() if not args.headless: - example.plot.plot() + example.plot.plot( + { + "velocity": {"streamlines": {"density": 2}}, + "pressure": {"contours": {}}, + } + ) diff --git a/warp/examples/fem/utils.py b/warp/examples/fem/utils.py index 6331c795f..e1ba23df4 100644 --- a/warp/examples/fem/utils.py +++ b/warp/examples/fem/utils.py @@ -1,4 +1,4 @@ -from typing import Any, Optional, Set, Tuple +from typing import Any, Dict, Optional, Tuple import numpy as np @@ -465,313 +465,444 @@ def _block_diagonal_invert(values: wp.array(dtype=Any)): # -def _plot_grid_surface(field, axes=None): - import matplotlib.pyplot as plt - from matplotlib import cm +class Plot: + def __init__(self, stage=None, default_point_radius=0.01): + self.default_point_radius = default_point_radius - if axes is None: - fig, axes = plt.subplots(subplot_kw={"projection": "3d"}) + self._fields = {} - node_positions = field.space.node_grid() + self._usd_renderer = None + if stage is not None: + try: + from warp.render import UsdRenderer - # Make data. - X = node_positions[0] - Y = node_positions[1] - Z = field.dof_values.numpy().reshape(X.shape) + self._usd_renderer = UsdRenderer(stage) + except Exception as err: + print(f"Could not initialize UsdRenderer for stage '{stage}': {err}.") - # Plot the surface. - return axes.plot_surface(X, Y, Z, cmap=cm.coolwarm, linewidth=0, antialiased=False) + def begin_frame(self, time): + if self._usd_renderer is not None: + self._usd_renderer.begin_frame(time=time) + def end_frame(self): + if self._usd_renderer is not None: + self._usd_renderer.end_frame() -def _plot_tri_surface(field, axes=None): - import matplotlib.pyplot as plt - from matplotlib import cm - from matplotlib.tri.triangulation import Triangulation + def add_field(self, name: str, field: fem.DiscreteField): + if self._usd_renderer is not None: + self._render_to_usd(field) - if axes is None: - fig, axes = plt.subplots(subplot_kw={"projection": "3d"}) + if name not in self._fields: + field_clone = field.space.make_field(space_partition=field.space_partition) + self._fields[name] = (field_clone, []) - node_positions = field.space.node_positions().numpy() + self._fields[name][1].append(field.dof_values.numpy()) - triangulation = Triangulation( - x=node_positions[:, 0], y=node_positions[:, 1], triangles=field.space.node_triangulation() - ) + def _render_to_usd(self, name: str, field: fem.DiscreteField): + points = field.space.node_positions().numpy() + values = field.dof_values.numpy() - Z = field.dof_values.numpy() + if values.ndim == 2: + if values.shape[1] == field.space.dimension: + # use values as displacement + points += values + else: + # use magnitude + values = np.linalg.norm(values, axis=1) - # Plot the surface. - return axes.plot_trisurf(triangulation, Z, cmap=cm.coolwarm, linewidth=0, antialiased=False) + if field.space.dimension == 2: + z = values if values.ndim == 1 else np.zeros((points.shape[0], 1)) + points = np.hstack((points, z)) + if hasattr(field.space, "node_triangulation"): + indices = field.space.node_triangulation() + self._usd_renderer.render_mesh(name, points=points, indices=indices) + else: + self._usd_renderer.render_points(name, points=points, radius=self.default_point_radius) + elif values.ndim == 1: + self._usd_renderer.render_points(name, points, radius=values) + else: + self._usd_renderer.render_points(name, points, radius=self.default_point_radius) -def _plot_tri_mesh(field, axes=None, **kwargs): - import matplotlib.pyplot as plt - from matplotlib.tri.triangulation import Triangulation + def plot(self, options: Dict[str, Any] = None, backend: str = "auto"): + if options is None: + options = {} - if axes is None: - fig, axes = plt.subplots() + if backend == "pyvista": + return self._plot_pyvista(options) + if backend == "matplotlib": + return self._plot_matplotlib(options) - vtx_positions = field.space.node_positions().numpy() - displacement = field.dof_values.numpy() + # try both + try: + return self._plot_pyvista(options) + except ModuleNotFoundError: + try: + return self._plot_matplotlib(options) + except ModuleNotFoundError: + wp.utils.warn("pyvista or matplotlib must be installed to visualize solution results") - X = vtx_positions[:, 0] + displacement[:, 0] - Y = vtx_positions[:, 1] + displacement[:, 1] + def _plot_pyvista(self, options: Dict[str, Any]): + import pyvista - triangulation = Triangulation(x=X, y=Y, triangles=field.space.node_triangulation()) + grids = {} + scales = {} + markers = {} - # Plot the surface. - return axes.triplot(triangulation, **kwargs)[0] + animate = False + for name, (field, values) in self._fields.items(): + cells, types = field.space.vtk_cells() + node_pos = field.space.node_positions().numpy() -def _plot_scatter_surface(field, axes=None): - import matplotlib.pyplot as plt - from matplotlib import cm + args = options.get(name, {}) - if axes is None: - fig, axes = plt.subplots(subplot_kw={"projection": "3d"}) + grid_scale = np.max(np.max(node_pos, axis=0) - np.min(node_pos, axis=0)) + value_range = self._get_field_value_range(values, args) + scales[name] = (grid_scale, value_range) - X, Y = field.space.node_positions().numpy().T + if node_pos.shape[1] == 2: + node_pos = np.hstack((node_pos, np.zeros((node_pos.shape[0], 1)))) - # Make data. - Z = field.dof_values.numpy().reshape(X.shape) + grid = pyvista.UnstructuredGrid(cells, types, node_pos) + grids[name] = grid - # Plot the surface. - return axes.scatter(X, Y, Z, c=Z, cmap=cm.coolwarm) + if len(values) > 1: + animate = True + + def set_frame_data(frame): + for name, (field, values) in self._fields.items(): + if frame > 0 and len(values) == 1: + continue + + v = values[frame % len(values)] + grid = grids[name] + grid_scale, value_range = scales[name] + field_args = options.get(name, {}) + + marker = None + + if field.space.dimension == 2 and v.ndim == 2 and v.shape[1] == 2: + grid.point_data[name] = np.hstack((v, np.zeros((v.shape[0], 1)))) + else: + grid.point_data[name] = v + + if v.ndim == 2: + grid.point_data[name + "_mag"] = np.linalg.norm(v, axis=1) + + if "arrows" in field_args: + glyph_scale = field_args["arrows"].get("glyph_scale", 1.0) + glyph_scale *= grid_scale / max(1.0e-8, value_range[1] - value_range[0]) + marker = grid.glyph(scale=name, orient=name, factor=glyph_scale) + elif "contours" in field_args: + levels = field_args["contours"].get("levels", 10) + if type(levels) == int: + levels = np.linspace(*value_range, levels) + marker = grid.contour(isosurfaces=levels, scalars=name + "_mag" if v.ndim == 2 else name) + elif field.space.dimension == 2: + z_scale = grid_scale / max(1.0e-8, value_range[1] - value_range[0]) + + if "streamlines" in field_args: + center = np.mean(grid.points, axis=0) + density = field_args["streamlines"].get("density", 1.0) + cell_size = 1.0 / np.sqrt(field.space.geometry.cell_count()) + + separating_distance = 0.5 / (30.0 * density * cell_size) + # Try with various sep distance until we get at least one line + while separating_distance * cell_size < 1.0: + lines = grid.streamlines_evenly_spaced_2D( + vectors=name, + start_position=center, + separating_distance=separating_distance, + separating_distance_ratio=0.5, + step_length=0.25, + compute_vorticity=False, + ) + if lines.n_lines > 0: + break + separating_distance *= 1.25 + marker = lines.tube(radius=0.0025 * grid_scale / density) + elif "arrows" in field_args: + glyph_scale = field_args["arrows"].get("glyph_scale", 1.0) + glyph_scale *= grid_scale / max(1.0e-8, value_range[1] - value_range[0]) + marker = grid.glyph(scale=name, orient=name, factor=glyph_scale) + elif "displacement" in field_args: + grid.points[:, 0:2] = field.space.node_positions().numpy() + v + else: + # Extrude surface + z = v if v.ndim == 1 else grid.point_data[name + "_mag"] + grid.points[:, 2] = z * z_scale + + elif field.space.dimension == 3: + if "streamlines" in field_args: + center = np.mean(grid.points, axis=0) + density = field_args["streamlines"].get("density", 1.0) + cell_size = 1.0 / np.sqrt(field.space.geometry.cell_count()) + lines = grid.streamlines(vectors=name, n_points=int(100 * density)) + marker = lines.tube(radius=0.0025 * grid_scale / density) + elif "displacement" in field_args: + grid.points = field.space.node_positions().numpy() + v + + if frame == 0: + if v.ndim == 1: + grid.set_active_scalars(name) + else: + grid.set_active_vectors(name) + grid.set_active_scalars(name + "_mag") + markers[name] = marker + elif marker: + markers[name].copy_from(marker) + + set_frame_data(0) + + subplot_rows = options.get("rows", 1) + subplot_shape = (subplot_rows, (len(grids) + subplot_rows - 1) // subplot_rows) + + plotter = pyvista.Plotter(shape=subplot_shape) + plotter.link_views() + plotter.add_camera_orientation_widget() + for index, (name, grid) in enumerate(grids.items()): + plotter.subplot(index // subplot_shape[1], index % subplot_shape[1]) + grid_scale, value_range = scales[name] + field = self._fields[name][0] + marker = markers[name] + if marker: + if field.space.dimension == 2: + plotter.add_mesh(marker, show_scalar_bar=False) + plotter.add_mesh(grid, opacity=0.25, clim=value_range) + plotter.view_xy() + else: + plotter.add_mesh(marker, opacity=0.25) + elif field.space.dimension == 3: + plotter.add_mesh_clip_plane(grid, show_edges=True, clim=value_range) + else: + plotter.add_mesh(grid, show_edges=True, clim=value_range) + plotter.show(interactive_update=animate) + frame = 0 + while animate and not plotter.iren.interactor.GetDone(): + frame += 1 + set_frame_data(frame) + plotter.update() -def _plot_surface(field, axes=None): - if hasattr(field.space, "node_grid"): - return _plot_grid_surface(field, axes) - elif hasattr(field.space, "node_triangulation"): - return _plot_tri_surface(field, axes) - else: - return _plot_scatter_surface(field, axes) + def _plot_matplotlib(self, options: Dict[str, Any]): + import matplotlib.animation as animation + import matplotlib.pyplot as plt + from matplotlib import cm + def make_animation(fig, ax, cax, values, draw_func): + def animate(i): + cs = draw_func(ax, values[i]) -def _plot_grid_color(field, axes=None): - import matplotlib.pyplot as plt - from matplotlib import cm + cax.cla() + fig.colorbar(cs, cax) - if axes is None: - fig, axes = plt.subplots() + return cs - node_positions = field.space.node_grid() + return animation.FuncAnimation( + ax.figure, + animate, + interval=30, + blit=False, + frames=len(values), + ) - # Make data. - X = node_positions[0] - Y = node_positions[1] - Z = field.dof_values.numpy().reshape(X.shape) + def make_draw_func(field, args, plot_func, plot_opts): + def draw_fn(axes, values): + axes.clear() - # Plot the surface. - return axes.pcolormesh(X, Y, Z, cmap=cm.coolwarm) + field.dof_values = values + cs = plot_func(field, axes=axes, **plot_opts) + if "xlim" in args: + axes.set_xlim(*args["xlim"]) + if "ylim" in args: + axes.set_ylim(*args["ylim"]) -def _plot_velocities(field, axes=None): - import matplotlib.pyplot as plt + return cs - if axes is None: - fig, axes = plt.subplots() + return draw_fn - node_positions = field.space.node_positions().numpy() + anims = [] - # Make data. - X = node_positions[:, 0] - Y = node_positions[:, 1] + field_count = len(self._fields) + subplot_rows = options.get("rows", 1) + subplot_shape = (subplot_rows, (field_count + subplot_rows - 1) // subplot_rows) - vel = field.dof_values.numpy() - u = np.ascontiguousarray(vel[:, 0]) - v = np.ascontiguousarray(vel[:, 1]) + for index, (name, (field, values)) in enumerate(self._fields.items()): + args = options.get(name, {}) + v = values[0] - u = u.reshape(X.shape) - v = v.reshape(X.shape) + plot_fn = None + plot_3d = False + plot_opts = {"cmap": cm.viridis} - return axes.quiver(X, Y, u, v) + plot_opts["clim"] = self._get_field_value_range(values, args) + if field.space.dimension == 2: + if "contours" in args: + plot_opts["levels"] = args["contours"].get("levels", None) + plot_fn = _plot_contours + elif v.ndim == 2 and v.shape[1] == 2: + if "displacement" in args: + plot_fn = _plot_displaced_tri_mesh + elif "streamlines" in args: + plot_opts["density"] = args["streamlines"].get("density", 1.0) + plot_fn = _plot_streamlines + elif "arrows" in args: + plot_opts["glyph_scale"] = args["arrows"].get("glyph_scale", 1.0) + plot_fn = _plot_quivers -def plot_grid_streamlines(field, axes=None): - import matplotlib.pyplot as plt + if plot_fn is None: + plot_fn = _plot_surface + plot_3d = True - if axes is None: - fig, axes = plt.subplots() + elif field.space.dimension == 3: + if "arrows" in args or "streamlines" in args: + plot_opts["glyph_scale"] = args.get("arrows", {}).get("glyph_scale", 1.0) + plot_fn = _plot_quivers_3d + else: + plot_fn = _plot_3d_scatter + plot_3d = True - node_positions = field.space.node_grid() + subplot_kw = {"projection": "3d"} if plot_3d else {} + axes = plt.subplot(*subplot_shape, index + 1, **subplot_kw) - # Make data. - X = node_positions[0][:, 0] - Y = node_positions[1][0, :] + if not plot_3d: + axes.set_aspect("equal") - vel = field.dof_values.numpy() - u = np.ascontiguousarray(vel[:, 0]) - v = np.ascontiguousarray(vel[:, 1]) + draw_fn = make_draw_func(field, args, plot_func=plot_fn, plot_opts=plot_opts) + cs = draw_fn(axes, values[0]) - u = np.transpose(u.reshape(node_positions[0].shape)) - v = np.transpose(v.reshape(node_positions[0].shape)) + fig = plt.gcf() + cax = fig.colorbar(cs).ax - splot = axes.streamplot(X, Y, u, v, density=2) - splot.axes = axes - return splot + if len(values) > 1: + anims.append(make_animation(fig, axes, cax, values, draw_func=draw_fn)) + plt.show() -def plot_3d_scatter(field, axes=None): - import matplotlib.pyplot as plt - from matplotlib import cm + @staticmethod + def _get_field_value_range(values, field_options: Dict[str, Any]): + value_range = field_options.get("clim", None) + if value_range is None: + value_range = ( + min((np.min(_value_or_magnitude(v)) for v in values)), + max((np.max(_value_or_magnitude(v)) for v in values)), + ) - if axes is None: - fig, axes = plt.subplots(subplot_kw={"projection": "3d"}) + return value_range - X, Y, Z = field.space.node_positions().numpy().T - # Make data. - f = field.dof_values.numpy().reshape(X.shape) +def _value_or_magnitude(values: np.ndarray): + if values.ndim == 1: + return values + return np.linalg.norm(values, axis=-1) - # Plot the surface. - return axes.scatter(X, Y, Z, c=f, cmap=cm.coolwarm) +def _field_triangulation(field): + from matplotlib.tri import Triangulation -def plot_3d_velocities(field, axes=None): - import matplotlib.pyplot as plt + node_positions = field.space.node_positions().numpy() + return Triangulation(x=node_positions[:, 0], y=node_positions[:, 1], triangles=field.space.node_triangulation()) - if axes is None: - fig, axes = plt.subplots(subplot_kw={"projection": "3d"}) - X, Y, Z = field.space.node_positions().numpy().T +def _plot_surface(field, axes, **kwargs): + Z = _value_or_magnitude(field.dof_values.numpy()) - vel = field.dof_values.numpy() - u = np.ascontiguousarray(vel[:, 0]) - v = np.ascontiguousarray(vel[:, 1]) - w = np.ascontiguousarray(vel[:, 2]) + if "clim" in kwargs: + axes.set_zlim(*kwargs["clim"]) - u = u.reshape(X.shape) - v = v.reshape(X.shape) - w = w.reshape(X.shape) + if hasattr(field.space, "node_grid"): + X, Y = field.space.node_grid() + Z = Z.reshape(X.shape) + return axes.plot_surface(X, Y, Z, linewidth=0.1, antialiased=False, **kwargs) - return axes.quiver(X, Y, Z, u, v, w, length=1.0 / X.shape[0], normalize=False) + if hasattr(field.space, "node_triangulation"): + triangulation = _field_triangulation(field) + return axes.plot_trisurf(triangulation, Z, linewidth=0.1, antialiased=False, **kwargs) + # scatter + X, Y = field.space.node_positions().numpy().T + return axes.scatter(X, Y, Z, c=Z, **kwargs) -class Plot: - def __init__(self, stage=None, default_point_radius=0.01): - self.default_point_radius = default_point_radius - self._surfaces = {} - self._surface_vectors = {} - self._volumes = {} +def _plot_displaced_tri_mesh(field, axes, **kwargs): + triangulation = _field_triangulation(field) - self._usd_renderer = None - if stage is not None: - try: - from warp.render import UsdRenderer + displacement = field.dof_values.numpy() + triangulation.x += displacement[:, 0] + triangulation.y += displacement[:, 1] - self._usd_renderer = UsdRenderer(stage) - except Exception as err: - print(f"Could not initialize UsdRenderer for stage '{stage}': {err}.") + Z = _value_or_magnitude(displacement) - def begin_frame(self, time): - if self._usd_renderer is not None: - self._usd_renderer.begin_frame(time=time) + # Plot the surface. + cs = axes.tripcolor(triangulation, Z, **kwargs) + axes.triplot(triangulation, lw=0.1) - def end_frame(self): - if self._usd_renderer is not None: - self._usd_renderer.end_frame() + return cs - def add_surface(self, name: str, field: fem.DiscreteField): - if self._usd_renderer is not None: - points_2d = field.space.node_positions().numpy() - values = field.dof_values.numpy() - points_3d = np.hstack((points_2d, values.reshape(-1, 1))) - if hasattr(field.space, "node_triangulation"): - indices = field.space.node_triangulation() - self._usd_renderer.render_mesh(name, points=points_3d, indices=indices) - else: - self._usd_renderer.render_points(name, points=points_3d, radius=self.default_point_radius) +def _plot_quivers(field, axes, clim=None, glyph_scale=1.0, **kwargs): + X, Y = field.space.node_positions().numpy().T - if name not in self._surfaces: - field_clone = field.space.make_field(space_partition=field.space_partition) - self._surfaces[name] = (field_clone, []) + vel = field.dof_values.numpy() + u = vel[:, 0].reshape(X.shape) + v = vel[:, 1].reshape(X.shape) - self._surfaces[name][1].append(field.dof_values.numpy()) + return axes.quiver(X, Y, u, v, _value_or_magnitude(vel), scale=1.0 / glyph_scale, **kwargs) - def add_surface_vector(self, name: str, field: fem.DiscreteField): - if self._usd_renderer is not None: - points_2d = field.space.node_positions().numpy() - values = field.dof_values.numpy() - points_3d = np.hstack((points_2d + values, np.zeros_like(points_2d[:, 0]).reshape(-1, 1))) - if hasattr(field.space, "node_triangulation"): - indices = field.space.node_triangulation() - self._usd_renderer.render_mesh(name, points=points_3d, indices=indices) - else: - self._usd_renderer.render_points(name, points=points_3d, radius=self.default_point_radius) +def _plot_quivers_3d(field, axes, clim=None, cmap=None, glyph_scale=1.0, **kwargs): + X, Y, Z = field.space.node_positions().numpy().T - if name not in self._surface_vectors: - field_clone = field.space.make_field(space_partition=field.space_partition) - self._surface_vectors[name] = (field_clone, []) + vel = field.dof_values.numpy() - self._surface_vectors[name][1].append(field.dof_values.numpy()) + colors = cmap((_value_or_magnitude(vel) - clim[0]) / (clim[1] - clim[0])) - def add_volume(self, name: str, field: fem.DiscreteField): - if self._usd_renderer is not None: - points_3d = field.space.node_positions().numpy() - values = field.dof_values.numpy() + u = vel[:, 0].reshape(X.shape) / (clim[1] - clim[0]) + v = vel[:, 1].reshape(X.shape) / (clim[1] - clim[0]) + w = vel[:, 2].reshape(X.shape) / (clim[1] - clim[0]) - self._usd_renderer.render_points(name, points_3d, radius=values) + return axes.quiver(X, Y, Z, u, v, w, colors=colors, length=glyph_scale, clim=clim, cmap=cmap, **kwargs) - if name not in self._volumes: - field_clone = field.space.make_field(space_partition=field.space_partition) - self._volumes[name] = (field_clone, []) - self._volumes[name][1].append(field.dof_values.numpy()) +def _plot_streamlines(field, axes, clim=None, **kwargs): + import matplotlib.tri as tr - def plot(self, streamlines: Set[str] = None, displacement: str = None): - if streamlines is None: - streamlines = set() - return self._plot_matplotlib(streamlines, displacement) + triangulation = _field_triangulation(field) - def _plot_matplotlib(self, streamlines: Set[str], displacement: str): - import matplotlib.animation as animation - import matplotlib.pyplot as plt + vel = field.dof_values.numpy() - if streamlines is None: - streamlines = [] + itp_vx = tr.CubicTriInterpolator(triangulation, vel[:, 0]) + itp_vy = tr.CubicTriInterpolator(triangulation, vel[:, 1]) - def make_animation(ax, field, values, plot_func, num_frames: int): - def animate(i): - ax.clear() - field.dof_values = values[i] - return plot_func(field, axes=ax) + X, Y = np.meshgrid( + np.linspace(np.min(triangulation.x), np.max(triangulation.x), 100), + np.linspace(np.min(triangulation.y), np.max(triangulation.y), 100), + ) - return animation.FuncAnimation( - ax.figure, - animate, - interval=30, - blit=False, - frames=len(values), - ) + u = itp_vx(X, Y) + v = itp_vy(X, Y) + C = np.sqrt(u * u + v * v) - for _name, (field, values) in self._surfaces.items(): - field.dof_values = values[0] - ax = _plot_surface(field).axes + plot = axes.streamplot(X, Y, u, v, color=C, **kwargs) + return plot.lines - if len(values) > 1: - _anim = make_animation(ax, field, values, plot_func=_plot_surface, num_frames=len(values)) - - for name, (field, values) in self._surface_vectors.items(): - field.dof_values = values[0] - if name == displacement: - ax = _plot_tri_mesh(field).axes - - if len(values) > 1: - _anim = make_animation(ax, field, values, plot_func=_plot_tri_mesh, num_frames=len(values)) - elif name in streamlines and hasattr(field.space, "node_grid"): - ax = plot_grid_streamlines(field).axes - ax.set_axis_off() - else: - ax = _plot_velocities(field).axes - if len(values) > 1: - _anim = make_animation(ax, field, values, plot_func=_plot_velocities, num_frames=len(values)) +def _plot_contours(field, axes, clim=None, **kwargs): + triangulation = _field_triangulation(field) - for _name, (field, values) in self._volumes.items(): - field.dof_values = values[0] - ax = plot_3d_scatter(field).axes + Z = _value_or_magnitude(field.dof_values.numpy()) - plt.show() + tc = axes.tricontourf(triangulation, Z, **kwargs) + axes.tricontour(triangulation, Z, **kwargs) + return tc + + +def _plot_3d_scatter(field, axes, **kwargs): + X, Y, Z = field.space.node_positions().numpy().T + + f = _value_or_magnitude(field.dof_values.numpy()).reshape(X.shape) + + return axes.scatter(X, Y, Z, c=f, **kwargs) diff --git a/warp/fem/__init__.py b/warp/fem/__init__.py index 38bad00f9..c4032de2c 100644 --- a/warp/fem/__init__.py +++ b/warp/fem/__init__.py @@ -1,7 +1,17 @@ from .cache import TemporaryStore, borrow_temporary, borrow_temporary_like, set_default_temporary_store from .dirichlet import normalize_dirichlet_projector, project_linear_system from .domain import BoundarySides, Cells, FrontierSides, GeometryDomain, Sides, Subdomain -from .field import DiscreteField, FieldLike, make_restriction, make_test, make_trial +from .field import ( + DiscreteField, + FieldLike, + ImplicitField, + NonconformingField, + UniformField, + make_discrete_field, + make_restriction, + make_test, + make_trial, +) from .geometry import ( ExplicitGeometryPartition, Geometry, diff --git a/warp/fem/cache.py b/warp/fem/cache.py index 366a0abb3..f5e684da4 100644 --- a/warp/fem/cache.py +++ b/warp/fem/cache.py @@ -1,3 +1,4 @@ +import ast import bisect import re import weakref @@ -18,7 +19,7 @@ def _make_key(obj, suffix: str, use_qualified_name): return _key_re.sub("", f"{base_name}_{suffix}") -def get_func(func, suffix: str, use_qualified_name: bool = False): +def get_func(func, suffix: str, use_qualified_name: bool = False, code_transformers=None): key = _make_key(func, suffix, use_qualified_name) if key not in _func_cache: @@ -29,14 +30,15 @@ def get_func(func, suffix: str, use_qualified_name: bool = False): module=wp.get_module( func.__module__, ), + code_transformers=code_transformers, ) return _func_cache[key] -def dynamic_func(suffix: str, use_qualified_name=False): +def dynamic_func(suffix: str, use_qualified_name=False, code_transformers=None): def wrap_func(func: Callable): - return get_func(func, suffix=suffix, use_qualified_name=use_qualified_name) + return get_func(func, suffix=suffix, use_qualified_name=use_qualified_name, code_transformers=code_transformers) return wrap_func @@ -97,6 +99,92 @@ def wrap_struct(struct: type): return wrap_struct +def get_argument_struct(arg_types: Dict[str, type]): + class Args: + pass + + annotations = wp.codegen.get_annotations(Args) + + for name, arg_type in arg_types.items(): + setattr(Args, name, None) + annotations[name] = arg_type + + def arg_type_name(arg_type): + return wp.types.get_type_code(wp.types.type_to_warp(arg_type)) + + try: + Args.__annotations__ = annotations + except AttributeError: + Args.__dict__.__annotations__ = annotations + + suffix = "_".join([f"{name}_{arg_type_name(arg_type)}" for name, arg_type in annotations.items()]) + + return get_struct(Args, suffix=suffix) + + +def populate_argument_struct(Args: wp.codegen.Struct, values: Dict[str, Any], func_name: str): + if values is None: + values = {} + + value_struct_values = Args() + for k, v in values.items(): + try: + setattr(value_struct_values, k, v) + except Exception as err: + if k not in Args.vars: + raise ValueError( + f"Passed value argument '{k}' does not match any of the function '{func_name}' parameters" + ) from err + raise ValueError( + f"Passed value argument '{k}' of type '{wp.types.type_repr(v)}' is incompatible with the function '{func_name}' parameter of type '{wp.types.type_repr(Args.vars[k].type)}'" + ) from err + + missing_values = Args.vars.keys() - values.keys() + if missing_values: + wp.utils.warn( + f"Missing values for parameter(s) '{', '.join(missing_values)}' of the function '{func_name}', will be zero-initialized" + ) + + return value_struct_values + + +class ExpandStarredArgumentStruct(ast.NodeTransformer): + def __init__( + self, + structs: Dict[str, wp.codegen.Struct], + ): + self._structs = structs + + @staticmethod + def _build_path(path, node): + if isinstance(node, ast.Attribute): + ExpandStarredArgumentStruct._build_path(path, node.value) + path.append(node.attr) + if isinstance(node, ast.Name): + path.append(node.id) + return path + + def _get_expanded_struct(self, arg_node): + if not isinstance(arg_node, ast.Starred): + return None + path = ".".join(ExpandStarredArgumentStruct._build_path([], arg_node.value)) + return self._structs.get(path, None) + + def visit_Call(self, call: ast.Call): + call = self.generic_visit(call) + + expanded_args = [] + for arg in call.args: + struct = self._get_expanded_struct(arg) + if struct is None: + expanded_args.append(arg) + else: + expanded_args += [ast.Attribute(value=arg.value, attr=field) for field in struct.vars.keys()] + call.args = expanded_args + + return call + + def get_integrand_function( integrand: "warp.fem.operator.Integrand", # noqa: F821 suffix: str, @@ -104,9 +192,6 @@ def get_integrand_function( annotations=None, code_transformers=None, ): - if code_transformers is None: - code_transformers = [] - key = _make_key(integrand.func, suffix, use_qualified_name=True) if key not in _func_cache: @@ -132,9 +217,6 @@ def get_integrand_kernel( if kernel_options is None: kernel_options = {} - if code_transformers is None: - code_transformers = [] - key = _make_key(integrand.func, suffix, use_qualified_name=True) if key not in _kernel_cache: diff --git a/warp/fem/domain.py b/warp/fem/domain.py index d842e9b04..93c0d65db 100644 --- a/warp/fem/domain.py +++ b/warp/fem/domain.py @@ -1,4 +1,3 @@ -from enum import Enum from typing import Optional, Union import warp as wp @@ -12,6 +11,7 @@ GeometryPartition, WholeGeometryPartition, ) +from warp.fem.types import ElementKind GeometryOrPartition = Union[Geometry, GeometryPartition] @@ -19,12 +19,6 @@ class GeometryDomain: """Interface class for domains, i.e. (partial) views of elements in a Geometry""" - class ElementKind(Enum): - """Possible kinds of elements contained in a domain""" - - CELL = 0 - SIDE = 1 - def __init__(self, geometry: GeometryOrPartition): if isinstance(geometry, GeometryPartition): self.geometry_partition = geometry @@ -108,8 +102,8 @@ def __init__(self, geometry: GeometryOrPartition): super().__init__(geometry) @property - def element_kind(self) -> GeometryDomain.ElementKind: - return GeometryDomain.ElementKind.CELL + def element_kind(self) -> ElementKind: + return ElementKind.CELL @property def dimension(self) -> int: @@ -177,8 +171,8 @@ def __init__(self, geometry: GeometryOrPartition): self.element_lookup = None @property - def element_kind(self) -> GeometryDomain.ElementKind: - return GeometryDomain.ElementKind.SIDE + def element_kind(self) -> ElementKind: + return ElementKind.SIDE @property def dimension(self) -> int: @@ -326,7 +320,7 @@ def __eq__(self, other) -> bool: ) @property - def element_kind(self) -> GeometryDomain.ElementKind: + def element_kind(self) -> ElementKind: return self._domain.element_kind @property diff --git a/warp/fem/field/__init__.py b/warp/fem/field/__init__.py index b4954549c..327bffb03 100644 --- a/warp/fem/field/__init__.py +++ b/warp/fem/field/__init__.py @@ -3,7 +3,7 @@ from warp.fem.domain import Cells, GeometryDomain from warp.fem.space import FunctionSpace, SpacePartition, SpaceRestriction, make_space_partition, make_space_restriction -from .field import DiscreteField, FieldLike, SpaceField +from .field import DiscreteField, FieldLike, GeometryField, ImplicitField, NonconformingField, SpaceField, UniformField from .nodal_field import NodalField from .restriction import FieldRestriction from .test import TestField @@ -85,8 +85,8 @@ def make_trial( """ if space_restriction is not None: - domain = space.domain - space_partition = space.space_partition + domain = space_restriction.domain + space_partition = space_restriction.space_partition if space_partition is None: if domain is None: @@ -98,3 +98,14 @@ def make_trial( domain = Cells(geometry=space_partition.geo_partition) return TrialField(space, space_partition, domain) + + +def make_discrete_field( + space: FunctionSpace, + space_partition: Optional[SpacePartition] = None, +) -> DiscreteField: + """Constructs a zero-initialized discrete field over a function space or partition + + See also: :meth:`warp.fem.FunctionSpace.make_field` + """ + return space.make_field(space_partition=space_partition) diff --git a/warp/fem/field/field.py b/warp/fem/field/field.py index c33864cc0..f10991002 100644 --- a/warp/fem/field/field.py +++ b/warp/fem/field/field.py @@ -1,9 +1,12 @@ -from typing import Any +from typing import Any, Dict, Optional import warp as wp +from warp.fem import cache +from warp.fem.domain import GeometryDomain, Sides from warp.fem.geometry import DeformedGeometry, Geometry +from warp.fem.operator import integrand from warp.fem.space import FunctionSpace, SpacePartition -from warp.fem.types import Sample +from warp.fem.types import NULL_ELEMENT_INDEX, ElementKind, Sample class FieldLike: @@ -36,6 +39,14 @@ def eval_arg_value(self, device): """Value of arguments to be passed to device functions""" raise NotImplementedError + def gradient_valid(self) -> bool: + """Whether the gradient operator is implemented for this field.""" + return False + + def divergence_valid(self) -> bool: + """Whether the divergence operator is implemented for this field.""" + return False + @staticmethod def eval_inner(args: "ElementEvalArg", s: "Sample"): # noqa: F821 """Device function evaluating the inner field value at a sample point""" @@ -66,14 +77,64 @@ def eval_div_outer(args: "ElementEvalArg", s: "Sample"): # noqa: F821 """Device function evaluating the outer field divergence at a sample point""" raise NotImplementedError + @staticmethod + def eval_degree(args: "ElementEvalArg"): # noqa: F821 + """Polynomial degree of the field is applicable, or hint for determination of interpolation order""" + raise NotImplementedError + -class SpaceField(FieldLike): +class GeometryField(FieldLike): + """Base class for fields defined over a geometry""" + + @property + def geometry(self) -> Geometry: + """Geometry over which the field is expressed""" + raise NotImplementedError + + @property + def element_kind(self) -> ElementKind: + """Kind of element over which the field is expressed""" + raise NotImplementedError + + @staticmethod + def eval_reference_grad_inner(args: "ElementEvalArg", s: "Sample"): # noqa: F821 + """Device function evaluating the inner field gradient with respect to reference element coordinates at a sample point""" + raise NotImplementedError + + @staticmethod + def eval_reference_grad_outer(args: "ElementEvalArg", s: "Sample"): # noqa: F821 + """Device function evaluating the outer field gradient with respect to reference element coordinates at a sample point""" + raise NotImplementedError + + def trace(self) -> FieldLike: + """Trace of this field over lower-dimensional elements""" + raise NotImplementedError + + def make_deformed_geometry(self, relative=True) -> Geometry: + """Returns a deformed version of the underlying geometry, with positions displaced according to this field's values. + + Args: + relative: If ``True``, the field is intepreted as a relative displacement over the original geometry. If ``False``, the field values are intepreted as absolute positions. + + """ + return DeformedGeometry(self, relative=relative) + + +class SpaceField(GeometryField): """Base class for fields defined over a function space""" def __init__(self, space: FunctionSpace, space_partition: SpacePartition): self._space = space self._space_partition = space_partition + @property + def geometry(self) -> Geometry: + return self._space.geometry + + @property + def element_kind(self) -> ElementKind: + return self._space.element_kind + @property def space(self) -> FunctionSpace: return self._space @@ -108,7 +169,6 @@ def divergence_valid(self) -> bool: def _make_eval_degree(self): ORDER = self.space.ORDER - from warp.fem import cache @cache.dynamic_func(suffix=self.name) def degree(args: self.ElementEvalArg): @@ -130,10 +190,6 @@ def dof_values(self, values: wp.array): """Sets degrees of freedom values from an array""" raise NotImplementedError - def trace(self) -> "DiscreteField": - """Trace of this field over a lower-dimensional function space""" - raise NotImplementedError - @staticmethod def set_node_value(args: "FieldLike.EvalArg", node_index: int, value: Any): """Device function setting the value at given node""" @@ -143,6 +199,393 @@ def set_node_value(args: "FieldLike.EvalArg", node_index: int, value: Any): def name(self) -> str: return f"{self.__class__.__qualname__}_{self.space.name}_{self.space_partition.name}" - def make_deformed_geometry(self) -> Geometry: - """Returns a deformed version of the underlying geometry using this field's values as displacement""" - return DeformedGeometry(self) + +class ImplicitField(GeometryField): + """Field defined from an arbitrary function over a domain. + Does not support autodiff yet, so if gradient/divergence evaluation is required corresponding functions must be provided. + + Args: + domain: Domain over which the field is defined + func: Warp function evaluating the field at a given position. Must accept at least one argument, with the first argument being the evaluation position (``wp.vec2`` or ``wp.vec3``). + values: Optional dictionary of additional argument values to be passed to the evaluation function. + grad_func: Optional gradient evaluation function; must take same arguments as `func` + div_func: Optional divergence evaluation function; must take same arguments as `func` + degree: Optional hint for automatic determination of quadrature orders when integrating this field + """ + + def __init__( + self, + domain: GeometryDomain, + func: wp.Function, + values: Optional[Dict[str, Any]] = None, + grad_func: Optional[wp.Function] = None, + div_func: Optional[wp.Function] = None, + degree=0, + ): + self.domain = domain + self._degree = degree + + if not isinstance(func, wp.Function): + raise ValueError("Implicit field function must be a warp Function (decorated with `wp.func`)") + + self._func = func + self._grad_func = grad_func + self._div_func = div_func + + argspec = integrand(func.func).argspec + arg_types = argspec.annotations + + pos_arg_type = arg_types.pop(argspec.args[0]) if arg_types else None + if not pos_arg_type or not wp.types.types_equal( + pos_arg_type, wp.vec(length=domain.geometry.dimension, dtype=float), match_generic=True + ): + raise ValueError( + f"Implicit field function '{func.func.__name__}' must accept a position as its first argument" + ) + + self.EvalArg = cache.get_argument_struct(arg_types) + self.values = values + + self.ElementEvalArg = self._make_element_eval_arg() + self.eval_degree = self._make_eval_degree() + + self.eval_inner = self._make_eval_func(func) + self.eval_grad_inner = self._make_eval_func(grad_func) + self.eval_div_inner = self._make_eval_func(div_func) + self.eval_reference_grad_inner = self._make_eval_reference_grad() + + self.eval_outer = self.eval_inner + self.eval_grad_outer = self.eval_grad_inner + self.eval_div_outer = self.eval_div_inner + self.eval_reference_grad_outer = self.eval_reference_grad_inner + + @property + def values(self): + return self._func_arg + + @values.setter + def values(self, v): + self._func_arg = cache.populate_argument_struct(self.EvalArg, v, self._func.func.__name__) + + @property + def geometry(self) -> Geometry: + return self.domain.geometry + + @property + def element_kind(self) -> ElementKind: + return self.domain.element_kind + + def eval_arg_value(self, device): + return self._func_arg + + @property + def degree(self) -> int: + return self._degree + + @property + def name(self) -> str: + return f"Implicit_{self.domain.name}_{self.degree}_{self.EvalArg.key}" + + def _make_eval_func(self, func): + if func is None: + return None + + @cache.dynamic_func( + suffix=f"{self.name}_{func.key}", + code_transformers=[cache.ExpandStarredArgumentStruct({"args.eval_arg": self.EvalArg})], + ) + def eval_inner(args: self.ElementEvalArg, s: Sample): + pos = self.domain.element_position(args.elt_arg, s) + return func(pos, *args.eval_arg) + + return eval_inner + + def _make_eval_reference_grad(self): + if self.eval_grad_inner is None: + return None + + @cache.dynamic_func(suffix=f"{self.eval_grad_inner.key}") + def eval_reference_grad_inner(args: self.ElementEvalArg, s: Sample): + return self.eval_grad_inner(args, s) * self.domain.element_deformation_gradient(args.elt_arg, s) + + return eval_reference_grad_inner + + def _make_element_eval_arg(self): + @cache.dynamic_struct(suffix=self.name) + class ElementEvalArg: + elt_arg: self.domain.ElementArg + eval_arg: self.EvalArg + + return ElementEvalArg + + def _make_eval_degree(self): + ORDER = wp.constant(self._degree) + + @cache.dynamic_func(suffix=self.name) + def degree(args: self.ElementEvalArg): + return ORDER + + return degree + + def trace(self): + if self.element_kind == ElementKind.SIDE: + raise RuntimeError("Trace only available for field defined on cell elements") + + return ImplicitField( + domain=Sides(self.domain.geometry_partition), + func=self._func, + values={name: getattr(self.values, name) for name in self.EvalArg.vars}, + grad_func=self._grad_func, + div_func=self._div_func, + degree=self._degree, + ) + + +class UniformField(GeometryField): + """Field defined as a constant value over a domain. + + Args: + domain: Domain over which the field is defined + value: Uniform value over the domain + """ + + def __init__(self, domain: GeometryDomain, value: Any): + self.domain = domain + + if not wp.types.is_value(value): + raise ValueError("value must be a Warp scalar, vector or matrix") + + self.dtype = wp.types.type_to_warp(type(value)) + self._value = self.dtype(value) + + scalar_type = wp.types.type_scalar_type(self.dtype) + if wp.types.type_is_vector(self.dtype): + grad_type = wp.mat(shape=(wp.types.type_length(self.dtype), self.geometry.dimension), dtype=scalar_type) + div_type = scalar_type + elif wp.types.type_is_matrix(self.dtype): + grad_type = None + div_type = wp.vec(length=(wp.types.type_length(self.dtype) // self.geometry.dimension), dtype=scalar_type) + else: + div_type = None + grad_type = wp.vec(length=self.geometry.dimension, dtype=scalar_type) + + self.EvalArg = self._make_eval_arg() + self.ElementEvalArg = self._make_element_eval_arg() + self.eval_degree = self._make_eval_degree() + + self.eval_inner = self._make_eval_inner() + self.eval_grad_inner = self._make_eval_zero(grad_type) + self.eval_div_inner = self._make_eval_zero(div_type) + self.eval_reference_grad_inner = self.eval_grad_inner + + self.eval_outer = self.eval_inner + self.eval_grad_outer = self.eval_grad_inner + self.eval_div_outer = self.eval_div_inner + self.eval_reference_grad_outer = self.eval_reference_grad_inner + + @property + def value(self): + return self._value + + @value.setter + def value(self, v): + value_type = wp.types.type_to_warp(type(v)) + assert wp.types.types_equal(value_type, self.dtype) + self._value = self.dtype(v) + + @property + def geometry(self) -> Geometry: + return self.domain.geometry + + @property + def element_kind(self) -> ElementKind: + return self.domain.element_kind + + def eval_arg_value(self, device): + arg = self.EvalArg() + arg.value = self.value + return arg + + @property + def degree(self) -> int: + return 0 + + @property + def name(self) -> str: + return f"Uniform{self.domain.name}_{wp.types.get_type_code(self.dtype)}" + + def _make_eval_inner(self): + @cache.dynamic_func(suffix=self.name) + def eval_inner(args: self.ElementEvalArg, s: Sample): + return args.eval_arg.value + + return eval_inner + + def _make_eval_zero(self, dtype): + if dtype is None: + return None + + scalar_type = wp.types.type_scalar_type(dtype) + + @cache.dynamic_func(suffix=f"{self.name}_{wp.types.get_type_code(dtype)}") + def eval_zero(args: self.ElementEvalArg, s: Sample): + return dtype(scalar_type(0.0)) + + return eval_zero + + def _make_eval_arg(self): + @cache.dynamic_struct(suffix=self.name) + class EvalArg: + value: self.dtype + + return EvalArg + + def _make_element_eval_arg(self): + @cache.dynamic_struct(suffix=self.name) + class ElementEvalArg: + elt_arg: self.domain.ElementArg + eval_arg: self.EvalArg + + return ElementEvalArg + + def _make_eval_degree(self): + @cache.dynamic_func(suffix=self.name) + def degree(args: self.ElementEvalArg): + return 0 + + return degree + + def trace(self): + if self.element_kind == ElementKind.SIDE: + raise RuntimeError("Trace only available for field defined on cell elements") + + return UniformField(domain=Sides(self.domain.geometry_partition), value=self.value) + + +class NonconformingField(GeometryField): + """Field defined as the map of a DiscreteField over a non-conforming geometry. + + Args: + domain: The new domain over which the nonconforming field will be evaluated + field: Nonconforming discrete field + background: Uniform value or domain-conforming field determining the value outside of the geometry of definition of `field` + """ + + _LOOKUP_EPS = wp.constant(1.0e-6) + + def __init__(self, domain: GeometryDomain, field: DiscreteField, background: Any = 0.0): + self.domain = domain + + self.field = field + self.dtype = field.dtype + + if not isinstance(background, GeometryField): + background = UniformField(domain, self.dtype(background)) + elif background.geometry != domain.geometry or background.element_kind != domain.element_kind: + raise ValueError("Background field must be conforming to the domain") + self.background = background + + self.EvalArg = self._make_eval_arg() + self.ElementEvalArg = self._make_element_eval_arg() + self.eval_degree = self._make_eval_degree() + + self.eval_inner = self._make_nonconforming_eval("eval_inner") + self.eval_grad_inner = self._make_nonconforming_eval("eval_grad_inner") + self.eval_div_inner = self._make_nonconforming_eval("eval_div_inner") + self.eval_reference_grad_inner = self._make_eval_reference_grad() + + # Nonconforming evaluation is position based, does not handle discontinuous fields + self.eval_outer = self.eval_inner + self.eval_grad_outer = self.eval_grad_inner + self.eval_div_outer = self.eval_div_inner + self.eval_reference_grad_outer = self.eval_reference_grad_inner + + @property + def geometry(self) -> Geometry: + return self.domain.geometry + + @property + def element_kind(self) -> ElementKind: + return self.domain.element_kind + + @cache.cached_arg_value + def eval_arg_value(self, device): + arg = self.EvalArg() + arg.field_cell_eval_arg = self.field.ElementEvalArg() + arg.field_cell_eval_arg.elt_arg = self.field.geometry.cell_arg_value(device) + arg.field_cell_eval_arg.eval_arg = self.field.eval_arg_value(device) + arg.background_arg = self.background.eval_arg_value(device) + return arg + + @property + def degree(self) -> int: + return self.field.degree + + @property + def name(self) -> str: + return f"{self.domain.name}_{self.field.name}_{self.background.name}" + + def _make_nonconforming_eval(self, eval_func_name): + field_eval = getattr(self.field, eval_func_name) + bg_eval = getattr(self.background, eval_func_name) + + if field_eval is None or bg_eval is None: + return None + + @cache.dynamic_func(suffix=f"{eval_func_name}_{self.name}") + def eval_nc(args: self.ElementEvalArg, s: Sample): + pos = self.domain.element_position(args.elt_arg, s) + cell_arg = args.eval_arg.field_cell_eval_arg.elt_arg + nonconforming_s = self.field.geometry.cell_lookup(cell_arg, pos) + if ( + nonconforming_s.element_index == NULL_ELEMENT_INDEX + or wp.length_sq(pos - self.field.geometry.cell_position(cell_arg, nonconforming_s)) + > NonconformingField._LOOKUP_EPS + ): + return bg_eval(self.background.ElementEvalArg(args.elt_arg, args.eval_arg.background_arg), s) + return field_eval( + self.field.ElementEvalArg(cell_arg, args.eval_arg.field_cell_eval_arg.eval_arg), nonconforming_s + ) + + return eval_nc + + def _make_eval_reference_grad(self): + if self.eval_grad_inner is None: + return None + + @cache.dynamic_func(suffix=f"{self.eval_grad_inner.key}") + def eval_reference_grad_inner(args: self.ElementEvalArg, s: Sample): + return self.eval_grad_inner(args, s) * self.domain.element_deformation_gradient(args.elt_arg, s) + + return eval_reference_grad_inner + + def _make_eval_arg(self): + @cache.dynamic_struct(suffix=self.name) + class EvalArg: + field_cell_eval_arg: self.field.ElementEvalArg + background_arg: self.background.EvalArg + + return EvalArg + + def _make_element_eval_arg(self): + @cache.dynamic_struct(suffix=self.name) + class ElementEvalArg: + elt_arg: self.domain.ElementArg + eval_arg: self.EvalArg + + return ElementEvalArg + + def _make_eval_degree(self): + @cache.dynamic_func(suffix=self.name) + def degree(args: self.ElementEvalArg): + return self.field.eval_degree(args.eval_arg.field_cell_eval_arg) + + return degree + + def trace(self): + if self.element_kind == ElementKind.SIDE: + raise RuntimeError("Trace only available for field defined on cell elements") + + return NonconformingField( + domain=Sides(self.domain.geometry_partition), field=self.field, background=self.background.trace() + ) diff --git a/warp/fem/geometry/deformed_geometry.py b/warp/fem/geometry/deformed_geometry.py index 49d66117b..deeb54ca9 100644 --- a/warp/fem/geometry/deformed_geometry.py +++ b/warp/fem/geometry/deformed_geometry.py @@ -10,18 +10,32 @@ class DeformedGeometry(Geometry): - def __init__(self, field): - """Constructs a Deformed Geometry from a displacement field defined over a base geometry""" - - from warp.fem.field import DiscreteField - - self.field: DiscreteField = field - self.base = self.field.space.geometry + def __init__(self, field: "wp.fem.field.GeometryField", relative: bool = True): + """Constructs a Deformed Geometry from a displacement or absolute position field defined over a base geometry. + The deformation field does not need to be isoparameteric. + + See also: :meth:`warp.fem.DiscreteField.make_deformed_geometry` + """ + + from warp.fem.field import DiscreteField, GeometryField + + if isinstance(field, DiscreteField): + if ( + not wp.types.type_is_vector(field.dtype) + or wp.types.type_length(field.dtype) != field.geometry.dimension + ): + raise ValueError( + "Invalid value type for position field, must be vector-valued with same dimension as underlying geometry" + ) + if field.eval_grad_inner is None: + raise ValueError("Gradient evaluation is not supported on the passed field") + + self._relative = relative + + self.field: GeometryField = field + self.base = self.field.geometry self.dimension = self.base.dimension - if not wp.types.type_is_vector(field.dtype) or wp.types.type_length(field.dtype) != self.dimension: - raise ValueError("Invalid value type for position field") - self.CellArg = self.field.ElementEvalArg self.field_trace = field.trace() @@ -60,7 +74,7 @@ def __init__(self, field): @property def name(self): - return f"DefGeo_{self.field.name}" + return f"DefGeo_{self.field.name}_{'rel' if self._relative else 'abs'}" # Geometry device interface @@ -74,20 +88,28 @@ def cell_arg_value(self, device) -> "DeformedGeometry.CellArg": return args def _make_cell_position(self): + @cache.dynamic_func(suffix=self.name) + def cell_position_absolute(cell_arg: self.CellArg, s: Sample): + return self.field.eval_inner(cell_arg, s) + @cache.dynamic_func(suffix=self.name) def cell_position(cell_arg: self.CellArg, s: Sample): return self.field.eval_inner(cell_arg, s) + self.base.cell_position(cell_arg.elt_arg, s) - return cell_position + return cell_position if self._relative else cell_position_absolute def _make_cell_deformation_gradient(self): + @cache.dynamic_func(suffix=self.name) + def cell_deformation_gradient_absolute(cell_arg: self.CellArg, s: Sample): + return self.field.eval_reference_grad_inner(cell_arg, s) + @cache.dynamic_func(suffix=self.name) def cell_deformation_gradient(cell_arg: self.CellArg, s: Sample): return self.field.eval_reference_grad_inner(cell_arg, s) + self.base.cell_deformation_gradient( cell_arg.elt_arg, s ) - return cell_deformation_gradient + return cell_deformation_gradient if self._relative else cell_deformation_gradient_absolute def _make_cell_inverse_deformation_gradient(self): @cache.dynamic_func(suffix=self.name) @@ -129,14 +151,27 @@ def side_arg_value(self, device) -> "DeformedGeometry.SideArg": return args def _make_side_position(self): + @cache.dynamic_func(suffix=self.name) + def side_position_absolute(args: self.SideArg, s: Sample): + trace_arg = self.field_trace.ElementEvalArg(args.base_arg, args.trace_arg) + return self.field_trace.eval_inner(trace_arg, s) + @cache.dynamic_func(suffix=self.name) def side_position(args: self.SideArg, s: Sample): trace_arg = self.field_trace.ElementEvalArg(args.base_arg, args.trace_arg) return self.field_trace.eval_inner(trace_arg, s) + self.base.side_position(args.base_arg, s) - return side_position + return side_position if self._relative else side_position_absolute def _make_side_deformation_gradient(self): + @cache.dynamic_func(suffix=self.name) + def side_deformation_gradient_absolute(args: self.SideArg, s: Sample): + base_def_grad = self.base.side_deformation_gradient(args.base_arg, s) + trace_arg = self.field_trace.ElementEvalArg(args.base_arg, args.trace_arg) + + Du = self.field_trace.eval_grad_inner(trace_arg, s) + return Du * base_def_grad + @cache.dynamic_func(suffix=self.name) def side_deformation_gradient(args: self.SideArg, s: Sample): base_def_grad = self.base.side_deformation_gradient(args.base_arg, s) @@ -145,7 +180,7 @@ def side_deformation_gradient(args: self.SideArg, s: Sample): Du = self.field_trace.eval_grad_inner(trace_arg, s) return base_def_grad + Du * base_def_grad - return side_deformation_gradient + return side_deformation_gradient if self._relative else side_deformation_gradient_absolute def _make_side_inner_inverse_deformation_gradient(self): @cache.dynamic_func(suffix=self.name) diff --git a/warp/fem/integrate.py b/warp/fem/integrate.py index 47cae3f06..570dfb721 100644 --- a/warp/fem/integrate.py +++ b/warp/fem/integrate.py @@ -9,14 +9,25 @@ DiscreteField, FieldLike, FieldRestriction, - SpaceField, + GeometryField, TestField, TrialField, make_restriction, ) -from warp.fem.operator import Integrand, Operator +from warp.fem.operator import Integrand, Operator, integrand from warp.fem.quadrature import Quadrature, RegularQuadrature -from warp.fem.types import NULL_DOF_INDEX, NULL_NODE_INDEX, OUTSIDE, DofIndex, Domain, Field, Sample, make_free_sample +from warp.fem.types import ( + NULL_DOF_INDEX, + NULL_ELEMENT_INDEX, + NULL_NODE_INDEX, + OUTSIDE, + Coords, + DofIndex, + Domain, + Field, + Sample, + make_free_sample, +) from warp.sparse import BsrMatrix, bsr_set_from_triplets, bsr_zeros from warp.types import type_length from warp.utils import array_cast @@ -108,7 +119,10 @@ def _replace_call_func(self, call: ast.Call, callee: Union[type, Operator], oper try: # Retrieve the function pointer corresponding to the operator implementation for the field type pointer = operator.resolver(field) - except AttributeError as e: + if pointer is None: + raise NotImplementedError(operator.resolver.__name__) + + except (AttributeError, NotImplementedError) as e: raise ValueError(f"Operator {operator.func.__name__} is not defined for field {field.name}") from e # Save the pointer as an attribute than can be accessed from the callee scope setattr(callee, pointer.key, pointer) @@ -209,39 +223,15 @@ def _check_field_compat( f"Passed field argument '{name}' does not match any parameter of integrand '{integrand.name}'" ) - if isinstance(field, SpaceField) and domain is not None: - space = field.space - if space.geometry != domain.geometry: + if isinstance(field, GeometryField) and domain is not None: + if field.geometry != domain.geometry: raise ValueError(f"Field '{name}' must be defined on the same geometry as the integration domain") - if space.dimension != domain.dimension: + if field.element_kind != domain.element_kind: raise ValueError( - f"Field '{name}' dimension ({space.dimension}) does not match that of the integration domain ({domain.dimension}). Maybe a forgotten `.trace()`?" + f"Field '{name}' is not defined on the same kind of elements (cells or sides) as the integration domain. Maybe a forgotten `.trace()`?" ) -def _populate_value_struct(ValueStruct: wp.codegen.Struct, values: Dict[str, Any], integrand_name: str): - value_struct_values = ValueStruct() - for k, v in values.items(): - try: - setattr(value_struct_values, k, v) - except Exception as err: - if k not in ValueStruct.vars: - raise ValueError( - f"Passed value argument '{k}' does not match any of the integrand '{integrand_name}' parameters" - ) from err - raise ValueError( - f"Passed value argument '{k}' of type '{wp.types.type_repr(v)}' is incompatible with the integrand '{integrand_name}' parameter of type '{wp.types.type_repr(ValueStruct.vars[k].type)}'" - ) from err - - missing_values = ValueStruct.vars.keys() - values.keys() - if missing_values: - wp.utils.warn( - f"Missing values for parameter(s) '{', '.join(missing_values)}' of the integrand '{integrand_name}', will be zero-initialized" - ) - - return value_struct_values - - def _get_test_and_trial_fields( fields: Dict[str, FieldLike], ): @@ -297,36 +287,6 @@ class Fields: return cache.get_struct(Fields, suffix=suffix) -def _gen_value_struct(value_args: Dict[str, type]): - class Values: - pass - - annotations = get_annotations(Values) - - for name, arg_type in value_args.items(): - setattr(Values, name, None) - annotations[name] = arg_type - - def arg_type_name(arg_type): - if isinstance(arg_type, wp.codegen.Struct): - return arg_type_name(arg_type.cls) - return getattr(arg_type, "__name__", str(arg_type)) - - def arg_type_name(arg_type): - if isinstance(arg_type, wp.codegen.Struct): - return arg_type_name(arg_type.cls) - return getattr(arg_type, "__name__", str(arg_type)) - - try: - Values.__annotations__ = annotations - except AttributeError: - Values.__dict__.__annotations__ = annotations - - suffix = "_".join([f"{name}_{arg_type_name(arg_type)}" for name, arg_type in annotations.items()]) - - return cache.get_struct(Values, suffix=suffix) - - def _get_trial_arg(): pass @@ -812,7 +772,7 @@ def _generate_integrate_kernel( ) FieldStruct = _gen_field_struct(field_args) - ValueStruct = _gen_value_struct(value_args) + ValueStruct = cache.get_argument_struct(value_args) # Check if kernel exist in cache kernel_suffix = f"_itg_{wp.types.type_typestr(output_dtype)}{wp.types.type_typestr(accumulate_dtype)}_{domain.name}_{FieldStruct.key}" @@ -949,7 +909,7 @@ def _launch_integrate_kernel( for k, v in fields.items(): setattr(field_arg_values, k, v.eval_arg_value(device=device)) - value_struct_values = _populate_value_struct(ValueStruct, values, integrand_name=integrand.name) + value_struct_values = cache.populate_argument_struct(ValueStruct, values, func_name=integrand.name) # Constant form if test is None and trial is None: @@ -1401,7 +1361,7 @@ def interpolate_to_field_kernel_fn( return interpolate_to_field_kernel_fn -def get_interpolate_to_array_kernel( +def get_interpolate_at_quadrature_kernel( integrand_func: wp.Function, domain: GeometryDomain, quadrature: Quadrature, @@ -1409,13 +1369,13 @@ def get_interpolate_to_array_kernel( ValueStruct: wp.codegen.Struct, value_type: type, ): - def interpolate_to_array_kernel_fn( + def interpolate_at_quadrature_nonvalued_kernel_fn( qp_arg: quadrature.Arg, domain_arg: quadrature.domain.ElementArg, domain_index_arg: quadrature.domain.ElementIndexArg, fields: FieldStruct, values: ValueStruct, - result: wp.array(dtype=value_type), + result: wp.array(dtype=float), ): domain_element_index = wp.tid() element_index = domain.element_index(domain_index_arg, domain_element_index) @@ -1430,25 +1390,15 @@ def interpolate_to_array_kernel_fn( qp_weight = quadrature.point_weight(domain_arg, qp_arg, domain_element_index, element_index, k) sample = Sample(element_index, coords, qp_index, qp_weight, test_dof_index, trial_dof_index) + integrand_func(sample, fields, values) - result[qp_index] = integrand_func(sample, fields, values) - - return interpolate_to_array_kernel_fn - - -def get_interpolate_nonvalued_kernel( - integrand_func: wp.Function, - domain: GeometryDomain, - quadrature: Quadrature, - FieldStruct: wp.codegen.Struct, - ValueStruct: wp.codegen.Struct, -): - def interpolate_nonvalued_kernel_fn( + def interpolate_at_quadrature_kernel_fn( qp_arg: quadrature.Arg, domain_arg: quadrature.domain.ElementArg, domain_index_arg: quadrature.domain.ElementIndexArg, fields: FieldStruct, values: ValueStruct, + result: wp.array(dtype=value_type), ): domain_element_index = wp.tid() element_index = domain.element_index(domain_index_arg, domain_element_index) @@ -1463,9 +1413,56 @@ def interpolate_nonvalued_kernel_fn( qp_weight = quadrature.point_weight(domain_arg, qp_arg, domain_element_index, element_index, k) sample = Sample(element_index, coords, qp_index, qp_weight, test_dof_index, trial_dof_index) - integrand_func(sample, fields, values) + result[qp_index] = integrand_func(sample, fields, values) - return interpolate_nonvalued_kernel_fn + return interpolate_at_quadrature_nonvalued_kernel_fn if value_type is None else interpolate_at_quadrature_kernel_fn + + +def get_interpolate_free_kernel( + integrand_func: wp.Function, + domain: GeometryDomain, + FieldStruct: wp.codegen.Struct, + ValueStruct: wp.codegen.Struct, + value_type: type, +): + def interpolate_free_nonvalued_kernel_fn( + dim: int, + domain_arg: domain.ElementArg, + fields: FieldStruct, + values: ValueStruct, + result: wp.array(dtype=float), + ): + qp_index = wp.tid() + qp_weight = 1.0 / float(dim) + element_index = NULL_ELEMENT_INDEX + coords = Coords(OUTSIDE) + + test_dof_index = NULL_DOF_INDEX + trial_dof_index = NULL_DOF_INDEX + + sample = Sample(element_index, coords, qp_index, qp_weight, test_dof_index, trial_dof_index) + integrand_func(sample, fields, values) + + def interpolate_free_kernel_fn( + dim: int, + domain_arg: domain.ElementArg, + fields: FieldStruct, + values: ValueStruct, + result: wp.array(dtype=value_type), + ): + qp_index = wp.tid() + qp_weight = 1.0 / float(dim) + element_index = NULL_ELEMENT_INDEX + coords = Coords(OUTSIDE) + + test_dof_index = NULL_DOF_INDEX + trial_dof_index = NULL_DOF_INDEX + + sample = Sample(element_index, coords, qp_index, qp_weight, test_dof_index, trial_dof_index) + + result[qp_index] = integrand_func(sample, fields, values) + + return interpolate_free_nonvalued_kernel_fn if value_type is None else interpolate_free_kernel_fn def _generate_interpolate_kernel( @@ -1493,17 +1490,20 @@ def _generate_interpolate_kernel( _register_integrand_field_wrappers(integrand_func, fields) FieldStruct = _gen_field_struct(field_args) - ValueStruct = _gen_value_struct(value_args) + ValueStruct = cache.get_argument_struct(value_args) # Check if kernel exist in cache if isinstance(dest, FieldRestriction): kernel_suffix = ( f"_itp_{FieldStruct.key}_{dest.domain.name}_{dest.space_restriction.space_partition.name}_{dest.space.name}" ) - elif wp.types.is_array(dest): - kernel_suffix = f"_itp_{FieldStruct.key}_{quadrature.name}_{wp.types.type_repr(dest.dtype)}" else: - kernel_suffix = f"_itp_{FieldStruct.key}_{quadrature.name}" + dest_dtype = dest.dtype if dest else None + type_str = wp.types.get_type_code(dest_dtype) if dest_dtype else "" + if quadrature is None: + kernel_suffix = f"_itp_{FieldStruct.key}_{type_str}" + else: + kernel_suffix = f"_itp_{FieldStruct.key}_{quadrature.name}_{type_str}" kernel = cache.get_integrand_kernel( integrand=integrand, @@ -1547,20 +1547,20 @@ def _generate_interpolate_kernel( FieldStruct=FieldStruct, ValueStruct=ValueStruct, ) - elif wp.types.is_array(dest): - interpolate_kernel_fn = get_interpolate_to_array_kernel( + elif quadrature is not None: + interpolate_kernel_fn = get_interpolate_at_quadrature_kernel( integrand_func, domain=domain, quadrature=quadrature, - value_type=dest.dtype, + value_type=dest_dtype, FieldStruct=FieldStruct, ValueStruct=ValueStruct, ) else: - interpolate_kernel_fn = get_interpolate_nonvalued_kernel( + interpolate_kernel_fn = get_interpolate_free_kernel( integrand_func, domain=domain, - quadrature=quadrature, + value_type=dest_dtype, FieldStruct=FieldStruct, ValueStruct=ValueStruct, ) @@ -1592,6 +1592,7 @@ def _launch_interpolate_kernel( domain: GeometryDomain, dest: Optional[Union[FieldRestriction, wp.array]], quadrature: Optional[Quadrature], + dim: int, fields: Dict[str, FieldLike], values: Dict[str, Any], device, @@ -1604,7 +1605,7 @@ def _launch_interpolate_kernel( for k, v in fields.items(): setattr(field_arg_values, k, v.eval_arg_value(device=device)) - value_struct_values = _populate_value_struct(ValueStruct, values, integrand_name=integrand.name) + value_struct_values = cache.populate_argument_struct(ValueStruct, values, func_name=integrand.name) if isinstance(dest, FieldRestriction): dest_node_arg = dest.space_restriction.node_arg(device=device) @@ -1623,7 +1624,7 @@ def _launch_interpolate_kernel( ], device=device, ) - elif wp.types.is_array(dest): + elif quadrature is not None: qp_arg = quadrature.arg_value(device) wp.launch( kernel=kernel, @@ -1632,19 +1633,25 @@ def _launch_interpolate_kernel( device=device, ) else: - qp_arg = quadrature.arg_value(device) wp.launch( kernel=kernel, - dim=domain.element_count(), - inputs=[qp_arg, elt_arg, elt_index_arg, field_arg_values, value_struct_values], + dim=dim, + inputs=[dim, elt_arg, field_arg_values, value_struct_values, dest], device=device, ) +@integrand +def _identity_field(field: Field, s: Sample): + return field(s) + + def interpolate( - integrand: Integrand, + integrand: Union[Integrand, FieldLike], dest: Optional[Union[DiscreteField, FieldRestriction, wp.array]] = None, quadrature: Optional[Quadrature] = None, + dim: int = 0, + domain: Optional[Domain] = None, fields: Optional[Dict[str, FieldLike]] = None, values: Optional[Dict[str, Any]] = None, device=None, @@ -1654,18 +1661,26 @@ def interpolate( Interpolates a function at a finite set of sample points and optionally assigns the result to a discrete field or a raw warp array. Args: - integrand: Function to be interpolated, must have :func:`integrand` decorator + integrand: Function to be interpolated: either a function with :func:`warp.fem.integrand` decorator or a field dest: Where to store the interpolation result. Can be either - a :class:`DiscreteField`, or restriction of a discrete field to a domain (from :func:`make_restriction`). In this case, interpolation will be performed at each node. - - a normal warp array. In this case, the `quadrature` argument defining the interpolation locations must be provided and the result of the `integrand` at each quadrature point will be assigned to the array. - - ``None``. In this case, the `quadrature` argument must also be provided and the `integrand` function is responsible for dealing with the interpolation result. + - a normal warp ``array``, or ``None``. In this case, the interpolation samples will determined by the `quadrature` or `dim` arguments, in that order. quadrature: Quadrature formula defining the interpolation samples if `dest` is not a discrete field or field restriction. + dim: Number of interpolation samples if `dest` is not a discrete field or restriction and `quadrature` is ``None``. + In this case, the ``Sample`` passed to the `integrand` will be invalid, but the sample point index ``s.qp_index`` can be used to define custom interpolation logic. + domain: Interpolation domain, only used if `dest` is not a field restriction and `quadrature` is ``None`` fields: Discrete fields to be passed to the integrand. Keys in the dictionary must match integrand parameters names. values: Additional variable values to be passed to the integrand, can be of any type accepted by warp kernel launches. Keys in the dictionary must match integrand parameter names. device: Device on which to perform the interpolation kernel_options: Overloaded options to be passed to the kernel builder (e.g, ``{"enable_backward": True}``) """ + + if isinstance(integrand, FieldLike): + fields = {"field": integrand} + values = {} + integrand = _identity_field + if fields is None: fields = {} @@ -1683,14 +1698,11 @@ def interpolate( raise ValueError("Test or Trial fields should not be used for interpolation") if isinstance(dest, DiscreteField): - dest = make_restriction(dest) + dest = make_restriction(dest, domain=domain) if isinstance(dest, FieldRestriction): domain = dest.domain - else: - if quadrature is None: - raise ValueError("When not interpolating to a field, a quadrature formula must be provided") - + elif quadrature is not None: domain = quadrature.domain kernel, FieldStruct, ValueStruct = _generate_interpolate_kernel( @@ -1710,6 +1722,7 @@ def interpolate( domain=domain, dest=dest, quadrature=quadrature, + dim=dim, fields=fields, values=values, device=device, diff --git a/warp/fem/space/basis_space.py b/warp/fem/space/basis_space.py index 7275b5c68..8810f2609 100644 --- a/warp/fem/space/basis_space.py +++ b/warp/fem/space/basis_space.py @@ -1,5 +1,7 @@ from typing import Optional +import numpy as np + import warp as wp from warp.fem import cache from warp.fem.geometry import Geometry @@ -136,6 +138,8 @@ def __init__(self, topology: SpaceTopology, shape: ShapeFunction): self.node_tets = self._node_tets if hasattr(shape, "element_node_hexes"): self.node_hexes = self._node_hexes + if hasattr(shape, "element_vtk_cells"): + self.vtk_cells = self._vtk_cells if hasattr(topology, "node_grid"): self.node_grid = topology.node_grid @@ -244,6 +248,16 @@ def _node_hexes(self): hex_indices = element_node_indices[:, element_hexes].reshape(-1, 8) return hex_indices + def _vtk_cells(self): + element_node_indices = self._topology.element_node_indices().numpy() + element_vtk_cells, element_vtk_cell_types = self._shape.element_vtk_cells() + + idx_per_cell = element_vtk_cells.shape[1] + cell_indices = element_node_indices[:, element_vtk_cells].reshape(-1, idx_per_cell) + cells = np.hstack((np.full((cell_indices.shape[0], 1), idx_per_cell), cell_indices)) + + return cells.flatten(), np.tile(element_vtk_cell_types, element_node_indices.shape[0]) + class TraceBasisSpace(BasisSpace): """Auto-generated trace space evaluating the cell-defined basis on the geometry sides""" diff --git a/warp/fem/space/collocated_function_space.py b/warp/fem/space/collocated_function_space.py index 250a5d675..486e58a22 100644 --- a/warp/fem/space/collocated_function_space.py +++ b/warp/fem/space/collocated_function_space.py @@ -44,6 +44,8 @@ def __init__(self, basis: BasisSpace, dtype: type = float, dof_mapper: DofMapper self.node_tets = basis.node_tets if hasattr(basis, "node_hexes"): self.node_hexes = basis.node_hexes + if hasattr(basis, "vtk_cells"): + self.vtk_cells = basis.vtk_cells def space_arg_value(self, device): return self._basis.basis_arg_value(device) diff --git a/warp/fem/space/function_space.py b/warp/fem/space/function_space.py index 91d8fcf96..a970621a7 100644 --- a/warp/fem/space/function_space.py +++ b/warp/fem/space/function_space.py @@ -1,6 +1,6 @@ import warp as wp from warp.fem.geometry import Geometry -from warp.fem.types import Coords, DofIndex, ElementIndex +from warp.fem.types import Coords, DofIndex, ElementIndex, ElementKind from .topology import SpaceTopology @@ -47,6 +47,11 @@ def geometry(self) -> Geometry: """Underlying geometry""" return self.topology.geometry + @property + def element_kind(self) -> ElementKind: + """Kind of element the function space is expressed over""" + return ElementKind.CELL if self.dimension == self.geometry.dimension else ElementKind.SIDE + @property def dimension(self) -> int: """Function space embedding dimension""" @@ -71,7 +76,7 @@ def trace(self) -> "FunctionSpace": def make_field(self, space_partition=None): """Creates a zero-initialized discrete field over the function space holding values for all degrees of freedom of nodes in a space partition - space_arg: + Args: space_partition: If provided, the subset of nodes to consider See also: :func:`make_space_partition` diff --git a/warp/fem/space/shape/cube_shape_function.py b/warp/fem/space/shape/cube_shape_function.py index d064ca031..7370c81d2 100644 --- a/warp/fem/space/shape/cube_shape_function.py +++ b/warp/fem/space/shape/cube_shape_function.py @@ -322,6 +322,87 @@ def element_node_tets(self): return grid_to_tets(self.ORDER, self.ORDER, self.ORDER) + def element_vtk_cells(self): + n = self.ORDER + 1 + + # vertices + cells = [ + [ + [0, 0, 0], + [n - 1, 0, 0], + [n - 1, n - 1, 0], + [0, n - 1, 0], + [0, 0, n - 1], + [n - 1, 0, n - 1], + [n - 1, n - 1, n - 1], + [0, n - 1, n - 1], + ] + ] + + if self.ORDER == 1: + cell_type = 12 # vtk_hexahedron + else: + middle = np.arange(1, n - 1) + front = np.zeros(n - 2, dtype=int) + back = np.full(n - 2, n - 1) + + # edges + cells.append(np.column_stack((middle, front, front))) + cells.append(np.column_stack((back, middle, front))) + cells.append(np.column_stack((middle, back, front))) + cells.append(np.column_stack((front, middle, front))) + + cells.append(np.column_stack((middle, front, back))) + cells.append(np.column_stack((back, middle, back))) + cells.append(np.column_stack((middle, back, back))) + cells.append(np.column_stack((front, middle, back))) + + cells.append(np.column_stack((front, front, middle))) + cells.append(np.column_stack((back, front, middle))) + cells.append(np.column_stack((back, back, middle))) + cells.append(np.column_stack((front, back, middle))) + + # faces + + face = np.meshgrid(middle, middle) + front = np.zeros((n - 2) ** 2, dtype=int) + back = np.full((n - 2) ** 2, n - 1) + + # YZ + cells.append( + np.column_stack((front, face[0].flatten(), face[1].flatten())), + ) + cells.append( + np.column_stack((back, face[0].flatten(), face[1].flatten())), + ) + # XZ + cells.append( + np.column_stack((face[0].flatten(), front, face[1].flatten())), + ) + cells.append( + np.column_stack((face[0].flatten(), back, face[1].flatten())), + ) + # XY + cells.append( + np.column_stack((face[0].flatten(), face[1].flatten(), front)), + ) + cells.append( + np.column_stack((face[0].flatten(), face[1].flatten(), back)), + ) + + # interior + interior = np.meshgrid(middle, middle, middle) + cells.append( + np.column_stack((interior[0].flatten(), interior[1].flatten(), interior[2].flatten())), + ) + + cell_type = 72 # vtk_lagrange_hexahedron + + cells = np.concatenate(cells) + cell_indices = cells[:, 0] * n * n + cells[:, 1] * n + cells[:, 2] + + return cell_indices[np.newaxis, :], np.array([cell_type], dtype=np.int8) + class CubeSerendipityShapeFunctions: """ @@ -632,6 +713,12 @@ def element_node_tets(self): raise NotImplementedError() + def element_vtk_cells(self): + tets = np.array(self.element_node_tets()) + cell_type = 10 # VTK_TETRA + + return tets, np.full(tets.shape[0], cell_type, dtype=np.int8) + class CubeNonConformingPolynomialShapeFunctions: # embeds the largest regular tet centered at (0.5, 0.5, 0.5) into the reference cube @@ -656,6 +743,7 @@ def __init__(self, degree: int): self.NODES_PER_ELEMENT = self._tet_shape.NODES_PER_ELEMENT self.element_node_tets = self._tet_shape.element_node_tets + self.element_vtk_cells = self._tet_shape.element_vtk_cells @property def name(self) -> str: diff --git a/warp/fem/space/shape/shape_function.py b/warp/fem/space/shape/shape_function.py index c53b19703..0ec8a29e4 100644 --- a/warp/fem/space/shape/shape_function.py +++ b/warp/fem/space/shape/shape_function.py @@ -1,3 +1,5 @@ +import numpy as np + import warp as wp from warp.fem import cache from warp.fem.geometry import Element @@ -100,3 +102,8 @@ def element_inner_weight_gradient( return grad_type(0.0) return element_inner_weight_gradient + + def element_vtk_cells(self): + cell_type = 1 # VTK_VERTEX + + return np.zeros((1, 1), dtype=int), np.full(1, cell_type, dtype=np.int8) diff --git a/warp/fem/space/shape/square_shape_function.py b/warp/fem/space/shape/square_shape_function.py index f34ef5367..9b9e62f28 100644 --- a/warp/fem/space/shape/square_shape_function.py +++ b/warp/fem/space/shape/square_shape_function.py @@ -206,6 +206,31 @@ def element_node_triangulation(self): return grid_to_tris(self.ORDER, self.ORDER) + def element_vtk_cells(self): + n = self.ORDER + 1 + + # vertices + cells = [[0, (n - 1) * n, n * n - 1, n - 1]] + + if self.ORDER == 1: + cell_type = 9 # VTK_QUAD + else: + middle = np.arange(1, n - 1) + + # edges + cells.append(middle * n) + cells.append(middle + (n - 1) * n) + cells.append(middle * n + n - 1) + cells.append(middle) + + # faces + interior = np.broadcast_to(middle, (n - 2, n - 2)) + cells.append((interior * n + interior.transpose()).flatten()) + + cell_type = 70 # VTK_LAGRANGE_QUADRILATERAL + + return np.concatenate(cells)[np.newaxis, :], np.array([cell_type], dtype=np.int8) + class SquareSerendipityShapeFunctions: """ @@ -501,6 +526,12 @@ def element_node_triangulation(self): return element_triangles + def element_vtk_cells(self): + tris = np.array(self.element_node_triangulation()) + cell_type = 5 # VTK_TRIANGLE + + return tris, np.full(tris.shape[0], cell_type, dtype=np.int8) + class SquareNonConformingPolynomialShapeFunctions: # embeds the largest equilateral triangle centered at (0.5, 0.5) into the reference square @@ -516,6 +547,7 @@ def __init__(self, degree: int): self.NODES_PER_ELEMENT = self._tri_shape.NODES_PER_ELEMENT self.element_node_triangulation = self._tri_shape.element_node_triangulation + self.element_vtk_cells = self._tri_shape.element_vtk_cells @property def name(self) -> str: diff --git a/warp/fem/space/shape/tet_shape_function.py b/warp/fem/space/shape/tet_shape_function.py index 607f6326a..32067f4a1 100644 --- a/warp/fem/space/shape/tet_shape_function.py +++ b/warp/fem/space/shape/tet_shape_function.py @@ -442,6 +442,14 @@ def element_node_tets(self): return np.array(element_tets) + def element_vtk_cells(self): + cells = np.arange(self.NODES_PER_ELEMENT) + if self.ORDER == 1: + cell_type = 10 # VTK_TETRA + else: + cell_type = 71 # VTK_LAGRANGE_TETRAHEDRON + return cells[np.newaxis, :], np.array([cell_type], dtype=np.int8) + class TetrahedronNonConformingPolynomialShapeFunctions: def __init__(self, degree: int): @@ -450,6 +458,7 @@ def __init__(self, degree: int): self.NODES_PER_ELEMENT = self._tet_shape.NODES_PER_ELEMENT self.element_node_tets = self._tet_shape.element_node_tets + self.element_vtk_cells = self._tet_shape.element_vtk_cells if self.ORDER == 1: self._TET_SCALE = 0.4472135955 # so v at 0.5854101966249680 (order 2) diff --git a/warp/fem/space/shape/triangle_shape_function.py b/warp/fem/space/shape/triangle_shape_function.py index eb3871c99..f28799612 100644 --- a/warp/fem/space/shape/triangle_shape_function.py +++ b/warp/fem/space/shape/triangle_shape_function.py @@ -309,6 +309,14 @@ def element_node_triangulation(self): return np.array(element_triangles) + def element_vtk_cells(self): + cells = np.arange(self.NODES_PER_ELEMENT) + if self.ORDER == 1: + cell_type = 5 # VTK_TRIANGLE + else: + cell_type = 69 # VTK_LAGRANGE_TRIANGLE + return cells[np.newaxis, :], np.array([cell_type], dtype=np.int8) + class Triangle2DNonConformingPolynomialShapeFunctions: def __init__(self, degree: int): @@ -317,6 +325,7 @@ def __init__(self, degree: int): self.NODES_PER_ELEMENT = self._tri_shape.NODES_PER_ELEMENT self.element_node_triangulation = self._tri_shape.element_node_triangulation + self.element_vtk_cells = self._tri_shape.element_vtk_cells # Coordinates (a, b, b) of embedded triangle if self.ORDER == 1: diff --git a/warp/fem/types.py b/warp/fem/types.py index 0a2b39ef8..b54c3f833 100644 --- a/warp/fem/types.py +++ b/warp/fem/types.py @@ -1,3 +1,5 @@ +from enum import Enum + import warp as wp # kept to avoid breaking existing example code, no longer used internally @@ -31,6 +33,11 @@ def get_node_coord(dof_idx: DofIndex): return dof_idx[1] +class ElementKind(Enum): + CELL = 0 + SIDE = 1 + + @wp.struct class NodeElementIndex: domain_element_index: ElementIndex diff --git a/warp/tests/test_examples.py b/warp/tests/test_examples.py index 225d4f9c2..74c2056fd 100644 --- a/warp/tests/test_examples.py +++ b/warp/tests/test_examples.py @@ -412,6 +412,18 @@ class TestFemDiffusionExamples(unittest.TestCase): devices=test_devices, test_options={"headless": True}, ) +add_example_test( + TestFemExamples, + name="fem.example_magnetostatics", + devices=test_devices, + test_options={"headless": True, "resolution": 16}, +) +add_example_test( + TestFemExamples, + name="fem.example_nonconforming_contact", + devices=test_devices, + test_options={"headless": True, "resolution": 16, "num_steps": 2}, +) if __name__ == "__main__": # force rebuild of all kernels diff --git a/warp/tests/test_fem.py b/warp/tests/test_fem.py index e0b8fb0db..ae0f76d7d 100644 --- a/warp/tests/test_fem.py +++ b/warp/tests/test_fem.py @@ -1305,6 +1305,93 @@ def test_particle_quadratures(test, device): assert_np_equal(measures.grad.numpy(), np.full(3, 4.0)) # == 1.0 / cell_area +@wp.func +def aniso_bicubic_fn(x: wp.vec2, scale: wp.vec2): + return wp.pow(x[0] * scale[0], 3.0) * wp.pow(x[1] * scale[1], 3.0) + + +@wp.func +def aniso_bicubic_grad(x: wp.vec2, scale: wp.vec2): + return wp.vec2( + 3.0 * scale[0] * wp.pow(x[0] * scale[0], 2.0) * wp.pow(x[1] * scale[1], 3.0), + 3.0 * scale[1] * wp.pow(x[0] * scale[0], 3.0) * wp.pow(x[1] * scale[1], 2.0), + ) + + +def test_implicit_fields(test, device): + geo = fem.Grid2D(res=wp.vec2i(2)) + domain = fem.Cells(geo) + boundary = fem.BoundarySides(geo) + + space = fem.make_polynomial_space(geo) + vec_space = fem.make_polynomial_space(geo, dtype=wp.vec2) + discrete_field = fem.make_discrete_field(space) + discrete_vec_field = fem.make_discrete_field(vec_space) + + # Uniform + + uniform = fem.UniformField(domain, 5.0) + fem.interpolate(uniform, dest=discrete_field) + assert_np_equal(discrete_field.dof_values.numpy(), np.full(9, 5.0)) + + fem.interpolate(grad_field, fields={"p": uniform}, dest=discrete_vec_field) + assert_np_equal(discrete_vec_field.dof_values.numpy(), np.zeros((9, 2))) + + uniform.value = 2.0 + fem.interpolate(uniform.trace(), dest=fem.make_restriction(discrete_field, domain=boundary)) + assert_np_equal(discrete_field.dof_values.numpy(), np.array([2.0] * 4 + [5.0] + [2.0] * 4)) + + # Implicit + + implicit = fem.ImplicitField( + domain, func=aniso_bicubic_fn, values={"scale": wp.vec2(2.0, 4.0)}, grad_func=aniso_bicubic_grad + ) + fem.interpolate(implicit, dest=discrete_field) + assert_np_equal( + discrete_field.dof_values.numpy(), + np.array([0.0, 0.0, 0.0, 0.0, 2.0**3, 4.0**3, 0.0, 2.0**3 * 2.0**3, 4.0**3 * 2.0**3]), + ) + + fem.interpolate(grad_field, fields={"p": implicit}, dest=discrete_vec_field) + assert_np_equal(discrete_vec_field.dof_values.numpy()[0], np.zeros(2)) + assert_np_equal(discrete_vec_field.dof_values.numpy()[-1], np.full(2, (2.0**9.0 * 3.0))) + + implicit.values.scale = wp.vec2(-2.0, -2.0) + fem.interpolate(implicit.trace(), dest=fem.make_restriction(discrete_field, domain=boundary)) + assert_np_equal( + discrete_field.dof_values.numpy(), + np.array([0.0, 0.0, 0.0, 0.0, 2.0**3, 2.0**3, 0.0, 2.0**3, 4.0**3]), + ) + + # Nonconforming + + geo2 = fem.Grid2D(res=wp.vec2i(1), bounds_lo=wp.vec2(0.25, 0.25), bounds_hi=wp.vec2(2.0, 2.0)) + domain2 = fem.Cells(geo2) + boundary2 = fem.BoundarySides(geo2) + space2 = fem.make_polynomial_space(geo2) + vec_space2 = fem.make_polynomial_space(geo2, dtype=wp.vec2) + discrete_field2 = fem.make_discrete_field(space2) + discrete_vec_field2 = fem.make_discrete_field(vec_space2) + + nonconforming = fem.NonconformingField(domain2, discrete_field, background=5.0) + fem.interpolate( + nonconforming, + dest=discrete_field2, + ) + assert_np_equal(discrete_field2.dof_values.numpy(), np.array([2.0] + [5.0] * 3)) + + fem.interpolate(grad_field, fields={"p": nonconforming}, dest=discrete_vec_field2) + assert_np_equal(discrete_vec_field2.dof_values.numpy()[0], np.full(2, 8.0)) + assert_np_equal(discrete_vec_field2.dof_values.numpy()[-1], np.zeros(2)) + + discrete_field2.dof_values.zero_() + fem.interpolate( + nonconforming.trace(), + dest=fem.make_restriction(discrete_field2, domain=boundary2), + ) + assert_np_equal(discrete_field2.dof_values.numpy(), np.array([2.0] + [5.0] * 3)) + + @wp.kernel def test_qr_eigenvalues(): tol = 1.0e-6 @@ -1398,6 +1485,7 @@ class TestFem(unittest.TestCase): add_function_test(TestFem, "test_dof_mapper", test_dof_mapper) add_function_test(TestFem, "test_point_basis", test_point_basis) add_function_test(TestFem, "test_particle_quadratures", test_particle_quadratures) +add_function_test(TestFem, "test_implicit_fields", test_implicit_fields) add_kernel_test(TestFem, test_qr_eigenvalues, dim=1, devices=devices) add_kernel_test(TestFem, test_qr_inverse, dim=100, devices=devices)