From 1d3feae72f4ff2c0909f9946757398ff01970fd2 Mon Sep 17 00:00:00 2001 From: Mark Messner Date: Wed, 10 Apr 2024 12:54:11 -0500 Subject: [PATCH] Okay have it both ways --- .vscode/settings.json | 6 + include/kernels/ADStokesPressure.h | 25 +++ include/kernels/ADStokesPressureByParts.h | 25 +++ include/kernels/ADStokesStressDivergence.h | 3 +- src/kernels/ADStokesPressure.C | 35 ++++ src/kernels/ADStokesPressureByParts.C | 35 ++++ src/kernels/ADStokesStressDivergence.C | 7 +- src/materials/StokesLinearViscous.C | 4 +- src/materials/StokesPowerLawCreep.C | 6 +- src/materials/StokesStrainRate.C | 3 +- tests/stokes/flow.i | 155 ------------------ tests/stokes/gold/flow_out.e | Bin 157288 -> 0 bytes tests/stokes/{manufactured.i => patch.i} | 4 + tests/stokes/patch_by.i | 182 +++++++++++++++++++++ tests/stokes/structural.i | 149 ++++++++++++++--- tests/stokes/tests | 20 ++- 16 files changed, 460 insertions(+), 199 deletions(-) create mode 100644 .vscode/settings.json create mode 100644 include/kernels/ADStokesPressure.h create mode 100644 include/kernels/ADStokesPressureByParts.h create mode 100644 src/kernels/ADStokesPressure.C create mode 100644 src/kernels/ADStokesPressureByParts.C delete mode 100644 tests/stokes/flow.i delete mode 100644 tests/stokes/gold/flow_out.e rename tests/stokes/{manufactured.i => patch.i} (98%) create mode 100644 tests/stokes/patch_by.i diff --git a/.vscode/settings.json b/.vscode/settings.json new file mode 100644 index 0000000..3fe7657 --- /dev/null +++ b/.vscode/settings.json @@ -0,0 +1,6 @@ +{ + "files.associations": { + "ADStokesStressDivergence.C": "cpp", + "ADStokesIncompressibility.C": "cpp" + } +} \ No newline at end of file diff --git a/include/kernels/ADStokesPressure.h b/include/kernels/ADStokesPressure.h new file mode 100644 index 0000000..e1fd12c --- /dev/null +++ b/include/kernels/ADStokesPressure.h @@ -0,0 +1,25 @@ +//* This file is part of the MOOSE framework +//* https://www.mooseframework.org +//* +//* All rights reserved, see COPYRIGHT for full restrictions +//* https://github.com/idaholab/moose/blob/master/COPYRIGHT +//* +//* Licensed under LGPL 2.1, please see LICENSE for details +//* https://www.gnu.org/licenses/lgpl-2.1.html + +#pragma once + +#include "ADVectorKernel.h" + +class ADStokesPressure : public ADVectorKernel +{ +public: + static InputParameters validParams(); + + ADStokesPressure(const InputParameters & parameters); + +protected: + virtual ADReal computeQpResidual() override; + + const ADVariableGradient & _grad_pressure; +}; diff --git a/include/kernels/ADStokesPressureByParts.h b/include/kernels/ADStokesPressureByParts.h new file mode 100644 index 0000000..5c589cf --- /dev/null +++ b/include/kernels/ADStokesPressureByParts.h @@ -0,0 +1,25 @@ +//* This file is part of the MOOSE framework +//* https://www.mooseframework.org +//* +//* All rights reserved, see COPYRIGHT for full restrictions +//* https://github.com/idaholab/moose/blob/master/COPYRIGHT +//* +//* Licensed under LGPL 2.1, please see LICENSE for details +//* https://www.gnu.org/licenses/lgpl-2.1.html + +#pragma once + +#include "ADVectorKernel.h" + +class ADStokesPressureByParts : public ADVectorKernel +{ +public: + static InputParameters validParams(); + + ADStokesPressureByParts(const InputParameters & parameters); + +protected: + virtual ADReal computeQpResidual() override; + + const ADVariableValue & _pressure; +}; diff --git a/include/kernels/ADStokesStressDivergence.h b/include/kernels/ADStokesStressDivergence.h index 6974b35..66f359c 100644 --- a/include/kernels/ADStokesStressDivergence.h +++ b/include/kernels/ADStokesStressDivergence.h @@ -21,6 +21,5 @@ class ADStokesStressDivergence : public ADVectorKernel protected: virtual ADReal computeQpResidual() override; - const ADMaterialProperty & _stress; - const ADVariableGradient & _pressure; + const ADMaterialProperty & _deviatoric_stress; }; diff --git a/src/kernels/ADStokesPressure.C b/src/kernels/ADStokesPressure.C new file mode 100644 index 0000000..b229f55 --- /dev/null +++ b/src/kernels/ADStokesPressure.C @@ -0,0 +1,35 @@ +//* This file is part of the MOOSE framework +//* https://www.mooseframework.org +//* +//* All rights reserved, see COPYRIGHT for full restrictions +//* https://github.com/idaholab/moose/blob/master/COPYRIGHT +//* +//* Licensed under LGPL 2.1, please see LICENSE for details +//* https://www.gnu.org/licenses/lgpl-2.1.html + +#include "ADStokesPressure.h" + +registerMooseObject("DeerApp", ADStokesPressure); + +InputParameters +ADStokesPressure::validParams() +{ + InputParameters params = ADVectorKernel::validParams(); + params.addClassDescription("The pressure part of the velocity kernel for Stokes flow, do not integrate by parts"); + + params.addRequiredCoupledVar("pressure", "The pressure"); + + return params; +} + +ADStokesPressure::ADStokesPressure(const InputParameters & parameters) + : ADVectorKernel(parameters), + _grad_pressure(adCoupledGradient("pressure")) +{ +} + +ADReal +ADStokesPressure::computeQpResidual() +{ + return _test[_i][_qp] * _grad_pressure[_qp]; +} diff --git a/src/kernels/ADStokesPressureByParts.C b/src/kernels/ADStokesPressureByParts.C new file mode 100644 index 0000000..845a0c6 --- /dev/null +++ b/src/kernels/ADStokesPressureByParts.C @@ -0,0 +1,35 @@ +//* This file is part of the MOOSE framework +//* https://www.mooseframework.org +//* +//* All rights reserved, see COPYRIGHT for full restrictions +//* https://github.com/idaholab/moose/blob/master/COPYRIGHT +//* +//* Licensed under LGPL 2.1, please see LICENSE for details +//* https://www.gnu.org/licenses/lgpl-2.1.html + +#include "ADStokesPressureByParts.h" + +registerMooseObject("DeerApp", ADStokesPressureByParts); + +InputParameters +ADStokesPressureByParts::validParams() +{ + InputParameters params = ADVectorKernel::validParams(); + params.addClassDescription("The pressure part of the velocity kernel for Stokes flow, do not integrate by parts"); + + params.addRequiredCoupledVar("pressure", "The pressure"); + + return params; +} + +ADStokesPressureByParts::ADStokesPressureByParts(const InputParameters & parameters) + : ADVectorKernel(parameters), + _pressure(adCoupledValue("pressure")) +{ +} + +ADReal +ADStokesPressureByParts::computeQpResidual() +{ + return -_grad_test[_i][_qp].tr() * _pressure[_qp]; +} diff --git a/src/kernels/ADStokesStressDivergence.C b/src/kernels/ADStokesStressDivergence.C index 258d884..4ff61e4 100644 --- a/src/kernels/ADStokesStressDivergence.C +++ b/src/kernels/ADStokesStressDivergence.C @@ -17,20 +17,17 @@ ADStokesStressDivergence::validParams() InputParameters params = ADVectorKernel::validParams(); params.addClassDescription("The unstabilized stress diveregence kernel for Stokes flow"); - params.addRequiredCoupledVar("pressure", "The pressure"); - return params; } ADStokesStressDivergence::ADStokesStressDivergence(const InputParameters & parameters) : ADVectorKernel(parameters), - _stress(getADMaterialPropertyByName("stress")), - _pressure(adCoupledGradient("pressure")) + _deviatoric_stress(getADMaterialPropertyByName("deviatoric_stress")) { } ADReal ADStokesStressDivergence::computeQpResidual() { - return _stress[_qp].contract(_grad_test[_i][_qp]) + _test[_i][_qp] * _pressure[_qp]; + return _deviatoric_stress[_qp].contract(_grad_test[_i][_qp]); } diff --git a/src/materials/StokesLinearViscous.C b/src/materials/StokesLinearViscous.C index 140c27d..3ae9a0a 100644 --- a/src/materials/StokesLinearViscous.C +++ b/src/materials/StokesLinearViscous.C @@ -24,7 +24,7 @@ StokesLinearViscous::validParams() StokesLinearViscous::StokesLinearViscous(const InputParameters & parameters) : DerivativeMaterialInterface(parameters), - _stress(declareADProperty("stress")), + _stress(declareADProperty("deviatoric_stress")), _strain_rate(getADMaterialPropertyByName("strain_rate")), _mu(getADMaterialProperty("mu")) { @@ -34,4 +34,4 @@ void StokesLinearViscous::computeQpProperties() { _stress[_qp] = 2.0 * _strain_rate[_qp] * _mu[_qp]; -} \ No newline at end of file +} diff --git a/src/materials/StokesPowerLawCreep.C b/src/materials/StokesPowerLawCreep.C index 3536e63..7093c49 100644 --- a/src/materials/StokesPowerLawCreep.C +++ b/src/materials/StokesPowerLawCreep.C @@ -25,7 +25,7 @@ StokesPowerLawCreep::validParams() StokesPowerLawCreep::StokesPowerLawCreep(const InputParameters & parameters) : DerivativeMaterialInterface(parameters), - _stress(declareADProperty("stress")), + _stress(declareADProperty("deviatoric_stress")), _strain_rate(getADMaterialPropertyByName("strain_rate")), _A(getADMaterialProperty("A")), _n(getADMaterialProperty("n")) @@ -41,6 +41,6 @@ StokesPowerLawCreep::computeQpProperties() auto effective_strain_rate = std::sqrt(2.0 / 3.0 * dev_strain_rate.contract(dev_strain_rate)); auto effective_stress = std::pow(effective_strain_rate / A, 1.0 / n); - + _stress[_qp] = 2.0 / 3.0 * dev_strain_rate / (A * std::pow(effective_stress, n - 1.0)); -} \ No newline at end of file +} diff --git a/src/materials/StokesStrainRate.C b/src/materials/StokesStrainRate.C index d095a32..1984ec3 100644 --- a/src/materials/StokesStrainRate.C +++ b/src/materials/StokesStrainRate.C @@ -32,4 +32,5 @@ void StokesStrainRate::computeQpProperties() { _strain_rate[_qp] = 0.5 * (_grad_vel[_qp] + _grad_vel[_qp].transpose()); -} \ No newline at end of file + _strain_rate[_qp] = _strain_rate[_qp].deviatoric(); +} diff --git a/tests/stokes/flow.i b/tests/stokes/flow.i deleted file mode 100644 index cf4a2f5..0000000 --- a/tests/stokes/flow.i +++ /dev/null @@ -1,155 +0,0 @@ -[Mesh] - [mesh] - type = GeneratedMeshGenerator - dim = 3 - xmin = 0 - xmax = 1 - ymin = 0 - ymax = 1 - zmin = 0 - zmax = 1 - nx = 5 - ny = 5 - nz = 5 - elem_type = HEX20 - [] -[] - -[Variables] - [p] - order = FIRST - family = LAGRANGE - [] - [u] - order = SECOND - family = LAGRANGE_VEC - [] -[] - -[AuxVariables] - [stress_xx] - family = MONOMIAL - order = CONSTANT - [] - [stress_yy] - family = MONOMIAL - order = CONSTANT - [] - [stress_zz] - family = MONOMIAL - order = CONSTANT - [] -[] - -[AuxKernels] - [stress_xx] - type = ADMaterialRankTwoTensorAux - variable = stress_xx - property = 'stress' - i = 0 - j = 0 - [] - [stress_yy] - type = ADMaterialRankTwoTensorAux - variable = stress_yy - property = 'stress' - i = 1 - j = 1 - [] - [stress_zz] - type = ADMaterialRankTwoTensorAux - variable = stress_zz - property = 'stress' - i = 2 - j = 2 - [] -[] - -[ICs] - [u] - type = VectorConstantIC - variable = u - x_value = 1e-15 - y_value = 1e-15 - z_value = 1e-15 - [] -[] - -[Kernels] - [equil] - type = ADStokesStressDivergence - variable = u - pressure = p - [] - [incompressible] - type = ADStokesIncompressibility - variable = p - velocity = u - [] -[] - -[Materials] - [strain] - type = StokesStrainRate - velocity = u - [] - [stress] - type = StokesLinearViscous - mu = 1.0 - [] -[] - -[BCs] - [no_slip] - type = ADVectorFunctionDirichletBC - variable = u - boundary = 'bottom top' - [] - [inflow] - type = ADVectorFunctionDirichletBC - variable = u - boundary = 'right' - function_x = '-sin(pi * y)' - [] -[] - -[Preconditioning] - [FSP] - type = FSP - topsplit = 'up' - [up] - splitting = 'u p' - splitting_type = schur - petsc_options_iname = '-pc_fieldsplit_schur_fact_type -pc_fieldsplit_schur_precondition -ksp_gmres_restart -ksp_type -ksp_pc_side' - petsc_options_value = 'full self 300 fgmres right' - [] - [u] - vars = 'u' - petsc_options_iname = '-pc_type' - petsc_options_value = 'lu' - [] - [p] - vars = 'p' - petsc_options = '-pc_lsc_scale_diag' - petsc_options_iname = '-ksp_type -ksp_gmres_restart -pc_type -ksp_pc_side -lsc_pc_type -lsc_ksp_type -lsc_ksp_pc_side -lsc_ksp_gmres_restart' - petsc_options_value = 'fgmres 300 lsc right lu gmres right 300' - [] - [] -[] - -[Executioner] - type = Steady - - solve_type = 'newton' - line_search = none - - l_max_its = 10 - l_tol = 1e-14 - nl_max_its = 15 - nl_rel_tol = 1e-8 - nl_abs_tol = 1e-6 -[] - -[Outputs] - exodus = true -[] diff --git a/tests/stokes/gold/flow_out.e b/tests/stokes/gold/flow_out.e deleted file mode 100644 index 9a084442405d5a963a8722d772b742a7e63c9dbb..0000000000000000000000000000000000000000 GIT binary patch literal 0 HcmV?d00001 literal 157288 zcmeF42YeL8`^FPMnqWbiih6?5l#YsEIJ%%TX=1_Sl3c<;l8ecOW>@UJgB2AE*s-8k z4p9-T2q>U*MC@3wi~gTycHS#Frk+;y7k&%gLMy;T~(B~-)`idE|J~fKU3O%|0LdL3sT8tM3a($(0dzI}7 zYBJ4N;Nwo`f!bEwlU*=D?L_gqlFw7>3tD+Is!MvQKNrW1N#93n_=_mhLFrI?U#;UQ z_Lbn6YD0ab*I~M)Q_92xjkJyl(mEwb>zp9%m;`CZ>a@M3JVolT8Do(imBz!GX*{f% z#>1LvJgk|1xaaQ!mI>GJ}A5qzF;3iE&(U~g}+ z-!suyQtA&BA)2ul$n2T9@Op$bY|O-_Qg8eMe{1=&BL`$EpiSh% zG7x;~+=V|iv43$%pqOn>_N8t4y_(g>RQak7=z|n1ujHSJc&(z+xZ|-*^5MRXtw^0& z2y0^h=;Lu;r{hWfm$7fH()i=C9mps8cq@=t6Z=OWkNYNjyl=2?^_&1&$7Cf>o3L%l z2W4+G(~8vPajku)sFZT{J%VdAWWlcb)wW-RsU zaIE)mkSp|ZS}^@l4P9<2#WON;R9$icMMb`xAo~H2tNK4eId&~IJ_V-~TP9ZfRPWz+ zbjSAC2_NXS*w6YJcotzDwSVmKu}}De@0r`J*llkhtf}@dRf93h%DGNH%$I9LroVu& zX4HI9&b~96a9^#<0iFADT#-+d?u-kOD7sGje#U*ZE++-*(oXV;(!HZU5=GZ(-^aMG z*5$-P-MTW8F!D1+E@IxzqL*MjoSAf&JUt1zNnXY-ek@^!xJ}sZot0TQ135&(_#FhIr({r$F4ik zTYyo^(qm`Ghnn=A3?NL4jP02>evf!eJ{lSyhkUbdz&5!DeELIkegS%4VY3!Iy2=t)|vy0J~MC-(b=2i}5!$ zZ^vZwZXX?&w^Ootf9rv~*7yq(c!J(+#A!Run_>riRp&u^+3n@myCVZ?u-VKeA+GxR=b}=MD)~V$>f= zmrypjUyRKsaoH}RXWxE--OPBq=*X6Sdzj>FR_uR7GjN7HleoL(KTStks6Ekl5JVW+-;IOBZ{Zq1e+dUF>9rVkf(FU2oG~W<1U=Hu>WAGwmbswEbq_ zbBCjy#a?EZwqM$Q7&qhEkI^rkwo8AN@lfL7d}g~d<7PXQxH-O=@eS4y;^r-Ld^6+b zI4E&*d^6+bI4E($USfxE{q^>w?c!8_skd1#i96LNU2fIy+}?5R*l6!g z+oj!0zmRx1pE(|wap@0|E^+aXT_}E&?Gl%MVi!vPk?j&U#{)Ai{Yug$ZrDrg5U#)8 zp0r(@>hDxfr~0JJogN=}ot$JliEGD3dw1F{?OytY#KZZd-`j=KA7s14rQh3y(obZ& z#HHWcg<==kE^+DicA@ku*)DOzUSfxE{q^>w?c!8_r+PZoCtdFJ_>dkS?g+1E$v7aP z^dB=W{oXE=ekI!_uKT|^9+>g;_+VMY(Jx|OGfdn2|7Bd0b}XUHe-cW%gpw{{IQ>Ny zKeE1W_J{C(x<93wKW=}sUmC|Zj~9+7yI)%0^5<#>@T*-Y=R309jN7FryPsWtNtZBH zz1@{h&Ts9)4OYvHTSqhZujqcFW9FB13GKGqrSGnMU(Vns_4(e+XLnxUw4ZD@^Gmvf zcH8aJcUL}ho;Bk#?@GEEw@XiUKfC;rE@7&AyD49W@sW8*LSwV-eQuc#C7;A4zogs7 zU)23zi51%JXKw$0&SyJstk(50$89rix1Xd-XqT^|+hu&S3pZFENO{~g^-P_e-7e#b zU1*L+W?b^wrSGck(x2@@b9^x4lFu%E7j2h*D53OY2_;=ZNtZC3{^AUN!q3N~#|Lc> zvtOI>)i(P{x`cN5D!Se5r)GSE<=OQ3p#3A`hJ>k(FVcVQLbKnSamien|Ihi-=MR>r!sCe3IJ2|crGMImX1_J# zlFu%ES8bR6YZsdR-i%8=yYyYOUHYMf((fgdbO|M0!f?9LPt)fQGx)i~ZO;27eLe0J%( zXuH)4?ElI*U>Dlww6@4NVH+;d{lg|S_04ja+a(@%ytKZUMuq=2`-op-3W-R#oQD5o zzyC@;v%GfY4)6D(?iUgo^*8gI@i-kLe|Wo2m$-!gHoJ$*Wm6B?E@8O*Y5R#uR4DC4 z!kCRwk+l7E|6A=4%6Xez=#<~wZt9rZ?c)Ed`J`S>_1;z6MOVrpp`_b|PW$cZ?NTl| zj)anK7dq{?tG7$J%v&*FkpS25V4t+l6-d?9xpgyZy}c|J8g_FSEW9 z-__gA<2W79++X5m=%ll=+vRwoBcY_*g--kJ>g`f4IgW&qZWlW3x2v{)IfEbP=SyV0 z!VI1Ao7*dz|9>@~)XTJ!#CP>}^Egh&GxwLc89M3g>~=Yx=twB(cA?XLyL!8nOO7L< zq}zo~`|YajvaV+r%6gz(XqV3}-PEz$&rJVc%_sFT>nrhHz1=*H)A7vxC2odJIy<{v zjwd=2O1fR>wBN4YF6ENrNGR!cq0@f5XnV#U{>k>{SZ^HNL+y|q2nGY{Rjqrb0ZRP& zBfucRx0O|2RJ^t^>?LR&Kr^o+2x&6y|KsroY zU=Otb-+gPo^Oi50=R4*3e&Jex??mVO%lXb;z89SD@aFrr_XB+AHQ!~;cU0F0e5W+u z6V3Pi$~%1dze5fNd@uST;81WFaDl@?Q_u_?0geRCK?~3lv;wWcQQ&CM2Jrtgv;*w{ zUy01uLh}{RodMs2&G%a$2fBdcL08ZX@V((#pgVAb9-t@a1$u)%pfBhLP6YkI0B{mG z84Lu2z+f;03}vOx~W1wN1m@&VtC z><8lk*C`4Buhxq|0Py1NY)}ID-eeBtWndzh1SW$iU@ABVoD0qa=YtEth2SD^F}MUw z1DAr!z~$fya3#13Tn(-P*MjT7_233@Be)6N3~m9pg4@9DU^#TwHh_&_6W9zs16#o7U@Q0nd}d0fgAJyJwY$f8}tEvK|gRJ=nn>flfcPf zAQ%J&gCSrjI0Xy?!@&qJ5{v?;g3;hKa5^{xoC(H&vw#PT1zwO1azHNdfjp28#sNPV z4<>*DPzZ`Z02G6>K?x`YK~M%Jf=OU9m;$DPbHKUaJa9g^09*(z0vCfzz%+0vxC~ql zt^ikptH9OZ8gMPR4qOjz05^i0z|G(ma4WbC+zzIL8Q>0ZC%6mT4ekN=f|=kxa6fne z%mNRBhrq+&5ilFf0gr;mz~kTv@FaK&JPqc8XTY=IIq*D~2VMX#g8ASj@G^J>yb3~K z0eB581h0cPz?)zZcnd5BOTgRU9q=xA54;aP087D#U>R5rR)CMd$6zH`1y+Mkz^7mh zSPRyH^Yuc}piLeT&3aWv}?XaWuahl0a^3mgubf@a_da3p9BT7Z_I6=)5P0!M>3pe<+z+Jg?D zBj^M=gJZz4;5g6)91psJZr}uv1-b(_=mC0yUZ6MV1Nwq~;6%_L3;-vAlfghR2n+^8 zz))}s7zT!e5nv=31x^K{!D--ha0WONi~(l>4;Tx)ARFX>T;KzFARmkaelQ+P00p2B z6oCLJ24{m3Pzr*e3`_)*z+^B5Oa9^3$K1UG@3!7boca2vQCOb0W-9pFxI7q}bT1MUSg!F}L<@Bo+v9t018 zhruIYHkbn*1&@Kp!4u#~@Dz9&%mvSYXTfvec`y&W0A2+1!Asy}@CtYpgunvu8dwNk z2XBBk!6NV$SPYhcx4}E$UGN@wAAA6of)BwmupF!aAAyg-O0Wv72A_aW!5Xj@tOM)8 z2Cxxq0-M2SU<>#hYz1F{FTq#fYw!*D7JLV`f$zZ&;79Ni*q(uo4XS{upc>cTs7b9xpfRB4f#4u;FrY5Ant($9H4g=c z0T-YywGIbO0X3=B3>*Qdc_e5KS^(-&t0ia!sM#7E1&#*PrB)lz7EqH~?Ld1#%?_X= z=me-st&+z)eW2gsF?-212>>9wR(V_fSSEPZ_o!&mp-7~ z`htFdKA_!B1pNX1FaVqcP6qS=?KTh$0`$XRFa!(*)S}%^0o0@x?KTWhlYSTuMgZzk zi*_3csObWuz^Q<`)S}%+18P!>b~_DFvjsREoB^myE!yo&Ku!8#3^)r=ms+%&2T+q* zwA)xfP5QwLvH^9eMZ4tyYIXs+zz3*HE&3u4P?K8pMLrk{=!bE@52#Bm`eHnwW*;yC z6aebd2eexuC<62W?G^yVfPOd|lz>t|AJA?=PzLCSiC_|#45&rBO##%T7VS0_P?LT* z2b>G2OD)>%JU~tQ;e2obpf0s&w+jI^sYSb81gJ?rTnsJ&)TI{fHVsgdez+7|2B=Fd z+U;_1A)wY+a0Q?y{ct6?3Q(6?wAiE!ypVKu!AL0Wb?tms+&jgMga!!$aU<{XJ z13-Q7CHM+_4ZZ>2g73gK@I9bbBhVO7lUhH39|1KF1&09_pf0t30^0#Ke+IvRUjcQg z)eIZ~s7bBg!0&*Xt-(>?Xh2?2)cgzV0DlAOQmY+k52#74f2`9QrFksk$AK<@ zy3oqZ1Xb`InweEWH9)&jms;Jx34of^f(_LA2BPMCt+Zh0%vu1NF>4>tjhnye1K3U0 zMdIdb`T%yz+!ySJ@9;zPnn~RJO&`E+YAvNMzPrFEa4HxK=mXd-vjKo6{SaPnVOyXN zXt&eAFhC!`7nw+>X54&DAJA?G0%+0?;q{t^h*OJp8w;pOEd-fO0FOmK|K>4CM!WN(&XbD;Yw$&2>?KTwzkVY-oO|3uC7WBjU-~s?$rA1$$ zeP>d$5oiP20_d9ax6-07+5xuTv%sa`G5}rWSK1ABQkn-N>;O6f=$iAl@+*DO2|&}F zuaysIx6Xi?=!4-kGUzJ5(ihi*E3l1PE^r)xra50LAJ7-@yV9f|jt5=AG1$Ktz-}^s z(+A2I=!0B;Gv{mN1KKSMbO-c<8}tA@wHAGGKe!!f^egP9*6g5Z&ezHZ^hIw#P5Pk^ z=nJ50&fiLlzUT*fVLScsICug)gfwc=ZfF-ulYZzA1_0=q^SAOVeQ^?ira50LE&AeQ zK+R*pKrjeE*POqVU+Ig%0Gj4}t$aXV3<1=nABKWc0Cdgyn?6u$p=je=voz;xs6}6#4yZ{#oB_@R&^708C`$GoClz3&euwd)Pggt=8OWx0}H&pcLE+ZUgKa(2AbFl@@(59o&d* z@I&-`t$aYc%>dLq2iyVf1kg3-Z{=6o?JfXKbG}wSpxy2U)TAHo0rvvvn)A2v0qr&u zP!oQLp0AY;Xt((l40J~*A0%n7Hpe1MpS_ArkcAEnp1@yyX;BoK-pjLfw6rd)x=!+)-H5-98 zpe>*-wP?4e05$1{r@>r6U1~K3?Ep2YMY}x%sCh8x06GHdQj2zb7EqIZcn&-ds7tLT zpc9}bwP?3_fSQMaW5BV1y40fGUI5gjA6^9W0d=Y60>=SrQfn4?2~d-MI39Eb)TI`E z@iL$${qPES6;PKxpxwHG699ccyM@34KtH%a56~0P2ejL3U?HF%UI%Z0HvzS1w_bpn z)S}%M0cz3@eL!D8U24&8Zvkr34~xMPKwWClZv6l?sYScJ4X8;!^ald~b*V+Wy#uI8 zKfDXx1JtD!?RFBNCbej{_W?EOhk;-apf0s&w+{d{>4&A@LqJ_>(HDaOHK|2kECbY} zABKWc0ClNFUn~dIq#ss*j{tS)1KMpE7!K$I+U;Yo63`E$z^Py~pbu!bRbVxsA3gz} zf;E6zwA*Qbn$)7*)&gqM4`+Zg0d=WGyR8G%q#xFU4S>4TqTR*-YEp}K+X$#hKX|}c zKwWClZkqr#>4(kWGeBKx(QaNqO={6@TL3lbha8X#s7o!{?Q=j)`e7^h0#KJ)^o0*l zlUnq}mw=k|Lp~S>s7o#S;wwN+`r&Kv4WOK_zh623`_*nq}K1?4?xW+U@ABVP?uVNg1-PYcYwdaKY+T_Iv1P=s0pnq z86Xqisd*u|2wV)P3#}^Hp-NSJhh~*(U@uS|aq3d*5-<%wQ)$&eSQAk5a&QH>5jsVa! z=WFEy+U-a{P5PlZXaS&W&fm(f^hHYmO>@3hKA2 z(+`h<$H5awgAbH1+9GTR=!f>81L%l-=mXmAN$?2L=vUZH-T#HAIbSO++O0F7CjD>> zI2J(HoWGS8eQ_M;4`{c(pdX+gP6YkI0PF)FM9<$!3wEn=5`d;TUn?!z?PNer`e7g#1fXlq z-^#D_#b5wUbG}wSpxuT5YSIrw!6^W`=KQVvN?!~E&@|_3H2Q#c zI}?ln^ut-e1IA(>YAph50W_5sec=VrH0NujMZ0AKYSIrmAQwQ_oWGS{=?fo#ra50L zAJA@jfSU9}J{Sj}YtG-wuk?i@C(wYbpf~#Q1c>iF}MUk*POqV)--Ip6i~Ao*b7v*bd}a+h+htvYfAEM`LLpLvzpqGzHY6-DU%7(hqaMqky{9qAywkYEp~7cnnaJerOGj z0@S4zeepP;CjIaPcoI;zA!q=O25kU*pnUNZ!lwcK&>nOE9RYnnyUhj90Q%us@Emv^ zP^$^(1gJ^B(r)trH4g*FfMWr5sYScJ0H{enya?t4>Qak#I}T8jTD03sfSS#~@t`Z9 zF12X4mjN{&2Cslu0d=WGyLAK9q!#TK0@Q2)vOsr0U24&83jj6ghu6SDKwWCl7j8gJ zYS9<318UL_JwY!(U24%6Zvbl24{w4+fV%Vn?baLg0rUaw_7+$S=!X+Qe=q>h2ejJ~ z@HU_y-U07|_W-qMx03)hsYScJ52#5$3&7^u;JZO={5>p8#so52L|pfV$Lr1AGdoNk6OsYXNoX1KRC$a0Z|cXt#A>J)j@X z0v<3H&M?dAp4q!#VA1yGZI$N{;4y40fGJ_ppKAGU%o0ClNF zyZHb$sYSbe38+aw*wAHD|P0P0eUcJl*jQj2!`7EqIZm;ed@b*V+W zeFvyXKWqcv1L{(Xz9YxUw32K4bU~f9>Cjs|T&ThI=)2OU61&^%0oy3sB-DA&@G5u>ECjEECEy+KE?5edf#qN&SPeb_>%a!E5o7||e=k5^ z)BbhAK45=9`!@myf+m1t*WsWUp#3?9wFb0*JJ11i2FHTqK{r7A_W-@X7VtUv5_|=| z1wVkFz|Y_}@CWz{{Db15e3?}N%9dFZ)CPVq9>5Rkp8nYg%K+xZ%qakKp}LoTA;L?* zr2u}&yb4?c;D_*92ac`oBjDJXzk@%)4n3U0E^1Ao213|Htp~s^Y7M>~LL67EY2esu z?E=THavrz9|U>rc5)wu%O z5p+~rf>y@OfbC))VB44vUIO&DHJ~Pq&)IFpX2;p;eFX?O4op#(8 zU<_9GGY&wAwxrh5TO)LX9sqVw=ftNX91Y-ewSE;qNIRnL;Wd=o@%=IIICuiQ4&DH7 z0OSHYgr9>ifY|>V)E~!F*JwDVT1)jH#PQU+C61@AXYvr@ zc3bZ{rQ8_Wc+f&~EWK>2eqLf9kY z3xKj^pj;U!PbS);dN%R^ghzrF0BtZ6{V}sU=n3dE`U`DPjXwnlVT(+ZBQprl7S;0@ zXoH!v!J`0eFmnM|2;Kx}Bg$WBgKGTYv0;l$*n#&y)%sFnga?B*pdIK4`T*ELtzn#j z5N)K&Sa3GLvDG^4WQ6B}JHTDwUhp({24HMZ*Wzf$YMo;}!j0fFK#dYm3aG^~Y$`Yh zQ13i&KDY>63@!oFz@^{{a3#13U_4OkYu6*Z0o({+6V*TAoAA2k418x9;Ul#k#&XUC z_k#z(EbtI`7(4>z0N7IXYuHn*0X>C~b$b>(2Ux#(;6*SWyaZkW@V{Eqcnu+ZpvD9E zL9JW7g%Ir_X8q@V#Fv5(!7{J{d;~rQtH5gT34s69we~uMwB;tiwoiL*0b2pum|6q) z8lklRZTJqKtM>mB!tKRqfP3<%%!6jT2Y<>u=x%%Pr_6()W)J?921|UUrJl)qYA?|B zpR%X+0$u;9dvGtvME6%TOYqsLp`n?pcZts%^yRv;r?>_U88W=DtJ%&8R(c7Wp^SmT zLT^#7YoNc#=jvIKUsmWV3YNNXs0x(T)ox^|ucWkHAx4!VUrD=SyyMxM?<;MW>qCrx z+60P&t~P#GUO`||TR)EC8ifr0K#^-eQC=YKu>DHBg}zCFk_l>aNuf80UAtDGWRk`5 zgFa?gSLBoL54v(-3k1h?@pb8t-7z=U>q5CYw(r=vP5Vx5I(8O~m~fy!dyubmoKzzs z8M9Xkk)eG@49|&|bGbUT?byC;rxcX*uEPiR9cnMwA$7stWhEtOB(4$uLSLl5Fy*Oe z>w~-{u8zmKI&|oIO#7}KIw|ux9is1KUrt%jn_b{jnx#Q+VX+d973_R`*G|VKCFpd( z|Nrf5`2&iI%Yv?cepI_+-EEgv;|yHoRqDs(@{|?%&o1+MygBT~o+9+^NC$MevRplL zbBFl~R3|ST=Pw@QEyC!5tHZ<}7dmWFuq06633~HOqm}gcEZ2a+qk0Y;&?j;JxKx+R zU&OLaj44)>aF(lOQNWYy%k!2M1Y5=Bk4vDSKj`-scyfbLGLf&Tmk{Houf&rbD8rbY z;FEd@y#V?aB;QiEp|8NlS=>V4T^%uao@fne#7 z?D4*wpk*VLIqH-8c<6GK+xg&>K7(+XP=aYg^Z z%-4$8P|YXak|{Al5~+w_&NyGr1Z$dttOqSASmMv&#SlIfv7xU>ogYe>V#5jvs*_HP zvPH34DwgGHF>t`(zCDMzMhqC#cle0DLp^;5_i2&HA1aNR@^I=9C`lm3q>7Z#Z= zvs~?MOOt2=GmWoo!9ZSKsW0)F4Xd1b35!bt#h!d$5Z%F3T9%jRpBz=>iu7d)`bvU! z4ci)i6WJ?K8p~4Zv)XV}xguXHQsuHzpC{K}T3q1G@#T7O3P|)L(y32c{$Q#_^l!ny~Dhm>IFEiNB<@B@(IhTK6%g9Q4M25 znkl(?S?6ou5K_=_UPd|Wbd2(iYDhuHc{yHBbmTY^Z+LFdFs?L70@?yX~W8#mq)4V zW!y{?dl|NB1-P<5IBPHBRFkJGTEO&{cv?*UJg=sdkvi4iv2T z^>S=^l|3n#XgPHv~5!@S%v!|<=@CJZU)I4@h}uxTfT zJt&yU_43wH4!R=fe7$Vc+s=*@OlDqAXs<-=lY&Xi%XWHZIU6N8FWc!QuqR8H+`No= zB0>I``I*O!o|mN@PWj5+j-~r@wj_+sDyt;-xfR_zW5{c*vjuJqIQ?&QdUyQxV5@X$e}lTT%5~e2ZP%jjUrE z?OFU`_+1G)o!dmWU(u1!v?m42y;Ydy>OExe@DV)+kI0~|nO}6nAq7iwyaoQEe147p zM5-KJPJ3M{IGlYA3Is}h(mi{5OMUVEKT4m1at|E;@wm{FhsPZ}K|DE98uS%oW0q@R z&k)lpU-yO36JYe`Axf-amgP&jdG>k+|TT1 zlF!^v;*DstN}5KoY=VAG!6c{ADEph`mhvH_V6xMwF|H);KNL)58jT(&qRJKN&z4|P z(`Y&MjOg8?1f8bQsQu&XX&i@w$xfr?>Pr17nB+8Cu3osD6ij&y5=t?L$Mm%~N<-)k_ykX=E*T+n!4J&CH#VQBv zM!_Vf(dhcy)K|~T$xfs3?M8H??9URWGL0r2C!+0W2_`j-hW!KNA3MM&?Ce|~{Nx_t-(WvoJ^dwaG_2Or}L|3CR{tf4k8#h?O zc)^6zXjHkx4~G5N4k+k2jY@tqOjJ$^COM6U+neZ$o`fo$f{xQ@tbb$j$Jv2`j?-w2 zKcsv%p%QeOMw1(7ta?z;aT<-QS7h9g`iDabCNYhg^^^Q&y0wm-)7SgPwE*J#?PC{=1){#l>3Zcy|H{9WpAmcldj$N zWY?MP^2_lh#Jr40(EUaHF&@|1o$@(vXR~oWbLiDO^7wgn4y>fhkf44PO!CYj(%+RD!Y2O(2*se_Z{UF;6$Jut@ZmNC_SNq1sin zpeOplZV|;(_ial3MJV50SVy^byWJBBfWQ>9o@AqA5gXRL99 zS^b5@1^yiTt~k+$OZRd9HOndO2O%|V8jK|m8)A_==qq&bw6;aeeF4#l)@QC_Jh)WO zWBeWL{E)C6QP8d%a~tlb#QcaduH8khd#rLYkKl=1`<hw=Gi=LM9Tg1F2n z?{@)ZaqIcz=_vU*YJwML1mgn*cjD|T>Uu)^nSx24eMvdZFt*$jbU6EJksS!)H`uOVptyzed_}=j&c4Fq zgjr6BTY}b4njFr);_a2F|54EK>?_=_9Lr~&A0~PBm8d-N?M9U_xwEh2%1OZ_&%Wa7 zYmT$3{X`4eoPEjtZ)tDkhcSY7XJ0Y?6+LMu64zmjpxxP5^f;5So>BgypySzB)c#T~ z3FGuBn9|u-;_*xLt#;^m_GP;+7tskHhb2tq?8_Wa~WoX0PeSI)j-+JSjoiAyN;HRDXja-Z$`RHFT({TpXLwpZSl zN~PUK`6qF{GpyGb;h{rc$(WeyuJY;?T~66ghYlxmbc*$cv8V%E3u{pc(d$7t5h)7> zJQ#3rQfECM8T;fQFMOivX_O~!p27NIxo3yb593AGGqL^Ht{lz|<1QY<<6cZVi}E)F z%R4)yV2AK|uawIS#SYYUI6I`^@Sx9|J0f@QT;7zMaWYzYvVD1h5+B<iO5*PzcmGP@U=N6ne{KL z++o4uL~AT+PLbDhIgG>Fj$y$PA71$tuBcJM_FWv;2gDA%?7>(V^b{9c>*5h9>ziIb z6hGL7;%BVeV9}YcN{H#j;kH;5C@J))UwXy%XDN^9M8b%m?!|Um3HynZH=>7mx!mih z*zr*E*_V&wVf28NpzWZisAiL3UP-8H2o~fAO8mibg%Q_exhz}6kB3rT2`&9tL5v1D zJ`X(+er}gujk6vEyJke_VgWG@1{a@%E9m=2z;0A@T00U_b`&1mlg=&kRv+m#e^&ps_@I*=jb) z+L3|@ToEHE*SBVQ%($`x1rrT2VMS%b4xOXUBgG$PXz4>RGW3U4qQ2k}beuazoi|9i zZ2S@<=+W<5jncy)Q|t?t=GcyhQcgJzFK#D$Cg=GIJ!QqYk(b7VHM3luJ09Qp_+yXj zc>FO@<+Ldue@9=KmzV3YUr)<&O)9qU$ML#&Y(jY%hN3aOC7xoO!Qz(@`MQj7O8fS9 z_ba5H5<*`+hrx!J>z|m=-SqNNnDwIDVSKqF{U|KR-?t=Ep_H@l49C1l74P zZa*f_V+q^ile8Vff(eUfRrwhEb|~eFgi4U(YI$A`kMrUOyAZlz!ScP_C8ixl$`jvy zAXrxD$uIHeCNiK^WheVddDU_7uwY&ux`{VGP~7<57ltz!7HY55uTvh5y5zSC2F6zGs^K(1X1OAJAARp6bhuCqOfUP{0cCbVC}#BWH_g1D{=pDjiEm|l*#1D;qVDwTq` z_34>d>KX41UovCH$J^1vPKnyz@v9cZy;Kk0{Z7$9w%=QX=?d@O8|N!gBetVT)DB88 z=q(9azwgY(&q8;;%3LyyUaLa6A|JNHR6ZZtXSaTNi$OL} z#tTS&W-Milj5f(BLU4k=*yEeb)gZnf+qzN|IX0VAtBjZMmlo#X?VGH5BB!MKYAHNK zh6l84&*!pSgZrL3V#wg+Dx8slD&yIhU@&$SGE(CT^b|hz5+z`MP2GlgM60M2qhX<1 zl`9S6l^Z$u_XYv|D(ckxa6Dc|IHuutVmfUmXXeMepY0F_4@kWa=Re%+^#)`C#w`3F4?d?0`{5<5MZp2R%a>=KB!8??zL0EHg2}zB zO0tw6Qn1_=#ERHZuPv4FUn>;*igJBLIg#fHQXOqr!c<<1RDON2m4ci<%YApBqe|2c z6qJFV-{U@6uBd($mC^1i1tuq6!M3#F6knST$)4F*Ri0AWf>SG3u##UR{2!R4{@wr1 z0F6^kQoggx`~{UXQL*(pPg3T8G_d}gfwopVRUYX%6iiQ2v_j;=7ab1|1}6AQhpT@F z*a!axrUWmC$?=KRVngP{azipdCyC8y7o%XM&w8mm&}dh%a5_tv`cU!p?Aq_KpKzpaSL;Six9PnM*&{3rjTLF#8;4$scvb`{qp9Be0> z<=Huj|NBBU(ff<|^BMLRKc0Blxf!3$W@jh*aC!Kzi&U8^({L zR7k<}xecul`S3A*y(R*44MQV~F@dXl0QA|IKg2C8SFNBK*00%fsQi!osabEYS$i2e92JxRs3kpwZCsj!THZ9X`2 zXq(EjX_QINda=f=$a4;r{8+*sc=3pTw_G{%oz-m|F3dRESU=hBBmePvXh{I8=Y{$z z4aJG=e-3kAMNWz3)=a`TvrsVA%ZhSe=u?prmU{z*>KORwmM)Qx{aTvd5nr(TMq$x>GDc^9KuI(FzFa$H;67nfi5#T-0y>JOll(O0aw+%Xv0q(MASDkcuSq*b3x-wVo*b`03Hb=R`rSgfcSL_l-OJ@Z63XGLAN~i{ptzyl&NaB9 zUgYNk=whqr((!xhRptPzp1e1uta$2z-rQWQ*ktDh3b9y*5h7*@gC$ILPEl`0K!x!m zwrpQX0xzdJ1$YaJ$9X*^fuMA@*s7;q!fao#9K8xux!eo0ir6q458~zW9TR#y3{3Lj zJs|kcYu0U=)V-X8|36j|zz>na-47kudW8Z9u+(iB8L`7dv^8Rf2XKfc?`fVaVdV1I zTCyS~q@Zn&&T@@D?R2}2n;;$Qn_e`V+!8!vhqDkLUKQlU^Wl-Un0g6w{Za2mhzQBI zRP4tBJ3N=hc9^_FEI|zLIO|GKS*whxS2^EXSXfTqtALP9Z!=itMmv`LnaU9W_%oa2Dwot(%uS+3KP7UPalf@(;N^LFa(FkHg$upcA3 z6M}fn3x2K^qXyTcE~tJ)X`NMiCXK@=8+)>o`tXI{c7`h(ZY+^6sg^J$4<@s-TDs?~ zURetntXmTImvO`sWde&VU$kRMjABzTy>w42L_S$AZ9HG@plC}PrVB@z6Ti_eP5A0o zZcM#lYr&`yi#p*?g#{D6x|OOG5{#@mc<}y-*w@sjeg=ywiy&7W`15kSzACE9sa~h> zcLC*`e`UEkg*}+e7nTIq%7r$QC~erk@92)n&r?d!>3<<5Q!?2{398=(u$;IcP```I z_vxQQS7gM-LW{Sov=pmz*q&yM$9JyM$9KyM$9Kw}k1tW->SYYr#khGV@5)}{67*ZC*6QPy#&=<<#aOuctw0EGeFi>P2*t&4kg-BQp{A8$@ zFNN?VdLboezSL*EC_%o)6h8%#8B~R-UeV*?cAPI~f_i0ljF3c9FJVrgAbM3;gd|M; z%z-QIfr135Qp(bRiNeqkj$CBLYwxF8VQfu(pAa>_V7I}jyceoeh$SNh{I z<4bX~&Hj^m!^ZirzYuQ=VoB7eph6~wys4z{yfzfu}EGN#Es@~ zOHaKW7W+!_u+ZQwj{eEo&Qv)fhzCz5#im9(T6yv_U_5Td;9%i^iy*4BlH|@QufmWc((Sq`(`>>NMFnn1hZ&V-oSwse5^HVbmPA^Me0@7l;jw z{l`z0b-ytANyuqWbw%HA_|Qr#8Qtn|D|5InSXLZu8DEKYQ5v0w2hiG~Tx{6sTbZ;$ z^al?#(lkkbn77h@m6=yoVBein9}iuw38lsA-^4`MH2fVu%rC;dLrnkHzftq$mS1~j zEcahwsDFe}zDCSe3KrqPzA5(ea_YA-TcRq*j5nxqo+^i{Q~UNV6xD@i_QnNW1!XQP*@|MLRX|HH@$;Z+8-Bmc9^zOs z*})RTyOZr_z0?L9OEB@nk-CJb{>F#d$~n8KSWhX<8CQn?b(QG7l(3y=4dM?=G0=Q@`*rt^>@8o@WjmDwAudR=r|QMPlE4pM?oCIVo!| zOEBSBEjFz5X*8t*W?02YPotI{#3|`%)S;DIhDuMP4rbcDo9Przqm_2SY}#K)Q(&~ z^P9eQHZPJbLzO!(w@|CdrM`kZb%m=hcGXp_sfJEz3GUu`xgzH_>3P}ml=b(&aCzIJ z%vGGA!P!g%R=vX4%ZV;;RX?)H9-nLprsw5&%Osz|NwjUdN?UZ!$oLQ3Umhj@0<{vW zoxH`9*t3%F1Xxo*IcGL}5(D=F#4gsI0Ev~#`hRr6+4agfv!P(>Q$U`#&|hFa7GrgC z3J⁡;^2BPfYSfz06ncBaIbVS;7N|K|fv+5t-(tcnB*}<#O&-P%zc!U~wz6BscMU z2v|_u0Ut1Ic(SiHu$qnj@4->^jQU!UDwlI2!V(S&;9t5*+v|mY&>Y?D;=Z#i;uQ4` ztvtN5DEzNwv60PY1fEHb%W0F4Sa6~@{7_Fts_ZQ(@!GX(lrWVWNS0I7sS=Bt=u?pr z@`c-hBK%A?(IRX`+M&;n7pWvF9n&gO<#O&-P%!nLl2_gn-W4gK!*dR)-zc=6CQjrL zWy915SoOj}{5SQ4eo4K`m3-!4(<3XMT*>R~J3EagSgWKX{=cSCl*|K^vLX9?5O zD6J6r;8tee-a`iWNw#G&XL=gNyqum!lWj86_0rR5q}%aZMb;1Zuxa$)e*vCKCrIy` zwm{kfX$zz+khVbD0%;4REs(ZA+5%|{q%DxPK-vOn3#2WOwm{kfX$zz+khZ`cZ-En* zRekcFt)m*I7&bEYZ=mD5Dzqz~T{&cbvpitMYh}N z*riK836s^8{UzMp`Q*4}Xx7u*F7af`D|!-2elsrV$?D4f67KGNa$GaCtCzW7vgH*$ z2}Q?@OM0@pvcH78JD(ia4DIS=?q?S_^GjSp$!EqT-7ZXZf6@83=Ckv+l-mrY9(JMZ z@1!gFWq%2GcRo3;U1*lmseGc7YCb1jr|oiFyD-`O>2h1=k;&T2&fdwUoAypN-|pHk zRXa$1>_REGldj~K{UzMp`Q*5Ep;=C+@`+BW`J8l}w##wt!esNO%bgw{qOY@)Z8y^H zr2W{1vY(T#$BavQvbwUrgu6SR9M=r(>SgX{7dP`u zTtdlb#wFb@Om%U1*oj zOqcv-x>NjrJ)cwkclGv)+ILr%%Sk_7?)3V6^tfiW2Z>84?a7Qwx*3}DtC>&Y5=wqE zF6m~ttM$#}ndOqWgtEUGmvjlua!b0IPvR2V?I-EU>dO8S?(TeYTr)K5X>ONzvgH*$ z2_?T7m-J+HWq%2GcRo3;8QRs$+%MVkik^g`W5y*tSzXy*!rh%uj%$W?^)mOfi<|i+ zE}`Tz*1%jr};(MdI*ldjWt zIj&upZ2ok))7RJ0=b_2A8)U@D-68X@bHD!TA*I#V+z@(i+2GzE-Q2_dch-j9RR=8&eK_HPd5v!z>R$Qd zq0NS$(>nCl+La9tdFfd9cYjsA^S2vkxj*Rc{o=lsABSeO^xk)O;aK;UFSpcqzIcZF z&mU@Bc}JIA_q}(Yv!K<0@7+tkxUFbHaDM3RRrMbG{o!uzFE{-&UV0+q~`{yR2)}|NJKIr_Xz5Va>amy4MG5Wp2J^x_fPd!3~Oj z_Pe*Gk zxo6j}?dbmbp{dy~5t=7*-cU-TYxz?!jb+}jQq`{+4eb#|}2YrsXfZhXo8 z=?gOsTwp($&YW?R`-jOR4y@j7t^1dmPdqqvp~wAus~TBF-B!DwyRh-!Kd)~PdUO6ufwPJ} zaz8ur?6UTvMvb@k$o%D_tY&|{)T8Q=!wzb&e7gJB zUuNH4Eq7Iqy&A3d9CFoY_dmZaZv9k^{oFgoKJi1nZcQ!wtyyN-uW*cIzfWG%_WRb_ zzIMB1za3v_J0Ga+_qMj*T~@x;&9(ie==>A4{f1liyXSHp?`+wx>C4)FqxJS>mi-!Z z(e^vYvfr2yI=-K_UmweUt-N}B6K%hNmi_j6K-({Grpe9ygBA<3%fLoLDng zx4*V}eCVUch3>ljRnz0a0zD3F)BXQAtNpcEVDHWz>wdo2>feXFX7%d^ zC+PO~p4Fe5)zkg>e5?Ic|3UZPu6kTQM33_>YaG8}rXJTX(&P3rJ&r%D$L(3xczyC0 zdi)-0jn9wt*W>kRdVJ2(hZO<9#4O>#?hKf^td@dkE3<1ar2Q4dYsI##?52iwZ_Swr&!}- z;0iq+9%zk^b3fPP;Z=J4%eKbD6I#0e8C9q2aM$?If^+?a{qNr9Ue)mlSIy)5h2Fn% z)vfm*@QVBICBJlk^tb~s0r*{Vu38r_N*U z&=)^0duzp=?x#P!WKPjp@43I8_F;>A*8b@J`qjexwRu(D8`}3f@7}zb?kz8O*`9w} zOZU3@-(E3l;|TZ4dxzJ#y8GqsEx(@jiGrCHF1*qP}~ZtMruTmQ8kV{c6Ddj}Lv`y`y2> zue+Sr!~NUFx??AN*vS2B?!BX)dA-#A>+N@DG(K=&_b(mW6kKsyEB7CNUi;v(7Du_a zpEhID5%aqD$hdXCAAX^J(9bpyJ(Uv9bi>_?B87+Ssbtu9?Iy4C%HyV0o8mt7cI*`w2y zhgN^v{lN*J|5aN3{?Nij*PZd+vt2?fu3G)tv^VFwH%xxw&FpOly4Rd@`Hd%bf7rcY z#MJvPYQ4t2e8)STrjK0g-gd=4XKrbV`Rw&mYI#C8xwqZc=BA--FLi&Ban_AT)O^hS z=ZROY9`*P__wRM@nK|RO{_a2OpM2LwcY*uoOBY}A*}Z4?s8xFHj?1#U_t@*X4tI}i z`&W4pIyNNnJL6`d}-986f?X&*f zx;u=NjD?|66Cv<;WsQb^?R{!zM zxBAVk8+1QutNV>x_mexU{*gUY_lK)>|7dOXhjTkv{b1pFR=aPG`Jeq@bV#@REmnIE z%+~FDk#6sQ>+)V}mH(ar+Wx0$`}ft`&(+(T=bI_%uHUEH zKhNs;HXWa%>p4KzyTA6=SY7W$+HX(l_(UClRo7qaF<<*}uC~WkJw9}?;^)=Z@fLb~ zI8u)f`)U97x9s#n8|~*8wcTFS@nx31ewnH5_mlQ}Yu(Ql>3+6A_p3v6e`>AUYi-@n z_SNn6Y~8PVg!`8@KK!Bk*^Kb`V2u;qb$gwp``H5BuW)^$#)&6%dmX3y*_pauy{W^C zb$fN|ezruH^GaRLi>-L4$8`Kq9e+fZd#p}BT$i_j4qJyqYg~Ism%E(~o9ggbUGDqC z<+aLvua4iZ%bl&uJ4=W6=i>5)jEEk@KS zeE3+mw}HA{-LAtcwY{@-ySiPspXxfy)b<{&+fU2z_+ZU14{Lij*X^W%?mv@tfAZ;g zM;$NG{b#i9Pq*pu;Va!=&eP#tx}Vgx`c2k)-7mW6(67gbe5-$K{Y3YN+B!T`j}JIM zQ{%(Ux}CqJ!;f|Qo~Or$clCMJ{rdWOKYjlCp+0{N>hZ0yK7Z(}$Lj`qygpN(Z=I#D zkE`kP*S7lnYlJ@U+N954r|9#C=6W1I+dAL+b(X&FZL7~;n_K5yb1&D|xxMuGxR)N^ zf49!J4trLgfAzM`UmwfR*RS8|^Vh@l`NMcU?jNnMQ?J+OUzh3g*A3QrSD{;9m#)y~ z5AW*vK(0S;)8}8c_4(@``g*C6zTR1*=iwqf55KR^>rd6!MUPqM^~db3uXldY^YDCq z{@z7j|D2<*YdYw8_yv8Q-c?`cjMCRH8}&S#Wu1=~XXxvhSM+t@G(8U=q0hU&bF1sV zuG_70uliJ%`#xRnxmLMr)YI21;Rg(RWbh-Vy-0$l6U|sG! zUGC|++)Z`-WnJ!vbh&@i$Fq#^=_`d-fgC@ zzs|C*$4=^^uWJ|T>)Lnqb=yb?rm?I&H4LZo9_1?)7Za*JpoN*SWLr)9YA4eO-I0 zzD{ecuiM7z>)!A5_1Qmq9c#9}9vh;sYhTstSl#t?TWjmOx9Prm{p3EqZqh`ro2=36 zCSU3Gsm^+R>UO<8FjKFSG}r4R)7;j&)46)x1AEt&)lNdM>^_tkjwRY z$8*+o{L#639i+%w-#D<5Ue_3`uh*Zq);oNedR?u*E_Vm3+-vITa?iK6uM6sRp@Z~( z&2_omy4)o?-chfYZMMpN`#N3j(K?Rx53Y+n{kK)_Y8`aBC+hebUGCAk+!?yu=jwPJ zy-u@Dx6k8r`y8&@;~RQC_e^c?PxX3C7j5r_y8WH0*Hcc@_Ab=z&8OE>@^$;VMX#sS z)9Wd>TJ7oK3A!Ckv+Uh!n$>>(Sfty{MY_EV)ApXJ+evNRzYfs-t*O$fN9_Wp=&-`#b4?`Ymi2z|Uff7;9!hlW1Ns(X6Q=^w>yU-9V+9a|J_3oT!E!{+Os zER*e_Mbkfg$v3gEA1{0J;r>^K)7OlA{<}Nt z2g31Be}26|orfd(??2f0#8Yav2z|05bMHGgUl-Z$)W?4>+PE+re?PeLiS6w(LaV!9 zy{_>q^~34ww;x&X+qJKT;~TDOc*MXbYlY(DT6RhJ2E zKi)eWU)G>r_hXO!F0}5c4R!Zj@kTg((+AC}m7RW;ZG8EnZO66ya!6?7v9JAg_1bEV z`)#WC?$9~4)`d2WzwO+{eU^nb&iv`|#e;tfZMbU5RnyPf5L$QR{6YI>j}5IKI=X*0 z)+si%yuiKh()UAa=k@E;<@P$EjR(y+|F?c8g*G(l_Iz2jdZ9JLhOQcNeO_q8VOxKG zZuCi^PbR#xU+1qwq0R5_XfgPpQ$w3DL)Xn+656mMYrp*(b_#9!_J`a5sCRB?(~JdG ztIb{&T7Sm3PrkXaYG`A}*8*o(A0OKCbFJa;}<^t!~vnNJl%fmKIEd% zmJhnV^igeZXltwcZrZlIZ|JL2A8Xp^=eD8GFK;ws*1V%aTfg6M?`?Np9{Mh@&jEc# z{1Xnp+&t;d+qXXw`X*z<*!@q}M)V^tCUq(Dpk;+iyYWvsaF}rcJ$* zwEd=N`+XSNeBeP-`^>&u+i!xl-^<#5b3>o+`0cnI?N@30Woi4>4}Jdh6;JodEYtRz ztL?X!w%?7RuRe2KG4s!2ZNC$={aS>+e(IRpvj(1{?RT+dpOG2iP~RWa_Uj&(ZrSf5 zJs$Pf0r&l=f%!=e^_7t^T`G`ng@G$L+mjyB@cP zOI$)dKK~VtTjTR2iAxx=#^DtbmoQ|FyQ?H_h9PTwJtRJz<7vAu&2))})VO)d?-G|V zWQ~)bOWX`Y)_8bIvgxv4XvMQRmk({#F0|Zt@3s@R#g})*oDWO(KdEYH`K3o?_vlk4 zobQdM!%v!U&f?JWp7+;#zQLKHw|<{pNo=m!-)&R z=^L}pxcB)H10(TTx30Y=`{ripAAoEILiWn7p3^zDs@zB2rzHM-m%>T>7nazCca{Ta$Jf5gpK_t)jF z8(O!a+TgnuU7^c;m@fCrq0N2XJt%+mDMqC$d-K@OBzeI#Mm0%9NDv#F{3PzZR{jlO`?RT zRI=}-&_a7Ey*e#Y2_;3@%D!jGAnWg)=bY#6^L_igKlA;3|GMYg&+{JFbDncw@42pG z_?)Qk%Q4b&CT6{+uUtMALPCf*R&|C1lu%lq<*`DYd@(>t%lt8@1Mb>P!6u}{x; zf!~WYKi^q#o!-xgig;JMgLS(N7co6|J$_#u#Xh~aiuv@mZs3RdGb-{@Js7_D${mu~g{dfFOlh0XJ#W|ZaV`JpqnRt6a6`iGhPjkP<@;R%e z#t%*Kp8bJAArW@mug&(_Wtb7{Cb!5@gng#nqMov_ld`Ait99g zW%B&GQS8(F`kRP}--T{RNzWRU@%&mt+(-JgLPrwsSLkAzzaI1a8Ye#YulYgJpQ`WS z|InLj_e{SV|Y{>oCkn6pg`|vB*JB#&`1J4ij zxZX9n-;S_OTEO$epQ7I6le8i|6kdzx2QRL7Pwt;g&HNzs$z8d>G_9mxBAUwA`TT5` zV|itr@>q`5S2_o%;T{`o0|hkRQl@8^?M#O#Ym_vaUTyy4JEnkCFeT%Q|=~ z>$r+}f&4Q^k%x4x(1)zs?y>H55%-bKT`pqMwH0&1fBIs+;y%*3;UXqot0Q9ayR=yM z=8FA)%MYY`wb);|LH>%ls(&QkX&n0|E!ZD0V1K2E{Q+I}O;oyF@|Ui$?(f6C$vPx| zNyQJ6kJOU&{jcngjA9=|m;H|JtmCcO2N|Z3AIL8jk({pv!CM0zDW+(djR_>N7*+y%l=4RuJ;`7pL-hl zf%-&z56XR}zOoO}THfdLGfR%;mCuQ0KZ4FfclmzxmHZ0vb*p~4$`AQrcW!dKJNw-+ zu1Vi=uYG#z3D6@yCt*R$Mi{rz`+0O=#oWI&yNB=ZajxiDcFDf0eF279)k$sG`grBf z>kXTNlw+|>*%z%B+$>=9ewkj=?!f#1%r&Yd%K@V5!Vk~SySWYsZSyv z_wH8xYn_EJgWJt>{&Ju8H|~1Xu@~-g+^T!iLFhC2`j1u&YIAtmr8iA}c;7y6ZQW*}YjOE?Bi*jef5wxu zkInw<&>M^gdVYDX^0|Z8=%n`QHyJ&&TWPy1f5q+l%)Ep2e!@6xH$UxYU0hy2zIe66 z7d$!W>frewYTz6Tvqe{mbvKBcSB3u-Shi{iG#nP=lic6R;`_Y8+3K2&0JrAlNS~4Zl@K6erc;ylK+lJ_v2+( ztLPPD#KnZ_qn6A-zd0dhhwjJXj(tmQYbu}NzAh&WdR%>l8=|sZFZqg}gQ(XH{hrM~ zdO63f(!S@}Xp0xFrXba?e)c`O9QEXkz*Cn-uDap6k35)Khygl)z4{xBB=NTZc^cq z>Q}ERQ2k6ns$Z+ON@X>lv9}%v0o5-9r1~XxRp#9B=F{u>S+|L~07X&|1_1F8QW3At3R_xu!sT;Rg~=G`$Yge}{qi;RQ&% zun|G(KbIho_#nC;5DzQ|>HP0b2I{|NK%(cNplLC(3nEUxCi^oJ#w2 zeua;#eqD}1;_Z0Bmq5H->CXe=XW`R=#Lq#xavpgFiHEIw$g#Yl#Jd?nKg;{_3dFDb zjG{r@#a`;f_@2c3K|dcE!Xf)&CSm-C|QU^Bm4&sSs2 zlS?*rmI!|imm3UD{%b)KT-R^fHcy>E46N(4N+2Swr`mle)dRzSGR6+Z6GrWVe$K?!$yF-o{2!W3Z_OY3At zM4*KGwE`sEO@lzf{re|HY|n*{!uDF>Ke#iKo`Zxt`58#KPvPHil>OJToNPSqlH&DHgpwxd^D^TjYr1vP*yQ3*e z{dUR-q&{=D0P3%AAl18nSCsndxB*E0w5=Sdk7|Kb?@eui`e&0y+|(yaHC>nYfzD67 zi0S-P($OfLleeP(<@@Cor1KCXeh)+5msgPbexXKvE#IeUWq$ab`5{gHK3rvKWoeq( z7o!^toDtOdrhJ|FVg7gKhfwAR@pam>W=3m<2Sn4bG7^Fzg4P5jVVz7H3E z9viQBNAMO(^@#j~`JtJ(PkNs|74z?wl;9if_gRSRQk|A0GCvHJ_y6t6{18q2u=ED= z!)jR4aFWr1Uf1#1%=H~&T23c^7+P22hyCAQGCz0`KV;2=*(sAZJ}dMge&|$A{P6A_ z@k5i}aQ?i*wdZGEXMRW_erWwC^FtAC9dPAP%Qd1`j8Gs1?QhaESRK2dVG0h#xvFV1Bqr{IK2=r2e{FO8n3_i1^{e z3g!nLkoss(72=1Di-;e-Rhs|g=KzTxPJfVNO)K+*(4#>7a7X=q=7*{3`*O_uQ0eFP z;=0uL3x3E~^8M;7>*pb?pT+mTd>_vbYk7X?DE4Xo8pZnfPu9oc_vcWWUpw`Z?+2P6 zPK%gy?-|y;BjtUh`BhiM#P1bB4XnG7b2$l%Y}7nrQi3E@6)s* z&98Q(H~*Z>dUGDl4+iEuKinf-S#p4NWe)4g1k#gT-DrL|<4OAQ<08_JwUw+Jt!aL+ z_)0p_xhw0$&#VuJ{W}O$#I&Tfl5A7_){1A2daOre=E7o&O z6kFezCpSq=CLMRPC(RGh)f92%DJ|BvuwvbIiS$~;7n&a?`Y3|V#{TXp>$5uTNSC$V z!}Ei)^75a)=2p@9ti!xC;by%xgY}h4Z%O`6KF<$rczyuVQTKbZj!I%3l|lMxU|XIa zo{(OGpI9&Hl1^%7Mf1ZY1?!>|)DmVDpH-Ew zBl&Nq*nf*-|4lqEl7Dtp+%M_c)x{zv|E!e#v*u!-bnbK!li%ekV$y9L*%wO?`=oP& zMNB?c2N9EQ3t+!X%t7+?AFb?TjbeZ59s5&l$#?o|G5bzc$)8GeXMd^-`%|sick*I? z>38y-W%r`$fyi7izJdeIXn2 zdmJ;_@9D4H=CsT($7C@3LJJj9w@nYtQ28^G-?Q#I`8xTF6w&)`jF_p?A>{WoU|%PW z{hMKm=og{q4yZmyzRq#>Z+>E*#zYyXw|@AERVHW)S-<9_+I; zWIx7<{g^oRW7@ONGKBn;_3hbb8OJ`07yBvR?3>uLpW?uN$^-UIj<7$nhJBOPb|l;2GBZ`u{3$KA;ov!_O- z&qs{&`>naFAMTu|Q~!3G0QoxZyWQ(n#MTAq_xN44PG)-Oo?x-L+pYKFdg`em?)`Lx zZo&lD7RzP@311WwO$&^h)$4-Fc5&6#ZxQvuq?=#oj&Z5L89e6wjCtdNF)Adkb}h5j zBLBrzBL_5UxJtySPn%Y)?Uo_pw3tCnPDX2EMtZVK*<|tiRhU_S;3mVb6+d^0Ug%Wx z_nEesT+e0u;8)`P;w8K`q5DqJCwR`^=c%uD#m^Tr^7~C|;~I~bQ&;}+`lFBdb2Ggs zz4W#a-y?V$yqoVIsfFpif38=MGDg1d->!JEd%QvQS#N)IKgRbzm{B(RFviw=mD_w~ z2prK#sr`L%eGE5k;W*EG58CPJ1z!E~0b?)d^pAO=nqT}~tF|z-0FOrBZoPEP!>CfT z$GsB5VB?lo=}Wuh;;DOc&%GV?6%Q&~+_7}sh#pJhi}kd^FjLu}By4q8OmUystY_k!uui1{bq+1>1b33?5-o-p5yc^@;oFKcdq7n%$S z{oBa_uZ(pnw2Y~bIXC)FsFywmbIXT^Z)w&TGXmS2PC7Rk&yO`U?wm6T^M~zmQqC<@ z{e8W>0d-H7V4l^giFJb3V1AG;T(z*lD-(uY(f7%~tg0(~I?Sw%7msTjq}>(1mU`VV zPb;uj!t6@=0@3C2#(U%2aQ)o5eh0aJ^QeBG-e62#x58&rJ8=D~!qMtJpB!8Uas2|Q ze%;n$_@#cg#|=s1`nAJpH9ltFY30K8`-$sUj9H7VYMB&-a{Zoj{Yo&kvrYT^?YnUO zR&o7Was5=kSHAg8&ePJST)(MYzZ#fbqfytnUS(Xr;atCLu3rx3mu*S-)gp=OH=gU4 zgcqB{_ARqq$Mu`6UN<@B`iZ$g>{I=WnU7i#A9ZiQ{1eCgc1O*GjB{o;tNY8o;XGQ7yX_2A)R?a z#TQcl{mOh0PCT&V0iXXZ+<$ZU+%H1uyqo7!|JjuCIoJI2SUS%S^Skub9 z{Z78lyxmTY<(2u_M8w3;mH%Ag;WqMp@`}W}>*ZK|MTuXdN-OOXpH})Fk?%v|P2p>c zeTgSWh@bDu`~T~TK2v%nwrDpL*Unqs`(n>ZzJY4~8mMS9I}1JcIeOG{TPOc{xbboB z!yV%nplhP-FYlkHqg_^&<1Oy*7uWq9dh}?ySNH}P*{7|xwt?{HaaTgw`fby%;ocsm zgFDO}C_Wct-?r@j-d7?{SnIqx{B@d${Tf`YHp9df1K$M1GIHhIp=tc$KQ9X$Gu_r+cFKE*da-43_LJTq!jw<$XR*}bLh&>#%%?C|~C30sVJ zHhkXPZXX_Rf2Om`&IWkYtp9@#M$<6L#M?v1$#62x!&uzf7Ytw zm1FJ`@%s((zWU1Nr@}Xq&d)tDAIkUhIoTk`>MPQD&^R}t)b}GVs((9p;8#QrO2G>X%pV+^hoEL~6MgZ}{WT590IDp=d__>@iGGzLcsC9wj^TZGS zAo0VT*J7XEoAA|?$sNB08#e4Gu2Ve*0KG?{pMl<|pM0NEd1i@Dx4U!X7^qHXSPT+B6h(r>4{v%Z59cFJd1l_-aSuQEfW!~6Ng(k zeyE~U*6%VZxcVrN_+dvp5I-~k;)hd8d>xiIG^7eh{9yS3h#yRV_~D%LT&Jg1w>g#q z@k2i#e%P*#SB`=B;V9#NK)By^QfwGJ|3!O;!HT&4jpC9nlmhX?V!~aw2@pR-67J*W zimlI|T>ZA)28bWr7`Hy-KBP>WTi8hd?`uH(u!V7VX57b=xknwIFI(ONBz{=>k#P(C zBIA`;#yyDp(~SDFXNkgbV!_e7zRrq}dYejTMs(nM7jYlnrg{%+#raX`{X3Xa8Yb?}6Mu9x_h#mHVWUk6r0L zIzPkYeL6q)ZRA*9fzFAq9ILND=V6fgz8rJkKak`9b_LSU;XwL%i~8qPKFB}wLl8*w z!{+Pqbs&AL2a-M>^%-dX68;yE9-blYr}-gQ#58}o1L@vLVxQ*M1Q8Rz?-eo4Uv+`> zY?Rn1{kl)YG(TJrG0k63fOM+h8~OT=Rv`Vk4oGi~2GX0YK$;&i0)Xa+zCgM%9!OXA z1Jadsfb?V`(EKn1B>i}!3`jpN0Md=qf#!z;Kyx|(>BPZ6`Y;Y?epmsL9<+G{qz6Z9 z!Vjcd8ZqvXg!|i3#+^*K&36Lnxvh+QIY@f0Z6V_>VBD5Kx~&Q0p3S%?Fm7KUeRhj+ z_h#JTjC(GS4x7TbjsJh#<6d%qW^sSM<$jE%esuHVdS9YGG}EVgS2f`NQ~mpF(%hQG z^?t(rW)0M5p+NoR#`W&OeH9MWPu@U%beZdYi~FamM%>gV8g+&IIX*w*L`>%=TaM)w z=$uqKm*To~9umIDvHA+s_rf<)zy6~b$bV}LlK<8{PyK!%UArF0KNCIf65spe`-l@Am7OsCEux=4Uj)&4U#|AJrBvB`V8bdbprC2 z9s&7Ien7s{6Ci)7Ban}z3*;}=1M-*t1oDwefczpQkdNdD4 zFz$ndyYovRU#Jb?HZ()>g-jTCQ^tJ_*zaN7j~RCx#$5x**J;MMr!ekJ#vKXd-<)OK znt6?IXL5gbCa~?voNt??c|F^OL3}Vmdz^bt8<2u&`QWuPe+VX<{8}{n zy9r+Kuetfj%ZpeDCu4i;Q~i4juWOWlxl{x1_OKpVw_y&R+wyRw;ix5edEw^7`_0nv zR@i|xmhWuv$-+6=AG>YDhh^)}T)gMX`Aj)q0Oy;>`C4$kFCgJNWzP5_8Q*foH-+&n zh0}Lt`+HetV8#m1j=3eH@q(ABPs^#V@zT*%CDuxJJX0K4*E9YIrZd@3JU!skW!g1<06ea7d+@8css?+||8JN&#M z{Jh7xe%H8uBe{NoT)#qoAFBM)`#8hzqlDkbRDK^${62oBe1%$^PoMK0?&qo8&o#N9i}`%0{8;IH z_^X`{Z$2N(IiDTpOXGagIG^h8aZBgz7N3ufj4y`qSu(ycjIX!q@Ao-ukJz{dkJWYB z8{_r`BX+IW+3%-VOdj7fE#r6?CN=KbEWkb%f}WK|FT8mhcKHT+ME`scw(Qvv+A*{a z9CZj@eD(f%*nj0}LH5+4kp0pzEA)T`1a6*mEcTTlBv!9#ax~i>GGuJDy}-Zb@-grBnu3pES(Bt|MlkE zUq)=#8U{x!`h>i38w=^b=@@3b55xW2wM#5_b%1rVm$mHm@i@e;DQY#L-eyQTrc>}o zsUGL+L;3Q*P(JOWlu!E-dPeCG&Xsp{wL|KPJLW_%umZ^bu$-W~it{ONg}Ht_R0 z^YeOg{fwx71-VqeTJ}`GD$)Et9Ql0&^ZT%-_wm7x--i>ukC%mJp zl+SY<<0~e7(;pJPO`eQz4dd&qc0TU$`H12EJIMW)!To%W&xbDca|rVJu%Ujop2+9J zn9fJ(Cq5t3_3k?S-xWR|O3pW)&&MgsH~lc<+r#JM37-!W#uvx_gY# z{_P%dO?!J#zAkks-{R_=&w%sAa=vqfFS$PBi)VbHjITfAn+~T&DIa~ZQ~iHq!i@{t zjLv|}eTELB);xwQ^<8&gZuSSHjyOB6+#>=Kl3Z*OJ)S|{1~-ot$t6&t|44hpJRNvw zdA;GF6AR&*o%^x-o900D!%r3^s{$c^R$1_*!eMZ4ghRyE1Eb+_*~u!66GEX#*Zq0- zR;M^$9_PEo`Koh1OU^f$@r_}8ml)qG#%Ip>Tp|3qZD_#uI7pq`DC)q@){tX8F=S9D}r zTpu%{>9N(JkYau>RX;-yt~+)M9=LN5JZKiMp?BSGs(%l0thI^TdAQ;B{N^P8Sk7m` z`BrnjE}YMk^Brbe=4a`B6mtFQ zas8HY{p|RC`1AYtmEVWz-$Rw&hmzmNSB2{^!h&xt4*)~ diff --git a/tests/stokes/manufactured.i b/tests/stokes/patch.i similarity index 98% rename from tests/stokes/manufactured.i rename to tests/stokes/patch.i index 9ab5ca0..31c574b 100644 --- a/tests/stokes/manufactured.i +++ b/tests/stokes/patch.i @@ -77,6 +77,10 @@ [equil] type = ADStokesStressDivergence variable = u + [] + [pressure] + type = ADStokesPressure + variable = u pressure = p [] [incompressible] diff --git a/tests/stokes/patch_by.i b/tests/stokes/patch_by.i new file mode 100644 index 0000000..5628db9 --- /dev/null +++ b/tests/stokes/patch_by.i @@ -0,0 +1,182 @@ +[Mesh] + [mesh] + type = GeneratedMeshGenerator + dim = 3 + xmin = 0 + xmax = 1 + ymin = 0 + ymax = 1 + zmin = 0 + zmax = 1 + nx = 5 + ny = 5 + nz = 5 + elem_type = HEX20 + [] +[] + +[Variables] + [p] + order = FIRST + family = LAGRANGE + [] + [u] + order = SECOND + family = LAGRANGE_VEC + [] +[] + +[AuxVariables] + [u_exact] + order = FIRST + family = LAGRANGE_VEC + [] + [p_exact] + order = FIRST + family = LAGRANGE + [] +[] + +[AuxKernels] + [u_exact] + type = VectorFunctionAux + variable = u_exact + function = u_exact + [] + [p_exact] + type = FunctionAux + variable = p_exact + function = p_exact + [] +[] + +[Functions] + [u_exact] + type = ParsedVectorFunction + expression_x = '2 * x^2 + y^2 + z^2' + expression_y = '2 * x^2 - 2 * x * y' + expression_z = '2 * x^2 - 2 * x * z' + [] + [p_exact] + type = ParsedFunction + expression = 'x + y + z - 3.0/2.0' + [] +[] + +[ICs] + [u] + type = VectorConstantIC + variable = u + x_value = 1e-15 + y_value = 1e-15 + z_value = 1e-15 + [] +[] + +[Kernels] + [equil] + type = ADStokesStressDivergence + variable = u + [] + [pressure] + type = ADStokesPressureByParts + variable = u + pressure = p + [] + [incompressible] + type = ADStokesIncompressibility + variable = p + velocity = u + [] + [body] + type = VectorBodyForce + variable = u + function_x = -7 + function_y = -3 + function_z = -3 + [] +[] + +[Materials] + [strain] + type = StokesStrainRate + velocity = u + [] + [stress] + type = StokesLinearViscous + mu = 1.0 + [] +[] + +[BCs] + [u_fixed] + type = ADVectorFunctionDirichletBC + variable = u + boundary = 'bottom top back front left right' + function_x = '2 * x^2 + y^2 + z^2' + function_y = '2 * x^2 - 2 * x * y' + function_z = '2 * x^2 - 2 * x * z' + [] + [p_fixed] + type = FunctionDirichletBC + variable = p + boundary = 'bottom top back front left right' + function = 'p_exact' + [] +[] + +# [Preconditioning] +# [FSP] +# type = FSP +# topsplit = 'up' +# [up] +# splitting = 'u p' +# splitting_type = schur +# petsc_options_iname = '-pc_fieldsplit_schur_fact_type -pc_fieldsplit_schur_precondition -ksp_gmres_restart -ksp_type -ksp_pc_side' +# petsc_options_value = 'full self 300 fgmres right' +# [] +# [u] +# vars = 'u' +# petsc_options_iname = '-pc_type' +# petsc_options_value = 'lu' +# [] +# [p] +# vars = 'p' +# petsc_options = '-pc_lsc_scale_diag' +# petsc_options_iname = '-ksp_type -ksp_gmres_restart -pc_type -ksp_pc_side -lsc_pc_type -lsc_ksp_type -lsc_ksp_pc_side -lsc_ksp_gmres_restart' +# petsc_options_value = 'fgmres 300 lsc right lu gmres right 300' +# [] +# [] +# [] + +[Executioner] + type = Steady + + solve_type = 'newton' + petsc_options_iname = '-pc_type' + petsc_options_value = 'lu' + line_search = none + + l_max_its = 10 + l_tol = 1e-8 + nl_max_its = 15 + nl_rel_tol = 1e-10 + nl_abs_tol = 1e-12 +[] + +[Postprocessors] + [u_error] + type = ElementVectorL2Error + function = u_exact + variable = u + [] + [p_error] + type = ElementL2Error + function = p_exact + variable = p + [] +[] + +[Outputs] + exodus = true +[] diff --git a/tests/stokes/structural.i b/tests/stokes/structural.i index 1648fcf..0a40b80 100644 --- a/tests/stokes/structural.i +++ b/tests/stokes/structural.i @@ -6,9 +6,9 @@ dx = '1 1' dy = '1 1' dz = '2' - ix = '1 1' - iy = '1 1' - iz = '2' + ix = '2 2' + iy = '2 2' + iz = '4' subdomain_id = '1 2 3 4' [] [] @@ -16,7 +16,7 @@ [Variables] [p] order = FIRST - family = MONOMIAL + family = LAGRANGE [] [u] order = SECOND @@ -25,6 +25,30 @@ [] [AuxVariables] + [deviatoric_stress_xx] + family = MONOMIAL + order = CONSTANT + [] + [deviatoric_stress_yy] + family = MONOMIAL + order = CONSTANT + [] + [deviatoric_stress_zz] + family = MONOMIAL + order = CONSTANT + [] + [deviatoric_stress_xy] + family = MONOMIAL + order = CONSTANT + [] + [deviatoric_stress_yz] + family = MONOMIAL + order = CONSTANT + [] + [deviatoric_stress_xz] + family = MONOMIAL + order = CONSTANT + [] [stress_xx] family = MONOMIAL order = CONSTANT @@ -37,39 +61,108 @@ family = MONOMIAL order = CONSTANT [] + [stress_xy] + family = MONOMIAL + order = CONSTANT + [] + [stress_yz] + family = MONOMIAL + order = CONSTANT + [] + [stress_xz] + family = MONOMIAL + order = CONSTANT + [] [] [AuxKernels] - [stress_xx] + [deviatoric_stress_xx] type = ADMaterialRankTwoTensorAux - variable = stress_xx - property = 'stress' + variable = deviatoric_stress_xx + property = 'deviatoric_stress' i = 0 j = 0 [] - [stress_yy] + [deviatoric_stress_yy] type = ADMaterialRankTwoTensorAux - variable = stress_yy - property = 'stress' + variable = deviatoric_stress_yy + property = 'deviatoric_stress' i = 1 j = 1 [] - [stress_zz] + [deviatoric_stress_zz] type = ADMaterialRankTwoTensorAux - variable = stress_zz - property = 'stress' + variable = deviatoric_stress_zz + property = 'deviatoric_stress' i = 2 j = 2 [] + [deviatoric_stress_xy] + type = ADMaterialRankTwoTensorAux + variable = deviatoric_stress_xy + property = 'deviatoric_stress' + i = 0 + j = 1 + [] + [deviatoric_stress_yz] + type = ADMaterialRankTwoTensorAux + variable = deviatoric_stress_yz + property = 'deviatoric_stress' + i = 1 + j = 2 + [] + [deviatoric_stress_xz] + type = ADMaterialRankTwoTensorAux + variable = deviatoric_stress_xz + property = 'deviatoric_stress' + i = 0 + j = 2 + [] + [stress_xx] + type = ParsedAux + variable = stress_xx + coupled_variables = 'p deviatoric_stress_xx' + expression = 'deviatoric_stress_xx - p' + [] + [stress_yy] + type = ParsedAux + variable = stress_yy + coupled_variables = 'p deviatoric_stress_yy' + expression = 'deviatoric_stress_yy - p' + [] + [stress_zz] + type = ParsedAux + variable = stress_zz + coupled_variables = 'p deviatoric_stress_zz' + expression = 'deviatoric_stress_zz - p' + [] + [stress_xy] + type = ParsedAux + variable = stress_xy + coupled_variables = 'deviatoric_stress_xy' + expression = 'deviatoric_stress_xy' + [] + [stress_xz] + type = ParsedAux + variable = stress_xz + coupled_variables = 'deviatoric_stress_xz' + expression = 'deviatoric_stress_xz' + [] + [stress_yz] + type = ParsedAux + variable = stress_yz + coupled_variables = 'deviatoric_stress_yz' + expression = 'deviatoric_stress_yz' + [] [] [ICs] [u] - type = VectorConstantIC + type = VectorFunctionIC variable = u - x_value = 1e-15 - y_value = 1e-15 - z_value = 1e-15 + function_x = '1.0 * x' + function_y = '1.0 * y' + function_z = '1.0 * z' [] [] @@ -77,6 +170,10 @@ [equil] type = ADStokesStressDivergence variable = u + [] + [pressure] + type = ADStokesPressure + variable = u pressure = p [] [incompressible] @@ -93,7 +190,7 @@ [] [stress] type = StokesLinearViscous - mu = 1.0 + mu = 1 [] [] @@ -129,13 +226,11 @@ value = -1.0 [] - [test] - type = ADVectorFunctionNeumannBC - boundary = 'front' - variable = u - function_x = '0' - function_y = '0' - function_z = '2.0/3.0' + [fix_zero] + type = ADDirichletBC + variable = p + boundary = 'right top' + value = 0.0 [] [] @@ -143,8 +238,8 @@ type = Steady solve_type = 'newton' - petsc_options_iname = '-pc_type' - petsc_options_value = 'svd' + petsc_options_iname = '-pc_type -pc_factor_shift_type' + petsc_options_value = 'lu NONZERO' line_search = none l_max_its = 10 diff --git a/tests/stokes/tests b/tests/stokes/tests index 3ca9200..5659e77 100644 --- a/tests/stokes/tests +++ b/tests/stokes/tests @@ -1,8 +1,20 @@ [Tests] - [./flow] + [./patch_no_by_parts] type = Exodiff - input = 'flow.i' - exodiff = 'flow_out.e' - requirement = "Stokes flow gives the right result for simple creeping flow" + input = 'patch.i' + exodiff = 'patch_out.e' + requirement = "Stokes flow gives the right result for a simple patch test, no integration by parts" + [../] + [./patch_by_parts] + type = Exodiff + input = 'patch_by.i' + exodiff = 'patch_by_out.e' + requirement = "Stokes flow gives the right result for a simple patch test, with integration by parts" + [../] + [./structural] + type = Exodiff + input = 'structural.i' + exodiff = 'structural_out.e' + requirement = "Stokes flow gives the right result for a structural problem (no integration by parts)" [../] []