{ "cells": [ { "cell_type": "markdown", "id": "3757f472", "metadata": {}, "source": [ ".. index:: Equations; Wave Equation\n", "\n", "# Wave Equation: double slit domain\n", "\n", "In the following we will consider the wave equation\n", "\\begin{align*}\n", "\\partial_{tt}\\psi - \\Delta\\psi &= 0, && \\text{in } \\Omega\\times(0,T), \\\\\n", "\\psi &= g, && \\text{on } \\Gamma_D\\times(0,T), \\\\\n", "\\nabla\\psi\\cdot n &= 0 && \\text{on } \\partial\\Omega\\setminus\\Gamma_D,\n", "\\end{align*}\n", "where $n$ is the outward unit normal and $g$ is given Dirichlet data\n", "given on part of the domain boundary. We repeat a double slit experiment,\n", "i.e., we have a wave channel which ends in two slits on the right and\n", "forcing the wave on the left end by defining $g=\\frac{1}{10\\pi}\\cos(10\\pi t)$.\n", "\n", "We discretize the wave equation by an explicit symplectic time stepping\n", "scheme for two unknowns $\\psi^n$ and $p^n$ where $p=-\\partial_t\\psi$ based\n", "on\n", "\\begin{align*}\n", "\\partial_t \\psi = -p, \\\\\n", "\\partial_t p = -\\triangle\\psi\n", "\\end{align*}\n", "The discretization first updates $\\psi$ by half a time step then updates\n", "$p$ and concludes by another update step for $\\psi$. So given\n", "$(\\psi^n,p^n)$ we compute\n", "\\begin{align*}\n", "\\psi^{n+\\frac{1}{2}} &= \\psi^n - \\frac{\\Delta t}{2}p^n, \\\\\n", "p^{n+1} &=\n", "\\begin{cases}\n", "p^n - \\triangle\\psi^{n+\\frac{1}{2}} & \\text{in the interior}, \\\\\n", "g(t^{n+1}) & \\text{on }\\Gamma_D\n", "\\end{cases} \\\\\n", "\\psi^{n+1} &= \\psi^{n+\\frac{1}{2}} - \\frac{\\Delta t}{2}p^{n+1}.\n", "\\end{align*}\n", "Note that the update for $\\psi$ does not involve any spatial derivatives\n", "so can be performed directly on the level of the degree of freedom\n", "vectors. We use *numpy* for these steps. For the update of $p$ we\n", "use a variational formulation\n", "\\begin{align*}\n", "\\vec{p} &= \\vec{p} + M^{-1}S\\vec{\\psi}\n", "\\end{align*}\n", "where $S$ is the stiffness matrix\n", "$S=(\\int_\\Omega\\nabla\\varphi_j\\nabla\\varphi_i)_{ij}$ and $M$ is the mass\n", "matrix, i.e., $M=\\int_\\Omega\\varphi_j\\varphi_i)_{ij}$. To simplify the\n", "computation we use a lumped version of this matrix, i.e,,\n", "$M\\approx{\\rm diag}(\\int_\\Omega\\varphi_i)$. In this form the algorithm\n", "only involves a matrix vector multiplication and no implicit solving step\n", "is required.\n", "\n", "The following test case is taken from\n", "https://firedrakeproject.org/demos/linear_wave_equation.py.html" ] }, { "cell_type": "code", "execution_count": 1, "id": "99984a12", "metadata": { "execution": { "iopub.execute_input": "2024-02-29T12:01:25.993610Z", "iopub.status.busy": "2024-02-29T12:01:25.993316Z", "iopub.status.idle": "2024-02-29T12:01:27.164778Z", "shell.execute_reply": "2024-02-29T12:01:27.163400Z" }, "lines_to_next_cell": 2 }, "outputs": [], "source": [ "from matplotlib import pyplot as plt\n", "import numpy\n", "import dune.fem as fem\n", "from dune.grid import reader\n", "from dune.alugrid import aluConformGrid as leafGridView\n", "from dune.fem.space import lagrange as solutionSpace\n", "from ufl import TrialFunction, TestFunction, grad, dot, dx\n", "from dune.ufl import Constant, BoxDirichletBC\n", "T = 3\n", "dt = 0.005\n", "t = 0" ] }, { "cell_type": "markdown", "id": "0fd305a5", "metadata": {}, "source": [ ".. index:: Grid construction; gmsh (from file)\n", "\n", "We use *gmsh* to define the domain and then set up a first order\n", "scalar Lagrange space which we will use both for $\\psi$ and $p$. We can\n", "construct general grids by either defining the grids directly in Python\n", "(as demonstrated in the following example) or by reading the grid from\n", "files using readers provided by Dune, e.g., Dune Grid Format (dgf) files\n", "or Gmsh files." ] }, { "cell_type": "code", "execution_count": 2, "id": "237cfb6e", "metadata": { "execution": { "iopub.execute_input": "2024-02-29T12:01:27.168879Z", "iopub.status.busy": "2024-02-29T12:01:27.168364Z", "iopub.status.idle": "2024-02-29T12:07:46.271837Z", "shell.execute_reply": "2024-02-29T12:07:46.270851Z" }, "lines_to_next_cell": 2 }, "outputs": [], "source": [ "domain = (reader.gmsh, \"wave_tank.msh\")\n", "gridView = leafGridView( domain, dimgrid=2 )\n", "gridView.hierarchicalGrid.loadBalance()\n", "V = solutionSpace(gridView, order=1, storage=\"numpy\")\n", "\n", "p = V.interpolate(0,name=\"p\")\n", "phi = V.interpolate(0,name=\"phi\")\n", "pVec = p.as_numpy\n", "phiVec = phi.as_numpy" ] }, { "cell_type": "markdown", "id": "76a664a1", "metadata": {}, "source": [ "Next we define an operator for the stiffness matrix including the\n", "boundary condition which are time dependent so we use a `Constant` for\n", "this. We use the `BoxDirichletBC` class which is derived from the more\n", "general `DirichletBC` class which takes a function space, the boundary\n", "function $g$ as a ufl expression, and finally a description of the part\n", "$\\Gamma_D$ of the boundary where this boundary condition is defined.\n", "This can be a ufl expression which evaluates to $0$ for\n", "$x\\not\\in\\Gamma_D$, e.g., an ufl `conditional`.\n", "\n", "Note that the stiffness matrix does not depend on time so we can\n", "assemble the matrix once and extract the corresponding `scipy`\n", "sparse matrix." ] }, { "cell_type": "code", "execution_count": 3, "id": "0ba5bb69", "metadata": { "execution": { "iopub.execute_input": "2024-02-29T12:07:46.275393Z", "iopub.status.busy": "2024-02-29T12:07:46.275086Z", "iopub.status.idle": "2024-02-29T12:09:37.355102Z", "shell.execute_reply": "2024-02-29T12:09:37.353977Z" } }, "outputs": [], "source": [ "u = TrialFunction(V)\n", "v = TestFunction(V)\n", "p_in = Constant(0.0, name=\"g\")\n", "# the following is equivalent to\n", "# x = SpatialCoordinate(V)\n", "# bc = DirichletBC(V, p_in, conditional(x[0]<1e-10,1,0))\n", "bc = BoxDirichletBC(V, p_in, [None,None],[0,None], eps=1e-10)\n", "\n", "from dune.fem.operator import galerkin\n", "op = galerkin([dot(grad(u),grad(v))*dx,bc])\n", "S = op.linear().as_numpy\n", "lapPhi = V.interpolate(1,name=\"e\")\n", "lapPhiVec = lapPhi.as_numpy" ] }, { "cell_type": "markdown", "id": "2b284047", "metadata": {}, "source": [ ".. index::\n", " pair: Quadratures; Mass lumping\n", "\n", "Next we construct a lumped (diagonal) approximation $M$ to the mass matrix.\n", "We then use this to form the matrix $\\Delta t M^{-1}S$ where $M^{-1}$ is\n", "easy to compute since $M$ is diagonal.\n", "To compute $M$ we can either use the *row sum* approach\n", "which is a diagonal matrix with entries $m_{ii} = \\sum_j\\int_\\Omega\\varphi_i\\varphi_j$.\n", "We could either setup a mass matrix using" ] }, { "cell_type": "code", "execution_count": 4, "id": "2bceeb67", "metadata": { "execution": { "iopub.execute_input": "2024-02-29T12:09:37.361929Z", "iopub.status.busy": "2024-02-29T12:09:37.361568Z", "iopub.status.idle": "2024-02-29T12:10:28.051457Z", "shell.execute_reply": "2024-02-29T12:10:28.050623Z" } }, "outputs": [], "source": [ "M = galerkin(u*v*dx).linear().as_numpy\n", "diag = M.sum(axis=1)" ] }, { "cell_type": "markdown", "id": "407bad85", "metadata": {}, "source": [ "Another approach is to construct a mass operator $ = \\int_\\Omega u\\varphi_i$\n", "and applying that to $u\\equiv 1$:" ] }, { "cell_type": "code", "execution_count": 5, "id": "92ac8945", "metadata": { "execution": { "iopub.execute_input": "2024-02-29T12:10:28.058188Z", "iopub.status.busy": "2024-02-29T12:10:28.057616Z", "iopub.status.idle": "2024-02-29T12:10:28.271634Z", "shell.execute_reply": "2024-02-29T12:10:28.270682Z" } }, "outputs": [], "source": [ "from scipy.sparse import dia_matrix\n", "lumping = galerkin(u*v*dx)\n", "lumped = lapPhi.copy()\n", "lumping(lapPhi,lumped) # note that lapPhi=1\n", "N = len(lumped.as_numpy)\n", "M = dia_matrix(([dt/lumped.as_numpy],[0]),shape=(N,N) )\n", "S1 = M*S" ] }, { "cell_type": "markdown", "id": "eff6669c", "metadata": {}, "source": [ "An alternative approach is to use an interpolation quadrature, where the\n", "integral for the mass matrix \\int_\\Omega \\varphi_i\\varphi_j$ is\n", "approximated by a quadrature based on the Lagrange points.\n", "We can change the quadrature by using the metadata argument when defining\n", "the ufl measure:\n", "\n", ".. todo:: currently mass lumping only for schemes\n" ] }, { "cell_type": "code", "execution_count": 6, "id": "90345318", "metadata": { "execution": { "iopub.execute_input": "2024-02-29T12:10:28.275953Z", "iopub.status.busy": "2024-02-29T12:10:28.275709Z", "iopub.status.idle": "2024-02-29T12:11:51.317298Z", "shell.execute_reply": "2024-02-29T12:11:51.315531Z" } }, "outputs": [], "source": [ "from dune.fem.scheme import galerkin\n", "lumping = galerkin(u*v*dx(metadata={\"quadrature_rule\":\"lumped\"})==0)\n", "lumped = lumping.linear().as_numpy.diagonal()\n", "M = dia_matrix(([dt/lumped],[0]),shape=(N,N) )\n", "S2 = M*S" ] }, { "cell_type": "markdown", "id": "795f6127", "metadata": {}, "source": [ "We can now set up the time loop and perform our explicit time stepping\n", "algorithm. Note that the computation is carried out completely using\n", "*numpy* and *scipy* algorithms we the exception of setting the\n", "boundary conditions. This is done using the `setConstraints` method on\n", "the stiffness operator which we constructed passing in the boundary\n", "conditions. At the time of writing it is not yet possible to extract a\n", "sparse matrix and vector encoding the boundary constraint." ] }, { "cell_type": "code", "execution_count": 7, "id": "8624a139", "metadata": { "execution": { "iopub.execute_input": "2024-02-29T12:11:51.323086Z", "iopub.status.busy": "2024-02-29T12:11:51.322830Z", "iopub.status.idle": "2024-02-29T12:12:08.034437Z", "shell.execute_reply": "2024-02-29T12:12:08.033280Z" } }, "outputs": [ { "data": { "image/jpeg": "/9j/4AAQSkZJRgABAQEAMgAyAAD/2wBDAAgGBgcGBQgHBwcJCQgKDBQNDAsLDBkSEw8UHRofHh0aHBwgJC4nICIsIxwcKDcpLDAxNDQ0Hyc5PTgyPC4zNDL/2wBDAQkJCQwLDBgNDRgyIRwhMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjL/wAARCAHuAeIDASIAAhEBAxEB/8QAHwAAAQUBAQEBAQEAAAAAAAAAAAECAwQFBgcICQoL/8QAtRAAAgEDAwIEAwUFBAQAAAF9AQIDAAQRBRIhMUEGE1FhByJxFDKBkaEII0KxwRVS0fAkM2JyggkKFhcYGRolJicoKSo0NTY3ODk6Q0RFRkdISUpTVFVWV1hZWmNkZWZnaGlqc3R1dnd4eXqDhIWGh4iJipKTlJWWl5iZmqKjpKWmp6ipqrKztLW2t7i5usLDxMXGx8jJytLT1NXW19jZ2uHi4+Tl5ufo6erx8vP09fb3+Pn6/8QAHwEAAwEBAQEBAQEBAQAAAAAAAAECAwQFBgcICQoL/8QAtREAAgECBAQDBAcFBAQAAQJ3AAECAxEEBSExBhJBUQdhcRMiMoEIFEKRobHBCSMzUvAVYnLRChYkNOEl8RcYGRomJygpKjU2Nzg5OkNERUZHSElKU1RVVldYWVpjZGVmZ2hpanN0dXZ3eHl6goOEhYaHiImKkpOUlZaXmJmaoqOkpaanqKmqsrO0tba3uLm6wsPExcbHyMnK0tPU1dbX2Nna4uPk5ebn6Onq8vP09fb3+Pn6/9oADAMBAAIRAxEAPwD3+qGoatFp1zZwyxSN9qlESMrJwxIA4LBj1z8oOACTgCr9VbrTra9lhkuBIxiZXVRK4TcpDAlQcNggEZB6UAMudY0uyn8i61Kzgmxu8uWdVbHrgnNRf8JDon/QY0//AMCU/wAaJ/8AkZbH/rzuP/Q4a0qAM3/hIdE/6DGn/wDgSn+NH/CQ6J/0GNP/APAlP8a0qKAM3/hIdE/6DGn/APgSn+NH/CQ6J/0GNP8A/AlP8a0qKAM3/hIdE/6DGn/+BKf40f8ACQ6J/wBBjT//AAJT/GtKigDN/wCEh0T/AKDGn/8AgSn+NH/CQ6J/0GNP/wDAlP8AGtKigDN/4SHRP+gxp/8A4Ep/jR/wkOif9BjT/wDwJT/GtKigDLbxJoSnDa1pwPvdJ/jSf8JNoH/Qc03/AMC4/wDGqWvf8fyf9ch/M1l151XHOnNx5djyq2YunUcOXbzOh/4SbQP+g5pv/gXH/jR/wk2gf9BzTf8AwLj/AMa56is/7Sf8v4mf9qv+X8Tof+Em0D/oOab/AOBcf+NH/CTaB/0HNN/8C4/8a56ij+0n/L+If2q/5fxOh/4SbQP+g5pv/gXH/jR/wk2gf9BzTf8AwLj/AMa56ij+0n/L+If2q/5fxOh/4SbQP+g5pv8A4Fx/40f8JNoH/Qc03/wLj/xrnqKP7Sf8v4h/ar/l/E6H/hJtA/6Dmm/+Bcf+NH/CTaB/0HNN/wDAuP8AxrnqKP7Sf8v4h/ar/l/E6H/hJtA/6Dmm/wDgXH/jR/wk2gf9BzTf/AuP/Gueoo/tJ/y/iH9qv+X8Tof+Em0D/oOab/4Fx/40f8JNoH/Qc03/AMC4/wDGueoo/tJ/y/iH9qv+X8Tof+Em0D/oOab/AOBcf+NH/CTaB/0HNN/8C4/8a56ij+0n/L+If2q/5fxOh/4SbQP+g5pv/gXH/jR/wk2gf9BzTf8AwLj/AMa56ij+0n/L+If2q/5fxOh/4SbQP+g5pv8A4Fx/40f8JNoH/Qc03/wLj/xrnqKP7Sf8v4h/ar/l/E6H/hJtA/6Dmm/+Bcf+NH/CTaB/0HNN/wDAuP8AxrnqKP7Sf8v4h/ar/l/E6H/hJtA/6Dmm/wDgXH/jR/wk2gf9BzTf/AuP/Gueoo/tJ/y/iH9qv+X8Tof+Em0D/oOab/4Fx/40f8JNoH/Qc03/AMC4/wDGueoo/tJ/y/iH9qv+X8Tof+Em0D/oOab/AOBcf+NH/CTaB/0HNN/8C4/8a56ij+0n/L+If2q/5fxOh/4SbQP+g5pv/gXH/jR/wk2gf9BzTf8AwLj/AMa56ij+0n/L+If2q/5fxOh/4SbQP+g5pv8A4Fx/40f8JNoH/Qc03/wLj/xrnqKP7Sf8v4h/ar/l/E6H/hJtA/6Dmm/+Bcf+NPi8RaJNKkUWs6fJI7BURLpCWJ4AAzya5uo5f9baf9flv/6OSrp5g5zUeXcunmbnNR5d33O6ooor0j1gooooAKKKKACiiigAooooAKKKKACiiigDNn/5GWx/687j/wBDhrSrNn/5GWx/687j/wBDhrSoAKKKKACiiigAooooAKKKKACiiigDnde/4/k/65D+ZrLrU17/AI/k/wCuQ/may6+fxX8aR8xjP48vUKKKK5zmCiiigAooooAKKKKACiiigAooooAKKKKACiiigAooooAKKKKACiiigAooooAKKKKACiiigAooooAKKKKACo5f9baf9flv/wCjkqSo5f8AW2n/AF+W/wD6OStaH8WPqjbD/wAaHqvzO6ooor6M+qCiiue0Xxno+veIdX0Sykc3elsEl8wBRIclW2DOSFYbScAZIxkEGi41FtNrodDRRRQIKKKKACiiigAooooAKKKKAM2f/kZbH/rzuP8A0OGtKs2f/kZbH/rzuP8A0OGtKgAooooAKKKKACiiigAooooAKKKKAOc184v0/wCuQ/mays1o+I226hGP+mQ/max/Mrza2FlKbkj4zH4iMcVNeZPmjNQ+ZS76weFkjlWJiybNGaiD07dWUqDRoqyZJmimbqXNZOm0aqpcdRSZpc1m4tFqSCiiikMKKKKACiiigAooooAKKKKACiiigAooooAKKKKACiiigAooooAKKKKACiiigAqOX/W2n/X5b/8Ao5KkqOX/AFtp/wBflv8A+jkrWh/Fj6o2w/8AGh6r8zuqKKK+jPqjjfiX4yXwd4XeWFv+JleZgs1BXKsQcyEHqqdeh5Kg4zmvm/wr4juvCHiC11ezTzDDlZYC5UTRkYZCR+BGcgMFODitj4k+MF8ZeK3ubWR20y1XyLMMCu5erSbSeCx74B2qmRkVyFcVWo3LTofS4HBxjQamtZbn2lZXlvqNhbX1pJ5ltcxLNE+CNyMAQcHkcEdanrxn4GeLonsZfCl5MizQs01gGIBdDlnQcclWy3JJIc9lr1+a8tbeaGGa5hjlmJESO4DSEddoPXqOldcZcyueBXoujUcH0J6Kpy6tpsLTrLqFpGYMCYPMo8vPTdzxn3pf7U08OU+3Wu5YvPK+cuRH139fu+/SqMS3RVRNU0+R4EjvrZ2uATCqzKTIB12884welW6ACiiigAooooAzZ/8AkZbH/rzuP/Q4a0qzZ/8AkZbH/rzuP/Q4a0qACiiigAooooAKKKKACiiigAooooA5TxPn+0o8f88R/M1i4NbviPH9ox5/55D+ZrIwKPb8ulj88zSjzYyo79SHmjcam2ik2ChYiD3R57oTWzIw5pQ9KY6QpVc1GQrVYjxJTg9QFSKMkVnLCQn8JccTOO5aDU4GqoepFeuKrgmtjsp4tMnzS5qIPTwa8+dFx3O2FVPYfRSZpc1g4tGylcKKKKkoKKKKACiiigAooooAKKKKACiiigAooooAKKKKACiiigAooooAKjl/1tp/1+W//o5KkqOX/W2n/X5b/wDo5K1ofxY+qNsP/Gh6r8zuqhuziynI/wCebfyqaobz/jyuP+ubfyr6M+qPixE/dr9BTtlSxr+6T6Cn7K8ly1P0CNHRHZfCG0juPiLYiUN8scrKVcqQdh6EEH1H0JFfQmraJNqF2lxBdpAwEe4SQmQExyCRMYYY+Yc+o9OteC/B5cfEaz/65S/+gGvpOu/Du8D5POI8uKa8kc7L4dvVed7PVEt3bzhC5tizRiaUSSZO8bjkYU8bc96ik8Ib4zCt3HFBgOkccBGyUQrCOrH5Nqj5Ov8AtV09FbnlnOweGJIrpZmvlYPIstwogxuZZpJl2Hd8o3SEHO7IHbk10VFFABRRRQAUUUUAZs//ACMtj/153H/ocNaVZs//ACMtj/153H/ocNaVABRRRQAUUUUAFFFFABRRRQAUUUUAcZ4tuTDq0Shc5gB/8easIX3qprovFKqdTiJAz5I/9CasMoh7CvTpOjyJSjdn5fnDksdVs+pGL9Se9Srdof4qjNvE38NNNqmOOKqVLCy6WPPVWotmWxcKeMjNPEgNZb2j5yjUwfaIfXFQ8tozXuSNI4uotzZ4NIUFZsV8QcOKtx3cb8BhXFVwFek7o6I4mnPSSJClN5FSBwaUqDWSrSjpNFOjGWsGRh6lV6jKU3kU5UqdVaExqTpvUtBqeDVRXqZXrza+EcT0aOKUibNLTA1OBrzZ02jvjO4tFFFZGgUUUUAFFFFABRRRQAUUUUAFFFFABRRRQAUUUUAFFFFABUcv+ttP+vy3/wDRyVJUcv8ArbT/AK/Lf/0cla0P4sfVG2H/AI0PVfmd1UN5/wAeVx/1zb+VTVDef8eVx/1zb+VfRn1R8dxJ+6T/AHRTtlSRJ+5T/dFP2V4TlqfqsKPuo7L4RLj4iWf/AFyl/wDQDX0fXzr8Jlx8QrL/AK5y/wDoBr6Kr1MI70z4XiCPLjWvJBRRRXSeIFFFFABRRRQAUUUUAZs//Iy2P/Xncf8AocNaVZs//Iy2P/Xncf8AocNaVABRRRQAUUUUAFFFFABRRRQAUUUUAeZ/ETUpLPX4I0zg2qt/489cl/bs3ckV3vjXTY7zWYZHxkW6r/483+Nc2dBgIr6nB4jDRoRU1rY/Ns2nRWNqKUdbmbDrz9yatx6+ueSDQ+gp0GKqTeH2x8tdX+x1Dzv9nl5GvFrMMmM4q4l1DIMhhXHyaZcwnjdSCW6tsZyaieXUZ6wY3h4y+CR2LwxyjIx+FVJbV0OUPHrWFBrMiMN2a2bbVkkHzMKxeHr0dtUc86M4bofHdTQnDcj3q/Bfxy8ZwfQ1XZ4Zhxj61VltmQ71/SuapQo19JqzFCo4vRm8HDUFc1hQ3zwsA/IrVt7pJQMGvIxOW1aHvQ1R2wxCl7tREhUigNipOCKYy1zQqqXuzHOk4+9AkV6lDVUBwalR65sThOqOnD4noywDTqiVqeDXjVKbTPVhUTHUUUVgbBRRRQAUUUUAFFFFABRRRQAUUUUAFFFFABRRRQAVHL/rbT/r8t//AEclSVHL/rbT/r8t/wD0cla0P4sfVG2H/jQ9V+Z3VQ3n/Hlcf9c2/lU1Q3n/AB5XH/XNv5V9GfVHyTDGfIj/AN0fyp/lmrkEP+jxcfwD+VSeT7V8y6mp+ywiuVHSfCpCvxAsif8AnnL/AOgGvoWvAPh0ssfjixaCNHfEgw7lRjY2eQD/ACr2DXL3UbTUbY24ufs42MywWxlEg8xRIGwpIxHkjGCT69K9rAu9G5+dcTK2Pfojforj9QvtdS31honvlkjguGiEVmH2OrgQCP5D5m9clvvY/wBmp7651V9Un+yaheR2s1qxQf2Wx+yNtBWTJXMhz/B156cGuw+fOporj7G/8RPqWnLP9oNttCvm02m4O+VS7Hb+6wqxOFOCd+OSK7CgAooooAKxtU1drLWdOsYrq0Et0/NvKMOyA/MwbcMYzwNpJPHqRs0UAZMsjHxZaRmF1VbKciQldrZeHgc549wK1qzZ/wDkZbH/AK87j/0OGtKgAooooAKKKKACiiigAooooAKKKKAPMviJrZ03xBbw5I3Wqv8A+PuP6Vyy+KuOTj610vxI0Y6j4it5h/DaKv8A4+5/rXHt4YcDgV9jgYYR4aHPvY/Pc0WFeMqc+9zQj8TqxGWFX4dfgfGSK5h/D0qdAagfS7mPO3NdUsHhJ7HnvD4aXwyO5W+tZscjPpTZYLeYZGPyrgi97bn+KrltrM8RAfNZPLXDWnIiWAktYO50c+jqRuUcVmy2UsB+TNWbfXUdcFquJcRTilGVanpPUw5q1PSSMuG7lhIDZrTh1AOoUtTJreNxkVQeBo2yvSqcadXdDtCr5M2GRZV7VVbzbZsqTUdtclflarpKypWLTg7PYxacHZ7Fqx1NZcK5w1aqsGFclNC8Tb0q/YatyElNeXjsrjUXtKO50U6jhqtUbrLTOVNOjkWRQQQQaVlzXiQm4PkqGsoKS54Co9Tq1VM4NSo9c2Kwya5kb4bENOzLINLUamng14lSFmexCd0LRRRWJqFFFFABRRRQAUUUUAFFFFABRRRQAUUUUAFRy/620/6/Lf8A9HJUlRy/620/6/Lf/wBHJWtD+LH1Rth/40PVfmd1UN5/x5XH/XNv5VNUN5/x5XH/AFzb+VfRn1R8yW8f+jRf7g/lUvl1NbR/6LDx/Av8ql8r2r4yU9WfrsKnuo3vh0m3xxYH2l/9FtXudeJ+AU2+NdPP/XT/ANFtXtlfR5Y74derPgOI3fHN+SCiiivQPBCiiigAooooAKKKKAM2f/kZbH/rzuP/AEOGtKs2f/kZbH/rzuP/AEOGtKgAooooAKKKKACiiigAooooAKKKKAOA8c6tFY63DE5AY26t/wCPN/hXOLr9sR1H50z4tWNzdeKrWSEnaLJBx673rgTp19H0LV9ngMvw9TDQlJ6tHwOZ4KjUxdSTnZtnoY1S2k7j86eJLeUdVrzoLfRdd1SLqN5DjIaup5Yl8EjzZZX/ACSO/ksLeZcgLWZc6JGwO0CsO18RyIQHJrattejmXDEVm6GJovR3RhLD4mjsZU2lyRN8oNEUk1uec10azQTL2qpcWqNnArSOIctJouOKcvdqIqw6gTw1XBKrisx7co3FPiYr3pyhF6oJ04vWJbfAbIqeCXBGTVEljTBOUfmpcLqxDp8ysdBgSrjrWXd2rRnenWp7O6DAZNaDRiZDjmuTmdGWuxypypSKel6mUYRyHj3rpI5A6gg5rjb20e3bzFq7pOr4YRSHiuTMMvjiI+1p7nTF8vvw26nSstNBwacjh1BByDSMtfOwbT9nMqav78SVGqYGqitip0avNxmHs7o78LXurMmFLTQadXjyVmerF3CiiipKCiiigAooooAKKKKACiiigAooooAKjl/1tp/1+W//AKOSpKjl/wBbaf8AX5b/APo5K1ofxY+qNsP/ABoeq/M7qobz/jyuP+ubfyqaobz/AI8rj/rm38q+jPqj57tY/wDRIf8Armv8ql8upbSP/Q4P+ua/yqby6+AlPVn6fCp7qNXwOu3xlpx95P8A0W9ey15D4NTb4v04/wC1J/6KevXq+ryd3wy9WfE5874xvyQUUUV6h4wUUUUAFFFFABRRRQBmz/8AIy2P/Xncf+hw1pVmz/8AIy2P/Xncf+hw1pUAFFFFABRRRQAUUUUAFFFFABRRRQBwnjR4F1mESY3fZ16/7zVzJ+yN2WqvxfvLq38W2qwZ2mxQnHr5klcENWvkPO6vtcvy6VTDQmpbo+BzPLJ1MXUmpbs9Ea0tpOwqpNo0EvRRXIQ+IrhPvBq0rbxPz8xrqeDxFP4WeY8BiaesWWbrw8vJUVkzabPbNlc4rpINbhnxkjmrJMM69qccRWp6TQRxVelpURy0F5PAQGzWxbagJFAai60+NgSorMaFoX4rR8lVeZs/Z1le1mbjsrrxiqjfI1QwTEgA1K4LCslHldjBQ5XYtwhZBilnscrkCs6O5MD4NbVrdpMoGazqKUHdGVVTpvmWxhmSS1l5zit/TL9JAATUN/p6zRblHNc2ZptOuOc7c1XJDEwstzRRjiYabnfz2yXEWBiuV1C1ks5t6jgVuaLqa3UYUkZq/fWSXUJBHNeVTrywlX2dTY5ablRk0/mZuiasJFEbt9K6IEMK88mhl027yMhc11+k34ubcAt8wFY5tgE17ekdEkqbvH4WaB4NSI1NIyKapwa8OUVVhZ7kRl7Ody4pp4qBGqYGvnsRS5We7QqcyHUUUVxHWFFFFABRRRQAUUUUAFFFFABRRRQAVHL/AK20/wCvy3/9HJUlRy/620/6/Lf/ANHJWtD+LH1Rth/40PVfmd1UN5/x5XH/AFzb+VTVDef8eVx/1zb+VfRn1R4hZx/6Fb/9c1/lU/l1JZR/6Bb8f8sl/lU/l+1fmsp+8z9BjU0RY8NSQ2vijTZJ5UiTzHXc7BRkxuAMmvSr3WYLG9htXimdpCgZ0A2x73CIWyQeWOOAehziuA8NJt8U6Yf+mkn/AKKevQr7R7LUpEkuY3LoMBkmeM9QRnaRnBAIz0PIwa+yyN3wi9WfJ5w74pvyRRufFFtbRXsn2S7kS1jmlBQJ+9WJgsuzLD7pPO7Ge2avvqQTUZLIWs7SrAZ0K7MSAHG1ct1zjrgc9agn8PaZcpdpJBJtu0ZJQs8i/Kxy4XDfJuPJ24z3zUv9jWf2pbr/AEnzlg+zhvtUv3P++uv+11zzmvYPKIrbXIboaWyW1wqalD50TsFwny79rYbOcemR71qVmQeH9Otl09YkuAunqVtgbqUhBjGDlvm44+bPFadABRRRQAUUUUAZs/8AyMtj/wBedx/6HDWlWbP/AMjLY/8AXncf+hw1pUAFFFFABRRRQAUUUUAFFFFABRRRQB554806C71uGSUAsLZV5/3mrkn0G1YcKKufFnxA2leK7W3BOGskf/x9x/SuKi8ZE8E19rl+FxTwsJwejR8BmeDxcsXUnB6Nm7N4chYHAFZdz4aIBKCrEHiuJyAWrVg1m3nHLCuzmxVLc8/nxlH4jj5NPubVvl3cVZtdRmhYK+a6qQQTg4waybvTUOSoFbRxKqaVEbRxcaq5aqHw6gJByaWXa4zWWIWhf2q9F8y1MoRjrEmVOMXeI1fker0JEgxWfcIyjIqK3vTFIAxxQ4OSuglTc43Ro3dkWUlRzWdDdPaThWJAzXS2ciXMfY5FUNX0nchdBzWNOsr+zmYUa65vZ1Da0u6iuowCwyah1nRUuIC6KMgVx9hqMlhdCNzjBr0PTrxL22ByCcciuHGU6mEmq0NjGtQnh6ilHY89trqXSrwI2QM16BpeoJfW4II3YrA8T6KJFM0S4PtWL4f1J7K6EUhI5xzXRiKNPH4f2kPiOqpFYml7WHxLc7HW9PE8DSKPmFc5pl49leeW3Az3rt4pEuYAwIIIrkNd05re4MqDjORXFltfnTw1XdHNRaa9nLZnYwSrNEGU5BFObrXO6BqO9BGxrojyK8jFYZ4es10MpXV4S3Q+NqsKaqKcGrCGvDx9Hqd+Cq9CcUtNFOrwZqzPai7oKKKKgsKKKKACiiigAooooAKKKKACo5f9baf9flv/AOjkqSo5f9baf9flv/6OStaH8WPqjbD/AMaHqvzO6qG8/wCPK4/65t/KpqhvP+PK4/65t/Kvoz6o8msY/wDiX23/AFyX+Qqx5VPsI/8AiXW3/XJP5CrHl1+VTn7zPs41NESaCm3xNph/6av/AOipK9Irz7R12+I9LP8A02f/ANEyV6DX3XD7vgl6s+czN3xDfkgooor2zzwooooAKKKKACiiigDNn/5GWx/687j/ANDhrSrNn/5GWx/687j/ANDhrSoAKKKKACiiigAooooAKKKKACiiigDyL4p+Hxqvie2nI+7Zqn/j7n+tcFJ4OIHAr0r4j63DpviK3hkIBNor/wDj7j+lcrH4ltJerCvuMurYqOFhyrSx8HmVfGwxdTkXu3OOuPDU8JyuaqiG9tG/iIFeirf2twOCtQ3NnbzKcAc16MMwntURywzSp8NWJydnrMyEK5P41tx6gJU5NUbzSVUkqKpxI0LYzVyjSqLmiazhRrLmias7gjNMgnUMATSRqZEqncxvCdwzxURin7pnGCfunQpEsy1lanpzoC6DGKTS9WXeEc11qQRXtvxg5Fc1SpLDSvLY5Jznhamuxxelas9tOI5CRzXeWjRX1t7kVwXiHSZLSUyxg8c1a8L68UkWKRuemDVYzDqvS9tR3NcXh416Xt6XzLPiPRDGxlQdO4qt4c1l7eYRSNgg4wa7y4hjvrToCCK821mwfT70yoMYNZYHERxVJ0Ku6JwlVYim6FTfoenAR3lt2IYVwPiHTGs7kyxjHPatzwxq4nhEbtWrrNkt3ak4ya8/DyngcU6cvhZyU5ywta0vmZHhfWfMjEMjfTNb+pWwurYjHIFeeQ7tN1HHQZr0TT7lbq0VsgnHNGZ0FRqrE0y8XT5J80NmcdbFrK+KngZrtLScTQhs9q53XrPypfMUe9WtEusoFJ+tXjILE0FVW5nWfPBVV8zf71NGarE8ZqWNulfN4mlzUzPD1OWZcU08VEhqUV8pXjZn0tGV0LRRRXMdAUUUUAFFFFABRRRQAUUUUAFRy/620/6/Lf8A9HJUlRy/620/6/Lf/wBHJWtD+LH1Rth/40PVfmd1UN5/x5XH/XNv5VNUN5/x5XH/AFzb+VfRn1R55p8f/EtteP8Alin8hVjy/an6dH/xLLT/AK4p/wCgirPlV+Pzn7zPpVU0ItNTb4g0o/8ATZ//AETJXdVxlom3XdKP/Td//RMldnX6Dw274BerPGxzvWuFFFFe+cYUUUUAFFFFABRRRQBmz/8AIy2P/Xncf+hw1pVmz/8AIy2P/Xncf+hw1pUAFFFFABRRRQAUUUUAFFFFABRRRQB4P8bNKur7xjZyw7to09FOPXzJP8a8xbStRgOQWr6B8ezWya7AkxUMbZSM+m5q5UxWc44CGvvsrzKdLCU4OOiR8Vj82qUcVOm43SZ5ZBfX9ocOG4rf0/xDIQA5rob7RbaVSVUVzF1pHkSFkr1FXoYhaqzMlicPilrGzN4XyzpmqFyQDmqtmSpCk1oyWhliyOax5Y05HPyRpSt0HWFwjMFJFa8tgtzDwM8VxFxJLYT7jkDNdd4f1iO5VULc1niaUoL2kNjLGUJwj7ansctqlnNp9x5ig4Bro/DHiBZNsTt7Vuaxo6XtoxVQTjNeZXEcujaju5C5qqM6ePouD3RtRlTzCi6cviR65qNgl9anABBHFeYalZTaTqHmKCADzXoHhnWEvLVUZsnFQeK9KWeBpFXmvOwNeeFrvD1Njz8FWlhazpVNnoSeGNYW5t1jZucVa8Qaat1AzY5rzzRbyTT9Q8tjgBq9Rtrhb2yHc4qMdQeFxCrU9mTjqLw1bnhs9TzmxlfTb/YSQueK9Hsblbq0GTniuI16x8q43qOhzWr4fvD5YQn2rfH0liKSqLc0xsVWpKtHch8QWWJCyjmrHhrUCuInPtV/VYxNGa5+1BtrrI4GamFq2H5JmdOXtcPyPdHW6tEs9uSMGsLTyYJ8ds1sJMJbfmst1CXOR61zYaLjTdNnLSb5ZQZ0KyboxzU8D5H0rNhf92KsW0mJMetefWoe60c6fKzYjPFTCq0RqwtfE42HLNn0mDnzRH0UCivLZ6KCiiigAooooAKKKKACiiigAqOX/W2n/X5b/wDo5KkqOX/W2n/X5b/+jkrWh/Fj6o2w/wDGh6r8zuqhvP8AjyuP+ubfyqaobz/jyuP+ubfyr6M+qOS02P8A4ldp/wBcU/8AQRVny6ZpoH9lWf8A1wT/ANBFWsCvxad+ZnpquUjCz6tpKpNJCftLHcgUn/UyccgjmtHXNGu77Uba7tVtpGh2FPPcoYysiudpCn74GwnjA9elVwANY0j/AK+X/wDRMtdPX6Rwv/yL16s4q8uadzkrjwveS3VxKFs3VpWkYOx/0sGZZAsvynARVKL97g9hxVT/AIQzUgjhbuISvZ/ZxP5rZRfJ2eUo28Lvw+/O7qMc13FFfQmJzNl4XETaV59tZqtjNLMuz52TJby0VioJUbyxJwcqvXk101FFABRRRQAVk6/c3Fpa20ts11v+1wqyW9uZt0ZkUOGAViBs3HIx061rUUAc7qt3PZ+KLB/MhMclpOqKUIIO+HOTu5/IfjXj/iT45eKdF8R3+mxWGkNHbylEMkMpYjtnEgru/ipqDaZJpNwjbTtmH1+aKvCfijGB42muABi6ghn46coB/SqyX97mVajW1jypryta/wB/MvuLxNCdOlSrp+7LmXzVmvwuei6L8avEup+H9avnsdJE9gIWRUik2lXYqxI8zPHHp1ok+NXiVPCM2rix0ozpepbhfKk27SjMT/rM5yB3rzr4cbLnVdS0uViEvtPljAHdlww/9BNaMlg8fwt1kMOYb6ByPrlf61vioxo5vCi/hlbTpreP5q/zN6UYyy2tVt70ZRV/JuP+bOk0b49eKdR1zT7Gaw0dYrm5jhcpDKCAzAHGZOvNaSfGnxI3iZdMNjpXlG8FvnypN23ft678Z/CvFdCl8jxDpk3/ADzu4m/Jwa73+zCPi+LYA7Bqfm/8B37qjiNLCuDpaXjJ/db/ADNcmpwrTqqqr2UX/wClX/Q2dX+Pfimw1q/s4rDRzHb3EkSF4ZckKxAz+868VpXfxs8R2+h6Pdiy0r7ReRSSSqYpMACRlXA8zPIXvmvF9ekE3iLU5RyHu5WH4ua6XxKvlwaBEBhU0iAj/gW5j+rGu7O6VPD0IOnGzcrfLll+tjkyWPt60IVdfdv+S/U9D1f41+JLDw5o2ow2WlNNemcSh4pCo2MANv7z0POc0zSvjd4lv9A1e/ksdJEtm0AjCxSbSHLA5/eewrznxICvgjwzn+KS7YfTco/oaZ4cwPBfiXPeS0A/76esvZQ/sdV7e/3/AO4lvyHyr677Nbc9vlc9O0b40+I9RmuxNZaUqQWc1xlIpOqISP8Alp0ziodI+OHiPUNG1e6ks9JW4slikjVYpAGVn2Nn952LL+defeHsix8REfeGkS4/77jB/Qmqnggb7nWY2AMbaTcl89toDL/48q1zYOCqZdXrS3i3Z+kVL8djozBQpY+FOKtH3Lr1k0/wL3izx3rHinV4r+7WCCSOAQhbZWVSAzHJ3EnPzHv6VQtfE17ARuYsK6fwf4ftdW0eaeZQWWcoM+m1T/Wreo+CbcIWjUfhX3mSY3DvA0qdTex8jm+OwMMfVw847OxQsPFfngK7c1oyXS3Eeeua4y60iTT7jgnANdJoymdApNenXoUor2lPY8rE4ehBe1pbFaS5EE+T610+jXEd0oGQa5/XtJkWAyIDwKwdC1uSwvVSQ8ZxSlQWIo80HqiZYZYrDuVN6o9B13QBcWrOi84zXn1rdzaNqexyQua9o0q5h1LTx0ORXnvjjw/5cjTRr0OciuHLMbepLC1jjyrF+88NX2Z3Hh/U01CzUEgnFYHjLRA8TSItc14M1d7acQyEjBxzXpV60d7YnPORXJWpTwOM5o/Czjr05YDF3jseY+GdRksb3ymJGDivUDKt7ZY65FeX6jZG01HzF45rsdEv82wVj2rrzKjGry1obnTmlJVUq8DnNZsTBemRRjmum8PXx8gKT2qtrEazAmqWmMYXxSm/bULSJqP2+GSlujb1iMTIT7VkaeTBL7ZrUllEkeDVJIvnyKypaU+VnPRdqbgzaMoki5NZksQL5HrT/NZFwafB+8as4rk1Mox9ndomtpTGuCaGw0gNEsRRc1BFJmUDNSle7RKV7yRsQj92KVH2y0+FMxCqsh2yiuRWk2jjSu2dFA25QfWra1m2T7olrRTpXxObU+WbPZy6d42JBS0gpa+fe57S2CiiikMKKKKACiiigAooooAKjl/1tp/1+W//AKOSpKjl/wBbaf8AX5b/APo5K1ofxY+qNsP/ABoeq/M7qobz/jyuP+ubfyqaobz/AI8rj/rm38q+jPqjl9NP/ErtP+uCf+girWap6af+JXaf9cU/9BFWc1+STpe8yfaiKf8Aic6T/wBfL/8AomWuorlUOda0n/r5f/0TLXVV99w7HlwKXmxp31CiiivdGFFFFABRRRQAUUUUAeT/ABu/49dI7f67n8Y68f8Ailj+3dLI6nSbfI9Dg17d8WLE6gdJhC5IWY4xnvEP614d8U5VPjM2oOfsdrDbn6hc/wBaeR+9nE7dIO/zcDux1SLy7D009eeT+Vmv1RD8L/8Ako2k+haTP08tq7PW0SH4c+JcDG65hH4765P4WxqviyS+kICWNpLPk+uNo/8AQq3PEuoif4XPIn/LzqKxt+Clq6M1Tnn+FjHok/uc5foYYdSWW4qfRumvmpXf4NHm+mDOq2Y/6bp/6EK9ojuIP+FvuhxvwVH+9sxXi2nNt1O0b0mQ/wDjwr0e7uHi+MgkXJK6ki49RkDFTxdG8aXpU/8AbDqyNXqV4/zQ5fvueb3ysuoXKt94SsD9cmu01yNrvwn4Y1ML1tGtG+sbkD9DXL+JUWPxVq6KMKt7MAP+Bmup1Kcp8OvC1tnqbiX/AMiYr0+JPewsJf31/wCkyOPIpXxVOXdO/pa/5pEHxAiNpY+GbH+GLThJ+LsWNU9FBXwPrLjo15aqfwWU1o+PG+1eGvCd+TmSS1kib/gDAVW0NA3w51k91vrc/ow/rXNKX/GPprtFfNTSf4pl0HbM2pf8/J/nIf4YyW1lf4W0m5Df98Z/mBVPwwgi8N+JLwEiTyYLZSD0DyAt+iY/GreiMbbQPE14Oq2AgHt5kiKf0zVTwmWn0jxJYqAd1klwCexjlU/yZq5sFGX9lYhr+b8EoN/gdOYNPNKb7ci+fM3+qF0zxVc6FGbaIZRm8w89yAP6V1mm+NE1BdknDehrjU8LajqWjzaxbxhraFzE5zyCAD/7MKzNKDrfxqDg7sV9lkUIVsHTi4pq266Psz5jO8BhMRiq0kveT1O316RZIzIKzfD2sJFdrG5wCe9dTNoElxpW/aTlc15hdxyWN8yHIKNXt0alKVGUX0/A8bARo4mlKinqj3lLKLUtPyMHIryTxVosmmX7OqkKTXceA/EH2i1EUrfMvFWfGVnFeWzNgZxXm4WpVwmLdKfws8fBVamAxjpT2ZgeBfEDx4gkbpxzXaa4kd9ZHjORXk2lQy2d/kZxmvQ4dQD2gDHtW2YYWMcQq0NzTNMKo4hVaXU4c2rWWp7l4Ga73TdQ32gVj2rnrm3W4n3Ad6txK1untWmJarRV9y8W1iIR5tx2qQrMxPemWJaAYponEsm01oLa/usgVk3yxUZGMpckFCRDLceZxUkFucZArMnk8q42n1rodNAkhB61FX3IXRnW/dwuinJIYzgmrVmBJVPVx5WTjFTaE/mEDNRJfuudGc1+550TXqGNc1Hps26XFampwAWxPtXPabKBe7c96zpNVKTZnS/eUW+x1NzHugJ9q55JNl5tJ711wh8yz6dq4q/zBqA+tYYKaqc0exlgvevE7eyAe296zL4bJfxrQ0d/Mswap6wu181w0ZWxMoM5Yq0y/pcm6PFbMZ4rntHfPHtXQR9K+ez6FptnoYB2m0TClpBS18jLc9+OwUUUVJQUUUUAFFFFABRRRQAVHL/rbT/r8t//AEclSVHL/rbT/r8t/wD0cla0P4sfVG2H/jQ9V+Z3VQ3n/Hlcf9c2/lU1Q3n/AB5XH/XNv5V9GfVHI6af+JXaf9cU/wDQRVrNUtOP/ErtP+uKf+girWa/OJUfeZ5jqixH/idaT/18P/6JkrrK40StHq+kskLzH7Sw2oVB/wBTJzyQOK09c16TStRtoi9tFA2wuZ8gyAyKrbTkAbFJc5B49OtfYZJHlwiXmztoS5oXN+iuSuPFF5FdXEQazRVlaNi6H/RAJljDS/MMh1Yuv3eB3HNNh8WO7L5uo6bAF08XVwHjObclQQT8+WB5O0AEDGTyN3rmx19FctYeJJrq/wBJtjeac8l2jySQIMP5YLhZAd5AyVUbMMch+TtJHU0AFFFFABRRRQBzXiFVbXtN3AHFrckZ9d0NfJ3i64e78YavM5yTdSD8AcD9BX1lr7Kuv6buON1rcKPqXgr5R8a232Txrq8OMf6Szf8AfXzf1ro4cdP6/i19q1P7rS/4ByVuf20L7Wlb743NjwQCPDfi6ROJFsowD/smQbv0ov8Ad/wqqPP3Tq4x/wB+mqz8KEjvNb1HSpv9Ve2LIR7ggitb4nadF4d8KaPosbZLXMk7fUDH8iKuumuIqdPuub5KMl+a/E7oYmkstq0L++6kdP8AwB3/APJWea6VEZ9XsogMl541A+rAV6hHZG4+M8rbcpDdtO/oAgz/AErgvBNv9p8baNF/09o35HP9K9DsJWbxR4m1A9Y7a6bPpwQKjij95WoUF1Ul/wCBOCOvKeaFPFVlvCF16+9/keT6hdG91K6uz1nmeQ/8CJP9a67VsjwZ4VB6+ROce3mmuIruPFhNvLpelkBTp+nQRMB/fZd7fq/6V6PEk17CEF1lf5KMl+qMskp/7Uktoxf6L9RvjmTyNE8LaaRhobFp2/7aOT/SotEfZ8OdYz/y0v7dB9drn+lM+IEnnaho8n8LaPa4/wC+SD+oNLZHZ8NnHeXVwf8AvmH/AOzrnrJRyGCXVQfzcoyf43Iwz58z9ak//bn+hPoCfatD8TWZ/i07zx9Y5Fb/ABqn4ZItPDfiC+5EkiRWUZ/32LN/47GR+NWvDkggtfEEzHAXSJh+LMiD9WFUvDZF1oGu6fgmRY4r2MD/AKZttb/x2Qn/AIDXNg3J5TiEv5v0hf8AC52Y7lWawvt7n38zt+hc03xbFoPhzUdNitC93fMN0xb5QgGAMeuS351zOnTlNRic/wB8VoweHNR1hfPs4S8anYT79f6itK38D6lARLPHtA5r63h6lGjhqbjaKer11b279tDwM3xOGhjK3PNcz038tEew6RcQS6KgLL93kGvIvGdih1NpIhnJ5xVm48Sz6VF9mBOQMdaz9Ov31fUlWbkk8V7GCwLw9SdRu6lrY+Uy/A1cLUliPsk3hrz7F/MIIBNdHf62JYtjNzWy3h9ItP8AMC4+WvMdXneLUWQEgKa2pOljKjkuhVF08xrOa3R1un2H2k7wKk1DdZJg8CtXwVGLq0VvUVB43g8iAkDHFc/tr4r2LONVnLGexkZ+jzi6mAznmuivbPbalgO1cP4OmMl8FJ716pf24/swn/ZrLMJewrxj3MszvQxKijzWK5K6hsPrXeWcHm2YPtXmVxKY9b25/ir1nQQJNMB77aM1fsqcZoebR5IQkupwOvE2916c10vhh/PgArm/GY8u5P1re8DHfGpp4r/clMeKXNgIzJ/EtvsiY1m+GZcz7c966DxYoFoTjtXJeGJD9vx/tVjhX7TBtmGHXPgpeR3uqR/6CSPSuDtJNmq4969Fvk3aefYV5q3ya0MdzWGVS56UkZZcrxnHyPT7P57NPcVwniP91fZ967nTTmxQ+1cT4vGLrOO9cmVv/a5xM8CrV0vU6rw4++wBB9KNbXCg1B4UO6w6+lW9bH7kGuVvlzFoxnG0/mVdFb58V00XSuU0U/va6qLpXl8Qx1OnC6V2WBS0gpa+JnufQQ2CiiioKCiiigAooooAKKKKACo5f9baf9flv/6OSpKjl/1tp/1+W/8A6OStaH8WPqjbD/xoeq/M7qobz/jyuP8Arm38qmqG8/48rj/rm38q+jPqji9OP/EstP8Arin/AKCKtZqlpx/4llp/1xT/ANBFWc18ZKjqz5t1dSS3P/E80n/r4f8A9EyV2Fcban/ie6V/18P/AOiZK7Kvo8sjy4dLzZ7GBlelcKKKK9A7AooooAKKKKACiiigDlvE6GTXdKUf88Lg59Pmhr5g+IsqzePtWZeglC/iFAr6j8QEjXdOI/59Ln/0KGvkjxLI8vijVXfO43cmc/7xro4cotZhiq1+kF+b/Q5a9XmrQp/yqX4uP+R1fwehMnjgOOkVs7MfQZA/rW/8cF8x9HuY2DRMJFBHr8v/ANeuc+Hlw2naZ4n1JDiSCzRFP+++2rHia/bV/hlp9zMczQak0RPsUJ/pW1avNcR0mvhUHB+rUpf5FrL1LDTxt9YzSt5WS/ORg/D5gnj7Rienn4/Q12+lyBrzxfEPvGyuAB7815z4XuPsvivSJs4C3kWT7bhmvRdMVIfihq+mOcJc/abce5KtiseI0oYyhVfTX/wGUX+p6mWO+GxsO8F+UzySu48bHzPFNxcKcx3EUM8Z9VaNSP8AD8K4ggqxBGCDg122vMLjQPDN5j5nsPIY+vluyj9MV3cSU2qdOp2lb703/wC2iyWVsXZ9Yv8ANfpcpeNxtfQFPUaPAfzZz/Wn2wz8OoCOi6rKD+MSf4U7x8u6TQJ1HySaPAoPuuQaNM/efDe8H/PLVYm/BonB/kKzxH/Iig/KH5xX6mODk1mab/nqf+3jNMjaTRvESqMkacG/AXEJNQeDcRXWqXjkiGDTZw+O5kXylH/fTj8q0/CyhoPEO77g0afd+a4/XFZ+iAr4M1t4/vNdWySf7mJD/wChBfyFc+CquOUYlf3mv/Aowj+p14+PPmkIf4X+L/yOz8A+JdN0fRJra8mVJWuWcA+hVR/Q10194p0u5tGWKVWJFeEXgbzQRnG3rS2Usv2lFVjyQOtfR5NgcNOjSlJNSa+R8jnmRU62Oq1uZ3buafiJWm1CSRR8p6Gl8LIV1eJj0Br0fTvCKXmmLLIgJK55rk9UtYtCvcpgYNfTUqtCtVag/eSscFHMadam8LDdKx6pLdRf2OVyM7a8O175tWlI6ZrrrXxA93b+WrE8YrFvtLeW480g9ajL8KsNKSn1OTKqDwdSXtOp3Pw+IjsV3HtTPiC4a2bHpUHhu4FpAFzjiq3iiY3cZA5rljSvmHtTgjBvMfadLnO+Ch/xNAe2a9i1Bh/ZTf7teVeFbX7PdBiO9ehajd/8S8rntUZvD2uKg10FnL9pi04nlV6C2vj/AHq9f8OcaYv0FeU+SZdY3e9eo6TJ5OngZ7VedLnoRgi85lelTj2OH8cEG4I962/AfEK/QVgeK/313j3rpfBsflQqaeK0y9RFiHy5dGJo+LiPsR+lcf4VXOo/8CrqPFku63I9q57wrHi7B96wwS5MC0zHDO2Cm+56LeYFi/8Au15jLzrX/Aq9J1J9tg3uK82UF9Y/GubJY2pzbM8t+2/I9N0v/jwj+lcZ4w/4+D9a7PThtsYx7VxXipt11+Nc2WL/AG2bMcF/HXzOg8I5+w1e1v8A49xVTwsu2x/KrOtn9yBXNU1zJszqu836mfov+uH1rrIulcroo/eiuqi6V5vEL1N8N/HZYFLSClr4ee59BDYKKKKgoKKKKACiiigAooooAKjl/wBbaf8AX5b/APo5KkqOX/W2n/X5b/8Ao5K1ofxY+qNsP/Gh6r8zuqhvP+PK4/65t/KpqhvP+PK4/wCubfyr6M+qOF04/wDEstP+uKfyFWd1U9PP/Ettf+uKfyFWc14Lo6nxMqurJrM517Sv+u7/APomSu0ribE/8T/Sv+u7/wDomSu2r1MHHlpWPpcqlzYdPzYUUUV1HohRRRQAUUUUAFFFFAHJ+K7yOx1rS5JcbGguEOfdof8ACvmDx/YHTvHOrQno03mqR3Djd/WvePjPdPawaS8Zwx84frHXjPxR+bxHYT45m0y3kJPckHn9KrIJyhm1aHScU/nHl/STNcXhIww1HFLeTlH9f0D4fqbyy8S6Ygy9xp/mKP8Arm4P9am1a1e3+FUW7q2r5/8AIbCqnwtkZfHtlCPu3EcsTD1BQn+ldZ4pgSb4bakqDizv0cfiSv8AWt8e/Y8QUO00n87Tj/kaUZ3yvEQ6qcH8m4r84nlWnEjU7Qr1EyY/76Fei3Mxt/jEZF7amufoWAP868703nVbPH/PdP8A0IV6obASfGfJUlFuhOfwGaXFrSjT/wANT/2w6cia9rVvsopv0V7nmfiGBbbxLqkCcJHdyqv0Dmul1Rf+KI8KsOnlXAz7+aa5fW5/tWvajcf89bmR/wA2Jrr9SKD4a+FkYfvN9yw/3d9enxHdYSF/51/6TM4Mjv8AWKSfZ/8ApLIvGif8Ub4PlYfOYJ1z7CTj+dQaFj/hXWt/9ftv/JqsfEJhFpvhe0Q5jj0/eMerHJ/lVPR32/D7VlH8WoW4P/fEh/pXJJP/AFfj6Rf31E1+DRrR97NH/wBfJ/nIm0djF4Y8Tyr9/wCxxx/8BaVAareDSJLLxHayH9y+mmUj/bSRCn6nH41a8NYkttftm5WXSZj9ChVwf/HaqeGlCeFfEs4HzlbaHPorSbj/AOgCsMFJLKcQmvtP8YwSfyZ05gn/AGpTX+B/Lmf+TOi8PSaWngLV49SmgZ3m/cQsg3qQoywbrzkDHbafWuG0pU/tWL03iq93vMoC5xt7Vf0XTbia4WQIQoOc4r7HhqiqeEja7TvJt7LS1l5HiZtBU8TXqSl8Xy6fn3Z7faatBbaQuSMhK8Z8Vap9v1WQqflBwK6PUbm4hsCmT0xXHWmny314FAPJ5Nerg8CqKnUg9X17Hy+UYOFCUq8mdT4I0qS7cyMp254rrtasI7S2JIAIFa3hXTItN05cgAgVy/jzXEUGGNhk8VyLETxWN5IbI8uVepjcfansYdtqGLkxoe9bZtWuItzDNcn4YtJL2/3YJGa9PktBbWfI6CuzHTjRmox3OrMpxw9VQjuczbKtrN6VoXF0ZYNoNc/fXoF3sU85ra0+3aaEMRWdWFkpyMa1OyVSZmwWuLveR3rq4rnZa4z2rHuUEDZ6URzmRdoNZVV7VJsxrr2yTZl6opnvM+9dXoH7q3HbisJ7bMm4itazlEUeKWJfPSUELFPnoqCG+I5fMQiq3huHZIGp1/8Avj61Z0tfKxWfw4fkRlfkw3Ije1WX/RCM9q4qzh3alu966PUbjdBjPasrTov9I3Ed6ywkfZUmjPCXp0pPudvakJZKT2FcNrv76/A9661p9lpt9BXKXC+bfZ965svp8lSc2Y4F2m5djq/D6COyApNaOVAqbSgFtgPaqeqOXeuKEebGORhzXn8xdGT5810sXSsLSI8Ln2rej6V4mfzvOx24LWq2TClpBS18bLc+gjsFFFFSUFFFFABRRRQAUUUUAFRy/wCttP8Ar8t//RyVJUcv+ttP+vy3/wDRyVrQ/ix9UbYf+ND1X5ndVDef8eVx/wBc2/lU1Q3n/Hlcf9c2/lX0Z9Uef6ef+Jba/wDXFP5CrOap6ef+Jba/9cU/kKsZrl9ifm8qvvMs6ef+J/pX/Xd//RMldxXA2cMNzrulRzxJKn2hjtdQwyIpCDg12V7pOm6mUN/p9pdFMhDPCr7c9cZHFb048sbH2OSS5sIn5sZBrWnXENzLHcfu7ZPMlZ0ZAE5+cZAyvythhkHB5qJfEWmM8KefIryv5YVoJFZW3KuHBXKZLKBuxnIxUMXhm0itri2+0XLQTwi2ZGK4EIDhYxhc4HmMQfvdMkjim/8ACMWxnWd7u7kmMgkmdigM5DIy78KBgbFxtxx1zmrPXLB8Q6csMkrPcKkT7JS1pKDHwDlgVyq4OdxwPerVtqNtd3VxbRGQTQYMiSROmASQCNwGQSrcjI4rLXwuot1hbVtQf/SftMjP5JMz9t/7vBAwMDtgdgAL9hpYsLq8uDd3Fw91JvbztnyDnCqVUHaAcAHP5kkgF+iiigAooooA8n+NsTzWmkqmf+WxOPrHXj3xQyuv6dEfvRaXbofY4J/rX0D49shqGo6TCRn93O30+aL/ABr53+KEyy+P7+NPuQCOJfoEH+NPIvfzir/dh+bh/kzpxteMsFQoLdSk/uVv/bg+GKN/wnVpMuQIIppSQOgEbD+tdLdT+b8M/Ezl9wa4g2jOcDfwK5vwOWs9P8Saopw0FgIF+srhf5A1P9qP/CstXiLctqFuD9MMf/Za6cy/e53Qa+xyfem5fqi8OrZViX1coL5Jxf8A7czk9JBbWbEDqbiMD/voV6z/AGkifGTywBg3H2cn3Ix/OvLfDq7/ABPpKf3ryEf+PiupluHf4mNcZ+Y6vu/8i1HFkeb2a/uz/wDbTXIVerWXeKXyblf8jj9Wi8jWb6H/AJ53Ei/kxFdZ4mfZp3hu1X/Vx6TE+P8AacsxP8vyrnPE/wDyNesY/wCf6b/0M1u61J9p0Xw3deuneSfrHLIv8sV6PEd5YenL+/8A+2yOPh9t4iEn/I/v93/gh4vf7T4V8K3J++IJ4G/4BJx+hFV9CHmeCteHaO5tZP8A0Yv9af4h/wCRH8O5/wCe13j80pvhw48HeJfdrQf+PvWP/MhX9bVNPyLXuZjp/P8Am9fzLPhxvLt9fl/uaTP/AOPFU/8AZqpeECbi11/T84WbT2mH+9E6yfyDfnVzw5h4tcgP/LbSZwPquH/9kqn4RDQWev6gACIdPMI57yuqf+glvyrlwFv7LxN97u3ryxt+J0Zmpf2nTt/c/wDS2dR4M0Cz1PSpbm4QF0nKAkdgqn+tdQ+n2VjEdoUYrzPTfEk+j2rW0WcM5f8AMAf0pl14ovbrI3YBr7TJsHXq4Kk+b3bHxGeZbi6+Z1Zc3u30Og1u6hdiikU7QLaJHDkCuWtY7q8mDNnFdXZxSW8Q+le5WpqlT9mmcmIpKjS9kpHSajrq2dkQrAcV5ZeTT6vqRxliTgVvao0tyfLGcVp+GdAVZVkdaihGlhKTm92Z4VUcBSdV/Ezo/Bfh9bO1WR1569Ks+LNRjs7ZlBAOK1Zb6LTrEgEDArybxNq8mpXphQkgnFeTg6NXGYt1qmyPJwVGpj8V7WexHpSvqep7uSM16la2At7IEjGBXOeCtC8tFldfzrqddvY7KzYAgYHSrzHE+0xCo0ys0xHtq6o0tkcVrl4BOVU96t6NAZUDEVzIkfUtU+Xkbq9F0qw+z2gJHIFdGKaoUlHqbY1rDUYw6mfeRrElUYJS74BqTXbsIxUGmaLbmchsVlFWpc8jGC5aPPIu+Txk06Nghq7cxiGPmshZN02BWUHzq5hB+0Vy3NmSnWkex6mjiOzOKTIRqjm0siHLTlRbnl/c4zWZHFun3Ed6s5Mhx2pyoFcVMfdViIe4mjZtG2Q/hVO6+eWp4W/dioyu6UVxwjyzcjjjpJs09NTbEK1U6VStE2xqKvp0r47OavNUZ6uXR6jxS0gpa+alue7HYKKKKQwooooAKKKKACiiigAqOX/W2n/X5b/+jkqSo5f9baf9flv/AOjkrWh/Fj6o2w/8aHqvzO6qG8/48rj/AK5t/KpqhvP+PK4/65t/Kvoz6o8508/8S21/64p/IVYzVSwP/Eutf+uSfyFWM11exPyadX3mW9MP/FQ6X/12f/0TJXeVwOlHPiLS/wDrs/8A6JkrvqwqR5XY+94dlzYJPzYUUUVB7gUUUUAFFFFABWFq+jXF/q9ndR+QUhMZ3SEh4dsquxTg5LgbTyOMdelbtFAGJe6fbXHiqxklV3P2Sc4MrbQQ8OMLnA/Ac1j6h8KPBWq3819e6L5tzM26R/tUwyfoHAroZ/8AkZbH/rzuP/Q4a0qKf7qbnT0k92tG/ViaTtfocfa/CzwZZ2V1Z2+jbLe72ecn2qY7thJXkvkYJPSkPws8GHT5LD+xv9GkkWV0+1TcsoIBzvz0Y/nXY0U3JuftG/e79fvLUpKDgno+nT7vkjiLX4Q+BbO7hurfQ9k0LrJG32uc4YHIOC+Ooqz/AMKw8HDUft/9j/6V5vnb/tM3385zjfjr26V11FFWTrfxfe9dQpzlTbdN2v20OIufhD4FvLqa5n0PfNM7SSN9rnGWJyTgP61PJ8LfBktnbWj6NmC1DCFftM3yhm3Hnfk8k9a7CiqqVJ1Fy1G2vPUmn+7acNLdtDj7n4WeDLuxtrKfRt9vbFzCn2qYbSxBbkPk5wOtNt/hV4LtLO5tINF2wXJQyr9qmO4rkryXyMZPSuyopc8vZ+yv7vbpvfb119R3fNz9d7nH23wt8GWbu8GjbGeN4mP2qY5VlKsOX7gmm23wq8FWlpc2sGi7Ibnb5q/apju2nI5L5HPpXZUUk2ouC2e66FSqTnJTk22uvU8B+IHw10uy1+CPRrU21sbVWZPMd8vubJyxJ6AflWJaeBoYyC4yfevR/iTrUGneIbeKQjcbRW/De4/pXDT+MIEB2kV91ltTF/U6cKe1j4bNMTmE8XUjC9rlhNEt7ReFHFUb1oowVGKybzxaZiQh/KqCXNxePnBwa9OnhavxVWc9LB1/jrM0IkWSfJ9a6O0uEtYu3FYdpauqgkVLcOwXaDU1Yqo+W5FeKqvlvoV9f1iSYGNCeeMVS8O6C91diWUZyc81attONxOGf1rs9OhisogcAU62IWHpezp7sdfFRw1H2VHdmvAsOm2I6Aha838Way1xOYY2zk4Arb8Ra9siaNG5xjiuV0jTZdSvxLICQTXPl2FVK+IrbmOW4ZU08TWN3wboZdllkXrzXb6pKllZYGBxUmmWken2YyAMDrXG+Ltb3ExRt7ACuHnnj8Zp8KOJuePxV+n6GFdXDX+pbRyM13uiWHk2oYjtXJ+E9Ia5nE0g6nNeh3JSzsjjAwK2zTEKMo4enubZnVV1Qp7I5rXboRgrnmqGjwtO+4is3Urpry/KKc812Ph6w8uBXYVdZrDYe73FUX1fDpdWLMnkRVitMXmwPWtjXZ1iQqDzWJpUTXFxux3rLD60vaSOfDxtTdSRtW1viMMRUUpAlCiteRBBa/h1rCR/NuvbNYUp+0vLoYU25tyZrQj92Klhj3Sj60sceIxVm2j+bdXHVrKKbOe92X4RgCrS9KgjFWBXw2Pqc02fQ4KHLEcKKKK8lnpoKKKKQwooooAKKKKACiiigAqOX/W2n/X5b/8Ao5KkqOX/AFtp/wBflv8A+jkrWh/Fj6o2w/8AGh6r8zuqhvP+PK4/65t/KpqhvP8AjyuP+ubfyr6M+qPMLCUf2dbf9ck/kKseaKybGT/QLb/rkv8AIVP5lfRLDaH43O/Mza0Z93iTSx/01f8A9EyV6FXmnh593ifTB/00f/0U9el15GNhyVbH6Jwz/uC9WFFFFch9AFFFFABRRRQAUUUUAZs//Iy2P/Xncf8AocNaVZs//Iy2P/Xncf8AocNaVABRRRQAUUUUAFFFFABRRRQAUUUUAeDfGzSrq+8Z2ckBbaNPRTj18yT/ABrz+Hwrdyfe3V7t47MA12Dzdu77MvX03NXMm6s4x1WvvMuzKrTwlOnFbI+HzLNq1LFTpwjszhbPweykFga3rfQ4rdRkDitC41q1iHDCsS88SIchDXV7XFV3qeb7bGYl6l24EcKkDFY8rq78VRk1Ga6f5c1ctbaRsM2a3jS9mryZ0xoukrzepoWhVADS32pFIyqmo2jKLVX7MZn56VkoxcuaRlGEJS5pGbHaS6jc5YEjNd5omlx2cQYrg1S02zjhAJAq3fakttCQCOlYYqtOt+7hsc2MxE67VKnsLr2srbQFFbnHrXA29vNq2oByCVzVi7km1O72jJBNdn4c0VLeNXZRnrWqdPAUL9Tpi4YChf7TNfRbBLG0HABxWJ4p1lUjaJG9q2dX1GO0tWVSAcV523m6tqHGSua8/L8O61R4mqceCo+1m6tTZal7w7p73l15rKTk5r0XCWdr6BRVHQ9NWytVOPmIqp4k1NLeAxq3OK58VWeNxSpQ2RliKjxVb3duhzetX5ubvy05ycV03h2wEdusjDmuW0KwfUb3zSCRnrXokcaWtuFHAUVtmldUaaoQ3Zri3GEVRj03MnXLoRRlB2rO0iEyPu9TVbVbk3d7sU5Ga6DSLXyoAxHaplbDYVJ7nPJezopdWXfK2oBU0KYAoIqWNa+fxFZqmzGjC87E6DipRTEFSCvkcRO8mfTUI2iLRRRXGdQUUUUAFFFFABRRRQAUUUUAFRy/620/6/Lf/wBHJUlRy/620/6/Lf8A9HJWtD+LH1Rth/40PVfmd1UN5/x5XH/XNv5VNUN5/wAeVx/1zb+VfRn1R4tZy/6Db/8AXNf5VP5tZtpJ/oUHP/LNf5VN5vvX20YaI/JZ0/eZu+HWmk8T6aIJER/Mc5dCwx5b54BH869LvY9SkKfYLu0gAzvE9s0ufTGJFx+teY+EX3eLtNGf4pP/AEU9etV81mitiH6I+94dVsCl5s5aztdeisNQjIuRczW/lxySTq6i42yFpVyx2xk+WAoHH90DJqBbDXPtVuyLfJbCcNFHJeBmhTzIy3mnefMBUSYGWxnHHbsKK84904E6X4s/swRme+aYOHkIuUDPL5bZ2ndgQb9vy8NgnjpXU2NjPFreoXcj3CwuEjije4Z0JwWZ1UsQuSwXGB9z0NatFABRRRQAUUUUAZs//Iy2P/Xncf8AocNaVZs//Iy2P/Xncf8AocNaVABRRRQAUUUUAFFFFABRRRQAUUUUAeF/Gi8u7fxhaJBnabBCcevmSV50brUpRwGr3Lx/pUN9r8EsgG4Wyrz/ALzf41zS6JaR/wAK193l2YUqWEpxcdUj4jMczo0sVODhdpnmKWOoXLfMW5rWs/Dc74Lg13Ys7WHkBaJLq2hXqtdc8znPSETzamb1J6U42MG20JIVBYc1aMUcS44pt5rUS5CkViy6jJOx25qYwq1NZEQhXq+9MvTOpbFLCADmqkMUrnJzV1Yyq1UkkrFySirXJpLvy04rGuHlvJNozir0ib2xVq0tFDAkURcaav1CMo0lzdR2jaQseHZea6Ga8S0g4OMCqJuEt4657Ur+S5fy0JNcjpSxM7z2OP2c8VUvLYg1a+k1C48tDnJroPDejCMLI61Q0XRiziR1JOa7LdHZW3YYFTjsTyQVCj1LxeIjGKoUthb68Szt2YkA44FebajdS6rqQRSSua0dd1SS5lMUZzk1f8NaF8wnlHvk08LRp4Ci6s92aUIrDU/az+J7G74e04WdmpIwSKXXr8W9qyg/Ma05ZEtoCTwAK4fUrl7+92KeM15WDhLF4l4ieyOWnF1J+96sfo9s11deY3c120aCOMKO1Z2jWAt7dSRya1G6VhmeK9rV5I7ImpL2knPp0GgZNWIxUSDJqygrwcdWsrHRgqV3ceop4pop1fN1HdnvQVkFFFFZmgUUUUAFFFFABRRRQAUUUUAFRy/620/6/Lf/ANHJUlRy/wCttP8Ar8t//RyVrQ/ix9UbYf8AjQ9V+Z3VQ3n/AB5XH/XNv5VNUN5/x5XH/XNv5V9GfVHz/aS/6HBz/wAs1/lU3m+9ZtrL/okP+4v8ql82v0KMNEfmkqfvM6vwU+7xjpw95P8A0W9exV4r4Dfd4008f9dP/RbV7VXyWcq2Kfoj7PI1bCJebCiiivKPYCiiigAooooAKKKKAM2f/kZbH/rzuP8A0OGtKs2f/kZbH/rzuP8A0OGtKgAooooAKKKKACiiigAooooAKKKKAPIvinrjaZ4ntoR/FZq//j7j+lcLJ4qlYYGa9C+JmhjU/ElvOR920VP/AB9z/WuQTwmnUgV9zl1TCxwtPn3sfC5lUwccXPnWtzAk1+6lPG6ojNfXH97muxi8NQR4yBV5NKtoh90V2PH0I/BE4JZhh4fw4nE22k3E5BkzW7aaKsYBYVtsbeBewxWdc6tDGCARWUsTVraRRhLF1q+kVoK0EcQwBVKWRQcA1Vl1RpmwtNiikmbJzVxpuKvMuNGUVebLEeGaray+WtMS22ikkQ9KltNkScZOxXuJXlO0ZqfTtL3yB2H51NbWwLZatVZI4I/pUVKriuWBFWu4x5IF2MRWkP0Fc5rWsbv3cbZJ44pNQ1J5Moh61W07SmuJhJJk896zoYeNP97V3JoUI017WqGj6S9zcCWQZye9d5DGlvCFAAAFV7O1S1iHGDVLVdSEUbIhrzcVUnjqqhH4UY1K0q1S/wBxQ13Vc/uozUOg6a00nmyDryc1Ts7SS/uwxBIJrt7O1S2hVFFbY3EQwND2cNzZrlj7KG73JlUIoApp5NOY0iDJr5iMrJzkZyV2oIkjWp1FNRakArwsXW5pHs4alyoUUtFFec3c70rBRRRSGFFFFABRRRQAUUUUAFFFFABUcv8ArbT/AK/Lf/0clSVHL/rbT/r8t/8A0cla0P4sfVG2H/jQ9V+Z3VQ3n/Hlcf8AXNv5VNUN5/x5XH/XNv5V9GfVHzPby/6LF/uD+VSebWfBL/o8X+4P5VJ5tfp0Kfuo+ElT1Z2nw8fd44sB7Sf+i2r3KvAvh3crF43sHcSEYkGERnPKMOigmvcb3U4NPKCaO7YvkjyLSWbp67FOPxr4rPlbGP0R9RlKthkvNlyisC28TCXT768mtNiW1mL5Akm4vEwfbnIG1v3ZyOQMjk1JDrly15a201lEjSXb2kxS4LBGERlBX5BuBAwc7cH1rxj0zborEOuXK2fntZRBk1BLORRcEgBpFjDg7OT8wOOO/NXbO+kuL+/tZIEj+zOoVlkLb1ZcgkYGD7c/WgC9RRRQAUUUUAZs/wDyMtj/ANedx/6HDWlWbP8A8jLY/wDXncf+hw1pUAFFFFABRRRQAUUUUAFFFFABRRRQB5v8QNVjsdegicjJtVb/AMeb/CuRfxJCO4rW+K2lzXvie2ljzgWSLx/vv/jXCr4euWODur7bLsPh5YWEpvWx8FmWGw0sXUlOWtzam8ToOhrOn8SyOcITUkPhaR+oNXofCoXBIP5V23wdM4b4Gn5mG1/dXXTPNLHp1xO2WzzXX22hwxDlatGCCFf4azlmEE+WmjKWYRjpSic9ZaLgAsOa0haJCvYVPNfQxjAwKyrnU95IU1mnVqu7MOatWd2TyyKvAqHduaqqb5TnmriQHHNaOKiauKhuSrKEXiqtxM8hwM1Y8qnpbgnOKlOK1IUox1K9pYeY4ZhXS2cCQJnAqnbxhBTri68lMg81x4iU6z5UctWpKrKxLf34iQgGue2S38/GSM09hLez4GcZrpdL01beMMy/NSnUpYGld7m0F7JWWsmP0zTltYgSPmrRJwKXoKjJya+VqVZ4qo5zNHamvNidTUyLTUWp1WuPGV0lyo3wtFt8zHKKcKBS14NSV2e1CNkFFFFZmgUUUUAFFFFABRRRQAUUUUAFFFFABUcv+ttP+vy3/wDRyVJUcv8ArbT/AK/Lf/0cla0P4sfVG2H/AI0PVfmd1UN5/wAeVx/1zb+VTVDef8eVx/1zb+VfRn1R8nwyfuI/90fyqTzKoRSfuU/3RT/Mr9ahD3UfJunqd78Ln3ePrEf7Ev8A6A1fQNfO3wnfd8QbIf8ATOX/ANANfRNfBcRK2Ofoj3suVqFvMqpplhGWKWVspZzIxWJRlyCCx46kEjPoT60waNpYFsBptmBauXtwIF/csepTj5T7irtFeEdxntoOjtC0LaTYmJ5fOZDbptaT++Rj73v1qxDYWdvdTXUNpBHcT482VIwHkx/eI5P41YooAKKKKACiiq8t9bwXlvaSM4muNwi/dsVJAJI3YwDgE4J7UAVZ/wDkZbH/AK87j/0OGtKs2f8A5GWx/wCvO4/9DhrSoAKKKKACiiigAooooAKKKKACiiigDiPGMsCatEJcZ8gfluauZa6tVPak+Kktwnii2EOdv2Jc49d71wxe+bsa+vwGAU8PCbluj88zTA8+NqSct2d3/altGOq1TuNdhXoRXIiC+k9akTSLuXrursjl9CDvJnAsFRj8UjYm8QA5Cms2bVppj8uat23hyQgb+TWrbeH0TlgKp1cLR2DnwtL4dTm0juLgjOa0bfSmbG4V0a2MEAyQKhmuYYhhccVk8a6mlNGU8ZKWkEVo7JYl5FRzOqHio5tQLnC0yOGSZskUKMt5kKMvimxUJduKvRRkDJFPgtlQAmieVUGBWcp8ztEylPmdohJKI1681QKyXcu0A4zViK3kupAMcVu2enpAM4y3rWFfF08LG8tyoKztHVkOnaYsChmHzfyrUwFFHCio2bNfLVq9TF1Ly2OtRjRV3uDNmhVzQq5NTItZ16saUeVCo0nUldjkWpAKQCnCvBrVXJnt0qfKrC0UUVynQFFFFABRRRQAUUUUAFFFFABRRRQAUUUUAFRy/wCttP8Ar8t//RyVJUcv+ttP+vy3/wDRyVrQ/ix9UbYf+ND1X5ndVDef8eVx/wBc2/lU1Q3n/Hlcf9c2/lX0Z9UfHEcv7pP90U7zaopL+7XnsKd5vvX7PCHuo8J09T0b4QPu+IlmP+mUv/oBr6Rr5n+DT7viRZj/AKYy/wDoBr6Yr854nVswfoj1MIrU7BRRRXzx0hRRRQAUUUUARzQRXMLwzxJLE42ujqGVh6EHrWaPDtjDc2ctkosI7aVpvItIo0SRipUlhtz0JHBFa1FAGTLBCniy0nWJBK9lOHkCjcwDw4BPetas2f8A5GWx/wCvO4/9DhrSoAKKKKACiiigAooooAKKKKACiiigDiPGNjFc6xE7gFhAo/Dc1YH9mWq/w1Y+I2oS2niC3RM4Nqp/8feuQOsXJHQ19Pg8LXnQi4zsrH5tm2HqyxtRra51H2S1j5+Wnh7WMfw1yB1C7foGxTd17J/erp/s+T+OZ531SX2mjqpdUt4h8uKzp9dUfdNZMemXUx+Ytir0OgscbgatYfDUvidw9nQh8UrleXU5pj8ueabHbz3DZbNb1vo0cYG4VdWCGEdBUTx1KHu00J14xVoIx7fSgACw5q+sccQwAM0+a5C5C1XSKW4boQKyc5zXNUdkc8pSnq2MlmJyEGT7U62095m3SdK0rfT0TkjcauqgUV52JzWFJclLc6KVCUvJEUFskKgAYqbOBSFsVGWJNeK/aV5c02dN4UlaIrNmhVzQq5qZUoq1o0o2iFKlKpK7BUqUCgCnAV4deu5M9ejRUUAFLRRXI3c6krBRRRSGFFFFABRRRQAUUUUAFFFFABRRRQAUUUUAFRy/620/6/Lf/wBHJUlRy/620/6/Lf8A9HJWtD+LH1Rth/40PVfmd1UN5/x5XH/XNv5VNUN5/wAeVx/1zb+VfRn1R8RKfkX6U7NMU/IPpS5r9shJcqPNaPRPgof+LlWn/XCb/wBBNfT9fLfwZiW4+I9ojlwPJlPyOyHhc9VINfTF7pkGoFDNJdqUyB5F3LD19djDP41+bcUu+YP0R2UPgJIb60uBMYLqCUQsVl2SBvLI6hsdD9aZFqunTm3EN/ayG4LCAJMp83b97bg84xzjpWRD4YeOxu7Nr1WiuLUWRIhwywqsgXnccvmTJbocfdGc1KuhXgurS5N/AZI71rycC2IEhMRi2qN/yfKep3ZPPtXzhsXl1vSXjSRdUsjHJJ5SMLhMM+M7Qc8nHarMV3bTzSww3EUksJAlRHBZCemQOn41zr+ES9pNAJbCPzZd6+VY7RAdoXdEN/ySd93Pbjjm/o+gDSb68uFnVluHd9qx7TlpHclzk7iN+AcDAGOaANmiiigAooooAzZ/+Rlsf+vO4/8AQ4a0qzZ/+Rlsf+vO4/8AQ4a0qACiiigAooooAKKKKACiiigAooooA4Xxnp0d3rUMj9Rbqv8A483+Nc//AGNBnt+VdN4uuDFq0SgZ/cA/+PNWD9qOM7DXt4eeJVKPK9D8vzibWOqq/UjTSrde1TLZQJ0UVGLmVvurSFrh+AMVo1Xb96djy+a5aCRp0UCmyXEaDrz7VALWZvvNUqWAPLEmsn7GLvUncajKWyIGvmbhFpnlXEx5BAPrWlHaonRRU4QCsJ5nRpaUonRDCzlq9ChDYKvLcmrqRKo6U/IFMZ686ria+IfkdMaVKlq9WPyBUbPTSSaULmlGjGGsxSqynpETk09Up6pUgWuevjElaJtRwrbuxFWngUoFLivIq1nJnq06SiGKWiiuZu5ulYKKKKQwooooAKKKKACiiigAooooAKKKKACiiigAooooAKjl/wBbaf8AX5b/APo5KkqOX/W2n/X5b/8Ao5K1ofxY+qNsP/Gh6r8zuqhvP+PK4/65t/KpqhvP+PK4/wCubfyr6M+qPh5T8o+lLmo1Pyj6Uua/W4VvdRwWPSPgif8Ai5dp/wBcJv8A0GvqGvlz4IH/AIuZa/8AXvN/6DX1HXwPEUubHN+SOqj8IUUUV4RqFFFFABRRRQAUUUUAZs//ACMtj/153H/ocNaVZs//ACMtj/153H/ocNaVABRRRQAUUUUAFFFFABRRRQAUUUUAcp4miV9SjYgZ8kfzNY/kL6CtrxK2NRjH/TIfzNY3mUc9faL0PzzNI0frlRy3uAhUdqURgUnmUhc0uWvLdnBejHZEmFFGQKi3Gk5NH1Zv4mH1hL4USlxTC9N2k04JVKFKG5LnUnsNyTShSakEdSBKyqYyMdImsMLKWsiJUqVUp4WlxXm1cXKR30sMoiBadijFLXDKo2dkYJBRRRWZoFFFFABRRRQAUUUUAFFFFABRRRQAUUUUAFFFFABRRRQAUUUUAFRy/wCttP8Ar8t//RyVJUcv+ttP+vy3/wDRyVrQ/ix9UbYf+ND1X5ndVDef8eVx/wBc2/lU1Q3n/Hlcf9c2/lX0Z9UfDKn5R9KXNRg/KPpS5r9GjW91HJY9J+Bx/wCLm2v/AF7zf+g19S18s/A0/wDFzrX/AK95v/Qa+pq+PzqXNi2/JG9P4QoooryTQKKKKACiiigAooooAzZ/+Rlsf+vO4/8AQ4a0qzZ/+Rlsf+vO4/8AQ4a0qACiiigAooooAKKKKACiiigAooooA5XxMudSj/64j+ZrF2Gt/wAQrnUI/wDrkP5msnZXPPG8j5ex8HmWE58XUl5lfZSiOp9lLtrN49nGsEiAR04R1NtpcVzzxsmbRwkURBKeFp2KXFc08RJnRGgkNApcU7FFYOo2bKCQmKWiiobuWkFFFFIYUUUUAFFFFABRRRQAUUUUAFFFFABRRRQAUUUUAFFFFABRRRQAUUUUAFRy/wCttP8Ar8t//RyVJUcv+ttP+vy3/wDRyVrQ/ix9UbYf+ND1X5ndUyaPzYJIwcb1K59Min0V9GfVHwteWdxp19cWN3H5dzbStDMmQdrqSGGRwcEHpUFe3/HzwS8N6njCzXMM+yC+QKxKuBhJCeQFICp2wQvUsceIE4Ga+wwmIVakp/ec8lZnpvwLsdQuPH/2yzt0khtLZzO8jsiqG+UDcFb5jyQDjIVueK+l729ntSgh0y7vA2cmBoht+u91/TNcr8LvBCeCfCccMy51O82z3rFVyrlRiIEdVTkdTyWIxnFdtXzWNrqtWc1sbRVkctZ+J7iWw1CeRbaSaG386GGPIbzdsjG3bk5kUR5OMHnoKZF4nme7s7YX+mSyTXRiTYhVbmMbd0iHzDjG4rj5ssMDGeOsorkKOF/4TLUP7N+0f6Du8vzN2xtu/wArf9m+9/rc8Z/8dzxW7pGs3F/q15ay+QUiaQbY1IeHbKyKHyTkuBuHA4z161u0UAFFFFABRRRQBmz/APIy2P8A153H/ocNaVZs/wDyMtj/ANedx/6HDWlQAUUUUAFFFFABRRRQAUUUUAFFFFAHO6+P9OT/AK5D+ZrKxWrr3/H8n/XIfzNZdeBipP20j5fGRXt5eomKMUtFc/Mzm5UGKKKKLhYKKKKQwooooAKKKKACiiigAooooAKKKKACiiigAooooAKKKKACiiigAooooAKKKKACiiigAooooAKjl/1tp/1+W/8A6OSpKjl/1tp/1+W//o5K1ofxY+qNsP8Axoeq/M7qiiivoz6ooa5o9n4g0O90i/TdbXcRifABK56MuQQGBwQccEA14B8MfhbqMfxGuH1uNFg8PTqzbdxWebG6LYwxwAVkOeeVBX5jj6OoreliJ04ShHZiaTCiiisBhRRRQAUUUUAFFFFABRUc08VtC808qRRINzu7BVUepJ6VmjxFYSzW62ji8hnl8kXFrLHIivjO04bdnAzwDgcnjNAEk/8AyMtj/wBedx/6HDWlWFf6jBa+KLJJI7piLOf/AFVrLIOXh7qpHarf9uWn/PHUP/BfP/8AEUAaVFZv9uWn/PHUP/BfP/8AEUf25af88dQ/8F8//wARQBpUVm/25af88dQ/8F8//wARR/blp/zx1D/wXz//ABFAGlRWWPENgZmiCXxkVQzL/Z8+QDkA/c74P5U7+3LT/njqH/gvn/8AiKANKis3+3LT/njqH/gvn/8AiKP7ctP+eOof+C+f/wCIoA0qKzf7ctP+eOof+C+f/wCIo/ty0/546h/4L5//AIigDP17/j+T/rkP5msurGralDd6htggvnZIlLL9hmBAJbBwU9j+VU/Nf/ny1D/wBm/+Jrw8TRqOrJqL+4+dxVCrKtJqL+4koqPzX/58tQ/8AZv/AImjzX/58tQ/8AZv/iaw9hV/lf3HP9XrfyP7mSUVH5r/APPlqH/gDN/8TR5r/wDPlqH/AIAzf/E0ewq/yv7g+r1v5H9zJKKgS68xpFS1vmaNtrgWU3ynAOD8voQfxp/mv/z5ah/4Azf/ABNHsKv8r+4Pq9b+R/cySio/Nf8A58tQ/wDAGb/4mjzX/wCfLUP/AABm/wDiaPYVf5X9wfV638j+5klFR+a//PlqH/gDN/8AE0ea/wDz5ah/4Azf/E0ewq/yv7g+r1v5H9zJKKgjuvOUtHa3zqGZSRZTHkEgj7vYgj8Kf5r/APPlqH/gDN/8TR7Cr/K/uD6vW/kf3MkoqPzX/wCfLUP/AABm/wDiaPNf/ny1D/wBm/8AiaPYVf5X9wfV638j+5klFR+a/wDz5ah/4Azf/E0ea/8Az5ah/wCAM3/xNHsKv8r+4Pq9b+R/cySioIrrz4Ulitb543UMrLZTEMDyCPlp/mv/AM+Wof8AgDN/8TR7Cr/K/uD6vW/kf3MkoqPzX/58tQ/8AZv/AImjzX/58tQ/8AZv/iaPYVf5X9wfV638j+5klFR+a/8Az5ah/wCAM3/xNMluvIheWW1vkjRSzM1lMAoHJJ+Wj2FX+V/cH1et/I/uZPRUfmv/AM+Wof8AgDN/8TR5r/8APlqH/gDN/wDE0ewq/wAr+4Pq9b+R/cySio/Nf/ny1D/wBm/+Jo81/wDny1D/AMAZv/iaPYVf5X9wfV638j+5klFR+a//AD5ah/4Azf8AxNMkuvJUNJa3yKWVQTZTDkkAD7vckD8aPYVf5X9wfV638j+5k9FR+a//AD5ah/4Azf8AxNHmv/z5ah/4Azf/ABNHsKv8r+4Pq9b+R/cySio/Nf8A58tQ/wDAGb/4mjzX/wCfLUP/AABm/wDiaPYVf5X9wfV638j+5klRy/620/6/Lf8A9HJR5r/8+Wof+AM3/wATUckrebZ7rW9QG8thue0lVR++TqSuBWlGjUVSLcXuuhrQoVVVi3F7ro+531FFFe+fShRRRQAUUUUAFFFFABRRRQAUUUUAFZ2oWVks41idJ2ms4nKGOZwAvUjaCFOcDqOcCtGigDnotfnmiQjT4ftxuRaohuP3eTCJ8+ZsyBt/2eo9Oa2bC7TUNOtr2NWVLiJJVVuoDAEA/nUA0TSVt3t10uyEL43xi3Ta2CSMjGDgkn8TV4AKAAAAOAB2oAWiiigAqpqV7/Z+ny3Qj8xkwFTdtBYkAZPYZIye1W6ZLFHPE8UsayRupV0cZDA9QR3FAGDNr91DZQX66Urx+dJBeyC4AFuI5ChYfLlxkMQAAa6GqI0bS1EYGm2YEb+YgEC/K3A3Djg/KvPsPSr1ABRRRQAUyUyCFzCiPKB8quxUE+5AOPyNPpksUc8TxSxrJG4KsjjIYHqCO9AGBN4l8vS7HUUsVZrmwN9KplwUiRVZgDt+Zh5gwOAeeRXRVRj0TSooI4I9MskijcSJGtuoVWAADAY4OABn2q9QAUUUUAFUdYvpNM0e7v4oEnNtE0pjaQpuVRk84POB6VeqG6tLa+tntry3iuIHGHimQOrfUHg0AUL3Uriw1a1h/s8GyuWCyXnmgbJDwq7MZOcKM9Oa1apx6TpsUsMsen2qSQ58p1hUGPJJO044zk9PU+tXKACiiigArP1fUJtNt45YbeOcs4Ty2lKO5PQIAp3N14OBxkkDJrQqtd6dZX5iN5Z29wYjuj86JX2H1GRwaAKMesY1ttPW1VYROYPND/MZfKEx+XH3drfeznPbvWvVaLTrKCZZobO3jlWMRK6RKGCDooIHTgcdKs0AFFFFABWffahNaX1nBHbxzLcPsI80iQDuwXaQVA5JJH4kgHQqpdaVp19Ms15YWtxKqlFeaFXIU9QCR09qAKunalcS6nd6dc6eLMQDdbkShxNFuK7sAfL0HB55rVqCCxtLaaaa3tYIpZjuleOMK0h55Yjr1PX1qegAooooAKzLvUXh1m006S1R7e7DKJSzcEKzEEFdp4U8bs98YBI06qy6bYT3S3UtlbSXC42yvEpcYORzjPBANAFPRdabVt2+2EOYIrmPEm/Mcm7bngYb5Dkc445Na1Q29na2nmfZraGHzGLv5aBdzHqTjqfepqACiiigArIl1Sf+159O+wJIy27XEB8wjzSpXC/MgUHLDlS2OM4yM69U20nTXuHuG0+1aZ87pDCpZsgg5OMngkfQmgCPRtRl1KyaS5tfslzHI0c1v5m/y2HbcAAeCDx61oVFbWtvZwLBawRQQrkiOJAqjJyeB71LQAUUUUAFc1B4ttr+PUDbQwzi1uYIIh54PmGR1VWbAOwbjnucDPXiulqtLp9lO8jzWdvI0ilHZ4lJZTjIORyOBx7D0oAbpt6dQsUuDH5bFmR0DbgGVipwcDIyDg4GRVumRRRwQpDDGkcSAKiIoAUDoAB0FPoAKKKKACiiigAooooAKKKKACiiigD/2Q==", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "step = 0\n", "while t <= T:\n", " step += 1\n", " phiVec[:] -= pVec[:]*dt/2\n", " pVec[:] += S2*phiVec[:]\n", " t += dt\n", " # set the values on the Dirichlet boundary\n", " p_in.value = numpy.sin(2*numpy.pi*5*t)\n", " op.setConstraints(p)\n", " phiVec[:] -= pVec[:]*dt/2\n", "phi.plot(gridLines=None, clim=[-0.02,0.02], figure=plt.figure(figsize=(10,10)))" ] } ], "metadata": { "jupytext": { "cell_metadata_filter": "-all" }, "kernelspec": { "display_name": "Python 3 (ipykernel)", "language": "python", "name": "python3" }, "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.10.12" } }, "nbformat": 4, "nbformat_minor": 5 }