From 8d383e31861c09769cc0a3edbedaab95b4e9de2b Mon Sep 17 00:00:00 2001 From: jtgrasb <87095491+jtgrasb@users.noreply.github.com> Date: Thu, 26 Oct 2023 18:58:53 -0400 Subject: [PATCH] Capytaine v2 Hydrostatics (#284) --- examples/data/bem.nc | Bin 36896 -> 39290 bytes examples/tutorial_1_WaveBot.ipynb | 37 +---- examples/tutorial_2_AquaHarmonics.ipynb | 67 ++++---- examples/tutorial_3_LUPA.ipynb | 56 +++++-- examples/tutorial_4_Pioneer.ipynb | 64 +++----- tests/test_core.py | 17 +- tests/test_hydrostatics.py | 204 ------------------------ tests/test_integration.py | 49 +++--- wecopttool/__init__.py | 1 - wecopttool/core.py | 163 ++++++++----------- wecopttool/hydrostatics.py | 169 -------------------- 11 files changed, 194 insertions(+), 633 deletions(-) delete mode 100644 tests/test_hydrostatics.py delete mode 100644 wecopttool/hydrostatics.py diff --git a/examples/data/bem.nc b/examples/data/bem.nc index 33d708163f9101d70e137f143dbee0f64e2ccbf7..bca7066d50a889a512b2855b0431960bc0dd06de 100644 GIT binary patch literal 39290 zcmeFa2|QL^*D!vW35gIz$&@J}Q;NL}nL?R|GG)j-D-AMM8b#(phE&p|D8fEeNFgMm zq(X*@%uW2yaa~u>efQk=`#kqMe82DgfA{)v?X%Zed#|12!?<%Xzq`iiakI%c%&4z#qw3rIR{$~^%u@`dGMa(GvP94?|I(GXWW z0wj7opkpUx(30=S(2+wkN=Ra)+#sWW`DGL$k(hqL{d1IY^}jBVXibP_?gClf@ z&XVK_#Cj5m1rH|f-p(Y&;TNQhM0Ju#gQi@~r27;Ep~ZTdo)*$L!3B^OD&3X+`v8fT zpuo#QfKX@S_f69AYUJ48zyG9xh5Bdtoz8w&_rI%v-zooh2Jkxr`kf{G&I*24S--37 z-`U3R?BaKpx722P>OYm?Z^VCI1LU6hEDg;8oViO|pba=zJkibOKc=G~1(TEZB{=ROzoENgC z!>JbI(&5=srPJ~B^>qEcTzC=!*zop;9#6oE7ltT~uBVTWm9tN%+aD_O&q_lg%|_sb zxTG9Pt|qzPR1;YfRE^-OTc{>}Tq{ITiGQ+Qm`M!0t1aR_^L*hPBT_n;uE;1&D z9Q<^8)?sz=xhs8s(?&SYlbw~7I zWZkVItP6)&a`01ktpAg`BZ~Uhbtn19y5qnzBQo^#bH@$T&&AD^qCBY!bX;WRN&ZH8 z)W?k){`K~OM9QNmZTuId{f`1;a`7R7aFmJr!{$@sV0$3$yxriIL7};izAwwzpgytKejlc!yq)_;_n;a;}(Wr zm5RH5U2uUJh#Q1E=qf?-#seIdEVO3C)n*&z+a%I!tx5@9bMwt4suxd#JH)(1?UEu+ za(LuCKPE)qymq1(G^IuZ}Lg%1vK4hnJfcXuS1hlb#d1zv50a`9o?f=~)_3hUR( z%dJ&VvXGZoQ(TY#$trD7P};Dbj>JLEV^@%yo1de5s2{n~iPS`cKom*G&Bw#v*VoNC zl-yh3brs^gGt_5+kE37E?|X|s=*0*x`VW)t` zBETPy{p;@Z+d^&<_3|?eg?lXv7F@FYxlgIXMS8 z2K)Q$b_-fyNjz@^OP|o+7nkr;7LvFGxsjEM86oKmX82s$H$ueL@Wqec_m`FSEK{^4sh&!X&C-hX{_-s zRp@7z%C?f2QTXZ6hvMxJ+4iYV2=4NcNJ-y48L1_vC-G9oF8`?#6FnV1*Kp*@9&$JC zMMI-VU})+7F*YL8{KNUC@FT_s3qTG3Grm!}9666}N?z2l6T_c(?*;_vUm80t^p*>W zC|vM>u=oX&(;<(c{_A5GHj2mjhsI0oF*RLX-CTuyor8l(B%x=k=N7v5YhGlzXz+2^ zS3EvA&GdEv(LHOmqa4+hBLQ&3PYyggf}lM?R*y<_KguUpYa!4Cqz*h_#z`~&6w z82G_z7rTdVoI*DV?zg}8crQ3jp`QiEvfYfWMZQwfeFxh-N95}y<|ye&@I<=m@fFKx zAUPje=yz~`cpE<*B=e<*rY=?f#|0Q5IbTMoq))R`dleHT%gGG+*6LV%Tg?K=d=|m) zoJVwTiLyenoNTaX`jM*l8g@wLyBOAVvq_6caX@l;mq7F0i>G{3mQdv2gb&Of9eI?> zNs)&ODnA^K6ier#$h8z+RAJi}m9dl}A2;MLTpRrAEH_0?9(YV(tIvi^9*VrY@cXF% zx8y8dirmYfA5Ege=j>$^`T5}M;9J(lIee6I@IyWAtqW6#R9B6tvflzRoF;%Wt4T5s z4cg^~$EkDIRQxv25El;|Pk4v~5S}0b!2x*3?5OXCHw_DMVt$YSVm6WhcfxJ>4FO?X zzXf89LjaWzl|RwJ6F|jN`4hed0e{AyGGP8w{zT{eTUeNr#1X{vMS!O-Ib22#e&U13 z#UoLM=s)$)Kv5Jwg`TP(DxPp&383NwDGH)a^-}RvdI2)|KhfiB6AM74ze-Ug-U$B= zPt9l5pWual!Bgc~{R^I&kMJ*eD!s@rc&a?26g)K$r(6)7G66=G+cX_LaX&mL*v%D> zuUWwULqDr1gc5&(#~ZE%pz47tr_?WaYCgDU@rU%d1MxfjFvXzoe(-mAYCf{R;HiAq z{eq|FL+LAi)+1Fv^1ski^QU%H)N)br>wlrA@?Gfp@O!_~Q}a z#L$x%=1}8=&rFQ3sc~X1;J0yNX#3mr#Mt_`>4~xUZ|RA#{BPsbF+7f57)Mk25Myg% zyiJW0WAIgnMojN!i9r04M1t{<7zR@BQ^$}b(pUs}dxAbBhpO;OApvrfxCA7UlSrZz zybu&INE9_l6rkMagQwo9;(xqwfsDir*Od?E9+%~!&@X`|A0FU;iz(^YDELK`^vv+^ z5dxSf`7pxHZ4TOd*%>JL(ZdIFZ1?5Z=qUNpLf+!jhh$u7DEW{mbR-I&Il$Gm#baI0 z4?vdlJ6Pddc9ApX8z7hKD`?vA{#pNQSuuD+fJtXaV_(wmW!;`Klq({7B(dQ+3jD_ z?Q~iumpB#sFX^mV6wfWmjevqZ-YxU*Dt<>FSt(PQ%+dU$Jb&{~eE8tCdoFXPhO{DXK*hiId)CSW|;M`F=c;N*H8NKeiS!r6OH{@y1(&%Yv6AU{H=k% zH9*z?u?EPzFgZ_0BSoFmq)eK#X_J?K+3?v+dSdg*KeiSU)Ahs@G~q7n658qSx_4on zS7;Bu@i2(I>p>{U!(Rwrd-wErW%D7g$#USU9r*r^e{3lo_wiZT*sBX$j2MYR5l$q&#mJXT^N%e? z>n?3OKy2dpcXUK8J2j;0h^jNT-DW*&uAhnewBH8?emhRPew#AZAsGHBVpke!HWV`u z(w+b!G(#dca^$etnsxaZbDpq%T~R{_Xab?aEhEw^h5@~}t!Hp9IKId-w zz(SVz{BnC`%vEE>O4cW-Xv)|>^PBesxLnN4C%~_OwXj>{mwxqx-4D7ovJ9HQB9fYm zV$~2B4M=e-&8|lVBPH);Hib~%L3U1|cQRVUq~9}j_<*ePlCsYgh9JC`w2dl}-Bb(m@94~7bKcY-a|?J$``FTLqS zX22Vdy)*&_Q(%u+S;8wGJM37mlt*Bx7mO^Yf2dyK3r$vT`Z!Rmg4vEgKb{ZMQ9yOp zB=3?*&@k*E|D;O+%U6E(VPwh^UYoLCs?pm79(IKCXI~fs{f|oK^=TVWYp?tH^_gX; zeeZjcqEs;4yFc@xv4kBKd9Nzzy4nn&pO;Abx?u|JFc`Rfhh>yy~{QR_+BuS9E^z zZuf@6xH_fEa;`VoQMSbhl?j|60 zpln}d_z*C^o3$$T%_H=9dZw-VKpB!Rq8p!A4u;JSG(($K*kMsw-ly&=%m9neoH15w zr+}H>WckKPTg>d8tO8;GGH%)Xec$>*I-eD!nG99ztIZl7_O3KE!Z*8KUT6}8(XmM< zjjqQOj?n?q7u^1D)jD~9X#%Z9v$wc+4*~6+#pA9 z#3f{j9cJ7;ooOXF1Li))(n*ULFSO&dk5n+29-(BSpuZTPjistM2@j#Jp_GXyLWzlRQW z*Q1Mjn;4fEm7xjm39n__gW(pL^ct>NTWtS`(!=7RX)v@}zmA}7Xq|W!MBXM{F^a!!@(s(Ih)dG(+uhV(K zYb#UJM`xPA*%{e%DdizxeASN8b!9yg<~cKZi@OZ*-;`Fk?i&mjuduJ%Ib@4fUX{(0 zIW!HjN+v$+W}O119P^4b_ieE|H&o-FKlFk_K5USwpc|9XG_r1X)v`T^~^omDNtM5zUBG_TP$;!k$1~uFL=C2 zuhml37y2!_Wa7f5hDEzE_jrzEpvY&p7^T7|z_zGUOXB5}v2W*&C|_prf}YouGK@N! zfXVW4vt2?%pow{yck$ujd9QCV|Bh`vdz~@aHeIth>j<7UMn5pP&593wrIg9dxzxfexcp zrB>`}SjpCiRS&PEqc@i9{qCLPpeRi0X6lSGCi2?DS(?KOhMm}vI5*n_>|@$vUeFGK zAV(XCbhrB`xbybCu-j#*65ctcIUEdkgum3^;9`p{$&o$m_Z$uEcHDr8mt4KJ zDmw#(RW;0w(v5?`J9FmWrj@b!(*Ai9E-#qh5u&>8RTH=}!|7(oJ_PuKugjlG!q8E# zIE#G_9A@I<{_oNjxxSp4JM0b%U0QD!+w5cI4 zsHjf+wIj|8j`!%%&~&L{z1Mv#Mmo|F&87Xkq|={4_NMmRh7l^5XB$6BnbQlhXP?4A3?g?)UUA;E4L=~GpF|h4IY&z<&ZHv?6 z{|u0$)!MOB%2+;wYwInXf6h+#af?q)z==-$rU`EUVOwv7=Dxd!j*n))Kg4orK5XkBo#mTnH29slK4YPYQ2AB$Cj|(4~0-QbC zGNZHBSk%kU^MP#Mux;<5J$iN?P;6-2u}o4Kv!^|o{}iPoo;2fz&edZ;a&D*m%rzzK z0E>lBBL4inSQB2Mmed5?bpj*ZWru*-N0`oJb_W?`l$>^Zau40wP!TTwF$8{07Eo$z zu*NzQS)av8j)5<4IMw#&PXT|nyUC9Vtg$}r#GDv@Z^&mpdQ&ic7d+at$?Mg;BF1R6 zd*Vev8tQF(F~)8+3OM`AxL+?<#HQ=~CdQfY`WoxpmW=y9S$csPl81%>Nua-o2t-85P(;D+@S**8A z#v6ulEQ6*JDRy_vP%s{ZEIu zeYZaHGy{hF2If))LqK|EWcb{}TZrZTuv(jU1)98K5|yG83d@#TtzGoh3hSxrDluRD z2}~^g(mmWe1ujdh-QNAo3hUf&s&*ZE!}tz~FkV4-Xf42v7&KlF0 zVY@$qv)eFw*Y~*p$9C(35YE4nwkPpDQ!}u-iBFEU41tmDeipop#VGr+aHZwl3gmsl zdc4Cc6o&geoIiEK3e!9`U%IyY1IQJUlzKmh`=6tq`Bl8Fuu!Lu-<_ zcq=R@b=_%stcG^7(3&X~)s~L=#h(2DG`PteDS~r5i3kuGi&vynGvv7_mk}C@@de@ zbWx9y&I-#^jyx{u=M8sw&|Z!#afOLRJZd41a#-%&+Fs+QsVH)@b_&h?5g>js%Bq)B z4qF1=Y@r?XgnbA67+xwg0h$NF!l~lJAkost@dft{#8@0R8mwB0mOT5KSe69kJgEgWT+XqNRNaSr`o}g(Z z3Tt>O6L=~VKHPNFP_oAY>j;d`9eXhdnz9Bg_6AOaS;vB9amWITl`$HzZ1si_d;G@6 z&2~bQY3{?i#*!FUw{K2(W(rEaJp1}&>mXSBJ*D4~UJ@G`=!-be5OUgm!0m)G?4|4}%2e+QtYuQd3aNyY*>jePzK89?)nRbXifPyG&O$WBasSlmRh7tM2YuAIS}5$y zx*@a_=P%S;b}J|AJurXxgZK3DX|PT9yWKUMe=6hbi?jPx#6xzkH zAE`-F{ShfBmg7lNf6jaG{V6A#o%b4Sj{TmC{WDJ}Ho3Nn!>bXPs12o=Ru2Q0uX*eV zHHGM|Yq8BYScy!vEG50R`fu30?&(S4_uy0A+*gz`4fY25^Ixp7z^;Zx zhiv`s4O220JcekTA^%N2J?<7UtZC+qZK6>(5>;I+iIof)gV zP(J|7GrBDU&rbvYFp+~w-z=~;Mv=$INj~uHyve2X%T7@CT$?V>TT!f7$N#y&p=6X> z@4n=G`2eUjP>T@UD~hd>Kl3iF*%Q(h>YO=Z-3U}fF6)WU4}-5->)G{E3Xz7)NDM1e zC2Cs1DpP9}0tJtpOxhG;iFF%Ztq?pk094%5C0j4!3pC3MfP}Rr2HYYIE%5U5j=7sz z+dDyl-76|owuxdxA0K7iyqAOy35M@HVmAOhp3EnSj*DQ^2Sbe%asFQrpWmkCjet(| zRY@Jk2zYkHZRL(JfR1|#k~^+V zgG+2G_m;n~#I7rRv$y}~4ZU>lGu&C_1ZBo=uVkDO!R!kI>zA!QgG!Yd{3lrl!0rbO zjU{;^7>~Hz{`$wBu(s`1E3Hi<7}%J3<}1#ByI#RjMUN{eN>x>8*+d1(QGKuScwaE2 zzw-IbXoeLwMN?GE^`#$_%%_}3IDdVAhjF@aD{QmL%Uh1^-Z1}EoW1fJN7yXTv@y6& z1Y5_>D43pu^B32#QhwAACgr({6i$d>YATyww>Ee}PIb*6j!kU?uroSpH3sX>64l@*1_Ge zD!0TNay-5g1acjrL4@<3>>LrSviW0>`MHxQy=u@x_GUjw8?g!8qbq{lxxC@@QoQ^x zeQMsv;`Q&cHFRo{Z3MJ_*6Jy6FG4;Y$HU*&R3JAcu9Y97gJ4rrW#ng38*Dh8ee0Fd zeo)zc%E065G*CUbLa3J61_M-dl+g4Nc1Ht1tKj<~J} zU&zYu2NQi!<(8}>*opByr^TB*p^|0Eu5GXpj30@4t;IJ2*al*Tgog37C{pd^iy1Hn)%@@o}~UgNC&2$9{0-;u-6y>(d~z zQfynvBU@~(mH4(SC2v@iQ!>K0%n9D?lTx-G6U8A5qtQQSu&~znN>!wOSk#DZFJW(1y?djIUw1uGaUxdB$J>LRHZ&x18L9+zda=wQp+K93k1 zc)|VptP@YgpMcWxenXdA{b1v(CxXguZOA_<^=j0|Tj)B}?Y5N-fGJLUM8ps8z-Bf+ zy1ix-3GUe=zPO@)0x+7FE%$KWfwg`tx+9$J1*OjE`byljh8A|}o3>xz$HKE*$Hb0Z zKo``cC+n`g205EIjOXp(!#<5SWXvN^_)*;J`eU9ZFmUA5jWxJ`^ktj1^K4=sx|?Ns zSMy9I3VQO&{ab=R6c2hbW?^lQ<-E28y(YcjK;j3sz?o?)wJ&*2i4}yH{d0-^C6ed*lgUp~Rks zo<^Yk@RTXTz7dc$X&;fa@e<-WRjSjYUx|+Azst7n@PlG@3apa?4%pJ_6VJaKdJCLy zafN7iO@sYxMZ!mZ*kdP)ChdnddBFoGSIk)FY==b=N9x?~3Sk8&xwapSKZ7I(51(l5 zcnfB#hxJ!W31Q^%XYe7NSF3i^4RFE76M+ zVTBWG{b0^P0af#74%nXfWlprZ@4#DeS&~ibG)Qtvj-I{kfcd8?cNpU1`I%OJS7oss zu;xPg*Uucn*be;k-!pzbyhc*~Dfk)F< zob`i8fMYD@$sYb9B-SA${86A1F@8B*P?vhbQ_wUR z*u_y_zrzs|yiA*CkNcBb2Uc&1jMxFsYVhnMaf@IZ*O=MpXeT0WM{%sp`5m|{?Okz7 zXf-COVE%LrZ~t-+3v9h*)C5uzih_<=jewN(Ah(d|DynYf@An(7Kvi;qyYg6kq1k#r z*7J*;up4tH+}3}22cEy+mkq$%pEyy+?7&G!%(Q6zk!6*h@UxNT6|8;-R1nr3+I|J^ zXS!=&yGO^NXMG;`ymQ|H^II1Sj+L**+_kNpvKu{Nn5fdWjIKstt+6p?kJboaO)=+L zy6P%w_zG;sQZ%93FQBh_RFJ&cVM@eNq3u|FqZPAttScZ ze>Vw4T#)?H2sS#-__eB!fc+9Q-b}`L`l)&!X)JY^ ztUwB{9Man>ykYf$E7nc*PMD1LxaY?+ec%Y!k){m1{Tp_D$yQzJgw1`=@l?0P>wn?o zFB#hZQ}57G7yszDQdWN`i+}%{-=f1`@*vK<_`|p8h|N*NdK2X>RDVAqS7+Y<`~Zyq zOIJHrt5R6 zxM_zT*P7=VclOdvm@Ao9J| z_Yzdx!Kw4tZf~S504055UDu%kXvul^!A~!!%KKuJs(6MR=(rZ!R}yN57^J#OpJew` zX|CjYazo;shVoO($|C8vRU$HSAB$CjLCgA@!}`7R8r1?HUgi{IRkI-tZMq>pG{R|a z=P6~*f=@3P^8Ca6!Ok}=dsgh|s$v>{k~b3>4Gz}??tSH14f@%SU8L7i20b;sB$=xj zRV|sTmPyDO0@_!LlNTNLL$&G~+#e*os|qtL7JIEG3UnQZZ_BW;qnmkbo0s;+gV|Lq zrtHr2;45~rMqo1yYT`P}?$*u^Xg+qF&Z=XB_XR`UQ%Vk@=gzL2CZumwQnTDl2jh}K z=kc{a9yYXs)4`2ib1p`R8J=q$`*FSMYv`H09>E`=WiGO&o<{tOwvRuI7X_z|*6r806$;o+csuxLnWIf!xxtmeH-QE=6?U|2Gg@J; z^pSCwDdH1iNsQvw0rv_ljkIjmK$9-Lr3^b8P{;N_qqiDsP|j8v4grmups^t|;#0zV z5LES~3FHNDFTZE7a9|#4$>xTJ1 z-pofjH8#A_2Teiha`pM1FNHwkCY{jeu2}SV{7k;hqFT^DmCELa1)+-Fx{UcLG3ba@ z{F9Z(-NDk;M@}AIqXui6^+Alv8Z$=w;cw9@FuL0GGUzN8X8Uy+CPvewJ z@=&{zseWfNA1JZQ8-3EZ2b|L1&-``8Rm7vpz+hS91mcqob@!dd2hOGYqv|wM(Zp>= zk~L@s^=C}>(uqf-H-;ax)4QEINd^(W>=pT zo$&$KT7PYqHLW0-y;RM^`65F12hOq2Jp)e?_HAg?i${m<_`~Z}1!(N42tZ(Ccki9yrqi zAlb0BL2b+nlzw`g%UE<5?YXvZa7ceQaG-;>DiUu2H{BPL^>QVsJEnehCv7)arq)%= zdnXOCxp~}I$i9toJRUJzD>(r^%&k4SRm%h}Ewj$M6|@A~8Pyl`@P#WnG8MEo%BCHh zI}mwmjlH*Rs%F8n2BlGPYk2h4pqItfsUB~PA(W)#WY3nHLH+gM z>!y@t*z8KRtsV6vRhXzlMX_iW=!w*>%Rs^`^pr$*cmVn zhJ7J^_>5xtr|MS@x=x@l?V9pA=acBxS*9fwRS(cZ!<^;IYH)qFR5z*bW`>ty_VgNw zwW2j@@7HX-_#W8qZmAxAc>|s36kaKm)q)ZqhVy2!X92gDnwPuBLHof3U@yD7(A9G(^ehw?-g^2uYUaNpAtW#YVnRsClDkXL<~&2XgYE4|DNX>K z&dLOF9TLqMGY-&)BktPeJ}Jz^Myp6nG#sU{sKz{}=?B#hM3yt{U4gw$;5*yc@*aHM zSspOq9SjWo@$HKHEdfr_2hM-^O61u4q=_S)n)!&CMy%{52mu&bPQxAGoEp$CUl zl$fRnLOnn`?m>T?Fu` z{jmCh*iE$PrYk*rZX@cnUx|XySgZ4@B%DR7e1jNW=Kl8@2O0>NCO0dzx3 zQ26o1mZb0wP-O5q`)`D(3wHh(+qiZpn)P%X6#RlwR+GEY37+ z_w3cM111{+jSP<=z0)fj%iKr7WNFc--AxME!F|eSxfAH1sxB4J=9|$-;Pxv`quF5lovt5#z3TAR0oJgW(GwsnIq&<0XOlp*J~1VHbse(y${ZA3 zISGVjx>btPF9Mh7Y`L6=p76E*`4qpy>e$nzH*Rv5o5ia#CAuszv5!ahDLVL!oo^=ZJjni<$)-7`4bMc;(y2b7 zhr?Y*Kfmgn0*kzoE?nSzgz9ytZ_mjv!0;amgZGv`1dHx`6kv+=hI=|TpQKq2u?nri z2+Of}w6!1EhYx-Rb<(nWa|No{a?w%FT`r3uqwSsjJK0mgGL=_?i?R{YcwJ)G?EVy= zS-BOh_VheBz}7$Ybh#l678WsTXyk$+BEu^=9!vuvw($JhHyhB0%h&y1(=3I~@q@`f zJiCD4*JWERmixjj>U>(4BQ!Bi%uY%@z>2BQCBS6XGnYm_+AIw}ksa~F`h2^|gwm2o0h*-1i zRv4E~fcmjB^rlxeFy)?)Usp;jh376w5BLb42O~3&Ukw=7pcb3Nx8DVxqNz&p_YMK) zf#2aXty?!(LeqD<_|{xl1JmCfzcD^N15CDvE^Xe}h+Yd`S1!N327akdkUiq}1+3Zo z{tY|cCoXyK`|A8bZS2ETw$4kzGw7v-qwQeuB#4_+mziOO*rrQ2dtSfghGy|0f@S5I zVEkN)Hc$V3^sKFm#_xDDIxvP_OCHYx`N69@O}yPH(AQlnvD0tr8oceTb6q8O(DHP3ZM8EuIVR=fJCt5Bd5o&q-mj(p_t)xtPv zz@7Q&Ww2z2O7rHsIiPuxS?uZAT9k0Kw^qdU2?|SoJAUOs8mQW#+^u213-+Zx=hW-o z43o;D1(O`#hB3P)RJuDy?#Z*w@LC_I)0r_1Gs(ELva4i|*1gczJ-!B{?P+sBDuz$jn@eT&|lJL5;^~ zSd8YL{A@Dl(7Rq;9OexR>Oq(cpAoE>cui{hJ_{}cJ4EPwXhcsj$hf%1 z85pv4mHQKiK-lEBc{xX&4t9)YbmD>eS)?PY*-7U#4H#pfq^qhH7S?DOas3K6yeD|l zRE#+Z+$~Rf)wHJ;J#_iZaZaoOiO|$WXO1TVJDap)n|l4>eu4Y+AWGxTog5wML?U*~3u7`)}eN*^d17X^>`#ocRI@mY8#~cSfoxXxQJ( z4lsX(m{C^e3@Fom7C-LRfV!XQPHhZyfJ%Jd<@29wL$0}(LHox8pqKX-7w$YQEGec& zOXo{2ve6B-eK|1&KATh(3YkF+tvXoDpveW>G_E<_N=X4LOY4EB@&g1tgj?}#JIGaA zv_D?^B)-#B{d?Z_y->&4H15hJXQ=QsAz&b@g%pY?En!YD}kvai&Hlp1GHY`a)iIK0Y@E;7f{~@k=ku(Z{vu$NtuhCzI|YHOZ0M>D#d&YN#*HIc+~&71C%O%Hj%F>xK5d=$Hnd zG$ctDS015DFJ$#kkzC=D+vD_KqK#mK@o5gsA^?W_dsPMuL+nU^+%@tqLN`elJ!@cz2_-iK8w;Mwi0%Fo{)Ail@Oweq|*h&SNXUB5l2@f|;M zx>{Gm;gtTjB10Q@xN4s*ox#d!5VPH0E_iVrVl8p4S2l2mi=?V2W5c(>Yf8$m1wQ%1 zv=x`xibes(vgC3&lTSWcukiF@?9)l0X4`J;=?1XNThQrbS58RkWGqy4i33hMoW~#b z-AANY55JQQ7@7?jy<8%494r<&^m*b`1eBa&?>9T{0WFhu*zUmfXp`nDxN{P0*cE8IEP4}W z*}DeJ7A}S|Yi;dHN_-KNQdODLu0fY+rT0YeBecG=Bkz9IVUT)qz_D0jypxm8mWKs(-f@tw{Mm`cc#B{lKP@SC_=b^7KVAak z(*LnC-<19ch1?v02dEo=i zwXdz<$=wc_hBdy>ihuPr8#!C~xQPM)D^5_PPv*&TzkO>W3 zJwIQ@>KKpO@9pouz*U2Si=T2z=anO!Qe9P0e**BI%Z*>&6$N?a%U0Z;@r9_+L&w^A z5;Uc_&FP@~Xmj_|#=CEPVeT%Lv^zI#prwnFe|w}aeB_}rh;ePi>NH!JYdXAu#CQ1* z>iEvmjvY=;Hco0-mQm|#UyFJ0%0FxQ65j|AoyrZTmeim;u81kmo>F9LR&~Snc0Bly zooZaWB^oMi1D__h_`{=oiwoU&CV`wFTe!c{eH0@8Ls&r3ADWqc^;u$Q2WOlQSGEuM z!NqP;VK)XgVaBng?;>Mv1D2FPP78b|pn5`TK>If}Y@k8i`=c5wTvg+%AsKNDIh^*o zlI>E1Qkks{>NnrO=Z*HZtJubXV&RQJce10QO+<*(VYvW!ZYp<)h|UzS<#;z6Fnu45 z(C{cm3F7U|=d1P+=eI+j(7UPufk&ZktUYD_h<2Wp&Oo$)0e^q`HWPF zTr2F#yHA7nVr5~IUJmKu=SL_no$;>F{Wpkn>${9=ryIf9SkM~lx*J^)aIPwDkA~-K zQ!Cms1K=}v>6|SGxL|6>ZQ~bT-lM41pUM>v1;Elw#|KLu?tsQExh9w6oZ#qq)`mz+ zAxxsIm6q451sxk(?Q{RxTW~OV|C?*_64S6mZVl&oq|xDhd;u zxnKMPL)WCOzjUui1Pmqx(!0l^VZg)2qSg-qFe;<+=!9iIV6MJ;uVilyy8m^3t$SMl zjK6c^hGmI8th(I%gehew%#152Pik3%g?~7gjJ)yoR5^C5viy6n$y2rc!!9Z8(*A>C z(=XnG3tuAC5pyI6=aCz4w!4XP40PYEt|~*5Y3m;&!_(lh-^5wQ$!G{>B})}o1;Uzr zR&0Z>2Em=}l1@DPt|FahVWQ&9f$#xWRC9xm1DuXstd+NGC#*B|OK>TWz~0@qjck^# zLB@LDz1zg!gZQ!b+3mBE*!@x6FX9PrfPBW@IK75g5U>xB`uE;M@XLXL{O_e`P0QP> z6R%GJj!%xso8CpkW2&Q{=T!os;g9#O;bns$L#vLbeANx4YIWzhxJV$p&nj(S`OX2B zxC!b?A!o?mEIUzsbPe{{d)u}yYYZ7Z7&J6;7yz-mp3S}Lmc-sRZQpY!<~_*Yp`O2@ zJOMFwuWVsRyoDsP*{*1+lpwCIW3uZH9R|ah+i%d-M8n+ze5L8H-h(~x z9pl5E>uA%rlELSqfl%^8hL?+&Bb?IgD||TM1hr$%gWOqBtg*i)TIOm6isGpYpN_y+ zj)lF>)h&|1u5N4?S7;vvMai8^%RXH}!w%E$)ig?x(o31TZE?3yq%mhvT2d5Hyi|Jd z{F!JNf7}c`mkxybyo`B6uLi)x$cc@tkFKC@gY@THIRoL&7q4neH5}p7+P6!aCLLip zoB5(gKm_};?vwwQi94u>86f-p_=?ojkHtCkYq0b>vDB)uA+TzyLB06Fda&s;ldjA4 z>nK#(vrvld2I4#YLDg~s?-$g^TfA+e;pQQ(nG*Iu*uYy`%(}N9+~ZdAde3wTB`j*L zU-vNp^0Mvd9+z{3CXUyHRGA$h@71h)tvjnR?A(B6)}Fga2k*;bl->c+iQHQ4Zc$AB z^+aXL)3?C&;-?#x9zKYfhhwn!`c+h(bRuGR(lx|;H$moPV-%(^RWPJ6- z;Ubi<_2Rcgr~N?a;)U7IHBnGz(%GD*FaQoRtWT~Vc?;5Kd(_fA^U>&lZvU%Tyj_Uo zqj!JoKv~%%(a_-U(vjDo{@Q7m)RnrjUsiXhD|Iv^^2_M|{7W3JPrSKtClJ?%~VPOawgZo&)ekcBc6x_AR7U7Et_Hp zl;b){lnVZNIw0A)D2GOpzx+geA%y~DD`Tfd$@aIH8YNp02Q^Bzm!HoS*({0e0}DX5 z5bCj}6q_Voi(UY-U7na0b)q2t9}gD3F+lxz9ifZh)xjUh9U&28SomrT^$bcHdg8Pi z0@#1|c@pZGkTmpEQ(64`PjFBVnWUko8qeb2e{O<$Dkcp*)ts0ppH!fpuSr8swX8+X z#5v23#5v2<_mpYqsdmOm5ez@bpR#@A=lk9CRO6%BMeji_1`R#c0_pIV;+J?*u~b83 z`S}7pDUf=PYL%?Nd8M8$h|3#w`3IwlU#cg8hAd_iD1iAq2ztu0j6z zXCgv8T^xf$aDn{Xf`hA5$aLi2KkM)R%I9}P@vNvm&?8E}Pc7rmv)X^l^#8)^-R!RW z|4bME?iZt|I(R|R0o9cMxA}pWLX48gx>!z2v*4l9($mmTof8F0?cp8q0^okj0{pA{ z|2dypaSu87NXmms)g}F(`_x1)h5OWh;rmhg}fGBHVQ1TU`HjH=WqaM(lw5 zub|8$ba&lZt#pW5aA$g220+nx?8YCJrQ=x7wKY=vjrnxK_T*T zHwN#ddr?VO=cU7kEivKo<(;jUw_#52`{$N-6eEFSU*hP5hENWNOp3RdBR2OYSvh0g z3d1JL+@H?hLBr2w>U;Y?qkc~4KB(@7-G6Sz5Y=RlInlN8bG^BPiVa_g7ck8tOt5h) zcda)D3Z?8D9d}}+B9~qTe4fGR)76w(nt4=k6}zB)ED$>r)Dmjht_Pq4Q zW&^~QRwRsqNTKI2l8Fs-4QbgFYjzO(uDiUC&&dyi32O$XZAq}rw%fbOX9*VZVO6rx zfhg?!1@~`GYXa~$SCn$DDJ+HsE%}2AJzQ8%N5`>@(HQJpXiwkiw*lB(iY`|EfFEXM z7)+_z@?e4ufzRf*AI1(c+nu1(55!2k8SxK#MByBF|LUFlmSH@1?y_&_j>9~{uRP8R z2*k#&*K~9U$w6Kty)8e4__3{)F3E>9k7ExHwcU!;3&dPq_if~7RfSINB||hv`LXv> zZHawXPGUum%S0vn0x;6CA7vqyny}yH%6SdQhu!JnJ(%(}9^dhC*$o;Mr`11_i*e!SldpR@MzI)OOQ@*rbon^5-EU5Jo1(rNe zXvMJ$mqE_%YzXzY$nY|v7T z-M0&t8l{CM%PPYS^5=(U;*+r!pcrfQ!vzbq4UtV>7|=FDS=j2SY|GABdk3`zFtkVGV@P@!@{2q#5RecS0cSMPg$@B4l4d!6gM zzCYge$Fr`rfBV^cueI*Io@cG!Z{2q;&x}6U?#5;&<5G|%%Vo(lV(0KDi;V2S3%;q!6_oh zNoVntE*aIpjAF&My?L%31W)(0 zL<;c}ymf;Czv`wW^ngEcF!m8MQV8zZeqSjN0_Le6&7H6Z<29Wxm)R0f9(208ri!9# zWF~W;C_fmjxID-cU=5?jo0Jz#5`e~Q+lN3VZ&b5YrI7B6Cm3(LKkplD3AjSR^Cag~ zq@O;rI47qL$NRbDig!DM?D^b9-xgztporaSf1iX5%3R;xtrg)CsRCB9Z^ z{|U&Zw;2%3OGKtt8f)810T_Gg8MpC}p)dF%KDW#RDR$mQz9|VvzuI_xRSyOyGbYoz z;`Go{ix!^^w5o_&KGKRQKOTiIsk3RT0dUAC6@)C>BHi5=C*3DbB6`mDGvjq}h$hxS zk%QPEZCRDorgYsMDZ1PY$=qs*nudz|IJ#m{L(45i7hVA4w80}Ep8KKc=5}ANgLV-WigeaJ>D3`ByXak<5d$h2JXla!bXlAMy77eD8QnSV$-P@-pqQy#VaV#Bx< zmLvBw30?StOR!Pv>B4QSz-IQI)Pe_|rp#9Kh3)`cn-Gc*IKC4ZWRzD$5IcpoKPoXf znSTlQj_$cJZm0yM;)h++Htt8Vo7PGcT7ARJHkUQDZOX=*B&gWgI<b5t5v{IGQl%pNJ9GT-1>{}-GS>ct#8BzBXp;1 z{ov`%$Dz1}H^ecMo^aVAV5aZn19^**U8yZ5NGb1(ZuO`xv~O@JC_cR%$&4jZ%A_tp z%5;%P;Fr^g%+Y5rPBMXqyyJ9Fzw#h|A^m6Hw_gPPXinwt$IhUsd$XVUB+S6@Yk1%e z%0ATa^w?Z{UNEdGu~5@UBd_?}yrV;x;A2d<=X|yq+I8rc-xM{RfrFdR%gzYUsH)hzA&)31 zT;GV)HkzU9oBOB=34Mf67w_BX4l=5J>ML)@5(BYwR)dlEC`iN7=ZlZ17tH3gSMgS= zAs^E;K7*&RK=ysN#p9P#=o)K>)Si#$VUbTZIE(oNa%o^XZq5=9%-)nC$~{9AS~fR4 zF-Ek9*I*S!yN z5pdve+8C0eAX}r>{jF@NP}es@za#o2?9`5V>2oj&;`9S=7&x6l9U39yeK%5}J|=sz zH_!|?cd^i^vR?-AGhS07=gg4u)(kr_)-;&XJ61I2ZUtWk&cBzPh=JfRow-q6bL4vc z7n%JEY0y?Ac(aq%9x51I2kv&n!g!s?B0sMMns;s5UM!FXTA9Taua+Ev`{DlMfw6Jm zxNQ&Lwpt4mTYD};_DL!jm8h6yW;uiI4QIw~rFgLV5F>n_&JtDs8u#0xVk$IrMbO3e zxez=yH?P82JQ(^vpZKg|iHvv${43m2K(VHq)j7=--d|1T+US-5Tyds}iawU;X@8?| zp?NZtYd6GDx?JJBO^yoNR05D*Fm3X>WQol44~!WzCc!cgzTObERPVpF0EXW)RlX}5-4#opG`)zoV;1Ff0X5owlG7an0-60VL zGyR&3@3U{!T-5 zrAjL4*?ALKA}7wikJm>Ta@)d?b1c+e7n7V_Hb6a#a})We4ucLx{rG9rX*4N(A=&vu z42W#HvG}o6A2Bvaa#M$>A>|E&-2BhYQRJokvl%Uy0bLvp{?XMCjsZBJ5+~02e+)W2%dD)-Q89$3m{oMCk+_putO(sH@ry{}n?i{>$poPl7 zSk|Y&0c9sHHuDENp#Lvx6N4c?W^eyln;>GT{_IusFI}5BL^SfxYZEJhmVPc16X3ro zuLM?NV-Pdi%haT8I?@C--}q8CifaLJnB-OJ(2S-YoGE{z1oIpc-9As6fGt1TNa3o+ z3J&tVWzCcejt!^Vlt7^HfskQK6Li>LiHLrsw1SVH^z3bhcuE&bnG#r=Y}Q;VZi4CY zbFHFTk}G)l`j)9?gc92KMkv8l)i8(Mr6!o+3u6$Ar$+?N=y!Mx&YiXsr!!T8_^kV} zm1mlOec#}ziA@SC_{k@pOt_}!>8>LJN)Yw=>Vv1^O)xJn?e0RaMc@JL3qf3XrL4b| z?-Us_>GTgxF5iNl?Oe6n3IB=2_Z#Q8?!!s6>f0@g$Z$*GXwp~w7Ce{t3|cpHhQLv| zmWTvCBCtH#c!~_NJ>CM`zPI4=OO$_W;4Fc|#n@~)Jk~_H7~4dKvJlIDE3sQ3I@G=X z;2k>xchbDlO~xzh14D1rAAwP&a*e6sMlh+f>nR`;vpf>(X!Cw`yiz*yPT#2`;OL3? z+QS+_h~esikC83}PEK3)*Tw^6W(BoKM^GpFI7g?BgjWxg>(F67cjot*#yz@z!CC=L1KI!kQ0F*79uj(8dAdC_h z-S*apz`H%8vu*I$_f^02uU7z;fjt(XztuxxewkTl=XnAz8aX}afQP?+_dzxnAed$t z?m5XkgZ0B*q(j`^T?Q+mML-TAr`3AIgOxk!| zS$X~kk^ZfIc&sH^c%8f~6#4Kge_N;lXNO1cEFTf|&!Mo~lncP4$OY$?%nn0)XmPKV zUJaDY@sGN`B-#%?h1q2X;>heQl|-lvXy{)T^sT6d^7ffdpLU}CDOTnDv_Uvsu1#0e zHEGb@zv*Z)O*PEEb?e2_Jc;{b?X>Q}I6J1fJlQ4%`xtVLZ8NR{WzyHmWBNq@f*NPg zL!10q@nH zf3PWYMI`Q^YVL_}JZr<1XS1x5a7H`W>`q!akTUm->1h-8*JF2?*&_l!{#kf~ujoN| zCH%>A@?aU<+f~lZ@tLq+&G8tXO;I?v%PW$DvIKa#MVnAPDh00ftafLQIU-W!p36!V z(Ks)?p4cYi17IGx>)|`Q61awud0*DqBhtvQqV=g5e7NSmNW(dCQ0!VO5lveRx^1k} ztTMKU#4LSxS5z#15%*0P3={*Zz`*dnf@{zfro7X6zzUICZ|ciC#o;;j7u{wPMWMaE z4l5%Uf{RL>m~g5kB6)q#;WLiMIjg6{XthLOrzO6A?RWw3dK7Wyd6^^9kVa0gdIJ92 zKKTdQKKlVAw8L6s@*$?ZGe^;vc>j&>>ai;&;=LWknrp5J!U?T!lh0)H;E)Sf#3hx} zh~!*dop~e)=R8pxY}zjXuUISM_-1p!{6I>h(~vPDxw0?SswCr?b9`6v&3i$_ZOxj^ zdfA{s_ee#l$&l#(y8;oKDY!-^<%a+9-Ee#|=f%x|OmGdHdnP=rhe$=+4d+c#aSgwj zk&j9IFj$`HI^mcO2hN%eEIYY$i^>?}e$-QFh1?=bs8sii2eKQM&3g$|&H3(&R+k6?JFUn09a>3e66UH+}|K z_zjQZG-H1J&EfIP`&}oYci&^*z&rbpuU2!uvYHt-dthq{+dfe|_~>MMe!T_o+hz=3 zb2xx1tJo;ahioxQ@Wv^tnnO7GWiC8)ae%##y`NollSXdgt@R4HGuFSw*Cg2q@UqSW zo2_oU!syAA&Z~{`NLGb|ns$RH)|7DC%leHHK3i|f?PTB$%H5=vw8IK$Uh!JD@<(4x zb-bBhPf8uH(7$<^#^VBve3x3to>o8|D!X%gn*uTYsA)Hy1}&VParSr!br6UL+NM(x zzh}E;r?i7qLNIQp3x@qIdU#eYwas?lP+;{N=X&TQi&_tidaJWs!g4nH<*t)HiL>p^ zovvn!fWCzzH7dg5NV9%Fha-P9CTp6Br!$S8f6 zTe{yEn;G#kB&a+gW{KS0bRWOkiNomUFK5|5#$%~cFO`*gEpS7vxY=#(@sN6``}J6t zGK>~~>vKJpgsI-++%TSFg;UW~yRqL+gn1u?a$0nNf4>7$Xk7}n_ry~19zR?BoZJ_l z-s)r!ez)^z2cg<5%8EFsM(j3mI@~uMY3zXW-jLC|SeOc5lFC{yW?De&PMNpq&6!w* zDjORss}p{IdW&d(Qaa(sGvAqp+YY2cuV)|8%f>E~WBZako$-!CCGP$HnXumUH);L- z&M?bMv$Xs<2jlq=zI49L1#j(H67N5q1;-}#i!O7zgEC9UaP_e~EJ{=Js`YnQ{FD8x zM*~|9XdeA8w8O+3D0`NpnX>Y+!moxVAsp^FLz{i9G*vDb0 zkq52C>>l`*V-H6sU+04Ej)FXUl|aZjd!%?ut`K{aRP5mX*#l3K_kCk{I}dVC>XoF{ z1%WT4D9?P@HSFn01_4TyCq91Ot$%Y|KA~#go$4nM3YG7N3z|BMu)&*bFOZ@auG467 zI{Q=sC<^b&X>1G!#`iBvF42}?7Ao}eo!7nat-?>!VtB8@S}r~98O2DTb1@RkmM+C2 z&f7YM%y{903v2aH4PS+1uJEBJZP6elX|_p`avgIFX(+1O=#4kVizs;{6vEvPN7IG4 zV_~mNfN6tw8FrR|-Mfd%8|So^(w3CJ21C?e+#l=4K_7jHo?37@79-su(D&2}7cA&x ze>-#y1a7lNWckMf|Bf0yIsXbwC`R_h!3Zy$afe=fqgxSJQD!*4q$hwn@3Q)P@k&hK zyy}}~E-&0IL@4RWw<2J9vXvp|S|YU9*yuH6R$@1t-Qh& zDyUHc9n(6)Dn=>rqJ;-HPO8Gra?ErUhP&gb>JO9yvPz&_cfmV*ODY5u>NdD^S7DhB z1rif=ZaBBh?gPt1CEy{crz%{N3Qg2=Mq;0k2V>;k30UjMdL^=caQ&^FQXaaUsOo+(>sygN9wAC!H-SG2Oa4 z7Z=~N_|n{qwEPc6a7@bB-(Wr!7+Wc?XmTsC``<2@ljAIK>$Ox{I~9xI3A$a}XPo+n zs@XL=?lg>43_*afn5e7-L=kMts<^P_2I`8{&)qo!v zY1KaeiDE-(<LRVA{jydf&{ds5;9M>SyW1rRJV!r>k9sU0~wi@dm;JhO%3vungyHs2TQP;ADiCl&ccYO$ArWDwd;~dHnfhto+`qW+ zb9cX+!eXak&8B0q(0Ts@D{rnufT1Npi$eB-&XGdDRelqdnPNr>CnO=yG)>DhZyJZJ zUl?w@th>>~kdV%VMvEX#<0=nQIu!^6$*u&sD{g~9H%4@2b5K`@ug7u%=w!rx26F!aqcp- zKcsx7aRHPmD{!6tt%o8^!yr0GD-zEhe0$P`u;=Le~kAZ z2K2|k|1s!)jLjco`-ds~VSRtt#2;4o^r-gLzen)DrvbWWmSA8Y90_q#F)%PFPZK0- zCPCT+!zTfzZYT-`QdBu9sy~U_jfIm4oiR+|`5xFl$BCSu~7x?ufkQ-djd z+MGqh_xE%2a3<=5!KzE?x>Ed2S9Ebmu?YF6bTyZOVTvD(CtXN33X3pZSEN&#JAJT# zz@h7Gt^{#N6mxp@OJ}VAn>wS3_t$ksxR!rUXEd)#(^-(StCw?#7uCbj)!&_HZ+Y=# z16_AY|D(E7B52cf=Pf}T5@ncP{n8!V|EBI}UH$92BUZZqp6+-EVl;wYRCnJ{H>!)9 z>$LLNFNNyRm1qB(@@S~ELe_uX@H~mDY06`!@JCOPqU6)7U&{MmH9VTP{_Dz<{P&ed zn}BGVaPjvG@O2BHBFseG)ZD~mLfoR{WTNlo=jHn6G(ucNlY~Oqb2L}^-WOY6x~gS= zQ}xVg=D#u@G6YXC`8W6Wtp82>q)kqLUB&;xKDRlCxjDLe1-ZF|c==Na4#&|!iz(&P z=fbnR9|P$+{K&vStRIPiVfu$pnK@sQK4vETzoO%}(`&k!^C!zc*YR4Ke1FO3(2a#O z8rtRbDV}D1T>5n1!$~+WnhP`Z^7XZLChq2c*Ll)h?M%~}B1f6n&80=Rb1u45a{0{( z(tH_h{Z3=Jjc|)Wj;?M2A)XRJp8gVB{aqs*sm^|G5*}Q>G#QyFJcJ`6Wci2JGg)i7 zkusOKnC5@x5Pwtm4k60Sp3!6@ofK`Des8PHXitqZWtsA;`D^qxOl>8m$WLPt)f5+6 zT|d3grio80H14!Q!bvGeI5%zGPZZKXNH~cC<~ z+hk3`r`i9aL(-on_X|$q)3!^rLee4eY5N&kA>pJvw9Ve{h2}zkFEqD6Df-2V8mj#wMZ` zQ-#DQ3;RvS6%<5z*WQht(- z#SDCsj^zwE+24&b;3R#k8E}%m^$a*k$7Tkcq+>e+PSV*l15VPJiYgN4$bOM@>}TMU zbT-d`ll^s=0VnBinE@y1JI;WU^qpqFN&3z+;3WO6GvFj0ml<%94$WuN3fW(h4t+NL z2b`os+dk6@iBHn;po4$GNjjd>aI*56f#*McKYP!Bll**Uz)60-GvFjYzZr0nK5b_} zE2P{ceSbRm7o4OYFau7~51fXRRSl@*m$*-HKJDfx5!CYd{LY%jT`cCn2qz3&N%>8#FYDPr&|J_Ts z|I14i2>wizh4kl9zjNK{=tT{2^9Z5~5#Ud>{`JGD$th|w^d7r=QN2QlFe)u_=NJ&= zPmLgsX&Uqo_M&Of(a$Z!)8BPUk+Y^%;qD(4668!}B7%2QTASKw)1I+R2|!SBbMcz8 zQJVgKT}0cCtE(~*fT@ezs2*Naw<)@`2oX)iy3WB~E~d^QL0;hoQ)ihd{PaHiO#vNU zoLxNK9D}{KyZuX*(QJ*T?I~Wq1ScX`5bW<8<`y)?l6LdZSo((k*F z34D8!3+AW0KVavE#; zN^Rw3)Tt+z+Nuhg3QEg7UES5Rl@yee6_g3Y{|n3iH?e$m8M>;-Z7=!!`0rT$e>;{> zi)_*EQ!=(si)j7crbn`RqmFI+-Ssw+fy3#*EX|uUPtUA`8~Bf(u+hfyG7gUWgP<-A zN|gU9!RwO^iFD(nFfdHHVcO05*B%^^3yY@da09tpCC!cDpFcX%r^(Jl{5cmT zhU{nr|BL0L+ULNjRRarc#-eTUiB-e&(%>JfhG5YQGVDEr&)&n{?JJf2pyQ4vA%@rL znBpdJ_aaj!63kt;%cFW2B+p{BUC5<`-EKzCn<(&x6jqye91Kk$v0O|O*7Sor^(P1Q z9^6G$>MXqsA8#NTk;EWNhhQjnNa)eNF*_W4R(ZQt!x+$yxZYGdJ_@K4b`v?}b~x%p z*`nE=KG3|p)(U$0!3D?pH@qGs`1f!t+q!0;`^9qWjrxXxNk7}qkNcIe;sGtS$y{&v zEa{1RC}$JkZr@Wny0srvWJ}1b*S?2jSM~~+y}yAj?%NJe_y$AS<82<_+3m5u>BTok z-i?75ZKYgHA4b8$`_TnmPwgjI%Vl_Jb4eFS+EryoXLoYsRN0 z-atvkMMZ2e!O-4ZqQc>u9e!$d|AU?KI2gL7xTmvy6bxEYzn^|g^uORvnvJ~=3|7EW zqO4SSxq+fsdY0fHceQij*$k8-c(zB=bOeO6@T@Q8P{GBcn=HN){LPePP2OH@0!Cd| z^d+78fj_G_;?k=}3P}Z$SAuV#%oE-+x(UJ1q||?G$A}#+vFKlBv3(rOGNMMjX&nXO zHqT!^tGB~7QEIpNw)j9z?fjT8tEn))Q~S8w7{PzCgDcb`3n}g{WIh)-0zfIh&H6HB ze6-6(;8ngil;FC0qszAm1PJdc1_Aw`STFSQmSgqk)bXFicH3_tuFBzC+(&|;p^K-B z{#!fTmS7h=ZH#F*L ztCh=Z0wU3_8l4jT;O84Pza@L?QGy4bVTbq)RB}r%=}Sp4T-1E;VNZ=6KGk<;i{hPe z(2yZ5@cPOqF#q=Wk#&|GZVqiN=OXN1NWo0=QXUoVo+suMpt%B5s*2S_BQnw8-1Dn` z=8k};@8bC7{8jNuKgH)Oa=f9fU&h%&)C8VgIq$i(2$hd-(2y1RJ$z#IDAs~2QZ z;S#ny)uC$y{{h8qX?KqysoV$7hxU&E3;zlIu;;2+mcz)E^Ncsdq3CgZOcNNiQ+_^& zrynTIv);|P;~vtgy8U?UeK`u6_n0!#8VoZiulx*8*x_0O8*Q7N<3K-g5g$j+DA;mS zh;lm24&SxO3{mv=fo2wIC-b*c;TlV|`u5c;G3D-J7q&yk(0$F9*Vw&BKu*Ymce97p zuuIh}ZTWNFu*~&|i)3{Z*m=`SFYf6(Fl$Uq0+-!I&xU3FHLsSV{>hx%to^~TpNfJs z58C0blP|HU&^WlS(owoNbreLk4)7Fu+TnzGS>nqBePGVVp8`X=ROtNp*lPVl1b<$m zx|kbT$b_5W(5sarApexLS67QV-XC|-e14%f^n9+|-n_IK)Ha43GP?5)sJ%;*PC0)E z-Lg;Xp2#joHRAKwyjepamk{H8=Rl&Lx$A^0bH_mF(3)*3`$oY&?urYK*V^GVk@6|G z!+fCrwteF^!+vmG!6O$?PVjG%&3{fNuE_8nx29M9l6{J?HV-eL|fu@&;`(J)r zY;krIc=C?-cXUpW~uv7BM5g|{X;wc9l(J+?+alaYV~d?5N5f78vVTP zXzGW+s|{?HLrQkIucs%MVf!fX9w|y)kTMF^ikz6Wi_;E^WKoWqMESrUrS^HxdHkS6 zwAl{+D+GVOox!|jSqPq&=VLM-2Ai!igpMX^;F|8hSD&vE_rLlhNln83*&g4(1NrZO z`0eCxm2YmL7*TgTX?p{WuG;E+&5Pi#B4mDr#SXtzvTp5rI|5#GmI*&UF$yNekD-ZI zn=lvC(zvd@K5*H@d*_XNeBpEd-sF+P1pk6yIfcn%=#tu4n0WVRaGROb6C1xkRl&nL7&IO!wSC9D7W1XCDrM`L~vG3O?S1 zMcLWvZ-|WmeVA&PaCQ{9Jy_Lw_rfOJFOv9q{{bIZ+0{2!vcm^jw_JbS#;_7cX#_b5@D7wVS#CQ21EV&r zAQ$af6)3vzPUOb&5IDP5LQpbq6P~?mTWZAiVKA4ftlv+}e`x}T{9Ge8VP>UWR%&TJ z&`N|y{+zN8jMI;KGCP8}|3?$l%_FnW;?IGR?v&5q-O`cNA|XwjXwbW#qs$vlT+dw^ zT-gXFLgM-ayx)O&ZM*NhRzs-o%&tZ2c2}TpiMn5P`a+;vvgTg3&6}|OiE^EjXFh|F zcg+WuJRSu%)EpAc)HdN0nTOZ&<@&%AUMq_~ZTE&YFF6knog?O7rgAX0Bnug71&Lpa z9|97R{C78qYvBDB!LdmM|CF7&2F6#Kz{Y!rL+2CrA2ZA$E|FG=_S|a@&MvD!i(PIS z@yrc{$HNpS=FZ-Ps~t9}YkeOAFVyUo{TLnv%-bB@OZ#lGzIw<@$s!-nt{0I7Rrs8g z**?`4YhINKK2Yuh&v`{ZmmKqi-xT(%9lEZL;U_)aY}YLGY2)UmIr}~Vp>`j~suk*Z z@uz&FwU@nN$3UjXN1G-P74ZRq9+_sTsN#6uZA&3 z_YG!eBB_)9Nv1rXfVb2$?xbU?cpi>jvA@_GembMy=0>c46gTxgjg;vJ?;d$~4}Z9d zC}-aDYWG(nyT+q)*XD;pA+s`WuCF#&XMS4fPNq-5ZpC4(5I! zQvy3)Hi74}Qtz4@_k)u*Ct;HQ6|};kKDumS6)J39)PDJXC`{)zST!fv1|RoP@wnLd z0VMC{ulVFU28_E31HQW3;7<$W8AU((KxtkF>hW{#@Z3`lUokZmoZA!98z_>AmQ@zM zH&^)pPFsX8aB@9Um%68UwG6DhZZ;wZ?+i+On-#e4(Oy%A=zn z-QdR$RE2bkGRFC=FUzyi5w-H#P?p7eP&4AiGkjkO?+{;B{N^~}e{a8$yt23n9O+rQ z#WTAfh^lt}d}&pJZm&yU9eKS9%}Xc~72ykm?q-hVy_c-9yt?ki_vZ$|2mjv9o@d5D zHp4Ujo`cr-tA^s%>+^l#*tr%%j+JiESi?F&`>rDH3&>b;%`*c%U!v((95M(*M8c)Z z!xgc>8u^jQBi=B|_)y%!kw#zy%wk3!^aF0ds8+Y*5-KS33Jg)BTayQu(+VIIeIF-5_ zv+>p%1Pf)LgQt#5e#;sFUwZ2gx2#sg{1unZYa|lupYAWqjfWoqW&gHoBNzKYT;Z|z z9gron^;_ZD5jiKt@Am_v>@rryhhspOMcR1qwiSNl zw{J<8g)cMWoGUAocBdoO z4Hk{<^ZJ2m->Ox{Tjem*2OHx@ao(``5^HSNp+->88`E5`HUNCFQlt2lBGePL(7Bqc z8dWk{t9pop!TxQT;@4fQFxx8D+26z80ZR{G2Qda>9h$)L)?D8TKkxYZ$~D#(8qd;U zT5)$PeAT&1SnZ=MR?O#hT&a+bRtl!F_yoQKMWgS#cNEBC_Z8hupLY}emkUdHA>9aK z#X`OmZ5#j+UTwLfxCn#Tid1ZGug6I;b+goi1fcU*ES;N@#$jp8;hSykO;pfcK_tpAf=~*+7~XaY+HBoi8K7@>z!}u zu>@c8dA6&YCmp%hEL|eKq7St5IGSiN$l#idE8>eHydfLqqomp5MxfhWp_S}A0IthT zUU=en9*O70JW>5th3@j+ez-0s6y7S07U&`LU-~>v=ESvLVCxWm#C*v(5XcMfh$ZCG zx60x0x#$Z+s#yYSsLpWD)z^2gY+j6Qjm<9&)*nS*o7a|psq6)}DwZZL3|x#iHsZ6( z!wCPMrvIHwv=P*vJNKp|X#m{z5E)b2bROMIe4EhMP=ykkI5daVLgC6e+iNV3TjA~B zc3up#>jg0)3=>dk90c2UuZ-Mdg@5!gMb=#Qg=svk|r9V86E*mqe zKEhQzLoG zSTGKp9`DNh;A)LGZ(nFSH`f=+$@IDP#XCWr@QBwaL<)c29U>oh>@c!T%BY^5*#i!V zZIbp_D22=MwvvtjV*S&1sHBPL|NXn71@ct`;9`<7%$`$(EVy$;w-i+&j;x;2u<>Ac z;cEJtY&jboccoS5avaf*!yXU&3I01X{GNT{vcauh_9D@HePQ3xk32cnPVfeAlo`9B z6z=Q){7Fmv2>Ma@Xjj3;9&mc@G!l`Gc`+ z?d2DX&~bQv!R)+O%^CbA^@6k!mYO=vK zFY=wfxcfqV&WZVU>P|3Ray(d1TncwQr|fIlbO1f8$-Bg;)dTpFoM%0`xd_Yd$%9LX z{=@e;cRj8DJ`Kf!HwHlKYYh$`u?uMbq0Zrv?N#W**+)9=UxMJ$C7TnsAzN&*D6c+$ zbq`=)^|L!sXBjWP`UwAY5$*5FRe*U>*5nf9Xd~uVQe`IQ}y{TQ;2#l*( zCF2?gz{3+vnw6c!#D47@8&3!^{&~?Cvd049iPaoW#)y6P!$abk#SuLqeK^HXX7xCD zykl4gQ|<8eU3LaFT}1rL@qv`UE~5Vq7OtDZq%lQ0wYUDvHl(2Ay759q4~Ph`-|}FN z6wa^NbY_U)&o(>3`87`?2%vsv%R>X;sq{N9rOlU6CY;wEOt;1+e)IraJwvS;A`W#{UN$e?)*e?sJ8Bwu$p^Mqc(|n>c7|E?my&Cj$zV#F z z6e#J}J0jpV0v3&azp_-!0XxlZKN_>v2lC%;)Ee7h3p0(RVD#iSg;*agI{zh2hqnNYpS8TT&#W5V z317g;tPlY64m=E-FR%q8?*$>?dIzwb%eCDvaU2wt$^R_>?0~OZX!xiK`@lz0PEt-* z_OMGeSY@PA9Ao3Nr}PgWK_~XS=omWL2^N%}z0JEq9CI!bbvZ&M{C}{u&Dilqz-FyH z`?c&KFbf>^Jl}c-nIz?{FlMes?G+cB%d`FAu_i;?3$M1|?XTbMi2v3J9JbtE6z4k* zSh2O)wyG`IFugakrP~|w*j==jx7Z8?8wBG#7A?fkCTig&t|O?DX8`(`bpbb%SrIXZ zB=CBdgY(P?|A)<4_6Q^p_OCru;VCf)dN?~pOH9tAJL6|KZm?A&r;>*2b3amHLAYp1 zUZx|yE_3FE!_6*0JsGoSz3VtIyfbItt_Vl`=*8+h{%eFkcTTLbY1<4BHe9)z6(fmr zsho9n{YOx4=0kzXoGuVE`;A!9R|(9$d3BD$HbVdJjk;TTiTKx~c*({^gCJqfr8iX% zijbGY=F?@*s*r)Zimj*x6|Q%Gta4D&34cHKNKaR>8_aSZJ@D9k98~{U)?35xgdg2x zSUiVVPx>A@6W*of0P{H*nY9iq!X3IF>($Fr(OW06u9vU6z}lX}nLGy;;{FO(Gu9w) zDDrtNUmQ;pIMKr)ah-J#XdXVjbl1%bh|RPmd1+x4(o}OD+kVRr8r~G%Y8>c<&!j9L z=8WhDYL}~@KZ4`n^|4^XNGB)!`N{U`7H4mmmbNc%;0; z=K#45uxr=)OW~--1G`r+CZm$zY*+ua-N4=Y)is|(lKAL}MHRio`ePy{B9(~pfX}A9 z6{kK70QP%LD>!*Bq7cDr0rN|%P_)QqZ1cz$y6;f$nEllWvpobqCwsd=_0m%AjV65AHMcIU@2*?zIOB?l4rF`@l>its^LF3yHI_hmVVzA ze^qDPFR{toTlg(VQL1rcB1u@_eK>rD0Subg8fZ!{ zLa*O)iNAb$wYps6>H5caf_2p0)gQ+1zEBMws+JA+=z*GZT8r9h~Y7RtP<_(R8Ly_<8@$}a#i40{bq1c0scTI`yb zUsTs^6Itat9t|>8ikmjH-2``w;wQE>tAj?PfU7+=$Ex#NqPnlP7=t(NTZ3ZP`=fh% z?({M8bXULYVSjH^DFs6116By$o{f}+6+76iQ^5O3p5X0W%wVAj^P7{G8PJ*Pq&qjZ zi-Lr~u#HW-XT$vk4@cu+JgVhyeJFNeviiX8(66TaX&|$B{U;8+$H44c`{9zgrU>p@ zy~5V*Qni`5KEuO@o58E(bH!zx(MWHy{);Kkq)wy%HQVC?mLQ$W{QTY+6(rwcHPHJu zA2j;;L@eARhK`@U<*iq`9F14(eDrRo65#sK5ZGfX2Jd&WZkaeyf(nJ7x9%GhK-Q}c z9`Tud3Rq6+8nY*N04~9!^LjEuP;<+vq|I;lRl9!TIB+T36SR#bS8FmHLUlIYB3lHc zL2yD5r-V)z2xOZSxL9S0PPGqyqjr^m1cBun4(l7D_||=MJzC8X2X$=|^HF^u@`JO;*);o~L9M$>L{SWyYMH^h1Nj!YlmcR9x6y)OwBWahA_TfEX? z<)u1bdK^gX0|r&Z7G@*I#!Z`*asW7<$N#j={V-rm?~Y47P! zc0yY^Km3%d^+P|I*be?ow+5kK_uOEPQHBETOI4^*uTiFr^rpda4dk&;{XBJH1!!rO zP2{=q0rYdK<+1Bzp)<9)7aiglz`^s6V;F+NfZ-|GCywf;kT7LlV)zGh;7RpKjXixH zJlJ~NG@(BhRqK?T@kqZ7`eXa<2WSSP!KZoU41%%faxedc+dL1@v*%p@yT+AJL(N5S z+rp1X``m5k``U(RPTlw^`G{M98ek_f=i?`^WO>V58@k>OzJ z#GH$UtryV9MrrAEd1sIZ+BQah#Xv?at|X^59Z8>5cR9<{3``O~jml<4qa^02Y#H8k z6jS=nYDZltDB$i6EBDuhn@^nYXNjCd^OTpjKbN;a=eu{MeoLqaLNYJeIC(yTxNRGS zxY|!6+{?>oSgT?2A!9$VT_!Fm_1GhD%i~T9$9N$s ztRDDlc|K zZAL?<9?WvHeh8RE>Z(pZ>jGDA&12-`%|~^aIbKii0&s-;%bEOaA7Gb%PUZdDa%5&T z$9~__FwpN+{J|%?6Kn|htk`({Dq_8!D|!9m3-EKtAg{hvCfbtKohQHODr!8t@QaUM zGAK^E`e~eL12n$Qx^UtJFJ2enuw28@4e_(xV_5K}1(=A<`Qe#$8OdL@E?l3#7^QYP zTzOt+3erH|MPl28C|7Fm7$ZW$HxHP8IvNEI<*(9zv~CD&O@Hhu)Oi#6?UszXDE|hC z6fz{AznG8Cooxp1{BELMy`zO}*Jot=EGF`8pB4Wdg+XI9>R@efylB~}1itXs@hP!5Uo;Aq6{7r|?`FPMC z2TE*K{s6icY%-+A>ro~D(!SM2-Qdi57T>*;bLjRn6Vpt?dNitiN#40T3s9pHHVG@* zKe6ebPnMYXgNnaNY^VCM5{%GCFy~R#|I`}O1 z^!x*8)G6;)lwKWbxw>wL)b?aBswc7JK0iA=6Se(5HS-aAt`;R(qSFVqzdg(P)Tjg{ z#aqd>tYYTtDlN_R4$E=E~w1^+7^x+ncF}AK%GE+{yNe2 zQ~Y?-IO}=UDE&0l>|7Lo@d0RkZqwt;n)XxKLvA2_V5+wY{Or5&qN_iNBK~ zhXr-F96#M2iE6caii{H8f?@84?ayVz@P)~`B;ef#v==UY5Bh?^C$5{eyn9Zf6Qv)j z6#VLu@B#Nvl1yn}(T>$gI<5;K)53`P!oqC`9g^E6xauP)dGbQmL9znvDhk|&M4ur$ z(N%m(1{olkufJ2Ce=A%R?D!^YSP>WXY5mx~DjFqLy6=}Z=>zkgeJ|(;UW9M8pGaBu z^fTBZq&WBW;0}K!(>YYtr**gV=zTPsA^1$cb0#>955(@JE`|<~^MZG5Y(b{W zzsd5zAuya|)AM){MuU}FpN_OXM!oNcKygbFc)5|{?PcH!-&TvXEL*3F7a11!IUR{b zZGm4q8aM_(#X*;hdnK~C!{L~UxYHLfq_Ox;&b{5>a&KLJ^x+b;(@(VevSTCqXZlcSfA><)ZRL@fH{kDkr> z1hm-#_dGUK!I6O_W$U?_;gD9~@5|srWezz8l>9IJRjN?08-&vG?j7q?-KX$=eenAXVT6f0g1HaOy2PORBdw z96Xht?5d}Omw%N|Z}Ce;yN+J;%`hASCJ9xn3y-T~eV)}S<2oGhd3bQ{0-i&_Wl6$I zeaUJx`rwhbR&g_WyM+677=>I;Px8C;9AAYP@D$LzQz1zGPZnJa4l8LYEBm|v^E z0_&71fG%ksSUbzxb|5z$9FknW{HiKO^E|#*-f3z^tgmWQTUO-(kJ}vII&T_7Zr0gY_!@Nwj@Iw)@;uVxHW4642n>jSp^@UT@9@Zsve8OfbCU# zxQYEsiI(N&eTM!}T&VYLyPzH}Ze7q-l9`J9I<@>mRY$;^R%4yhfdE^1Kd>v$7l4`# z;yDLnbAi3R%}K6XHK^xaexq1yGm?IHBjQL~E^uKgV3xjV1=Bj44R#R0Z~=wQ_7*7< z;9$dv;NGzZ=w8Wbz917BDEOe_pfTT9aCw$|Khx>}*v2*9={l*0wQ>``sAe2SF7rKQ z8plRJV)f_6=Yk+k+WqoXlBW=~FHx-i9)1GsZteKt=~|1#Y=@oBENw<@pKam_*5`l( z;g3wNg7$ET$8MQmxe_!{3|Z&R_XT*w7@w{@-H7yBl()GfC7ApE^&YR?tngJCn{Y@i z!9R{ST)0>tL$L(!1uRFA=J%nCMwz4Fx`N0tW0Tdme^Mx&AzldH77Z&5e3%2`O5t3` zeRW7tljBO!!6syvU9aXPb{uSWsOtWj=>$XLm5#3ZtO<7~h0T2#to5FO{%~!GVV#c5dAr2s6yb`HGwjaD#%b z0O#IxWK?6q{yKUL1V27!$ke5a7dZP}bxaX~_K|KQ7i9CnQHC+;=w)}%drz*{Pv<^B zgD&V4uou9;y?Mj8{jUDc32yli6SoN8r`}#22%Mlvf`X|@v=G2JzhWrPA)7@)prV=zc!da3;WXK^6W3* z;zO)<`+Orh?&iv2RA&O&G?vF*YF!2kx{Lb89|S@Z&y=8@9) zI0lBrJM}0Rbg}*G*G6BL3qjX~X~CNkj(~tm1O8dF?jULDA=%;k_2}B8r$KLg4g=yo<-wkBU%&PS+I=A1r`^0nC8=YJgob5DEj7^=|2JHj^2VXEhYvEOfQRAM*{(4=sy z!mYcgR#TMKU*R4KWL{fCN;=tGn}0t$_JCdQ0nX2112ZM&}bZ>SHexO=GL59CYePR{-~iF)$W%?_xW@ z9xlJgd8xH#4qW7(%$@e~2oO$rllH9PHo936qGdw-{)M$|GOlt?3TQa(6M0!O6t)a* z{uZ8U2TjH`Wvc5we3XLf%6098q&efd8$DXamku6Od6BLA@@|KLnrCW{xt2D*&(YJMI@v}kR z9=F?QZ;Ftvx$!};ZSQI(skCiy^Y^&u;SG*(g;v}1Re=*A(wWuoEq4R@>aKQ%SJe?R z-DitBW~C4PVkIkwwgke`!q+DDqI#HR(yz8O=p?FiUl(d0H3kNGnJyl=0rAEQOE}uD z^202@-DmZeq=Ve}uiNXc*P`Ux_l`ettwRU5hIhC9NCva2rFVY#9s#{2B_3UzyA>W; zxg2&an*ehw^RimE)}woOm*qd3y%oxTRq;!|Yy=t4_oSH%2Eu*4P&O@07r(P`IQwcv z0gBpUR6TzQ@w{-Dz5I*IYMi*l++$}T58ThrQ(edy3*gVeFkgo|s4+P<6tde^j7_`C^ zieJMo^{WCP&+g#^m#tUhxw^KyjbENbbxj}C`)`i|BRlCc``53=t|wNR)%Nm3JEh(_=~xF&7k?+hKyAX0kCFU*D|&(5YJ0^&~f8K0V-E!TvlN< z3KmNlZ&@}6;z@y8(VoYAFkWzt!rGq+;HNT2=^SGIw@R;F%B+49Vdkw29uJeic;43& zVauZ6oW46tJ!QOL-W=u7^NYtoUv^!E(vLb67IW+VHFht!<#yPX3(}VGr+qGS@Q^>` zm780nVgzts;Q5_%4&|eo!e=VZ#P5cNv;-xSG_`S`!N*`7e;#;UnW3vKn2JPXRZ_dl zYLM6k+u8d%5OPqg+Gu6A7krK6?0v+z2g+U)yZKtn2QEIFO?ie#!6}}jq9?m+k@L#| z0c$bhcV|l^EcQECL3cjNb|m8ulODQpRP}3PAy7ZxY*ij|ds}g+Eo1~dQMO@|+_e%% zJDj*ayo?=M%;9!V{JI&$72gYtG^kz^UJQSwVg!elphVf3b=IHd&Rcy%R`5i<(XZn#I?kTUI>ilxYQe zA}2%%iAn|@hPf;3O`_o*mWPbHuT!D0P#L)0Mf{HR{mAvx+O6NCN4@9!9v-P^i^eHFeeAH5-*n*tj- zpFxp%Q9x#$(d3cX8e}y~OYL^Rb=1;+WN98q0Wrrvgtp^oXv!seLxLj!T7h}TJ{%nZ zxyKy^953BM$G7WoG*9}&NqNy$v(I)g_duA?u_AwXiN8+pOpZ2A?W4rhvt9*>RKdFC zHKTxC-t+L@Gpq3Ou`)ZAcbw48F+HM6CJB}CiiVlB)u6&Osdeu!TtW|dzxZ*r#DH3L z;g&w(7?{M!XT2yR0M73`ov3|(6ugP=YfuTRMYlgpQVU`OV9=bLO&iQNLznPHD-8t` z;iiEujmtz_w8XY8S9us-_o~o8CxI>zwpyM>oY$*h`_#ehx?Fq zC*}0V4*@Wwv3R*sqyr3&q}Xi#>;w}m&pu(VmcYrkyKDBfJw&FI*7ZiFU7*rxAUOWM z3^wBz6P#z#1$306c`N6o0~=TUT#vd7$X@qcOMv%HbmHu?h?@sf!B-wB%ObBBB2GM5 zCAEoI*HqqRjcR`j;=?B*lHcVKzbmMb5(a_rbHcoqZNXchL&yiENiG-oRKy&xc+24T zN0OHXeyBx*alP9L>Use;OS0V{uNI0_R*>fEHqrm1%Plw*zD?z2hODa3J z-at};#c9l!Q$WjguEUmoG4Sle>nrc=4}?2|B>nN#OZOfih2SH|olh}gIg#MS`WKlUC2(zhD zgOxV#z=vCkgwMEy-n_gHo}>iAe9nbj0cRcIr)IABYQ?S4A(<(5n}iJZdt~$`l=u-A z`>|5ykoP^na`utH&3l&ORYQ#1lD_r=mWONA{1sA>heeKn2=5gnzP6^516@Tj=Zg8~ ziR=Zf0Vkb#SH?h<6$_ZP&jv!n_i1}4RTx>6OtwIBSM@3TDdqCWz49YK%!7nx2r(V?h;ip)T0)2ThI8_)|j@-0jl zjV?kO9%5(q?+S$NnxBj%n4F+vYIDw6CnvZnpo`jnV-bG(I?L7|xeS>!^EHS3=mwjf z-AVkGz69TC>C{m&=?8b;`E*94s)ETOZN5#bFQFw1jH;gpT|#n&%W|I{+XpI2SB$qj zh=!K}&8%Cu2f`ZLnloXHJ>ZU0Pf^c>GiWh;R6q=2mq!%ie6uGU;esoT$oK9RC_6Dy zv0p|Kvy>hjR~0Nr2R3MwhF<%kR>M-qP5{&o@pV%c;g%B* z!Yq8-+MKFR|!-K@U_~a9FV~!&fj~g&AJFy5(y%xfp0%sEkXy4?XfK4ph3pVvJN?*r7>Gk|CD;Pm-V;5(;X z>^70l(|9A`R?BP8-uYwG_!A&8$IaAc(bwtgz5z3@0i`=KKc=sz!2Ow;N%Q6~5c468 z4kM)Q+7bD9E)%5lWrpT1)&5ChERar@6{;CAY}FH=1?h6KK@oZVjgt%6A)OBg{KmJ} z@Uj#qq|3<#+sE!{`Yf3Z>3q3i%?mEMMN4@gy}!KBvZL^TUmEYU9DK0W;@;kS>3q}j z@I&>w{wV28{%N@cV4()r_B~kw)A9+zx##4AUmX*imQx5Oif!;+aa?FxUSas{K!DrP z6T;JS&w*5i)Gb3fbEf4NfvW#0I`9ZVwqUH`3tzcui;2L9H-|8W}l zpL&n>>Hp>L(VDNby)pFx33~nnTE2qmLZ&PCoNkkeH>O-24MuZ$zZKdF>OWJMf4h?Q z7B+<|ylKnKsX`y~*<{m{I9N!e)|f7|1U6Dr zId7%_bi?~~4<<|@^O4AO16(1wwUdPWUH_*wFqPDY{Pq&pWR^e%W^#DB|I8~+rUhhRCXJNC z+0~WEiSOqe98Bg|VqhjsmTy`xB6E(bvtIy_{)NoD#K272F@wim8M=ZPm`P)1BC=*T zxsXtH`T>RgR~~KTL!KiIob$Jw+YP?tS<;00T)o`ggPf=G7&#INdtBVed`}F_?W6@){wEX5k?Tpa{x#`MclOr=H{Jdi zI;V$j$^wiGgyxBLB{Ks9c{Q1Qitw;gg;;h@m49{qucgzWjT4azD2PtmoL>FlJ)=!? zS-+pBJE`|A?d&uIre#26`yVw_Qbw{S4V4gMs*r|CwtsoVX#nZ_NYmtBcK?seQ^C{w zn5gMh$Y07lMKg;PF57nee(zq1rgQ!;6Vw1Ookoyu?y3jr^IGFc=IHt1 zg68V3CNp3BF`EgKH27rFE1=Db!*v**J-)$wewfq*9)PHR7BLP6)AK`6Ui9XI;1B(TlQ&6%YF2rhD0GUGks6(=}&8yKZpd z7~@YiM+A4{t&)xJc1!wW`_c9M$7?7sAwOB+eG4yE7cYFPI=BaKVG{Pdbszu_u|56q zd6*ljRD~*xSPNhS?fJ)){A2ORiSK54>jSaH$eV|!yhLGF!3pZiH$vDUf5dn0Ks??f zbXM#5wLqL;?P_T4Dh+q34IM4u6~SQ%yv3&5_u)^s4@R!;2*itS+6C|?C_+_1S-<9+ zbMYx}@$%7+NjRpdY{4Aj`c1RPUw$djgh%e+8tO7pJi%GIw%+RimU&uGr)L(3p9?$< zbgf+t9~WJ84SpelhZa1t%ymh@A?KHiHBD>4&kFQfT9sTp!H5IN|)aI}T7L=>38_uY{mVO~=cks3Uj}&*~fd^}X<@ z(waCyF&F56MRwEhVOiLp^TX0<(^0(Kwc%dUOLzQTa9jKxUQfty)yi#ur8?xk|LW@a z+%&9l;OcgRNLOqU>YW`{>j(9v*3WwIdNu5|5c_$pDh;#L)K9drY{kQ%x-3B?27ppDHTpJq&A98tZbzip?dT4yOw;($MpY!D%NV%~A3;ldiYHt+{ z*DLRJS#!$`-Y)j2w78Oi57cyMBo?p5AL^Dobf1iYB5&T-7N~ha(K|8cScvmhDOsE? z(fYV9ma1W69|z06*SiaN`a<=VlP2uPGBCgXsnrQu0Bbpf?RB(^hh^Wg%U4PG!=A)e z&pX~3_>AFtEX1XO9ks&_U+;;Bmepy-cMb)>-M6krzLw0uqG`)`pDQZjxE3M%rlEMq z*|$r~VM!ni=Hza;h|)3Ox%yR6_Wx<`JcFuex&(Yhg1#auSyVuhA_hP}Ky(h05s4x> zUCBA;qzHl_T;P%=sDK2?NQP-J@roe42r5d_3n-$JgMjRHcdhzre}7xGKX$67y3Rax zx@UTNX6l?ir$-3A6%MK>_lSZk5jW)CD|o=5sSjo4Y#I{pVl+F%#)D>r1)cP4BO&xY zxo&XS9nwbLe)%4siZ&FTbcD~bp>0teoi}?TU|>g)PvBQKsHiK9-A+zHe%N`Io-1tV zag$`^<=Sf?{o1WqfX)TP%VSd2Qj<}C^rx6&8ZjhZGq)>EEflo9>B4^;bO6PRcH@+{ z$#6C+VQa1`0GT`(*c+Ve4`v57$DKKBVcj&^j-Dw4tvK$_t(w<>!};LvLiPYPZfxSX zk0~tB1Q;C&N<|-@IN{*;S$QeqpW1a@KeQ(MH9e^enIZtfh)+2Ff-qqr;3b?AD=hmN=8b1 zMhg6@07McAS$qXfXv1OMuR6s9t^ecyQsrn8GJP!lO?L!rn6ZCP&p@JELDKy&Q)SQ%z%cJqHJIsK- zwD6$e4FxpT5O(0edul%>4#qEdzr+VsCvu(&+JO0I5vAf)bwtCP*Jkul6l68Ou`qXx z;Y?opjVCWT!?SHtIilBf&^J}hi_H~h!Q8$$Hw;S15o@tBS9|T^8 zhW}IqnW4ihlnn<*eGs^4!Y04Ofixq+DN1t z>7lpv5;FJQuY;jogw59{`lzxcXd{1-4AZ{X#-9IFM$Y~#xIz34hdMwD z_1%{Wsrz~K>ve^YT$A>$t5s3Znl?W(+-8F8@AMPbIMQIUzll$6_gPrau#rGE(9M8)Bkf;`t3=p6skUSDLp)zqoXDnfmdQ{<=2|}9Yl8@+oS^Q%(xJuA zrE!ei8^)!c+s4i&!;$7PkNWMl=zV>4%j|qA_4&|FKB(ahTDfPIvz*9KU7HS?EVju1 z=1GD`S_*i|mg?xRLE_6V-^->hU#O5-@DdHnvk- zfme}-qp&CiWInJ~J5sFC?K9KM8K2`Ic3)KXHw8!Vc&AY4IF$l3_ntg|)^3eLCI4pl z_#+yWvK?v?4%vX7RoXr_lT`TXcRp&Yjznna;7JV)|IftCZ^t#AA)j5%!J6xfsmc(^7UAd8Zd5N zz46#N9zyl+2$U!opuasvODMBmkgc5%H)AM>v`%@dGp~?rj+bh(yCFV;E0y`X` z+4ItudUbW={9`=s8&4u^FucZk9E{N2;S?i&aWlxfz_ZPs-vH75VEFNChy>qDp*uF! z5S{R!6Z@M{4&n~o&9M@(MAm1gww5hOaL%kTLiC&gaz9OPs((}+T`$OLls2|TcX+eU zheapAQb7MnN$Q$Ip3_6=oPC5Jfp$oW zqrW4{-U)sBZZ@l!8wUqyQ2Ra&V|bVcP{*oA3gqC zoeM6%jZNzQ$Jy|EmgQe%xBnr#Q)L2yNyuR1UoGrC=Kpz^!24eXf$@LcJNCz^@ZSv* z{x}u>$mG2RSsN>-gHB|$f@mmrIcJ^@{AQOOT-&r6Jq?_#l!1#bL76tb6>v-IlntWh zZ;lqF%iFYr3+kG6l;P`R9kNSME9jbz2d23w|K_VP-~MUqC4O1$JwxSl7i?UusQPNL z8*?Ym{N`f^gTA$)U1E0IS(Rb6!l`RFw8DT2$H!f2yAhSAJ}cga?;U1c>6}!8A4NSs zFLSiQJ6hY*dxGSDbGExv`|!Qq2Q;cHl|VF#WUe*U0t#iuU*iLGs2p7#5aPvq(^#i% zJe6ReQ>Xf`$`%+=EZk2wVfpnb!H65TOdFU`@ZI6N9B3CN6gOSW89y;$hQ^2X^;II_t9n;9|W^O zPiHD86+G8d!Yxfp2PY&IVW(^EtFN)mAky zM#Qzv@j| zl7gtbBlDJmBVPI3luYrJhX#)e=vl0zPN4k|5$TJI%6W9L#bPj{zVv-WuMRZ)4A&nT|JFB_A3GF) z+s`b%XRJ66io$q&VB;h3;$T(0^W*pV@sLTKLDVD#$$by{&cX8M)5YankHA$;gkDmU z%EK!jZHoxS8?HHPJ)M^W9_tB)F27pnEN1uH5kkFw)EZGz5RQvJ>#3J#mV?`*CR%~M z8mRHZL^rQf`w!xCB=%m%(_=4aEuA_G7YaxL%tvZq!qUx;!!CeYf8y5H@Ccli*@&n0 zqAU#8h8HXZR70J!PM?Y>_5Osvd$zc81GgMh;`5y3rKJ3cxFXEL`_Fe{XT0a+`j46_QD-1^R#{@jyVe@dwMlZYcyM`m)2U7g=-k}<&yI+P5Tp}a`DQ0|eDD2C0hf4O zdhK!Kx53k(%*H#uovjqS79)%Axw;^N_2^k7ES+gU#U)8IG=j&Veis_p^*jnq#c7%UNC z9?Pp|qsqaYwy|LFx& z9!xWB?@cx~LImbBeP|*LXLV5~sZ%uc{A?t)}Bx+0nd5 z`voBUL3#nHMNs80)R0WN7be z9U5^`M|sXfb`h0qoVE1sl-|}s2)*>SjhQ|PuJwLlo1Ilgd9;gW!df}_4d%DTZ&&w2 z_|r~3q)4Skh>ne&u3-!!dXk!uG)dL{V2_T%sPY?Q3XX>n;xW5+=iE=yW`u2w4L z_UKvot|I2(&c=l98xFoe3Sz(?oD!)7SV#1S9v2})Ope8Ao^2stcDmW zxHR=?SO6!UySrOIzyuCrf_iPvBFL~zu1B+th?$BE1tuhl;#SASCS>JoU@J^qY3s{r zB*6DY_5Fw=*6BU3!*Nj>@4k-RQ#s`dvPudO&K$BRcr77eP}CD+9N*OGeg}9`$ofFb zS#LPn%4(n;hM@*Nv8#dEewh9GN!3s{WjrSMc~kCM0MNPwZle^-qe^bN=UL?;80CY) z=}qi7cgjf%1s` zw77lEt4OSNcf=6^tt+^X?m3%ZyKaJ9__a%~h0dXw7=F?BburkT(~YZd$Bgkh;%df= zj$2@@h;!V(Cyr`IG%0f{n6-d_ zN_PSgk5MPx$)p;v4Kk*hlyPwgwOZ72cTL9FCFE#s8rb0LaVIO)AQ@=a9tI6rsKAp0 z4ZN-WDHu(u#b1Y|?QzmWhbM#bDF9U>7TgyugQ0`3lbU8a)}(TlsY1vJKUUpOh5+8Af zx#Q7ckAoY3(vDw6ksqw>Xghr|uo;n?ST zo4W{OwhtoFsQKe|Hhyo^nDU|e#%=w~jwlG1-#Y5*Sd0bfJCxIn_~V_s$nndg`EciY zocfj_)uZLjufan5_(>jponW_Amg?6t3@|=YIghG7Hm|H%UR0lP~z=y#PGw%vjmVKMx>jxKR*)oCx8jMJs746v=0Q}+ih~ZMdA_y@t z;@VhD0zQ+8<_Z2vEbVaHHpL?WcqqNK9Ja3*B!xdFdU29L!<2eMGApqurj;W8CV!lG z?-EZ>PBFaVy?eMuBn2$n=UaC&RAKD8qZ;92{&F`#r^>2Hq#$vIeFoAJ56V12TNh1JU7)*Efb1TzfJLztFg{5{9vq~C%$SDW%+2N6k6HTi&(QV z!BVF^%RRmtliFc$P%_j5cjU|C*xpwPsV%v;LSJRVldBiYU45%D3F*}yoonv+iP}u* z)rY0faYwU9XfhKx54pdQxmb zKP&`vU8=&OoabkK7#wibQB(epr%GVsc8G-&OC}s$`YYAAyAo>{COyP}osU`0}0LBJmRR910 diff --git a/examples/tutorial_1_WaveBot.ipynb b/examples/tutorial_1_WaveBot.ipynb index 34c0cd99..01dcc0b9 100644 --- a/examples/tutorial_1_WaveBot.ipynb +++ b/examples/tutorial_1_WaveBot.ipynb @@ -97,25 +97,6 @@ "ndof = fb.nb_dofs" ] }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "Next we will add the mass and hydrostatic stiffness properties. \n", - "If these values are known they can be added directly.\n", - "Here we will use the fact that the WaveBot is free floating and assume constant density to calculate these properties, which Capytaine can natively perform with the `FloatingBody` created above. For convenience, this functionality has been wrapped in `wecopttool.hydrostatics`." - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "stiffness = wot.hydrostatics.stiffness_matrix(fb).values\n", - "mass = wot.hydrostatics.inertia_matrix(fb).values" - ] - }, { "cell_type": "markdown", "metadata": {}, @@ -145,7 +126,9 @@ "\n", "We will analyze 50 frequencies with a spacing of 0.05 Hz. These frequencies will be used for the Fourier representation of both the wave and the desired PTO force in the pseudo-spectral problem. See the Theory section of the Documentation for more details on the pseudo-spectral problem formulation.\n", "\n", - "If you would like to save our BEM data to a NetCDF file for future use, see the `wecopttool.write_netcdf` function." + "If you would like to save our BEM data to a NetCDF file for future use, see the `wecopttool.write_netcdf` function.\n", + "\n", + "The run_bem() function now calculates the hydrostatics and stores them as part of bem_data. The inertia and stiffness can still be defined manually as part of the `FloatingBody` if desired." ] }, { @@ -236,8 +219,6 @@ "source": [ "wec = wot.WEC.from_bem(\n", " bem_data,\n", - " inertia_matrix=mass,\n", - " hydrostatic_stiffness=stiffness,\n", " constraints=constraints,\n", " friction=None,\n", " f_add=f_add,\n", @@ -512,8 +493,6 @@ "# Update WEC\n", "\n", "wec_2 = wot.WEC.from_bem(bem_data,\n", - " inertia_matrix=mass,\n", - " hydrostatic_stiffness=stiffness,\n", " constraints=constraints_2,\n", " friction=None,\n", " f_add=f_add_2\n", @@ -558,8 +537,6 @@ "source": [ "wec_2_nocon = wot.WEC.from_bem(\n", " bem_data,\n", - " inertia_matrix=mass,\n", - " hydrostatic_stiffness=stiffness,\n", " constraints=None,\n", " friction=None,\n", " f_add=f_add_2)\n", @@ -712,10 +689,6 @@ " nfreq = 50\n", " bem_data = wot.run_bem(fb, freq)\n", "\n", - " # Mass & hydrostatic stiffness\n", - " stiffness_3 = wot.hydrostatics.stiffness_matrix(fb).values\n", - " mass_3 = wot.hydrostatics.inertia_matrix(fb).values\n", - "\n", " # Impedance definition\n", " omega = bem_data.omega.values\n", " gear_ratio = 12.0\n", @@ -764,8 +737,6 @@ "\n", " # Create WEC\n", " wec = wot.WEC.from_bem(bem_data,\n", - " inertia_matrix=mass_3,\n", - " hydrostatic_stiffness=stiffness_3,\n", " constraints=constraints,\n", " friction=None, \n", " f_add=f_add,\n", @@ -874,7 +845,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.10.10" + "version": "3.10.9" }, "vscode": { "interpreter": { diff --git a/examples/tutorial_2_AquaHarmonics.ipynb b/examples/tutorial_2_AquaHarmonics.ipynb index dc0c1425..a1548f14 100644 --- a/examples/tutorial_2_AquaHarmonics.ipynb +++ b/examples/tutorial_2_AquaHarmonics.ipynb @@ -1,7 +1,6 @@ { "cells": [ { - "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ @@ -29,12 +28,15 @@ "import matplotlib.pyplot as plt\n", "from matplotlib import cm\n", "from scipy.optimize import brute\n", + "import xarray as xr\n", + "import logging\n", + "logging.basicConfig()\n", + "logging.getLogger().setLevel(logging.DEBUG)\n", "\n", "import wecopttool as wot" ] }, { - "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ @@ -62,7 +64,6 @@ ] }, { - "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ @@ -83,14 +84,12 @@ ] }, { - "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "#### Hydrostatics and mass\n", "The AquaHarmonics device is positively buoyant (i.e., the buoyancy force is greater than the force due to gravity).\n", - "Therefore, we'll calculate the hydrostatic stiffness from the mesh, but set the rigid-body mass manually. \n", - "We will also calculate the displaced volume and mass from the mesh, which we will need later for defining forces and constraints." + "Therefore, we'll set the rigid-body mass manually, but allow the hydrostatic stiffness to be set automatically by run_bem() based on the mesh. We will also calculate the displaced volume and mass from the mesh (before manually defining the mass of the FloatingBody), which we will need later for defining forces and constraints." ] }, { @@ -99,17 +98,17 @@ "metadata": {}, "outputs": [], "source": [ - "mass = np.atleast_2d(5e3) # [kg]\n", - "stiffness = wot.hydrostatics.stiffness_matrix(fb).values\n", - "\n", "g = 9.81\n", "rho = 1025\n", - "displaced_mass = wot.hydrostatics.inertia_matrix(fb, rho).values # [kg]\n", - "displacement = displaced_mass/rho # [m^3] " + "fb.center_of_mass = [0, 0, 0]\n", + "fb.rotation_center = fb.center_of_mass\n", + "displaced_mass = fb.compute_rigid_body_inertia(rho=rho).values # [kg]\n", + "displacement = displaced_mass/rho # [m^3] \n", + "\n", + "fb.mass = np.atleast_2d(5e3) # [kg]" ] }, { - "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ @@ -134,7 +133,6 @@ ] }, { - "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ @@ -167,7 +165,6 @@ ] }, { - "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ @@ -233,7 +230,6 @@ ] }, { - "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ @@ -257,7 +253,6 @@ ] }, { - "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ @@ -281,7 +276,6 @@ ] }, { - "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ @@ -321,7 +315,6 @@ ] }, { - "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ @@ -349,7 +342,6 @@ ] }, { - "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ @@ -410,7 +402,6 @@ ] }, { - "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ @@ -474,15 +465,13 @@ ] }, { - "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "#### WEC object\n", "Finally, we can use all the different components we've developed thus far to construct a `WEC` object:\n", "\n", - " * **BEM data** - defines the hydrodynamics of the hull\n", - " * **Inertia and hydrostatics** - rigid-body inertia and hydrostatic stiffness\n", + " * **BEM data** - defines the hydrodynamics of the hull and includes hydrostatics\n", " * **Constraints** - limitations on the hardware (max power, max torque, etc.) and the constraint to prevent the tether from going slack\n", " * **Additional forces** - this captures all of the forces on the WEC that are not due to the interaction between the hull and water (PTO, etc.)" ] @@ -490,20 +479,19 @@ { "cell_type": "code", "execution_count": null, - "metadata": {}, + "metadata": { + "scrolled": false + }, "outputs": [], "source": [ "wec = wot.WEC.from_bem(\n", " bem_data,\n", - " inertia_matrix = mass,\n", - " hydrostatic_stiffness = stiffness,\n", " constraints = constraints,\n", " f_add = f_add,\n", ")" ] }, { - "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ @@ -526,7 +514,6 @@ ] }, { - "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ @@ -544,7 +531,6 @@ ] }, { - "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ @@ -558,7 +544,9 @@ { "cell_type": "code", "execution_count": null, - "metadata": {}, + "metadata": { + "scrolled": true + }, "outputs": [], "source": [ "scale_x_wec = 1e1\n", @@ -583,7 +571,6 @@ ] }, { - "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ @@ -603,7 +590,6 @@ ] }, { - "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ @@ -621,7 +607,6 @@ ] }, { - "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ @@ -719,7 +704,6 @@ ] }, { - "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ @@ -747,7 +731,6 @@ ] }, { - "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ @@ -783,12 +766,11 @@ " # Unpack geometry variables\n", " mass_ratio = x[0]\n", " mass_var = mass_ratio * max_mass\n", - "\n", + " bem_data['inertia_matrix'] = mass_var\n", + " \n", " # update WEC \n", " wec_mass = wot.WEC.from_bem(\n", " bem_data,\n", - " inertia_matrix = mass_var,\n", - " hydrostatic_stiffness = stiffness,\n", " constraints = constraints,\n", " friction = None,\n", " f_add = f_add,\n", @@ -819,7 +801,6 @@ ] }, { - "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ @@ -853,7 +834,6 @@ ] }, { - "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ @@ -912,6 +892,13 @@ "ax3.set_yticks([1, 5, 10, 15, 20, 25])\n", "ax3.legend(loc=9)" ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] } ], "metadata": { diff --git a/examples/tutorial_3_LUPA.ipynb b/examples/tutorial_3_LUPA.ipynb index 52a89035..a57602dc 100644 --- a/examples/tutorial_3_LUPA.ipynb +++ b/examples/tutorial_3_LUPA.ipynb @@ -232,7 +232,7 @@ "source": [ "#### Combined `FloatingBody`\n", "\n", - "With both WEC bodies defined separately, we can now create a union of the bodies and define the properties for the overall LUPA device. At the equilibrium position the float is neutrally buoyant while the spar is positively buoyant and requires mooring pre-tension. The combined center of mass and moment of inertia can be found by using the given values weighted by the mass (via the parallel axis theorem for the moment of inertia). This is also when we can specify the surge and pitch degrees of freedom.\n", + "With both WEC bodies defined separately, we should first define the respective centers of mass and rotation centers for the `FloatingBody`s. Then, we can create a union of the bodies and define the properties for the overall LUPA device. At the equilibrium position the float is neutrally buoyant while the spar is positively buoyant and requires mooring pre-tension. The combined center of mass and moment of inertia can be found by using the given values weighted by the mass (via the parallel axis theorem for the moment of inertia). This is also when we can specify the surge and pitch degrees of freedom.\n", "\n", "We are using the density of *fresh* water, $\\rho = 1000 kg/m^3$, since we are modeling LUPA in a wave flume." ] @@ -246,19 +246,23 @@ "# density of fresh water\n", "rho = 1000\n", "\n", - "# floating body\n", - "lupa_fb = float_fb + spar_fb\n", - "lupa_fb.name = 'LUPA'\n", - "\n", "# mass properties float\n", "mass_float = float_mass_properties['mass'] \n", "cm_float = np.array(float_mass_properties['CG'])\n", "pitch_inertia_float = float_mass_properties['MOI'][1] \n", + "float_fb.center_of_mass = cm_float\n", + "float_fb.rotation_center = float_fb.center_of_mass\n", "\n", "# mass properties spar\n", "mass_spar = spar_mass_properties['mass'] \n", "cm_spar = np.array(spar_mass_properties['CG'])\n", "pitch_inertia_spar = spar_mass_properties['MOI'][1] \n", + "spar_fb.center_of_mass = cm_spar\n", + "spar_fb.rotation_center = spar_fb.center_of_mass\n", + "\n", + "# floating body\n", + "lupa_fb = float_fb + spar_fb\n", + "lupa_fb.name = 'LUPA'\n", "\n", " # mass properties LUPA\n", "lupa_fb.center_of_mass = ((mass_float*cm_float + mass_spar*cm_spar)\n", @@ -279,6 +283,36 @@ "lupa_fb.add_rotation_dof(name='Pitch')" ] }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "#### Define Hydrostatics Manually\n", + "\n", + "When combining multiple bodies, we need to be careful to set up the hydrostatics correctly. The bodies move separately in heave but move together in surge and pitch. Therefore, we should define the individual heave inertia values for each body and define the total surge and pitch inertia values. The inertia then needs to be reformatted as a DataArray for Capytaine. Lastly, the hydrostatic stiffness can be calculated for the total immersed body." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# reorganize inertia values into DataArray for Capytaine\n", + "rigid_inertia_matrix_xr = xr.DataArray(data=np.asarray((inertia)),\n", + " dims=['influenced_dof', 'radiating_dof'],\n", + " coords={'influenced_dof': list(lupa_fb.dofs),\n", + " 'radiating_dof': list(lupa_fb.dofs)},\n", + " name=\"inertia_matrix\")\n", + "\n", + "# Set FloatingBody inertia matrix\n", + "lupa_fb.inertia_matrix = rigid_inertia_matrix_xr\n", + "\n", + "# Calculate hydrostatic stiffness after keeping immersed value\n", + "lupa_fb = lupa_fb.keep_immersed_part()\n", + "lupa_fb.hydrostatic_stiffness = lupa_fb.compute_hydrostatic_stiffness(rho=rho)" + ] + }, { "cell_type": "markdown", "metadata": {}, @@ -329,12 +363,8 @@ "try:\n", " bem_data = wot.read_netcdf(filename)\n", "except:\n", - " bem_data = wot.run_bem(lupa_fb, freq)\n", - " wot.write_netcdf(filename, bem_data)\n", - "\n", - "# hydrostatics\n", - "_ = lupa_fb.compute_hydrostatics(rho=rho)\n", - "hydrostatic_stiffness = lupa_fb.hydrostatic_stiffness" + " bem_data = wot.run_bem(lupa_fb, freq, rho=rho)\n", + " wot.write_netcdf(filename, bem_data)" ] }, { @@ -910,8 +940,6 @@ "\n", "# WEC\n", "wec = wot.WEC.from_bem(bem_data,\n", - " inertia_matrix=inertia,\n", - " hydrostatic_stiffness=hydrostatic_stiffness,\n", " constraints=constraints,\n", " friction=friction,\n", " f_add=f_add,\n", @@ -1296,8 +1324,6 @@ "\n", " # create WEC object\n", " wec = wot.WEC.from_bem(bem_data,\n", - " inertia_matrix=inertia,\n", - " hydrostatic_stiffness=hydrostatic_stiffness,\n", " constraints=constraints,\n", " friction=friction,\n", " f_add=f_add,\n", diff --git a/examples/tutorial_4_Pioneer.ipynb b/examples/tutorial_4_Pioneer.ipynb index d0a8fde6..a87d36b6 100644 --- a/examples/tutorial_4_Pioneer.ipynb +++ b/examples/tutorial_4_Pioneer.ipynb @@ -1,7 +1,6 @@ { "cells": [ { - "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ @@ -34,6 +33,7 @@ "import autograd.numpy as np\n", "import matplotlib.pyplot as plt\n", "from scipy.linalg import block_diag\n", + "import xarray as xr\n", "\n", "import wecopttool as wot\n", "\n", @@ -42,7 +42,6 @@ ] }, { - "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ @@ -85,7 +84,6 @@ ] }, { - "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ @@ -108,7 +106,6 @@ ] }, { - "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ @@ -137,7 +134,6 @@ ] }, { - "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ @@ -178,7 +174,6 @@ ] }, { - "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ @@ -197,12 +192,11 @@ ] }, { - "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "#### Hydrodynamics and hydrostatics\n", - "As mentioned above, the `FloatingBody` object in Capytaine only needs to model the buoy, since no other components are being excited by the waves." + "As mentioned above, the `FloatingBody` object in Capytaine only needs to model the buoy, since no other components are being excited by the waves. The inertia matrix needs to be defined manually here since it is not based on the displaced mass." ] }, { @@ -216,29 +210,36 @@ "pnr_fb.center_of_mass = np.array([0., 0., buoy_props['CG']])\n", "pnr_fb.rotation_center = pnr_fb.center_of_mass\n", "ndof = pnr_fb.nb_dofs\n", - "pnr_fb.show_matplotlib()" + "pnr_fb.show_matplotlib()\n", + "\n", + "pnr_fb.inertia_matrix = xr.DataArray([[buoy_props['MOI']]],\n", + " dims=['influenced_dof', 'radiating_dof'],\n", + " coords={'influenced_dof' : ['Pitch'],\n", + " 'radiating_dof' : ['Pitch']},\n", + " name=\"inertia_matrix\"\n", + " )" ] }, { "cell_type": "code", "execution_count": null, - "metadata": {}, + "metadata": { + "scrolled": true + }, "outputs": [], "source": [ "rho = 1025. # kg/m^3\n", "freq = wot.frequency(f1, nfreq, False) # False -> no zero frequency\n", - "bem_data = wot.run_bem(pnr_fb, freq)\n", + "bem_data = wot.run_bem(pnr_fb, freq, rho=rho)\n", "omega = bem_data.omega.values\n", "\n", "pnr_fb.keep_immersed_part()\n", - "k_buoy = pnr_fb.compute_hydrostatic_stiffness(rho=rho).values.squeeze()\n", "k_spring = spring_props['Max torque'] / spring_props['Max displacement']\n", - "print(f'Hydrostatic stiffness from Capytaine: {k_buoy} N-m/rad')\n", + "print(f\"Hydrostatic stiffness from Capytaine: {bem_data['hydrostatic_stiffness'].values.squeeze()} N-m/rad\")\n", "print('Hydrostatic stiffness from experiment: 37204 N-m/rad')" ] }, { - "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ @@ -280,7 +281,6 @@ ] }, { - "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ @@ -292,7 +292,6 @@ ] }, { - "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ @@ -316,7 +315,6 @@ ] }, { - "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ @@ -351,7 +349,6 @@ ] }, { - "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ @@ -372,7 +369,6 @@ ] }, { - "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ @@ -397,7 +393,6 @@ ] }, { - "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ @@ -439,7 +434,6 @@ ] }, { - "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ @@ -484,7 +478,6 @@ ] }, { - "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ @@ -506,7 +499,6 @@ ] }, { - "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ @@ -537,7 +529,6 @@ ] }, { - "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ @@ -569,7 +560,6 @@ ] }, { - "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ @@ -595,7 +585,6 @@ ] }, { - "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ @@ -610,8 +599,6 @@ "outputs": [], "source": [ "wec = wot.WEC.from_bem(bem_data,\n", - " inertia_matrix=np.array([[buoy_props['MOI']]]),\n", - " hydrostatic_stiffness=k_buoy,\n", " f_add=f_add,\n", " constraints=constraints,\n", " uniform_shift=False,\n", @@ -619,7 +606,6 @@ ] }, { - "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ @@ -647,7 +633,6 @@ ] }, { - "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ @@ -678,7 +663,6 @@ ] }, { - "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ @@ -745,7 +729,6 @@ ] }, { - "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ @@ -753,7 +736,6 @@ ] }, { - "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ @@ -767,7 +749,7 @@ "metadata": {}, "outputs": [], "source": [ - "hydro_data = wot.linear_hydrodynamics(bem_data, np.array([[buoy_props['MOI']]]), k_buoy, friction=None) \n", + "hydro_data = wot.add_linear_friction(bem_data, friction=None) \n", "hydro_data = wot.check_linear_damping(hydro_data, uniform_shift=False) \n", "\n", "Zi = wot.hydrodynamic_impedance(hydro_data)\n", @@ -784,7 +766,6 @@ ] }, { - "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ @@ -856,7 +837,6 @@ ] }, { - "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ @@ -886,7 +866,6 @@ ] }, { - "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ @@ -968,11 +947,18 @@ " axi.label_outer()\n", " axi.autoscale(axis='x', tight=True)" ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] } ], "metadata": { "kernelspec": { - "display_name": "conda-wot", + "display_name": "wot_dev", "language": "python", "name": "python3" }, @@ -986,7 +972,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.10.10" + "version": "3.10.9" } }, "nbformat": 4, diff --git a/tests/test_core.py b/tests/test_core.py index 6695357d..f372d760 100644 --- a/tests/test_core.py +++ b/tests/test_core.py @@ -55,17 +55,22 @@ def bem_data(f1, nfreq): ndof = 2; ndir = 3; radiation_dims = ['radiating_dof', 'influenced_dof', 'omega'] excitation_dims = ['influenced_dof', 'wave_direction', 'omega'] + hydrostatics_dims = ['radiating_dof', 'influenced_dof'] added_mass = np.ones([ndof, ndof, nfreq]) radiation_damping = np.ones([ndof, ndof, nfreq]) diffraction_force = np.ones([ndof, ndir, nfreq], dtype=complex) + 1j Froude_Krylov_force = np.ones([ndof, ndir, nfreq], dtype=complex) + 1j + inertia_matrix = np.ones([ndof, ndof]) + hydrostatic_stiffness = np.ones([ndof, ndof]) data_vars = { 'added_mass': (radiation_dims, added_mass), 'radiation_damping': (radiation_dims, radiation_damping), 'diffraction_force': (excitation_dims, diffraction_force), - 'Froude_Krylov_force': (excitation_dims, Froude_Krylov_force) + 'Froude_Krylov_force': (excitation_dims, Froude_Krylov_force), + 'inertia_matrix': (hydrostatics_dims, inertia_matrix), + 'hydrostatic_stiffness': (hydrostatics_dims, hydrostatic_stiffness) } return xr.Dataset(data_vars=data_vars, coords=coords) @@ -75,11 +80,9 @@ def hydro_data(bem_data): """Synthetic hydro-data containing inertia, stiffness, and friction in addition to the coefficients in `bem_data`.""" ndof = len(bem_data.influenced_dof) - inertia_matrix = np.ones([ndof, ndof]) - stiffness = np.ones([ndof, ndof]) friction = np.ones([ndof, ndof]) - data = wot.linear_hydrodynamics( - bem_data, inertia_matrix, stiffness, friction + data = wot.add_linear_friction( + bem_data, friction ) return data @@ -1300,12 +1303,12 @@ def test_round_trip(self, bem_data, data_back): class TestLinearHydrodynamics: - """Test function :python:`linear_hydrodynamics`.""" + """Test function :python:`add_linear_friction`.""" def test_values(self, bem_data, hydro_data): """Test the function returns expected values.""" mat = np.array([[1, 1], [1, 1]]) - calculated = wot.linear_hydrodynamics(bem_data, mat, mat, mat) + calculated = wot.add_linear_friction(bem_data, mat) xr.testing.assert_allclose(calculated, hydro_data) diff --git a/tests/test_hydrostatics.py b/tests/test_hydrostatics.py deleted file mode 100644 index d4178cf3..00000000 --- a/tests/test_hydrostatics.py +++ /dev/null @@ -1,204 +0,0 @@ -""" Unit tests for functions in the :python:`hydrostatics.py` module. -""" - -import pytest -import capytaine as cpy -import numpy as np - -import wecopttool as wot -from wecopttool.core import _default_parameters - - -# setup simple constant density rectangular barge -@pytest.fixture() -def rho(): - """Water density [kg/m^3].""" - return _default_parameters['rho'] - - -@pytest.fixture() -def g(): - """Gravitational acceleration [m/s^2].""" - return _default_parameters['g'] - - -@pytest.fixture() -def lx(): - """Barge length [m].""" - return 4.0 - - -@pytest.fixture() -def ly(): - """Barge width [m].""" - return 5.0 - - -@pytest.fixture() -def lz(): - """Barge height [m].""" - return 2.0 - - -@pytest.fixture() -def cog(): - """Center of gravity of the barge.""" - return (0.0, 0.0, 0.0) - - -@pytest.fixture() -def mass(lx, ly, lz, rho): - """Mass of the barge.""" - return lx*ly*lz/2 * rho - - -@pytest.fixture() -def fb(lx, ly, lz): - """Simple constant density rectangular barge.""" - rect = cpy.RectangularParallelepiped( - (lx, ly, lz), resolution=(100, 100, 10), center=(0.0, 0.0, 0.0,)) - rect.add_all_rigid_body_dofs() - rect.center_of_mass = None - fb.mass = None - return rect - - -@pytest.fixture() -def fb_mass(fb, cog, mass): - """Simple constant density rectangular barge with the mass - properties specified. - """ - fbc = fb.copy() - fbc.center_of_mass = cog - fbc.mass = mass - return fbc - - -@pytest.fixture() -def stiffness(lx, ly, lz, rho, g): - """True (theoretical/calculated) value of the stiffness matrix for - the barge. - """ - stiffness = np.zeros([6, 6]) - stiffness[2, 2] = rho*g * lx*ly - stiffness[2, 3] = 0.0 - stiffness[2, 4] = 0.0 - stiffness[3, 3] = (rho*g * lx*ly**3 /12) + rho*g*(lx*ly*lz/2)*(-lz/4) - stiffness[3, 4] = 0.0 - stiffness[3, 5] = 0.0 - stiffness[4, 4] = (rho*g * ly*lx**3 /12) + rho*g*(lx*ly*lz/2)*(-lz/4) - stiffness[4, 5] = 0.0 - stiffness[3, 2] = stiffness[2, 3] - stiffness[4, 2] = stiffness[2, 4] - stiffness[4, 3] = stiffness[3, 4] - return stiffness - - -@pytest.fixture() -def inertia(lx, ly, lz, mass): - """True (theoretical/calculated) value of the inertia matrix for - the barge. - """ - inertia = np.zeros([6, 6]) - inertia[0, 0] = mass - inertia[1, 1] = mass - inertia[2, 2] = mass - inertia[3, 3] = mass/12 * (ly**2 + lz**2) - inertia[3, 4] = 0.0 - inertia[3, 5] = 0.0 - inertia[4, 3] = 0.0 - inertia[4, 4] = mass/12 * (lx**2 + lz**2) - inertia[4, 5] = 0.0 - inertia[5, 3] = 0.0 - inertia[5, 4] = 0.0 - inertia[5, 5] = mass/12 * (lx**2 + ly**2) - return inertia - - -class TestStiffnessMatrix: - """Test function :python:`hydrostatics.stiffness_matrix`.""" - - def test_inferred_cog(self, fb, rho, g, stiffness): - """Test the function with the center of gravity not provided, - but inferred. - """ - stiffness_calc = wot.hydrostatics.stiffness_matrix(fb, rho, g) - assert np.allclose(stiffness, stiffness_calc, rtol=0.01) - - - def test_given_cog(self, fb, rho, g, cog, stiffness): - """Test the function with the center of gravity provided as a - function input and not in the floating body. - """ - stiffness_calc = wot.hydrostatics.stiffness_matrix(fb, rho, g, cog) - assert np.allclose(stiffness, stiffness_calc, rtol=0.01) - - - def test_given_cog_fb(self, fb_mass, rho, g, stiffness): - """Test the function with the center of gravity provided in the - floating body and not as a function input. - """ - stiffness_calc = wot.hydrostatics.stiffness_matrix(fb_mass, rho, g) - assert np.allclose(stiffness, stiffness_calc, rtol=0.01) - - - def test_given_cog_redundant(self, fb_mass, rho, g, cog, stiffness): - """Test the function with the center of gravity provided in both - the floating body and the function inputs. - """ - stiffness_calc = wot.hydrostatics.stiffness_matrix(fb_mass, rho, g, cog) - assert np.allclose(stiffness, stiffness_calc, rtol=0.01) - - - def test_cog_mismatch(self, fb_mass, rho, g): - """Test that the function fails if the center of gravity is - provided in both the floating body and the function inputs but - have different values. - """ - cog_wrong = (0, 0, -0.1) - with pytest.raises(ValueError): - wot.hydrostatics.stiffness_matrix(fb_mass, rho, g, cog_wrong) - - -class TestInertiaMatrix: - """Test function :python:`hydrostatics.test_inertia_matrix`.""" - - def test_inferred_mass(self, fb, rho, inertia): - """Test the function with the mass not provided, but inferred. - """ - inertia_calc = wot.hydrostatics.inertia_matrix(fb, rho) - assert np.allclose(inertia, inertia_calc, rtol=0.01) - - - def test_given_mass(self, fb, rho, cog, mass, inertia): - """Test the function with the mass provided as a function input - and not in the floating body. - """ - inertia_calc = wot.hydrostatics.inertia_matrix(fb, rho, cog, mass) - assert np.allclose(inertia, inertia_calc, rtol=0.01) - - - def test_given_mass_fb(self, fb_mass, rho, cog, mass, inertia): - """Test the function with the mass provided in the floating body - and not as a function input. - """ - inertia_calc = wot.hydrostatics.inertia_matrix(fb_mass, rho) - assert np.allclose(inertia, inertia_calc, rtol=0.01) - - - def test_given_mass_redundant(self, fb_mass, rho, cog, mass, inertia): - """Test the function with the mass provided in both the floating - body and the function inputs. - """ - inertia_calc = wot.hydrostatics.inertia_matrix(fb_mass, rho, cog, mass) - assert np.allclose(inertia, inertia_calc, rtol=0.01) - - - def test_mass_mismatch(self, fb_mass, rho, cog, mass): - """Test that the function fails if the mass is provided in both - the floating body and the function inputs but have different - values. - """ - mass_wrong = 0.9*mass - with pytest.raises(ValueError): - wot.hydrostatics.inertia_matrix(fb_mass, rho, cog, mass_wrong) diff --git a/tests/test_integration.py b/tests/test_integration.py index 46ed235e..a93c178e 100644 --- a/tests/test_integration.py +++ b/tests/test_integration.py @@ -104,20 +104,16 @@ def irregular_wave(f1, nfreq): @pytest.fixture(scope='module') def wec_from_bem(f1, nfreq, bem, fb, pto): """Simple WEC: 1 DOF, no constraints.""" - mass = wot.hydrostatics.inertia_matrix(fb).values - hstiff = wot.hydrostatics.stiffness_matrix(fb).values f_add = {"PTO": pto.force_on_wec} - wec = wot.WEC.from_bem(bem, mass, hstiff, f_add=f_add) + wec = wot.WEC.from_bem(bem, f_add=f_add) return wec @pytest.fixture(scope='module') def wec_from_floatingbody(f1, nfreq, fb, pto): """Simple WEC: 1 DOF, no constraints.""" - mass = wot.hydrostatics.inertia_matrix(fb).values - hstiff = wot.hydrostatics.stiffness_matrix(fb).values f_add = {"PTO": pto.force_on_wec} - wec = wot.WEC.from_floating_body(fb, f1, nfreq, mass, hstiff, f_add=f_add) + wec = wot.WEC.from_floating_body(fb, f1, nfreq, f_add=f_add) return wec @@ -130,10 +126,13 @@ def wec_from_impedance(bem, pto, fb): w = np.expand_dims(omega, [0, 1]) A = bemc['added_mass'].values B = bemc['radiation_damping'].values - mass = wot.hydrostatics.inertia_matrix(fb).values - hstiff = wot.hydrostatics.stiffness_matrix(fb).values + fb.center_of_mass = [0, 0, 0] + fb.rotation_center = fb.center_of_mass + fb = fb.copy(name=f"{fb.name}_immersed").keep_immersed_part() + mass = bemc['inertia_matrix'].values + hstiff = bemc['hydrostatic_stiffness'].values K = np.expand_dims(hstiff, 2) - + freqs = omega / (2 * np.pi) impedance = (A + mass)*(1j*w) + B + K/(1j*w) exc_coeff = bem['Froude_Krylov_force'] + bem['diffraction_force'] @@ -146,9 +145,7 @@ def wec_from_impedance(bem, pto, fb): @pytest.fixture(scope='module') def resonant_wave(f1, nfreq, fb, bem): """Regular wave at natural frequency of the WEC""" - mass = wot.hydrostatics.inertia_matrix(fb).values - hstiff = wot.hydrostatics.stiffness_matrix(fb).values - hd = wot.linear_hydrodynamics(bem, mass, hstiff) + hd = wot.add_linear_friction(bem) Zi = wot.hydrodynamic_impedance(hd) wn = Zi['omega'][np.abs(Zi).argmin(dim='omega')].item() waves = wot.waves.regular_wave(f1, nfreq, freq=wn/2/np.pi, amplitude=0.1) @@ -215,7 +212,7 @@ def test_same_wec_init(wec_from_bem, bem_res = wec_from_bem.residual(x_wec_0, x_opt_0, waves) fb_res = wec_from_floatingbody.residual(x_wec_0, x_opt_0, waves) imp_res = wec_from_impedance.residual(x_wec_0, x_opt_0, waves) - + assert fb_res == approx(bem_res, rel=0.01) assert imp_res == approx(bem_res, rel=0.01) @@ -227,17 +224,17 @@ class TestTheoreticalPowerLimits: @pytest.fixture(scope='class') def mass(self, fb): """Rigid-body mass""" - return wot.hydrostatics.inertia_matrix(fb).values + return fb.compute_rigid_body_inertia() @pytest.fixture(scope='class') def hstiff(self, fb): """Hydrostatic stiffness""" - return wot.hydrostatics.stiffness_matrix(fb).values + return fb.compute_hydrostatic_stiffness() @pytest.fixture(scope='class') - def hydro_impedance(self, bem, mass, hstiff): + def hydro_impedance(self, bem): """Intrinsic hydrodynamic impedance""" - hd = wot.linear_hydrodynamics(bem, mass, hstiff) + hd = wot.add_linear_friction(bem) hd = wot.check_linear_damping(hd) Zi = wot.hydrodynamic_impedance(hd) return Zi @@ -246,14 +243,12 @@ def test_p_controller_resonant_wave(self, bem, resonant_wave, p_controller_pto, - mass, - hstiff, hydro_impedance): """Proportional controller should match optimum for natural resonant wave""" - + f_add = {"PTO": p_controller_pto.force_on_wec} - wec = wot.WEC.from_bem(bem, mass, hstiff, f_add=f_add) + wec = wot.WEC.from_bem(bem, f_add=f_add) res = wec.solve(waves=resonant_wave, obj_fun=p_controller_pto.average_power, @@ -275,19 +270,17 @@ def test_p_controller_resonant_wave(self, power_optimal = (np.abs(Fex)**2/8 / np.real(hydro_impedance.squeeze()) ).squeeze().sum('omega').item() - assert power_sol == approx(power_optimal, rel=0.02) + assert power_sol == approx(power_optimal, rel=0.03) def test_pi_controller_regular_wave(self, bem, regular_wave, pi_controller_pto, - mass, - hstiff, hydro_impedance): """PI controller matches optimal for any regular wave""" f_add = {"PTO": pi_controller_pto.force_on_wec} - wec = wot.WEC.from_bem(bem, mass, hstiff, f_add=f_add) + wec = wot.WEC.from_bem(bem, f_add=f_add) res = wec.solve(waves=regular_wave, obj_fun=pi_controller_pto.average_power, @@ -317,14 +310,12 @@ def test_unstructured_controller_irregular_wave(self, regular_wave, pto, nfreq, - mass, - hstiff, hydro_impedance): """Unstructured (numerical optimal) controller matches optimal for any irregular wave when unconstrained""" f_add = {"PTO": pto.force_on_wec} - wec = wot.WEC.from_bem(bem, mass, hstiff, f_add=f_add) + wec = wot.WEC.from_bem(bem, f_add=f_add) res = wec.solve(waves=regular_wave, obj_fun=pto.average_power, @@ -343,4 +334,4 @@ def test_unstructured_controller_irregular_wave(self, power_optimal = (np.abs(Fex)**2/8 / np.real(hydro_impedance.squeeze()) ).squeeze().sum('omega').item() - assert power_sol == approx(power_optimal, rel=1e-3) + assert power_sol == approx(power_optimal, rel=1e-2) diff --git a/wecopttool/__init__.py b/wecopttool/__init__.py index d36607f3..9b3d6b66 100644 --- a/wecopttool/__init__.py +++ b/wecopttool/__init__.py @@ -31,7 +31,6 @@ from wecopttool.core import * from wecopttool import waves -from wecopttool import hydrostatics from wecopttool import pto try: diff --git a/wecopttool/core.py b/wecopttool/core.py index e60f14a2..1ea44943 100644 --- a/wecopttool/core.py +++ b/wecopttool/core.py @@ -40,7 +40,7 @@ "standard_forces", "run_bem", "change_bem_convention", - "linear_hydrodynamics", + "add_linear_friction", "wave_excitation", "hydrodynamic_impedance", "atleast_2d", @@ -50,6 +50,7 @@ "decompose_state", "frequency_parameters", "time_results", + "set_fb_centers", ] @@ -284,8 +285,6 @@ def __repr__(self) -> str: @staticmethod def from_bem( bem_data: Union[Dataset, Union[str, Path]], - inertia_matrix: Optional[ndarray] = None, - hydrostatic_stiffness: Optional[ndarray] = None, friction: Optional[ndarray] = None, f_add: Optional[TIForceDict] = None, constraints: Optional[Iterable[Mapping]] = None, @@ -312,26 +311,16 @@ def from_bem( :py:func:`wecopttool.run_bem`, rather than running Capytaine directly, which outputs the results in the correct convention. The results can be saved - using :py:func:`wecopttool.write_netcdf`. - - In addition to the Capytaine results, if the dataset contains - the :python:`inertia_matrix`, :python:`hydrostatic_stiffness`, - or :python:`friction` these do not need to be provided - separately. + using :py:func:`wecopttool.write_netcdf`. + :py:func:`wecopttool.run_bem` also computes the inertia and + hydrostatic stiffness which should be included in bem_data. Parameters ---------- bem_data Linear hydrodynamic coefficients obtained using the boundary element method (BEM) code Capytaine, with sign convention - corrected. - inertia_matrix - Inertia matrix of size :python:`(ndof, ndof)`. - :python:`None` if included in :python:`bem_data`. - hydrostatic_stiffness - Linear hydrostatic restoring coefficient of size - :python:`(nodf, ndof)`. - :python:`None` if included in :python:`bem_data`. + corrected. Also includes inertia and hydrostatic stiffness. friction Linear friction, in addition to radiation damping, of size :python:`(nodf, ndof)`. @@ -362,32 +351,16 @@ def from_bem( If :python:`None` the names :python:`['DOF_0', ..., 'DOF_N']` are used. - Raises - ------ - ValueError - If either :python:`inertia_matrix` or - :python:`hydrostatic_stiffness` are :python:`None` and is - not included in :python:`bem_data`. - See :py:func:`wecopttool.linear_hydrodynamics`. - ValueError - If any of :python:`inertia_matrix`, - :python:`hydrostatic_stiffness`, or :python:`stiffness` are - both provided and included in :python:`bem_data` but have - different values. - See :py:func:`wecopttool.linear_hydrodynamics`. - See Also -------- - run_bem, linear_hydrodynamics, change_bem_convention, + run_bem, add_linear_friction, change_bem_convention, write_netcdf, check_linear_damping """ if isinstance(bem_data, (str, Path)): bem_data = read_netcdf(bem_data) - # add inertia_matrix, hydrostatic stiffness, and friction - hydro_data = linear_hydrodynamics( - bem_data, inertia_matrix, hydrostatic_stiffness, friction) - if inertia_matrix is None: - inertia_matrix = hydro_data['inertia_matrix'].values + # add friction + hydro_data = add_linear_friction(bem_data, friction) + inertia_matrix = hydro_data['inertia_matrix'].values # frequency array f1, nfreq = frequency_parameters( @@ -411,8 +384,6 @@ def from_floating_body( fb: cpy.FloatingBody, f1: float, nfreq: int, - inertia_matrix: ndarray, - hydrostatic_stiffness: ndarray, friction: Optional[ndarray] = None, f_add: Optional[TIForceDict] = None, constraints: Optional[Iterable[Mapping]] = None, @@ -454,11 +425,6 @@ def from_floating_body( nfreq Number of frequencies (not including zero frequency), i.e., :python:`freqs = [0, f1, 2*f1, ..., nfreq*f1]`. - inertia_matrix - Inertia matrix of size :python:`(ndof, ndof)`. - hydrostatic_stiffness - Linear hydrostatic restoring coefficient of size - :python:`(ndof, ndof)`. friction Linear friction, in addition to radiation damping, of size :python:`(ndof, ndof)`. @@ -504,7 +470,7 @@ def from_floating_body( bem_data = run_bem( fb, freq, wave_directions, rho=rho, g=g, depth=depth) wec = WEC.from_bem( - bem_data, inertia_matrix, hydrostatic_stiffness, friction, f_add, + bem_data, friction, f_add, constraints, min_damping=min_damping) return wec @@ -2157,12 +2123,17 @@ def run_bem( # radiation only problem, no diffraction or excitation test_matrix = test_matrix.drop_vars('wave_direction') if write_info is None: - write_info = {'hydrostatics': False, + write_info = {'hydrostatics': True, 'mesh': False, 'wavelength': False, 'wavenumber': False, } wec_im = fb.copy(name=f"{fb.name}_immersed").keep_immersed_part() + wec_im = set_fb_centers(wec_im, rho=rho) + if not hasattr(wec_im, 'inertia_matrix'): + wec_im.inertia_matrix = wec_im.compute_rigid_body_inertia(rho=rho) + if not hasattr(wec_im, 'hydrostatic_stiffness'): + wec_im.hydrostatic_stiffness = wec_im.compute_hydrostatic_stiffness(rho=rho, g=g) bem_data = solver.fill_dataset( test_matrix, wec_im, n_jobs=njobs, **write_info) return change_bem_convention(bem_data) @@ -2189,14 +2160,11 @@ def change_bem_convention(bem_data: Dataset) -> Dataset: return bem_data -def linear_hydrodynamics( +def add_linear_friction( bem_data: Dataset, - inertia_matrix: Optional[ArrayLike] = None, - hydrostatic_stiffness: Optional[ArrayLike] = None, friction: Optional[ArrayLike] = None ) -> Dataset: - """Add rigid body inertia_matrix, hydrostatic stiffness, and linear - friction to BEM data. + """Add linear friction to BEM data. Returns the Dataset with the additional information added. @@ -2205,62 +2173,32 @@ def linear_hydrodynamics( bem_data Linear hydrodynamic coefficients obtained using the boundary element method (BEM) code Capytaine, with sign convention - corrected. - inertia_matrix - Inertia matrix of size :python:`(ndof, ndof)`. - :python:`None` if included in :python:`bem_data`. - hydrostatic_stiffness - Linear hydrostatic restoring coefficient of size - :python:`(nodf, ndof)`. - :python:`None` if included in :python:`bem_data`. + corrected. Also includes inertia and hydrostatic stiffness. friction Linear friction, in addition to radiation damping, of size :python:`(nodf, ndof)`. :python:`None` if included in :python:`bem_data` or to set to zero. - - Raises - ------ - ValueError - If either :python:`inertia_matrix` or - :python:`hydrostatic_stiffness` are :python:`None` and is not - included in :python:`bem_data`. - ValueError - If any of :python:`inertia_matrix`, - :python:`hydrostatic_stiffness`, or :python:`friction` are both - provided and included in :python:`bem_data` but have different - values. """ - vars = {'inertia_matrix': inertia_matrix, 'friction': friction, - 'hydrostatic_stiffness': hydrostatic_stiffness} - dims = ['radiating_dof', 'influenced_dof'] - hydro_data = bem_data.copy(deep=True) - - for name, data in vars.items(): - org = name in hydro_data.variables.keys() - new = data is not None - if new and org: + + if friction is not None: + if 'friction' in hydro_data.variables.keys(): if not np.allclose(data, hydro_data.variables[name]): raise ValueError( - f'BEM data already has variable "{name}" ' + - 'with diferent values') - else : + f'Variable "friction" is already in BEM data ' + + f'with different values.') + else: _log.warning( f'Variable "{name}" is already in BEM data ' + 'with same value.') - elif (not new) and (not org): - if name=='friction': - ndof = len(hydro_data["influenced_dof"]) - hydro_data[name] = (dims, np.zeros([ndof, ndof])) - else: - raise ValueError( - f'Variable "{name}" is not in BEM data and ' + - 'was not provided.') - elif new: - data = atleast_2d(data) - hydro_data[name] = (dims, data) - + else: + data = atleast_2d(friction) + hydro_data['friction'] = (dims, friction) + elif friction is None: + ndof = len(hydro_data["influenced_dof"]) + hydro_data['friction'] = (dims, np.zeros([ndof, ndof])) + return hydro_data @@ -2320,7 +2258,7 @@ def hydrodynamic_impedance(hydro_data: Dataset) -> Dataset: ---------- hydro_data Dataset with linear hydrodynamic coefficients produced by - :py:func:`wecopttool.linear_hydrodynamics`. + :py:func:`wecopttool.add_linear_friction`. """ Zi = (hydro_data['inertia_matrix'] \ @@ -2569,3 +2507,36 @@ def time_results(fd: DataArray, time: DataArray) -> ndarray: out = out + \ np.real(mag)*np.cos(w*time) - np.imag(mag)*np.sin(w*time) return out + + +def set_fb_centers( + fb: FloatingBody, + rho: float = _default_parameters["rho"], +) -> FloatingBody: + """Sets default properties if not provided by the user: + - `center_of_mass` is set to the geometric centroid + - `rotation_center` is set to the center of mass + """ + valid_properties = ['center_of_mass', 'rotation_center'] + + for property in valid_properties: + if not hasattr(fb, property): + setattr(fb, property, None) + if getattr(fb, property) is None: + if property == 'center_of_mass': + def_val = fb.center_of_buoyancy + log_str = ( + "Using the geometric centroid as the center of gravity (COG).") + elif property == 'rotation_center': + def_val = fb.center_of_mass + log_str = ( + "Using the center of gravity (COG) as the rotation center " + + "for hydrostatics.") + setattr(fb, property, def_val) + _log.warning(log_str) + elif getattr(fb, property) is not None: + _log.warning( + f'{property} already defined as {getattr(fb, property)}.') + + return fb + diff --git a/wecopttool/hydrostatics.py b/wecopttool/hydrostatics.py deleted file mode 100644 index 29454ddd..00000000 --- a/wecopttool/hydrostatics.py +++ /dev/null @@ -1,169 +0,0 @@ -"""Functions for calculating hydrostatic and mass properties for -floating bodies. -""" - - -from __future__ import annotations - - -__all__ = [ - "stiffness_matrix", - "inertia_matrix" -] - - -from typing import Iterable, Optional -import logging - -import numpy as np -from capytaine import FloatingBody -from xarray import DataArray - -from wecopttool.core import _default_parameters - - -# logger -_log = logging.getLogger(__name__) - -def stiffness_matrix( - fb: FloatingBody, - rho: float = _default_parameters['rho'], - g: float = _default_parameters['g'], - center_of_mass: Optional[Iterable[float]] = None, - rotation_center: Optional[Iterable[float]] = None -) -> DataArray: - """Compute the hydrostatic stiffness of a Capytaine floating body. - - .. note:: Only works for rigid body DOFs which must be named - according to the Capytaine convention (e.g., - :python:`"Heave"`). - - Uses - :py:meth:`capytaine.bodies.bodies.FloatingBody.compute_hydrostatic_stiffness` - on the immersed part of the mesh. - - Parameters - ---------- - fb - A capytaine floating body (mesh + DOFs, and optionally center of - mass). - rho - Water density in :math:`kg/m^3`. - g - Gravitational acceleration in :math:`m/s^2`. - center_of_mass - Center of gravity/mass :python:`(cx, cy, cz)`. - rotation_center - Center of rotation :python:`(rx, ry, rz)` for hydrostatics - calculations. - - Raises - ------ - ValueError - If :python:`fb.center_of_mass is not None` and - :python:`center_of_mass` is provided with a different value. - """ - fb = _set_property(fb, 'center_of_mass', center_of_mass) - fb = _set_property(fb, 'rotation_center', rotation_center) - fb_im = fb.copy(name=f"{fb.name}_immersed").keep_immersed_part() - return fb_im.compute_hydrostatic_stiffness(rho=rho, g=g) - - -def inertia_matrix( - fb: FloatingBody, - rho: Optional[float] = _default_parameters['rho'], - center_of_mass: Optional[Iterable[float]] = None, - mass: Optional[float] = None, - rotation_center: Optional[Iterable[float]] = None, -) -> DataArray: - """Compute the inertia (mass) matrix assuming a constant density for - the WEC. - - .. note:: This function assumes a constant density WEC. - - Uses - :py:meth:`capytaine.bodies.bodies.FloatingBody.compute_rigid_body_inertia` - on the full mesh. - - Parameters - ---------- - fb - A capytaine floating body (mesh + DOFs, and optionally center of - mass and mass). - rho - Water density in :math:`kg/m^3`. - center_of_mass - Center of gravity/mass. - mass - Rigid body mass. - rotation_center - Center of rotation :python:`(rx, ry, rz)` for hydrostatics - calculations. - - Raises - ------ - ValueError - If :python:`fb.center_of_mass is not None` and - :python:`center_of_mass` is provided with a different value. - ValueError - If :python:`fb.mass is not None` and :python:`mass` is provided - with a different value. - """ - fb = _set_property(fb, 'center_of_mass', center_of_mass) - fb = _set_property(fb, 'rotation_center', rotation_center) - fb = _set_property(fb, 'mass', mass, rho) - return fb.compute_rigid_body_inertia(rho=rho) - - -def _set_property( - fb: FloatingBody, - property: str, - value: Optional[Iterable[float]], - rho: float = _default_parameters["rho"], -) -> FloatingBody: - """Sets default properties if not provided by the user: - - `center_of_mass` is set to the geometric centroid - - `mass` is set to the displaced mass - - `rotation_center` is set to the center of mass - """ - valid_properties = ['mass', 'center_of_mass', 'rotation_center'] - if property not in valid_properties: - raise ValueError( - "`property` is not a recognized property. Valid properties are " + - f"{valid_properties}." - ) - if not hasattr(fb, property): - setattr(fb, property, None) - prop_org = getattr(fb, property) is not None - prop_new = value is not None - - if not prop_org and not prop_new: - if property == 'mass': - vol = fb.copy( - name=f"{fb.name}_immersed" - ).keep_immersed_part().volume - def_val = rho * vol - log_str = ( - "Setting the mass to the displaced mass.") - elif property == 'center_of_mass': - def_val = fb.center_of_buoyancy - log_str = ( - "Using the geometric centroid as the center of gravity (COG).") - elif property == 'rotation_center': - def_val = fb.center_of_mass - log_str = ( - "Using the center of gravity (COG) as the rotation center " + - "for hydrostatics.") - setattr(fb, property, def_val) - _log.info(log_str) - elif prop_org and prop_new: - if not np.allclose(getattr(fb, property), value): - raise ValueError( - f"Both :python:`fb.{property}` and " + - f":python:`{property}` were provided but have " + - "different values." - ) - elif prop_new: - setattr(fb, property, value) - - return fb