From 42a47dae312642eb8527ce14b2f9886c4d60df29 Mon Sep 17 00:00:00 2001 From: mleoni-pf <160495781+mleoni-pf@users.noreply.github.com> Date: Tue, 4 Feb 2025 09:03:52 +0100 Subject: [PATCH 1/5] Update elasticity_scaling.md (#242) MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit * Improve workflow handling (#118) * Create common action for all workflows and split them into multiple files * Dokken/update fspace and wmtgs (#127) * fix runs-on * Minor improvements in linearelasticity.md (#124) * Fix various things on release branch * Try adding trame * Revert fundamentals * Membrane code * More updates due to autoformatting * Lagrange -> CG back conversion * Change to github image --------- Co-authored-by: Julius Herb <43179176+juliusgh@users.noreply.github.com> * Update petsc arches * Fix bounding boxes and more (#135) * Start fixing code for nightly build * rerun all files * Try fixing subdomains * Fix meshio * PETSc python API update (#137) * Resolve https://github.com/jorgensd/dolfinx-tutorial/issues/136 * Fixes related to: https://github.com/FEniCS/dolfinx/pull/2703 * Update changelog * Make sure all notebooks run. Change to pathlib in some examples * Api changes related to: https://github.com/FEniCS/dolfinx/pull/2763 (#142) * Bump version numbers * Bump version numbers (#149) * Fix petsc arch * Prepare v0.7.1 (#154) * Merge branch 'main' into release and bump versions * Delete obsolete file * Temporary shift path in nightly test * Make sure there are no warnigns in build * Tabulation * Updates compatible with nightly branch of DOLFINx (#156) * Replace VectorElement and FiniteElement with `basix.ufl.element`. Replace `dolfinx.fem.FunctionSpace` with `dolfinx.fem.functionspace` * More updates * And more * Last vector element updates * Bump pyvista and dolfinx to v0.7.2 (#159) * Bump pyvista and dolfinx to v0.7.2 * remove -r * Add further pyvista deps * Bump versions * Update petsc solver and pc links. Resolves #143 (#160) * Fix typo and issue Issue on page /fem.html #122 (#161) * Compute entity-cell connectivity before calling locate_dofs_topological. Remove soon to be deprecated pyvista syntax for updating time dependent fields. Fix range->np.arange conversion * Update backend in workflow * Fix broken link (#140) + http:// to https:// (#162) * Run book build prior to parallel run * Add back write frame * Fix wrong link * Fix typo in changelog * Merge main into release (#189) * Merge main into release * Apply suggestions from code review * Various fixes when reading through the diff * Further fixes * add missing checkout (#191) * add missing checkout * Add test docker on release PR and release (only push on tag * Remove unused import * Update text to resolve #194 (#195) * Dokken/update nonlin options (#203) * Update options for non-linear solver to sensible choice, ref #200 * Update docker path * Update python file as well * Update config ptr (#207) * Pressure correction equation fix (#196) * Update fundamentals.md (#199) on line 110, updated a typesetting error in equation (3) * Change from vector to petsc_vec (#206) * dolfinx.fem.Form changed to dolfinx.fem.form (#213) * Update compiler_parameters.ipynb In the text it was mentioned that `dolfinx.fem.Form` was used to compile the form, instead of `dolfinx.fem.form`. Not matching the code cell. * jupytext sync * Form -> form in both py and ipynb and LinearProblem -> fem.petsc.LinearProblem * Update book_stable.yml * Update test_stable.yml * Fix python path in test * Fix CI (#227) * Bump docker/build-push-action from 5 to 6 (#221) Bumps [docker/build-push-action](https://github.com/docker/build-push-action) from 5 to 6. - [Release notes](https://github.com/docker/build-push-action/releases) - [Commits](https://github.com/docker/build-push-action/compare/v5...v6) --- updated-dependencies: - dependency-name: docker/build-push-action dependency-type: direct:production update-type: version-update:semver-major ... Signed-off-by: dependabot[bot] Co-authored-by: dependabot[bot] <49699333+dependabot[bot]@users.noreply.github.com> * Bump actions/configure-pages from 4 to 5 (#222) Bumps [actions/configure-pages](https://github.com/actions/configure-pages) from 4 to 5. - [Release notes](https://github.com/actions/configure-pages/releases) - [Commits](https://github.com/actions/configure-pages/compare/v4...v5) --- updated-dependencies: - dependency-name: actions/configure-pages dependency-type: direct:production update-type: version-update:semver-major ... Signed-off-by: dependabot[bot] Co-authored-by: dependabot[bot] <49699333+dependabot[bot]@users.noreply.github.com> * Update book_stable.yml (#225) * Update book_stable.yml * Update test_stable.yml * resolve #224 (#226) * Update test_stable.yml --------- Signed-off-by: dependabot[bot] Co-authored-by: dependabot[bot] <49699333+dependabot[bot]@users.noreply.github.com> * Update test_stable.yml (#235) * Add libgl flag to ci (#238) * Update elasticity_scaling.md --------- Signed-off-by: dependabot[bot] Co-authored-by: Jørgen Schartum Dokken Co-authored-by: Jørgen S. Dokken Co-authored-by: Julius Herb <43179176+juliusgh@users.noreply.github.com> Co-authored-by: rossbm1 <120818149+rossbm1@users.noreply.github.com> Co-authored-by: Manuel Pena Co-authored-by: dependabot[bot] <49699333+dependabot[bot]@users.noreply.github.com> --- chapter2/elasticity_scaling.md | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/chapter2/elasticity_scaling.md b/chapter2/elasticity_scaling.md index da15c47d..e6a0b2f2 100644 --- a/chapter2/elasticity_scaling.md +++ b/chapter2/elasticity_scaling.md @@ -18,6 +18,6 @@ where $\beta = 1+\frac{\lambda}{\mu}$ is a dimensionless elasticity parameter an ``` is a dimensionless variable reflecting the ratio of the load $\rho g$ and the shear stress term $\mu \nabla^2u \sim \mu \frac{U}{L^2}$ in the PDE. -One option for the scaling is to chose $U$ such that $\gamma$ is of unit size ($U=\frac{\rho g L^2}{\mu}$). However, in elasticity, this leads to displacements of the size of the geometry. This can be achieved by choosing $U$ equal to the maximum deflection of a clamped beam, for which there actually exists a formula: $U=\frac{3}{2} \rho g L^2\frac{\delta^2}{E}$ where $\delta=\frac{L}{W}$ is a parameter reflecting how slender the beam is, and $E$ is the modulus of elasticity. Thus the dimensionless parameter $\delta$ is very important in the problem (as expected $\delta\gg 1$ is what gives beam theory!). Taking $E$ to be of the same order as $\mu$, which in this case and for many materials, we realize that $\gamma \sim \delta^{-2}$ is an appropriate choice. Experimenting with the code to find a displacement that "looks right" in the plots of the deformed geometry, points to $\gamma=0.4\delta^{-2}$ as our final choice of $\gamma$. +One option for the scaling is to choose $U$ such that $\gamma$ is of unit size ($U=\frac{\rho g L^2}{\mu}$). However, in elasticity, this leads to displacements of the size of the geometry. This can be achieved by choosing $U$ equal to the maximum deflection of a clamped beam, for which there actually exists a formula: $U=\frac{3}{2} \rho g L^2\frac{\delta^2}{E}$ where $\delta=\frac{L}{W}$ is a parameter reflecting how slender the beam is, and $E$ is the modulus of elasticity. Thus the dimensionless parameter $\delta$ is very important in the problem (as expected $\delta\gg 1$ is what gives beam theory!). Taking $E$ to be of the same order as $\mu$, which in this case and for many materials, we realize that $\gamma \sim \delta^{-2}$ is an appropriate choice. Experimenting with the code to find a displacement that "looks right" in the plots of the deformed geometry, points to $\gamma=0.4\delta^{-2}$ as our final choice of $\gamma$. -The simulation code implements the problem with dimensions and physical parameters $\lambda, \mu, \rho, g, L$ and $W$. However, we can easily reuse this code for a scaled problem: Just set $\mu=\rho=L=1$, $W$ as $W/L(\delta^{-1})$, $g=\gamma$ and $\lambda=\beta$. \ No newline at end of file +The simulation code implements the problem with dimensions and physical parameters $\lambda, \mu, \rho, g, L$ and $W$. However, we can easily reuse this code for a scaled problem: Just set $\mu=\rho=L=1$, $W$ as $W/L(\delta^{-1})$, $g=\gamma$ and $\lambda=\beta$. From ce1216597d2c94f773a233383f19f5dcf1122c09 Mon Sep 17 00:00:00 2001 From: mleoni-pf <160495781+mleoni-pf@users.noreply.github.com> Date: Tue, 4 Feb 2025 09:01:44 +0100 Subject: [PATCH 2/5] Update navierstokes.md (#247) MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit * Improve workflow handling (#118) * Create common action for all workflows and split them into multiple files * Dokken/update fspace and wmtgs (#127) * fix runs-on * Minor improvements in linearelasticity.md (#124) * Fix various things on release branch * Try adding trame * Revert fundamentals * Membrane code * More updates due to autoformatting * Lagrange -> CG back conversion * Change to github image --------- Co-authored-by: Julius Herb <43179176+juliusgh@users.noreply.github.com> * Update petsc arches * Fix bounding boxes and more (#135) * Start fixing code for nightly build * rerun all files * Try fixing subdomains * Fix meshio * PETSc python API update (#137) * Resolve https://github.com/jorgensd/dolfinx-tutorial/issues/136 * Fixes related to: https://github.com/FEniCS/dolfinx/pull/2703 * Update changelog * Make sure all notebooks run. Change to pathlib in some examples * Api changes related to: https://github.com/FEniCS/dolfinx/pull/2763 (#142) * Bump version numbers * Bump version numbers (#149) * Fix petsc arch * Prepare v0.7.1 (#154) * Merge branch 'main' into release and bump versions * Delete obsolete file * Temporary shift path in nightly test * Make sure there are no warnigns in build * Tabulation * Updates compatible with nightly branch of DOLFINx (#156) * Replace VectorElement and FiniteElement with `basix.ufl.element`. Replace `dolfinx.fem.FunctionSpace` with `dolfinx.fem.functionspace` * More updates * And more * Last vector element updates * Bump pyvista and dolfinx to v0.7.2 (#159) * Bump pyvista and dolfinx to v0.7.2 * remove -r * Add further pyvista deps * Bump versions * Update petsc solver and pc links. Resolves #143 (#160) * Fix typo and issue Issue on page /fem.html #122 (#161) * Compute entity-cell connectivity before calling locate_dofs_topological. Remove soon to be deprecated pyvista syntax for updating time dependent fields. Fix range->np.arange conversion * Update backend in workflow * Fix broken link (#140) + http:// to https:// (#162) * Run book build prior to parallel run * Add back write frame * Fix wrong link * Fix typo in changelog * Merge main into release (#189) * Merge main into release * Apply suggestions from code review * Various fixes when reading through the diff * Further fixes * add missing checkout (#191) * add missing checkout * Add test docker on release PR and release (only push on tag * Remove unused import * Update text to resolve #194 (#195) * Dokken/update nonlin options (#203) * Update options for non-linear solver to sensible choice, ref #200 * Update docker path * Update python file as well * Update config ptr (#207) * Pressure correction equation fix (#196) * Update fundamentals.md (#199) on line 110, updated a typesetting error in equation (3) * Change from vector to petsc_vec (#206) * dolfinx.fem.Form changed to dolfinx.fem.form (#213) * Update compiler_parameters.ipynb In the text it was mentioned that `dolfinx.fem.Form` was used to compile the form, instead of `dolfinx.fem.form`. Not matching the code cell. * jupytext sync * Form -> form in both py and ipynb and LinearProblem -> fem.petsc.LinearProblem * Update book_stable.yml * Update test_stable.yml * Fix python path in test * Fix CI (#227) * Bump docker/build-push-action from 5 to 6 (#221) Bumps [docker/build-push-action](https://github.com/docker/build-push-action) from 5 to 6. - [Release notes](https://github.com/docker/build-push-action/releases) - [Commits](https://github.com/docker/build-push-action/compare/v5...v6) --- updated-dependencies: - dependency-name: docker/build-push-action dependency-type: direct:production update-type: version-update:semver-major ... Signed-off-by: dependabot[bot] Co-authored-by: dependabot[bot] <49699333+dependabot[bot]@users.noreply.github.com> * Bump actions/configure-pages from 4 to 5 (#222) Bumps [actions/configure-pages](https://github.com/actions/configure-pages) from 4 to 5. - [Release notes](https://github.com/actions/configure-pages/releases) - [Commits](https://github.com/actions/configure-pages/compare/v4...v5) --- updated-dependencies: - dependency-name: actions/configure-pages dependency-type: direct:production update-type: version-update:semver-major ... Signed-off-by: dependabot[bot] Co-authored-by: dependabot[bot] <49699333+dependabot[bot]@users.noreply.github.com> * Update book_stable.yml (#225) * Update book_stable.yml * Update test_stable.yml * resolve #224 (#226) * Update test_stable.yml --------- Signed-off-by: dependabot[bot] Co-authored-by: dependabot[bot] <49699333+dependabot[bot]@users.noreply.github.com> * Update test_stable.yml (#235) * Add libgl flag to ci (#238) * Update navierstokes.md --------- Signed-off-by: dependabot[bot] Co-authored-by: Jørgen Schartum Dokken Co-authored-by: Jørgen S. Dokken Co-authored-by: Julius Herb <43179176+juliusgh@users.noreply.github.com> Co-authored-by: rossbm1 <120818149+rossbm1@users.noreply.github.com> Co-authored-by: Manuel Pena Co-authored-by: dependabot[bot] <49699333+dependabot[bot]@users.noreply.github.com> --- chapter2/navierstokes.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/chapter2/navierstokes.md b/chapter2/navierstokes.md index dd23acc3..336562cc 100644 --- a/chapter2/navierstokes.md +++ b/chapter2/navierstokes.md @@ -59,7 +59,7 @@ $\langle -\nabla \cdot \sigma, v\rangle$. Just as for the [linear elasticity pro where $T=\sigma \cdot n$ is the boundary traction. If we solve a problem with a free boundary, we can take $T=0$ on the boundary. However, if we compute the flow through a channel or a pipe and want to model flow that continues into an "imaginary channel" at the outflow, we need to treat this term with some care. The assumption we then can make is that the derivative of the velocity in the direction of the channel is zero at the outflow, corresponding to that the flow is "fully developed" or doesn't change significantly downstream at the outflow. Doing so, the remaining boundary term at the outflow becomes -$pn - \mu \nabla u \cdot n$, which is the term appearing in the variational problem [](ipcs-one). Note that this argument and the implementation depends exact on the definition of $\nabla u$, as either the matrix with components $\frac{\partial u_i}{\partial x_j}$ or $\frac{\partial u_j}{\partial x_i}$. +$pn - \mu \nabla u \cdot n$, which is the term appearing in the variational problem [](ipcs-one). Note that this argument and the implementation depend exactly on the definition of $\nabla u$, as either the matrix with components $\frac{\partial u_i}{\partial x_j}$ or $\frac{\partial u_j}{\partial x_i}$. We here choose the latter, $\frac{\partial u_j}{\partial x_i}$, which means that we must use the UFL-operator `nabla_grad`. If we use the operator `grad` and the definition $\frac{\partial u_i}{\partial x_j}$, we must instead keep the terms $pn-\mu(\nabla u)^T \cdot n$. From 459ee5ba5191a9932ec38b64efa95843cd20ff81 Mon Sep 17 00:00:00 2001 From: Antonio Baiano Svizzero <107617271+bayswiss@users.noreply.github.com> Date: Mon, 3 Feb 2025 21:29:34 +0100 Subject: [PATCH 3/5] Adding Acoustic Helmholtz equation to the tutorial (#232) * Acoustic Helmholtz scripts helmholtz.md contains the mathematical derivation of the variational formulation, while the helmholtz_code scripts contain the implementation. * Fix * equations rendering bug fixed * Update helmholtz.md Equations rendering second fix, missing spaces * Major updates to helmholtz * Add to toc * Updates to helmholtz * Update to work in parallel --------- Co-authored-by: jorgensd --- _toc.yml | 6 +- chapter2/helmholtz.md | 125 ++++++++ chapter2/helmholtz_code.ipynb | 519 ++++++++++++++++++++++++++++++++++ chapter2/helmholtz_code.py | 290 +++++++++++++++++++ 4 files changed, 939 insertions(+), 1 deletion(-) create mode 100644 chapter2/helmholtz.md create mode 100644 chapter2/helmholtz_code.ipynb create mode 100644 chapter2/helmholtz_code.py diff --git a/_toc.yml b/_toc.yml index 3f46be09..42eb9934 100644 --- a/_toc.yml +++ b/_toc.yml @@ -35,6 +35,10 @@ parts: - file: chapter2/ns_code1 - file: chapter2/ns_code2 - file: chapter2/hyperelasticity + - file: chapter2/helmholtz + sections: + - file: chapter2/helmholtz_code + - caption: Subdomains and boundary conditions chapters: - file: chapter3/neumann_dirichlet_code @@ -48,4 +52,4 @@ parts: - file: chapter4/solvers - file: chapter4/compiler_parameters - file: chapter4/convergence - - file: chapter4/newton-solver \ No newline at end of file + - file: chapter4/newton-solver diff --git a/chapter2/helmholtz.md b/chapter2/helmholtz.md new file mode 100644 index 00000000..00d6d1ac --- /dev/null +++ b/chapter2/helmholtz.md @@ -0,0 +1,125 @@ +# The Helmholtz equation +Author: Antonio Baiano Svizzero + +The study of computational acoustics is fundamental in fields such as noise, vibration, and harshness (NVH), noise control, and acoustic design. In this chapter, we focus on the theoretical foundations of the Helmholtz equation - valid for noise problems with harmonic time dependency - and its implementation in FEniCSx to compute the sound pressure for any acoustic system. + +## The PDE problem +The acoustic Helmholtz equation in its general form reads + +$$ +\begin{align} +\nabla^2 p + k^2 p = -j \omega \rho_0 q \qquad\text{in } \Omega, +\end{align} +$$ + +where $k$ is the acoustic wavenumber, $\omega$ is the angular frequency, $j$ the imaginary unit and $q$ is the volume velocity ($m^3/s$) of a generic source field. +In case of a monopole source, we can write $q=Q \delta(x_s,y_s,z_s)$, where $\delta(x_s,y_s,z_s)$ is the 3D Dirac Delta centered at the monopole location. + +This equation is coupled with the following boundary conditions: + +- Dirichlet BC: + + $$ + \begin{align} + p = \bar{p} \qquad \text{on } \partial\Omega_p, + \end{align} + $$ + +- Neumann BC: + + $$ + \begin{align} + \frac{\partial p}{\partial n} = - j \omega \rho_0 \bar{v}_n\qquad \text{on } \partial\Omega_v, + \end{align} + $$ + +- Robin BC: + + $$ + \begin{align} + \frac{\partial p}{\partial n} = - \frac{j \omega \rho_0 }{\bar{Z}} p \qquad \text{on } \partial\Omega_Z, + \end{align} + $$ + +where we prescribe, respectively, an acoustic pressure $\bar{p}$ on the boundary $\partial\Omega_p$, +a sound particle velocity $\bar{v}_n$ on the boundary $\partial\Omega_v$ and +an acoustic impedance $\bar{Z}$ on the boundary $\partial\Omega_Z$ where $n$ is the outward normal. +In general, any BC can also be frequency dependant, as it happens in real-world applications. + +## The variational formulation +Now we have to turn the equation in its weak formulation. +The first step is to multiplicate the equation by a *test function* $v\in \hat V$, +where $\hat V$ is the *test function space*, after which we integrate over the whole domain, $\Omega$: + +$$ +\begin{align} +\int_{\Omega}\left(\nabla^2 p + k^2 p \right) \bar v ~\mathrm{d}x = -\int_{\Omega} j \omega \rho_0 q \bar v ~\mathrm{d}x. +\end{align} +$$ + +Here, the unknown function $p$ is referred to as *trial function* and the $\bar{\cdot}$ is the complex conjugate operator. + +In order to keep the order of derivatives as low as possible, we use integration by parts on the Laplacian term: + +$$ +\begin{align} +\int_{\Omega}(\nabla^2 p) \bar v ~\mathrm{d}x = +-\int_{\Omega} \nabla p \cdot \nabla \bar v ~\mathrm{d}x ++ \int_{\partial \Omega} \frac{\partial p}{\partial n} \bar v ~\mathrm{d}s. +\end{align} +$$ + +Substituting in the original version and rearranging we get: + +$$ +\begin{align} +\int_{\Omega} \nabla p \cdot \nabla \bar v ~\mathrm{d}x +- k^2 \int_{\Omega} p \bar v ~\mathrm{d} x = \int_{\Omega} j \omega \rho_0 q \bar v ~\mathrm{d}x ++ \int_{\partial \Omega} \frac{\partial p}{\partial n} \bar v ~\mathrm{d}s. +\end{align} +$$ + +Since we are dealing with complex values, the inner product in the first equation is *sesquilinear*, +meaning it is linear in one argument and conjugate-linear in the other, +as explained in [The Poisson problem with complex numbers](../chapter1/complex_mode). + +The last term can be written using the Neumann and Robin BCs, that is: + +$$ +\begin{align} +\int_{\partial \Omega} \frac{\partial p}{\partial n} \bar v ~\mathrm{d}s = +-\int_{\partial \Omega_v} j \omega \rho_0 \bar v ~\mathrm{d}s +- \int_{\partial \Omega_Z} \frac{j \omega \rho_0 \bar{v}_n}{\bar{Z}} p \bar v ~\mathrm{d}s. +\end{align} +$$ + +Substituting, rearranging and taking out of integrals the terms with $j$ and $\omega$ we get the variational formulation of the Helmholtz. +Find $u \in V$ such that: + +$$ +\begin{align} +\int_{\Omega} \nabla p \cdot \nabla \bar v ~\mathrm{d}x ++ \frac{j \omega }{\bar{Z}} \int_{\partial \Omega_Z} \rho_0 p \bar v ~\mathrm{d}s +- k^2 \int_{\Omega} p \bar v ~\mathrm{d}x += j \omega \int_{\Omega} \rho_0 q \bar v ~\mathrm{d}x +-j \omega\int_{\partial \Omega_v} \rho_0 \bar{v}_n \bar v ~\mathrm{d}s \qquad \forall v \in \hat{V}. +\end{align} +$$ + +We define the sesquilinear form $a(p,v)$ is + +$$ +\begin{align} +a(p,v) = \int_{\Omega} \nabla p \cdot \nabla \bar v ~\mathrm{d}x ++ \frac{j \omega }{\bar{Z}} \int_{\partial \Omega_Z} \rho_0 p \bar v ~\mathrm{d}s, +- k^2 \int_{\Omega} p \bar v ~\mathrm{d}x +\end{align} +$$ + +and the linear form $L(v)$ reads + +$$ +\begin{align} +L(v) = j \omega \int_{\Omega}\rho_0 q \bar v ~\mathrm{d}x - j \omega \int_{\partial \Omega_v} \rho_0 \bar{v}_n \bar v ~\mathrm{d}s. +\end{align} +$$ diff --git a/chapter2/helmholtz_code.ipynb b/chapter2/helmholtz_code.ipynb new file mode 100644 index 00000000..f85949bb --- /dev/null +++ b/chapter2/helmholtz_code.ipynb @@ -0,0 +1,519 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Implementation\n", + "Author: Antonio Baiano Svizzero and Jørgen S. Dokken\n", + "\n", + "In this tutorial, you will learn how to:\n", + "- Define acoustic velocity and impedance boundary conditions\n", + "- Compute acoustic sound pressure for multiple frequencies\n", + "- Compute the Sound Pressure Level (SPL) at a given microphone position\n", + "\n", + "## Test problem\n", + "As an example, we will model a plane wave propagating in a tube.\n", + "While it is a basic test case, the code can be adapted to way more complex problems where velocity and impedance boundary conditions are needed.\n", + "We will apply a velocity boundary condition $v_n = 0.001$ to one end of the tube and an impedance $Z$ computed with the Delaney-Bazley model,\n", + "supposing that a layer of thickness $d = 0.02$ and flow resistivity $\\sigma = 1e4$ is placed at the second end of the tube.\n", + "The choice of such impedance (the one of a plane wave propagating in free field) will give, as a result, a solution with no reflections.\n", + "\n", + "First, we create the mesh with gmsh, also setting the physical group for velocity and impedance boundary conditions and the respective tags." + ] + }, + { + "cell_type": "code", + "execution_count": 1, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Info : Meshing 1D...\n", + "Info : [ 0%] Meshing curve 1 (Line)\n", + "Info : [ 10%] Meshing curve 2 (Line)\n", + "Info : [ 20%] Meshing curve 3 (Line)\n", + "Info : [ 30%] Meshing curve 4 (Line)\n", + "Info : [ 40%] Meshing curve 5 (Line)\n", + "Info : [ 50%] Meshing curve 6 (Line)\n", + "Info : [ 60%] Meshing curve 7 (Line)\n", + "Info : [ 60%] Meshing curve 8 (Line)\n", + "Info : [ 70%] Meshing curve 9 (Line)\n", + "Info : [ 80%] Meshing curve 10 (Line)\n", + "Info : [ 90%] Meshing curve 11 (Line)\n", + "Info : [100%] Meshing curve 12 (Line)\n", + "Info : Done meshing 1D (Wall 0.000926869s, CPU 0.000589s)\n", + "Info : Meshing 2D...\n", + "Info : [ 0%] Meshing surface 1 (Plane, Frontal-Delaunay)\n", + "Info : [ 20%] Meshing surface 2 (Plane, Frontal-Delaunay)\n", + "Info : [ 40%] Meshing surface 3 (Plane, Frontal-Delaunay)\n", + "Info : [ 60%] Meshing surface 4 (Plane, Frontal-Delaunay)\n", + "Info : [ 70%] Meshing surface 5 (Plane, Frontal-Delaunay)\n", + "Info : [ 90%] Meshing surface 6 (Plane, Frontal-Delaunay)\n", + "Info : Done meshing 2D (Wall 0.0311526s, CPU 0.030342s)\n", + "Info : Meshing 3D...\n", + "Info : 3D Meshing 1 volume with 1 connected component\n", + "Info : Tetrahedrizing 1285 nodes...\n", + "Info : Done tetrahedrizing 1293 nodes (Wall 0.0127086s, CPU 0.011753s)\n", + "Info : Reconstructing mesh...\n", + "Info : - Creating surface mesh\n", + "Info : - Identifying boundary edges\n", + "Info : - Recovering boundary\n", + "Info : Done reconstructing mesh (Wall 0.0264354s, CPU 0.02449s)\n", + "Info : Found volume 1\n", + "Info : It. 0 - 0 nodes created - worst tet radius 2.85527 (nodes removed 0 0)\n", + "Info : 3D refinement terminated (1746 nodes total):\n", + "Info : - 0 Delaunay cavities modified for star shapeness\n", + "Info : - 0 nodes could not be inserted\n", + "Info : - 6588 tetrahedra created in 0.0199287 sec. (330577 tets/s)\n", + "Info : 0 node relocations\n", + "Info : Done meshing 3D (Wall 0.0772706s, CPU 0.07553s)\n", + "Info : Optimizing mesh...\n", + "Info : Optimizing volume 1\n", + "Info : Optimization starts (volume = 0.01) with worst = 0.0138257 / average = 0.769909:\n", + "Info : 0.00 < quality < 0.10 : 18 elements\n", + "Info : 0.10 < quality < 0.20 : 69 elements\n", + "Info : 0.20 < quality < 0.30 : 66 elements\n", + "Info : 0.30 < quality < 0.40 : 67 elements\n", + "Info : 0.40 < quality < 0.50 : 150 elements\n", + "Info : 0.50 < quality < 0.60 : 273 elements\n", + "Info : 0.60 < quality < 0.70 : 913 elements\n", + "Info : 0.70 < quality < 0.80 : 1842 elements\n", + "Info : 0.80 < quality < 0.90 : 2090 elements\n", + "Info : 0.90 < quality < 1.00 : 1095 elements\n", + "Info : 152 edge swaps, 0 node relocations (volume = 0.01): worst = 0.300511 / average = 0.784453 (Wall 0.00193693s, CPU 0.00201s)\n", + "Info : No ill-shaped tets in the mesh :-)\n", + "Info : 0.00 < quality < 0.10 : 0 elements\n", + "Info : 0.10 < quality < 0.20 : 0 elements\n", + "Info : 0.20 < quality < 0.30 : 0 elements\n", + "Info : 0.30 < quality < 0.40 : 67 elements\n", + "Info : 0.40 < quality < 0.50 : 135 elements\n", + "Info : 0.50 < quality < 0.60 : 271 elements\n", + "Info : 0.60 < quality < 0.70 : 907 elements\n", + "Info : 0.70 < quality < 0.80 : 1872 elements\n", + "Info : 0.80 < quality < 0.90 : 2113 elements\n", + "Info : 0.90 < quality < 1.00 : 1081 elements\n", + "Info : Done optimizing mesh (Wall 0.00743052s, CPU 0.007641s)\n", + "Info : 1746 nodes 9265 elements\n" + ] + } + ], + "source": [ + "import gmsh\n", + "\n", + "gmsh.initialize()\n", + "\n", + "# meshsize settings\n", + "meshsize = 0.02\n", + "gmsh.option.setNumber(\"Mesh.MeshSizeMax\", meshsize)\n", + "gmsh.option.setNumber(\"Mesh.MeshSizeMax\", meshsize)\n", + "\n", + "\n", + "# create geometry\n", + "L = 1\n", + "W = 0.1\n", + "\n", + "gmsh.model.occ.addBox(0, 0, 0, L, W, W)\n", + "gmsh.model.occ.synchronize()\n", + "\n", + "# setup physical groups\n", + "v_bc_tag = 2\n", + "Z_bc_tag = 3\n", + "gmsh.model.addPhysicalGroup(3, [1], 1, \"air_volume\")\n", + "gmsh.model.addPhysicalGroup(2, [1], v_bc_tag, \"velocity_BC\")\n", + "gmsh.model.addPhysicalGroup(2, [2], Z_bc_tag, \"impedance\")\n", + "\n", + "# mesh generation\n", + "gmsh.model.mesh.generate(3)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Then we import the gmsh mesh with the ```dolfinx.io.gmshio``` function." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "from mpi4py import MPI\n", + "from dolfinx import fem, io, default_scalar_type, geometry\n", + "from dolfinx.fem.petsc import LinearProblem\n", + "import ufl\n", + "import numpy as np\n", + "import numpy.typing as npt\n", + "\n", + "mesh_data = io.gmshio.model_to_mesh(gmsh.model, MPI.COMM_WORLD, 0, gdim=3)\n", + "domain = mesh_data.mesh\n", + "assert mesh_data.facet_tags is not None\n", + "facet_tags = mesh_data.facet_tags" + ] + }, + { + "cell_type": "markdown", + "id": "92d8cd12", + "metadata": {}, + "source": [ + "We define the function space for our unknown $p$ and define the range of frequencies we want to solve the Helmholtz equation for." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "06cebc38", + "metadata": {}, + "outputs": [], + "source": [ + "\n", + "V = fem.functionspace(domain, (\"Lagrange\", 1))\n", + "\n", + "# Discrete frequency range\n", + "freq = np.arange(10, 1000, 5) # Hz\n", + "\n", + "# Air parameters\n", + "rho0 = 1.225 # kg/m^3\n", + "c = 340 # m/s\n" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "lines_to_next_cell": 2 + }, + "source": [ + "## Boundary conditions\n", + "\n", + "The Delaney-Bazley model is used to compute the characteristic impedance and wavenumber of the porous layer,\n", + "treated as an equivalent fluid with complex valued properties\n", + "\n", + "$$\n", + "\\begin{align}\n", + "Z_c(\\omega) &= \\rho_0 c_0 \\left[1 + 0.0571 X^{-0.754} - j 0.087 X^{-0.732}\\right],\\\\\n", + "k_c(\\omega) &= \\frac{\\omega}{c_0} \\left[1 + 0.0978 X^{-0.700} - j 0.189 X^{-0.595}\\right],\\\\\n", + "\\end{align}\n", + "$$\n", + "\n", + "where $X = \\frac{\\rho_0 f}{\\sigma}$.\n", + "\n", + "With these, we can compute the surface impedance, that in the case of a rigid passive absorber placed on a rigid wall is given by the formula\n", + "$$\n", + "\\begin{align}\n", + "Z_s = -j Z_c cot(k_c d).\n", + "\\end{align}\n", + "$$\n", + "\n", + "Let's create a function to compute it.\n", + "\n" + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "metadata": {}, + "outputs": [], + "source": [ + "# Impedance calculation\n", + "def delany_bazley_layer(f, rho0, c, sigma):\n", + " X = rho0 * f / sigma\n", + " Zc = rho0 * c * (1 + 0.0571 * X**-0.754 - 1j * 0.087 * X**-0.732)\n", + " kc = 2 * np.pi * f / c * (1 + 0.0978 * (X**-0.700) - 1j * 0.189 * (X**-0.595))\n", + " Z_s = -1j * Zc * (1 / np.tan(kc * d))\n", + " return Z_s\n", + "\n", + "\n", + "sigma = 1.5e4\n", + "d = 0.01\n", + "Z_s = delany_bazley_layer(freq, rho0, c, sigma)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Since we are going to compute a sound pressure spectrum, all the variables that depend on frequency (that are $\\omega$, $k$ and $Z$) need to be updated in the frequency loop.\n", + "To make this possible, we will initialize them as dolfinx constants.\n", + "Then, we define the value for the normal velocity on the first end of the tube" + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "metadata": {}, + "outputs": [], + "source": [ + "omega = fem.Constant(domain, default_scalar_type(0))\n", + "k = fem.Constant(domain, default_scalar_type(0))\n", + "Z = fem.Constant(domain, default_scalar_type(0))\n", + "v_n = 1e-5" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "We also need to specify the integration measure $ds$, by using ```ufl```, and its built in integration measures" + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "metadata": {}, + "outputs": [], + "source": [ + "ds = ufl.Measure(\"ds\", domain=domain, subdomain_data=facet_tags)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Variational Formulation\n", + "We can now write the variational formulation." + ] + }, + { + "cell_type": "code", + "execution_count": 6, + "metadata": {}, + "outputs": [], + "source": [ + "p = ufl.TrialFunction(V)\n", + "v = ufl.TestFunction(V)\n", + "\n", + "a = (\n", + " ufl.inner(ufl.grad(p), ufl.grad(v)) * ufl.dx\n", + " + 1j * rho0 * omega / Z * ufl.inner(p, v) * ds(Z_bc_tag)\n", + " - k**2 * ufl.inner(p, v) * ufl.dx\n", + ")\n", + "L = -1j * omega * rho0 * ufl.inner(v_n, v) * ds(v_bc_tag)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "The class ```LinearProblem``` is used to setup the PETSc backend and assemble the system vector and matrices.\n", + "The solution will be stored in a `dolfinx.fem.Function`, ```p_a```." + ] + }, + { + "cell_type": "code", + "execution_count": 7, + "metadata": { + "lines_to_end_of_cell_marker": 2 + }, + "outputs": [], + "source": [ + "p_a = fem.Function(V)\n", + "p_a.name = \"pressure\"\n", + "\n", + "problem = LinearProblem(\n", + " a,\n", + " L,\n", + " u=p_a,\n", + " petsc_options={\n", + " \"ksp_type\": \"preonly\",\n", + " \"pc_type\": \"lu\",\n", + " \"pc_factor_mat_solver_type\": \"mumps\",\n", + " },\n", + ")" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "lines_to_next_cell": 2 + }, + "source": [ + "## Computing the pressure at a given location\n", + "Before starting our frequency loop, we can build a function that, given a microphone position,\n", + "computes the sound pressure at its location.\n", + "We will use the a similar method as in [Deflection of a membrane](../chapter1/membrane_code).\n", + "However, as the domain doesn't deform in time, we cache the collision detection" + ] + }, + { + "cell_type": "code", + "execution_count": 8, + "metadata": {}, + "outputs": [], + "source": [ + "class MicrophonePressure:\n", + " def __init__(self, domain, microphone_position):\n", + " \"\"\"Initialize microphone(s).\n", + "\n", + " Args:\n", + " domain: The domain to insert microphones on\n", + " microphone_position: Position of the microphone(s).\n", + " Assumed to be ordered as ``(mic0_x, mic1_x, ..., mic0_y, mic1_y, ..., mic0_z, mic1_z, ...)``\n", + "\n", + " \"\"\"\n", + " self._domain = domain\n", + " self._position = np.asarray(\n", + " microphone_position, dtype=domain.geometry.x.dtype\n", + " ).reshape(3, -1)\n", + " self._local_cells, self._local_position = self.compute_local_microphones()\n", + "\n", + " def compute_local_microphones(\n", + " self,\n", + " ) -> tuple[npt.NDArray[np.int32], npt.NDArray[np.floating]]:\n", + " \"\"\"\n", + " Compute the local microphone positions for a distributed mesh\n", + "\n", + " Returns:\n", + " Two lists (local_cells, local_points) containing the local cell indices and the local points\n", + " \"\"\"\n", + " points = self._position.T\n", + " bb_tree = geometry.bb_tree(self._domain, self._domain.topology.dim)\n", + "\n", + " cells = []\n", + " points_on_proc = []\n", + "\n", + " cell_candidates = geometry.compute_collisions_points(bb_tree, points)\n", + " colliding_cells = geometry.compute_colliding_cells(\n", + " domain, cell_candidates, points\n", + " )\n", + "\n", + " for i, point in enumerate(points):\n", + " if len(colliding_cells.links(i)) > 0:\n", + " points_on_proc.append(point)\n", + " cells.append(colliding_cells.links(i)[0])\n", + "\n", + " return np.asarray(cells, dtype=np.int32), np.asarray(\n", + " points_on_proc, dtype=domain.geometry.x.dtype\n", + " )\n", + "\n", + " def listen(\n", + " self, recompute_collisions: bool = False\n", + " ) -> npt.NDArray[np.complexfloating]:\n", + " if recompute_collisions:\n", + " self._local_cells, self._local_position = self.compute_local_microphones()\n", + " if len(self._local_cells) > 0:\n", + " return p_a.eval(self._local_position, self._local_cells)\n", + " else:\n", + " return np.zeros(0, dtype=default_scalar_type)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "The pressure spectrum is initialized as a numpy array and the microphone location is assigned" + ] + }, + { + "cell_type": "code", + "execution_count": 9, + "metadata": {}, + "outputs": [], + "source": [ + "p_mic = np.zeros((len(freq), 1), dtype=complex)\n", + "\n", + "mic = np.array([0.5, 0.05, 0.05])\n", + "microphone = MicrophonePressure(domain, mic)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Frequency loop\n", + "\n", + "Finally, we can write the frequency loop, where we update the values of the frequency-dependent variables and solve the system for each frequency" + ] + }, + { + "cell_type": "code", + "execution_count": 10, + "metadata": {}, + "outputs": [], + "source": [ + "for nf in range(0, len(freq)):\n", + " k.value = 2 * np.pi * freq[nf] / c\n", + " omega.value = 2 * np.pi * freq[nf]\n", + " Z.value = Z_s[nf]\n", + "\n", + " problem.solve()\n", + " p_a.x.scatter_forward()\n", + "\n", + " p_f = microphone.listen()\n", + " p_f = domain.comm.gather(p_f, root=0)\n", + "\n", + " if domain.comm.rank == 0:\n", + " assert p_f is not None\n", + " p_mic[nf] = np.hstack(p_f)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## SPL spectrum\n", + "After the computation, the pressure spectrum at the prescribed location is available.\n", + "Such a spectrum is usually shown using the decibel (dB) scale to obtain the SPL, with the RMS pressure as input,\n", + "defined as $p_{rms} = \\frac{p}{\\sqrt{2}}$." + ] + }, + { + "cell_type": "code", + "execution_count": 11, + "metadata": {}, + "outputs": [ + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAB9QAAAKsCAYAAAC9C4YeAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjkuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8hTgPZAAAACXBIWXMAAA9hAAAPYQGoP6dpAAD6DElEQVR4nOzdd3hUddrG8Xtm0jsJpEFC772DgBVR7L279vKqq4iuK7u6rrrq6hZ7XRXsvWEvKChKk947pBfSeyYz8/4xyWQCARJSzpTv57q8rjP9Acmck3Of3/OYHA6HQwAAAAAAAAAAAAAAoAmz0QUAAAAAAAAAAAAAAOCJCNQBAAAAAAAAAAAAAGgGgToAAAAAAAAAAAAAAM0gUAcAAAAAAAAAAAAAoBkE6gAAAAAAAAAAAAAANINAHQAAAAAAAAAAAACAZhCoAwAAAAAAAAAAAADQDAJ1AAAAAAAAAAAAAACaQaAOAAAAAAAAAAAAAEAzCNQBAAAAAAAAAAAAAGiGoYF6WVmZZs2apZ49eyo0NFRHHXWUVqxY4Xrc4XDob3/7m5KSkhQaGqrp06dr+/btBlYMAAAAAAAAAAAAAPAXhgbq1157rb7//nu98cYbWr9+vWbMmKHp06crMzNTkvTYY4/pqaee0gsvvKBly5YpPDxcJ510kqqrq40sGwAAAAAAAAAAAADgB0wOh8NhxAdXVVUpMjJSn332mU499VTX/WPHjtXMmTP14IMPKjk5WXfccYfuvPNOSVJJSYkSEhI0b948XXTRRUaUDQAAAAAAAAAAAADwEwFGfXBdXZ1sNptCQkKa3B8aGqrFixdr9+7dysnJ0fTp012PRUdHa+LEiVqyZMlBA/WamhrV1NS4btvtdhUWFiouLk4mk6lj/jAAAAAAAAAAAAAAAK/gcDhUVlam5ORkmc2HbupuWKAeGRmpyZMn68EHH9TgwYOVkJCgd955R0uWLFG/fv2Uk5MjSUpISGjyuoSEBNdjzXnkkUd0//33d2jtAAAAAAAAAAAAAADvlp6erh49ehzyOYYF6pL0xhtv6Oqrr1b37t1lsVg0ZswYXXzxxVq5cuURv+ecOXM0e/Zs1+2SkhKlpqZq9+7dioyMbI+yAQAGslqt+umnn3TccccpMDDQ6HIAAB6EfQQA4GDYRwAADoX9BAD4n7KyMvXu3btF+bGhgXrfvn21aNEiVVRUqLS0VElJSbrwwgvVp08fJSYmSpJyc3OVlJTkek1ubq5GjRp10PcMDg5WcHDwAffHxsYqKiqq3f8MAIDOZbVaFRYWpri4OH7BAQA0wT4CAHAw7CMAAIfCfgIA/E/D931LRoYfuiF8JwkPD1dSUpKKior07bff6swzz1Tv3r2VmJioBQsWuJ5XWlqqZcuWafLkyQZWCwAAAAAAAAAAAADwB4auUP/222/lcDg0cOBA7dixQ3/60580aNAgXXXVVTKZTJo1a5b+8Y9/qH///urdu7fuvfdeJScn66yzzjKybAAAAAAAAAAAAACAHzA0UC8pKdGcOXOUkZGh2NhYnXvuuXrooYdcS+zvuusuVVRU6Prrr1dxcbGmTp2qb775RiEhIUaWDQAAAAAAAAAAAADwA4YG6hdccIEuuOCCgz5uMpn0wAMP6IEHHujEqgAAAAAAAAAAAAAAnsrhcKiurk42m+2gzwkMDJTFYmnzZxkaqAMAAAAAAAAAAAAA0FK1tbXKzs5WZWXlIZ9nMpnUo0cPRUREtOnzCNQBAAAAAAAAAAAAAB7Pbrdr9+7dslgsSk5OVlBQkEwm0wHPczgcys/PV0ZGhvr379+mleoE6gAAAAAAAAAAAAAAj1dbWyu73a6UlBSFhYUd8rndunXTnj17ZLVa2xSom4/4lQAAAAAAAAAAAAAAdDKz+fAxd3Mr14/os9rlXQAAAAAAAAAAAAAA8DEE6gAAAAAAAAAAAAAANINAHQAAAAAAAAAAAACAZhCoAwAAAAAAAAAAAADQDAJ1AAAAAAAAAAAAAIDXcDgc7fKcliBQBwAAAAAAAAAAAAB4vMDAQElSZWXlYZ9bW1srSbJYLG36zIA2vRoAAAAAAAAAAAAAgE5gsVgUExOjvLw8SVJYWJhMJtMBz7Pb7crPz1dYWJgCAtoWiROoAwAAAAAAAAAAAAC8QmJioiS5QvWDMZvNSk1NbTZwbw0CdQAAAAAAAAAAAACAVzCZTEpKSlJ8fLysVutBnxcUFCSzue0T0AnUAQAAAAAAAAAAAABexWKxtHk+eku0PZIHAAAAAAAAAAAAAMAHEagDAAAAAAAAAAAAANAMAnUAAAAAAAAAAAAAAJpBoA4AAAAAAAAAAAAAQDMI1AEAAAAAAAAAAAAAaAaBOgAAAAAAAAAAAAAAzSBQBwAAAAAAAAAAAACgGQTqAAAAAAAAAAAAAAA0g0AdAAAAAAAAAAAAAIBmEKgDAAAAAAAAAAAAANAMAnUAAAAAAAAAAAAAAJpBoA4AAAAAAAAAAAAAQDMI1AEAAAAAAAAAAAAAaAaBOgAAAAAAAAAAAAAAzSBQBwAAAAAAAAAAAACgGQTqAAAAAAAAAAAAAAA0g0AdAAAAAAAAAAAAAIBmEKgDAAAAAAAAAAAAANAMAnUAAAAAAAAAAAAAAJpBoA4AAAAAAAAAAAAAQDMI1AEAAAAAAAAAAAAAaAaBOgAAAAAAAAAAAAAAzSBQBwAAAAAAAAAAAACgGQTqAAAAAAAAAAAAAAA0g0AdAAAAAAAAAAAAAIBmEKgDAAAAAAAAAAAAANAMAnUAAAAAAAAAAAAAAJpBoA4AAAAAAAAAAAAAQDMI1AEAAAAAAAAAAAAAaAaBOgAAAAAAAAAAAAAAzSBQBwAAAAAAAAAAAACgGQTqAAAAAAAAAAAAAAA0g0AdAAAAAAAAAAAAAIBmEKgDAAAAAAAAAAAAANAMAnUAAAAAAAAAAAAAAJpBoA4AAAAAAAAAAAAAQDMI1AEAAAAAAAAAAAAAaAaBOgAAAAAAAAAAAAAAzSBQBwAAAAAAAAAAAACgGQTqAAAAAAAAAAAAAAA0g0AdAAAAAAAAAAAAAIBmEKgDAAAAAAAAAAAAANAMAnUAAAAAAAAAAAAAAJphaKBus9l07733qnfv3goNDVXfvn314IMPyuFwuJ7jcDj0t7/9TUlJSQoNDdX06dO1fft2A6sGAAAAAAAAAAAAAPgDQwP1Rx99VM8//7yeeeYZbd68WY8++qgee+wxPf30067nPPbYY3rqqaf0wgsvaNmyZQoPD9dJJ52k6upqAysHAAAAAAAAAAAAAPi6ACM//LffftOZZ56pU089VZLUq1cvvfPOO1q+fLkk5+r0J554Qvfcc4/OPPNMSdLrr7+uhIQEffrpp7rooosMqx0AAAAAAAAAAAAA4NsMDdSPOuoovfTSS9q2bZsGDBigtWvXavHixfrvf/8rSdq9e7dycnI0ffp012uio6M1ceJELVmypNlAvaamRjU1Na7bpaWlkiSr1Sqr1drBfyIAQEdr+C7nOx0AsD/2EQCAg2EfAQA4FPYTAOB/WvOdb2igfvfdd6u0tFSDBg2SxWKRzWbTQw89pEsvvVSSlJOTI0lKSEho8rqEhATXY/t75JFHdP/99x9w/3fffaewsLB2/hMAAIzy/fffG10CAMBDsY8AABwM+wgAwKGwnwAA/1FZWdni5xoaqL///vt666239Pbbb2vo0KFas2aNZs2apeTkZF1xxRVH9J5z5szR7NmzXbdLS0uVkpKiGTNmKCoqqr1KBwAYxGq16vvvv9eJJ56owMBAo8sBAHgQ9hEAgINhHwEAOBT2EwDgfxq6nLeEoYH6n/70J919992u1u3Dhw/X3r179cgjj+iKK65QYmKiJCk3N1dJSUmu1+Xm5mrUqFHNvmdwcLCCg4MPuD8wMJAdIQD4EL7XAQAHwz4CAHAw7CMAAIfCfgIA/Edrvu/NHVjHYVVWVspsblqCxWKR3W6XJPXu3VuJiYlasGCB6/HS0lItW7ZMkydP7tRaAQAAAAAAAAAAAAD+xdAV6qeffroeeughpaamaujQoVq9erX++9//6uqrr5YkmUwmzZo1S//4xz/Uv39/9e7dW/fee6+Sk5N11llnGVk6AAAAAAAAAAAAAMDHGRqoP/3007r33nt10003KS8vT8nJybrhhhv0t7/9zfWcu+66SxUVFbr++utVXFysqVOn6ptvvlFISIiBlQMAAAAAAAAAAAAAfJ2hgXpkZKSeeOIJPfHEEwd9jslk0gMPPKAHHnig8woDAAAAAAAAAAAAAPg9Q2eoAwAAAAAAAAAAAADgqQjUAQAAAAAAAAAAAABoBoE6AAAAAAAAAAAAAADNIFAHAAAAAAAAAAAAAKAZBOoAAAAAAAAAAAAAADSDQB0AAAAAAAAAAAAAgGYQqAMAAAAAAAAAAAAA0AwCdQAAAAAAAAAAAAAAmkGgDgAAAAAAAAAAAABAMwjUAQAAAAAAAAAAAABoBoE6AAAAAAAAAAAAAADNIFAHAAAAAAAAAAAAAKAZBOoAAAAAAAAAAAAAADSDQB0AAAA+weFw6NZ3VuvkJ37Wjrwyo8sBAAAAAAAA4AMI1AEAAOAT1mWUaP7aLG3JKdMbS/YaXQ4AAAAAAAAAH0CgDgAAAJ+QVVzl2t5dUGlgJQAAAAAAAAB8BYE6AAAAfEJeWY1rO72QQB0AAAAAAABA2xGoAwAAwCfklVW7tjOKKmWzOwysBgAAAAAAAIAvIFAHAACAT8grbVyhbrU5lFNafYhnAwAAAAAAAMDhEagDAADAJ+S6tXyXpDTmqAMAAAAAAABoIwJ1AAAA+IS8/VakM0cdAAAAAAAAQFsRqAMAAMAn5O+/Qp1AHQAAAAAAAEAbEagDAADA61ltdhVU1Da5j0AdAAAAAAAAQFsRqAMAAMDr7SuvOeC+vQTqAAAAAAAAANqIQB0AAABeL6/0wECdGeoAAAAAAAAA2opAHQAAAF4vr+zAQL2wolZl1VYDqgEAAAAAAADgKwjUAQAA4PXyyqpd24EWk2s7vbDKiHIAAAAAAAAA+AgCdQAAAHg995bvQ5OjXdtptH0HAAAAAAAA0AYE6gAAAPB67ivUx/Xs4tpmjjoAAAAAAACAtiBQBwAAgNdzX6E+rlesa5sV6gAAAAAAAADagkAdAAAAXi+vzBmom03SmNQY1/0E6gAAAAAAAADagkAdAAAAXq+h5XvXiGB1iwxWWJBFEi3fAQAAAAAAALQNgToAAAC8ms3u0L7yWklSfFSwTCaTUmPDJEkZRVWy2R1GlgcAAAAfUlNn098+26BHvtrMcSYAAICfIFAHAACAVyusqHWdzIyPDJEkpdQH6rU2u3JLqw2rDQAAAL7lo5WZen3JXr348y4t3JpndDkAAADoBATqAAAA8GoN7d4lKT4yWJJcK9Ql5qgDAACg/axJL3Jtb88rN7ASAAAAdBYCdQAAAHi1vLIa1zaBOgAAADrS5uwy13ZWcZWBlQAAAKCzEKgDAADAq+W5tXTvFuVs+e4eqKcTqAMAAKAd1Nns2ppLoA4AAOBvCNQBAADg1fJKD1yhnuIWqO8tIFAHAABA2+3aV6HaOrvrdkYRgToAAIA/IFAHAACAV3Nv+Z5Qv0K9R5dQ1320fAcAAEB72Jxd2uQ2K9QBAAD8A4E6AAAAvFpeWWPL94YV6iGBFiXWh+u0fAcAAEB72LRfoF5aXaeyaqtB1QAAAKCzEKgDAADAq7mvUO8aEezabpijXlBRq/Kauk6vCwAAAL5lc3bZAfdlFVc380wAAAD4EgJ1AAAAeLWGGeqx4UEKCmg8vHWfo84qdQAAALTV/i3fJdq+AwAA+AMCdQAAAHgth8Oh/PoV6g3t3hukugXqzFEHAABAW+wrr3Edd7rLJFAHAADweQTqAAAA8FolVVbV2uySpG77B+pxoa5tVqgDAACgLdxXp/ePj3BtE6gDAAD4PgJ1AAAAeC33+enxkSFNHkuNDXdts0IdAAAAbeEeqE8fkuDapuU7AACA7yNQBwAAgNfKLa12bcdH0fIdAAAAHWNzdplr+4RB8a5tAnUAAADfR6AOAAAAr5VX6r5CvWmg3jUiSKGBFkkE6gAAAGibhhXqAWaThveIVmx4kCQpq7j6UC8DAACADyBQBwAAgNdyb/meENW05bvJZHKtUs8orJLd7ujU2gAAAOAbaups2pFXLknq2y1CwQEWdY8JlSTllFarzmY3sjwAAAB0MAJ1AAAAeK28MreW7/utUJeklPpAvdZmV24Zq4cAAADQettzy1VXf3Hm4KRISVJyjPNiTpvdoVy3izwBAADgewjUAQAA4LXcV6jHR4Yc8HiTOeoFtH0HAABA6zW0e5ekwUlRkqTk+hXqEnPUAQAAfB2BOgAAALxWvvsM9agDV6inxjae6GSOOgAAAI7E5uwy13ZDoN7dLVDPLCJQBwAA8GUE6gAAAPBaDS3fI0MCFBJoOeDx1Di3FeoE6gAAADgCza1QbxKos0IdAADApxGoAwAAwCs5HA5Xy/fm5qdL+7V8J1AHAABAKzkcDm3OcQbqXSOC1a3+uJOW7wAAAP6DQB0AAABeqbymTpW1NknNz0+XpB5dCNQBAABw5HJKq1VcaZUkDUmOct3fvQsr1AEAAPwFgToAAAC8UsPqdKn5+emSFBJoUUL9Y+kE6gAAAGilpu3eI13bceFBCgpwnlplhToAAIBvI1AHAACAV8ordQvUD9LyXWps+76vvFYVNXUdXhcAAAB8x+bsMtf2kKTGFeomk8k1Rz2zqEoOh6PTawMAAEDnIFAHAACAV8orq3ZtJ0Q13/JdklLc5qinF7FKHQAAAC23qckK9agmjyXHOI9BK2ptKq3iwk0AAABfRaAOAAAAr5Tv1vK9WwtWqEtSWgGBOgAAAFquoeV7UIBZfbqGN3msYYW6xBx1AAAAX2ZooN6rVy+ZTKYD/rv55pslSdXV1br55psVFxeniIgInXvuucrNzTWyZAAAAHiIJjPUIw++Qr1JoM4cdQAAALRQVa1Ne/ZVSJIGJEQowNL0VGqyW6DOHHUAAADfZWigvmLFCmVnZ7v++/777yVJ559/viTp9ttv1+eff64PPvhAixYtUlZWls455xwjSwYAAICHyCttbPkeH3XwFeo949xavhOoAwAAoIW25pbJXj8afXBi1AGPJ7NCHQAAwC8EGPnh3bp1a3L7n//8p/r27atjjjlGJSUleuWVV/T222/r+OOPlyTNnTtXgwcP1tKlSzVp0iQjSgYAAICHaLpC/eCBegor1AEAAHAENmUdfH66JPVghToAAIBfMDRQd1dbW6s333xTs2fPlslk0sqVK2W1WjV9+nTXcwYNGqTU1FQtWbLkoIF6TU2NamoaT66WljoPfK1Wq6xWa8f+IQAAHa7hu5zvdAC59SvUQwPNCjY7Dvq9EBNsVkigWdVWu/YWVPL94cPYRwAADoZ9BI7Exsxi1/aA+LAD/v3ERwS6tjMKOc4EvBn7CQDwP635zveYQP3TTz9VcXGxrrzySklSTk6OgoKCFBMT0+R5CQkJysnJOej7PPLII7r//vsPuP+7775TWFhYM68AAHijhjEhAPxXVqFFkknhZpu+/vrrQz43JsCiHKtJaQXl+uLLr2Q2dU6NMAb7CADAwbCPQGv8ttl5vClJ6euXqmBz08etdqnh9OrGPdn66quMTq0PQPtjPwEA/qOysuWdLD0mUH/llVc0c+ZMJScnt+l95syZo9mzZ7tul5aWKiUlRTNmzFBU1IGtmQAA3sVqter777/XiSeeqMDAwMO/AIBPqrbaVLVkgSSpV2IXnXLKhEM+/7PC1crZmq86h0njph2vxKiQzigTnYx9BADgYNhHoLXsdof+supHSTYlR4fovDOObvZ5j25cqPzyWlWZQnTKKcd0bpEA2g37CQDwPw1dzlvCIwL1vXv36ocfftDHH3/sui8xMVG1tbUqLi5usko9NzdXiYmJB32v4OBgBQcfOEMzMDCQHSEA+BC+1wH/ll3a2JIpISr0sN8HPbuGS1vzXa9NiYvs0PpgLPYRAICDYR+BlkorqFRFjU2SNCQ56qD/bpK7hCm/vFZ55TVymCwKCjB3ZpkA2hn7CQDwH635vveII7y5c+cqPj5ep556quu+sWPHKjAwUAsWLHDdt3XrVqWlpWny5MlGlAkAAAAPkVdW7dqOjzrwYsr9pcY2jv7ZW1DRITUBAADAd2zKblyxNDjp4F0vu8c4Ox85HFJOSfVBnwcAAADvZfgKdbvdrrlz5+qKK65QQEBjOdHR0brmmms0e/ZsxcbGKioqSn/84x81efJkTZo0ycCKAQAAYLS8shrXdnzk4du3uwfq6YUtn48EAAAA/7S5xYF6qGs7s7hKqXFhB30uAAAAvJPhgfoPP/ygtLQ0XX311Qc89vjjj8tsNuvcc89VTU2NTjrpJD333HMGVAkAAABPklfqtkI9snUr1NMI1AEAAHAYLQ3Uk90C9aziqg6tCQAAAMYwPFCfMWOGHA5Hs4+FhITo2Wef1bPPPtvJVQEAAMCTNVmh3oKW7z26EKgDAACg5TbnOAP1sCCLesYefNU5gToAAIDv84gZ6gAAAEBrtLble2iQxbWSPa2QE50AAAA4uLJqq9LrjxkHJkbKbDYd9Ln7t3wHAACA7yFQBwAAgNfJbWXLd6mx7fu+8hpV1tZ1SF0AAADwfltyylzbh2r3LhGoAwAA+AMCdQAAAHid/PoV6kEWs2LCAlv0Gvc56umsUgcAAMBBbMpq2fx0SYoJC1RooEUSLd8BAAB8FYE6AAAAvE5Dy/dukcEymQ7egtNdSixz1AEAAHB4m7MbA/UhhwnUTSaTkmOcI4gyi6vkcDg6tDYAAAB0PgJ1AAAAeJXaOrsKK2olOQP1luoZR6AOAACAw2sI1E0maVBi5GGf372L8ziz2mpXUaW1Q2sDAABA5yNQBwAAgFfZV17j2k6Ianmg3rTlO4E6AAAADmSzO7Q11zlDvWdsmMKDAw77mu71K9Ql2r4DAAD4IgJ1AAAAeJWGdu+SFB8ZcohnNpVKy3cAAAAcxu59Faq22iUdfn56g+ToUNd2RhGBOgAAgK8hUAcAAIBXySutdm3Ht6Lle7fIYAUHOA9/CdQBAADQHPf56S0N1Lt3aQzUWaEOAADgewjUAQAA4FWarFBvRct3k8nkWqWeXlgpu93R7rUBAADAux1JoJ4cQ6AOAADgywjUAQAA4FWOtOW71Nj2vabOrny3WewAAACAtH+gHtmi13R3C9QzCdQBAAB8DoE6AAAAvIp7y/durWj5LkkpbnPU9xbQ9h0AAABNbc4ukyRFhQQ0CcoPJTE6RCaTc5sV6gAAAL6HQB0AAABe5UhbvkuNK9Ql5qgDAACgqaKKWuXUX7w5KClKpoaU/DACLWYl1HdOyiyuPsyzAQAA4G0I1AEAAOBV8sqcJynNJikunEAdAAAA7cO93fuQFs5Pb5Ac4wzU95XXqNpqa9e6AAAAYCwCdQAAAHiVvFLnCvWuEcGymFu2aqhBalxjoJ5OoA4AAAA3m9oQqHfv0nicmV3CKnUAAABfQqAOAAAAr2GzO7Sv3BmoJ0SFtPr1KV1YoQ4AAIDmuQfqg49whbrEHHUAAABfQ6AOAAAAr1FQUSO7w7kdH9m6du+SFBpkUbf61xGoAwAAwN3m7DJJksVsUv+EiFa9tntMqGs7s4hAHQAAwJcQqAMAAMBrNLR7l6T4qNYH6lLjHPX8shpV1TLfEgAAAFJtnV078pyBep+u4QoJtLTq9cnRboE6K9QBAAB8CoE6AAAAvEZ+WWOg3i2y9S3fpcZAXZLSi1ilDgAAAGlnfrmsNmcrpNa2e5ek7l0aA3VavgMAAPgWAnUAAAB4jdzSatf2kbR8l6QUt0A9rYBAHQAAANKGzBLX9pDk1gfqyW4t37NKCNQBAAB8CYE6AAAAvEae2wr1Iw3Ue8U1Buo78svbXBMAAAC8n3ugPqJ7dKtfHxUSoIjgAEnMUAcAAPA1BOoAAADwGnllbivUo46s5bv7iqP1bidOAQAA4L/WuR0XDj2CQN1kMql7/Sr1rJJq2e2OdqsNAAAAxiJQBwAAgNfIK237CvV+3SIUEug8DF6fQaAOAADg76w2uzZllUpydjOKDg08ovdJjnFe8FlbZ1dBRW271QcAAABjEagDAADAa7i3fO8acWSBeoDFrCFJzlXqaYWVKqm0tkttAAAA8E7bc8tVU2eXJA3vEXPE7+M+Rz2zmLbvAAAAvoJAHQAAAF4jvz5QjwsPUlDAkR/KDndr47khi1XqAAAA/qyt89MbdO/SGKhnEagDAAD4DAJ1AAAAeAWHw+EK1LsdYbv3Bu4rj9bR9h0AAMCvrcssdm0P79GGQD2GQB0AAMAXEagDAADAKxRXWlVrc7bijI8KadN7NVmhnkmgDgAA4M/Wu11gOTQ56ojfx73le0YRgToAAICvIFAHAACAV3Cfnx7fxhXqfbuFKyTQeSi8nkAdAADAb9XW2bU5u0yS1KdbuCJDAo/4vVihDgAA4JsI1AEAAOAVckurXdttDdQDLGYNTXauUk8rrFRxZW2b3g8AAADeaVtumasLUlvmp0vOY1SL2SRJyiohUAcAAPAVBOoAAADwCu25Ql3av+17aZvfDwAAAN7HvVvR8B4xbXqvAItZifWjiTJp+Q4AAOAzCNQBAADgFfLK3Faot3GGuiQNcwvU12UWt/n9AAAA4H3Wuc1PH97GFeqSlBzjPE4tqrSqsrauze8HAAAA4xGoAwAAwCvklbbvCvURPdxXqDNHHQAAwB+tr7+w0mSShiZHtfn9ms5Rrz7EMwEAAOAtCNQBAADgFfKbtHxv+wr1vt0iFBpokdS01ScAAAD8Q02dTVtzyiRJ/bpFKDw4oM3vmewWqGcW0/YdAADAFxCoAwAAwCs0bfne9hXqFrNJQ+pXIaUXVqmoorbN7wkAAADvsTWnTFabQ5I0vEfb271LTQP1LAJ1AAAAn0CgDgAAAK+QV79CPSokQCH1K8vbyn1O5oYsVqkDAAD4E/f56SPaYX66JHXvQqAOAADgawjUAQAA4PEcDodrhnp8VNvbvTdwD9Rp+w4AAOBf1rsF6u21Qt19hnpmEYE6AACALyBQBwAAgMcrq6lTldUmSYqPbHu79wYj3E6cup9QBQAAgO9bV39BpdkkDUlq/5bvzFAHAADwDQTqAAAA8HgNq9Ol9g3U+3SLUFiQs308K9QBAAD8R7XVpu25ZZKkAQmRCg1qn5FCEcEBig4NlCRllRCoAwAA+AICdQAAAHi8vLJq13Z7tny3mE0akhQlScooqlJRRW27vTcAAAA81+bsUtXZHZKajgFqDw2r1LOLq2Wr/wwAAAB4LwJ1AAAAeLz8so5ZoS41nZfJKnUAAAD/4H7c117z0xt0j3FeAFpndzQ5jgUAAIB3IlAHAACAx3Nv+d6tvQP17gTqAAAA/mZdhlug3s4r1LszRx0AAMCnEKgDAADA4zVp+R7Zfi3fpf0C9QwCdQAAAH+wof5CygCzSYPrRwC1l2S3QD2LQB0AAMDrEagDAADA4+W5tcpMiGrfFep9ukUoLMgiiRXqAAAA/qCq1qZtuWWSpAEJkQoJtLTr+yezQh0AAMCnEKgDAADA47m3fI+Pat8V6hazSUOTnauSMourVFhR267vDwAAAM+yKbtEdodzu73bvUtS9y6sUAcAAPAlBOoAAADweA0t38OCLIoIDmj39x/GHHUAAAC/0WR+eo8OCNRp+Q4AAOBTCNQBAADg8RpWqMdHtm+79wYj3E6kbiBQBwAA8GnuF1CO6IBAvVtEsAItJklSRhGBOgAAgLcjUAcAAIBHK6m0qqymTpIUH9m+7d4buLf6XJ9BoA4AAODLGo73Ai0mDUyMbPf3N5tNSox2HreyQh0AAMD7EagDAADAo/26c59re1gHzLiUpN5dIxQeZJFEy3cAAABfVlFTpx355ZKkQYlRCg6wdMjnNLR9L62uU1m1tUM+AwAAAJ2DQB0AAAAe7edt+a7towd07ZDPsJhNGprsDOszi6tUWFHbIZ8DAAAAY23MKpXD4dzuqIs1JSnZbY76nn2VHfY5AAAA6HgE6gAAAPBYDofDFagHBZg1sXdch32W+wlVVqkDAAD4po6en+56b7djy8U79h3imQAAAPB0BOoAAADwWDvyypVVUi1Jmtg7VqFBHdOSU2p6QnV9RnGHfQ4AAACM436cN7wDV6gfOzDetb1wa16HfQ4AAAA6HoE6AAAAPNYi93bv/bt16GexQh0AAMD3ras/zgsKMGtAQmSHfU6vruHqFRcmSVq5t4g56gAAAF6MQB0AAAAe6+ftje0xjxnYsYF6n67hCq9fAb8+g0AdAADA15RVW7Urv0KSNDgxUkEBHXtq9JgBzuPXOrtDv+4o6NDPAgAAQMchUAcAAIBHqrbatGyX88RjYlSI+sdHdOjnmc0mDa1fpZ5VUq2C8poO/TwAAAB0ro1Zpa7t4R04P72Be9v3Rdto+w4AAOCtCNQBAADgkZbvLlRNnV2SdPSArjKZTB3+mcNp+w4AAOCz3LsQjege0+GfN6lPnGsV/MKt+XI4HB3+mQAAAGh/BOoAAADwSD+7z08f0LHt3hs0CdRp+w4AAOBT1rldMNkZK9RDgyya1CdOkpRdUq1tueUd/pkAAABofwTqAAAA8EiL6gN1s0ma2q9rp3ym+4lVVqgDAAD4lvUZxZKk4ABzh48TanCM24WhtH0HAADwTgTqAAAA8DhZxVXanudcwTMyJUYxYUGd8rm948IVERwgSdpAoA4AAOAzSqqs2lNQKUkakhylAEvnnBY9dmBjoL5wa/4hngkAAABPZXignpmZqcsuu0xxcXEKDQ3V8OHD9fvvv7sedzgc+tvf/qakpCSFhoZq+vTp2r59u4EVAwAAoKP9st2t3Xv/zmn3Lklms0lDk6MkSVkl1dpXXtNpnw0AAICOszHTfX56x7d7b9Cna7hSYkMlSSv2FKq8pq7TPhsAAADtw9BAvaioSFOmTFFgYKC+/vprbdq0Sf/5z3/UpUsX13Mee+wxPfXUU3rhhRe0bNkyhYeH66STTlJ1dbWBlQMAAKAj/bxtn2u7s+anN2gyR51V6gAAAD6h6fz0mE77XJPJpGMHxEuSrDaHluws6LTPBgAAQPswNFB/9NFHlZKSorlz52rChAnq3bu3ZsyYob59+0pyrk5/4okndM899+jMM8/UiBEj9PrrrysrK0uffvqpkaUDAACgg9jsDi3e4QzUo0ICNLJH560gkprOUd+QQaAOAADgC9a7HdeN6OTjS/c56gu3MkcdAADA2wQY+eHz58/XSSedpPPPP1+LFi1S9+7dddNNN+m6666TJO3evVs5OTmaPn266zXR0dGaOHGilixZoosuuuiA96ypqVFNTWNrztLSUkmS1WqV1Wrt4D8RAKCjNXyX850O+K7V6cUqqXL+jB/VN04Ou01Wu63TPn9wQrhre216Ed83XoR9BADgYNhHYG1GsSQpNNCs1JjgTv23ML5nlAItJlltDi3cmqfa2lqZTKZO+3wAh8d+AgD8T2u+8w0N1Hft2qXnn39es2fP1l/+8hetWLFCt956q4KCgnTFFVcoJydHkpSQkNDkdQkJCa7H9vfII4/o/vvvP+D+7777TmFhYe3/hwAAGOL77783ugQAHeTrdJMkiyQppipLX32V2amfb3dIwRaLamwm/b4rT1999VWnfj7ajn0EAOBg2Ef4pwqrlFHkPA2aGGLTt9983ek19I4wa1uJWZnF1Zr70ddK5DQl4JHYTwCA/6isrGzxcw0N1O12u8aNG6eHH35YkjR69Ght2LBBL7zwgq644oojes85c+Zo9uzZrtulpaVKSUnRjBkzFBUV1S51AwCMY7Va9f333+vEE09UYGCg0eUA6ABzX1omydmS86ZzjlNSdEin1/BOzgot31Ok4lqTJhx9grpGBHd6DWg99hEAgINhH+HfFu8okH5fKUk6elhPnXLKoE6vITt6j/75zTZJkil5qE45qmen1wDg4NhPAID/aehy3hKGBupJSUkaMmRIk/sGDx6sjz76SJKUmJgoScrNzVVSUpLrObm5uRo1alSz7xkcHKzg4ANPeAYGBrIjBAAfwvc64JtKKq1aVz/fsn98hFK7RhpSx8iUGC3fUyRJ2pJbqeO6RBhSB44M+wgAwMGwj/BPm3LKXdujUrsY8m/ghMGJrkD9lx0Fuv6Yfp1eA4DDYz8BAP6jNd/35g6s47CmTJmirVu3Nrlv27Zt6tnTeYVm7969lZiYqAULFrgeLy0t1bJlyzR58uROrRUAAAAdb/GOfbI7nNtHD+hmWB3Duke7ttdnlhhWBwAAANpufUbj8dxwt+O8ztQvPkLdY0IlSct2Faqyts6QOgAAANB6hgbqt99+u5YuXaqHH35YO3bs0Ntvv62XXnpJN998syTJZDJp1qxZ+sc//qH58+dr/fr1+sMf/qDk5GSdddZZRpYOAACADvDztnzXtpGBuvuJ1nUZBOoAAADerOECyfAgi3p3NabzkMlk0jEDnce3tTa7lu4qMKQOAAAAtJ6hgfr48eP1ySef6J133tGwYcP04IMP6oknntCll17qes5dd92lP/7xj7r++us1fvx4lZeX65tvvlFISOfP0gQAAEDHcTgc+nm7M1APDjBrYu9Yw2rpFReuyGDndKQNrFAHAADwWpnFVcosrpIkDe0eLYvZZFgtx7hdMLpwa/4hngkAAABPYugMdUk67bTTdNpppx30cZPJpAceeEAPPPBAJ1YFAACAzrYjr1zZJdWSpAm9YxUSaDGsFrPZpKHdo7R0V6FySquVV1at+Egu6AQAAPA2n67OdG1P7dfVwEqkKf26KtBiktXm0MKt+XI4HDKZjAv4AQAA0DKGrlAHAAAAGixya/d+jIHt3huM7BHj2mYFEQAAgPdxOBz6cGWG6/bZo7sbWI0UERygcT2dXZjSCiu1p6DS0HoAAADQMgTqAAAA8AieFqifNCzRtf2R24lYAAAAeIdVacXava9CkjS5T5xSYsMMrkg6dqB72/c8AysBAABASxGoAwAAwHDVVpuW7y6UJCVFh6hffITBFUmjU2LUp2u4JGnZ7kKlF7KCCAAAwJu4r04/d2wPAytpdMxA5qgDAAB4GwJ1AAAAGG7Z7kLV1NklSUf37+YRsyRNJlOTE68fr8o8xLMBAADgSaqtNn2xLkuSFBZk0Uy37kNGGpgQqcSoEEnS0l0FqrbaDK4IAAAAh0OgDgAAAMP97Nbu/WgPaPfe4OzR3dWQ7X+8OkMOh8PYggAAANAi323KVVl1nSRp5rAkhQcHGFyRk8lkcrV9r6mza+muAoMrAgAAwOEQqAMAAMBwDYG62SRN7dfV4GoaJceE6qi+cZKkvQWV+n1vkcEVAQAAoCU+cmv3fp6HtHtvcCxt3wEAALyK3wTqW3NKjS4BAAAAzcgqrtL2vHJJ0qiUGEWHBRpcUVPnjmk8Aet+YhYAAACeKbe0Wr9sdwbV3WNCNbF3rMEVNXVUv64KMDvbIC3aRqAOAADg6fwmUH9vRbrRJQAAAKAZntruvcHJwxIVHmSRJH25Lps5lwAAAB7uk9WZstdP6jl3bA+Z68NrTxEVEqgxPbtIknbvq9DeggqDKwIAAMCh+E2g/sW6LJXX1BldBgAAAPbz83bPDtTDggI0c3iSJKmspk7fbswxuCIAAAAcjMPh0IduXYXOHdPdwGoOzr3tO6vUAQAAPJvfBOqVtXZ9ujrT6DIAAADgps5m1+Lt+yRJ0aGBGtkjxtiCDsK97fvHqzimBAAA8FTrMkq0o36c0IReseoZF25wRc07ZgBz1AEAALyF3wTqkvTm0r1yOBxGlwEAAIB6azNKVFrt7CI0tV9XWTysHWeDib1j1T0mVJL0y/Z85ZZWG1wRALSPqlqbvtuYo6KKWqNLAYB20WR1+ljPXJ0uSUOSohQfGSxJWrKzgLFCAAAAHsyvAvUtOWValVZkdBkAAACo13R+elcDKzk0s9nkahdqd4jORwB8xl8/Wa/r31ipa15bYXQpANBmNXU2zV+bJUkKCTTrlPqxPZ7IZDK5VqlXWW1asafQ4IoAAABwMH4VqEvSm0vTjC4BAAAA9X7amufa9sT56e7OcWv7/tGqDDofAfAJq9OLJUmr0opVU8fqSADebcHmPJVUWSVJJw9NVGRIoMEVHdqxA+Nd27R9BwAA8Fx+E6hHhwZIkr5cl61CWtkBAAAYbnVakdZllEiSBidFKSk61OCKDq1X13CN69lFkrQtt1wbMksNrggA2q68ps61nVdaY2AlANB2HzVp997jEM/0DFP7dVXDxKOftuZxwSYAAICH8ptA/azRzhadtTa7Pvg93eBqAAAA8PIvu13bVx7V08BKWs79xOxHqzIO8UwA8A4VboF6dkm1gZUAQNvkl9VoYf04oaToEB3V13PHCTWIDgvU2PoLNnflV2jlXkZVAgAAeCK/CdTPH5fi2n57eZrsdq74BAAAMEp6YaW+3pAtSeoaEaQzR3U3uKKWOXVEkoIDnIfQn63JVG2d3eCKAODI2e0OVdY2tnnPLqkysBoAaJvP1mTKVn++7+zR3WVpWPrt4S4an+ra/t8vuwysBAAAAAfjN4F6r7hwTe3nvDJ1b0GlFu/YZ3BFAAAA/uvVX3er4frGP0zupZBAi7EFtVBUSKBmDE2UJBVVWpvMgAcAb1NRW9fkdg4r1AF4KYfDoQ+9rN17g9NHJishKliS9N2mXO3ZV2FwRQAAANif3wTqknTZpMYrPt9cutfASgAAAPxXSZVV769wjuAJCTTrskne0e69wbljGlfTu8/pBABvU1Fja3Kblu8AvNXGrFJtySmTJI1OjVHfbhEGV9RyQQFmXXFUL0mSwyG9snj3oV8AAACATudXgfr0wQmuKz5/2JxLOzsAAAADvLM8TRX1LYbPHdNDseFBBlfUOtP6d1N8pPOY8qeteSqsqDW4IgA4MuU1rFAH4BvcV6ef50Wr0xtcOqGnwoKcHZs+WJmuIo4vAQAAPIpfBeoBFrNrLpHdIb2zPN3gigAAAPxLbZ1d837dI0kymaRrpvY2tqAjYDGbdPZo5yp1q82h+WsyDa4IAI5Mxf6BeimBOgDvU1tn1/y1WZKcq71PG5FscEWtFx0WqAvGpUiSqq12vbWMzpoAAACexK8CdUm6aEKKLGaTJOnd5Wmy2uwGVwQAAOA/vlyf5Qpspg9OUB8vasfpzn0u50erCNQBeKcDAnVWqAPwQu4dg2YMSVB0aKDBFR2Za6b2Vv0pS837ba9q6myHfgEAAAA6jd8F6knRoTphULwkKa+sRgs25xpcEQAAgH9wOBz638+NMyGvm9bHwGraZkBCpIZ3j5Ykrc8s0bbcMoMrAoDWK9svUM8rq1YdF50D8DIfubV7P9cL2703SIkN08nDEiVJ+8pr9NnqLIMrAgAAQAO/C9Ql6bJJPV3bby5NM7ASAAAA/7FkZ4E2ZZdKkkb2iNb4Xl0Mrqhtzh3T3bXtfiIXALzF/ivU7Q4pv7zGoGoAoPUKymv045Y8SVJ8ZLCm9etqcEVtc63bBacvL94lh8NhYDUAAABo4JeB+tR+XdUzLkyStHjHPu3eV2FwRQAAAL7vf7/scm1fO62PTCaTgdW03RmjuivQ4vwzfLI6k1WdALzO/oG6JGXT9h2AF5m/Nkt1dmfofPbo7gqwePepzjGpXTSup/Oi02255Vq0Ld/gigAAACD5aaBuNpt06cRU1+23lu41sBoAAADftyOvTD9tdZ4Q7B4Tqpn17Sy9WWx4kI4b2DhKaPGOfQZXBACtU15z4Hxe5qgD8BZ2u0PvrUh33fbmdu/umqxS/2X3IZ4JAACAzuKXgboknTc2RUEBzj/+ByszVG098EQCAAAA2of7ycCrpvTy+tVDDdxP3H60KtPASgCg9VihDsCbfbE+W1tyyiRJI1NiNCAh0uCK2seJQxKadNbcmFVicEUAAADwjTOZRyA2PEinDk+SJJVUWfXFumyDKwIAAPBN+WU1+ni1M2yODA7QheNTDK6o/Rw3MF5dwgIlSd9uzFFeKUEUAO9R3kygnlNSZUAlANA6NXU2/evbLa7bs08cYGA17ctiNumaqb1dt19hlToAAIDh/DZQl6TLJrm1fV9G23cAAICO8MbSvaqtc84Xv3hiqiJDAg2uqP0EBZh1Xv0q9do6u55csN3gigCg5VihDsBbvbU0TemFzguApvSL09H9uxpcUfs6b2wPRYc6j5nnr81iHAcAAIDB/DpQH5PaRYMSne2gVqcV00IJAACgnVXV2vTmUueFiwFmk648qpexBXWAG4/pq4jgAEnSuyvStTO/3OCKAKBlKmqbW6FOaAPAs5VWW/X0j40XMc6ZOVgmk8nAitpfWFCALp/UU5JUZ3do3m97jC0IAADAz/l1oG4ymXRZ/cGpJL25NM3AagAAAHzPR6syVFhRK0k6dUSSkmNCDa6o/cVFBOuGo/tIkmx2h/797VaDKwKAlimvsbm2gyzO0wM5jK4A4OFeWLhTRZVWSdKZo5I1rHu0wRV1jD8c1dP13fz2sr3NjukAAABA5/DrQF2SzhrdXeFBFknSZ2syVVBeY3BFAAAAvsFud+jVxY0zH6+b1sfAajrWNdN6q2tEsCTp6w05WpVWZHBFAHB47i3fe3UNkyTlllbLbncYVRIAHFJ2SZVeqT++DLKYdeeMgQZX1HHiI0N05qhkSVJpdZ3eX5FucEUAAAD+y+8D9YjgAJ1bP/eystam/3y/zeCKAAAAfMOCLXnata9CkjSpT6zPrh6SnG05Z03v77r9z6+3yOEgkALg2cqrnYF6SKBZPbo4A3WrzaGC+s4iAOBpHv9+m2rq7JKkP0zuqZTYMIMr6ljXul2Q+uqvu1VnsxtYDQAAgP/y+0Bdkm45rp9rlfo7y9OYpQ4AANAO/vfLLte2L69Ob3Dh+BT17houSVq+u1ALt+YbXBEAHFpD++CI4AAlRoe47meOOgBPtDWnTB+uzJAkRYYE6Obj+hlcUccbmBipowd0kyRlFFXp2425BlcEAADgnwjUJcVHheiPJzhXFDkc0gOfb2JFEQAAQBusyyjW8t2FkqQ+3cJ13MB4gyvqeIEWs/50UmPb0Ue/2SIbbZMBeLCKWmegHh4coKSoxkA9u6TKqJIA4KAe/WaLGg6tbj6un7qEBxlbUCe53u3C1Jd+2cU5SwAAAAMQqNe7akov9YxztolatrtQX63PMbgiAAAA7/XUgu2u7Wun9pHZbDKwms4zc1iiRqbESJK25JTp09WZxhYEAIfQMEM9PGi/FeqlrFAH4FmW7CzQj1vyJElJ0SG68qhexhbUiab0i9OgxEhJ0tr0Yv2+t8jgigAAAPwPgXq94ACL7jl1iOv2w19tVrXVZmBFAAAA3um7jTn6YbPzhGe3yGCdM6a7wRV1HpPJpLtPHuS6/d/vt3FMCcAj1dTZZLU5VzlGBAcoKTrU9Vg2Ld8BeBC73aFHvt7sun3HjIEKCbQYWFHnMplMTcYnPb9wp4HVAAAA+CcCdTfTB8drWv+ukqTM4iq99POuw7wCAAAA7ipq6vT3+Rtdt+89bYhfnfCUpMl943TsQOesy8ziKr25dK/BFQHAgSpqGi/2CQ+2MEMdgMf6cn221mWUSJIGJUbq7NH+c7Fmg9NHJiuxfjTHj1vytGAzs9QBAAA6E4G6G5PJpL+dNkSW+pakzy3coaxiZscBAAC01OPfb1NWfRAzrX9XnT4iyeCKjHHXSYNkqu9y/8xPO1RSZTW2IADYT0O7d8k5Q909UGeGOgBPUVtn17++3eq6fffMQa7zdv4kKMCsOac0dkG6b/5GVdXSBQkAAKCzEKjvp39CpC6f1FOSVG21659fbzG4IgAAAO+wIbNEr/66W5IUHGDWP84aJpPJ/054StKQ5CidPcq5eqq40qoXF9GaE4BnKXcL1COCAxQRHKDI4ABJUm5pjVFlAUATby3bq7TCSknOWeLHDOhmcEXGOWNksqb0i5MkZRRV6ekftxtcEQC03Psr0nXl3OXakFlidCkAcEQI1Jtx+/QB6hIWKEmavzZLK/YUGlwRAACAZ7PZHfrrJ+tld47j1R+P76eeceHGFmWw208coCCL83D71V9300IZgEfZf4W6JNcq9eySKjkcDkPqAoAGpdVWPf3jDtftu08e7LcXa0rOzpoPnDnMdXz50s+7tC23zOCqAODwKmvrdM9nG7Rwa74e/YYFjAC8E4F6M6LDAnXHjIGu2/d/vlE2OycTAAAADuatZXu1tn62Zb/4CF1/dF+DKzJeSmyYLp/c2PnoyQXbDK4IABrtv0JdagzUq612RlUAMNyLi3aqsKJWknTmqGQN7xFtcEXG69stQjce6zzOrrM7dM+nG7gACoDHSy+sUm2dXZK0Nr2Y7y0AXolA/SAunpCqQYmRkqQNmaX6cGW6wRUBAAB4ptzSav3rm8bZlg+fPVxBARxmStLNx/VztVB+b0W6duSVG1wRADg1F6gnNZmjTlcNAMbJKanWK4udo4SCLGbd6bbwxd/ddGxf9YwLkyQt312oj1ZlGlwRABxaw+gOSSqtrlNGUZWB1QDAkeFM50FYzCbdd/pQ1+1/fbtVpdVcoQ8AALC/B77YpLL6YOaCcT00oXeswRV5jtjwINcqIrtDeoz2dgA8RPMt30Nd9zGmAoCR/vn1ZlVbnasZL5/cUymxYQZX5DlCAi164MxhrtsPf7VZRfUr+QHAE+0tqGhye2NWqUGVAEBTVbW2Fj+XQP0QJveN0ynDEyVJ+8pr9Yzb3CYAAABIP23N05frsiU5w+M5MwcbXJHnuWpKL8VHBkuSvtuUq5V7Cw2uCACk8prGEwfhwRZJrFAH4Bk+W5OpT9dkSZIiQwJ0y3H9DK7I8xwzoJtOHZEkSSqsqNVj33LRJgDPle62Ql2SNmUTqAPwDP/4clOLn0ugfhhzZg52tSyd++tu7cqnTScAAIDkvIrz3k83uG7/5ZTB6hIeZGBFniksKECzpg9w3X7wi82qs9kNrAgAmq5Q33+GuiTllNCKE0DnSy+s1D2fNB5fPnjmMI4vD+Jvpw1xfX+/szydizYBeKy0/QN1VqgD8AAf/J6uz+ov4mwJAvXDSIkN0w1H95EkWW0O/ePLzQZXBAAA4Bme+nG7a/bZpD6xOndMd4Mr8lwXjOuhPt3CJUlr0ov13MKdBlcEwN811/KdFeoAjFRns+u2d1e7RgmdM7q7zhrN8eXBJESF6M4ZjRdt/vWTDbJy0SYAD3RgoF5iUCUA4LQ1p0z3frbh8E90Q6DeAv93bF8lRjlPLPy4JU8/bc0zuCIAAABjbc0p0/9+3iVJCrKY9dDZw2UymQyuynMFWMz613kjZK7/K3pywXatTisytigAfq28mRXqSVFuM9RLCdQBdK6nFmzXqrRiSVJqbJjuP3OosQV5gcsn99Kw7lGSpC05ZZr36x5jCwKA/djtDqUXNe18lFVSraKKWoMqAuDvKmrqdNNbK1Vtbd2FiATqLRAWFKA5pwxy3f77/I1NTj4AAAD4E7vdob9+sl51dock6cZj+6pvtwiDq/J8Y3vG6pbj+0uSbHaHbn9vTZMVogDQmZpboR4VGqCQQOdpAlaoA+hMy3YV6JmfdkiSAswmPXXxaEWGBBpcleezmE166Kzhariu9fEftimzmJEdADxHXlmNausODK2Yow7ACA6HQ/d8ukE78yskSQMSWn4+k0C9hc4YmaxxPbtIkvYWVGrOx+vlcDgMrgoAAKDzvf97un7f61xd3btruG46tq/BFXmPW4/vp1EpMZKkPQWVeuDzTcYWBMBvldfYXNvhwRZJkslkUlK0c5V6LoE6gE5SXFmrWe+tUf21mrr9xAGu4yUc3siUGF02sackqbLWpvvnbzS4IgBotLegwrUdE9Z4oRRz1AEY4b0V6fpkdaYkKTzIov9cMLLFryVQbyGTyaT/XDBSkfVX7n++NktvL08zuCoAAIDOta+8Ro98vcV1+8Ezhykk0GJgRd4lwGLWExeOUliQ8+/svd/T9c2GbIOrAuCPKppp+S7JNe6srKZOZdXWTq8LgH9xOBy6+6P1rq4Yk/vE6cZjuFizte48aaC6RgRLkr7blKsfNuUaXBEAOLnPTz9xcIJreyNz1AF0ss3ZpbrP7cLDf547Qr27skK9Q/SMC9ej541w3b7/80188QMAAL9hszs0+/21KqlyBixnj+6uqf27GlyV9+nVNVx/P6NxJujdH69XDitBAXSyilpnoG42SaFuF0YlRYe4tnOZow6gg727Il3fbMyR5Fy5+N8LR8piNhlclfeJDg3UvacNdt2+b/5GVdYyWgiA8dLdAvUTBscr0OL8jqflO4DOVFZt1U1vrVJN/QiKyyf11Okjk1v1HgTqrXTK8CRdMdnZRqm2zq5b3l7NVfsAAMAv/Oe7rfp5W74kKS48SH89dfBhXoGDOX9sD80clihJKq606s4P1spuZ5wQgM5TXu0MWsKDA2QyNYZXiW6BOnPUAXSkHXnluv/zxlVCj547wjV2Aq13xshkTekXJ0nKLK7SY99sNbgiAGi6Qr1vtwj1j4+UJO3Mr1C11XawlwFAu3E4HJrz8Xrt3uccQTGse5TuOa315zQJ1I/AX04drOHdoyVJu/dVME8dAAD4vK/WZ+u5hTslSRazSc9cMsbVVhKtZzKZ9PDZw5UQ5fw7XLxjn179dbfBVQHwJ+X1Ld/d271LTVeoE6gD6Cg1dTbd+s5qVVudq4QunZiqk4YmGlyVdzOZTHrwzGEKsjhP9877bY8+X5tlcFUA/J17oJ4SG6YhyVGSnB3wtuaUGVUWAD/y5rI0fbHOOW4xMjhAz14yRsEBrR9fSaB+BIIDLHr2kjGueepfrMvWW8uYpw4AAHzTttwy3fnBWtftv54yWJP7xhlYkW/oEh6k/5w/ynX7sW+2ajNt7wB0koYZ6uH7BeqJbqtDGUcBoKM89s1WV7vf/vERuufUIQZX5Bv6dIvQvac3/l3++aN12pZLYAXAOGmFVZKkhKhghQRaNLQ+UJdo+w6g423ILNGDn29y3X7svBHqGRd+RO9FoH6EUuPC9JjbPPUHvtikDZnMUwcAAL6lpMqq61//XZW1zlZsZ4/urqum9DK2KB8ytX9XXTettySp1mbXrHfX0PYOQIez2x2qqP9e3z9QZ4U6gI62cGueXlns7MwTFGDWUxePVmhQ61cJoXmXTUzVOWO6S5Iqa2268Y2VjKsEYIiKmjrtK6+RJKXGhkmShiQ1Buobs8hTAHSc0vq56bU2Z0ekK4/qpZnDk474/QjU22Dm8CRdeVQvSQ3z1FdxgAoAAHyG3e7QrHdXa0+Bs0Xb0OQoPXz28CazdtF2d540UIPrTypszS3TP7/eYnBFAHxdpduFOxHBTUMs9xnqOSVVnVYTAP+QX1bTpPPRnJmDXMdBaB8No4UaQqtd+yp05wdrGVcJoNOlFzVt9y5Jg91XqGexQh1Ax3A4HPrzh+tcYydG9ojWX05p/dx0dwTqbTTnlEEa0cM5T31PQaXuZp46AADwEU/8sE0/bc2XJHUJC9QLl41l9VAHCA6w6MmLRik4oHHe5aJt+QZXBcCXNbR7l6TwoKYr1GPDghRocV44xQp1AO2psrZO177+u/aV10qSjhvYzbVQBe0rJNCiFy4bq6gQ53f8txtz9cKiXQZXBcDfpBU0BuoNK9SjQgJd25uzy2Szk6UAaH+vLN6trzfkSJKiQgL0zCVjFBTQtkicQL2NggMseubiMYqsP0D9cl223mSeOgAA8HLfbszRUz/ukCSZTdIzl4xxXVGO9jcgIbLJlbJ3frBWBfWt8QCgvZW7BeoR+7V8N5tNSohyrlLPKSVQB9A+rDa7bnprldamF0tyztL91/kj6XzUgVLjwvTkRaPV8Ff8r2+36Ncd+4wtCoBfaVgZKjUG6lJj2/cqq017Cio6vS4Avu3T1Zl66KvNrtv/Pn9ku5zTNDRQ//vf/y6TydTkv0GDBrker66u1s0336y4uDhFRETo3HPPVW5uroEVNy81Lkz/cpun/uDnzFMHAADea0deme54370V52BN6dfVwIr8wx8m99SxA7tJcrZDve3dNbLWz3kCgPbUZIX6foG61DhHvbjSqmq39vAAcCQcDofu/mi9FtZ3PooMCdC8qyaoa0SwwZX5vuMGxeu2E/pLkuwO6Y/vrFZWMeM8AHSOdLdAvWdcY5g1NNl9jjpt3wG0nx825eqOD9aqoZH4Lcf104yhie3y3oavUB86dKiys7Nd/y1evNj12O23367PP/9cH3zwgRYtWqSsrCydc845BlZ7cCcPc5unbrPrZuapAwAAL1RabdX1b6x0rV48fWSyrp3W2+Cq/IPJZNJj541QXHiQJGnxjn26+yPGCQFof01WqIccGKgnRoe6tnNo+w6gjR77dqs+WpUhSQqymPW/P4xjbnonuvX4/jqu/qLNwopa/d9bq1RTx8VSADre3sIDZ6hL0hDmqAPoAL/t3Keb3l7lGiVx6cRU3TFjQLu9/4G/OXeygIAAJSYeeHVASUmJXnnlFb399ts6/vjjJUlz587V4MGDtXTpUk2aNKnZ96upqVFNTWN7zNJS5xey1WqV1dqxAfefTuynVXsLtS6zVHsLKnXrO6v03MWjFGAx/LoFAPAZDd/lHf2dDvgju92h299do135zpZrgxIi9I8zBqmuru4wr0R76RJi0bMXj9Qf5q1UbZ1dH63KUFJUkG47oZ/RpXkF9hFAy5RWNP7OHBpgOuBnJj4i0LWdUViu7tFBnVYb0FHYRxjjtSV79fzCnZIkk0n693nDNDYliv8Pneyxc4bp7BeWKqOoSmvTi3XfZxv04BlDjC4L8CjsJ9pfWn0795BAs2KCza6/2wHxjeH6hsxi/s4BtNnajBJd99rvqq1zdno8bXii7j1l4GHPabbm+8fkMHDJy9///nf961//UnR0tEJCQjR58mQ98sgjSk1N1Y8//qgTTjhBRUVFiomJcb2mZ8+emjVrlm6//faDvuf9999/wP1vv/22wsI6fu5nQbX0r3UWVdmcA4rGd7Prkr52mRkJBQAAPNw36SZ9nWGRJIVZHLpjhE1dQwwuyk+tKTBp3jazHHIeRF7Ux6bJCaxUB9A+VuSb9OYO5/f9ub1sOjqp6ffLwmyTPtnjfPyyfjaN78b3D4DWW73PpNe2Nx7PnNfbpmmJfJ8YJaNCemK9RVaH8//HxX1tmhTP/w8AHcPukO5cZpHNYVJiqENzRjV2xnA4pL/+blFFnUkRAQ79Y5xNJvITAEcou1J6aqNFlXXOL5IhMXZdO9Culqx1rqys1CWXXKKSkhJFRR26g5KhK9QnTpyoefPmaeDAgcrOztb999+vadOmacOGDcrJyVFQUFCTMF2SEhISlJOTc9D3nDNnjmbPnu26XVpaqpSUFM2YMeOwfxntpe/IAl37xipZbQ6tyDdraL9e+svMgTKxVwCANrNarfr+++914oknKjAw8PAvANAin63N1tdL1kuSzCbpmcvGahpz0w1ziqSk3/bq4a+3SpI+2BOgE6eO1tH9+X9yKOwjgJYpWp4u7dgsSRo/eoROGdO9yePmjbn6ZM9aSVJi70E65WhGf8D7sY/oXEt2Feit11fJIWdg+3/H9Nbs6f0Nrgrd+mXpro83SJI+2huoC2dMaDLLGPBn7CfaV3ZJtWxLf5YkDUmN1ymnjG7y+Af5K/XrzgKV15k0/ugTFB8ZbESZALxcWmGlHnp5hSrrnF3YJvTqolf+MEYhgZYWvb6hy3lLGBqoz5w507U9YsQITZw4UT179tT777+v0NDQQ7zy4IKDgxUcfOCXb2BgYKftCI8ZlKinLx6tm95aJbtDmrckTXERIfrjCfziAADtpTO/1wFf98W6LN310XrX7TtPGqjjBycZWBEk6fpj+imntFav/rpbNrtDt767Vu/dMFnDukcbXZrHYx8BHFp1XeOKxOiw4AN+XnrEhru288pr+XmCT2Ef0fE2ZpXoprfXympzftdcMK6H7jp5MAtNPMAFE3pqXVap3lyapto6u255d60+v2WquoQz2gNowH6ifWS7hVQ9u4Yf8Hc6rHu0ft1ZIEnallep7rERnVofAO+XW1qtK19bqbwyZ5g+oke0XrlyvCJDWv4d3prve48a7h0TE6MBAwZox44dSkxMVG1trYqLi5s8Jzc3t9mZ657m5GFJeuSc4a7b//l+m95YutfAigAAAA70zYYc3fbuGtnrs5XLJqXq/47pa2xRcLnn1MGaOcx57FtRa9NV81Yoo6jS4KoAeLuKmsY5cuHBB15nnxTdeIF7dkl1p9QEwDekF1bqyrkrVF7/PXPCoHg9fPZwwnQPcu9pQzQqJUaSlFFUpSvnNf7/AoD2klbY+Htrz9gDR/EOceuOsSm75StEAUCSiipqddnLy5ReWCVJ6h8foXlXTWhVmN5aHhWol5eXa+fOnUpKStLYsWMVGBioBQsWuB7funWr0tLSNHnyZAOrbLkLx6dqzsxBrtt/+2yD5q/NMrAiAACARj9sytUf31klW32aftH4FD1wxjBOeHoQs9mkxy8cpbE9u0iS8stqdOXcFSqptBpcGQBvVn6YQL1rRJDM9buCHAJ1AC1UUF6jP7y6XPn1q4RGp8bomUvGKKAlAyzRaYIDLHr+sjHqVt9eeW16sa6et0JVtbbDvBIAWs49UE+NOzBQdx83sTGrpFNqAuAbyqqtumLucm3PK5ckpcSG6o1rJiq2gzvuGHpEe+edd2rRokXas2ePfvvtN5199tmyWCy6+OKLFR0drWuuuUazZ8/WTz/9pJUrV+qqq67S5MmTNWnSJCPLbpUbjumrG+tXeTkc0uz31mjh1jyDqwIAAP5u4dY83fTWKlcrznPH9NDDZw+X2UyY7mlCAi363x/GqXdXZwvmHXnluv6N31VTx0lPAEfGfYV6RDOBeoDFrPjIEEmsUAfQMuU1dbp63grt3lchSerbLVyvXjFeoUEtm1+JzpUUHao3r5momDDnKq7luwt145srOb4E0G6aBOrNrFDv3TVCIYHOeGpTFivUAbRMtdWm617/XesynBfixEcG661rJikxOqTDP9vQQD0jI0MXX3yxBg4cqAsuuEBxcXFaunSpunXrJkl6/PHHddppp+ncc8/V0UcfrcTERH388cdGlnxE/nzyQF08IUWSVGd36MY3V2rl3kKDqwIAAP7q1x37dP0bK1Vrs0uSzhiZrMfOG0GY7sFiw4M076rxiqu/2nbZ7kLd+cE62e2Ow7wSAA5UUdMYmIQHNx92NZyQKKioUW2dvVPqAuCdCsprdMn/lmpt/YnNhKhgvXb1BOZye7iBiZF6/eoJrgurFm3L123vrFGdje98AG3nHqj36HJgoG4xmzQo0blKfU9Bpcqq6cIG4NBKqqy64tXlWrrLma/GhAXqjWsmNtsFoyMYGqi/++67ysrKUk1NjTIyMvTuu++qb9/GmZ0hISF69tlnVVhYqIqKCn388cdeMT99fyaTSf84a7hOGe6svdpq11VzV2gzs0EAAEAnW7qrQNe8tsIVjswclqj/XjBSFsJ0j9czLlyvXDnedRX/52uz9Oi3WwyuCoA3Kj/MCnVJSqoP1B0OKa+MVeoAmpdZXKXzX1ziWiUUHRqo166e0Gx4As8zokeMXnU7vvxmY47u+pCLNgG0XXp9oJ4QFayQwOYv4HSfo74lp6xT6gLgnbKKq3T+C79p2W5nmB4eZNG8qyZoYGJkp9XQ/G/O+5k9e3ar3/iee+5RbGxsq1/nqyz18y/Lqn/XL9v3qbS6Tn94dbk+uvGoTrt6AgAA+Lff9xTq6nkrVG11huknDknQUxePZq6lFxmVEqNnLh6j69/4XXaH9OKiXYoJDdL/Hdv38C8GgHoVh5mhLqlJy7yckmrCMQAH2JFXpstfWe4aDZEQFaw3rpmoAQmdd2ITbTehd6xeunycrn3td9Xa7Pp4dabCgi168MxhMpm46BZA61XU1Glfea2k5tu9N2gyRz2zRON7kScBONCWnFJd+eoK5ZQ6jzljw4P06pXjNSolplPraFGg/sQTT2jy5MkKCmpZq6bFixfrlltuIVDfT3CARS9cNlaXvrxMa9KLlV9Wo8teWaYPbpyshKiO7+8PAAD815r0Yl05d4Uqa51tfo8b2E3PXDJagYTpXmf6kATdf+Yw3fvpBknSo99sUWVtnWafOICTngBapGGFenCA+aD7gSS3QJ056gD2tya9WFfNXa6iSmeL3t5dw/X61ROUcojgBJ7r6AHd9PQlo3XTW6tkszv05tI0hQcF6O6Zgzi+BNBq6UXu89PDD/q8IUmNgfomuvkCaMZvO/fphtdXqqz+d9iecWF67aoJ6tX14N8tHaVFgbokffLJJ4qPj2/RcyMjuRL1YMKDAzT3yvG64MUl2p5XrrTCSp37/G9645qJ6m3APwAAAOD7NmSW6PJXlrkClGn9u+r5y8YqOKD5tmvwfJdP6qmyaqse+2arJOnpH3eoosame08bzElPAIfVsD84WLt3SUqMDnVt5xCoA3CzePs+Xf/G764LNYcmR+m1qyeoa0SwwZWhLU4amqj/nD9St7+/Rg6H9OLPuxQeHKBbT+hvdGkAvMzeAvdA/eAXWg1KjJLZJNkdBOoADjR/bZbueH+NrDbnKJqRPaL1ypXjDTvmbNGSpLlz5yo6OrrFb/riiy8qISHhiIvydV3Cg/TGNROVEus8QZFRVKXznv9N6zKKjS0MAAD4nOW7C3Xpy8tUVu0MTyb3idNLl4876AwzeI+bju2nv58+xHX71V93a87H62Vj5iWAw2ho+X6wdu8SK9QBNO+r9dm6at5yV5g+sXes3r1+EmG6jzhrdHc9dNZw1+3/fr9NryzebWBFALxRw/x0SUqNCz3o80KDLOrTLUKStC2nXFabvcNrA+D5HA6HXvp5p259Z7UrTD9+ULzeMfiYs0WB+hVXXKHg4JYXeckllyg8nNXWh5IYHaIPbzxKgxKdq/kLKmp10UtL9fO2fIMrAwAAvuLT1Zm67OVlKqlytuKc0CtWr1w5TqFBhOm+4sopvfXYuSNkrl+U/u6KdN3+3hpORAA4pIoaZxB2qEA90W0sWU5pVYfXBMDzvbVsr25+e5XrxOaMIQl67eoJigwJNLgytKdLJqbqnlMHu24/+MUmvbs8zcCKAHibtMKWrVCXGueo19rs2pFX3qF1AfB8NrtD93++SQ9/tcV138UTUvTS5WMVFtTipusdos1DM3ft2qWNGzfKbuekXWslRIXovRsma0Jv56z5ylqbrp63Qp+tyTS4MgAA4M0cDoee/GG7Zr23RrX1werRA7rp1avGG37wifZ3wfgUPXnRaAXUp+rz12bpprdWqdpqM7gyAJ6ots7u2jdEBB/8AquEKFaoA3ByOBx69qcd+usnG+Sob4Rz/tgeeu7SMXQ98lHXTuujWdMbW73P+WS9Xvttj3EFAfAq7oF6ymECdfc56huzaPsO+LNqq023vL1K89yOOW6fPkAPnz1cAZY2x9lt1uIKrFar7rvvPp1++ul66KGHZLPZdPHFF6t///4aMWKEhg0bpj179nRgqb4pOjRQr189QScNdbbIr7M7dNu7a2inBAAAjkhtnV13frBOj/+wzXXfxRNS9coV4w45Kxfe7fSRyXrhsrEKCnAe3n+/KVfXvf67KmvrDK4MgKdpaPcuHXqFelCAWV0jgiQxQx3wZza7Qw9+sVn/+nar674bju6jx84b4REnNtFxbjuhv66b1luS5HBI983fqEe/2SKHg/FCAA6tIVAPDbSo22HaMw9JbgzUNxGoA36roLxGf3hlub7ekCNJsphNeuzcEbpten+ZTCaDq3Nq8ZHv3Xffreeff16JiYl69dVXdc4552j16tV6++239e677yogIEB//etfO7JWnxUSaNFzl47VxRNSXfc9+MUm/fNrDlIBAEDLlVRa9YdXl+mjVRmu++bMHKSHzx6mQE54+rzpQxI098rxCqtv6f/L9n264tXlKq22GlwZAE9S3sJAXXKOKpOkvLIa2ez8bgr4m6KKWl05d7le/bVx0cfdMwdpzimDPebEJjqOyWTSX04ZrJuP6+u67/mFO3XH+2sZLwTgoOx2hzIKneOCUmPDDru/aLpCvaRDawPgmVbsKdSpTy3W8j2FkqSwIItevmKcLhifYnBlTbV4mdKHH36oefPm6ZRTTtG2bds0aNAgffnll5o5c6YkKT4+XpdeemmHFerrLGaTHj57mOIjg/Xkgu2SpBcW7VR+WY3+ee5wToIDAIBDSiuo1FXzlmtnfoUkKTjArMcvHKVThicZXBk605R+XfXGNRN05dwVKquu04o9Rbr0f8v02tUTFBseZHR5ADxAhVvniojDjAFJjArVhsxS2ewO7SuvadIGHoBvW59RohvfXKnMYmco0nDe6sLxqYd5JXyJyWTSn04apMSoEP1t/kY5HNLHqzOVX16j5y8bSwcsAAfIKa12jRc6XLt3SYqLCFZiVIhySqu1KbtUDoeDi7YAP2G3O/TSL7v0r2+3ui7g7hoRrFevHKcRPWKMLa4ZLU5ps7KyNHLkSEnSgAEDFBwcrH79+rkeHzBggHJyctq/Qj9iMpl0+4kD9I+zhqlhn/HRqgzd8MZKVdUyAxMAADRvVVqRzn7uV1eYHhcepHeun0SY7qfG9ozVO9dNUpewQEnS+swSXfTSEtcJcQD+raUt3yUpKZo56oA/em9Fms594TfXsUPXiCC9de1EwnQ/dvnkXnr+0sbxQr9s36eLXlqivDL2DQCacp+fntqCQF2Shta3fS+rrlNGEb+3Av6gqKJW177+u/759RZXmD6xd6y+unWqR4bpUisCdZvNpsDAQNftgIAAWSyWxjcym2lP3k4um9RTz10yRkH1q9J/3JKnS15eqvyyGoMrAwAAnuar9dm6+KWlKqiolST17RauT26aojGpXQyuDEYa1j1a798wWfGRznl123LLdcbTi7V8d6HBlQEwWnlN48XaESEta/kuSTklnNwEfF211aa7P1qnP3+0XrV1ztWFY1Jj9MUfp2lSnziDq4PRTh6WqLeunajoUOf54Q2ZpTr3+d+0K7/c4MoAeJKmgXpoi17jPkd9I3PUAZ+3Kq1Ipz29WD9uyZMkmUzSLcf101vXTlS8B3dFa1Uf8W+//Vbz58/X/PnzZbfbtWDBAtftb7/9tqNq9EszhyfptasnKLJ+xcDqtGKd9vQvWrm3yODKAACAJ7DbHXpu4Q7d9NYq1dSf8JzUJ1Yf/98Upca17Cpw+Lb+CZH64MbJ6ln/76GgolaXvrxUby9LM7gyAEZyX6EeEWw5xDNZoQ74k4yiSl3w4hK9uyLddd8Vk3vq3esnN7m4Bv5tfK9YfXjjZCXX/5tIL6zSeS8s0eo0zlcCcEp3D9RbeG5iqFugvok56oDPcjgcevmXXbrghcYuirHhQZp31QTdedJABXj46OtWDbq54oormty+4YYbmtxmtkX7mtw3Tu/dMFlXzl2uvLIa5ZbW6MIXl+ieUwfriqN68fcNAICfKqyo1R3vr9FPW/Nd950zprv+ec4IVxtGQJJ6xoXrs5un6Ja3V2vxjn2y2hz6yyfrtSm7RPedPlSBHv7LCoD2V17d8pbvTVeoE6gDvuqX7fm69Z3VKqq0SpJCAs165JzhOnt0D4MrgyfqnxCpj2+aoivnLteWnDIVVtTqkv8t07OXjtbxgxKMLg+AwY6k5fuQpGjX9qZsVqgDvqikyqq7Plyrbzfmuu4b36uLnrp4tJKiW9bNwmgtPoNmt9sP+5/Nxpzv9jYkOUpf3DpVE3vHSpLq7A79/fNNuu3dNaqsrTvMqwEAgK9ZuqtAM5/8uUmYfvv0AfrP+SMJ09GsmLAgzbtqvK6Z2tt135tL03Tpy8tUUM5IIcDflDdZoX64GeqNJzZYoQ74HrvdoWd/2qErXl3uCtNTY8P08f9NIUzHISVGh+j9GydrUh/n+coqq03Xvb5S7yynExLg79wD9R5dWhao9+gS6urUu4mW74DPWZfh7MDtHqbfcEwfvX3dJK8J06VWtnyHMeIjQ/TWtRN1/dF9XPfNX5uls579lTlFAAD4CZvdoacWbNcl/1uq3FJnCBoXHqTXrp6g26b3p3MNDinAYta9pw3Rv88fqaD6VenLdxfqjGd+1UZa6gF+xb3le3jQYVaoR7FCHfBV+WU1uu713/Wvb7fK7nDed8KgeH1+y9Qms2yBg4kKCdRrV0/QaSOSJDl/X5nz8Xr99ZP1qqlj0RXgr9IKnIF6YlSIQgIPPV6ogdls0uD6fU9WSbWKKmo7rD4Anae2zq4nftimc5//TemFzhbv0aGBeuWKcZozc7DXdU1sUcv3+fPnt/gNzzjjjCMuBgcXYDHrL6cM1uiUGP3pw3Uqr6nTttxynfHMr/r3+SN08rAko0sEAAAdJK+0WrPeW6Pfdha47juqb5yeuHCU4qOYaYmWO29sD/XtFq4b3lipvLIaZRZX6bznl+jf54/UqSM4ngT8QXlty1u+hwZZFB0aqJIqq7JLqzq6NACd5It1Wbr30w2uVekmkzR7+gDdfFw/mc1cpImWCw6w6KmLRishKkSvLN4tSXprWZo2ZJbo2UvHtHh1KgDfUF5Tp4L6MLyl7d4bDEmK0vLdhZKcbd+n9Ova7vUB6Dxr04v154/WaUtOmeu+0akxeuaSMeoe4z2r0t21KFA/66yzmtw2mUxyOBxNbjeg7XvHmjk8SQMSI3XjGyu1Pa9c5TV1uvHNVbrh6D7600kDFeBlV3QAAIBDW7QtX7PfW+P6pdRskmbVn/C0cMITR2B0ahd9/sepuv6NlVqbXqwqq003v71Km7P7afaJAziRDvi4ila0fJekpOgQlVRZlVtSI7vdwXcE4MUKK2p172cb9OW6bNd9XSOC9O/zR+rYgfEGVgZvZjabdO9pQzQwMVL3frpBNXV2rc0o0WlPL9aTF43WMQO6GV0igE6S7tbuPaWVgfpQt+4oG7NKCNQBL1Vttenx77fpf7/scnVBsphNuvGYPrrthAFePa6yRZW7z0n/7rvvNGrUKH399dcqLi5WcXGxvvrqK40ZM0bffPNNR9cLSX27RejTm6fo9JHJrvte/HmXLntlmfLLmIMJAIAvsNrsevSbLbri1eWuMD0hKlhvXzdJt57QnzAdbZIQFaL3rp+kc8Z0d933zE87dN3rv6uQ9nqAT6uoabwIPjz48G04E6OdnVBqbXYVVvL9AHir7zbmaMbji5qE6acOT9K3s44mTEe7uGBcij6+6SjXqtTiSquunLtcT/6wXXa74zCvBuAL3Oent3qFulugzhx1wDst312omU/+ohd/bgzThyRF6bObp+hPJw3y6jBdauEKdXezZs3SCy+8oKlTp7ruO+mkkxQWFqbrr79emzdvbtcC0bzw4AA9ddEojUmN0UNfblad3aGluwp16lO/6NFzR+i4QfwyBACAt8ooqtRt767Ryr1FrvuOG9hN/z5/pOIigg2sDL4kJNCi/5w/UkOTo/XQl5tkd0gLtuTp5Cd+1n8uGKlp/VlNBPii8iNYod4gp6RaXdkPAV6lpNKq+z/fqI9XZ7ruiwkL1INnDmuyUANoD0OTo/X5LVN1xwdr9MPmPDkc0uM/bNPq9CI9fsEodQkPMrpEAB3IfYV6alzrWjr3j49UoMUkq82hjQTqgFcpr6nTo19v0RtL97ruC7KYddv0/rr+6D5eNyv9YFr9p9i5c6diYmIOuD86Olp79uxph5LQUiaTSVdN6a13r5+khCjnSY28shpdNW+F7v5oncqqrQZXCAAAWsNud+j1JXt00uM/u8L0ALNJfz1lsF65YjxhOtqdyWTSNVN767WrJ6hLWKAk5/Hk5a8s1z++2KSaOsY5Ab6mScv3kMMH6olRjSdDc0qqO6QmAB1j4dY8zXhiUZMwffrgBH13+9GE6egw0WGBeunycfrTSQPV0FRr4dZ8nfb0Yq3LKDa0NgAdqy0r1IMCzOofHylJ2plfrmorv4sC3mDRtnyd9PjPTcL0Makx+uq2qbr5uH4+E6ZLRxCojx8/XrNnz1Zubq7rvtzcXP3pT3/ShAkT2rU4tMy4XrH64o/TdLTbTKJ3V6Tr5Cd+0W879hlYGQAAaKkdeeW64MUl+ttnG1VR6/zFsUeXUH1w42Rdd3QfZtaiQ03r303fzDpa0/o3zql7efFunfnMr9qWW2ZgZQDaW0OgbjZJoYGHb/nuvkI9u5RAHfAGZdVW3f3ROl05d4VyS52jASNDAvTv80fqf38Yq/jIkMO8A9A2ZrNJNx/XT29cM1Fx9avSM4urdN7zS/T2sjQ5HLSAB3zR3gL3QD281a9vmKNud0hbcvg9FPBkuaXVmv3eGl3x6nJlFldJcv5++bfThuiDG49Sv/oLZHxJqwP1V199VdnZ2UpNTVW/fv3Ur18/paamKjMzU6+88kpH1IgW6BYZrNeuGq9Hzhmu8CDnSZHM4ipd8vIy/X3+RlXW1h3mHQAAgBGsNrue+XG7TnnyF/3u1uL94gmp+uq2aRqd2sXA6uBPEqJC9NpVE3TvaUMUVH8F8ZacMp3+9GK99tseTnwCPqKh5Xt4UIBMpsNfrJXYpOV7VYfVBaDtHA6HPl2dqRP+s0jvrkh33X/0gG767vajdd7YHi36uQfay5R+XfXFrVM1OjVGklRrs+svn6zX7e+tUUkVnTUBX9PQ8j000KKuEa0f8cAcdcDzVVtteubH7Tru3wubdEE6qm+cvp11tK6e2lsWH10U1OoZ6v369dO6dev0/fffa8uWLZKkwYMHa/r06RyUG8xkMuniCama2q+r/vThWi3dVShJmvfbHi3cmqf/XDBSY3vGGlwlAABosC6jWHd9uK7Jlde94sL0yDkjNLlvnIGVwV+Zzc4W8Ef1jdOsd9doa26Zaursum/+Ri3cmqfHzhupbpGMHgC8mStQb8H8dGm/Feq0fAc81qasUv19/kYt31Poui88yKJ7Thuii8ancM4OhkmKDtV710/Ww19t1rzf9kiSPl2TpWW7C/XYeSM0rX+3Q78BAK9gszuUUeS8+DI1NuyI9jtDk6Nd21+sy9LFE9h/AZ7C4XDoi3XZ+ufXW1wr0iUpKiRAfzllsC70g+PNVgfqkjO4nTFjhmbMmNHe9aAdpMSG6e1rJ+m1JXv06DdbVG21a09Bpc57YYmun9ZHt584QCEtaO0HAAA6RlWtTY//sE0v/7JL9vpFv2aTdN3RfXT7dPbTMN7gpCh9dssU/fPrLa4Tnz9tzdfJT/ysf50/QscPSjC2QABHrKLGOVYkPLhl+5qmK9QJ1AFPU1xZq/9+v01vLt3rOq6UnLPS7zt9iFJaOcMW6AhBAWb9/YyhGtOzi/768XqV1dQpu6Ral7+yXJdNStVfThmssKAjOk0NwEPkllar1maXpCPe94xOjVFqbJjSCiv1284CLdqWr2MHxrdnmQCOwNr0Yj3wxSatdOusaTGbdOnEVN0+fYC6hLe+I4U3alHL96eeekrV1S3/xfmFF15QWRkzLoxkNpt01ZTe+urWaa62Sg6H9OLPu3T604u1LqPY0PoAAPBXv+3Yp5Of/Fkv/dwYpg9OitJnN0/VnJmDCdPhMUICLfr7GUM176rx6hrhXJVeUFGrq+f9rjkfr6dNJ+CFHA6HKurHgUW0cIV6ZEiga6wYgTrgOWx2h95Znqbj/r1Qry9pDNN7dw3X3KvG6+UrxhGmw+OcMTJZ39x+tKb0a+zG9ebSNM188hf97tZdAYD3SSt0n59+ZPufQItZfzppoOv2P7/eIpud0WOAUXJKqjX7/TU689lfm4Tp0/p31de3TdMDZw7zmzBdamGgfvvtt7cqIL/rrruUn59/xEWh/fTpFqEPbzxKd88c5JqFuT2vXGc++6vu+XS9iitrDa4QAAD/kF1SpdvfW6NLXl6mvQXOXzSDApy/LM6/ZYqG94g+zDsAxjh2YLy+nTVN0wc3rgx4Z3mapv93kb5cl81sdcCLVNba1PAj29KW71LjKvXskmp+5gEPsCqtSGc9+6vmfLxeRZXOC9zCgiz688mD9M2saTqO1XzwYN1jQvXG1RP1wJlDFRLoPFe5t6BS57+4RI98tVnVVpvBFQI4EmkF7oF66BG/z6nDkzSy/vzIlpwyfbwqo821AWidaqtNTy2on5O+qnFOep9u4Zp75Xi9fvUEDUiINLBCY7ToN2iHw6ETTjhBAQEt+4W7qqrq8E9Cp7GYTbrxmL46bmC87vhgjTZklsrhcF4B+tX6HP355IE6f2yKzGbfnm8AAIARKmvr9OKiXXrx552qttpd94/v1UX/PHeE+naLMLA6oGXiIoL1vz+M01vL0vTwV5tVWWtTflmNbn57lY4fFK8HzhyqHl1YBQd4uor6+elS6wL1pOhQ7cyvUJXVptKqOkWHBXZEeQAOI6+0Wo99u1UfrmwaLpw+Mll/OWWQkqKPPMAAOpPZbNIfJvfStP7ddMf7a7QqrdjVWfOnrXn67wWjNKw7FxwD3sR9hXrPuPAjfh+z2aS7Zw7Wxf9bKkn6z3fbdPrIZLr5AZ2g2mrTO8vT9PzCncorq3HdHx0aqFnT++uyST0VaGnROm2f1KLfoO+7775WvemZZ56p2NjYIyoIHWdgYqQ+uWmKXl28W08u2K7KWpsKK2r154/W653l6XrwzGGsjgMAoJ3Y7Q59tjZTj369VTmljS1yY8ICdceMgbp0QioXs8GrmEwmXTapp44bFK/7PtugHzbnSZJ+3JKnJTsLNPvEAbpqSi8F+PEvV4CnK3cL1Fva8l3ab456aTWBOtDJiipq9cLPO/Xab3uaXKA5MCFS9585VJP6xB3i1YDn6t01XB/ceJT+98su/fe7baq12bUtt1xnPfurbjm+n24+rp9fn7gHvIl7oN7WkSOT+8bphEHxWrAlTzml1Xr119266dh+bS0RwEFUW216b0W6nlu4Q7mljUG6xWzSZRNTNcuP5qQfSocE6vBcgRazbjimr84c1V0PfbVZn6/NkiStSS/WGc8u1sUTUvWnGQP54QAAoA1W7i3SA19s0tr0Ytd9AfWrMG47oT9BBLxa95hQ/e8P4/TtxhzdN3+jcktrVGW16aGvNuuT1Zl65JzhGpkSY3SZAJpRUdPYRjc8uOWrfJLcAvXskioNTPS/9n6AEcpr6vTq4t3638+7VOZ2QUxkSIDuOHGALpvUkwvZ4PUaOmseO7CbZr+3VpuyS1Vnd+iJH7brmw05+sdZwzSuFwu3AE/nHqj36NL2jil/njlIP23Nk90hPf/TTl00PlWxZBZAu6qps+n9Fel69qedTRYDSdLJQxN1x4wB6u+Hrd0PpuWXpMOnJEaH6OmLR+viCSm677ON2p5XLodDentZmr5en627Th6kC8fRBh4AgNbILK7So19v0fz6C9YanDAoXn85dTDt3eEzTCaTTh6WpKP6ddV/vt2q15fulcMhbcou1VnP/aorJvfSHTMGKDKEi0cAT9J0hXrLfz6brFAvqT7EMwG0h2qrTW8u3avnFu5UYUWt6/6gALMum9hTNx3XV10jgg2sEGh/gxKj9OnNU/TMj9v17MKdstkd2pJTpvNeWKLzx/bQ3TMHKY5/94DHSq8P1BOjQtqlPfuAhEhdMC5F765IV1lNnZ5asF1/P2Nom98XgDNI/+D3DD330w5l7ff73YwhCbpten8NTaab9f4I1P3cUX276qvbpmner3v0xA/bVFFrU1GlVXM+Xq93l6fpvjOGakxqF6PLBADAo5VVW/XSz7v00s+7VFPX2IZzQEKE7j1tiKb172ZgdUDHiQoJ1P1nDtNZo7trzsfrtSWnTA6HNO+3PfpmQ47uOW2wTh2eJJOJizQBT1DRJFA/0hXqBOpAR7Ha7Prg9ww9tWB7k1VCFrNJF4zroT8e31/JMcxJh+8KCjBr9oyBmj4kQXM+Xq+NWaWSpA9WZui7Tbm66+SBumh8qiwsAAI8SnlNnQrqLwBLbWO7d3e3nzhAn63JUlX9hWZXHtVLvboe+Xx2wN/V1tn1wcp0PfvjgUH69MEJmjW9v4Z1J0g/GAJ1KNBi1nVH99EZo5L10JebXavq1maU6JznftNJQxP0p5MGql88rR0AAHBXXlOn137bo5d+3qWSKqvr/tjwIM0+cYAuGp9CG074hdGpXfT5H6fq1cW79fgP21RttSuntFq3vL1ac3vu0T2nDtZoLtIEDFdR2xioh7dmhnpUY4DHCnWg/dnsDn2+NkuP/7BNewsqmzx2xshk3X7iAPUmQIAfGdEjRvNvmaq3lu3Vv77dqrLqOpVUWfXXTzbo/RXp+sdZwzW8Byf8AU+RVtB+89PdJUSF6LppvfXUjztUZ3foX99t1bOXjGm39wf8RXFlrd5alqbXftujvLKaJo9NHxyv204YwH61BQjU4ZIQFaKnLh6tiyek6r75G7Qtt1yS9O3GXH2/KVfnje2hWdMHcDU0AMDvVdbW6fUle/Xiop0qqmwM0gMtJl15VC/dcnx/RYfS6hr+JdBi1g3H9NUpw5N0z6cbtGhbviRp5d4inf3cbzp9ZLLuOmlgu55gAdA6ZdVHFqi7r1DPLK5q15oAf1ZttenDlRn63y+7DgjSpw9O0B0zBmhwUpRB1QHGsphN+sPkXjp5WKIe+WqLPlmdKcm5AOiMZxfr8kk9dceMgfzeBXgA9/npPePa9/e964/pq7eWpamgolZfrsvWtVOLuFgbaKHd+yr06uLd+nBlhqqstiaPHT8oXrOm99eIHjHGFOeF2i1Q37Vrl2688UZ999137fWWMMjkvnH68tZpev/3dD35w3blldXI7pDe/z1Dn67J0hWTe+qmY/upS3iQ0aUCANCpqmqdbcZeWLTT1c5Mkswm6ezRPXTrCf3UM47VQ/BvKbFhmnfVeP20NU8PfblZO/MrJEmfr83StxtzdM3U3rrp2L7MVwcM0LTle8tPB8SEBSomLFDFlVYt2VWg7bll6p9ABzPgSJVUWfXm0r2a++tu7SuvbfLYUX3jdOdJAxm/B9SLjwzR4xeO0gXjUvS3zzZoe165HA7p9SV79dX6bM2ZOVjnjOnOiCHAQOlugXp7tnyXnMess6b3172fbZQkPfL1Fr13/SR+5oGDcDgcWr67UC8v3q0fNufK4Wh8zGySThqaqBuO6atRKTGG1eit2i1QLysr04IFC9rr7WCwQItZl07sqXNG99Dc33br+YU7VVZdp9o6u/73y269uzxdNxzTR1dP7a2wIBodAAB8W7XVpreXpem5hTu1r7yxNZLJJJ01qrv+eHw/9ekWYWCFgGcxmUw6flCCpvXvpneXp+nxH7arsKJWtXV2Pb9wp95fka7bGYsAdDr3QL01K9RNJpOundpb//5um2x2hx7+arPmXjWhI0oEfFpOSbVe/XW33l6WpnK3n0dJmtIvTjcf209H9etqUHWAZ5vcN05f3TZNry7erSd+2K4qq037ymt1xwdr9frSvfrrKYM1oXes0WUCfsl9hXpHdCS7aEKqXv11j3bvq9Dy3YVasDlP04cktPvnAN7MarPrq/XZevmX3VqfWdLksbAgiy4Yl6Krp/RWajt3kfAnJKE4pNAgi246tp8umZCq5xft1Lxf96imzq6ymjr9+7ttmvfbXt12Qj9dOD5VQQGcDAUA+JaqWpveW5Gm5xftVG5p0yD9tBHJuu2EfuoXzwo94GACLWZdPrmXzhjVXc/9tENzf92jWptdBRW1uufTDXrttz36yymDdezAbqwwADpBeU1jm7+IYEurXnvttD56e1maskqq9dPWfP28LV9HD+jW3iUCPmlHXrle+nmnPlmdKautcZmQ2STNHJakG4/py9xKoAUaRgydPjJZD36xSV9vyJEkrU0v1gUvLtH0wQm6e+Yg9YvnYmegM6V14Ap1yfmz/+eTB+rGN1dJkh75erOOHdiNi7MBOS/YfP/3dL2zPE3ZJdVNHkuMCtGVU3rp4gmpjEhpBwTqaJGYsCDNmTlYVx7VS0/+sF3v/54uu0PaV16jez/bqOcX7tT1R/fRRRNSFRLYuhMzANpm0bZ8PfzlZp0xKlk3H9fP6HIAn5BfVqM3luzRG0v3NpmRLkmnDk/SbdP7awCtboEWiw4N1JxTBuuyST316Ddb9MW6bEnS9rxyXTVvhSb0jtXsEwdoUp84gysFfNuRrlCXpJBAi/48c5Bue3eNJOmhLzfrqL5xnMgEDsLhcOi3nQWa99ueA9ptBgWYdd7YHrp+Wh/16sq4IKC1kmNC9fxlY13nQ7bmlkmSfticq5+25unC8SmaNb2/4iNDDK4U8A8NLd9DAy3qGtExY2JPGpqosT27aOXeIu3Mr9D7v2fokompHfJZgKez2R1auDVP7yxP049b8mR3NH18WPcoXTetj04ZnqRAfl9rNwTqaJWk6FD989wRunZaH/3nu62uK0GzSqr198836Zmfdujqqb11+aSezMUEOsmzP+3Q1twy/evbrTpjZHKHtFYC/MWOvHK9sniXPlqVqdo6e5PHThqaoFnTB2hwUpRB1QHeLyU2TM9cMkZXTSnSP77cpNVpxZKk5bsLddFLSzW5T5xuP3EA7TqBDlJe6xaoH8HortNHJOvVX/dobXqxtuaWcSITaEZJlVUfr8rQG0v3ald+RZPHIkMCdPmknrpySi+CPqAdHDOgm6b266qPVmboP99vVW5pjWx2h95elqZPV2fquml9dP3RfVp9ERmAlrPZHUovcgbqqbFhHdZ5zGQy6S+nDNK5zy+RJD3+wzadOSqZn2/4lcziKr2/Il3v/55+wGp0s0k6flC8rp3WRxN7x9IFsAO0+Ntm9OjRh/wfUFlZedDH4Hv6xUfo+cvGak16sZ5esF0LtuRJkvaV1+qxb7bq+YU7deVRvXTVlN6KDe+Yq9IAOGWXVLm2P1qVoVnTBxhYDeB9HA6HVuwp0ks/79QPm/OaPBZgNumMkcm6dlofDUkmSAfay9ieXfTx/x2lL9Zl6/EftrkChyW7CrTkxSWa0i9Ot08foHG9CNaB9uS+Qj3iCE4+ms0m3XvqYJ33gvNE5n+/36rTRyZxMTUgaVNWqd5Yukefrs5SldXW5LGEqGBdM7W3Lp6Qys8L0M4sZpMuGJ+i00cm69Vfd+v5hTtVXlOnylqbnlywXW8tS9PtJ/bXheNS6KoCdICc0mrXOJOOns08tmesTh6aqG825ii/rEYv/7Jbt03v36GfCRitzmbXj1ucq9EXbcs/YDV6UnSILhiXogvGp6h7TKgxRfqJFv8GfdZZZ3VgGfBWo1Ji9MqV47Upq1TPLdyhL9dny+GQyqrr9PSPO/TyL7t16cRUXXd0HyVEcfU10N4cDofyyxrnOn+4MkO3Ht9fZjNXoAGHU2ez69uNuXrpl11am17c5LHI4ABdMjFVV07ppaRoDkaBjmAymXT6yGSdMjxJ89dm6skftmtPgfMi3V93FOjXHUs0rX9X3X7iAI1J7WJwtYBvaEvL9wbjesXq1BFJ+nJdtvaV1+q5hTv155MHtVeJgFepqbPpmw05en3JXq3cW3TA45P6xOrySb00Y2gC7TaBDhYaZNHNx/XTReNT9FR9kF5nd2hfeY3++skGvbJ4t2ZNH6BThyfJwjkToN2kFXTs/PT93XXyQH2/OVc2u0Mv/rxT54zpTrdO+ByHw6ENmaX6bE2m5q/NUp7b+X+pYTV6gi6ZmKJjBsSzX+skLf4N+r777uvIOuDlhiRH6ZlLxmh2frmeX7hTn6zOVJ3doSqrTS8v3q3Xl+zVeeN66JqpvdW3W4TR5QI+o7ymTtXWxrbUGUVVWrq7QEf17WpgVYBnyy+r0fu/p+vtZWnKLK5q8lhSdIiuntJbF05IURSrh4BOYTGbdPboHjp9RLI+XZOlp3/crr31J2V+2b5Pv2zfp2MGdNPtJw7QqJQYY4sFvFx5jXPVbFCAWUEBRx7u3X3yIH2/MVe1NrteWbxbl0xI5UQm/MqOvHJ9tCpD769IV0FFbZPHIoIDdM6Y7rpsUk8NSIg0qELAf8VFBOv+M4fpyim99a9vt+ir9c5xlbvyK3TrO6v15A/bdOsJ/XXaiGQCCKAdNMxPlzonUO/TLUKXTEjVG0v3qrLWpvNfWKLXrp6ggYnsc+H99hZU6LM1Wfp0TeYBo4MkqXtMqC4an6Lzx6UoMZoFrJ2tVZek/z979x3eZn3v//8lyfLeeyfOdvZOTEJCBgkJm0BLS1vK4bQ95xdoIW1PS09bSlsOpf2eltIyTiktXSkUWvYIIYGEQPZetpM4ie3YlrdlyUPz94dsxSbKIrHl8Xxcly9Jt27JH0Gse7zu9/uzZcsWvf7663I4HFq0aJGuueaanhoX+qlhKdH6xW2TdN/Vo/S7Dcf0/PYytbs8crg9Wr21VKu3luqq0Sn68hVDNW9kClW0wCWq+cTVaZL00o5yAnXgE7xer7aU1OuvW09qzYEquT7RH2lsRqy+Om+Yrp2YQfUQECQhJqNunZatGydn6uXdp/Sb9UdUVu+76GVDcY02FNfoiuFJ+tr84Zo3Mpn5wIBPobNC/dO0e+8qJzFSd80dqv/bUCKHy6NH3ynUbz8/9XIMEeizmlqcem1fhf65s1x7PtHdSJJGpUXriwVDdfOUrEv+GwNw6fKSo/TkHdO0q7RBP3urUNtO1EuSjtXY9Y3n9+jxdUd078KRun4SwTpwKUp7OVCXpPuvHqUPj9ToRF2Lqqxtuu3pj/Xsl2doBlOGoR+qtbXrjb0VenVvhXaXNp7xvNlk0MIxqfrczFxdOTKFbVYQXfAe/ksvvaTPfvazioiIkNls1i9/+Us9+uij+ta3vtWT40M/lRUfoYduHK97Fo7Us5uO669bTsrWcfLmg6IafVBUo+EpUfryFUN1y9TsT91uEBjsAgXqbx2o1EM3jmNuPkC+E5//3FWuv209qWOfuLLTYJCuGpWif79ymK4YnkQ4B/QRZpNRn5meo5unZOmfO8v1m/VH/d0kPj5Wp4+P1WlMeoy+Om+Yrp+UyUUwwEXoPCaLCjNd8nutXDBCL+0oV53doTf2VequOfWaNoSTmBhYXG6PNh6p0T93ntLaQ76uDF2FGA26Zny6vjh7iGbmJbI/CfRBU3MT9MLXZuvjY3X69XtHugXr973QEawvGqHrJ2YyxzrwKXQN1HurY1FiVKhe+s8r9G/Pbde+8iZZ21z6wu+36vHPTdHScem9MgbgUjS1OPXeYYte21uhTUdr5f7kxOiSZuUl6qYpWVo+PkNxkZzn7wsMXq/3zP9TAUybNk0zZszQE088IZPJpEceeUS/+MUvVF9f39NjvCRWq1VxcXFqampSbGxssIczaDW1OvWP7WX60+YTKm/o3l43JjxEt8/I0ZcKhtImELhIb+yr0D2rd0uSws1Gf/v3R1dM0Gdn5AZzaD3G6XTqrbfe0vLly2U2szOBM3m9Xu0tb9LftpzU6/squk2LIEnJ0aH6zPQcfY72tEC/4HB59NLOcv1u4zH/HOudMuLCdffcPN0+M1fRYSFsI4DzGPX9t+VweTQmPUbv3Dfvkt/vr1tO6vuvHJAkTcqJ18v/eQVdyNBnXcw2orDKqn/uLNfLuytUazvzIub8jFitmJqlGydnKSUmrKeGDOAy83q92lxSp8feO6Jtx7uf085LjtI9C0boxskE64MVxxKfzo1PfKS9ZY0yGKTDP75G4eZLv3DzQtnbXfqPv+7Uh0dqJfnmlf7pTRP0+VkD85wo+jeLtU3vHqzSmoMWbSmpO6N7piSNSY/RTVOydMOkTGXGRwRhlIPPxWTIF1wWXFRUpBdeeEEmk+8L8Zvf/KZ++MMfqrq6WqmpqZc2Ygx4cRFmfWXeMP3b3DytPWTRHz86rq0dO67NbS498+FxPbvpuBbnp+nLc4aqYBiVgsCF6FqhfvuMXD338QlJ0os7ygdsoA6cTXVzm17bU6GXdparsKr5jOdn5SXqjtlDdM249EuaNxZA7woNMerzs3L12Rk5WnuoSk9vKPG32q1satNP3zysX687ojtmDdEXZmYFd7BAH+Z0e+Rw+S4yu1ztqG+fkaM/bz6hYotNe8sa9fq+Ct04mb9D9E8lNTa9tb9Sb+6v0uFK6xnPJ0WF6sbJWVoxLUvjMuOCMEIAl8pgMOiK4cm6YniyNh+r02PvFfvPTx6vteubL+7V4+uP6N+vHKZbp2YrIrT3gkGgvzrVUTyXGhPWq2G6JEWFhejZO2fov17aq1f2VMjjlb738n7VNLfr64tGkC8g6I7X2rXmYJXeOVAVcMogydft+YbJmbppcpZGp8f07gBxUS74KLqlpaVbOh8aGqrw8HDZbDYCdVwwU0c7tGvGp+tgRZP+9PEJvbKnQg6XRx6v9O4hi949ZNGw5Ch9dkaOVkzLVnI0V3sDZ9O1WmLhmFR9fKxWxRabdpxsUEmNTcNSooM4OqDntTndeveQRf/aVa6NxTX65MWdMeEhWjE1W3fMytXINHZKgf7Mtx+ZoaXj0rX9RIN+t/GY3jtcLcl3gebTG47p2U0lmpJo1JAKqyYPSQryiIG+pXP+dEmXbcqtEJNR/33tWN35h22SpEffLtTScem9fjIV+LSO19r11v5KvbGvMmCIbjYZtGhMmlZMy9ZVo1OYZgQYQAqGJ6lgeIG2lPhawW8uqZMknaxr0Q9eOaBfrS3WF2cP0ZcKhiiJc5NAQC63R3V237nJtNjwoIwhNMSoX35mslJiwvTMh8clSb96r1jVzW368Y3jmW8avcrj8Wr/qSa9d9iiNQerVGyxBVwvOyFCS8f5crJpuQl0+eonLuoo+ve//72io0+HMy6XS88995ySk5P9y77+9a9fvtFhQBuXGaef3zpJ37lmjP6+rVR/2XJSFqtvA1xSa9cjbxfqF2uKdPXYNN0+M1dzRySzAQQ+oWuFempsmG6blqOH3zosSXppZ7n+65oxwRoa0GO8Xq+2n2jQv3aV6819lWruEhB0mpwTr8/PzNV1kzIUGXp5QgMAfYPBYNDMvETNzEvUEUuznvmwRK/srpDD7ZHT7dW2GqNuemqLpubG60sFQ7VsQrrCQgj3AFuX7eXlqlCXpPmjUjR/VIo2FNeooqlNz246rpULRly29wcut5N1LVpzuEZv7qvUoQAhuuSbwuCWjnabCVGhvTxCAL1p9rAkzf5qkraW1Ok3649q01Ff6+h6u0O/XndET284plunZevfrxymvOSoII8W6FtqbQ51TiicGhOcQF2SjEaD/vvasUqNCfefF/3b1lLV2Rx67PbJXOyJHlVra9eHR2q0oahGG4/Uqt7uCLjemPQYLRmXrqXj0jQ2I5YOCv3QBc+hPnTo0PP+DzYYDCopKbksA7tcmEO9/3C6PXr7QJX+vrXUf1VoV1nxEfrM9BzdNj2b+SOADnf9cZveL6qRJO38/mJ5vNLsR9bJ7fEqPTZcH3134YC7EIU5rQavYzU2vbanQv/aXa6y+tYzns+Kj9DNU7J089QsDac7AzCoVFvb9MePT+ivW06qua37RTZJUaG6fWaOPj9riLLYh8QgVlTVrKWPbZQkfWZ6tn5+66TL9t5HLM265tcfyu3xKirUpPe/fVVQT6oCXXm9Xh2ssGrtwUr9c+tRldsDHx9NyonXtRPStWx8hnISI3t5lAD6igOnmvTMhyV6Y1+l3F1aoBkM0pKxafrqvGGaNiQxiCNET+F808XbV96oG377kSTpjlm5evjmCUEekfTy7nJ9+8V9/vmpZ+Yl6pkvTVdcBP9PcXm43B7tLmvUhqIabSiu0f5TTQHXMxikqbkJWjouTUvGpmsoF2X1ST0yh/qJEycudVzAOZlNRt0wKVM3TMrUiVq7XthRphd3lPtbWp9qbNWv3ivWr9cVa/6oFH12Ro4WjEml4giDWk3H34fJaFBCZKiMRoMWjE7Ve4ctqrK26cMjNbpqNNNyoP86VmPTW/sq9eb+yoDzokeFmrRsQoZumZql2XlJtEgCBqnU2HB955ox+trcIXr4b2u11x6n4mpfa7U6u0NPvH9MT31wTIvz0/SlgqGaMyKJq8Ex6HSvUL+8JxRHpsXoczNz9NctpbI73Prlu8X62YqJl/V3ABejzenW5mN1eu+wResLq1XZ1NbxTPfv/knZcbp2YgYhOgC/8Vlx+vXtU/Rf14zRHzYd1/PbfNs2r1dac9CiNQctmjYkQf8+N09Xj01TCFNBYBDr7DYrBbdCvaubp2QrMSpM//nXnWpxuLXteL1ufvIj/ffyfC0ck8pxIC6a1+tVSa1dW0rq9NHRWn14pPaMC/k7RYeFaM6IJM0flarFY1P7zN8FLg96oKJPGpocpe9cM0arrh6l9YXVen5bqTZ0zI3r8UrvF9Xo/aIaxYaH6NqJGbpxcpZmDk0kSMGg09nyPTk61P/v/7bp2XrvsEWSr+07gTr6m/OF6AaDNHdEsm6ZmqWl49Jp6Q7ALyosRHPTvXp4WYF2lzfrz1tOas2BKrk8Xnm80ruHLHr3kEXDUqL0hVlDdPOULFr5YtCwdwvUL/9FyfcvHqVXd1eoud2lF3aU6faZuZqcE3/Zfw9wNtXWNq0vrNZ7h6v10dFatTrdAdebmBWraydmavkEQnQAZ5cVH6EfXDdWX180Uqu3luqPHx1Xdcc5mJ0nG7TzZIMy4sJ1x6xcfXZGrlJimGcdg4/F2ua/nxbbd/4G5o9K0d+/Mlt3Pbdd9XaHSmrsuvtPOzQ1N17fXjpGBcOTgj1E9GFdA/QtJfXaUlLXbdrVTxqbEav5o1N01agUTR2SIDMXWg1YF3wGevPmzaqrq9N1113nX/bnP/9ZDz74oOx2u2666Sb95je/UVhY3/niRP9nNhm1dFy6lo5LV0Vjq17cUa5/7CjTqUZfq19rm0t/31amv28rU2ZcuG6YnKWbp2RpdHpMkEcO9DyPx6tam29Olq4HbgvHpCopKlR1dofePWRRU4tTcZG0NULfdr4QXfLNi37thAxdNylDGXG0bQZwdgaDQbOGJWnWsCRZrG36+7ZSrd5a6j8JWlJj14/fOKSfvV2oq8em6bbp2bpyZMqAmyYF6KproB51GedQ75QUHaZ7Fo7QI28XyuuVPve7LfrpTeO1Ylr2Zf9dgCQ5XB7tKm3QR0drtbG4RnvLA7fbDA0x6orhSbpqZJK8FQf0hZtn08oXwAWLizDrP68arrvn5unVPaf0zIclKrb4OiFVNrXp/71brF+vO6LlEzL0pYIhmpqbQAUsBo3qLiFjah8K1CXfVC7//M8rdN/zu/37CLtKG/W5Z7boypHJ+taS0ZrExZ/QxQfo8ZFmXTkyRfNHpWjeyGSlxlKFPlhc8FH0j3/8Y1111VX+QH3//v26++679eUvf1n5+fn6xS9+oczMTP3oRz/qqbFikMuMj9A3Fo/UPQtH6KOjtXpl9ym9c7BKLQ7fVecVTW16esMxPb3hmPIzYnXT5EzdMDmT0AUDVkOLwz+fV0r06Z1Ws8moGydn6Q8fHZfD5dFre0/piwVDgzRKIDCX26NdpY1ad9ii9w5bdKzGHnC9zhB92YR0ZSdQQQTg4qXFhuu+xaO0csEIvXvQoj9vPqGtx+slSQ63R2/u913MkxEXrhVTs3Xb9GwNSWJuMww8th4O1CXpy3OG6s39ldpX3qRWp1vffHGvth6v00M3jFdEKFN14dJ4vV4VW2zadLRWm47UaOvxev/5gE9Kjg7TwjEpWpSfprkjkhUVFtIxN+6BXh41gIEiNMSo26bn6NZp2dp4pFZ/2XxC6wqr5fVKTrdXr+6p0Kt7KjQ2I1Z3XjFEN0zKYtuHAa+6S4V6X2xtnZccpVdWztGagxb977tFOtIxLdiHR3xtu68Zl65vLhmlkWkU5w0mzW1O7Stv0q6TDdpd1qjdpQ1qaHGedf2oUJNm5CVq9rAkzR6WpAlZcVyMP0hd8FH0nj179JOf/MT/+Pnnn9esWbP0zDPPSJJycnL04IMPEqijx5mMBs0blaJ5o1L0U4dLaw9Z9MruU9p4pNYfLh6utOpwpVU/e6dQM4cmavmEDF0zPl1pXC2EAaRz/nRJZ7QWu216tv7w0XFJ0os7ywnU0SdY25zaWFyjdYer9X5RtRrPsrNKiA6gJ5hNRl07MUPXTszQEUuzXtxZrn/tKvd3e6lsatNv3z+q375/VLOHJeq2aTlaNoFpJTBwdJ9DvWf+XYeFmPTCVwv04GsH9I8d5ZKkf+wo156yRj3x+amcrMRFq7a2dQTotdp0tLZbJdwn5WfEanF+qhblp2liVhxTwgHoEQaDQfNH+SoTy+pb9LetpXphe6k/jDlUadV3/rlf//NWoW6blq07Zg9RXjIXa2Jg6rpd7qvn3Q0Gg64Zn66rx6bp1T2n9Kv3ilVW7+t++87BKq05VKWbp2Tp/sWjmApmAPJ4vCqptWnXyUbtLmvQ7tJGFVma5fWe/TWfDNDHZ8YqhDbu0EUE6g0NDUpLS/M/3rBhg5YtW+Z/PGPGDJWVlV3e0QHnERkaohsnZ+nGyVmqtbXrzX2Venn3Ke0pa5Qkeb3S1uP12nq8Xg++dlDThiRo2fh0LZuQoax4KtfRv3VtPZMc3T1Qz8+I1fisWB04ZdW+8iYVVTUzFQKCoqy+Re8dtmjd4WptPV4np/vMPVajQZo2JEFLxqYTogPoFSPTYvS95fn69tLRWl9YrRd3lOn9ohr/xZm+Nm++/cdrJ2ToxsmZmjUsiavQ0a/1dMv3ThGhJv381kmalZek779yQK1Ot4otNt3w249oAY9z8nq9Kqtv1dbjddp+ol7bTzToeG3gLkaS76LiK0cka86IZM0dmdxnT+QDGLhyEiP13WVjdN/ikXpjX6X+svmEv7V0U6tTv990XL/fdFwz8xJ1+4wcLRufQdU6BpTOOdRNRoOSokKDPJpzMxkNumVqtq6bmKkXdpTpN+uOqLq5XV6v9K9dp/T63gpdOyFDyydkaN6oFIWb+Vvtb1xuj47V2HWwokkHK6z+2+Y21zlflxBp1uSceM3MS9LsYYkanxXHPOgI6IKPotPS0nT8+HHl5OTI4XBo165deuihh/zPNzc3MwcVgio5Okx3XjFUd14xVMdr7Xp1zym9tqdCJV0OwHeebNDOkw366ZuHNSknXsvHp2vZ+AzlJhHeoP+pPUeFuiTdNi1HB04dlCS9uKNM379ubK+NDYNXc5tTW0rq9eGRGn14pPasJ0Gjw0I0b1SyFuen6arRqUrs4wdeAAYms8mopePStXRcuqqtbfrX7lP6x44ylXRMQ2Frd+mFHWV6YUeZ0mLDdP3ETN00JUvjMmOZGxP9jq39dGvsqLCeP0G4Ylq2JuXE6f/72y4VW2y0gMcZPB6viqubtb3jIvjtJ+plsZ69Aj3CbNLsYYmaMyJZV45M0ai0aL6LAfQJ4WaTbp2WrVunZWtvWaP+vPmkXt9XIYfLI0nadrxe247X68FXD+qGyZm6fUauxmexP4n+r7NCPSU6rN90hgkNMeqLs4fo1qnZ+vPmE3pqwzE1tjjldHv1yp4KvbKnQpGhJi0ck6pl4zO0YEwKXcv6oDanW4VVzV3Cc6sKK61q7/jePRuT0aAx6TGakhuvKTkJmjokQUOTIvk+xgW54G+C5cuX67vf/a4effRRvfLKK4qMjNSVV17pf37fvn0aPnx4jwwSuFh5yVG6b/EofWPRSBVbbHprf6XePlCpYovNv87eskbtLWvUI28XalxmrJaMTdei/FROkKLf6FqhHihQv3Fyph5+87Acbo9e2XNK31k2hqvrcNm5PV7tK2/smH+qRrtLG+XyBO6blJ0QocX5aVqUn6pZeUkKDeHfI4C+IzU2XP8xf7i+Nm+YdpU26B/by/XGvgrZO+bntVjb/VVGw1KidNPkLN04OZP51tFv2Huh5fsnjUiN0asr5+pHrx3UCzt8He1oAT94NbU4tbe8UXvKfD87TzaoqfXs81WaTQZNyo7X7GFJmjsyWVNzE9h/BNDnTcqJ1//mxOu/r83XSzvL9ML2Mh3ruFizud2lv20t1d+2lio/I1a3z8jRTZOzFBdJkRr6H5fb4y/2SY0987xkXxcRatLX5g/X52bl6vcfHtcfPzrur2Rucbj1xr5KvbGvUuFmo+aPStHyCRlaOCZVMeH8vfamVodbx2psOlLdrGKLTUcsNh2tblZpfYvOcvqxm7TYME3MjteU3HhNzU3QxOw4LpDAp3bB/3J+8pOf6JZbbtH8+fMVHR2tP/3pTwoNPV1N9oc//EFLlizpkUECn5bBYNDo9BiNTo/R/VeP0tFqm945UKm39lfpUKXVv17nVUy/eq9Y6bHhWpifqkVjUjVnRDLtXdBndQvUo8/ccY2PDNXVY9P05v5K1doc+qCoRlePTTtjPeBieL1eldTatbWkXpuO1uijo3VnPREaYjRoam6C5o9O0aL8VI1Oi+GCJQB9nsFg0LQhiZo2JFE/umGc3jts0at7KrShuNo/bUVJjV2/XFusX64t1uSceN04OVPXTshQKu2G0Yf1Vsv3T4oINenRWydq1rBE/ffL3VvAP3TjON06NbvfVDThwrW73Dpc2aw9pQ3aW96kvWWN3brHBRIZatK0IQmaOTRRM/ISNTknnuNxAP1WYlSovjpvuL5ype9izee3lemNfZVqdfou1jxcadWDrx3Uw28d1jXj0nXrtGzNGZHMFEPoN2ptDv881Kkx/fc4KDbcrFVXj9I9C0boo2O1ent/pd49ZFFji+9cV5vTozUHLVpz0KJQk1FzRyarYFiSJmbHaXxWXK/uVw9UHo9X1c3tOlFnV2ldS0eA7gvRyxtazznfeVdDkyI1LjNO47JifbeZsWdMkwpcigv+a09OTtbGjRvV1NSk6OhomUzdD2pefPFFRUdHf+qB/OxnP9MDDzygb3zjG3rsscckSW1tbfrmN7+p559/Xu3t7Vq6dKmefPLJbnO5AxdjRGq07lk4UvcsHKmTdXa9faBKb++v9M9vJElV1jat3lqq1VtLFW42au6IZC0c46uoZE429CXnq1CXpFunZ+vN/ZWSfG3fCdRxsbxer45U27S1pE5bjtdra0l9t+kGPmlYcpSuHOlrwzl7eFKvVcABQE+ICDXp+kmZun5SphrsDr19oEqv7jmlrcfr/et0Vlr++I1Dmj4kQcvGZ+ia8enKjI8I4siBM9mCUKHe1S1TszUxO04r/7ZbRZZmtTrd+q+X9unxdUd0+4wc3TY9h+OtfqrV4VaRpVmHK606VGHVvlNNOlxhlcN97pabCZFmzRiaqJl5iZoxNFHjMmMVQkctAANM14s1f3j9WL2xr1IvbC/TnrJGSZLD5dFreyv02t4KpcaE6YZJmbp5apbGZtBBE31bdXOb/35aP6xQ/6TQEKMWjE7VgtGpetjt0daSer11oFLvHqxSrc0hSXK4PVpfWK31hdWSJINBGpESrYnZ8ZqUE6eJ2fHKz4hRWAgXBH5Su8utysY2X2he36KTdS06WWfXyboWlda3nLdVe1cRZpNGpEZrdHqMxmX6wvP8jBi6B6DHXfRRdFxcXMDliYmJn3oQ27dv1//93/9p4sSJ3Zbff//9evPNN/Xiiy8qLi5O99xzj2655RZ99NFHn/p3AZ2GJEXpP+YP13/MH67KplatO1ytdYct+uhYnX+OozanR+8drtZ7h6ull6UJWXGaPypFV45M1hTazSHIas4zh7okXTkiWakxYapubtf6wmrV2tq5Mg/n5PF4VWRp1taSOm3tmM+y3u446/pxEWbNHZGsK0cma+7IZGUnRPbiaAGg9yREherzs3L1+Vm5qmhs1et7ffPrHe7oeuT1SttPNGj7iQb9+I1DmpwTr+UT0rVsfIZyEvluRPDZHcEN1CVfC/hXVs7p1gK+vKFV/+/dYv3qvSNaOCZVn5+Zq3mjUqjQ66Nqmtt1qCM4P1xp1aFKq0pqbOdtuRlqMmpsZqwm58Rrck68JuXEM18lgEEnJtysz83M1edm5qqoqlkvbC/Ty7vL1dBRCVvdfHqKoVFp0bppSpZumpzFhZrokyzW0+cl+3OFeiDmjkr0uSOT9ZMbx2v7iXq9c6BKbx+o7Pa5vV51VFLb9M9d5R2vNWhMeqzGZ8VpaFKkshMilZ0QoeyECCVGhQ7IfZ92l1uWpnZVNLWqqqlNlU1tqmpqVUVTm//xuYpzziYq1KQRaTEamRrt+0mL1sjUGGXFR9DhCkER9LIxm82mO+64Q88884x++tOf+pc3NTXp2Wef1erVq7Vw4UJJ0h//+Efl5+dry5Ytmj17drCGjAEoIy5CX5g9RF+YPUQtDpc+OlqndYctWldY3a0KeP+pJu0/1aTfvn9UUaEmFQxP1rxRvkpMTgagt3X+2ww3G896UjTEZNQtU7P19IZjcnm8emX3Kf37lcN6c5jo45pandpT1qhdJxu0u6xRe0obZG1znXX96LAQTR+aoNnDkjR7WJImZMVxwhvAoJMZH6GvzR+ur80frmJLs97YW6G3D1TpSLXNv05n5fr/vFWo8VmxWjY+Q8snZCgvmTnXERy2dl+LWYPB11o7WDpbwF8zPl3PfXxCG4/UyOuV3B6v1h6yaO0hizLjwvWZGTn6zPQcQoQg8Hq9qrM7dLTjBPGxapuOVttUZGnudnx8LsNSojQ5O16Tc+M1KTte+RmxXJAOAF2MTo/RD68fq+8sG633C6v1r12n9H7R6SmGii02/fydIv1iTZFm5SXqlinZumZCumKpwEQfMdAq1M/GZDT4z4H98LqxKqxq1r7yRu0tb9K+8kYVVTXL1eXKQqfb688QPinCbFJ2QoRyEk+H7OlxEYqPMCuuy09shDko59q8Xq9anW7Z2l2yt7vV3OZUvd3h/6mzO1Rv67i1t/uXNZ/jPOL5hJqMykmM0JCkKOUmRmpoUqSGJkdpZFqMMuPCyVvQpwQ9UF+5cqWuvfZaLV68uFugvnPnTjmdTi1evNi/bMyYMcrNzdXmzZvPGqi3t7ervf30AZ7V6qsYcTqdcjoDz/EKdGU2SFeNTNRVIxP10HVjdLDSqveLarS+qEYHK5r969kdbr132KL3DlskSdnx4b4r14YnqWBYomIj2MFFz+o8mZUcHSaX6+w7LjdPStfTG45J8rV9/9Ks7H69M9L5Xc53+sXzeLw6VmP3BeflTdpd2qhjtfZzzkUUGx6i6UMSNDMvQTOHJig/PaZbK06P2yWPuxcGDwAXIBjbiLzEcN27YJjuXTBMR6pt/vn1Ci2nw/UDp6w6cMqqX6wp0oiUKC0ck6JFY1I1KZuLktB7bG2+v4vIUNM59x17y9zhCZo7PEHlDa16adcpvbTzlCwd+7cVTW167L0jenzdEc0bmawbJmVo+pAEZcQNrOqnYHO5PTrV1KaTdS06VmPXsRqbjtXYdbTarsbWC/seNZsMGpkarfyMGI1Jj1F+x88Zx8Net5zO4O40chwBoC8ySlo0OlmLRierscWptw5U6bW9ldpZ2ijJVwG7paReW0rq9YNXD2j+qGQtH5+uBaOTFRka9FP7AwrbiYtT2dDiv58YGTJo/ruNTInQyJQIrZiSIUlqd7p1uKpZ+09Ztf9Uk/adsqrkLOfaWp1uf0X7+USHhSguIkSx4WbFRYQoMjREISaDzEaj79bUcWs0KMRkVIjRoBCTQfJKTo9Xbo9XLren+313x32PV3aHLzT3hecu2R1u2dtd5+06dLGMBik1JkzpceFKjw3XkMRI5SZGKDcxUkOSIpUWE3bWavO+cMyCge9ivrsMXu+5TqP3rOeff14PP/ywtm/frvDwcF111VWaPHmyHnvsMa1evVp33XVXt3BckmbOnKkFCxbo0UcfDfieP/rRj/TQQw+dsXz16tWKjKTVIi5Nk0MqbjKosNGgwiaDbM7AX/YGeZUdJY2M82pErFfDY70KZ+oUXEYuj/TNrb4Dp6HRXt0/4dwnp36136QTNt+/129NcCknuseHiCDzeqW6dqnMZlCZ3aAyu+9+q/vcwU2M2au8GN/31ohYrzIjfTu/AICLU90q7a03aG+dUWX2wF+k0SFejUvwanyiV6PjvApjfxE96MGdJjU6DIoze/Xj6X3vaji3VzrcYNDH1QYdajDIqzP/buJDffspnT9ZkRJTbp9bu1uqbZNq2wyqa/fddj5uaJc8Af47n01kiFdZkV5lRUnZUV5lRnqVFiFReA4Al19tm7SjxqAdtUbVtJ35XR1q9O1HTknyKj/eqyA2n8Eg9fwxozZX+3YCvj3RpWwacfm1uaWqFqm+3bf/Vd9uUH1bx2275PIOjBNtESavos1StNm3nx4fKsWHddx2PI4JlUwD4+NigGppadHnP/95NTU1KTY29pzrBu0ytrKyMn3jG9/Q2rVrFR5++a4yf+CBB7Rq1Sr/Y6vVqpycHC1ZsuS8/zGAi+HxeFVoadamo3XadLROO042+NsyedURXtkNWl/haw0zLjNGs/MSNTsvUVNz4xUVpHkLMTBUNrVJWzdKkkbmpGn58snnXL85tVzff/WQJGm3M0OfWzC+33ZRcDqdWrt2ra6++mqZzf3zM1xuXq9XpfWtOlBh1YEKqw52/JyrdbskhRgNys+I8c1lmR2nKblxyo6P6NcdDAAMbn11G1HW0KI1B6v13uFq7Spr9Fcr2FwGba0xaGuNFBpiVMGwRC0ak6IFo1OUHkslLi6v/961XpJLSXFRWr58brCHE9D1HbeVTW16adcpvbjzlG+/t0Ojw6DddQbtrvM9DjcbNTErTlNz4zUlN16j06KVHhs+aDo/tLs8qm7unKeyXVVW3zyVVdZ2VXbMV1lnd1z0+6bFhml4SpSGp0RreEqURnT89Pd5P/vqNgIAzuZL8h3v7y1v0qt7K/X2AYv/e93hOb1NjAo1adGYVC2fkKa5I5IVxpVOnwrbiYvz8l92SdW1kqRbli1ScvTAbft+OXk8vql1yhtbdaqhVdXN7Wpqdamp1ammVqesbU41tbpkbXWqqc0pa6urW0v5yy0q1KSosBD/bXRY5+MQRYWZlBAZqsQos5KiQpXY5Sch0iwzV7ZiAOjscn4hgpbo7dy5U9XV1Zo6dap/mdvt1saNG/Xb3/5Wa9askcPhUGNjo+Lj4/3rWCwWpaenn/V9w8LCFBZ25pe32WxmQ4jLblJukiblJmnlQqnF4dLWknptKK7R5mN1KrKcbg/v9ni1r9yqfeVW/e7DEwoxGjQxO06zhiVpxtAETctNVFwk/z5x4ZraTrdVSo0NP+/3241TsvXTtwrV5vRoXWGNFj+2SfctHqXPz8rttzs/g/V7vbnNqWJLswqrmlVY2ayiqmYdrrJe0HxFKTFhmpITr6lDEjQ1N0ETsuIUwWXsAAagvraNGJYap/9MjdN/LhipOlu71hdWa93ham08UqMWh69S2OHyaENxrTYU10o6rPyMWM0flaKrRqdo2pCEfru9Rt/g9XrV4vDtK0SH962/j0Byk81atWSMvrF4tDYfq9PW43XaebJBe8oa/X8zktTm9GjbiQZtO9HgXxZiNCgjPlw5CZG+n8QIZXfc5iREKjn67G0l+4JWh1u1tnb/3JS1NofqbL77dR1zVtba2mWxtqnWdvFheaeYsBDlJvlabeYmRvmC89RoDU+NHvDz8/a1bQQAnM+MYSmaMSxFP7rBo63H6/XGvgq9c6BKDS2+NrV2h1uv7avUa/sqFRMeoiVj07VsfLrmjkxWuJlj/ovFduLC1HTsh5iMBqXFRfXp/au+JjMsVJmJF9Y+1Lcf71aLwy2XxyOX2yun2yOXp+PW7ZXL42vn7nJ7ZTDI3/49pKM9fOet2WiUqaNNfESoSVGhIfx/w6B3Md/3QQvUFy1apP3793dbdtddd2nMmDH6zne+o5ycHJnNZq1bt04rVqyQJBUVFam0tFQFBQXBGDJwTpGhIVowJlULxqRKkups7dp6vF5bSuq0+Vhdt7lRXB6vdpU2aldpo57qWDYqLVrThiRq+pAEzRiaqJxEqkRxdjW205U6KTHnvwI0Jtys/7l5gr77r/1yuDxqaHHqwdcO6k+bT+iBZflanJ/Kv7c+xuHy6GSdXYVVvtC8sMqqwqpmlTe0XtDrU2PCNCErTuOz4jQhK04TsuOURrUjAARdUnSYbpueo9um56jN6dbmkjq9d8ii9w5bZLGenu7qcKVVhyutenrDMUWHhWjOiCRdNTpV80elKDM+IoifAP1Rq9Ptnw8xqh/Nt2oyGjR3ZLLmjkyW5Jvzu7CqWbtKG7TzpO/nk/tGLo9XZfWtKqtvlVR3xnuGhhiVEGlWbLhZsRFmxYaHKDbCrJjwkC7LzIoOD5HZaJCp44Sk0eA7GWnqWGYyGhTScQLS5fHI4fLK4fbI6fLI6fb47nec8HS6PWp1uNXc5pKt3SVrm1O2jvvN/lunmttcand5Lst/O998leHKSYxQbmKUhvjD80gNSYpSQqSZ/X8A6GdCTEbNGZGsOSOS9eMbx+vjY3V6Y2+F1hys8neoa25z6Z+7yvXPXeWKCjXpqjGpWjouXQtGpyhmgF8whd5V3ew7dknp4xcr9ncGg8FXMU63WyDogvZXGBMTo/Hjx3dbFhUVpaSkJP/yu+++W6tWrVJiYqJiY2N17733qqCgQLNnzw7GkIGLkhQdpuUTMrR8QoYkqaa5XVuP1/kD9mM19m7rF1tsKrbY9PdtpZJ8Ien0IQmaPjRR04YkKD8jRmEhXFUKn5rm0yfcLyRQl6RbpmZr1rAk/eKdQr2yp0KSVFJj11f+vEMFw5L039fma3xWXI+MF2fX1OLU0RqbjnX+VNtVUmPTyfoWuS+wpVN6bLjGZ8WeDs+z4pRKeA4AfV642aQFo1O1YHSqfnrTeB2ssOrdQxZ9UFStfeVN/vVs7S6tOWjRmoMWSb4LMTvD9WlDEqg8wnnZ2k93sunPJ+NCTEaN77hg8EsFQyVJFmubdp1s0O6yRp2ss/vC9IaWs3bvcbg8sljbu13A0t+YjAalxoQpPS5cGXHhyoiLUEZceMdj3/3UmDCF0NkCAAYss8mo+aNSNH9Uih6+eYI2Ha3RG3sr9e4hi3+7b3e49ea+Sr25r1KhJqPmjkzWNePStXhsmhKjQoP8CdCfudwe1dp8+1KpsbR6BzA49Okj6V/96lcyGo1asWKF2tvbtXTpUj355JPBHhbwqaTEhOm6iZm6bmKmJKm6uU07TjT4fk7W62CFtVt4VtPcrrcPVOntA1WSpFCTUfmZsZqSE6/JOfGalBOvoUmRVBUMUt0C9YuYoygrPkKP3T5FX56Tp4ffPKTtHe0xN5fU6frfbtItU7L17aWjlR5HGHs5WducKq1r0cm6FpXWt6i03q5jNb7g/GJadUaFmjQ6PUaj02OVnxGj0WkxGpMey5QRADAAGAwGf1C46upRqrW168MjNfqgqEYbi2v8LT2l0xdi/m5jicJCjJqZl6g5I5I1d0SyxmbEUiGCM9jbT7dJjw4bWBdgpMWGa9mEDC3ruJC5U1OrU2X1LSpvaFF5Q6vK6ltU1tCqisZW3/yUrU7Zu7SPDxajQYoOC1FMeEeVfIRZydGhSooKU1J0qJKiQpUUHdbtNi7CzN85AMAvNMSohWPStHBMmtqcbn18rFbvHKjS2kMW/z6kw+3R+sJqrS+slvFf0sy8RC0dl67F+WnKSYwM8idAf1Nrc8jbcRo7NYZziAAGhz4VqH/wwQfdHoeHh+uJJ57QE088EZwBAT0oNSa8WwV7i8OlPaWN2nGyQdtP1Gt3aWO3ShKH26O9ZY3aW9boXxYfadakbF/APjknXhOz45R0EeEq+q9PU6He1eSceP3jawV650CVfvZOoU7Wtcjrlf65q1xv7q/Qv88dpmUT0jU6LYbKlgvQ5nSrqqlNFU2tKq9v1cl6u0rrW1VaZ9fJ+hY1dglBLkS42ajhKdEalhKt0WnRGp0eqzHpMcqKj+DkKQAMEsnRYbp5SrZunpItt8er/aea9EFRtT4oqtHe8kb/Cax2l0cfHqnVh0dqJUkJkWZd0RGuzx2RzAlSSJLsA6RC/WLERZgV13GRytk43R7Z2nwt2K2tnbdOWTvar7s9Xrm9XrndXrk8Xnm8vlt3lx+P16tQk1HmEKPMJqPCQowymwwym3yPfc8ZFB5iUkxHK/nosBDFhocoOjxEEWYTF0kDAC6bcLPJH6673B5tO1GvNQeqtOagRVVW3/SBHq+0paReW0rq9dDrhzQmPUZXj03Tovw0TcyK47wDzqu6+fRUlGlUqAMYJAbHkTTQD0SGhuiKEcm6YoRvfkC3x6vCKqt2nmzQntJG7SlrVElt9zbxjS1ObSiu0YbiGv+yzLhwjeto+Tw+K1bjM2n9PBDV2C4tUJd8lXDLJmRoYX6q/rL5pB5fd0TWNpfanB799v2j+u37RxUZatKk7HhNyY3X1NwETcmNH3QXbbQ53aq1tctibVNFY5sqm1pV0dimisZWVTb5Hl9MlXlXKTFhGp4SpeEp0b6f1GiNSI1WRmw4B7AAAD+T0eC/gPK+xaPUYHdo45EabTpSq01Ha1XZdPqEVkOL09/aU5KGJEXqiuFJmj0sSbPykuhCM0h1vVA3OpzTAJ3MJqMSokKVQNtbAMAAFGIy6orhybpieLIevH6c9pY36p2DVVpzoEon6lr86xVWNauwqlm/WX9UqTFhWpSfpqvHpuqK4clMLYSAuk6dQ4U6gMGCI2mgjzIZDRqXGadxmXH6UoFvWVOLU3vLfeH63jLfbZ29e5BX0dSmiqY2rT1k8S9LjQnztxAdnxmrsZmxyoqPoBKiH+taoZ58iQF3WIhJ/37lMK2Ymq3H1x/RXzaflKtj+oEWh1ubS+q0uaTOv/6QpEh/uD4uM07pceFKjg5VWEj/OMjyeLxqbnOpsdWhxhan6lscqmluV62tXTXNXX5s7aptbpf1LPNvXgiDQcqMi1BOYoSGJEYpNylSuYmRGpIUqSFJUYqLoFU7AODiJUSF6sbJWbpxcpa8Xq9Kau366GitNh2p1eZjdWruEp6e7Jhy5O/byiT5tuOz8hI1Ky9Js4YlKjuBCvbBoGuFenQopwEAABhsjEaDpuQmaEpugr57zRgVW2x677BFaw9ZtKdLN8zq5nb9fVup/r6tVBFmk+aOTNbi/FRdNTpVaRTsoAMV6gAGI46kgX4kLtKseaNSNG9UiiTJ6/WqvKHVH7DvP9WkgxXWbhUokm9nuHOepE4x4SHK75iHOT8jVmMyYjU6LUYRof0jFB3sOgP1mPCQy3a1cEJUqB68fpzuuiJP7x6q0u6yRu0+2aCKLlVv0ukT8y/vPtVteVyEWakxYUqJCfPf+u6HKy7CrNAQXwtM362py/3TyyTJ4/W1z3R5vPJ0tNnsvHV7vGpzOFVhl3aXNcrhNqjF4VKr0y17u9t33+GW3eG739TqVGOLU42tTjW1OHy3rU5/m9xLZTT45u3MjI9QRpzvNis+QrlJkRqSGKmshIh+c6EBAKB/MhgM/k4nXyoYKpfbo32nmvRRR/X6rtIGOd2nN3yd2/F/7CiXJGXFR2jWsETNzvNVseckctHlQGQbhC3fAQBAYAaDQaPTYzQ6PUYrF4xQtbVN6wur9d5hiz48Uqt2l0eS1Op0a+0hi79oZ2xGrK4anaIFY1I1JSeeKQIHsW4V6gTqAAYJjqSBfsxgMCgnMVI5iZG6flKmJF/17cn6Fl+4fqpJ+0816cCppjOqbJvbXNp2ol7bTtT7lxkN0tDkKH/QPjItRiNTo5WbGMlOch/TGah/2nbv55KbFKl/v3KY/3FVU5t2lzZoV2mDdpX6LtxwdBxcddXUEVYfqbZd9jGdKUTat61Hf0NUqMl/UUDnhQEZceHKiI9QVny4MuIilBoTxt8GAKBPCTEZNTU3QVNzE3TvopFqcbi062Sjth6v09aSeu0pa5TDfXo7fqqxVf/adUr/2uW7UC4jLtxXwT4sSbPyEpWXHEXAPgDY293++9EE6gAAoIvU2HDdPjNXt8/MVavDrU1Ha/XeIYvWFVq6TXF3qNKqQ5VWPfnBMcWGh2jeqBRdNTpV80el9Mj5KfRdNV0q1Gn5DmCw4EgaGGCMRoPykqOUlxylGzpCdq/Xq7L6Vh2o8IXrhyutOlzZrCpr98pjj1cqqbGrpMauN/dX+peHmowalhLlD9hHpkZrZFq0hiRFyUyY2OtaHC7ZHb6Toim9MJ95ely4lk3I0LIJGZIkh8ujQ5VW7TrZoBN1dtU0t6u6o016dXOb2pxnhu3BZjBIseFmxUeaFR9hVlxkqOIjTj/uGpynRIcrOSZUkbRDBQAMAJGhIZo7MllzRyZLktqcbu0qbdDWknptPV6n3aWN/iokSapsatMreyr0yp4KSb6pg2Z2BOyz8xI1IjWagL0fsrU7/fepUAcAAGcTEWrS1WPTdPXYNHk8Xu0pb9QHRTX6oKha+8qb/OtZ21x6Y1+l3tjnO384MTtOV41O1VWjUzQpO14mI/uLA1nXCnWmAgAwWHAkDQwCBoPBN3dzUqSWd4SiktRgd+hwlS9cL6y06nCVVcUW2xnVxw63R4VVzSqsau62PKQjvB+ZFq0RqTH+oD0vOYo21z2otvn01cHBuAI4NMSoyTnxmpwTf8ZzXq9XtnbXJ0L2dtnbXWp3udXu9Mjh9py+dbnlcHnU7vItkyST0SCT0SCj0SCTQTIZjTIZfcuNBoMM8qq2qkKjhw9VdHioIkJNigw1KSo0xH8/ouOxLzAPVUx4iIwczAEAoHCzSVcMT9YVw30Be7vLrb1lTdpaUqetx+u182SDWp2nq5mrm9u7nSxNigrVjKGJmj40QTOGJmpsZiwXWPYDti4V6lFh7KcDAIDzMxoN/s5Hq64epZrmdm0ortH7RdX6sLimWzfMfeVN2lfepMfXHVFCpFnzR/law88bmaKEqNAgfgr0BEtHkZbJaFAS/38BDBIE6sAglhAV2u2EqiS53B4dr7XrcFWzjlbbdLS6WcUWm07U2uXydJ942uXx6ki1raPFd5V/uclo0JCkyI5q9hiNTPPN6zksJYqq38ugxna6s0Bfa6llMBgUE25WTLhZw1Kie+R3OJ1OvfVWuZYvHyOz2dwjvwMAgMEiLMSkmXmJmpmXqHvl60Sz/1STtnQG7Cfq/Z1xJKnO7tA7B6v0zkHfvl+E2aQpufGaPiRB04cmakpuvGLC2T73NfYuc6jT8h0AAHwaKTFhunVatm6dli2X26PdZY16v7Ba7xfV6HCl1b9eQ4vT3/HIYJAm58RrwehULRidqnGZsRQ8DADVnVNRRofx/xPAoMGRNIBuQkxGX2v3tJhuyx0uj07W2X0BusWmI9W+wL2kxt5tHk5Jcnu8/tbxaw5auj2XEReuvOQoDUuJ0rBkX8g+PCVamfERtIO6QJ3zp0t9L1AHAAD9W2iIUdOGJGjakAStXOC72PJAhdVfwb79eL2au4SzrU63Pj5Wp4+P1UmSjAYpPyPWX8U+fUii0uNoAxlsXQN1Wr4DAIBLFWIyasbQRM0Ymqj/umaMqpratKG4Wu8X1mjT0VrZOvY9vF5pd2mjdpc26pdri5UcHaZ5o5I1f1SKrhyZokSqm/sdl9ujWpvv3GRqLOclAQweHEkDuCChIV2C9gmnl7vcHpXWt3QE7c3+wP1Yja3bfJydKpvaVNnU5j/p2vX9hyZF+kP2YSm+1vHDU6IUH8nOdVfdAvVemEMdAAAMXiGm01O9fG3+cLk9XhVbmrXjRL22n2jQzpMNOtXY6l/f45UOVlh1sMKq5z4+IUnKTojo1iZ+REo0lSy9zEaFOgAA6EHpceH67IxcfXZGrhwuj3aebNAHRdV6v6haxRabf71aW7v+teuU/rXrlAwGaWJ2vK4alaL5zL3eb9TZHfJ2NDFNjeHCWQCDB0fSAC5JiMmoYSnRGpYSraXj0v3L3R6vyhtadMRiU3F1s45V23W81qaSWrsaW5xnvI/D5VGxxdZtJ7tTYlSohnVUtef5q9qjlJsYpdCQwTdnJxXqAAAgWExGg/IzYpWfEasvFgyVJJ1qbNWOE/XacaJB20/Uq8jS7D/JJknlDa0qbzill3efkiTFRZj9LeJnDE3Q+Kw4hZuZ17snUaEOAAB6S2iIUQXDk1QwPEkPLM/XqcZWX7heWKOPjtaq1embTsjrlfaWNWpvWaN+ve6I4iLMunKkr3p9/qgUpcYS1vZFnfOnS1IaFeoABhGOpAH0CN886lEakhSlxWPTuj1Xb3eopMYXrvtaw/vun6yzy+n2nvFe9XaH6u0O7TjZcMbvyEmI6GghH+1vIz88JUopMWEyGAbmVa01NgJ1AADQd2TFRyhrcpZunJwlSWpqdWpXaYN2dgTse8oau3Uuamp1al1htdYVVkuSQk1GTcyO8wfs04Yk0KHoMrO3u/33o8K4eAEAAPSerPgI3TFriO6YNUTtLrd2nmjQhuIabSiuUWFVs3+9plan3thXqTf2VUryTSPUGa5PG5IwKItq+iKL9fR5SSrUAQwmBOoAel1iVKgSoxI1fWhit+Uut0enGltVUmPXsY6Q/XiNXSW1tm47a53cHq9O1LXoRF2L3i+q6fZcdFiI8pKjNDItWqPSYjQ6LUaj0mOUGRfe74N2Wr4DAIC+LC7CrAWjU7VgdKokXyeiAxVN/jbxO07Uq6FLxyKH26MdJxu042SDnt7gWzYqLVrThyZq+hBfm/jshIh+vw8XTJ0t30NNRoWFEKgDAIDgCAsx6YoRybpiRLIeWJ6vyqZWbewI1z88UqvmttNddQ5XWnW40qqnNxxTdFiIrhiepPmjUzRvZIpyEiOD+CkGt+pmKtQBDE4E6gD6jBCT0V/VvmBMarfnbO0uf7heUmPvqG636XitXS0O9xnvZWt3af+pJu0/1dRteXRYiEalRWt0eoxGpsZodHqMRqXF9KtK785A3WDwXZwAAADQl4WGGDU1N0FTcxP01XmS1+vVsRq7dp48HbCfqGvp9prOqYBWby2V5DtZN31oomZ0tIrPz4hljs2LYHf4Tk5TnQ4AAPqSjLgI/9zrLrdHe8oa/dXr+8pPn9Oztbv07iGL3j1kkSQNT4nS/FGpmj86RbPyEpk+qBd1q1AnUAcwiBCoA+gXosNCNCE7ThOy47ot93q9sljbVVJj07GOkL2kxq7jtXaVN7TI84kO8rZ2l3aVNmpXaWO35Zlx4ZqcG68pOQmanBuvCX14Ls9am0OSlBQVqhAT7a4AAED/YjAYNCI1WiNSo/XZGbmSfJUuvhbxDdpxsl4HK6xyd9mRs1jb9ea+Sr3Z0QI0OixEU3LjNX2Ir0385Nx4RYZyeHs2nXOoM386AADoq0JMRl+HoqGJ+uaS0aq1tWvTkVptKK7RxuIa1dkd/nWP1dh1rOa4/vDRcUWYTZozIkkLxvg6JGXGRwTxUwx8NV0q1Gn5DmAw4WgaQL9mMBiUHheu9LhwXTEiudtzbU63jtXYVGxpVlGV77bY0qzyhtYz3qeiqU0V+6v01v4qSVKI0aAxGTG+gD0nXlNy45WXHBX0VqNer9dfoZ5Mu3cAADBApMaEa9mEDC2bkCHJFwDvLWv0B+y7TjbI3qUrka3dpQ+P1OrDI7WSJJPRoPGZsZqZl6g5I5I1My+RgL2Lzvap0QTqAACgn0iODtNNU7J005QseTxeHaywakNxtTYU12hXaaP/4stWp1vvHa7We4erJUlj0mO0YEyqFo5J1ZSceIpRLrOuFeppsQTqAAYPjqYBDFjhZpPGZcZpXGb3qnZbu0tHOsL1oiqbDlU2aV95U7fW8S6PVwdOWXXglFV/2XJSkpQQadbVY9N03cRMFQxPkjkIO+TWVpccbo8k9as29QAAABcjKizEP7+mJLncHhVWNfvmYT/ZoO3H61XdfPpkntvj1d7yJu0tb9IzHx6X2WTQlNwEzR2RrDkjkjUpO27Qnkx1uT1qd/n2H6lQBwAA/ZHRaPB3rrxn4Ug1tTr10dFavV9YrfeLalRrO71fWFjVrMKqZj31wTHFRZg1f1SKFuWn6qrRqYqLMAfxUwwMFquvQt1kNCiJqSgBDCIcTQMYdHwtQhM0JTfBv8zt8epIdbN2lzZqT2mjdpc16Ei1Td4uLeMbWpz6x45y/WNHuRIizbpmfIaum5ihWXmJvXaCtsZ2uq0SgToAABgsQkxGjc+K0/isOH15Tp68Xq/KG1q1/cTpediPVNv86zvdXm07Xq9tx+v1y7XFigkL0axhSZozIklzRyRrRGp00DsP9RZ7++mLRgnUAQDAQBAXYdbyCRlaPiHDX72+vrBa64uqta+80X8+r6nVqdf2Vui1vRUKMRpUMDxJS8am6eqx6UqPo7r60+i8qDUlOkxG4+DYnwYAiUAdACT5rqockx6rMemx+txM31yezW1O7Stv0p6yRu0ubdDmY3X+VqMNLU79fVup/r6tVMnRobpmfLqum5ipGUMTZerBncmulVgE6gAAYLAyGAzKSYxUTmKkbpmaLUlqsDu0uaROm47W6uOjtTpR1+Jfv7ndpfcOW/TeYYskKSs+QlePTdPVY9M0My8xKJ2HeovN4fLfjw4zBXEkAAAAl1/X6vVvLB6pWlu7Piiq0fuF1dpYXKPmdt++kMvj9U8Z9INXD2pSTryWjE3T0nFpGp4SHeRP0T+43B5/N4DUWM5LAhhcCNQB4Cxiws2a09EmVPLNyf5BUbVe31ep9Yer1er0heu1Nof+uqVUf91SqtSYMN0wKVP3LhrZI22karoG6syhDgAA4JcQFeqvVJKksvoWfXysVpuO1unjo7Wqszv8655qbNVzH5/Qcx+fUGx4iBaOSdWScemaNyplwM0zbm8/HahHMa88AAAY4JKjw3TrtGzdOi1bTrdHO040aO0hi9YcrNKpxlb/envLGrW3rFG/WFOkYclRWjQmRTHNkrdru0p0U2d3+Kv/U2Oo8AcwuHA0DQAXKNxs0jXjM3TN+Ay1OFxaX1itN/dVan1htX9eyurmdv1+03GtOVSl335uqiblxF/WMdRQoQ4AAHBBchIj9dnEXH12Rq48Hq8Kq5r10dFabTxSoy0ldXK6fWcDrW0uvbKnQq/sqVCoyag5I5J09dh0LR6bOiBOFNq6BOrR4ZwCAAAAg4fZZFTB8CQVDE/SD67L16FKq949aNG7hyw6XGn1r1dSa1fJJrukEL1UsUk3Tc7SDZMzNSI1JniD74M650+XpDQq1AEMMhxNA8CnEBkaousmZuq6iZmytbu07rBFb+yr1IaiGjncHpXVt+rWpz/Wfy/P151XDL1sc3TW2KhQBwAAuFhGo0FjM2M1NjNWX5k3TNY2pz4oqtHaQxZ9UFjtbwXqcHv0flGN3i+q0X+/Il0xPEmfm5mrJWPTFRrSP9vCd61QH2jV9wAAABfKYDBoXGacxmXG6f6rR6msvkXvdlSu7zhRL09H5XVpfaseX39Uj68/qrEZsbpxcqaun5SpzPiI4H6APqDaevq85EC48BQALgZH0wBwiaLDQnTj5CzdODlLpxpbde/qXdpV2iin26sfvX5IW4/X69FbJyo2/NJbwFOhDgAAcOliw826YVKmbpiUKYfLoy0ldVp7yKK1hyyq6qi88Xqlj47W6aOjdUqKCtWt07J1+8xc5SVHBXn0F6dby3cCdQAAAEm+bkZ3z83T3XPzVGdr15oDFXru/QM6ajX6w/VDlVYdqrTqkbcLNTMvUTdOztTy8RlKiAoN7uCDxNJMhTqAwYujaQC4jLLiI/TC1wr0izVF+t3GEknS2weqdLDCqic+P1UTsuMu6f0J1AEAAC6v0BCj5o1K0bxRKfrxjeO0/1ST1h6y6LW9FTpZ1yLJN1/k/20s0f9tLFHBsCR9blaulo5LU1iIKcijPz9bu9t/n0AdAADgTEnRYbptWraiLPs0/coFWnOoRq/urdDeskb/OtuO12vb8Xo9+OpBLcpP1RdnD9WcEUmXrStlf2DpWqFOoA5gkOFoGgAuM7PJqO8tz9fMoYn65ot71dTqVGl9i1Y89bG+f12+vjh7yKfe2a61OTp+h0FxEZde8Q4AAIDTDAaDJmbHa2J2vO5fPEpbSuq0elup1hys8s+5vrmkTptL6pQQafZXrQ9PiQ7yyM+ue8v3vn8BAAAAQDClxoTp3+bm6d/m5ulErV2v7a3Qq3tO6ViNXZLk8ni15qBFaw5aNCwlSl+cPUQrpmVfls6UfV1Nlwp1Wr4DGGz65yRwANAPLB6bpje/PleTc+Il+ebk/OGrB3XP6t2ytjk/1Xt2VqinRIcNqitgAQAAepvRaNAVI5L1289P1ZYHFul7y8d0a/fe0OLUMx8e16L/3aCVq3epvKEliKM9O1vXlu+hXFMPAABwoYYmR+nri0bqvVXz9ebX5+pr84YptUvHyJIaux56/ZBmPbxOD/xrvw5XWoM42p7XtUI9LZZAHcDgQqAOAD0oOyFS//hagf59bp5/2Zv7K3X9bzbpwKmmi3ovt8erentHoE67dwAAgF6TFB2mr84brvXfnK+/f2W2bpiUqVDT6cPpN/dVatH/btCv1har1eE+xzv1Plu3CnUCdQAAgItlMBg0LjNODyzP10ffXagnPj9Vs/IS/c+3Ot36+7ZSLfv1h7rt6Y/12t4KOVyeII64Z1R3VKibjAYlDdJ55AEMXhxNA0APCw0x6vvXjdXMvER968W9sra5dLKuRZ/5v81au2q+suIjLuh96uzt8vg6jRKoAwAABIHBYFDB8CQVDE9Svd2hf+4s11Mbjqne7lC7y6Nfrzuil3aW64HlY3TthIw+0VGoa8t35lAHAAC4NGaTUddOzNC1EzNUbGnWXzaf1L92lcvecVHl9hMN2n6iQcnRYbp7bp7umjNU4eaBMe1OZ4V6SnSYjMbg7+cCQG+iQh0AesmScel68+tXalJ2nCSpxeHWOweqLvj1ne3eJQJ1AACAYEuMCtVX5g3T+9+6SnfPzVNIx0nFU42tumf1bn32d1t0sOLiOhL1BBuBOgAAQI8YlRajn9w0Xlu+t0g/vnGcRqRG+5+rtbXr0XcKNf8X7+vv20rlcvfvinWX26Nam+/cZGos5yUBDD4E6gDQi3ISI/XwzRP8j/eVN17wa7sG6snR7LgCAAD0BXERZv3gurF6574rNW9Uin/5tuP1uv43m/S9l/erztZ+jnfoWXZavgMAAPSomHCzvlQwVGvvn6e/f2W2lk9IV2cBt8Xargf+tV9LHtuodw5Uyuv1Bnewn1Kd3aHOoafGMH86gMGHQB0Aetno9BiFhvi+fveVX3jVEhXqAAAAfdeI1Bj96a4ZevbO6RqaFClJ8nil1VtLteD/faA/bDout6f3T6Da20/P6R4dTqAOAADQUzqnB3ryjml65755Wpyf5n+upMau//jrLt385MfafKwuiKP8dCzWNv/9NCrUAQxCBOoA0MvMJqPGZsRKko7X2tXU6ryg19V0qWxKoUIdAACgzzEYDFqUn6Y198/TA8vG+CvCrW0u/fiNQ/rG87t7vd1n15bvkQNk/k4AAIC+blRajH5/53S99B8FmjE0wb98T1mjPvfMFt35h219YnqgC1VtPX1ekgp1AIMRgToABEHnPOqStP8Cq9SpUAcAAOgfwkJM+tr84Vr/rfm6bVq2f/kb+yr1Xy/t69VK9c6W71GhJhk7e48CAACgV0wfmqh/fK1Az945XaPTYvzLNxTX6NrHN+kbz+/uVv3dV1maqVAHMLgRqANAEEzMjvff33eq8YJeU2tz+O8TqAMAAPR9qTHh+sVtk/TsndNlNvnC7H/tPqUH/rVPnl4K1f2BOvOnAwAABEVnF6O3vnGl/ve2ScqKj/A/9+qeCi3/9Yf68EhNEEd4fpauFeoE6gAGIQJ1AAiCSTmnK9T3lV1ohfrpK0GTafkOAADQbyzKT9MTn5+qkI4K8X/sKNf3Xz0gr7fnQ/XOlu/RBOoAAABBZTIatGJattZ/a75+cN1YJUSaJUl1doe+9Idt+tXa4l7tZHQxup6XpOU7gMGIQB0AgiAvOVpRob45LPeVN17QazpbvkeFmqgwAgAA6GeWjEvX45+bIlNHqL56a6keev1Qj4bqXq/XH6iz/wgAANA3hIWYdPfcPK3/5lVaMDpFkuT1Sr9ed0R3/mGbam3t53mH3te1Qj0tlkAdwOBDoA4AQWAyGjQ+y1elXtHU1m1+9LPpXId27wAAAP3T8gkZ+uVnJqlzKvPnPj6hh9883GOhepvTo84ip6gwU4/8DgAAAHw6CVGhevbOGfr20tH+/cNNR2t17eMfatvx+uAO7hOqOyrUTUaDkqJCgzwaAOh9BOoAECSTcuL99/efZx71Nqdb1jZfdRGBOgAAQP914+Qs/eLWSTJ0nDT9/abjevSdoh4J1Tur0yVavgMAAPRFRqNBKxeM0OqvzPaf87NY2/W5Z7bo6Q3HemWKoAvRWaGeEh0mY2f6DwCDCIE6AATJxOzT86jvPc886l1bPTF/OgAAQP+2Ylq2fnbLBP/jpzcc06/eO3LZf4+9S6BOy3cAAIC+a/awJL319St1xfAkSZLb49XP3i7UV/68Q00tzqCOzeX2qK7j3GRqLOclAQxOBOoAECQTs+L99883j3rXlvBUqAMAAPR/n52Rq5/eNN7/+PF1R/SbdZc3VLcRqAMAAPQbKTFh+svds3TvwhH+Ze8drta1v/lQe8sagzauOrvDP41QagzzpwMYnAjUASBIchIjlBBpliTtK286ZwunboE6FeoAAAADwhdmD9GD14/1P/7ftcV6esOxy/b+XSvUYwjUAQAA+jyT0aBvLhmt5+6a4T9vWN7Qqluf/ljrDluCMiaLtc1/P40KdQCDFIE6AASJwWDQhOx4Sb4rPSua2s66bo2NCnUAAICB6K45efrv5fn+xz97u1DvF1Vflve2O6hQBwAA6I+uGp2qN79+pabmxkuSnG6vvvH8Hh2tbu71sVRbT5+XpEIdwGBFoA4AQTSpyzzq+87Ruqm22eG/T6AOAAAwsHxl3jB9a8ko/+PnPjpxWd7X1u723ydQBwAA6F8y4yP0wtcKdO2EDEm+6Xy+8uedamrt3TnVLc1UqAMAgToABNHEjgp1Sdpb3nTW9Wpsp3dcCdQBAAAGnv/vqhHKToiQJG08UqOKxtZLfs+uLd+jw0yX/H4AAADoXWaTUb+4baLyM2IlScdr7brv+d1ye84+deTl1q1CnUAdwCBFoA4AQTSxa4V6eeNZ1+s2hzqBOgAAwIBjNBr0mek5kiSvV3pxR/klv2fXQJ0KdQAAgP4pMjREv/viNP+c6u8X1eiXa4t67fdXd6lQp+U7gMGKQB0AgigtNtzfKmn/qSZ5znJ1addAPSmKQB0AAGAgunVatgwG3/0Xd5addd/wQjW3EagDAAAMBDmJkfrt56fKZPTtLD7x/jG9ua+yV363pUuFelosgTqAwYlAHQCCrLPte3ObSyfq7AHXqbH5dlzjI80KDeGrGwAAYCDKjI/QvJEpkqTyhlZ9fKzukt6ve8t3AnUAAID+bM6IZH1veb7/8bde3KvDldYe/72dFeomo0FJUaE9/vsAoC8ilQGAIJvUre37mfOoe71ef4V6SjTV6QAAAAPZ7TNy/Pef3156Se9ld3SpUA8lUAcAAOjv/m3OUN0yNUuS1Op066t/2aEGu6NHf2dnhXpKdJiMHRXyADDYEKgDQJBN6KhQl6S9AeZRt7W71Ob0SGL+dAAAgIFuUX6av/Ln3YOWSzpBamt3++9ToQ4AAND/GQwG/c/NEzSxo0CnrL5V9/59t1xuT4/8Ppfbo7qOzpmpsZyXBDB4EagDQJBNzDp3hXrX+dMJ1AEAAAa20BCjbp7iqzpyuD16Zc+pT/1e3Vq+hxOoAwAADAThZpOe/sI0JUf7LsLcdLRWj75T2CO/q87ukMfru58aw/zpAAYvAnUACLKEqFDlJkZKkg5WNJ1xRWm3QJ2W7wAAAAPeZ7u0fX9he5m8Xu+neh9bl0A9Ksx0yeMCAABA35AZH6En75imkI4W7M98eFyv7P70F2KejcXa5r+fRoU6gEGMQB0A+oDONk1tTo+OVNu6PVdrO93mkwp1AACAgW9kWoym5sZLkgqrmgN2MboQnRXqZpNBYSEE6gAAAAPJzLxEPXjDOP/j7/xznw6c+nT7jWdTbT1d6EOFOoDBjEAdAPqAzkBdkvZ9Yh71mubTV4ISqAMAAAwO3arUd5R9qvfoDNSjmD8dAABgQPrCrFzd3rHf2O7y6Kt/3qEGu+M8r7pwlmYq1AFAIlAHgD5hYna8//7eT1Qg1diYQx0AAGCwuW5ipqJCfVXlr+2pUIvDdZ5XnMnW7pYkRYUSqAMAAAxEBoNBD904zt/dqKKpTX/8+MRle/9uFeoE6gAGMQJ1AOgDxmfFyeCb8kj7PxmoNxOoAwAADDZRYSG6bmKmJN9c6G/tr7ro9+isUI+mQh0AAGDACgsx6fHPTZGpYz7157eVyun2XJb3ru5SoU7LdwCDGYE6APQB0WEhGpESLUkqrLKq3eX2P9c1UE+OJlAHAAAYLD7Tpe37P7ZfXNt3l9ujVmdHhXoY86cDAAAMZNkJkVqcnypJqm5u19pDlsvyvl0r1NNiCdQBDF4E6gDQR0zomEfd6fbqcGWzf3lny3eT0aCEyNCgjA0AAAC9b2puvEak+i663HaiXsdqbBf8Wrvj9AWazKEOAAAw8H1x9lD//b9sPnlZ3rNzDnWT0aCkKM5LAhi8CNQBoI+Y1GUe9X3ljf77nRXqSVGh/tZNAAAAGPgMBoNu71qlvuPCq9Q7271LtHwHAAAYDK4YnqRhyVGSpM0ldTpa3XyeV5yfpaNCPSU6TEbOSwIYxAjUAaCPmNhRoS5J+zrmUfd4vKq1OSQxfzoAAMBgdPOULJlNvpOX/9x56oLnw+waqFOhDgAAMPAZjQbdMXuI//Fft5Re0vu53B7VdXTOTI3lvCSAwS2ogfpTTz2liRMnKjY2VrGxsSooKNDbb7/tf76trU0rV65UUlKSoqOjtWLFClksl2fuDwDoa/IzYhXScaVnZ4V6Q4tDbo9XEoE6AADAYJQUHabF+WmSpFpbu9YXVl/Q62xUqAMAAAw6t07NVrjZF/v8c2d5t4ssL1ad3aGO05JKjWH+dACDW1AD9ezsbP3sZz/Tzp07tWPHDi1cuFA33nijDh48KEm6//779frrr+vFF1/Uhg0bVFFRoVtuuSWYQwaAHhNuNmlMRowk6Wi1TfZ2l786XfK1VgIAAMDg89mubd+3X1jbd3v76TnUCdQBAAAGh7hIs26YlClJam536bW9FZ/6vao72r1LUhoV6gAGuaAeVV9//fXdHj/88MN66qmntGXLFmVnZ+vZZ5/V6tWrtXDhQknSH//4R+Xn52vLli2aPXt2wPdsb29Xe/vpL3qr1SpJcjqdcjqdPfRJAODyGJcRqwOnrPJ4pT2ldXK6vf7nkqLMfI9J/v8G/LcAAHwS2wgMVLOHxisjLlyVTW16v6haZXXNSo89d5VQU0ub/354iIG/Cwx6bCMAAOcykLYTn5uerX/sKJck/fnjE1oxOV0Gw8XPf36qwea/z3lJAAPRxXyv9ZnL1N1ut1588UXZ7XYVFBRo586dcjqdWrx4sX+dMWPGKDc3V5s3bz5roP7II4/ooYceOmP5u+++q8jIyB4bPwBcDoZ6gySTJOnF97Yq2iz/4+rSo3rrrSNBG1tfs3bt2mAPAQDQR7GNwEA0MdqoyiajPF7pZ8+/ryXZ3nOuv63m9H7liaOFest2uBdGCfR9bCMAAOcyULYTQ6JNOmkz6HBVs578x9vKi7n49/jIcnp/svpEsd56q+jyDhIAgqylpeWC1w16oL5//34VFBSora1N0dHRevnllzV27Fjt2bNHoaGhio+P77Z+Wlqaqqqqzvp+DzzwgFatWuV/bLValZOToyVLlig2NranPgYAXBZ5lc16/snNkiRnbJayMmOlo8WSpPkzp2j5hPRgDq9PcDqdWrt2ra6++mqZzeZgDwcA0IewjcBANrGhVe/+6kN5vdJ+e7R+ec1cGY1nrzSq31oqHS2UJM2aOknLJ2f21lCBPoltBADgXAbadqIt45S+8y/f1LrHTTlauXzCRb/H0fVHpZISSdLCK6ZrweiUyzpGAAi2zi7nFyLogfro0aO1Z88eNTU16aWXXtKdd96pDRs2fOr3CwsLU1jYmfN5mM3mAbEhBDCwjc2KV1iIUe0ujw5UWJUZH+F/Li0+ku+xLvheBwCcDdsIDER5qWbNGZ6sTUdrVVrfql3lzSoYnnTW9VtdpyvYYyPD+JsAOrCNAACcy0DZTtw4JUePvFOsxhan3j5g0Q+vH6ek6IubB73WfroVcmZC1ID47wIAXV3M95qxB8dxQUJDQzVixAhNmzZNjzzyiCZNmqRf//rXSk9Pl8PhUGNjY7f1LRaL0tOp0AQwMIWYjBqX6eumcbKuRUeqT89VlBJzcTu9AAAAGFg+MyPHf/+F7aXnXNfW5vLfjw4L+rX0AAAA6EXhZpM+M9237+hwe/xzql+Mamu7/35abPhlGxsA9Ed97qja4/Govb1d06ZNk9ls1rp167RixQpJUlFRkUpLS1VQUBDkUQJAz5mYHa9dpY2SpM3H6vzLCdQBAAAGtyVj0xQXYVZTq1Ov76vU/lNNSogMVXxkqBIizUqIClV8pFkJkaE6VHm6dV0UgToAAMCg8/mZufrdRl/L9tXbTuqr84bJdI4pgz7J0twmSTIZDUqKCu2RMQJAfxHUo+oHHnhAy5YtU25urpqbm7V69Wp98MEHWrNmjeLi4nT33Xdr1apVSkxMVGxsrO69914VFBRo9uzZwRw2APSoSTlx/vvtLo8kKSzEqBhOhAIAAAxq4WaTbp6Spec+PiG3x6tjNXZJ9vO+LjrM1PODAwAAQJ8yNDlK80alaGNxjcrqW7WxuEYLxqRe8Os7K9RTosNkvIggHgAGoqCmM9XV1frSl76kyspKxcXFaeLEiVqzZo2uvvpqSdKvfvUrGY1GrVixQu3t7Vq6dKmefPLJYA4ZAHrcxOz4M5alxITJYGDHFQAAYLC7d+EIlda3qNjSrMYWp2ztrnOuHxMeooy4iF4aHQAAAPqSL84eoo3FNZKkv2w5ecGBusvtUa3NF6inxtI1EwCCGqg/++yz53w+PDxcTzzxhJ544oleGhEABF9eUpRiwkLU3OXkKO3eAQAAIElJ0WH6w5dn+B87XB41tjrU2OJUg92hhhanGlt8t/Z2lxbmp9LyHQAAYJBaOCZVWfEROtXYqveLqlVW36KcxMjzvq7O7pDH67ufGsP86QDAUTUA9DFGo0Hjs+K0uaTL/OnRBOoAAAA4U2iIUakx4ZzoBAAAwBlMRoM+PytXv1hTJK9X+tvWUn132Zjzvq6z3bskpVGhDgAyBnsAAIAzTewyj7pEhToAAAAAAACAi/eZ6Tkym3xTSf5jR5naXe7zvsZibfPf58JNACBQB4A+adIn5lFPpkIdAAAAAAAAwEVKiQnTNeMzJEn1dofe3l913tdUN1OhDgBdEagDQB80MZsKdQAAAAAAAACX7ouzh/jv/2XLyfOu361CnUAdAAjUAaAvyoqPUGJUqP8xgToAAAAAAACAT2PG0ASNTouRJO082aCDFU1nrNPc5tSu0ga9sL1UHxRV+5fT8h0ApJBgDwAAcCaDwaCJ2XH6oKhGEoE6AAAAAAAAgE/HYDDoCwVD9INXDkiSnnz/mK4anaIj1TYVW5pVXNWsiqa2gK9NjyNQBwACdQDoo26anKUPimqUFhumsRmxwR4OAAAAAAAAgH7q5ilZ+tlbh2V3uPXm/kq9ub/yvK+5ZUqWkqMp9AEAAnUA6KNumpKlCdlxSo0JU7jZFOzhAAAAAAAAAOinosNCtGJatv68+cw51GPCQzQqLUaj0qI1MjXGfz81lup0AJAI1AGgTxueEh3sIQAAAAAAAAAYAL69dLQMklqdbo1Ki9HItBiNTotRWmyYDAZDsIcHAH0WgToAAAAAAAAAAMAAFxNu1kM3jg/2MACg3zEGewAAAAAAAAAAAAAAAPRFBOoAAAAAAAAAAAAAAARAoA4AAAAAAAAAAAAAQAAE6gAAAAAAAAAAAAAABECgDgAAAAAAAAAAAABAAATqAAAAAAAAAAAAAAAEQKAOAAAAAAAAAAAAAEAABOoAAAAAAAAAAAAAAARAoA4AAAAAAAAAAAAAQAAE6gAAAAAAAAAAAAAABECgDgAAAAAAAAAAAABAAATqAAAAAAAAAAAAAAAEQKAOAAAAAAAAAAAAAEAABOoAAAAAAAAAAAAAAARAoA4AAAAAAAAAAAAAQAAE6gAAAAAAAAAAAAAABECgDgAAAAAAAAAAAABAAATqAAAAAAAAAAAAAAAEQKAOAAAAAAAAAAAAAEAABOoAAAAAAAAAAAAAAARAoA4AAAAAAAAAAAAAQAAE6gAAAAAAAAAAAAAABECgDgAAAAAAAAAAAABAAATqAAAAAAAAAAAAAAAEQKAOAAAAAAAAAAAAAEAABOoAAAAAAAAAAAAAAARAoA4AAAAAAAAAAAAAQAAE6gAAAAAAAAAAAAAABECgDgAAAAAAAAAAAABAAATqAAAAAAAAAAAAAAAEQKAOAAAAAAAAAAAAAEAABOoAAAAAAAAAAAAAAARAoA4AAAAAAAAAAAAAQAAE6gAAAAAAAAAAAAAABECgDgAAAAAAAAAAAABAAATqAAAAAAAAAAAAAAAEQKAOAAAAAAAAAAAAAEAABOoAAAAAAAAAAAAAAARAoA4AAAAAAAAAAAAAQAAE6gAAAAAAAAAAAAAABECgDgAAAAAAAAAAAABAAATqAAAAAAAAAAAAAAAEQKAOAAAAAAAAAAAAAEAABOoAAAAAAAAAAAAAAARAoA4AAAAAAAAAAAAAQAAE6gAAAAAAAAAAAAAABECgDgAAAAAAAAAAAABAAEEN1B955BHNmDFDMTExSk1N1U033aSioqJu67S1tWnlypVKSkpSdHS0VqxYIYvFEqQRAwAAAAAAAAAAAAAGi6AG6hs2bNDKlSu1ZcsWrV27Vk6nU0uWLJHdbvevc//99+v111/Xiy++qA0bNqiiokK33HJLEEcNAAAAAAAAAAAAABgMQoL5y995551uj5977jmlpqZq586dmjdvnpqamvTss89q9erVWrhwoSTpj3/8o/Lz87VlyxbNnj07GMMGAAAAAAAAAAAAAAwCQQ3UP6mpqUmSlJiYKEnauXOnnE6nFi9e7F9nzJgxys3N1ebNmwMG6u3t7Wpvb/c/tlqtkiSn0ymn09mTwwcA9ILO73K+0wEAn8Q2AgBwNmwjAADnwnYCAAafi/nO7zOBusfj0X333ac5c+Zo/PjxkqSqqiqFhoYqPj6+27ppaWmqqqoK+D6PPPKIHnrooTOWv/vuu4qMjLzs4wYABMfatWuDPQQAQB/FNgIAcDZsIwAA58J2AgAGj5aWlgtet88E6itXrtSBAwe0adOmS3qfBx54QKtWrfI/tlqtysnJ0ZIlSxQbG3upwwQABJnT6dTatWt19dVXy2w2B3s4AIA+hG0EAOBs2EYAAM6F7QQADD6dXc4vRJ8I1O+55x698cYb2rhxo7Kzs/3L09PT5XA41NjY2K1K3WKxKD09PeB7hYWFKSws7IzlZrOZDSEADCB8rwMAzoZtBADgbNhGAADOhe0EAAweF/N9b+zBcZyX1+vVPffco5dfflnr169XXl5et+enTZsms9msdevW+ZcVFRWptLRUBQUFvT1cAAAAAAAAAAAAAMAgEtQK9ZUrV2r16tV69dVXFRMT458XPS4uThEREYqLi9Pdd9+tVatWKTExUbGxsbr33ntVUFCg2bNnB3PoAAAAAAAAAAAAAIABLqiB+lNPPSVJuuqqq7ot/+Mf/6gvf/nLkqRf/epXMhqNWrFihdrb27V06VI9+eSTvTxSAAAAAAAAAAAAAMBgE9RA3ev1nned8PBwPfHEE3riiSd6YUQAAAAAAAAAAAAAAPgEdQ51AAAAAAAAAAAAAAD6KgJ1AAAAAAAAAAAAAAACIFAHAAAAAAAAAAAAACAAAnUAAAAAAAAAAAAAAAIgUAcAAAAAAAAAAAAAIAACdQAAAAAAAAAAAAAAAiBQBwAAAAAAAAAAAAAgAAJ1AAAAAAAAAAAAAAACIFAHAAAAAAAAAAAAACAAAnUAAAAAAAAAAAAAAAIgUAcAAAAAAAAAAAAAIAACdQAAAAAAAAAAAAAAAiBQBwAAAAAAAAAAAAAgAAJ1AAAAAAAAAAAAAAACIFAHAAAAAAAAAAAAACAAAnUAAAAAAAAAAAAAAAIgUAcAAAAAAAAAAAAAIAACdQAAAAAAAAAAAAAAAiBQBwAAAAAAAAAAAAAgAAJ1AAAAAAAAAAAAAAACIFAHAAAAAAAAAAAAACAAAnUAAAAAAAAAAAAAAAIgUAcAAAAAAAAAAAAAIAACdQAAAAAAAAAAAAAAAiBQBwAAAAAAAAAAAAAgAAJ1AAAAAAAAAAAAAAACIFAHAAAAAAAAAAAAACAAAnUAAAAAAAAAAAAAAAIgUAcAAAAAAAAAAAAAIAACdQAAAAAAAAAAAAAAAiBQBwAAAAAAAAAAAAAgAAJ1AAAAAAAAAAAAAAACIFAHAAAAAAAAAAAAACAAAnUAAAAAAAAAAAAAAAIgUAcAAAAAAAAAAAAAIAACdQAAAAAAAAAAAAAAAiBQBwAAAAAAAAAAAAAgAAJ1AAAAAAAAAAAAAAACIFAHAAAAAAAAAAAAACAAAnUAAAAAAAAAAAAAAAIgUAcAAAAAAAAAAAAAIAACdQAAAAAAAAAAAAAAAiBQBwAAAAAAAAAAAAAgAAJ1AAAAAAAAAAAAAAACIFAHAAAAAAAAAAAAACAAAnUAAAAAAAAAAAAAAAIgUAcAAAAAAAAAAAAAIAACdQAAAAAAAAAAAAAAAiBQBwAAAAAAAAAAAAAgAAJ1AAAAAAAAAAAAAAACIFAHAAAAAAAAAAAAACAAAnUAAAAAAAAAAAAAAAIgUAcAAAAAAAAAAAAAIAACdQAAAAAAAAAAAAAAAiBQBwAAAAAAAAAAAAAgAAJ1AAAAAAAAAAAAAAACIFAHAAAAAAAAAAAAACAAAnUAAAAAAAAAAAAAAAIIaqC+ceNGXX/99crMzJTBYNArr7zS7Xmv16sf/vCHysjIUEREhBYvXqwjR44EZ7AAAAAAAAAAAAAAgEElqIG63W7XpEmT9MQTTwR8/uc//7kef/xxPf3009q6dauioqK0dOlStbW19fJIAQAAAAAAAAAAAACDTUgwf/myZcu0bNmygM95vV499thj+v73v68bb7xRkvTnP/9ZaWlpeuWVV3T77bcHfF17e7va29v9j61WqyTJ6XTK6XRe5k8AAOhtnd/lfKcDAD6JbQQA4GzYRgAAzoXtBAAMPhfznR/UQP1cjh8/rqqqKi1evNi/LC4uTrNmzdLmzZvPGqg/8sgjeuihh85Y/u677yoyMrLHxgsA6F1r164N9hAAAH0U2wgAwNmwjQAAnAvbCQAYPFpaWi543T4bqFdVVUmS0tLSui1PS0vzPxfIAw88oFWrVvkfW61W5eTkaMmSJYqNje2ZwQIAeo3T6dTatWt19dVXy2w2B3s4AIA+hG0EAOBs2EYAAM6F7QQADD6dXc4vRJ8N1D+tsLAwhYWFnbHcbDazIQSAAYTvdQDA2bCNAACcDdsIAMC5sJ0AgMHjYr7vjT04jkuSnp4uSbJYLN2WWywW/3MAAAAAAAAAAAAAAPSUPhuo5+XlKT09XevWrfMvs1qt2rp1qwoKCoI4MgAAAAAAAAAAAADAYBDUlu82m01Hjx71Pz5+/Lj27NmjxMRE5ebm6r777tNPf/pTjRw5Unl5efrBD36gzMxM3XTTTcEbNAAAAAAAAAAAAABgUAhqoL5jxw4tWLDA/3jVqlWSpDvvvFPPPfec/uu//kt2u11f/epX1djYqLlz5+qdd95ReHh4sIYMAAAAAAAAAAAAABgkghqoX3XVVfJ6vWd93mAw6Mc//rF+/OMf9+KoAAAAAAAAAAAAAADow3OoAwAAAAAAAAAAAAAQTATqAAAAAAAAAAAAAAAEQKAOAAAAAAAAAAAAAEAABOoAAAAAAAAAAAAAAARAoA4AAAAAAAAAAAAAQAAE6gAAAAAAAAAAAAAABECgDgAAAAAAAAAAAABAAATqAAAAAAAAAAAAAAAEQKAOAAAAAAAAAAAAAEAABOoAAAAAAAAAAAAAAARAoA4AAAAAAAAAAAAAQAAE6gAAAAAAAAAAAAAABECgDgAAAAAAAAAAAABAAATqAAAAAAAAAAAAAAAEQKAOAAAAAAAAAAAAAEAABOoAAAAAAAAAAAAAAARAoA4AAAAAAAAAAAAAQAAE6gAAAAAAAAAAAAAABECgDgAAAAAAAAAAAABAAATqAAAAAAAAAAAAAAAEQKAOAAAAAAAAAAAAAEAABOoAAAAAAAAAAAAAAARAoA4AAAAAAAAAAAAAQAAE6gAAAAAAAAAAAAAABECgDgAAAAAAAAAAAABAAATqAAAAAAAAAAAAAAAEQKAOAAAAAAAAAAAAAEAABOoAAAAAAAAAAAAAAARAoA4AAAAAAAAAAAAAQAAE6gAAAAAAAAAAAAAABECgDgAAAAAAAAAAAABAAATqAAAAAAAAAAAAAAAEQKAOAAAAAAAAAAAAAEAABOoAAAAAAAAAAAAAAARAoA4AAAAAAAAAAAAAQAAE6gAAAAAAAAAAAAAABECgDgAAAAAAAAAAAABAAATqAAAAAAAAAAAAAAAEQKAOAAAAAAAAAAAAAEAABOoAAAAAAAAAAAAAAARAoA4AAAAAAAAAAAAAQAAE6gAAAAAAAAAAAAAABECgDgAAAAAAAAAAAABAAATqAAAAAAAAAAAAAAAEQKAOAAAAAAAAAAAAAEAABOoAAAAAAAAAAAAAAARAoA4AAAAAAAAAAAAAQAAE6gAAAAAAAAAAAAAABECgDgAAAAAAAAAAAABAAATqAAAAAAAAAAAAAAAEQKAOAAAAAAAAAAAAAEAABOoAAAAAAAAAAAAAAARAoA4AAAAAAAAAAAAAQAAE6gAAAAAAAAAAAAAABECgDgAAAAAAAAAAAABAAATqAAAAAAAAAAAAAAAEQKAOAAAAAAAAAAAAAEAA/SJQf+KJJzR06FCFh4dr1qxZ2rZtW7CHBAAAAAAAAAAAAAAY4Pp8oP7CCy9o1apVevDBB7Vr1y5NmjRJS5cuVXV1dbCHBgAAAAAAAAAAAAAYwPp8oP7LX/5SX/nKV3TXXXdp7NixevrppxUZGak//OEPwR4aAAAAAAAAAAAAAGAACwn2AM7F4XBo586deuCBB/zLjEajFi9erM2bNwd8TXt7u9rb2/2Pm5qaJEn19fVyOp09O2AAQI9zOp1qaWlRXV2dzGZzsIcDAOhD2EYAAM6GbQQA4FzYTgDA4NPc3CxJ8nq95123TwfqtbW1crvdSktL67Y8LS1NhYWFAV/zyCOP6KGHHjpjeV5eXo+MEQAAAAAAAAAAAADQ/zQ3NysuLu6c6/TpQP3TeOCBB7Rq1Sr/Y4/Ho/r6eiUlJclgMARxZACAy8FqtSonJ0dlZWWKjY0N9nAAAH0I2wgAwNmwjQAAnAvbCQAYfLxer5qbm5WZmXnedft0oJ6cnCyTySSLxdJtucViUXp6esDXhIWFKSwsrNuy+Pj4nhoiACBIYmNjOcABAATENgIAcDZsIwAA58J2AgAGl/NVpncy9vA4LkloaKimTZumdevW+Zd5PB6tW7dOBQUFQRwZAAAAAAAAAAAAAGCg69MV6pK0atUq3XnnnZo+fbpmzpypxx57THa7XXfddVewhwYAAAAAAAAAAAAAGMD6fKD+2c9+VjU1NfrhD3+oqqoqTZ48We+8847S0tKCPTQAQBCEhYXpwQcfPGN6DwAA2EYAAM6GbQQA4FzYTgAAzsXg9Xq9wR4EAAAAAAAAAAAAAAB9TZ+eQx0AAAAAAAAAAAAAgGAhUAcAAAAAAAAAAAAAIAACdQAAAAAAAAAAAAAAAiBQBwAAAAAAAAAAAAAgAAJ1AEDQPfLII5oxY4ZiYmKUmpqqm266SUVFRd3WaWtr08qVK5WUlKTo6GitWLFCFoul2zqlpaW69tprFRkZqdTUVH3729+Wy+XqzY8CAOhhP/vZz2QwGHTffff5l7GNAIDB69SpU/rCF76gpKQkRUREaMKECdqxY4f/ea/Xqx/+8IfKyMhQRESEFi9erCNHjnR7j/r6et1xxx2KjY1VfHy87r77btlstt7+KACAy8ztdusHP/iB8vLyFBERoeHDh+snP/mJvF6vfx22EwCAC0GgDgAIug0bNmjlypXasmWL1q5dK6fTqSVLlshut/vXuf/++/X666/rxRdf1IYNG1RRUaFbbrnF/7zb7da1114rh8Ohjz/+WH/605/03HPP6Yc//GEwPhIAoAds375d//d//6eJEyd2W842AgAGp4aGBs2ZM0dms1lvv/22Dh06pP/93/9VQkKCf52f//znevzxx/X0009r69atioqK0tKlS9XW1uZf54477tDBgwe1du1avfHGG9q4caO++tWvBuMjAQAuo0cffVRPPfWUfvvb3+rw4cN69NFH9fOf/1y/+c1v/OuwnQAAXAiDt+vlWAAA9AE1NTVKTU3Vhg0bNG/ePDU1NSklJUWrV6/WrbfeKkkqLCxUfn6+Nm/erNmzZ+vtt9/Wddddp4qKCqWlpUmSnn76aX3nO99RTU2NQkNDg/mRAACXyGazaerUqXryySf105/+VJMnT9Zjjz3GNgIABrHvfve7+uijj/Thhx8GfN7r9SozM1Pf/OY39a1vfUuS1NTUpLS0ND333HO6/fbbdfjwYY0dO1bbt2/X9OnTJUnvvPOOli9frvLycmVmZvba5wEAXF7XXXed0tLS9Oyzz/qXrVixQhEREfrrX//KdgIAcMGoUAcA9DlNTU2SpMTEREnSzp075XQ6tXjxYv86Y8aMUW5urjZv3ixJ2rx5syZMmOAPSiRp6dKlslqtOnjwYC+OHgDQE1auXKlrr72227ZAYhsBAIPZa6+9punTp+u2225TamqqpkyZomeeecb//PHjx1VVVdVtGxEXF6dZs2Z120bEx8f7QxJJWrx4sYxGo7Zu3dp7HwYAcNldccUVWrdunYqLiyVJe/fu1aZNm7Rs2TJJbCcAABcuJNgDAACgK4/Ho/vuu09z5szR+PHjJUlVVVUKDQ1VfHx8t3XT0tJUVVXlX6drUNL5fOdzAID+6/nnn9euXbu0ffv2M55jGwEAg1dJSYmeeuoprVq1St/73ve0fft2ff3rX1doaKjuvPNO/3d8oG1A121Eampqt+dDQkKUmJjINgIA+rnvfve7slqtGjNmjEwmk9xutx5++GHdcccdksR2AgBwwQjUAQB9ysqVK3XgwAFt2rQp2EMBAPQBZWVl+sY3vqG1a9cqPDw82MMBAPQhHo9H06dP1//8z/9IkqZMmaIDBw7o6aef1p133hnk0QEAgu0f//iH/va3v2n16tUaN26c9uzZo/vuu0+ZmZlsJwAAF4WW7wCAPuOee+7RG2+8offff1/Z2dn+5enp6XI4HGpsbOy2vsViUXp6un8di8VyxvOdzwEA+qedO3equrpaU6dOVUhIiEJCQrRhwwY9/vjjCgkJUVpaGtsIABikMjIyNHbs2G7L8vPzVVpaKun0d3ygbUDXbUR1dXW3510ul+rr69lGAEA/9+1vf1vf/e53dfvtt2vChAn64he/qPvvv1+PPPKIJLYTAIALR6AOAAg6r9ere+65Ry+//LLWr1+vvLy8bs9PmzZNZrNZ69at8y8rKipSaWmpCgoKJEkFBQXav39/t4OctWvXKjY29oyTbACA/mPRokXav3+/9uzZ4/+ZPn267rjjDv99thEAMDjNmTNHRUVF3ZYVFxdryJAhkqS8vDylp6d320ZYrVZt3bq12zaisbFRO3fu9K+zfv16eTwezZo1qxc+BQCgp7S0tMho7B6BmEwmeTweSWwnAAAXjpbvAICgW7lypVavXq1XX31VMTEx/jmo4uLiFBERobi4ON19991atWqVEhMTFRsbq3vvvVcFBQWaPXu2JGnJkiUaO3asvvjFL+rnP/+5qqqq9P3vf18rV65UWFhYMD8eAOASxMTEaPz48d2WRUVFKSkpyb+cbQQADE7333+/rrjiCv3P//yPPvOZz2jbtm363e9+p9/97neSJIPBoPvuu08//elPNXLkSOXl5ekHP/iBMjMzddNNN0nyVbRfc801+spXvqKnn35aTqdT99xzj26//XZlZmYG8dMBAC7V9ddfr4cffli5ubkaN26cdu/+/9u705Covj+O459pyjCvYYvlSEaLbYRgC0XRMi2kJlFUFhbZpAlthLQgQRsh0V5EZEVlRSUVkaHQRjklPqgki4TAkooCK8KENkfR+T8IL01zS6sp4/97v+DC3HPPnPOd+2QYPnPOLdWuXbuUmpoqie8JAEDz2bxer7eliwAA/LfZbDbL9pycHLlcLklSTU2NVq5cqdzcXHk8HsXFxWn//v0+22s9f/5cixcvltvtVkhIiObPn68tW7aodWv+PwYA/0+cTqdiY2O1Z88eSXxHAMB/WUFBgdasWaPHjx+rZ8+eWrFihdLT083rXq9XGzZs0KFDh1RdXa1Ro0Zp//796tu3r9mnqqpKy5YtU35+vlq1aqUZM2Zo7969MgyjJT4SACBA3r9/r3Xr1unChQt68+aNIiMjlZycrPXr1ysoKEgS3xMAgOYhUAcAAAAAAAAAAAAAwALPUAcAAAAAAAAAAAAAwAKBOgAAAAAAAAAAAAAAFgjUAQAAAAAAAAAAAACwQKAOAAAAAAAAAAAAAIAFAnUAAAAAAAAAAAAAACwQqAMAAAAAAAAAAAAAYIFAHQAAAAAAAAAAAAAACwTqAAAAAAAAAAAAAABYIFAHAAAAAAAB43K5ZLPZZLPZlJeXF9Cx3W63Ofa0adMCOjYAAAAAAFYI1AEAAAAA+IGvA+KvjydPnrR0af+s+Ph4VVZWKiEhwWz7XsDucrmaHY6PHDlSlZWVmjVrVoAqBQAAAADgx1q3dAEAAAAAAPzr4uPjlZOT49MWHh7u16+2tlZBQUF/q6x/Vtu2bRURERHwcYOCghQREaHg4GB5PJ6Ajw8AAAAAwLdYoQ4AAAAAQBMaA+KvD7vdLqfTqWXLlikjI0OdO3dWXFycJKmsrEwJCQkyDENdu3bVvHnz9PbtW3O8jx8/KiUlRYZhyOFwaOfOnXI6ncrIyDD7WK3oDgsL07Fjx8zzFy9eaNasWQoLC1PHjh01depUPXv2zLzeuPp7x44dcjgc6tSpk5YuXaq6ujqzj8fjUWZmpqKiotS2bVtFR0fryJEj8nq9io6O1o4dO3xquH///h9bof/s2TPL3QCcTmfA5wIAAAAAoDkI1AEAAAAA+A3Hjx9XUFCQiouLdeDAAVVXV2v8+PEaNGiQSkpKdPnyZb1+/dpnm/LVq1fr5s2bunjxoq5evSq326179+791Lx1dXWKi4tTaGioioqKVFxcLMMwFB8fr9raWrNfYWGhKioqVFhYqOPHj+vYsWM+oXxKSopyc3O1d+9ePXr0SAcPHpRhGLLZbEpNTfVbmZ+Tk6MxY8YoOjr6127YD0RFRamystI8SktL1alTJ40ZMybgcwEAAAAA0Bxs+Q4AAAAAQBMKCgpkGIZ5npCQoHPnzkmS+vTpo23btpnXsrKyNGjQIG3evNlsO3r0qKKiolReXq7IyEgdOXJEJ0+e1IQJEyR9CeW7dev2UzWdOXNGDQ0NOnz4sGw2m6QvYXdYWJjcbrcmTZokSerQoYP27dsnu92u/v37KzExUdevX1d6errKy8t19uxZXbt2TRMnTpQk9erVy5zD5XJp/fr1unPnjoYNG6a6ujqdPn3ab9V6cyUnJ8tut/u0eTweJSYmSpLsdru5VXxNTY2mTZumESNGaOPGjb80HwAAAAAAv4tAHQAAAACAJowbN07Z2dnmeUhIiPl6yJAhPn0fPHigwsJCnwC+UUVFhT5//qza2loNHz7cbO/YsaP69ev3UzU9ePBAT548UWhoqE97TU2NKioqzPOBAwf6hNgOh0MPHz6U9GX7drvdrrFjx1rOERkZqcTERB09elTDhg1Tfn6+PB6PkpKSfqrWRrt37zaD+0aZmZmqr6/365uamqr379/r2rVratWKDfYAAAAAAC2DQB0AAAAAgCaEhIR8d4vzr8N1Sfrw4YOmTJmirVu3+vV1OBzNfva4zWaT1+v1afv62ecfPnzQkCFDdOrUKb/3hoeHm6/btGnjN25DQ4MkKTg4uMk6Fi5cqHnz5mn37t3KycnR7Nmz1a5du2Z9hm9FRET43cfQ0FBVV1f7tGVlZenKlSu6c+eO3x8GAAAAAAD4mwjUAQAAAAAIoMGDB+v8+fPq0aOHWrf2/9ndu3dvtWnTRrdv31b37t0lSe/evVN5ebnPSvHw8HBVVlaa548fP9anT5985jlz5oy6dOmi9u3b/1KtMTExamho0M2bN/1WjjeaPHmyQkJClJ2drcuXL+vWrVu/NFdznT9/Xps2bdKlS5fUu3fvPzoXAAAAAABNYc80AAAAAAACaOnSpaqqqlJycrLu3r2riooKXblyRQsWLFB9fb0Mw1BaWppWr16tGzduqKysTC6Xy29b8/Hjx2vfvn0qLS1VSUmJFi1a5LPafO7cuercubOmTp2qoqIiPX36VG63W8uXL9fLly+bVWuPHj00f/58paamKi8vzxzj7NmzZh+73S6Xy6U1a9aoT58+GjFiRGBulIWysjKlpKQoMzNTAwcO1KtXr/Tq1StVVVX9sTkBAAAAAPgRAnUAAAAAAAIoMjJSxcXFqq+v16RJkxQTE6OMjAyFhYWZofn27ds1evRoTZkyRRMnTtSoUaP8nsW+c+dORUVFafTo0ZozZ45WrVrls9V6u3btdOvWLXXv3l3Tp0/XgAEDlJaWppqamp9asZ6dna2ZM2dqyZIl6t+/v9LT0/Xx40efPmlpaaqtrdWCBQt+4840raSkRJ8+fVJWVpYcDod5TJ8+/Y/OCwAAAADA99i83z6QDQAAAAAA/HVOp1OxsbHas2dPS5fip6ioSBMmTNCLFy/UtWvXH/Z1uVyqrq5WXl7eH6vnb8wBAAAAAIDECnUAAAAAAPAdHo9HL1++1MaNG5WUlNRkmN6ooKBAhmGooKAgoPUUFRXJMAydOnUqoOMCAAAAAPA9rVu6AAAAAAAA8G/Kzc1VWlqaYmNjdeLEiWa9Z9u2bVq7dq0kyeFwBLSeoUOH6v79+5IkwzACOjYAAAAAAFbY8h0AAAAAAAAAAAAAAAts+Q4AAAAAAAAAAAAAgAUCdQAAAAAAAAAAAAAALBCoAwAAAAAAAAAAAABggUAdAAAAAAAAAAAAAAALBOoAAAAAAAAAAAAAAFggUAcAAAAAAAAAAAAAwAKBOgAAAAAAAAAAAAAAFgjUAQAAAAAAAAAAAACw8D83rQHH5VkF4AAAAABJRU5ErkJggg==", + "text/plain": [ + "
" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "if domain.comm.rank == 0:\n", + " import matplotlib.pyplot as plt\n", + "\n", + " fig = plt.figure(figsize=(25, 8))\n", + " plt.plot(freq, 20 * np.log10(np.abs(p_mic) / np.sqrt(2) / 2e-5), linewidth=2)\n", + " plt.grid(True)\n", + " plt.xlabel(\"Frequency [Hz]\")\n", + " plt.ylabel(\"SPL [dB]\")\n", + " plt.xlim([freq[0], freq[-1]])\n", + " plt.ylim([0, 90])\n", + " plt.legend()\n", + " plt.show()" + ] + } + ], + "metadata": { + "jupytext": { + "formats": "ipynb,py:light" + }, + "kernelspec": { + "display_name": "Python 3 (DOLFINx complex)", + "language": "python", + "name": "python3-complex" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.12.7" + } + }, + "nbformat": 4, + "nbformat_minor": 2 +} diff --git a/chapter2/helmholtz_code.py b/chapter2/helmholtz_code.py new file mode 100644 index 00000000..09694ebc --- /dev/null +++ b/chapter2/helmholtz_code.py @@ -0,0 +1,290 @@ +# --- +# jupyter: +# jupytext: +# formats: ipynb,py:light +# text_representation: +# extension: .py +# format_name: light +# format_version: '1.5' +# jupytext_version: 1.16.6 +# kernelspec: +# display_name: Python 3 (DOLFINx complex) +# language: python +# name: python3-complex +# --- + +# # Implementation +# Author: Antonio Baiano Svizzero and Jørgen S. Dokken +# +# In this tutorial, you will learn how to: +# - Define acoustic velocity and impedance boundary conditions +# - Compute acoustic sound pressure for multiple frequencies +# - Compute the Sound Pressure Level (SPL) at a given microphone position +# +# ## Test problem +# As an example, we will model a plane wave propagating in a tube. +# While it is a basic test case, the code can be adapted to way more complex problems where velocity and impedance boundary conditions are needed. +# We will apply a velocity boundary condition $v_n = 0.001$ to one end of the tube and an impedance $Z$ computed with the Delaney-Bazley model, +# supposing that a layer of thickness $d = 0.02$ and flow resistivity $\sigma = 1e4$ is placed at the second end of the tube. +# The choice of such impedance (the one of a plane wave propagating in free field) will give, as a result, a solution with no reflections. +# +# First, we create the mesh with gmsh, also setting the physical group for velocity and impedance boundary conditions and the respective tags. + +# + +import gmsh + +gmsh.initialize() + +# meshsize settings +meshsize = 0.02 +gmsh.option.setNumber("Mesh.MeshSizeMax", meshsize) +gmsh.option.setNumber("Mesh.MeshSizeMax", meshsize) + + +# create geometry +L = 1 +W = 0.1 + +gmsh.model.occ.addBox(0, 0, 0, L, W, W) +gmsh.model.occ.synchronize() + +# setup physical groups +v_bc_tag = 2 +Z_bc_tag = 3 +gmsh.model.addPhysicalGroup(3, [1], 1, "air_volume") +gmsh.model.addPhysicalGroup(2, [1], v_bc_tag, "velocity_BC") +gmsh.model.addPhysicalGroup(2, [2], Z_bc_tag, "impedance") + +# mesh generation +gmsh.model.mesh.generate(3) +# - + +# Then we import the gmsh mesh with the ```dolfinx.io.gmshio``` function. + +# + +from mpi4py import MPI +from dolfinx import fem, io, default_scalar_type, geometry +from dolfinx.fem.petsc import LinearProblem +import ufl +import numpy as np +import numpy.typing as npt + +mesh_data = io.gmshio.model_to_mesh(gmsh.model, MPI.COMM_WORLD, 0, gdim=3) +domain = mesh_data.mesh +assert mesh_data.facet_tags is not None +facet_tags = mesh_data.facet_tags +# - + +# We define the function space for our unknown $p$ and define the range of frequencies we want to solve the Helmholtz equation for. + +# + + +V = fem.functionspace(domain, ("Lagrange", 1)) + +# Discrete frequency range +freq = np.arange(10, 1000, 5) # Hz + +# Air parameters +rho0 = 1.225 # kg/m^3 +c = 340 # m/s + +# - + +# ## Boundary conditions +# +# The Delaney-Bazley model is used to compute the characteristic impedance and wavenumber of the porous layer, +# treated as an equivalent fluid with complex valued properties +# +# $$ +# \begin{align} +# Z_c(\omega) &= \rho_0 c_0 \left[1 + 0.0571 X^{-0.754} - j 0.087 X^{-0.732}\right],\\ +# k_c(\omega) &= \frac{\omega}{c_0} \left[1 + 0.0978 X^{-0.700} - j 0.189 X^{-0.595}\right],\\ +# \end{align} +# $$ +# +# where $X = \frac{\rho_0 f}{\sigma}$. +# +# With these, we can compute the surface impedance, that in the case of a rigid passive absorber placed on a rigid wall is given by the formula +# $$ +# \begin{align} +# Z_s = -j Z_c cot(k_c d). +# \end{align} +# $$ +# +# Let's create a function to compute it. +# +# + + +# + +# Impedance calculation +def delany_bazley_layer(f, rho0, c, sigma): + X = rho0 * f / sigma + Zc = rho0 * c * (1 + 0.0571 * X**-0.754 - 1j * 0.087 * X**-0.732) + kc = 2 * np.pi * f / c * (1 + 0.0978 * (X**-0.700) - 1j * 0.189 * (X**-0.595)) + Z_s = -1j * Zc * (1 / np.tan(kc * d)) + return Z_s + + +sigma = 1.5e4 +d = 0.01 +Z_s = delany_bazley_layer(freq, rho0, c, sigma) +# - + +# Since we are going to compute a sound pressure spectrum, all the variables that depend on frequency (that are $\omega$, $k$ and $Z$) need to be updated in the frequency loop. +# To make this possible, we will initialize them as dolfinx constants. +# Then, we define the value for the normal velocity on the first end of the tube + +omega = fem.Constant(domain, default_scalar_type(0)) +k = fem.Constant(domain, default_scalar_type(0)) +Z = fem.Constant(domain, default_scalar_type(0)) +v_n = 1e-5 + +# We also need to specify the integration measure $ds$, by using ```ufl```, and its built in integration measures + +ds = ufl.Measure("ds", domain=domain, subdomain_data=facet_tags) + +# ## Variational Formulation +# We can now write the variational formulation. + +# + +p = ufl.TrialFunction(V) +v = ufl.TestFunction(V) + +a = ( + ufl.inner(ufl.grad(p), ufl.grad(v)) * ufl.dx + + 1j * rho0 * omega / Z * ufl.inner(p, v) * ds(Z_bc_tag) + - k**2 * ufl.inner(p, v) * ufl.dx +) +L = -1j * omega * rho0 * ufl.inner(v_n, v) * ds(v_bc_tag) +# - + +# The class ```LinearProblem``` is used to setup the PETSc backend and assemble the system vector and matrices. +# The solution will be stored in a `dolfinx.fem.Function`, ```p_a```. + +# + +p_a = fem.Function(V) +p_a.name = "pressure" + +problem = LinearProblem( + a, + L, + u=p_a, + petsc_options={ + "ksp_type": "preonly", + "pc_type": "lu", + "pc_factor_mat_solver_type": "mumps", + }, +) + + +# - + +# ## Computing the pressure at a given location +# Before starting our frequency loop, we can build a function that, given a microphone position, +# computes the sound pressure at its location. +# We will use the a similar method as in [Deflection of a membrane](../chapter1/membrane_code). +# However, as the domain doesn't deform in time, we cache the collision detection + + +class MicrophonePressure: + def __init__(self, domain, microphone_position): + """Initialize microphone(s). + + Args: + domain: The domain to insert microphones on + microphone_position: Position of the microphone(s). + Assumed to be ordered as ``(mic0_x, mic1_x, ..., mic0_y, mic1_y, ..., mic0_z, mic1_z, ...)`` + + """ + self._domain = domain + self._position = np.asarray( + microphone_position, dtype=domain.geometry.x.dtype + ).reshape(3, -1) + self._local_cells, self._local_position = self.compute_local_microphones() + + def compute_local_microphones( + self, + ) -> tuple[npt.NDArray[np.int32], npt.NDArray[np.floating]]: + """ + Compute the local microphone positions for a distributed mesh + + Returns: + Two lists (local_cells, local_points) containing the local cell indices and the local points + """ + points = self._position.T + bb_tree = geometry.bb_tree(self._domain, self._domain.topology.dim) + + cells = [] + points_on_proc = [] + + cell_candidates = geometry.compute_collisions_points(bb_tree, points) + colliding_cells = geometry.compute_colliding_cells( + domain, cell_candidates, points + ) + + for i, point in enumerate(points): + if len(colliding_cells.links(i)) > 0: + points_on_proc.append(point) + cells.append(colliding_cells.links(i)[0]) + + return np.asarray(cells, dtype=np.int32), np.asarray( + points_on_proc, dtype=domain.geometry.x.dtype + ) + + def listen( + self, recompute_collisions: bool = False + ) -> npt.NDArray[np.complexfloating]: + if recompute_collisions: + self._local_cells, self._local_position = self.compute_local_microphones() + if len(self._local_cells) > 0: + return p_a.eval(self._local_position, self._local_cells) + else: + return np.zeros(0, dtype=default_scalar_type) + + +# The pressure spectrum is initialized as a numpy array and the microphone location is assigned + +# + +p_mic = np.zeros((len(freq), 1), dtype=complex) + +mic = np.array([0.5, 0.05, 0.05]) +microphone = MicrophonePressure(domain, mic) +# - + +# ## Frequency loop +# +# Finally, we can write the frequency loop, where we update the values of the frequency-dependent variables and solve the system for each frequency + +for nf in range(0, len(freq)): + k.value = 2 * np.pi * freq[nf] / c + omega.value = 2 * np.pi * freq[nf] + Z.value = Z_s[nf] + + problem.solve() + p_a.x.scatter_forward() + + p_f = microphone.listen() + p_f = domain.comm.gather(p_f, root=0) + + if domain.comm.rank == 0: + assert p_f is not None + p_mic[nf] = np.hstack(p_f) + +# ## SPL spectrum +# After the computation, the pressure spectrum at the prescribed location is available. +# Such a spectrum is usually shown using the decibel (dB) scale to obtain the SPL, with the RMS pressure as input, +# defined as $p_{rms} = \frac{p}{\sqrt{2}}$. + +if domain.comm.rank == 0: + import matplotlib.pyplot as plt + + fig = plt.figure(figsize=(25, 8)) + plt.plot(freq, 20 * np.log10(np.abs(p_mic) / np.sqrt(2) / 2e-5), linewidth=2) + plt.grid(True) + plt.xlabel("Frequency [Hz]") + plt.ylabel("SPL [dB]") + plt.xlim([freq[0], freq[-1]]) + plt.ylim([0, 90]) + plt.legend() + plt.show() From 3876b4c44ded9bdebaf8c30c82e17a99d6be2144 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?J=C3=B8rgen=20Schartum=20Dokken?= Date: Mon, 3 Feb 2025 22:45:54 +0100 Subject: [PATCH 4/5] Make helmholtz backwards compatible (#248) --- chapter2/helmholtz_code.ipynb | 18 ++++++++++++++---- chapter2/helmholtz_code.py | 19 +++++++++++++++---- 2 files changed, 29 insertions(+), 8 deletions(-) diff --git a/chapter2/helmholtz_code.ipynb b/chapter2/helmholtz_code.ipynb index f85949bb..c33bcace 100644 --- a/chapter2/helmholtz_code.ipynb +++ b/chapter2/helmholtz_code.ipynb @@ -143,16 +143,26 @@ "outputs": [], "source": [ "from mpi4py import MPI\n", - "from dolfinx import fem, io, default_scalar_type, geometry\n", + "from dolfinx import (\n", + " fem,\n", + " io,\n", + " default_scalar_type,\n", + " geometry,\n", + " __version__ as dolfinx_version,\n", + ")\n", "from dolfinx.fem.petsc import LinearProblem\n", "import ufl\n", "import numpy as np\n", "import numpy.typing as npt\n", + "from packaging.version import Version\n", "\n", "mesh_data = io.gmshio.model_to_mesh(gmsh.model, MPI.COMM_WORLD, 0, gdim=3)\n", - "domain = mesh_data.mesh\n", - "assert mesh_data.facet_tags is not None\n", - "facet_tags = mesh_data.facet_tags" + "if Version(dolfinx_version) > Version(\"0.9.0\"):\n", + " domain = mesh_data.mesh\n", + " assert mesh_data.facet_tags is not None\n", + " facet_tags = mesh_data.facet_tags\n", + "else:\n", + " domain, _, facet_tags = mesh_data\n" ] }, { diff --git a/chapter2/helmholtz_code.py b/chapter2/helmholtz_code.py index 09694ebc..2eee4715 100644 --- a/chapter2/helmholtz_code.py +++ b/chapter2/helmholtz_code.py @@ -63,16 +63,27 @@ # + from mpi4py import MPI -from dolfinx import fem, io, default_scalar_type, geometry +from dolfinx import ( + fem, + io, + default_scalar_type, + geometry, + __version__ as dolfinx_version, +) from dolfinx.fem.petsc import LinearProblem import ufl import numpy as np import numpy.typing as npt +from packaging.version import Version mesh_data = io.gmshio.model_to_mesh(gmsh.model, MPI.COMM_WORLD, 0, gdim=3) -domain = mesh_data.mesh -assert mesh_data.facet_tags is not None -facet_tags = mesh_data.facet_tags +if Version(dolfinx_version) > Version("0.9.0"): + domain = mesh_data.mesh + assert mesh_data.facet_tags is not None + facet_tags = mesh_data.facet_tags +else: + domain, _, facet_tags = mesh_data + # - # We define the function space for our unknown $p$ and define the range of frequencies we want to solve the Helmholtz equation for. From c87c45793e67f7093e661ce3bd90c2812bfae13b Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?J=C3=B8rgen=20Schartum=20Dokken?= Date: Sun, 19 Jan 2025 18:47:57 +0100 Subject: [PATCH 5/5] Native arm build (#240) * Try new publish workflow * Fix platform choice * Update workflow to merge artifacts to one place * Set retention days to minimum * Fix name * change name to pattern * Minor improvements * Try new strategy * Remove duplicate tag * Revert platform tags * Use correct dockerfile * remove first login * Add merge skip if not pushing on release or tagging * Set push to false as we don't need the images before merging. * Remove push * Remove conditional to check that simpler code runs * Add back pushtrue * Revert comment --- .github/workflows/publish_docker.yml | 96 +++++++++++++++++++++------- 1 file changed, 73 insertions(+), 23 deletions(-) diff --git a/.github/workflows/publish_docker.yml b/.github/workflows/publish_docker.yml index 753bf033..bd09b935 100644 --- a/.github/workflows/publish_docker.yml +++ b/.github/workflows/publish_docker.yml @@ -1,4 +1,5 @@ -name: Publish tutorial docker image +# Recipe based on: https://docs.docker.com/build/ci/github-actions/multi-platform/#distribute-build-across-multiple-runners +name: Build and publish platform dependent docker image on: push: branches: @@ -15,8 +16,12 @@ env: IMAGE_NAME: ${{ github.repository }} jobs: - build-and-push-image: - runs-on: ubuntu-latest + build: + strategy: + matrix: + os: ["ubuntu-24.04", "ubuntu-24.04-arm"] + runs-on: ${{ matrix.os }} + permissions: contents: read packages: write @@ -25,12 +30,6 @@ jobs: - name: Checkout repository uses: actions/checkout@v4 - - name: Set up QEMU - uses: docker/setup-qemu-action@v3 - - - name: Set up Docker Buildx - uses: docker/setup-buildx-action@v3 - - name: Log in to the Container registry uses: docker/login-action@v3 with: @@ -38,30 +37,81 @@ jobs: username: ${{ github.actor }} password: ${{ secrets.GITHUB_TOKEN }} + - name: Set up Docker Buildx + uses: docker/setup-buildx-action@v3 + - name: Extract metadata (tags, labels) for Docker id: meta uses: docker/metadata-action@v5 with: images: ${{ env.REGISTRY }}/${{ env.IMAGE_NAME }} - - name: Build Docker image + - name: Set architecture tag (amd64) + if: ${{ matrix.os == 'ubuntu-24.04' }} + run: echo "ARCH_TAG=amd64" >> $GITHUB_ENV + + - name: Set architecture tag (arm) + if: ${{ contains(matrix.os, 'arm') }} + run: echo "ARCH_TAG=arm64" >> $GITHUB_ENV + + - name: Build and push by digest + id: build uses: docker/build-push-action@v6 with: - context: . - load: true - push: false file: docker/Dockerfile - platforms: linux/amd64 - tags: ${{ steps.meta.outputs.tags }} + platforms: ${{ env.ARCH_TAG }} labels: ${{ steps.meta.outputs.labels }} + outputs: type=image,"name=${{ env.REGISTRY }}/${{ env.IMAGE_NAME}}",push-by-digest=true,name-canonical=true,push=true - - name: Build (arm) and push (amd/arm) Docker image - uses: docker/build-push-action@v6 + - name: Export digest + run: | + mkdir -p ${{ runner.temp }}/digests + digest="${{ steps.build.outputs.digest }}" + touch "${{ runner.temp }}/digests/${digest#sha256:}" + + - name: Upload digest if: github.event_name == 'push' + uses: actions/upload-artifact@v4 with: - context: . - push: true - file: docker/Dockerfile - platforms: linux/amd64,linux/arm64 - tags: ${{ steps.meta.outputs.tags }} - labels: ${{ steps.meta.outputs.labels }} + name: digests-${{ env.ARCH_TAG }} + path: ${{ runner.temp }}/digests/* + if-no-files-found: error + retention-days: 1 + + merge-and-publish: + if: github.event_name == 'push' + runs-on: ubuntu-latest + needs: + - build + steps: + - name: Download digests + uses: actions/download-artifact@v4 + with: + path: ${{ runner.temp }}/digests + pattern: digests-* + merge-multiple: true + + - name: Log in to the Container registry + uses: docker/login-action@v3 + with: + registry: ${{ env.REGISTRY }} + username: ${{ github.actor }} + password: ${{ secrets.GITHUB_TOKEN }} + + - name: Extract metadata (tags, labels) for Docker + id: meta + uses: docker/metadata-action@v5 + with: + images: ${{ env.REGISTRY }}/${{ env.IMAGE_NAME }} + + + - name: Create manifest list and push + working-directory: ${{ runner.temp }}/digests + run: | + docker buildx imagetools create $(jq -cr '.tags | map("-t " + .) | join(" ")' <<< "$DOCKER_METADATA_OUTPUT_JSON") \ + $(printf '${{ env.REGISTRY }}/${{ env.IMAGE_NAME}}@sha256:%s ' *) + + - name: Inspect image + run: | + docker buildx imagetools inspect ${{ env.REGISTRY }}/${{ env.IMAGE_NAME}}:${{ steps.meta.outputs.version }} +