From acf5367556f93824fa17b0a352bd63cc7b2387be Mon Sep 17 00:00:00 2001 From: Emily Boudreaux Date: Thu, 3 Apr 2025 11:14:50 -0400 Subject: [PATCH] fix(poly): bug fixing in block form currently derivitive constraint is not working --- 4DSSEConsole.sh | 205 ------------------ .../laneEmdenVariationalBlockForm.pdf | Bin 204112 -> 204176 bytes .../laneEmdenVariationalBlockForm.tex | 2 +- src/poly/solver/private/polySolver.cpp | 118 ++++++---- src/poly/solver/public/polySolver.h | 5 +- src/poly/utils/private/integrators.cpp | 6 +- src/poly/utils/private/operator.cpp | 20 +- src/poly/utils/public/operator.h | 10 +- src/probe/public/probe.h | 2 - tests/testsConfig.yaml | 5 +- 10 files changed, 112 insertions(+), 261 deletions(-) delete mode 100755 4DSSEConsole.sh diff --git a/4DSSEConsole.sh b/4DSSEConsole.sh deleted file mode 100755 index a9e2b9f..0000000 --- a/4DSSEConsole.sh +++ /dev/null @@ -1,205 +0,0 @@ -#!/bin/bash - -# --- Functions --- - -# Display the main menu -display_menu() { - clear # Clear the screen for a clean menu - echo "--- Project Install Console ---" - echo "1. Reconfigure (with tests)" - echo "2. Reconfigure (without tests)" - echo "3. Build" - echo "4. Run Tests" - echo "5. Run Test with gdb (Specify test name)" - echo "6. Wipe (Clean build directory)" - echo "7. Generate Documentation" # Added documentation option - echo "8. Exit" - echo "-------------------------------" - echo -n "Enter your choice [1-8]: " -} - -# Reconfigure the project (meson setup) -reconfigure() { - local test_flag="$1" # Use local to avoid global variable pollution - - if [[ "$test_flag" == "true" ]]; then - meson setup build -Dbuild_tests=true --buildtype=debug - else - meson setup build -Dbuild_tests=false --buildtype=release - fi - - if [[ $? -eq 0 ]]; then # Check exit code for success - echo "Reconfiguration successful." - else - echo "Reconfiguration failed." - return 1 # Return a non-zero exit code to signal failure - fi -} - -# Build the project -build_project() { - meson compile -C build - - if [[ $? -eq 0 ]]; then - echo "Build successful." - else - echo "Build failed." - return 1 - fi -} - -# Run tests -run_tests() { - meson test -C build - - if [[ $? -eq 0 ]]; then - echo "Tests passed." # More specific message. - else - echo "Tests failed." - return 1 - fi -} - -# Run a specific test with GDB -run_test_gdb() { - local test_dir="build/tests" # Store the base test directory - local selected_test="" - local test_options=() - local test_index=1 - local choice - - # --- Find test executables --- - # Use find to get a list of potential test executables, *excluding* .p and .cpp.o files - # We sort it for consistent ordering. The '-type f' ensures we only find files. - # The '-executable' ensures the files are actually executable. - find "$test_dir" -type f -executable \ - ! -name "*.p" ! -name "*.cpp.o" ! -name "*.h.gch" | sort | while read -r file; do - # Remove the 'build/tests/' prefix and the file extension - test_name=$(basename "$file") - test_options+=("$test_name") - echo "$test_index. $test_name" - test_index=$((test_index + 1)) - done - - # Check if any tests were found - if [[ ${#test_options[@]} -eq 0 ]]; then - echo "No test executables found in build/tests." - return 1 - fi - - # --- Prompt for test selection --- - echo "-------------------------------" - echo -n "Enter the number of the test to debug (or 0 to cancel): " - read -r choice - - # --- Input Validation --- - if ! [[ "$choice" =~ ^[0-9]+$ ]]; then # Check if input is a number - echo "Invalid input. Please enter a number." - return 1 - fi - - if [[ "$choice" -eq 0 ]]; then - echo "Test selection cancelled." - return 0 - fi - - if [[ "$choice" -lt 1 || "$choice" -gt ${#test_options[@]} ]]; then - echo "Invalid test number. Please enter a number within the valid range." - return 1 - fi - - # --- Get the selected test name --- - selected_test="${test_options[$((choice - 1))]}" # Adjust for 0-based array indexing - - # --- Run GDB --- - # Construct full path. More reliable than changing directories. - local full_test_path="$test_dir/$selected_test" - - if [[ ! -x "$full_test_path" ]]; then - echo "Error: Test executable '$full_test_path' not found or not executable." - return 1 - fi - - gdb --args "$full_test_path" - -} - -# Wipe the build directory -wipe_build() { - echo -n "Are you sure you want to wipe the build directory? (y/N): " - read -r confirm - - if [[ "$confirm" == "y" || "$confirm" == "Y" ]]; then - rm -rf build - echo "Build directory wiped." - else - echo "Wipe cancelled." - fi -} - -generate_docs() { - echo "Generating documentation..." - doxygen - if [[ $? -ne 0 ]]; then # Check if doxygen succeeded - echo "Error: Doxygen failed." - return 1 - fi - - if [[ -d docs/latex ]]; then #Check if directory exist - cd docs/latex || return #cd into docs/latex and exit the function if unsuccesful - make - if [[ $? -ne 0 ]]; then #Check make - echo "Error: make in docs/latex failed." - cd ../.. || exit #go back, or exit the script if it fails - return 1 - fi - cd ../.. || exit # Go back to the project root - else - echo "Warning: docs/latex directory not found. Skipping LaTeX build." - fi - - echo "Documentation generated." -} - -# --- Main Loop --- - -while true; do - display_menu - read -r choice - - case "$choice" in - 1) - reconfigure "true" - ;; - 2) - reconfigure "false" - ;; - 3) - build_project - ;; - 4) - run_tests - ;; - 5) - run_test_gdb - ;; - 6) - wipe_build - ;; - 7) - generate_docs - ;; - 8) - echo "Exiting..." - break # Exit the while loop - ;; - *) - echo "Invalid choice. Please enter a number between 1 and 8." - ;; - esac - - echo -n "Press Enter to continue..." - read -r # Wait for Enter key press -done - -exit 0 #Good practice \ No newline at end of file diff --git a/docs/derivations/laneEmdenVariationalBlockForm.pdf b/docs/derivations/laneEmdenVariationalBlockForm.pdf index 71df7ad376ad1d9856b9712a020b21529aa54f4c..a26eb454eed104b99539b2f6cf5e36e62e0ffb66 100644 GIT binary patch delta 17707 zcmV(zK<2;Dx(txJ43H!OF*%np>jEf$?OIE7+_(|G`&Y~@oQjCx{W!QPyICpO)Y^w# z)ov5Gi$-^wo7&)nO8+T!pn1Y*G>-wO=deoD*rqDwmu%YQsdf z&Kct>`(4|t1J-qE&4>q(vqHO=V{TO-kF0^;p-q%BQUEc+#fT167 z;$3L>?XC<%`>UcZy0~P2*$+|dzPa)Tl>!Zi-lPy?4)(@Yf@r44vq~O`A%^o+J{fqs zt*%h%ZP28yZt$r*BH+~uh6Nozv>yOqo3(X$6L)l7T$PLJ{vj%1Ch9W9h=%4w3zpDT z5~fv`kNxdH%`W)A}(pA3M03II9{ zz)TndU{3%@a)1Qx7gK>LhfqXvF_X3cmZEe*56lQhy)E$xJy=EME5UBr=ZfG>!XcFs4)G>~O==~-nAz0lL|yN?)ZeQ!)n%1@ooHsSTBHjZz@ z^Dku=dx}2vL`fVVRXQu`eSC!MF&VxEdmyv3*;$m{Y!2eP64IiDWWcYR%JwV~8enBa zI`wGygu38=IBwKnIN9F>DHe0%367g~SA+~0V2sdQZ#`r86(nW#T|6QrX~RG>tVAMN zUqOaEbTGf%#2mUkX4rM;q9oyx5Ou%@xIb{I4hc1)qAFOPNO*M?Hga+%y+| zfN|I?AAdgMsAa}c&qL=K8J>8QaLfq|i-aQ+HO>-$;fOCG9L;1%aO5Pxaey1@<_RY{ z(AyTx{&t;9YGZv9+M5y+oF+aM;ioC%NZ=3UtKAQE+Szf?dn$A?*h(r86^63-quIi# zRR|3=R8V8sqUVFFpz3y1{i__`~7a*US>@EBFPK$Y5V zoK-}94hIOm(;7Ow9Sdk|hXNYMh2y;%T940vDWI{WHuO+ZpQYmxxXRzx67uMj-)lo4 zk-`%65Sm7-KDmv68KG8H_v>m<;>nF%iaZJhAt z;v7(Y6Qxfi9#@JE&>(uF@eA$13gXZr{Cr>FxSOizhWv61(XpSbd=}`#Kmp*y9AYSc zciN%cEQEWJeZS6ac9|k)ykV15Sn0#{zq()lFWdnz?ufN4Tmh(;i4RiE%BW>)0ij0B zqePMe4MYiXKUk&$fCBRa9l0CV_NPIKlCa0g3tVXo4yP0%2QSyg>_Het}oLt9x&MfM)@B~331QPFZd+^V2K@2lO z%00ax#>?EhaD7nj!t;2IJKomDxj-c#>I!9UcC+OK6wr~^# zc{G%9*+td$!x`2fEc>n;(m9xawTs`_AWfpwg){9NtE&yFyB>7L0h7D*7!ETMR7GR-@usWk+^8jb^VY7@R%VH|B1bJApz z@BDS?EYS?=pf3wBx*THWNax$67eH>z!B+UeS$I137@LlBiWr(y;ob*-Jj>@valfB} z0XE0*LM-t;fN8>=^Uh)ndn8I|V`f)!+oUcIbeRl#o!S|&A=Bu#ibr;O8)ac!cNRko z%A!C`_-WSMKLR8DC9ZSP#W$lR&4UXp!H^7q-{Lv^DWMRGG)VF{~2u1QEXFg7qbmvQw0D3{BI z0SSNI(lSVQNOzZ%bi)7x3^5GD&>hl^bT`u7NQVd_ozfv9(y7A5z4zIBpY#7~-MeP3 znR(uh=lvLHHMKdUt(`4F3eFHX2M;H=2tZa%i-#A$&HaLto0}Jtfk6ihcLe<<#$?a~ z!Cb-45Rrcw$ihHC_@hh?2!FIwbA|v^+#G)aJTCw|0wO$uBHY{nUT$vT{}?*ML;!L? zcd#`;jT4~a3<0@fGRQhZy!XJo#2N$xJSGQdE2{!D zpdiSf#;SiBumk>@4Sa z03gKrw;|Bc)%npM=ne!s0xcg6{-}Q(2vCrI1pq!a_*Z+bRxmIW?#k&3cKqEU*Y7Zo zQ}H8xdDFv9Yqre{;LY^e_Vf+AvVqc;XjLgoYH>^cK>Vt%zxbm3*f(FsXITC76f4a zhvO#P{M=TLUp)WMO#j2>|Bv9mqx|1E{=Wq&xH&rhp=SO=|Nl_~oxqM>e;GV-)(!s1 ze>LaFD}em3sXpk>n^gl@gWa6|SE~#MKHh;e#MbfeF@jwcz@8v$O)!7l%I;6H{3+M{ zePoVc2uRb}75w{E0dVkebN`R-@v^KO9$yaEM@Igk0zDqif2Wj(SUFq&zA#>X0RRvN z1A1X{KN5+TpC91E^LUEZAkRNI4B+B~IKv-Z0FUkY0c@OMn7=nFzz^V({w?|w2>`fc z{z1Y3F8RNapdf%t?Qef1{D^>mqu-TtS^o?20JuPZL;e>4E*r2r=wJHBtl)n^Apn=- zzu;qDr+?t%wQxEA3-SWEpnpTY$BLnkG;;n|b{;+e7wn(xkEOvq{vF`4bhm#X@1x%H zU+^)T*FWJOEAs}y{`CH@!CSe(V2@AFAE)uSmH*&B-Z>D+6J&pdxis%=B^qkq5ZZN9 zCr#``IJcX1QK)hr~s z@fcL}Y#gds84rJ(-+i@LEASHQXT>x}tZs}+P}wX)fo0Ad*>f~F2lW$H!V^!t*~_aL z!o0}`8kIygOuuijFG_rjHhAN1=e$2x>F~PtlGBqLkx?UG5{^$X`pA4qRUuvQDKA-v zji=2$?Vi|ldc))Rh1q*oUF!yRaAwBHFe$=JPR4~ytQ>zGv)-YU68Q>46M$Btu#EmH0 zwH^5xOS)zb$fThXA0u@Db0kN z28i!P%3gM`qWXSUzzhz3@EM2%EQbeH%*(T{ljr_$=|Uov=(guJOY=NTZFn_w&tx!g za1?OxI0E&s4Ta&YSCF#_go2o&sSJ>caVK)!SNd+fAaPpiqD_K05?zZC@c?Endpf>Y%c#)*nlRxp%W5FoDY)`4x zDkn`k)zPcYaZjJsOX_8oPm_91XBmLLTTw#6rDHF#`!etj(gaugC10g8KF#7o8Z>_z zoMlP7FJ_~Nbb2+&L?kS*xu`6NRtOZ@Ws=KAntD!Sz+H4Y)|Pq8_I1+wv#HEMD)BUq zQyo zJ%l6Ojlge`w2q~*YD-m9ytYIdUsgDYU+WACCX?N@cuXb-ut>Hem&iU@&rj8}Bo;K= za>&Ykd#NH?Wrc&O-us629A_z3k*c?}cw@;oD`|G%S%hx!gVjXMdf^w`_vMbW0h`GJ z({`BxsKyw)OY6;_@)DlW7;S&J>t1ATId*CKGO*~4wBb6PbSzZnP!kf0Gj;(n_p6>> zsKgDpriYqS7mYV0WuSV+1uRkc`e+3g#tIBJ<&LGvq{R!`D_>Kf2a$#cj;8abQT%EI z5=Z_lW59k$mlEW|l6pWO`O3$fMAgM9j`$rgL1U3SZN~(yn%-vfOmTmoFz^d&kc=KD zP@I23AV)2d7HrlM-&`=0m%z%{&UXr3gl&1o{@mmTAwSEHH?zaba*#UJHXo1MrIB?s zGeaBqTW;ZS5tizhS_eI>VMx_?-j|~1xB5nTs71k&t%!x!HA$Fva>X2kI?-!IQz`^U zl(>?MTKUqnmX*jl`QLw4&jg~)@9H*_EU+t3`A)E|hQx#-6IDe4X?CHJxMU&f+Ss3r zX&I|gi4v%Ghv)-*=dPU?kR)`vCmrfGpB%v4Ue_5)sX?Mz-9GADFi2A1y%M6n@KAhV zOC5MG+J$9Zg4VX`IOdEoLDA0p3p-hN^Ze4SS70k$c0QVcV#t3`=aNb5B4&ibei|Xq zd7g(_tpehwdLPe*aFQ&dao@}}ZgDLYvOkkI6Sqr(K=PiiEv~6j1rs+^20MCul)lms zI9I1cf+3w4V9Whf`PCcpbFnp1w6q`YzG6oL=(CDTJ5umi7lW`EiOu)%P-W6;hC1eF z1{dxrV|#*rJWYSXh>n%gwJ(=z0ugBhN(i2$rkW=A;j?AlZPYs^UYHAzoHM18AAWkV zYT>}DgWbYh^VT$X08fBlYmvjak6FI@o3sk;3s(_T5UX_I*VUqNY*eHBicxvjCiSTp zrzWdU6{nn>E4LEZ;r$UaXW0_3dcqmoDN{#0YjzLE+r57w>9{=+VGgq8bv7>g97`#` zLIhkQozYp-p=_L2w)T!j`LJhpaLV|Pk}pJfOW;eVOJ0nW3K_H>VplFS$GAszoeG?4 zGq^rvZ#P9VrMV$@mPs{|?t?^z>L1xu(uTKdiG~Ci^Xb23{YVhQ@*@072gSRdtV4!S zgk0Fyaq55Uq+lUjhLe3^tRx5*+MUhwt{=)#bdC`29iM$}Mp}IvsD)ol(KrPSw0ovC;t!cZpun^y5 z1w>>vv!3Ro1@-I*@8y$f*9^&76FDpnP$XH>*WKxzj)Ya1!BhGC?D#>6Yf9=7utBA2 z z4Ay@!y31)zYrPj%r4~3_0^i<}VDleEB203`C(Efe@tebzly(oQU2KEPYPp`l_RE6t zzt~RHYBD58djhJgzxJf|dgjkc|FSw9N+0W(EGYbGX#JTvlvO16#yQs7BBHQRTJ28L zGBhTP15%tAeozZXNiwBtE{Pm=v}Jrh11f({tH%41YRPX+G{veRkHMdSQ?ILIcSru& z(Y;5wlhT|8@phJcoE|Ehl{tJ!jn6jD+@=H&{T?h)JOlohxGO>Ix=UacG ze&_enu%MQ$*?L~W%xX<9T40!vTi$3G&OZ8!n24OE3VekO3IJ%(G*-{#@k0noDZxTA z{P(PT!=Y)O8<1`(nm+j-NOGDF5*_&J4$v1&%6B+ zLg0gPz^(O7x2EUiX$%!R<6O}40VPeH3B_=A!s^baLLOC_y2?pHI+dLU#O9OU^ZS9* zHfks-DL2UkfEeMU!D|&lg2b1CMyhhpfki17HhaoEk}X#FWVnCxH^kVoD`_Jcf`Sbk79&gaD7DpU?(67ZslPt;ZG6q> zIvDd2yGFo(cF>j%Vm%pdoa{lc6B-0z4z*s2OLLxernxRNfBo%^GcWG8*yqj)45Wy9 zsu7$%Jfvsgo@xooJ!Y8Q7uCk})^a*ujpp%hv8M;@2U|CSlfk|kocn+B!XB2|NFA<2 zGU(n3oAa5YO|;FKytXGu#q&Rj1640n*RcAqIuxr##Jw zoqVq$i9s)^`!r-Zo~lc~-4qulwdEH-(l_&*7-vTH3~CyDSk9FXcvjc93>PBB_fyZi zm{*RI8w-7tlpfhrnOT2A(?c7SK=vMK4K?t&sUKpC0&&k!Z{m`?4iWFDlNMMEqqHCP z0{E1B0#mfwQ3C9DG(Y1wRTJI!SCz%DsQMwm=TFpU3{G6ix&F z_9t($lb$HJQ|S!@Xd+*pEOO$bn=ly`UJb`x7zIb=5yetRM$G%}5ehS+XLk?buGzo@ z@itomUv>FIS~&o!{E@E+caOM?qma6($c&SN12weVKxL#qC3Y1DC$UPr5$+tK6n27Ok7zyh@6jZM$xR8WeqMUjY)6MlW6RG0)!;+g~+9%eM0y(i@!rk=J=8tlXP_1ik5^dRmZxTnctmhn6;{mTN|@S z`sffNU`~Hgg%$63_MLVtw=eM*Sr&G(|+dwrOEdai6dmbG*?E?pdV?3FZ$?Zu;FbCEmVSI^#{- zddh0swKlBlMWF)yB8W{TBdyuc9P4+POqVobl*IS(zw}CJoY4=y&KQ6T34SaZ@H<WBOTJy zuONw`H4vBuWg2vza|QUh-gr&qWFCJpulNb`L>_f>_wKkZrAG{+UAv)1?p=`X*SU@q zDdKoVm=u%Bo~a<@vjj1&_d2OXGEkyy=a2WGXPTf$_w;=lKNEz?R^=C|mJyby$h31` zUJ+boY5Ct4IXwu7hQw84xNB~)5GmIF+}_nBW58>$F=dF?fCmM{x9@I~Iof}8Njso# z6%b0ipU-=}*-YctA{Aji+1jxEfco}apO|sKoCTHc_TxfeU1rSu54HkQ`u7REdtmd7 zw=y~2?kL9=RHm{ABc@htUjl=sVmd;bP+3F9Zkua%)hM~!yTYyu7EC5a@6nuxX;#kT z&Ye2+FyFG5=oojMqZKX914Dn0nTB9-%kP$2KZ>*Ea~{CWo2S{C0!#ZfS(eW~mF2=% zX9za=^j||4=<~6tbbY8fVp__lf@LBp&jEQcVH7UQ`5A?nKnsm=UzhDj>vYtDwat$& z4k(ak20G|BzjB@XPOE9%wQ^X*P;S&>zZDhOQ@isW${VotKMTAz1h{`y7}q|bIB=Wb}E zi{=$_79E-i6d8`soAJBgM|i`R*H zz}UB%bZnY$IN73#riCQi{q*x?Co2j^$BY&bwCfjM_59>pIYBn3%47kLa4fF#y$$L96VO zZEVPQ5y1>b78ZZ@1X-b(YR&j~%hn@@cSw3^bZF0W?y+s?&{_Pnc)kb-!Wec)%y!_F z&XU+=qzW->`22-j%X>=Ph9L4FmF33PUixz~lV`Z4jHb0%?RJUTbR7E(E#5$Z``oV6 zB9;nmVm+%^O+A-(1JY%-?;?g$D381k1%#>bp^ z1s`XO)YXIKD5a4+~-A>-3tFDkoV( z&58RYGxy!tkfyEc%_T+^HHXDdf|!)ZH;JbciL@nU9KLTawe-|r0vOix1opAll^!)9 zioA~M*du>t22Kmkt1}#UZ!)F@+39FNMDCZXrg#X?t{ngm;ytwj(#qKgq=Vc#?GUeN z)FIt3g^_vmduX_^s*l{@Iv{vjMdf5X=oQWO30Eogy9m|szAv66cyU_z@$76NU@q&T zFoZKmI}?plK6q_QRP?zJ|!J_-ID=WW`8~K$DB#Nu76g@8n*?gSD~^g=P*n%N`oFw(In$ z6+?ecDb2|u$+^c}w!c$YmR^!#xl_hX?UN}*lNOMR=uUaq>;=LQ<+g-0Soa8VyBCt& z!`|Pd`>r#>)-A7)TitZ5psmLrn(PoMi_-2(i?`cnVlwTv1Oy7q+t@eIxGz1Dmd|Fz zGY$@|l0u#G>nIPCfm`K*1JbBJvF49ZLM?xsw53Ni$1mb?M5@x#D}tY3VTt9}U4KJb z_$Bq~$Pe743azr46hHM7yNwnVXjT$4be3&pcED9W#3Y}E!=-+&nY!Cd(Q&p-_ z?9_uiUg@UdesBJLpwNeGSFzQ*?lVBN*xzQmEt%=W;hh1BlU`(EZSR1$tm?O7(T9IB z*Q~cppLzyXo6ehLp3sCngA_`bE6${|QfXi%YeCR?l`5{MSK3}IWQoRhABEn8-PxmZ zk2&Aia1A%BB)oPj%a6We%`#DCMN6azXoRex^;QZUo|{9y`qlfUaHVHQxlb0Yjc|s5 z@FtMQWpnbPr^TXid?y-{ArU6TS08_9?wYnGg%d_MS&V_^f~Ug!?=0%y4=pDyyCp8u z<;5?E)RBf;SG|8nB{srQQNq*jRZ5EOn>!Qm9Z$Sfr!jr95#VaLA4aE_$sTtT>oH6z z7k+}`tU=Fc8c`%})ZklM@Wdok9yDqhKS$x|oeD>xzvo!PZev0Ak9*kPyOe*W{WXM# zuw~^DIQBK@T4U(EC)4e|mpuZ3hl9(80U_d4j#5iJ~jNr z6KDF7JEl-P%*8t)J3lDnhoG|Nk zmP;kdbjlkfT#ltWDPoxu$_ta|!C4nu2OIi_U7(WZmtS>EuUaLaseV=_Y~<8fT(*go z6-{RkOdf}SXG=)Pz9=bbH~p$z+ElBAnR-ihb~@Ij?)a8Nho(M8G1Y%QxFW6N>PIaC zy|`bBwOhldEt!_2pl^(H3#(vH9VRufg|hxZ z7o&dJV#vEUoAj*BDEEI2ou4`(6I49P?yZ${)YC1_8>uHP9kMUDtKN;&$G&tb+n|K1 zXx{S6uN`MZ4x3aA_ryIxdXe9+7wQs58Fh&<8zPX-hGbf;Jn3LKv>`;9@zDA_tViwj zZMsW^wmXQ3bXO=8e~I$8 zjG~8d5s;q8v%}9uU2|2#EXlTpRcRrcj^sRf%Z%~pb>Dx7%*0vSti~MH2ydEklb^Mn zzrvX(6$?Vj}!1^}9Dw3i1v%i>OiSKI-u*sZZ|AmxewSzAz{`?Re5g z;~jrgoO0$%*$v$EvI>)y#q!FFxi@?s=`eUBK3?iy_qGG04@#uZAQIqco{`&OLsItx zfPRgS0$r@%MCiT~pj-Wc_H6b{{@RWO2PxQi{_F6OIfVPm?OSI-iGx?w9pD8bOZoSDAlwfv3q~Ghm%>{*&GMH?Y zlshXMwN1Q1&v+uRFu-dd^>$J?qw9E!WpCBrWBo_hC07rnssX#EVr%nC%xTWqhd94m z{^5VEper?sVJ{h_fnr2L@dX^3ugI}CblA)s3Y_pCmK2LQ`xUl3f^4E>oUQ&fL$`VZ z!K$>9tzzE5M%jK#T;+kI=PauGSR0=anNBcK7@eh zv)f;Z$*a8lJR*K{6|ejF7z`UCuN?K<#58{p-I>7>4&zKPV(2iaN`P{PR>h9unXP0? zBDWqmud!;Ea~bJ}JP9`3^iIoE$(Q#={U5R)P%b_+#OIcbp7Xlx8f z4?WSET9Jv+`jFq+_rzdtzxD_2WaBr!WU{$unja^n0y0t39(;?$;;CEXZRbov-%Wp3 z1M#XB^Vn6E`n7OAUU!sBH2wG*+K*#oIUpwN$K*lHDhkhE0kqb&$JA_#B->tcpe;7Z zkF}@WV_~8)SKQ6*&1obdQl7u1q$1-hui$7{m#SlLskO9|Jo(_(TV%5}Zpxruz>#DT z`vXHka?|Y0bGGWUq`P&OzK$R(-^YJG5oBL-2sJzL#~Lv~H)9}noB`V>b4KX_UYCUX z7^~vM7BE@8uHLGdh#(YTX@B3;u($i-R8L(%fqJ?k@|FwLW)u3q}_j{0;?gi zA`ep&=?ZojqF$@)_fC-4`+TKdi7gv(e$@eNdvj`+e8o%a!)Vz_dF`%nr{uOEW*@*- zG^2DjQ*E&h5n9{UQMH;Ua#a+YgR#15l-JHF7ktRuhshX|;=`!jn+_%#Af@`pIJ3yT z^R`tOjc4C^j5t3Tz2fY*VIY6xrOmJ}iP&m?R#U|KC58#T9S0jR-)Y8%r5R{{g2&_N zoIoC~u8*6gV0a*LtA; zZmaIl??SSBovuRbd7)(uE?i}IzREG8>=XK7w4F;I$D2D96ExQe=A?hu##u@#gvfwx zpSO&sQ~CooC0ZC)l@}jrO5rO?qjy$n>+up;S?zRRB%L4hv!zvqa?+|Ws$~|UU^6TV z^>W2?JdD+ipM-a2f6xh_x~L6mO4NEM5_G_`Dv!$J_uzi5{gYn?HE(i#R||je5U&CbJV7OZ?BW1k)qfq4vu6v z4K;E1La?$#DaWG$f8d1Hf%av1HDPy%%oJ6=&%Qp#j}srDhWfGMFj4H`bYI#XlC?|Q z@q&5z4xPkX+%gE&JkIxpue(<*sr)Ul3rU+JD}!zYzzh=0^(lX0B!!4a#ish%qQAs7 z3rrKgdvx598+@Z3;jeF`_aN8##avaq&%9o+Acws;b1z>M(=)@Tj_ zSmqc9nt@CuG!=h2tw!eVEUOfCqkY`gRr!_?Nr559#Fe+@Qc0DO8pr4A5IIlHv=2@~ zqJyBWa+|u7$53HW;=tq3A^y$s2FaqSV37;_tz))4zp8731edzly4x7B-=Q4W4Z+>}SZ?!Sf^ zb5*pVg=15yq-Zfky7o{cffGURi&3DPekhY6G9b1K+TL~rEy!pq)OykaT0I-W0 z-a>T0lsA9W|FBk##GL}Chmv+m1kx=nX4Y*SN^7;TR@n3-p)&YG_?h!px!TIPjr6X zTolb4OONKJBAPX76 z2hJST#W5Xzxkp!ivvsj9IxEITr}TfkgF{WmF^K`8t-GFE7jX!WM+#0l6&l*!_z?sU zt^_vg)J=s>5NUiD&!3IPd$0Jza1S9omt5bw$;8&|Al@9MAh2K0qPkT1jX0}zjfaE^(3DAqjs0)LdEoG_6yphdu|yoz_CAV#0x})iSvytacsvOCJb3EU@Pv!q z=v#K!Eprh9TC%UOg@SlZh5ipMqjj5?QHB8%hYYX*w+yfX6Im%UF)%PQFf%PMG$$}H zCn*XqO>bmGVRU66C`39kFfuhTG%zzUFfcPPmq-Z%!~rswG3x>-f89IVa?{AR@BRus zo>RMs`hKaJ+EoO|A$unSgiBy39-KrZKHIS~w!^@D{aKfm)RL?alRTY*v|6pztJnP^ z+o}+SEFuiPP{tL~6L0EVQ1BMawY z;(;!}Dn$f>yeyTH!vw>aVRN(}}(*GLU0Kwg z)2I?pIwzQfnJ}bzkhm)flmwlU&ZQBmDA`$XCN@mu3=Nq}OcP?`QUCc|D}(H2FuW_R ze><@kYVibQg<*1dhB4CYVd)uZ)O!%@KRXl?TulTU!>|J^f&eVqr&lBrVbDII?2oaFz;MidUEz=E;X#D|eWyGG%E!3hs zOZQOWGEz&5PIKi^E>@-;yqfmWsuQkRe~A{2Iqhtq83G!BJ%{q~h|!uz27y5!7?fzB z6%U!qIfo2sHaB<dBO5zw5RCsA0Ah{hm7KmheT`jW&E3S&8T zF+@I(o>cSjF~=!I&<*O&jW-xl(>fHt5i4%60*bT(Sw-KF!l8<#pqPn5@;;03f01kq zTLo-HEZGqRJ@7|wX*UH=zz;TK9Ytu2Nzss1DGZFzRAB6cGbmnzJMeQfQcq6Qj$PUD zn)hhRP9taNtYCf+(;#|7zk?lx7o;pwyVd9n77P~lS=5aQXb4h46%>M?Pz&}7Hsd`K zr)cE41;g*p?lXck1#Q`F+dYudf3SVEKCfNLfZKb=en;b0P=)OyL6V0gKVgrj5!10H zE8vkKh*gEUPh&kT0vGL zC@Y{byNF~2NraLuIzi+abX}(?s7nk=G4N4jIgUJ#f^p93mN65VF=%_kf3cRnCK2(d z8q=|j;Bm(5f$Y~qb5Aocz})V0?`8p|k^Ff$BA)yYoYIgx(uD&WkNK_?ItODQ`3*NI+a06&xDz?kkH6MVwmoRR5+r3j&fl*2igBNQ)5&NFx33EB|J+!zI89UG+~ z&?JTKad%wemR;m~jBdEXf3(nuDPG#-4z5!TM~`(ak6LL;nV9;m!%BB-42G;gcubuZGsZPzI-$i843+;ariY9R%YGmK#awW8fd zql>DcJ?CS{6R#j*)>k{rHI1ve{+(zof^|dt8d!dwv}Cbb=6U-pSsczWIvFZM6~gvQ;1@U{9RF~hc%74-UYNJjK5g*J8F?oN0hxor1Cj1pMe;X=4j!`e@W0FHHnF9dn&**mXnrL zYBbi^j><(#6#%zs+&f!OeS)EOFQO2Q!E}hiCc$IZ=}eFrM}u$akr5BZ zBz-ys*8r0G0Lz`8!ur*5hoR6NjdS#rslsv*_iIN+Pj1z#nP^a1^^_$Csij( zE~d1OPV6Ohe;vAf*PwHR(}v;)8t2 za1J`;Kmk3}0gD<@D>M~j9K0m1lX<2#!ILft($gy>e@kv*xI1^;g;V$q81S$b1zBT~ z<)VzMR_d6!1R?q7b?soD%Brj`&PxW?ujaCJixnXLwvf!R<9(!#Iwj=1O=(Szk7;K) zRIVPI)%o-+BkxzuHaFg^W~s%4qT{z&%oYU)Hb`fYG0%zuxTL;}PjFG_sT&pW8tqNk zxk64aJU$42Qvs{ycw3j>Eu;NeCIUj1w2{p-?a^vSr>ez~RAIZ&456t6*JN2qA=%Sp zNm;p1CxcyTsi@p}=$Bm+15*K^m$4HA8Gp+^twtPqB~sR)eGM!>Hx}tI1`8L{LeVpp z5YN*{u|#=Q$F~QJh^7AX4l^gFa+5}0@!HVK(pop3GATwb$PStXuL%lsCVqCUX{REu z03D~6v&P3e6K;yaQ2f=n5eh&25yCY$mV$K4aCg?;AsVACrQv#}BBh4rIqje|D}Nn1 z6;@VhZv(~8G#Q)@lsahMK4WruP9whJf${|32vFJ-ieLRvEly@aZ;T9AQ0et0*dzuI zdhWTprQK;Bs9utSdR3p^S=aa?CE$Do3nK_ZhSX~@F%IYX3S&0ZdE|5Dv3$v4(q3Mt zc*-mI?oQTwbDRW&odUjijj3yp`G3HK20zL`&dl>Pe~{D8ea zp}9HRwDanFJKx;!yg9SItGCciy@zgK)B5)Mdb*h{7wOoJg48YwI)8rp`uO-0DCqcb zOhM3wvS$hteqy-#F2j9-4#`ULYE(&2!ERqRHv5jO_6$T141?`{MZ29|pIuDLIkN*M zoJ{xTYF6H?XJEY5W(tP;t^8J<1G^pJn4P@(=wELR|4qWNTV4sA!P= zC9z2?cKC9;@b&TK%YWBT4~Q?nUFWh6#$G-WUqO!Y^*gA!tK=&PQ2)D4)l2k)%?$nR z>(^EknYMe;YNsgv`~Is}rzaqnlQCZ~qV@&xAcZ&K%Edc8VIMHefjf-j{{voeSRRpQ zyeeN^SC>=pkGCMY_vI;xFN8WE~`Ipr*jbe z*RQjZ{(q~=uglvNq>8KZ>c`DhwWw7BS>t+oz5*ohZw9%e$sl!9%VJgIoAcXs`Dgj( z?Q&D`157mQ>UzfXtgCNu!g}_-mwdiK`uSe|DE}`1u2#!j0TKo6`m=7QJoM@9i-T8( zq}h@5tPQ8<7{cQy*4}@q*cYHkRshL*Msi)|d6)%2%=cN}DZ)QHJ>5THqaJZ~t=Xuf zvH0=ty;-*-tXrM4F1m#m+VT-t?P*hfU3IFF6>7c1?M@;4^zh{Qvm@lzjg%X;2g|4* zMbSZzDjNKv!^uuT`sCH=;rrLf$%xuioqVKHkHU!`k40EbZi|0FuTu>AK>F+jnkK}g z;Ig&syd!7~Yte7u4+C%5Da7v|9PfW1`v@b>Q1>Cii^q|w83`R3{^F` zDkTN5l-alGyjp*pSLO9=0XgIPcD|Y2%zu=N>Dhdmq=Ph6s|wP}Y6&ar+p~4GDQ_1L zNmuLXyd$ukV&&=kHz#jcU}J@m1$ki?hq=NB^%(LRce{JYpu_A=q4MP)N1ygFC%$_( zTJC@^X#BdF2kUp4D@fVl6>l3}l~TJGJI_;*gY9Y^Oxk}|l@kL-1Bd_}j&};F$M28+ zb^P%kZx2q6@xzA^4F_-ujo)ZC7QfTwTeipA(Xn3CiqS-Pb~~R}UDw_c$6hSXmlw0e zB^(9yE~tBfU-JGewMz)T3|WW6Z_hS#i;E!T{pq@*?|2X1!*nEtKR}@`j~~od>kX0$ zMm#Bxrs;o%UZ|YTE;d(d{0$6&zFWRuKnZtI6=5VY-9=9Stx1<0EBoXamy?5LmdG)7 zK#tskVZP8d)2%eemoznXS!m!b3m z6PGg^16ff+G&D9sMM5(;LN_%xGD0{vL^nn@K}0e`L_tC~I5#;yAUs1fG&Vv-LNho* zH#IjhLO3@>H%2x=L^4A}K|(h;H#t6C3NK7$ZfA68ATl{KmofVSDSsW!%`1da6vy%N zyN|&$m}k6<2V(}q42ChrJAZ<`GG)VNHj0&(OBuZMy7-=PKBt`-fC+(zzbdoNTAl)Q{l*+cF2eNT% z)eAYQTJ=FL_N|hT*LSOaz>8U>fa11Q8aUXp8UW5$tpdoS qm3l$oc~BxX`*a1KAJ=1YSN;KKBr!6VCmjPH3N|w`3MC~)PeuxquSRnK delta 17655 zcmV($K;yrVyA05}43H!OGBA@N2`PVCOLOG55x(!Qm|Hj%5yAU$*ectxlsIMW!B*L+ z?Ale_8SW0tGo(Utq?Hf(?db-<8FB!Ja{Mx)W_{`x_Eb@l2G7s4! zxcabQN)fH)LR&!;Ws9rz;%)YI*;Y4I{c)KK%CoCY83q?iNtOYzFU`d z82-34HY`q0;j~-#}3QeeXpJ zZ>xGBZ_Bt~UG2)c1B#9W4RTJn@lrY#mGF*I7O#kn&VnyEPaN`X&J};nexF*B1OL%7 zK^K}h12+#Z1tr4jg(QrLuyLJWZV7X6M@OR9%bbHzQTSn0Aa}4+_~S)}gJ?=D*aF%F zrSe4YAkFn2wa!s1RToB}CvLO82?V!jJ;7$-DJmbR-Xqan81Hbr-}cq@_AWREkiN%- zcVXN$yE2e=t72QUam#m%EJz48~8QW}!pq%;O5Yb9F>lu!O9V zNbM%F;+4Lmj1RX83&uO<+!>-RfCWRpcVW@r5)vq)1Y!#JazuaLj3-WEK(Pv*sWE{r zrGG`3HQr$5gW8t>ohBHqe2Pn zV8RR)U{4f~!~hBSFUAa04xxzTV!~}kS&GsTJTN01^|r)E@L(BXuqEIGY#rHsd3|tW zfXF`qSk}Mm0-%3_wq(0x89YKbh_H%qi-1w}0TS+C4->-U7UK}0qE&yBROmee{UULV zZ*$G&E)KR(+`W|PYQprX=_(o-%8BW=kHm;6xsMS?H#iL-hG##g_$gSG3JYGYq)X*+ zfNcF`S_U7DC{5fLgbi(;=1}v9A|G<7w7yc8g|KGrjGuqW$!&R$6RG{jrJk zO?dvbY~xJPg_+2SgG!aoitRo=LiU&p--11m+1czYa`3&E3Y`qLk_yxcLtgyR zY~j=@goYX_h|zD*^TAb6bUTQ?D(WyOZ{jS=0$A?7Z-*0MW*5o2Zb;kEfpIF9IQx2C zw2^-b{C-_wVb;aHn+Yk@>GOI>9)(6+i8p~Lnp75dpQ|E}hb(ggKs3Y?WU|BC-D8=F zPsP?NQBZ`%I?VYSUms?}J(SdE>9_>0@{gs2JUZp~ z+7L*jumnAXropODx)HD<)T-)!T@6CqR=ytwqoKkMPd*UAZ|n=0$+5zl>k!IJaOrR2 zg3F6@(CTI6KH_-XDLO!n+8c~tXa`mhNsI9FU4e8rRnhkO)dr$tH`@6uz=weXz=?mU zilN+T`f@!J?nU<9GPl`P3Y^Cv5;}ls73gijqbVct3xQlHWIn?rMF~j`f&XY|p|I1c?wxyvyyuKSP5U zR)my$+#tqg?pQ<2uDVc zM?)ExT~uw?pJ5HcvTw^ior8Z_JNu0d(j-b<=(@y<+az=VRA_-mqEPQdiDS|V+jvRK zK0ZNp)%wcvQ?(mDI|BJ!pUKV^Nf6E~q1ZJU=bhu!8Uo-Jt^;#w6TVtu9Ia<_(qxkF ze7kg#Xa;f6mjxJI4l#3p^G(tVAUEb0pFtKy z{XZxer9M-KRf3-&f|)w84&NkLffqWK$^wu*wHobH89hDB_@+38v9(~LQj23PWk-px zfK4>284A6@wEVkQ=U1=(2}jQZAcOY^xAzDEk9`9)H!_!TumKkWG&PrT>jEc}z%c}W z6zUo-Od}!P2n^CaLwC1yO2;t35CaUI(hbtxEuE4INF#`ZbV!4Ubc%3s?|t^(=luU# z_pVuMW}dg>c|QhP4Q&o7D=N{sNf6EC7&@kO2E1?f@xA2pkHwfB{r3 z+-xC^j~T%h4ghT@FcjkE{a+!NMQz>OoQ1i#JUu-*EgW4ro!~a&EbIVJsGBW+Knvmu zfqOu#0KXdss988d{_KnslL4S(3w8ZduI*&)=4kfvx2|@kI4br zN-6+#X9(<1W0gM**a3gd1_0s&{T=Qv@85x-us@tFz+fjwXA77&6lMdkhB`n1>hda_ zZeDKe01KGaZ$k?QSEomR3l9r_sDp*&qro4QTL9#xGyxWm4gS@hD;N%Sc5~%)g*yCh zk?VJu$0^IftYn-V9U(9`SIpn_$wJ`}@Z;FMx&GX&JY*}?#!hN zgSxmwlw|)hc@$y(BeQ|H0r-GGpnwoB0OA6Gc!6!Xeuvlbc82_+1pO9&KX%~n>+IwV zuzqX;;t#ckJpN$%x>|TZ0B&%1h`;Z@75|MeK_Gw?6zm4DgxEk~nEyn76ho~4#E;8| zL%je-z(?+Z0Kng$zyFy$(#y&T=HUI0`5*h`da0+csjtoUXUBi*q@|s_0KObNyZ{bv zJ`g}iP!J%%2L$;4cN7hO3+P`}0RM4Sf>}EOg#Ik{aZ3Lw*yFGLGyioPEP(%xrRMZV zS_pvoAC4OX`GDZZFVO!p)BkY!|0DSCDE~K(|8GI^?hX!rsG0xJ|9{jLj!*~hzYHEZ z>+bf*e^sZ)D}ep4sXpY-n^lEaLERnySF7Y^@puPPFdK)z#|U+Qm4|vktTdo*VB0^* z@~2$)_mMe3VGs=`SLp9o1;7CU0{=(%cv)ck$CtzPk&%C>AdiRh-znu_U?;2J7sk!U z53qp4Exa*-k3{0;;{*7D9#7E<;`IlI0bHCgC$~ozz+-#<0Ba{W=I@R2^8>h~e~bP^ zLI5tgzmb3dfJ^m%ZzS}HEdD`{mf*k9@5;HX{sln*F38`Ij|aeI4fTNhOV117vi%q2 z18_nA1qA_I4*!CWsT}`-kH^F1^f%;vtl#;OMo$0A4tk6Z|0nxnX>OiQ|586T;{Ffh ze$;#Y3qFST{wMrnDLxSRpWgp9c(6Mh{`mC#aTxW?#vRkCbv5==RcDSPLv!-H4RW$xN_`E#a>w6%@mZ`6_Ssg`x`ywcnlds>2s4v^lW`#vDaFLDbton|meVGo@#Y>c zE3bJeoeH#Dy0WV$K_Ao|1e*A?u!OdIZ}Xdc_RyJ$8(*?UXxS2{HqpTP_#&tf^P3oP zB8#_ehVJ4>SB;0jE2Aq!Ep0+L=5kc+?*cb$avD8@_E&UhN0z;Ym1a(qeWoHk_BUA_ zh>p*H!P!oO-82UslGb$guyo|xB>1K&1@Ds?L#r0Qx41at69Fci^W$KSLfMXZmzlPG*% z7nWTmh*(#SjK60NFFMn&`EWM={_M!4m;)$(%2{cnEscw<=8!^Yn5Ukgax`5S&qrx0 z)Z9mOFI@hjgB8{9hdgF*=!0)x6kstTuyR(8eT_Wtr%M+SsrWlPplQ0-L0Y5cz&+E; zzWu{j`(aUCP;%@M5pr+C5seGCHk8dGBogt5hWTaZU20H}wq8A=|LY;)vYypSi!&>K zpl$!0@2gj$j%Z)FJrBq7({96w!)*8l1CWE$hK~?1Z$h3?3MmBJ2F3aIdk?n9U*RiSoajk)M%pbXm-^<)ASxCWh>s>ym3zm3hq%i&b$Z`Ou&c^Mx4 zC`lsmna=*nnX4U>ta?tv622rl(i;1J!{e96v~3`rC`gG;-E7-+rKGceSpLe!y}obr zki_z<-4dG{TKM%ta#G>(&!_p08dWS_Ot+h_8M*m=&7G}hZE>SWuO{yB(Bq&ezidvZ zRI4V8JJryu&u~tt8zl6yDyB%irm|l`et;>Rp=F~luzNFc_tOQIKS*e1Fp^|{^CI;d z_Rp}S-kA>0EI!VCXewJs?0aQYuR3hu-r}9*i6aAlt=$Qm zp>3gny3hM=`wqId#6Gl@CQK23aEYH1tkyJZ*M13KE+?6beWLy)cqRxU4qu3EKYTydXn3QPcTpKC^_8CCIEV^Suncl@gb3O`@1;G#JG{^q>Vbm{ch0(MGQ6zD;u5rM-Q-02j* znke+d@3YF8|DLoOL-QHWZn#I@*F}+c^q{EJRHk z`?C=(V+|@{!n5rG`d5Bmt{fSV#C6_H*w?Q=*@wG_)f-5v!lGNTEvVOf=;wJkf0I$?}cv~&N$PS#yNyKwK}-^h@EnT=tf7%Od&rLT@^t~|LNf;ddQDHqp+|g=@#emGCWp%edM*X5@`)X zJu~&obC1-2(H#MQP_q!CLzPtBi^bYNL>m55{3mH?Cds{cY+1Ko8yphP&G?DWn9|4( zKJzS_+q3Fmzh$m{W0Kd0%g?7Z&tcTdELZbgN|~0&RTve*DpmAtxp)j4)$qP@SkARs zZ8Fxe8T`5Ogp+gWMjSihL*(>nj<{xb1YmpfFK-jHiQmuqXKcRuzCpP8u!Oc3t0e;2;`tR946GXAR3BJ)e<6ceF zBf}_vLeA~#Id!&Du@EjI$o3ek@FN7bXYzd-267dgB87U!W}cgp*4zYY;gt}2Mfe%g znhxN;_ms1kvu>FBG13Zke8`F8<)!4ZtW|sX&P0)UiTq;e1JCw5{FBrDhM+2S+X-c< zq0FqNE^*Ux0wIt5=>%_9N^~YAN4K`=T45`HHM*np;NldOY04c7Ew6$;+ORt-+IL!5 zh$-2xBC}dpPYTe2y0?UO3P`nU2V|`X?dSU_k}T=#Z}m=w!YfVP(s=!C`5=j_ifWPY ze#IK4Qg-!Lbbh9fgu~>AWluFhLwmez}J6u7cL2&S!2}a*p$D;F}v_Y`()Ngb9w<$+9ZVd}iaB)lXlytKtRyjFzs?OM$%L+@G-}X)!G>{U8d8y7gS(ez$!A44$*L}g!I$u~L08B2 zmVCv*qg$wx(u@W1W`=x>-dQF)Yw&;yk8O;(O*J0lOV^^tPraBVb?3s-Avds$dB>N6 z4PbVBe`8>+AYZkC6!QSneo`bzISY#go1%^|Wm;?$CL!k>AHQzlf?pJETzB(-jfv%> z8Setsd;b>(g>`H#R(<}EANmk=QyrTd@Po#9Q^xd=#w$NjNWl$+jd|g;%PxrOf zQ=-?oZ@u|pJh|`Jx?4*w3Dr`653Lyp&^6s}1#k2%y|cnj+GJ0jk~}qGn}3`5E?Q;@ z*QfmghVP5A57hduTigBOB=#9Q> z@)Qd57tzrWMqBu!hfv5PGMlEtV_UTi4Pm5AE z_!f>2QBdoCspc#dmBx#v>uO>tVv~qS9mk;}yc>JPMO?jDJLXJiF`xDXS9tR?$m?b( z<#c$^8R_L>M=W+{CrQ$V1R!+dCg9;L2E)ng6)8`U=w zzosxo*Z$Z~*tPsGY5Q&HU{(_mM#-M|TcJTP=1{AJ_;jZkCz`8rv#{?ePTV+~qAQ)1 z7)X%~&xW4%;v!K;c&R2Rb(>gSr_Q-6C**0V zjnv^fAdT*mus)kL+)Ue&#cgwpR5E)<#9x8;b!prZ9SqN)X+K>`)TaR9ypnpp276gF zb^>Za?BorDB?i5q>Q$HFAW@Tgvo0n?YQraXsBh{uKE{md71Z4Su$U+Jin_jc(M^yP z&tEP7d{!x5b~H49B`G7SyDF=crkge>foue66*cg=i9h07d7|!tp2P(?9YXG5M=hu* zM%f2Aj|B;E=wN`B+Ik~tXb-hBT04hPpV(KAau%T$ZOQE6X(i6r{_+;C^h|v)X3*1Z zQ1U!qzIUQv6D1dR$1Iv&$7Wsm`9gVn-)T^AE*`6?ygvnh4V?tG2IWqCTObG8-BCaq zg=1fU-Em4z(i3@)XL^GGny44Y^PG6-#!Lo9mxJ-=hQZPKgmF|+k+Xg~1VW7HIqwE= zR;}Rzxa)5NHM;^}tsDRqz9>zC?L#i3Xry<~$c&PM1J$+MA?2iZ;@b-S6Ii7_2)A{1 z8XxP}?euqlgHlMe7Hi!Fnq8iL@4U8ibifRc&f>|ZT|Ookf+vL@P{cto_VLW!`vqwm z`RFN{z$g_ry1JQn>Z2k7{!10OZTQz#;-S&Q6i-HR#?o${(%5voM@L0;(7ki=DE$?k z&xYUhT$8~^te?I@7Xi^-w3S$yUrQnVa+Xy64t&5q(MPoE(x_%DeXG;9Z=J{m#ZW35+nL24f)zpb{>RQ zd&7@AK3V>=NZf1!j5v@-_(H)hA>xW9z&5m5Q`` z8?$GB#_#|m;ERGXEAG+Ed+j)&AJLwS*Pv5(DO4>s@P(wBWJc@(QO1>ZUH{g~O$hgD zeaQ#tEWA))5;dIRSrMUE^yD~isX)x3;K$a((7=X>;UYf!fG@)L?{Z~5ou5~U*qzi% zq1@JBd@|S8GJb1-g44h`<*D|S^p{LSqW!jiH(^!4s=FfHKC>XsgqlDTzyo)UONo>k z%SL1`8K`Pa;X9x7Q=W$kBRMWa+fxzd>@V8@JVd`eIQPETval`*39QJG;*pTMjOWzC z%a_JMTtnRsN=F4x==!OjRS2>9`>EX^|zTwkmxWk#i8O3pN=8w(p`tPPoef+X@ z#+o*8VtuS5~mi zh2(sIqrJ69g^b|J=n!FO_|vCCI?jfYuTAko#s-ngLkpk+&aGA4Z&QOcRjQf<8cUpe zp0<)$N$#>}ne;x`N-P>VhB-^nefKI(2ePAh?npTKG@WL*Ey^KM>$UfIZq_05lD`lr>hK+~XUlaVBeEau${jUo=bUu6 z-gT&0;i-3|aS5r+sWL(VOAzB)kE3c710~94!B{tXmNAM{ckk!1Qvqk0>VjgGa)MH2 z>2~17CH_UWR=`cM;{(4)NPG>3hsFjAp+eo==C%eI18$?W3Bzl3x1d*luiLjb$sDY^ zr0mf*3JJtVX7j_=TWI{>N=BMZv^H)&puRcNCt}>KU_qt3`7{?;pA|d%ldX`Hek6f= z2Wpo2MmpEW1Les4nTgE)kO`P=FED5_wj;C|l{IAarlodUl@i$A6@FDXXFNW9kLEN; zvvd}J=GdW!`G&nz$EfRn46S%!)*|$XX#gI-_elD9)RvD;; ztLvY5_9>93`#R{?zj2-UO{r?#wsM%qQhu$&ej~!aqk8K%kl$y26L1=MWdLxmG-_}- zWp?UXctD#8vZ2Uyki-4%ncUjeI~&Rt82ifk+-kAur32Ax*_ZA#L`8WDSnbr7+^R^B zny1@QtdP06>ous>(n~MFcB#E|IhZy_SZcz3I(052x37V8Z%FK;_*lIbxDZ4~jNqocD8*S@; zfx!5VS_(0o9xpaHZieADUI=Wgf&E;8`S~jLz((cGuZx!O_3x-7J-tM!Q+NnIVJOdH z^V8X$CvOE(wh=A1g;eToC`M-{T`ZEf#DNqR^RE=WqBjh1yY_V&$ET`B5M2{E`kb|6 z(W-i7ni>m#T!b-$k%ff3AYe2Tt!ZB$nFi#D4he7d4((aa9k#DJbmn*F&*uR_7=!kS zIriLA+2Y%b&q7QaSN6!Yd?v-L@uT+BSgvjCq*judydo^6HLOJIHcL&X;@PKZar+8A zzHB?rV=2=nHn57;HgIV-A`O2)5>p)$@_$KXw4}v}2>=WZvFo^J^Gr+ZKQe^O;Pz-{Q3z*{b`~1jFsURX?z`wj(^ZluUBB5 z)0lFyaoEky=zB77;gk9^D*Wx4@0xfhV#QksiF-sFW4Zb*l*lB4zSN11GP?7vBCVQk z_mrEtMzzs~Yn24}c7s{Xeg?>zyfni*GNQzWkcs&pq)uDgx3Xbw!CE;6g3|}<ert^IHOounR(Bn+bL-K^W?Mwc;`IBnlFjz%*eu%(e*QwUHukS*zzff$ z#nTzF%>4s!QmA7=J>@~N#YTldpA_nU9oFn2N~pP`w$!l3*m-=eaCLe{W$+U$EYX7c ztM5p2za%ve{h`e&&ehfvVkiEhH!)%dvl9C@(WB&5593A##erjzp%V%nAJ(ap?JT~5 zs>{?%9J`Uns@#=5?#+Jm75S2FE42F5uK+|!0<1UNl9`U}-@in0)Qd{2>*@16Pj>pSdqAy!gL1fGj*(FEf_ktV&&D;QX9`)wn*H&!_e#S zTRT+XsMEDI*I$3*e! z5N8MkcN}?4CO1E3N;KxF-*{7hGAz=V=<*}YZS$stP{Qy!i;;zyz=_cAd-H~ofyKl{ z_ryiI{MU2B^`sG2)g$kpi4Jj8mV!Qbmyu%oL3xdY!?92=NhEztd2SrjKjUI!Z%zNOZK3G3_p6>svsHpxWkrdgiBo-k z(K=2>B!fLLdCctxTS7w4d1-OG$v5q?<~l{pv>US1lhH0Uhc_I5Iy4Qj3Tbx1mFXRq zKkE?a#r#vP+#5e{NWV=A`p!rHSPFqhh)TB9@cyqubAVJ zw*Ki#mG=?LCezTi?43KC(7XAi4M2`@KdTb8mLH$gfi&##Gp%XYQw_ZNtiD`eVCRwRPimJ%{c+ zyYUg}^}u4AiCNq*j(l&gmVkS!e1^~F!dq5%;WD&|UtMcOl%qd$&kR|NBc77$&9_kD zOKp4&p0ziBa(ksWV^o%s_!b35%j`%zZ;iw6@9rUTQj1JdR!Pf&%ha>{gYLz4sHu;G z*FF~=HKU5q%rkSyk<8c6dr>mHsTMOeyX`+&ioBvRKD4)KQt_QoH6EPXFU4fX^wBo9 zS7NKBdNOR$ku1x94VtUfRGN*gcN^BCESKiwKYP9j zBeIH(1a9;%e}QBneZvK_dW7%} zSUOT7^q0n>Ny3%6Q!zBgD?zm@Kx&HFVVP5MwLWf;>62r^pbg?3~uai5{sQbZYnVk zaYW0JVBGmp;bA)LeR``gz+mskT3mdSmzQz=8DgUpMc-e0+x~?vOKN23wB#9f$^3Wl zJ#N;6UeZI9kmQ7K?acmG5&?>bUybHRIMK*|&$HF5_(JPEP1H(e{Ofhi5r}?EXlbao zzp04$v5o>n>Rj+=UwF?G;PgHgbzs)#vp|9rkG+rn+WtWDxk;We4^@^jhf`@QJjRM# zqJ_kqveTC8N~)GNY++&by}`0(;Dx~F>|b-3i}`KM!9OAB^YMJ@`A(9twbWaTcmYy> z(G*$P5qLG9$UsOhPF7Qt^d}upK%@CiPX$@xl$O~WV{O5p&mgm`sVFfC=w6iy@dPA| z-!~|;_?TbG5}MRmmzghWq(&|B9l}d}7tWI=q0&Fqz$j19KxKY46R5_F*N#rC$|Q=t zzUl#_QybZW$SiE5n9euqZ=buMs)k~JC+vZgT$lR-E{FVTg6;I<2&eAKI513nP?>8a zXqITOJ9E@DU$*W&<>;hZ9w)i{-a&wv+^6_3e$PSoMft}XgsJ~+ZRV8*B}dDTABbEQ zN=khPU3j4f8PP40NWt}p{DP(?QMpXn1rf2}ptRtvAMhV??0> z(GBxRM)Eg9dTZeX6Em_4w_|b4IFd`(HI~Hgt$P?vUmMMhKdm*mmg}vE3fW{C=^h_{^G!i5v)<&j0;28 zRi;+sEc@QgV0zusz*aNanSYqHKhUaj^dl|&>J>wdO{o;RGC#0vRjR9hdD%ryI;&r! zWZ5__<}-2kTr%ZpydayMcS)~k#GyQ4i-9h^26@6I8M`EbJP9{iSITAT>W_CGxD4DR z)0!Xi5t=evb9_HX4x-L|Jj-L&mD?RWC8a^5qY5Ij)DCD|ygeO!g5`-_F_h#ntQMy! zZ(b^+9H-%$(#6ojt|O*@mX%SsGxMD4Kp6CUT?%EJi+l?wYp8lP+%W_={Oy zE=z=+3F-19w3%k(Kw&YRAGEn3?N2*Tqo3X(>8QiIv-DgVBx4({!{>vZiGru>5kQNg zy2ttR+h|<%8~?jceNU?Pl;8W3JKf+Z_ldXNaG8nv_}hJ>$o1KOs3r)C?-4ER9Uvt; zSeu{Y&|kqzGgZ_S#w8klElWvZ2TLl)$J(unrz~O2yB;`=;B&9a1>xe}>bb^OyO;T% zWf`T`E|yH55*&%9X#26(tIJ}Zu$b*A-PjDs??irTU6GG4+BKu`CC9?X3zK?9XL&~kLwy)>GW0Y01-Xx&^#j-EKh>OGofYdo%XL6QLO*17K|JRfIv62S2zg4q zarAP!+j}iTSJ~%@yI*IxRiU4GOPjke(Q*DwW;(jkSdosIED!sjGb8GKZM!ETe#^J4FbP0A6J9chreJCnb-RBvg`&6y%)Z*^b2?txdmfCgTCDg&ycPLFqnCnYd>kcs z`tU2Fq>*|&n}Jb;z+ruLv2fRTNy64xrj|7&RiKoP$R}h>u`dr{&d1+tUTYRTIrI#7 zeZaK13a1^8N+>z4!2PlU(2AcYG$(X@&E*X_fprpp@J-mCb-##r92=t`COsYd;E6mw zs=mEPaI6Th;9K6$V0OSEjVanB_=u_&`i}R9yAa{2;s;t5_t(=gEhjIwA(-_AFZx|I z*=;>H4S8SsRoTb#b&W5kND#_7uag~Dh&|=+!j;|Pa>L&gm(%84X&TwNlQuIqT*x_R z=W`%`2$0F04tN<&Yo7+~OCAGD;x>l=VD{Sio-Z7n!8ljc#q}D`aqSc4##>&|qQ*fVh5>iiv%M45HtumA1rxU1Ebj-Ao956bn3%Y^ zor62uTDX0cjXAU-_3i zg}qa#mOg_MFsmUtEpDTzN&2a{|F5}&Qw7Zce!ri6fh zUv*?xYgnO4lG!KE5;N)LFik2ud#+Es^B|>&Nyou z)Cswq;`?~z=N>jLLtf`6u`|(vS2{I=%e2uT1W%qI9JrO`w>``G8L57OF5theWo_zQ zGt!!>HjAn_G(2Q)2`3C;Ru@qwcQ{vn;pd%P`j~^v`@GcSv|Cq}qxYf4`Os`(z}awZ zDDbs_r@#U-b-{G!(2AxTBU&NElU8k@{2hkcR@x@6Gqh@=z7LRlIWWVg;^j=oAvr$ zx%QVlGTn+rR;MG0Ss+=pgZI7#28AMJ9HLvF>dgu>4|tFEUl_^MnD>(Xi-aDgI+9vIKL@S3Vbeye+(yG%k8nzFJXdM>^ zPCE+eg?te7#NVG6V?&Mp)GYRPOTXdJ^Nj-k+sZdMhacX*r>C-R7)`-&zP~RgF%b_P zClbB4!p>e69c=7ToQCS2zwp_R%=Nf9AAs@2qt^z=f?938#>mQXS)v<%OWUI|1PtwQ zEUjo$M%tNGur3k&ulZ61K6Wz)QK!7bV`(XO1TUf0SM}fbvN?5hSfb76`GCcAU0aGJ zettXrR;wzbCPZ9j)${Vu(x7De-1R9`Zmi})Q==FClbw{M54qEN+z)Ee@pO?2T=FmD z28m9HkOs>C14gP@Tnc4}I^MTEf@%D6(hLhDc{ zSs1SiB?{|op?u*DJav?&yB)ef$JNy7$hO@Y`PY^8sl zaF7wO1UL-**KmeWg%uL68qRWNWK03zt}wQ6Q3JdyyvBPV!n;_QSOfw4tSo|xp|H9L z*2C2nAx2wRUqlsv@K{6>-N*sCu(5DJr%+ZDz>pMhWZ^taJkSMLrHDX~mxYunXgCoj z7TEUyEQ7LeHn1z4Dm(xyEzC&6Bfx(L>$qa!d=Q^JPTvLH=BVvGPW0*tY6i~{t=;W-E# zKvEGfm?*{Q@#9hX?#ESCl&{C<)u?`XRo}(VtAO=}a>nJ~(0|Djh#r*VbwJJUXbI_;!b}^qZHC2Fki^d zje~U_877HCOF)svj0rwFO1}fC*ydeZHHwab9(A;3!XVaF4GNSi?+ zq?fSbBaV>x6kMXVfG!FP#|IGs2v;Cy2S~#LH#y^xK*?idhHW|vc{E6;(y)^AyJ5^E zXaPQK8Z*@xr|OZdFMofU@X(x}MvYtt8-bE&fhw?F6_HX4vefq)8Bab2<$!OfUqFaF z@`{frN^cd(!=r4Z#JN%+Xegy{uhys&g46-mt+7bEvZy_$Q6-*qPA~~GVMy~JaaR;5 z2|6X6OCwZKvbEq$teMCe>NA&^Cd9_0{`0w32HDMEcvo6?VkdvB#S@SfhRNX>#z?b+ zrDvp3??ABs>`+W_H4$tK!w#?r0Jb)#EpRltcL_ng0V}}j2OgB8D^#`n#5qmbbP>b>`-9v@TNG&Nk&6P*FSebV4 zYT8GuPPlp{S~P#=w6lh02xtIy9LmQdMr$G&1O|a%P@;iWJY+8C95STY+~7Hc8dxP$ zVW)>;m{>L_O)U^dKyz}QL}5818iNP{0o3#8OA<#YjOEb95cxcMQqBFx9H$gP*Qi@J z-e5>g>rnhgthm7nDAEdK6@5brhboqWVkQd7`z*povN3;b6|fPpWJeJ6z#qM(-4r|l zKiG_Q6rnLDMMGAlFfc+>fw2?Lpm+`Lz|YZ0JvmV`wq?g_-lHiyjhvygg84yAgXj(Y z4t5k?kg`baR--dmFj&}UQ8y-_AxH&PPzZuTE!ZpAjQ2>KqLJqo48J?O%?Q#Iv}L<( zw?In6=GA}odCf`&+}u0%I~q5FDr_GKk~}2&341(^n2s%30gntptSZ!f8Z#pKG#@eL zIO`~ZQKI$a@rtu0CSZ(Ga1J;X6jU@i#^ca8QOVHI3bGnOSpkjNMI<9gB9v^=2_ny+ z>pDe2U1Cs*fsZ1~apZv%jB{2ujG4%cLEGz(Yw3S$5)qH8HXU0F9%sBB$bLOE_cQ|o z%rKODGj+JT{xienD0uVvo{8kUw@QLdb)6~PNSy?@H0vFjOpev z!6)4H8JRv;(EeP zBg1eNna>+IV_ZF^6B-=BQ2F0ry3d%?#UpW{DWK4Z_slT{A-Cog$}6rrq4v?93Ngx3 z1ZkwIy(M_eDvlhA+>k*KEpMj-sGD)vMl*kUhB0iQRy5mabWzo`=X~^e;uS>9`f6sm zrg1gbzY~o`ux@By4a?7ymMm7wJa3*Qi^DlaCqreZ!Z?btGpubfVwANW*R72F6rz|S ze^(UhVNGMMcL8k)<1beIj#?zt5oPZXsk{(ruhEfmwrpC&jTi1``aREm#vJER7w>;0 zS`T+nSTm!%s65q%;u6{&W*196r%BDD6^`OAc&v?7b89N07(`I^bJQw=o#BCzYEUd_ zJk;)%s!z5o#+iCV|ETu484UIa+vf5;T8EO=9BOo(eFH<)kH*8jUr!qjJ$w1;A|@_s-T+ zpJ1rnizq~6Fdd?>N$}Y9IPMJS6i+8|Ium5Z(coKpWW<9pNuN%^HGrf(z;dUjuzq#i zVJLJ*;~YI@s<2$d{o0YylUsFaCK^;$IohlgyUm~aN!5vxiz%(66MISBhvI*rMGBtv7$gIMEYZPmqFjDDe9J+%ENU6baO?uRv_#j_0oP!QIP(TmsfJKd{6`G1M z4qg)1$vjh=;7OMR>FE`cCAWVt+?_k_!YTX)40u?Jf~+yga#6-rD|O6Vf{^_4y0$k@ zWmQ%e=OqK{S94jq!3q$6Q%L66@jg;Vof2~1rZgtU$27AXDp!xq>U?^Zk#|?K&5bv! zS=Qn~(ec|XW{ZLY8>F+ym}f-+TvA`gC%7o|)U^tDjrJz&Tp_0yAA~x;u7Fj2yj_>y z4Wr#zCIUj1w2{p-?a^vWr>en`RAJlD456t6*JN2qA=%SpNm;p1CxcyTsi@q!@0Xqv z15TSM_PO$Mg}r4E`m&zM}E(}=Hlpgh4h0+cp|;#YrEi<6nq8zaLNRC;|0 zHi^N5o_nrtX?L0js+XjoUe%{})@yu`5^%nPg%JcHL+Z7d82j^lg)tlIJo35nSiXPc zFljHZQ#|Dre0L|;dvlxwgPj7tc#WxRkonUu*{75EIpVWwdAgWft>z2| zXXRn}rW}`NWmTS+lk&2hmGg4(Z>DH(0(eZwKO+4ZgRkEFad5oP45|SJH8O~ezK##W zpooJEJ}I99lh4b&avxZI3CtdqAAsT0n&m&sX<~amtIDhLs#;9uXXTDgp?H68ybyZ2Z4 z0sH!d`sQTS%&YIse09b1`pkdkZoP$W)_dp%HZ89&FUPCNe3p)FDM)Rjpu;Dx4-Y?q zf)4kG6a;N3d!{hqCx)wUGTcY#kgO!H29@*#?Dj=%v+u}iPeAm*Fxc)_wA=CJ$=SG^ zGCN?x$#kzSCgs&~0>)dd#$dSL%5T*vu-g)j*~+UAetWb3I|;{bl#PG(Ni^ET2NF#b zgQ7w9=fozl*y77=!`H_bFJ3>{Bfk7*oy&DF_VR)F3UZLIm!Rghk}n`Y{dbe9=jaEk z3HsaDuZ<`&ZFi#8R#E)r?yFbFMuualDgf2+!`^Xmnqii`5%$JIqOTdM@J#^v~Q0Z8EA1ae26LDo?%i&eF6POq2c zU*%ud^Hs$UFwrcl%L&u7tiHhs%gOgn^7#ho=X?31{HOe}NlZROV#C(_ats?vZ+XN*0)8pgaBR1** zXV;pIIvR@~|K6H)JHWctIqRZbc%dyHfYlz?<<~{48d;#$Tik9HvXA$Vo;^K4Zrwn+ zvG!mY^@Aun=s`t;U$i)X*(yjMy*l21{~9?NP@7sOAF0%XaN>t!5mu9%BGBs;gFcWx zJAtMSF)6rgEIV%q8vR=I8~8)d8@3AZyL*SbAILt!fHSQ75aGqcNYxEWm1AE;qPbx2 z&v{=?$IFXyl{4$NihA$=*#F}I>zB8m2CMH_Jfh+o`YtHGZNn~q@P4;UZ-y(F@nSLm zn{|vCs;YBUN(x{plW*f`H9M`!%gGFK#^v>NHMyGpC}-o7={QLTX{Z(zq?N@SR+iT% z%W73#&mfX6mg8wlU|Yq?AZ3eJysddvN^M{4JWE9mwyW!4(zL3a7%&<@1ZZ)* zRY*O2fAHJk$6wy=9UbC_4+9zw;1U|Y(QGV!r^`2N54EF1y{Hv~iSXolI<4BSy&;Z0 zpPkOnCbM%m3hG@@_YA+}{aI?45Pa#g4u@Z#tmqaOLCU*-<7GwP@$S8cX-NuyfI?p$ z-G{fHB`r>5N=bcpHMFlOeO7`$1Va9>Yxds^BB zeutL61&fxNL?am4B_`e>6u0Z*HZlGufmk-~q8{IdXJA~nWMb#mwZ^yc2J36zS6mQl=PoaBR zJ#L82_?x&T{wVGmKQS1UzfaDfMB-naN&X3&;@_qzy6cooPsdQX&Cf@SC;MgG{{RRL zcIgV2UmF8!QZP0`HZ?XhML9t+Lqan}Mno_;HbOH+I7LJ>L_|b3H$EUdFg8LqH8wOw zIYBW)LNi51L@+ouLNi4;MMN}2L_{_>K3xhgOl59obZ8(lIWZuUAqgmdCD6G_gFzIA z;q%RJG>azgF^TK%7S}B9UW2WL5-V#{5g>*h(IJJtQ^UbS&||1WRWb9IkG|);DySzZG(8kK|*1hn;`+x0*Rm5 z+&amVWzt4kA^w%Eu97u>5+orKCJ~Y*?WBW5NsM%oF49eUNH6Il{Ui>_?b;m!kV5m# z3CQWYRT5Ge1!ulXfIK<2<|7&zIt8UZdMR--_rYBdIYWvzU;T-s_JZhOaS z0`9PGH3?Tdx0-^x_N}Jj%4vbC>?sdLrFyG8-77DTO6^stCzQq`mhy3|e3q2&n(~uW qej5VM-6(-eJ9HCX@J@%Ry!-=H+cl1tXB`6{3N<${3MC~)PeuwCM=WLl diff --git a/docs/derivations/laneEmdenVariationalBlockForm.tex b/docs/derivations/laneEmdenVariationalBlockForm.tex index 589cbaf..77e0b36 100644 --- a/docs/derivations/laneEmdenVariationalBlockForm.tex +++ b/docs/derivations/laneEmdenVariationalBlockForm.tex @@ -14,7 +14,7 @@ \title{Deriving The Full Lane Emden Weak Form} \author{Emily M. Boudreaux} -\date{March 2025} +\date{April 2025} \begin{document} diff --git a/src/poly/solver/private/polySolver.cpp b/src/poly/solver/private/polySolver.cpp index 4b24351..8f290ab 100644 --- a/src/poly/solver/private/polySolver.cpp +++ b/src/poly/solver/private/polySolver.cpp @@ -34,6 +34,8 @@ #include "resourceManagerTypes.h" #include "operator.h" +#include "debug.h" + #include "quill/LogMacros.h" @@ -149,14 +151,14 @@ void PolySolver::assembleBlockSystem() { // A full derivation of the weak form can be found in the 4DSSE documentation // --- Assemble the MixedBilinear and Bilinear forms (M, D, and Q) --- - auto Mform = std::make_unique(m_feTheta.get(), m_fePhi.get()); - auto Qform = std::make_unique(m_fePhi.get(), m_feTheta.get()); + auto Mform = std::make_unique(m_fePhi.get(), m_feTheta.get()); + auto Qform = std::make_unique(m_feTheta.get(), m_fePhi.get()); auto Dform = std::make_unique(m_fePhi.get()); // TODO: Check the sign on all of the integrators - Mform->AddDomainIntegrator(new mfem::MixedVectorWeakDivergenceIntegrator(negOneCoeff)); - Qform->AddDomainIntegrator(new mfem::MixedVectorGradientIntegrator(negOneVCoeff)); - Dform->AddDomainIntegrator(new mfem::VectorFEMassIntegrator(oneCoeff)); + Mform->AddDomainIntegrator(new mfem::MixedVectorWeakDivergenceIntegrator()); + Qform->AddDomainIntegrator(new mfem::MixedVectorGradientIntegrator()); + Dform->AddDomainIntegrator(new mfem::VectorFEMassIntegrator()); Mform->Assemble(); Mform->Finalize(); @@ -187,47 +189,24 @@ void PolySolver::solve(){ // --- Set the initial guess for the solution --- setInitialGuess(); - // --- Set the essential true dofs for the operator --- - mfem::Array theta_ess_tdof_list, phi_ess_tdof_list; - std::tie(theta_ess_tdof_list, phi_ess_tdof_list) = getEssentialTrueDof(); - m_polytropOperator->SetEssentialTrueDofs(theta_ess_tdof_list, phi_ess_tdof_list); + setupOperator(); - // --- Load configuration parameters --- - double newtonRelTol = m_config.get("Poly:Solver:Newton:RelTol", 1e-7); - double newtonAbsTol = m_config.get("Poly:Solver:Newton:AbsTol", 1e-7); - int newtonMaxIter = m_config.get("Poly:Solver:Newton:MaxIter", 200); - int newtonPrintLevel = m_config.get("Poly:Solver:Newton:PrintLevel", 1); + // It's safer to get the offsets directly from the operator after finalization + const mfem::Array& block_offsets = m_polytropOperator->GetBlockOffsets(); // Assuming a getter exists or accessing member if public/friend + mfem::BlockVector state_vector(block_offsets); + state_vector.GetBlock(0) = *m_theta; + state_vector.GetBlock(1) = *m_phi; - double gmresRelTol = m_config.get("Poly:Solver:GMRES:RelTol", 1e-10); - double gmresAbsTol = m_config.get("Poly:Solver:GMRES:AbsTol", 1e-12); - int gmresMaxIter = m_config.get("Poly:Solver:GMRES:MaxIter", 2000); - int gmresPrintLevel = m_config.get("Poly:Solver:GMRES:PrintLevel", 0); + mfem::Vector zero_rhs(block_offsets.Last()); + zero_rhs = 0.0; - LOG_DEBUG(m_logger, "Newton Solver (relTol: {:0.2E}, absTol: {:0.2E}, maxIter: {}, printLevel: {})", newtonRelTol, newtonAbsTol, newtonMaxIter, newtonPrintLevel); - LOG_DEBUG(m_logger, "GMRES Solver (relTol: {:0.2E}, absTol: {:0.2E}, maxIter: {}, printLevel: {})", gmresRelTol, gmresAbsTol, gmresMaxIter, gmresPrintLevel); - // --- Set up the Newton solver --- - mfem::NewtonSolver newtonSolver; - newtonSolver.SetRelTol(newtonRelTol); - newtonSolver.SetAbsTol(newtonAbsTol); - newtonSolver.SetMaxIter(newtonMaxIter); - newtonSolver.SetPrintLevel(newtonPrintLevel); - newtonSolver.SetOperator(*m_polytropOperator); - mfem::GMRESSolver gmresSolver; - gmresSolver.SetRelTol(gmresRelTol); - gmresSolver.SetAbsTol(gmresAbsTol); - gmresSolver.SetMaxIter(gmresMaxIter); - gmresSolver.SetPrintLevel(gmresPrintLevel); - newtonSolver.SetSolver(gmresSolver); - // newtonSolver.SetAdaptiveLinRtol(); + mfem::NewtonSolver newtonSolver = setupNewtonSolver(); - mfem::Vector B(m_feTheta->GetTrueVSize()); - B = 0.0; - - newtonSolver.Mult(B, *m_theta); + newtonSolver.Mult(zero_rhs, state_vector); // --- Save and view the solution --- - saveAndViewSolution(); + saveAndViewSolution(state_vector); } @@ -295,22 +274,77 @@ void PolySolver::setInitialGuess() { } -void PolySolver::saveAndViewSolution() { +void PolySolver::saveAndViewSolution(const mfem::BlockVector& state_vector) { + mfem::BlockVector x_block(const_cast(state_vector), m_polytropOperator->GetBlockOffsets()); + mfem::Vector& x_theta = x_block.GetBlock(0); + mfem::Vector& x_phi = x_block.GetBlock(1); + bool doView = m_config.get("Poly:Output:View", false); if (doView) { - Probe::glVisView(*m_theta, *m_mesh, "solution"); + Probe::glVisView(x_theta, *m_feTheta, "θ Solution"); + Probe::glVisView(x_phi, *m_fePhi, "ɸ Solution"); } // --- Extract the Solution --- bool write11DSolution = m_config.get("Poly:Output:1D:Save", true); if (write11DSolution) { std::string solutionPath = m_config.get("Poly:Output:1D:Path", "polytropeSolution_1D.csv"); + std::string derivSolPath = "d" + solutionPath; + double rayCoLatitude = m_config.get("Poly:Output:1D:RayCoLatitude", 0.0); double rayLongitude = m_config.get("Poly:Output:1D:RayLongitude", 0.0); int raySamples = m_config.get("Poly:Output:1D:RaySamples", 100); std::vector rayDirection = {rayCoLatitude, rayLongitude}; - Probe::getRaySolution(*m_theta, *m_feTheta, rayDirection, raySamples, solutionPath); + Probe::getRaySolution(x_theta, *m_feTheta, rayDirection, raySamples, solutionPath); + // Probe::getRaySolution(x_phi, *m_fePhi, rayDirection, raySamples, derivSolPath); } +} + +void PolySolver::setupOperator() { + mfem::Array theta_ess_tdof_list, phi_ess_tdof_list; + std::tie(theta_ess_tdof_list, phi_ess_tdof_list) = getEssentialTrueDof(); + m_polytropOperator->SetEssentialTrueDofs(theta_ess_tdof_list, phi_ess_tdof_list); + + // -- Finalize the operator -- + m_polytropOperator->finalize(); + + if (!m_polytropOperator->isFinalized()) { + LOG_ERROR(m_logger, "PolytropeOperator is not finalized. Cannot solve."); + throw std::runtime_error("PolytropeOperator is not finalized. Cannot solve."); + } +} + +mfem::NewtonSolver PolySolver::setupNewtonSolver(){ + // --- Load configuration parameters --- + double newtonRelTol = m_config.get("Poly:Solver:Newton:RelTol", 1e-7); + double newtonAbsTol = m_config.get("Poly:Solver:Newton:AbsTol", 1e-7); + int newtonMaxIter = m_config.get("Poly:Solver:Newton:MaxIter", 200); + int newtonPrintLevel = m_config.get("Poly:Solver:Newton:PrintLevel", 1); + + double gmresRelTol = m_config.get("Poly:Solver:GMRES:RelTol", 1e-10); + double gmresAbsTol = m_config.get("Poly:Solver:GMRES:AbsTol", 1e-12); + int gmresMaxIter = m_config.get("Poly:Solver:GMRES:MaxIter", 2000); + int gmresPrintLevel = m_config.get("Poly:Solver:GMRES:PrintLevel", 0); + + LOG_DEBUG(m_logger, "Newton Solver (relTol: {:0.2E}, absTol: {:0.2E}, maxIter: {}, printLevel: {})", newtonRelTol, newtonAbsTol, newtonMaxIter, newtonPrintLevel); + LOG_DEBUG(m_logger, "GMRES Solver (relTol: {:0.2E}, absTol: {:0.2E}, maxIter: {}, printLevel: {})", gmresRelTol, gmresAbsTol, gmresMaxIter, gmresPrintLevel); + + // --- Set up the Newton solver --- + mfem::NewtonSolver newtonSolver; + newtonSolver.SetRelTol(newtonRelTol); + newtonSolver.SetAbsTol(newtonAbsTol); + newtonSolver.SetMaxIter(newtonMaxIter); + newtonSolver.SetPrintLevel(newtonPrintLevel); + newtonSolver.SetOperator(*m_polytropOperator); + mfem::GMRESSolver gmresSolver; + gmresSolver.SetRelTol(gmresRelTol); + gmresSolver.SetAbsTol(gmresAbsTol); + gmresSolver.SetMaxIter(gmresMaxIter); + gmresSolver.SetPrintLevel(gmresPrintLevel); + newtonSolver.SetSolver(gmresSolver); + // newtonSolver.SetAdaptiveLinRtol(); + + return newtonSolver; } \ No newline at end of file diff --git a/src/poly/solver/public/polySolver.h b/src/poly/solver/public/polySolver.h index 92ceec9..eb89cf1 100644 --- a/src/poly/solver/public/polySolver.h +++ b/src/poly/solver/public/polySolver.h @@ -1,6 +1,7 @@ #ifndef POLYSOLVER_H #define POLYSOLVER_H +#include "linalg/solvers.hpp" #include "mfem.hpp" #include #include @@ -53,7 +54,9 @@ private: // Private methods std::pair, mfem::Array> getEssentialTrueDof(); mfem::Array findCenterElement(); void setInitialGuess(); - void saveAndViewSolution(); + void saveAndViewSolution(const mfem::BlockVector& state_vector); + mfem::NewtonSolver setupNewtonSolver(); + void setupOperator(); }; diff --git a/src/poly/utils/private/integrators.cpp b/src/poly/utils/private/integrators.cpp index 8c7e1a5..7590604 100644 --- a/src/poly/utils/private/integrators.cpp +++ b/src/poly/utils/private/integrators.cpp @@ -60,7 +60,8 @@ namespace polyMFEMUtils { double u_safe = std::max(u_val, 0.0); double u_nl = std::pow(u_safe, m_polytropicIndex); - double coeff_val = m_coeff.Eval(Trans, ip); + // double coeff_val = m_coeff.Eval(Trans, ip); + double coeff_val = 1.0; double x2_u_nl = coeff_val * u_nl; for (int i = 0; i < dof; i++){ @@ -94,7 +95,8 @@ namespace polyMFEMUtils { for (int j = 0; j < dof; j++) { u_val += elfun(j) * shape(j); } - double coeff_val = m_coeff.Eval(Trans, ip); + // double coeff_val = m_coeff.Eval(Trans, ip); + double coeff_val = 1.0; // Calculate the Jacobian diff --git a/src/poly/utils/private/operator.cpp b/src/poly/utils/private/operator.cpp index 4de7639..7dd3815 100644 --- a/src/poly/utils/private/operator.cpp +++ b/src/poly/utils/private/operator.cpp @@ -3,6 +3,8 @@ #include "linalg/vector.hpp" #include +#include "debug.h" + PolytropeOperator::PolytropeOperator( std::unique_ptr M, @@ -19,8 +21,12 @@ PolytropeOperator::PolytropeOperator( m_Q = std::move(Q); m_D = std::move(D); m_f = std::move(f); +} - +void PolytropeOperator::finalize() { + if (m_isFinalized) { + return; + } m_Mmat = std::make_unique(m_M->SpMat()); m_Qmat = std::make_unique(m_Q->SpMat()); m_Dmat = std::make_unique(m_D->SpMat()); @@ -28,14 +34,14 @@ PolytropeOperator::PolytropeOperator( m_negM_op = std::make_unique(m_Mmat.get(), -1.0); m_negQ_op = std::make_unique(m_Qmat.get(), -1.0); - MFEM_ASSERT(m_Mmat.get() != nullptr, "Matrix m_Mmat is null in PolytropeOperator constructor"); - MFEM_ASSERT(m_Qmat.get() != nullptr, "Matrix m_Qmat is null in PolytropeOperator constructor"); - MFEM_ASSERT(m_Dmat.get() != nullptr, "Matrix m_Dmat is null in PolytropeOperator constructor"); - MFEM_ASSERT(m_f.get() != nullptr, "NonlinearForm m_f is null in PolytropeOperator constructor"); + m_isFinalized = true; } void PolytropeOperator::Mult(const mfem::Vector &x, mfem::Vector &y) const { + if (!m_isFinalized) { + MFEM_ABORT("PolytropeOperator::Mult called before finalize"); + } // -- Create BlockVector views for input x and output y mfem::BlockVector x_block(const_cast(x), m_blockOffsets); mfem::BlockVector y_block(y, m_blockOffsets); @@ -89,6 +95,9 @@ void PolytropeOperator::Mult(const mfem::Vector &x, mfem::Vector &y) const { } mfem::Operator& PolytropeOperator::GetGradient(const mfem::Vector &x) const { + if (!m_isFinalized) { + MFEM_ABORT("PolytropeOperator::GetGradient called before finalize"); + } // -- Get the gradient of f -- mfem::BlockVector x_block(const_cast(x), m_blockOffsets); const mfem::Vector& x_theta = x_block.GetBlock(0); @@ -111,6 +120,7 @@ mfem::Operator& PolytropeOperator::GetGradient(const mfem::Vector &x) const { void PolytropeOperator::SetEssentialTrueDofs(const mfem::Array &theta_ess_tofs, const mfem::Array &phi_ess_tofs) { + m_isFinalized = false; m_theta_ess_tofs = theta_ess_tofs; m_phi_ess_tofs = phi_ess_tofs; diff --git a/src/poly/utils/public/operator.h b/src/poly/utils/public/operator.h index 00d02e9..b809280 100644 --- a/src/poly/utils/public/operator.h +++ b/src/poly/utils/public/operator.h @@ -20,13 +20,19 @@ public: void SetEssentialTrueDofs(const mfem::Array &theta_ess_tofs, const mfem::Array &phi_ess_tofs); + bool isFinalized() const { return m_isFinalized; } + + void finalize(); + + const mfem::Array& GetBlockOffsets() const { return m_blockOffsets; } + private: std::unique_ptr m_M; std::unique_ptr m_Q; std::unique_ptr m_D; std::unique_ptr m_f; - const mfem::Array &m_blockOffsets; + const mfem::Array m_blockOffsets; mfem::Array m_theta_ess_tofs; mfem::Array m_phi_ess_tofs; @@ -39,6 +45,8 @@ private: std::unique_ptr m_negM_op; std::unique_ptr m_negQ_op; mutable std::unique_ptr m_jacobian; + + bool m_isFinalized = false; }; diff --git a/src/probe/public/probe.h b/src/probe/public/probe.h index d9157f6..dfe3427 100644 --- a/src/probe/public/probe.h +++ b/src/probe/public/probe.h @@ -25,7 +25,6 @@ #include #include #include -#include #include #include "mfem.hpp" @@ -65,7 +64,6 @@ namespace Probe { */ void glVisView(mfem::Vector &vec, mfem::FiniteElementSpace &fes, const std::string &windowTitle = "vector", const std::string& keyset=""); - double getMeshRadius(mfem::Mesh& mesh); diff --git a/tests/testsConfig.yaml b/tests/testsConfig.yaml index 0032269..83cbc2a 100644 --- a/tests/testsConfig.yaml +++ b/tests/testsConfig.yaml @@ -3,8 +3,8 @@ Debug: true Probe: GLVis: Visualization: true - Host: "10.28.92.45" - # Host: "localhost" + # Host: "10.28.92.45" + Host: "localhost" Port: 19916 DefaultKeyset: "iimmMcaa" GetRaySolution: @@ -25,6 +25,7 @@ Poly: Constraint: Gamma: 1e2 Output: + View: true 1D: Save: true Path: "output/Poly/1D/poly.csv"