commit 4f3ecfeed4a8277de6b1be71e183afd75dc72a40 Author: Laureηt Date: Tue Mar 15 08:04:08 2022 +0100 ptdr diff --git a/TP-EDP.ipynb b/TP-EDP.ipynb new file mode 100644 index 0000000..79866aa --- /dev/null +++ b/TP-EDP.ipynb @@ -0,0 +1,420 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": 1, + "metadata": {}, + "outputs": [], + "source": [ + "import matplotlib.pyplot as plt\n", + "import mpl_toolkits.mplot3d\n", + "import numpy as np\n", + "import scipy.sparse as scps\n", + "import scipy.sparse.linalg as ssl\n", + "import math" + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "metadata": {}, + "outputs": [], + "source": [ + "def maillage_carre(n: int):\n", + " \"\"\"\n", + " Une discrétisation possible d'une EDP elliptique sur le domaine ]0,1[ x ]0,1[.\n", + " Le carre [0,1]x[0,1] est maille uniquement avec des triangles.\n", + " Les conditions limites sont de type Dirichlet uniquement -> `neumann=[]`.\n", + "\n", + " Args:\n", + " n (int): nombre de points par cote du care => Npts points de discretisation au total\n", + "\n", + " Returns:\n", + " coordinates : matrice a deux colonnes. Chaque ligne contient les coordonnes 2D d'un des points de la discretisation. Ces sommets seront identifies a l'indice de la ligne correspondante dans la matrice coordinates.\n", + " elements3 : matrice a trois colonnes. Chaque ligne contient les indices des sommets d'un element triangle, dans le sens antihoraire.\n", + " dirichlet : vecteur colonne des indices des sommets de la frontiere de Dirichlet.\n", + " neumann : matrice a deux colonnes. Chaque ligne contient les indices des deux sommets d'une arete de la frontiere de Neumann. (neumann est vide sur cet exemple)\n", + " \"\"\"\n", + "\n", + " h = 1 / (n - 1)\n", + " n_pts = n * n\n", + " n_elm = 2 * (n - 1) * (n - 1)\n", + " coordinates = np.zeros((n_pts, 2))\n", + " elements3 = np.zeros((n_elm, 3), dtype=int)\n", + " neumann = []\n", + " dirichlet = np.zeros((4 * n - 4, 1), dtype=int)\n", + "\n", + " # Coordonnees et connectivites :\n", + " e = -1\n", + " p = -1\n", + " x = np.zeros((n + 1, 1))\n", + " x[n, 0] = 1.0\n", + "\n", + " for l in range(n + 1):\n", + " x[l, 0] = l * h\n", + "\n", + " for j in range(n):\n", + " for i in range(n):\n", + " p = p + 1\n", + " coordinates[p, 0] = x[i, 0]\n", + " coordinates[p, 1] = x[j, 0]\n", + " if (i != n - 1) & (j != n - 1):\n", + " p1 = p\n", + " p2 = p1 + 1\n", + " p3 = p1 + n\n", + " p4 = p2 + n\n", + " e = e + 1\n", + " elements3[e, 0] = p1\n", + " elements3[e, 1] = p2\n", + " elements3[e, 2] = p3\n", + " e = e + 1\n", + " elements3[e, 0] = p4\n", + " elements3[e, 1] = p3\n", + " elements3[e, 2] = p2\n", + "\n", + " # Liste des sommets de la frontiere de Dirichlet:\n", + " p = -1\n", + " for j in range(n):\n", + " p = p + 1\n", + " dirichlet[p, 0] = j\n", + "\n", + " for j in range(n * 2 - 1, n * (n - 1), n):\n", + " p = p + 1\n", + " dirichlet[p, 0] = j\n", + "\n", + " for j in range(n * n - 1, n * n - n - 1, -1):\n", + " p = p + 1\n", + " dirichlet[p, 0] = j\n", + "\n", + " for j in range(n * n - 2 * n, n - 1, -n):\n", + " p = p + 1\n", + " dirichlet[p, 0] = j\n", + "\n", + " return coordinates, elements3, dirichlet, neumann\n" + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "metadata": {}, + "outputs": [], + "source": [ + "def show(coordinates, u) -> None:\n", + " \"\"\"Fonction d'affichage de la solution u sur le maillage defini par elements3, coordinates.\n", + "\n", + " Args:\n", + " elements3 : matrice a trois colonnes contenant les elements triangles de la discretisation, identifies par les indices de leurs trois sommets.\n", + " coordinates : matrice a deux colonnes contenant les coordonnes 2D des points de la discretisation.\n", + " u : vecteur colonne de longueur egale au nombre de lignes de coordinates contenant les valeurs de la solution a afficher aux points de la discretisation.\n", + "\n", + " Returns:\n", + " None, plots a figure\n", + " \"\"\"\n", + "\n", + " ax = plt.figure().add_subplot(projection=\"3d\")\n", + " ax.plot_trisurf(\n", + " coordinates[:, 0], coordinates[:, 1], u, linewidth=0.2, antialiased=True\n", + " )\n", + " plt.show()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "**Partie I : maillage triangulaire et conditions de Dirichlet**\n" + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "metadata": {}, + "outputs": [], + "source": [ + "def f(x, y):\n", + " # return 2 * np.pi ** 2 * np.sin(np.pi * x) * np.sin(np.pi * y)\n", + " return np.ones(x.shape[0])\n", + "\n", + "\n", + "def u_ex(x, y):\n", + " return np.sin(np.pi * x) * np.sin(np.pi * y)\n", + "\n", + "\n", + "def u_d(x, y):\n", + " return np.zeros(x.shape[0])\n", + "\n", + "\n", + "def g(x):\n", + " return np.cos(x)" + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "metadata": {}, + "outputs": [ + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAPsAAADyCAYAAACLfbNuAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAACEsElEQVR4nO29d3hkd3X//7rTm0Yzo96llbb36oZtSgBj3Kg2JAZTY0I3JIGEJ3EK5QGSfPPFlN8XnJhQjLEN2GCzAQw2GBt7d70rrXZXq97LjKb3en9/jO7dGWmatJJ21zvv5/EDq7lz752Zez7nfM55n/cRRFGkjDLKePlDcaFvoIwyylgflI29jDIuE5SNvYwyLhOUjb2MMi4TlI29jDIuE5SNvYwyLhOoirxersuVUcbaQ1iPi5Q9exllXCYoG3sZZVwmKBt7GWVcJigbexllXCYoG3sZZVwmKBt7GWVcJigbexllXCYoG3sZZVwmKBt7GWVcJigbexllXCYoG3sZZVwmKBt7GWVcJigbexllXCYoG3sZZVwmKBt7GWVcJigb+wWAKIrEYjESiQRlKe8y1gvFxCvKWGWkUilisRiRSET+m1KpRK1Wo1KpUCqVCMK6aBmUcZlBKOJZym5nlSCKIolEgkQigSAIxONx+e+iKJJKpWQjj0ajVFRUoNFoysZ/eWBdfuCyZ18HSGF7pkFLEAQBQRBQKBTysUNDQ7S3t2MwGICy5y9jdVA29jVGIpFgcnKSZDJJU1MTgiDI3jyX0UrGr1QqUSqVstcPh8Py8SqVSv6vbPxllIqysa8RMsP2VColh+/LRS7Pn0wmSSQS8jEqlUr2/AqFomz8ZeRE2djXAKlUing8LoftkjcvFYWOl84nYbHxC4KQ5fnLxl+GhLKxryIkw5OSb5I3zme8+UL55SCX8ScSCfkepKhApVKh0WjKxn8Zo2zsqwRRFInH4ySTySUGuNjYi3n75UYCi9+72PgnJycBaGhoKHv+yxhlY18FSLVzyVPnyrhfKPJM5v1ICb94PJ7l+aU9v1KpLBv/yxhlYz8PLK6dS2H7Yqzmnv18IWX6JeQyfinZp1Kpci5eZVyaKBv7CrG4dl7IIBYb7+zsLP39/ajVaiwWC1arlcrKyiVGuJrId3+5jD8WixGNRoF03kGtVsuev2z8ly7Kxr4CSEm4fGH7YkjGnkwm6evrIxaLsX//fkRRxOv1Mj8/z9DQEEqlEqvVKmfyLwQKGX9msi8z7C/j0kCZLrsMSGF7T08PXV1d6HS6kt7ndDqZnp4mEAjQ1NRES0tL1mIhIRqN4vF4GBkZQRRFDAYDVqsVq9WKyWRasUednJxEoVDQ2Ni4ovdLkJ6VVCrF6dOn2bJlixzul43/vFCmy15MyKydS8m4UuF0OnE4HBw4cACz2Zz3OK1WS11dHR6Ph/r6ejQaDR6Ph4mJCfx+f5bxGwyGdQ+nMxN9kUhEjgBisRixWAyg7PkvYpSNvQgW186lULYUY08kEpw5c4ZIJEJtbW1BQ8+EFPbr9Xr0ej0NDQ2Iokg4HMbtdjMyMkIwGMRkMsnGr9frz+tzrgSZ1F4419Sz2Pgzef1l479wKBt7AeSrnQuCUHRP7ff7OXnyJK2trRiNRqanp8/rXgRBwGAwYDAYaGpqQhRFgsEgbrebgYEBIpEIFRUVsvFrtdrzul6p97T434tr/KIoEo1GlyT8ysa//igbex4Uqp0XKo2JosjU1BTj4+Ps3LmTiooKvF7vqpfeBEHAZDJhMploaWkhlUrh9/txu92cPn2aeDxOZWUlVquVZDJ5QYyqmPGLopgV8kulvjLWBmVjX4TMsD1f7VyhUOT07IlEglOnTqFQKDh06BAqVfrrXQ9SjUKhoLKyksrKStrb20mlUni9XtxuN3Nzc3IkYLVasVgs8r2tJ3IZfyqVkoU8pqenaW5uRqPRlDv61gBlY89AqbXzXMbr8/no7e2lvb19SdY71/FrRZeVoFAo5JBeo9EAoNfrcbvdjI6OIghC3hr/emHxdzw7O0tjY2NZxWeNUDb2BRSjvGYiM0EniiITExNMTU2xe/dujEbjkuMvJF1WgkKhoKqqiqqqKgDi8Tgej2dJjd9qtWI2my/YXnpxjX9xL3/Z+FeOy97YS6W8ZkJK0MXjcXp7e9FqtRw6dCivd7yY6LIS1Go1NTU11NTUAOnymdvtltl9Wq12VWr854Ncvfxl4185LmtjT6VS2O12RFHEarWW/KAIgoDf76evr4/Ozk7q6+uLHn+hPXuxz6bRaKirq6Ourg6ASCSC2+3OW+O/ECjF+MsqPvlxWRp7ZhLO7/cjiiI2m63k9/p8PiKRCPv37y/pwb/Qnn0l59LpdDQ0NOSt8UejUaanp7Fareh0uovG82cmVyVoNBq0Wu1l39F32Rn74rBdqVTKhJliiMVinDx5kkQiQVdXV8ke7mLw7OeDXDX+F154gWQyecFq/Pnuc7Hxj4+Po1arqa2tzWrnvRx7+S8rY88lF5WvjLYYUv1648aNBAKBZT0kF9qzrzakRbKlpUWu8QcCAVwu15Iav8VikasBF+I+pVq+SqW67CW8LgtjL1Q7L0Z9FUWR4eFh5ufn2bdvH3q9nmAwuKyutJebsS+GQqHAbDZjNpuX1PglZV2pzLfeNf5UKpUlD1ZMwuvlbPwve2MvJBcFhamv0WiUkydPUlFRwcGDB4tqyuXDpWa854vMGj+kW4I9Hs8FqfFnGvti5DL+l7OKz8va2EupnecL451OJ319fWzatEkuTxV7Tz4sNvZgMEhPTw+CIGC1WrHZbFRUVKx4MbnYoVQqS6rxS7Lbq1njX875lqPicyka/8vS2JdTO19suKIoMjg4iMfjYf/+/Tl71s/Hs8/OzjI8PMyWLVtQq9V4PB6mp6fx+/3odDqsViuxWOyCJbnWA7lq/C6Xi0QiwdGjR9FoNHJkUFFRcV4GdT6LRyEhj+npaerq6jAYDJeMhNfLztiXIxcF2cYeiUTo6enBZrNx4MCBgnTZZDJZ8j1JWwWp3VXaEiQSCbm2nVnecrlc2O12XC6X7Plfzsav0WioqalhamqK/fv3yzX+yclJAoGAvAharVaMRuOyDGo1I4VM43e5XNTV1WWp+Eie/2Lt5X9ZGbsUcklffCkPhWTsDoeD/v5+tmzZIoebhd6zHM8eiUQIBoM0NjayZcsWBEHImugC2eWtZDKJWq3GZDLJGe5EIiFnuK1W6wVpZFlLZBplvhr/6OgowWAQo9GY1cdf6Hde7W2BhGQyKY/ognNcBqmX/yMf+Qif+9zn2LJly6pfe6V4WTwxUtje19eHxWKhtrZ2We/3eDxEo1EOHDhQkgddThgvLSJarZb29vaS3iM9vBUVFVRUVNDW1kYymZQz3GNjY/J+X0pyFXugL+bwEvIPzCjUxz84OEgkEpFFPHJFQGtl7KlUKivEz6TvQtrzXwhBkUK45I09s3a+3MRZOBzm1KlTCILA/v37SzaIUq4j7f29Xi8HDx7k6NGjJd9XLiiVSmw2m8z0i8fjuN1u7HY7AwMDaDQabDbbunDZ1yJ5KIpiyX0Ji/v4A4EAbrebM2fOEIvFsmr8a2Xs0r3kg6QkdDHhkjX2XKOWlEplycYuGUlnZyfT09OrSpKJRqP09PRgtVoLLiKFvFkxg5JYYVIUE4lEcLlcjI+PEwgEMBqNsvGvtnGuxtiqxcg1zroUZNb429raSKVS+Hw+ec/v9/sZGRmhqqpqVWv8xb7TUChUNvbVQL7aeSkeN5VKcfbsWcLhMAcPHswaj1QqCtXmJabd5s2bqa6uXtZ5zwc6nY7GxkYaGxvlUNflctHf34/f78doNMr1b7VavW73VSpWawFRKBRYLBYsFgsdHR0cOXKEqqoqvF4vY2NjAOtS44/H4xeMOZgPl5yxF6qdK5XKglnyUChET08P9fX1cqJsJRrtuRJ0oigyOjqK3W7PW7IrFedbZ88MdVtbWxkbGyORSBAIBJiYmEAURSwWCzabbUUP/Fp59rUKt6urq+WFN5FI4PF4cDqd59XHf7HnQHLhkjH2UmrnCoUib1PL7OwsQ0ND7Nixg8rKyqz3LNfYFxtjPB7n5MmTGAyGLKbdxQJBEDAajXIrrvTAS6QWlUqVRe4p9iCvhbGvxTlzQaVSZRl/LBbD4/EwNzdHf39/STX+Ygvxen2W5eKSMPZSa+e5DDdzCsuhQ4eWhLArMfbM93i9Xk6dOsWGDRuK9rWXikLbhNXA4gc+Go1m7XENBoO8389V2rqQCbrVhkajWZL7KFbjl8puhXAxGvxFb+zLGbW0OEEXCAQ4efKkPIUlXzJsuZCMcWJigsnJybxyVCvFej8kWq2W+vp66uvrEUWRUCiUVdqS2ldtNpu8D71YEnSrjVJq/BUVFUB+g04kEhdE068YLlpjX4lclEKhkPfs09PTjI6OsmPHjpKHMyzn3lwuF4IgFJSjKoRSQuULASnkNxqNNDc3Z0lU9/b2kkwmqaioIB6Pk0gkVjW7fTEYeyby1fjn5uYIhUIcOXIka1CHlKeRFoVieO9738svfvELHA5HryiKO3JcXwD+E7gRCAF3iaL4UsbrZuA08DNRFD9S7HoXpbEnk0kCgQAajWZZfGOJgnry5ElSqVSWnPNqIRAI0Nvbi0ajYceOJb/PquBiaoRZLFGdTCZxOBy4XC5OnDghZ/gXN/MsF2uZoFstSIlPSG99tm7dKi+E0lZRrVbz7LPPlmTsd911Fx/5yEfYv39/vkPeAGxc+O8K4JsL/yvhX4Dfl3r/F9W3K3nzSCTC8ePHl91YEIlEmJubw2q1smvXrlU39JmZGXp6eti0adNFx45aL0jZa5PJxIEDB9ixYwcGg4Hp6WmOHj1KT08PExMTBIPBZS1Ya5X0WwtIe3ZBEOT6/p49ezhw4ABWq5WBgQG6u7s5ePAg9913X97zXHfddcXk0G4F/kdM40+ARRCEBgBBEPYDdcCvSr3vi8azZ9bOpdJWqT++VCsfGxujsrKS5ubmVb23VCqVleRb65HKa+HZV9OQMu8tU6gyl1Zd5n6/EBX5Uirn5UvQKRQK2tvbueuuuwD4+te/zuzs7PlcqgmYyPj3JNAkCMIc8G/AXwB/VurJLgpjX45m+2JIU1iUSiW7d+9mcHBwVe8tsza/detWuYnlfI0xmUyu21im9WLQ5drjZspVFWrmWQvDXEtjL3ReicFoMBjYsGHDql8f+CvgSVEUJ5djKxfU2AvJRZWCxVNYIpHIqnpciVK7fft2LBaL/PeVlOsy4ff76enpAdJlMInzLnHaL6Y9ey6U6oUFQSi5mUdiQ64m1tuzS1hFquwU0JLx7+aFv10FXCsIwl8BJkAjCEJAFMXPFDrZBTP25fadL37v+Pg409PTWWWvYgy6YueU7iGVSjE4OIjP5+PgwYNLaI8rNUZRFJmZmZGrBFqtlng8nsVpN5lMshLKxYqVLkSFmnnm5+flz7xazTzr1fG2GKVm40vA48BHBEH4EenEnFcUxRngz6UDBEG4CzhQzNDhAhn7+YTthaawrNTjSsYrCILcxGKz2fI2sazE2EVR5MyZM0SjUQ4dOiRTdbVabVZdV6K0Op1O3G63TGu1WCwX1QKwGl44s5nHZDIRj8dRqVQ5m3lWkhBdS89eKPlbasfbO97xDp5++mmAzYIgTAL/CKgBRFH8FvAk6bLbIOnS23vO577X1diXWztfHC56PB5OnTqVdwrLSo1dep/H4+HMmTNFm1iWe51oNEooFMra9+eKQKSwt6amBr1eT2trKx6PB5fLxfDwcM6Q/0JgrTLnGo0mbzNPNBrN2u+X0syzlsZeKNkYDAZLGjry4IMPSv8354cR0x7lw4XOIYriA8ADRS/GOhp7Ls32QpBCcknve3R0lLm5Ofbu3Zt3OMNKH0BBEBgZGcHtdpfUxLIczy51wel0Otrb20ve64qiuESoMRqNLgn5JeNfT9mqtcqcZ3rLxc08Uuuqy+UquZmnWLi9UhRL0AWDQVpaWvK+fqGw5sa+uO+81LBdor5KU1iMRiOHDh1a9ZU6Ho/LIeOBAwdKFlAoBimvMDMzw759++jp6cm7QAzNuAnHEmxvrS547lwhf2am+2IN+UtBsQUks3UVSmvmWatqRyl79outlx3W2NgXh+3L8QZKpRKXy8XQ0BAbN25cttRUKfB6vfT29qLX6+no6Fi1ByOZTHLq1CkUCgUHDx6UCRiZxi59F4+/2M/nvv97uhosTLkCXL2lmT0tFrbU6ilUtMmV6c4X8q+2J74YauKlNPNIidXVvt9i2fjLztglr/ziiy/KCalSITVjDA8Py1NYVhOieG6m+t69e+nv71+1kl0oFKK7u5uWlpYscs9iY08kU3zhkef4n9/1AhCMJnAHIjxxdJAnFhSsNjb0cs22Fl6xtZlDGxvRafL/XIVCfrfbjV6vJ5lMrkrIfzG2uOZq5hkfH8fr9XLkyJGczTwrxTqW3lYVq27si2vny62fSlNYALZt27bqhp5JwpGy+ctVi80HSVxycc88ZBv7vC/Eh7/5JC/0T8uvj8x52NtRx/GROflvAzNuBmbcPPBUDxqVkv2d9bx+bwf7OhvY0lxYATcz5B8ZGZEHWJ45c4Z4PH5eIf/FaOyZkJp5pLbU5ubmJc08mfv95dKqSyXVXGxYVWPPJxdVKjKnsMzPz6/Y2+Z7cAKBAD09PbS1tdHU1CT//XxJMqIoMjQ0hNvtzlmXh3PG3j0yx91ff4IZd2DJMbOeICqFQCK1dOGJJZIEIzG++rMX2FBnYcYd5OqtzbxiazOv2NpClTn/oigIAjqdjtraWlpbW2VyS2bIb7VaqaqqKinLvxaEn7Vk0OVq5pHGUY2MjCy7maeUPbvUBnsxYdWMXRRFotHoimrnkhJr5hQWl8t13jXzTEgtrzt37lzyQ5yPscfjcXp6eqioqCg6WOJnL5zlH374LLFEbuLPjDvAwa4GjgzOLHnt0Mb030UR3MEoTn+Ix17o57EX+hEE2NxUxbXbWrhmazMHOhvQqPM/jIvJLdJ+d2JiQtark17PV5m4mD27hHwLyOItTywWw+12MzMzw9mzZ2XBCpvNhsFgWHJfpYTxL2tjlww8n4xPvh8y3xSWlbLhpCy+9CNLk1ji8XjelteVGrvf7+fkyZN0dnZSV1eX97hYIsn/98wgvzk9R1eDDYtJh8sXZnDGxWIfOTjrRq9WEo6nP7teo2RzUxUvDpxbACbmfRzoauDowqIgitA36aRv0sm3f3UCvUbFwY0NvHpXO4c2NhT9kRfvd6X6dl9fnzx+WSK3KJXKiyJBV+o5SwnRM5t5IC0x7nK58jbzFLvXVWTQrSpWNYzPZTTS33KthIWmsKzU2CUBC5VKJTexNDQ00Nrauiw5q2KQdOd27dpVMBlj9wT50Def5Nhguvvp7JRTfs1s0LChzopSITDm8DLvC+MORNjXUcNLIw7qLUZ0GhUnRuxLzjs868agURGKJZa8Fo4lGJ3z8v8dPs69D/6B6godV3TV8Wf7NnL1lmaspvw8gsX17cyQf3R0FIVCIXv71TT6i2kB0ev1NDU15W3miUQizM/P553MU4xhd6Gw5nckGW2msadSKQYGBvD7/XmnsJyPsadSqbxNLPneU+o+VGp3jcfjXHPNNQV/1GODM3zoG09i9wZzvu4LxTiRkZBrqzFTazESicbYXGdk1p9g1pP7va5AhIMbGzgysDTk39VWw9Csh2A0zW2Y90d44vgYTxwfQyEIbGut5tqtLbxiWzN7NtShLhCSLg75Y7EYY2NjuFwuXnzxxZJC/lJwsXa95Spxvvjii/h8PsbHxwGyJvOUumBJKjW1tbX09vbmum5OlRpBEPaQFrEwA0ng86IoPlTKNdfc2FUqFYlEQk5ahcNhenp6qKmpKThAQcoeLxcKhYKhoSGi0WjeZNliCCUKPEYiEbq7u6mtrUWv1xc09O/9rod/fvD3xJO5z1uh11BtNmA2aNColCBCOJ5gzhOktkLP4EyEDQ1WNmltzLqDTMz7lpyjd8xBVYUepz8s/+3gxgaODc6SyrN4pUSR3jEHvWMOvnn4JYw6NddsaeaVO9s4tLGe1hpLwe9AUl9VKpV0dHTkDfmXO5DhYvLshaBUKlGpVHR2dgLnxk/b7Xb6+/v54he/SDKZpKenh507d+a9vqRS8653vSvfpfKp1ISAd4miOCAIQiNwTBCE/xVF0VPs3lfV2HP9WJkeem5ujsHBQbZt24bVai14LqVSSTQaXdb1JWXQuro69u3btyw5q2LGLtFepS3H7Oxszgc0Gk/wue89zRNH+mmqrqDSoEWrUhEKBRGUagLRBA5vCH84hj8cy3qvRqVgR1stR4fSIf+p8Xn5tdpKAy3VZpKpFEMzbvyRuMy6c/rD6fe21uT09AJgq9BSU2nEoFWjUihIpESCkTiuQJjfnxrH4Qvx6HN9zHmCXLO1mVdsa+HqLU1U6JdGXZlJ2GIhv+T1i0lUXyrGDtnViMXjp7/yla/wjne8gy9/+csEAgEee+yxnOe47rrrGB0dLXQZWaUG+JMgCBZBEBpEUezPuI9pQRDsQA3gKXbf6xLGS/VdaQpLKd42UzyyFEhlO7PZTGNj47IenELGLooiY2NjzM3NZfHmpWggc3sy7fLzHz/7E9MuP83VZgamlybg8qHarMdi1PHSUG5lE7s3hN0bAkCpENjUaMNi1OINRtjebEOpUuILx9i7oQ61UkFKFAlFE3iCEeyeIE5/FKd/6eJpM+loqTZzfHgOhSDQXlvJQ8+e4aFnz6BUCOxqr+Warc1cu62FXe21KAuoCOUK+V0ul8xqKxTyX6xhfK5zFnq2WltbMZvN/OAHPzjfS+VUqQHk1VwQhEOABhgq5YTrwo0/ffo0zc3N8hSWUlDqnl0URYaHh3E6nezfv5+RkZFl7/XzGbtEwFGpVEuGPyxmxD3fN8lHvvVLvMEotRYD064AlQYtnfVWFAqBwal5POHc25Jmqw5/JMHgjDvvPVbo1AthvzZdVhPTiTizQYtSqSASS1Jp0DI+78XhDec9TybaaiuJxOIMLFw3JYqoVec+YzIlcnx4juPDc9z3xDHMeg1Xbm5iV3Mlu5vMtBc5v0ajKSnLb7FY1syzr3aPwDr2sheEkNai+x7wblEUS8our2kYPzMzg8PhoL29veRxxRJKMXapSUYSP1zucEcJuSbJBINBenp6ltBeM98jXef+Xx/niw//kcTC/ryxqoJpVwBvKMpLw+c8dWt1BWatgmA0wYQrRCIlsr+znp4xOwpBoMFiwKAWkOxNpdYSiieZ94XTYX/Em3UPezfUcXLMQSSeZGtzFWcm05n+1hozdRYj4Wic/mkXscTS72PnQhIvFM3+3GenXEuYfBJ84Ri/OjHCyIyZHz+fRPvzkxza2MgrtrZw5eYmjLr8bafFQv5QKMTk5CTV1dUlTaUpBWvRCLOOvPh8KjUIaQnpJ4C/XxCiLAlr4tmlKSzxeJyWlpYVcbGLGa3UxLK4SeZ8J7zAOTmqXLRXCYIgEIrG+Jvv/o6f/els1msnhuZospmYcmWz5Mbn/fL/r9BrOLSxnkRSpKPWQv+0ixlPaNFVcmfiIU2yyay9RzOIOuMOH+OOdEJPq1aysc6MUa/BFYgxPu8rmsSbdPnRqZVE4ksX2/2d9XSPzJFIiRzaWMkPnjnFD545hVqpYM+GOl6xtZlrtrawo7UGhaLwQI/MkP/IkSPodLqSQv5SsRZhfCnGvpYqNYIgaICfkt7PP7KcE666sS+ewjIxMXFe9fLFyGwdzdXbvhLPLoXkEpNPmqleKLfgCMT43Fceoy+jbi4hkUpRb6tYYuwSrEYd9VYjT/WM01lvYWjWQ5OtgkabiXA8wcCUk2gObwygVQq0VZuyDB1geNbDvg11vDSc7ZGj8SQDc2nDVykErtrcRCyRZFd7LcNzbnyh7CQhgMMb4tDGRl4cmM76++JSX8/oHHWVBua8IeLJFEcGZjgyMMN/PH4Ei1HH1VuaeMVCI0+9tbC3EwSB+vr6JcIVZ8+eJRqNZnH5S83yr5WxrwahRlKpmZ+fZ5kqNW8HrgOqFiSpIF2WO1Hsmqtq7C6Xi97e3qwpLCqVilhs6QNVDLnC+EQiQW9vL2q1Wm4dXYzlJvak98TjcV566SUqKioKlgQBnu+b4PM/P8WkJ//e+MTwLM1VFUw6/Vl/31BvIRiJyyG3Tp3+CaZcfqZc6WM1KgUbqo1UWczMes6V3eotRrQaFf1z2eG8hClXALVSkbPcV6FT01xt5vmzU+c+tyCwscGK1aTDE4wwOONGouV3ZxiyRqVge8vSTH8knmST1cScd3FEAp5ghCePDfHLY0Mc3NiA0x/h2m1pr39oUwN6TXbIn7lnzyVcIbXvLifLv1YJutUI4zNUaiAdomchn0qNKIrfB75fyr0uxqoae2Vl5RJKqlKpJJFYyvIqhsXGLlFTJSXZfFhJGB+JRJienmb79u0Faa8A33jyKF/9yfOkRJG2GjP11gr84Sj9U86sBpZkSqRCk/2g7WytYmDGSyR+7vs4NTGftd8GiCVSDM8HGZ5Ph/E1lQa2NVcRjiU4PTFPPsx5gmxvNHNqOrsmX1epR0DIugakE3IDGUnBCr2GznoLSqWCCYePBpuJaCJJbaUx5x4eoGfUzvaWak7luC+NSsG2lmo5ChmadfPAb0/KHXyvWDD+rQsdfIUYjqVk+Rdr1V2IMP5i7XiDVTZ2SUIqE+fDcZfeNzU1xdjYWFFqKhQe25wL0rmrq6sLGnowEuNv/vs3PHH0nC79mMPH2MLe2KhTs73BhkalZMzuwe4NcXbWS2u1mYl5HztbbPSMLw35Ib0wFEJrjZk/9k2xu72WSCwhl92cvjBDc56sY8dcYUxaFYFoekHprK1gxhMiFCv+GyRTKdzBKJUGDW21lRi1arY2q/GHo2hVyqy8QCa8oeiSbr1Ko5a6SmNOqm8skeT5s1M8f3aKr/z0BarNejbX6Lk1VckrtjZTbc4tOyZhcZY/FApladVJIf9aZPgvVeEKWCcG3UqMXapjS/3Hpc5tUyqVRCKRosdlTnnZtm1bwckdo3MePnjfL+ifduU9JhiJ053h/ZptJlpqKtFpVNjMek4M5/aMAP3TLna01dA75sj6u0KA/Z3nuuDG530oFYqs+7AadbTXVS7cpxd3MCJ3zu3vrOfEyJy8mCgEqK7QYaswYNKn5+jFEykCkRgObwhvKMqYPb1F2NpcRd+kUyb+aFRKtrVUY9KpmXZ6mXSdC90nnf6s/XyjzSR/rkJQKQRqLUYsRh3zgTAP//EMDz97hgarkTqriau3NHOwq3AHn9S7bjQaaWlpIZVKyVn+UCjESy+9JLfvrkaWv2zsC8hHtFhJGB8OhwmFQrS0tOQdt5wLpYTxEu21rq5OHs6X7z2/7R7hE9/5Fb5QbjZfVYWearMBk05NKpkgFI4gqDR4QlGm3QFSKRG9RsXmhkr0GjV2f4TpHIm7UCQ7GjFoFHTU2bLaXdOJs+wsvDsYwT2cXtwEATbUWVAoBF6zu41AOM6ejjrcviC+SBx3MIrdF8HuK7wY7u+sp3vULpcSIe2NM7cQNqOWjnoryVSK4TkPp8bnqarQYTPpcfhCeILp76vKpMdm1lGh06BSKkimRILROO5AGIc3zLQrIH8fVaEUZr2GI4Mz7O2o49u/OoFOne7ge8W2dN/+xsbCqq1Sb7rVasXlcrFjxw7cbjdTU1P4fD559rzNZluRMEopCbq1kFBbDaxbI8xyINFqdTodra2ty3pvMWN3uVycOXOGrVu3ynvAXI0woijytV8c4f5fH6fGrKez3oJmYZsSjidwuP24gzGc/nAWN13C7o5aBmfcBCNx9m6o4/jwOS/XXG2m0WYiGIkxsFAHH57zsKejlhMjdppsRuLxRM59cN+kkwqdGn9k6VZFFNP7doNWzQv905h0ajbUW9FrlPgi8aLbBYD9G2o5lofFl/U9BqO4hmYxaFTUWY00V5nRaxTEEikMWjVmfQS7L4QzEMYZKE7yqTFpUCgVjCxEFmMOL2a9Bl84xh9OT/CH0xPA89RZjHJ575qthTv4YOksusUhf2b7bimRY7G22Yu1vRUusjA+lUrR399PMBjk0KFDHDlyZNnXy7e45KO9wtJGGH84yud//AdGZj2011TSP+1iaNZT8j1cubmJF/qnkNaPEyNzNNuMTLrSCbfJeR+TCxl2nVrJjtYa9FoVsXiSnW3VjMx5CeQwZkgTWxZ7dwn1FiNatZLe8fR2IBCJ0zN6bs/cUm2m3mIkFI0zMJNNttGoFGysr+TYcPYeWyFAbaURq0mPQatCqRAIhMKEYkm8oTjuYISROS/VFQb+eGaWTY1WJp0BOust1FoMTDh8ObP1mdjUaGNq3kswcK5q4wpE2N9Zv2ThmfMEefT5szz6/FkUgsAVmxvZ11HPVVua2NtZuIOvUMg/NjaWpVhjNptzRpOlaMZf1mF8KcYuiVhUV1ezefPmFe+tcnl2qWSn0WiW0F4Xv2dw2sUHv/4Ew7NudrXX0jNqR6NSsr21BpNOzZQrIBvqYhg0KjY3V/OnjPIWpD1uRR52WSSelI3z4MYGpp1hNjfa8PgDTHkiOYktPaN2qs165n3nPOamRhsObyhvSyykRS8mFi8yGhVOf4g6q4lwNMHBjQ2IKRF/OIrLH8YViDLrCeY9rwAcyFDXicSThKNxujMWmeaqChqsaQ7B4LQr6zPtaq/h7JSLaI7P+dLwLJsabXn3/rWVBsbsXp7vm+LrvzyGUavm0KZGubZfDJkhP5xTrJmenubs2bPo9folIf+lKjYJ6+DZS6l7S00smaG1hOVmVBcbu0R7bW1tzdKdy/Wewy8N8en7fyV7ValEFkskOTV+LnnWYDNh1SoIx+LM+GJE4kkabSY0KiXHh3OHwGemPbRWmxifX7pfVykFdrfXyQmuc/V2JdtbqjHo1FltrpF4ku2tlbKx791Qx6nx+bxyVxqlQLVZj81sRL9Q148mkvhCUSLxOAgCdk8Qq0lHNJ7m6IeixfMsUmktM68wavdysKueI4PnvodJp1/mG2Qm+tRKBX/qn867vRDFdBdhLl2+ttpKAuEYzgyuQzAa53cnx/jdyTEgvTV4VX+IV2xt5uotzZgNhZmcpYT8kUgkL6sSLrMwPp8+ei5ITSwulyuniEXmVJhSkcmgk2ivO3fulEk+ee6aH/xplJ+deI7MrXv/lIutLdWcWbR3nnEF5NYjrVrJK3e0Ek+KjNk9Be/NqF3q3a1GHbUWQ859ciyRzNq311mMtFSbiSUS9E85aa5Ks+6ODMxQXanHZtJj0mlQKgQSyRSBSBynP4zLH2baHWLanR1Ob22uYsrpxxeOcaCrnqMLBqpSKtjSXIVZr8HhDcn76EyYtCqsBnXO0tqZSRdWow53cGkiUEr0HdrYgCsQwWrS0Vpjxu8PYA8m8AazE6FjDt8S5t6WZhuT8wECkcJkLUcgxo+fPcOPnz2DQaPiFdtaONjVwK6OWrmDLx/yhfz9/f0MDQ0xPj6eU6SylDB+pcIVC6+9G/jcwqH/KoridwteLAMXTDsnFovR09OD2Wxm//79eYUBl2vsCoWCRCLBwMBASbRXbzDCR/+/w/z+1HTO1xVFoordHXX84fSE7J2abBU0VVcQjqUNMjM8PTPlpqvBKne3ddRWEorFOTtVuEQF6W2AQasmnkiiU6vY1lqNSachGk+yYYFyW2q3GyzNuB8fnqO12sz4vI9EMkVfBgGnqkJPW00lKVFkeM6DTpme0zfhzr0XD0Ri7GmvyWnsSkFgz4Y6Oeewd0MdLw2ly5IKQaCz3kJVhR5vKMbgtIukKHJy1E6D1cSMO1Aw7K+q0FFVYcCo06AUwBcIgFKTXvACYY4NztI/7eJXJ0bom3Jy5eamtDrvthZaqgs5g3Mhv8lkor29HbVaLYf8fr8fvV7P0NBQSWKTKxWuEATBRppWewAQSQtXPC6KYv52yQxcEGOXBjQWm/Sykky+1E1VWVlZlPZ6dmqeu+97IqfnknBq3MHGRhsDi/aNaqXA1iYbL/ZnLxKZtFedWsXO9lr0ahVjdjdz3rBME93TUUff5Ly8f1UpBeoq0zVnnVpJKBxCo9PjD8fOiV1EPED6obYYdbzQP0NXg5WhWQ81lQZaq80kFsQt8iX4YGkTDaSJPfnC3MyKw+ZGG4lEHLNBS4VRy1AGxTYT3aMOmi06Jj3nDF6vUdJZb8uKYsbsXkw6NYFInJQoMjTrkZOhJp2aznorapUCpUJBU5UJpy/MpkabTDOOxBN4glEc3iBOfwSnP3dZsd5iRKEQGLV7GXf46Gqw8usTI/z6xAiQJi69YkGq68rNTZh0uR2EtGfPFfL/4Q9/4OzZs7zhDW/g2muv5atf/WrO8t5KhSuAVwK/FkXRBSAIwq+BG4AH854pA2sexkuQ/jY2Nsbs7GzBAY0SltvU4vP56OnpQaPRsHHjxoLHPv5CP3/7wG9Qq5Ts66zH7w8wH0zk9EaaRfmYmko9BrUyLyNOQiSe4GRGoqqzzoytQsdrdrUxM++io9pIMBrHH0ngDceZcgUWNc/4l5xzQ50FfzgmG4RamY6IHN4QjoWst0ohsLnJRqVBmxWGSxz3XJl8gN5xB9tbq7MUcjKxu72Ws1POhQUqfW9mvYYN9RaUCoUsmglpt6PR6hCIIAImrZJKvVpORkpIZ93rODa0lHQUiJxL9B3a2MCsO0iD1UQknqB33E64BFYgQGu1mVAszvzC/j4t7BFHo1LIFYlxh48fOk7xw9+fQqVQsLujlldsa+H67a3saKuRz5UrQSeF/H/1V3/FQw89xLPPPstLL710Ppp8+YQr8v29JKyLZ1coFMRiMc6cOYNGoyl5QONymlqmpqYYHx9nz5499PT05D0umUrxxYf/yHd+dRxIC0BI6jAKQWBjow2bSY/TH5LD7dMTLlprzIw7fGxuqmLeF2IsR6LNoFFRazFSadShXWB9RWJptZg5d5BJZwCjTosvEGTWF2Vjo41ao5HovJdUqDjFd3d7LWenXUQyFGXPTDqXcNMTKTFra2AxaOisq8Ro0BVk8gEEwjEUAku89cGuBo4NLW2L9YVjWfv2thoztZVGApEYAzMu9nfVM+UKIEDeLsDjw3YaLTqmPbm98sGuc5GIlOhTqxRsba7CtJBXGM0TnXU1WLF7Q0tIUZNOf94SZiKV4tjQLCfHHGxusi0x9kLPbiqVQq/X84pXvCLvMRcK6xbGHz16lA0bNtDQ0FDye0oJ4yVd+EQiwcGDBwvu713+MB/51i95rm8y97lEMStcN+vVtNiMaLRatGoV9RYj0+4AdRYjjVY9SkFAoVTjX6CbeoKRhYdu6YNXXaFFq1LQs0CJ3dxUxckMemxLdQX11jTRpn/KtST7vH9DLS+N2MkRNC0RoFgMnVrJpCvA3IgjS9IqJ7fe4eNAZ72sg7e4tFYMYw4foVgcm0nPvg0NVBo0aJRKRhz5t0opUYSUiABZMl4qhcDO9tqc144nUlmNPVUVetpqKxFFkZE5D55glA01RqZdgbzfz0tDs2yoszC86DsA0GtUfOPu13PN1uzRy4Waa1ZxUk4+4Yop0qF85t+fLvWkaxLGZ2JqaopAIMCePXvkiZulopixZ9Je29raCu7Pe8fs/OV9T8j76XywGLVUVegRknEqjHpUajWRWBK9RsAbidNkMzPnCeT1JLmwqcnGrMvPfIYGnFKZfa8T834mFsQt9BoVrWYNNouZaacfm1G7hOySiRG7N2cvO6Qz7uMOL8GFUloyJWZz6006OuosiCmRoTkPvlCUoVkPJp2aWCLJ1ubqJcamUymoNhuwmfVo1SpE8VwEY/eGcHjTVNgdrdWcHLMTXohEGm0mGm0VRGIJhmZdWWH4tC/K7rYqusfSBqxVCrTVpLXxSoEoigTCUYw6TXpPrxLw+EO01mgZnHbnHKmVSKU1DBaX9kw6Df/vw2/gQNdSx1Rs2tFKxp7lQD7hiv8FviAIgqTW+jrgs6WedM08ezKZ5MyZM6RSKaqrq1ddraZQbX4xHvnjGf7+e78lGk+iVSvTiTCTDp1ahSAIROIJ5pwefNEUnmBU5nXjCC4ITZiWRAPVFToaLAYUSjVnp+Zzkl8gnfXuGZ0jnsx+2E6Pz7OlyUZfjkx8OJZgeD6BM5ykuaoCkTThJhCO0T/tylmXnnYHUCkFEhnXycVxXwx3III7kLGNWehv12lUJFMpwrEkhzY2EE+mCIRj6UEWwQiTrgCTecJy+doLijbyPWbw4KV6u1GnZtYVYMLpZ8QRwGLUkhLT32//7LkFVaUQqDHrqTIb0GnUCKS5AtIC4wpEcAUi8rX/eHZW/p6MWjXbGqxoVQomnYGsOXsjdm9WOG8xarn/o29kZ9vy+e2SAEoxZApXNDc3MzU19T5KEK4QRdElCMK/ABK19J+lZF0pEIrc3LLjkkQiISfKmpubaW5u5syZMzQ2NhYd1rAYIyMjaLXarP51URQZHR3Fbreze/funEmQ5557jquvvpp4Isl3fnWc589OYvcEcXhDuALhnKFwLnTWWwhE4swVYKVB+sHd1GTDoFUzOe9jemGPemARuWQxtjRV5VS6Aait0KBUqZlxB7Nq4Eatmq4GKyqVgjG7h/mMppbtTZWcmkobSL79qASDVkVNpZFKgwatSoVI2ju7gxGUCgEBAbVKQYU+nZWWOupKQbFr54LVoGFzSzVGjYZYMkk0nkg3zUTizPtDOP2l/W4HN6bHYhU6ttFmoslWIROI4qkUDRYToVic//7YTWxuyj8hVxpBnguRSISbbrqJF154ofiNZmN1+3DzYNU9u8PhoK+vj+3bt8tMo9XoaYfitNes+/CmxyLPeoKMO7zUW4101ltpjFVwdmo+pwhjJvZ11nN6fD5LaCIfYokk084A1ZV6GqwmNjVVoRCTDE87805lBeibci4RrgDY3lrN8IyLcDxNGjkxYqfeYmTWEyS4iIratiAu6Q1FmfEGMWmVNFTqeHFghiqTliqzARUptBoNCAr8kRjzC11pYzm2IhsbrHLX2p6OWjmMFgTkGvicy8uEK7QkiScA+7vyG7pcXjTp0GvUCEJaNssTCGP3hph2BojEE9gq9KRSIhajjngyidMXLsnrHOiqz6mbvxiZEUaFTs3ujgZqzQY+fNN+2mstJVwpNy5m9hysgbEbjUYOHjyIWn2OLXY+ajWSEIU0brmYUg3AgD3APT95kFl32iPvaqulZ8wu/1uvUbG7vQaNWsmo3Z1FRhGAQ5ubeGERv12nVlK3kGnXqVXE43GCkSjRVJpq6gqkiRtNNhMzLj+RaAxfJIFGrWJ7ow1BTDHm8OFexBBbbDAHuho4PjybFaonkikabaac/PRMAY0Gq5H9nQ2EYwl8UZE5bwhnoPRBG9mltfQis7nRxtlpF6JIVg28Qqems8GGQiEwZvcSCMfY0lxF/7SLznoLFfr0pJtkKl3mci20tC4tL6bRVmXAHYriC0UxaNVMOHwkF9yzxailo84CCIzZPXK4nomdzRY5+lmMKpOeKrMOk06LUimQTIlpqu0Cf2Da5ecLd76SpqrCZJhiyjcXs0oNrJGxLzbs8/HskUhEbnktTnuFHz7Ty788cTZbImpRTBeOJegePZf4aa2uwKgSEZRaaqxGgpE4V2xqWqCcxnD4Qrj84bRROXI3wQBsb61hbM4tq8RIk1YzRS066izUVhrwhKIMTDk5O+VkW0sVZydd7O2slyezLsaJEbvMcMuFpqoKRFHkj31TaFVK3MEI7bWV1FYasLv9jDsDOckvEvKV1hJ5cib+SFyeU2c2aNnRVoNapaDRZmRgyr3kOy+ErU1Whue8sshmml9/rgLgCUazIowNdRaqzHp8oSgjsx52ddQy6wqwvaUanSbtWAKhCP5wHG8kUbDNdmODlf/++E3UVhY30ktZuALWIRsP5zek0eFw4PF4OHToUFa0sBixRJJ//MHTPPj7U0teOzXuoKPOwkiOEgukJZ7rzFpUaoGp4QBdDTYSyRQDM668ohWLcbCrgZcWeeRplx+lQsj628icR76PCr2GDXWVVOq16Dtq8xo6pMtTVpMup7FnctwBti80p4zavXLVwKhVsbHRhlKhYGQu2zvuaLLkLa0NzXrYt6E+S/8+E7WVBrRqVRYrTmK+qZSKokMr9m2op3t0bknS8dTE/JLOPkg3xwzPeRie86RJQm01xBMp6m2mktppM7GjtYb7P/rGoj3xEi6WARErxbrU2VUq1bLntsViMQYHBxEEoejctll3gA9948m8HWeQrsPmM/bdHbWcnZwn4kvf47Gh9IOvVAhsaa6i0qhLl9tyqLqqFAK7O+pyGsu0K8D+rgaO5TEkfziG3RfG6Y8sTI+xYDPpcfmX1r8BukftsvS0hH2ddfSMOrIy7idG5+Q9voRgNCGTXyTvWFNpQKtS8OyZ3LwDCRPzvpw68u21lQQisSVDJzOZbwBtNekIIxCJMzDjlCsGhZJpoWicTY3WJcYuwahV05KjNCe103r8QUbn/UuqIBI21Rr5m9d2EPTMo1XYirI5oTSVmsvKs+fCcj27NACisbGRUChU0NBf7J/iG08eJZanFVJC98hcTk9xxeYmXswQmshEMiVmNYTUWYwLI5MSnJ2cT4etVlNBZZc5TzAnIw3SHlmaHKPTqNJ74oX5fBU6FV2NVYgpkcHZc1z3TAnmfFnveCL/Hh/S3tEZCKNUCAzMuDHrNXTWWxGEtP68Z1E04/CFssJqgE2NVqZdgYIcfAljDi9jDi8GjYrOOisNNhOCkC4/For2T4zY2dFaTe8i+m6lUUt1hT7rt5GQ3U6rYHtrFQaNmhlPgMkFHsPVW5r45oduQEzEcblcDAwMyC2sVVVVebXpL2VlWbgIw/jJyUkmJibYu3cvyWSSkZGRvMc+8FQ3n3/oD7JOeoVeQ1eDjXAoyIwvhjfjoY0nU9QYlEhOKC00UbUkEZcLygVhRJtJjyimE3z7O+txenwIYkLWV8/5eeZ97Ots4KWhbKM80NVA98icfO+DM252ttXIrDp/JCF7LYnrXqHXYPeG2N5ajUapLFjeOj4yR0dtZc4mH0kQUpKR9oVjslT0uay7AVcgzNCMG5H0QIhKnQpvJMGmOhMjc96c+vRS55lJr0alUBJPJrPq84Mzbkx6jbxAtlRXUGXUEohEGZsPLjmnKxDJ0sJPbxuUJSkHxRIpTo3PYzXqqDbruWZLExsbbXz6tivTIpYaNQaDgebm5ixt+pGREVQqFTabjaqqKoxGI4IglPfsJV2kBGmqVCrF6dOnZSVZpVJJKBTKSaqJxBP83f/8lp8815f1d384JofyCkFgc1MVFqOWyXkvU64gk94YBo0Ki0m3IDRxLgQ0aZXUVpqoNOnRqBSkUiKhWAKXP4zdG0z3sC9kkXe319I/5SScEda211ZSZzHhC0U4O+XM8uTzvmAWFXRXqy3n/tyfJz+QyXWvNGixGLUoFAI7WqsZnHHnJPSIIhhzdG5takxzxT3B3NdanHWvNGjZUG8BQK2AaDzOtNNPs0WLUadFpVaRSKUXjPS02PydZ9oFPkJmJJTJHNRpVGxprkKrVjK18H1PuwJyP3uzzUQ8lZKPl6AUpMVYh16rQkDAFwgSTog4/ZG0KGcwwqamKv7mzVehUi4NxRdr00ejUXkoRTAYpKKiouj04YtZpQbWgFQjiuKSCTB+v5+RkRF27dqV8z3hcJienh7q6+tpbW2Vo4NIJMKpU6fYv3+/fOzkvI+7v/HEEtnlYmi0meist6HXKAlFE8STKWKJJN5QRG4hLQUHNzZwdGCm4BdjNmjparAiAIOzbrzBKHs66uifctJZb+Fknq4yIMu7L0ZzVQUpUcxSp9WoFHQ1WDFo0pJZmewwIKuOv7i0Vuj+aysNmHQaNCqFXD6zmnSkRBGfL0AkCSN2X8kPiFmvocFmKql3X0JTVQWNVhOxRBKFkI7OlAoFWrWSlCgSjiZwBSI4vMG82zcJb716C//659cXnD+XD6Io4vP5mJiYwOv1otPpsuSppX381772NRobG7nrrrsKnu/w4cN8/OMfJ5lM8v73v5/PfOYzWTclCEIb8F+k5667gL8QRXFy4bUvA28EFMCvgY+LJZLyL/ieXaK9btu2TdYCy/e+585M8JFv/TJnnRXSe7naSiNKMYleqyYUiZBESSCSYM4TpKW6kpE5N65AmI2NVWiUCrzBaEmGrlEp2NZcXRJpwxeKyp10Bq2KKzY1YjZoUKtqOTKQWyRDgidPiWhrcxVTLv+S2WyxRIrTE+f2rk02Ew3WNBusf9olJ+72tFfTM+YgJYqyJ7SatBi0GgQB4okk3lAMhzeIb6HeLUGlFNjZWstzfVNp1t9Mei9kNeroqKtMRwN55sZBeva8UaspaujpeviC2o5SIJlM4Q1FMWrVaSNNimhUSqad/oJU3cV496t28ndvu3rFnHVBEKisrCQcDmM2m2loaJCFK3w+H0ajkWAwiNPpLNpanUwm+fCHP8yvf/1rmpubOXjwIJ/97Ge3iaJ4OuOwr5LuZ/+uIAivBr4I3CkIwtXANYDkNZ8FrqfEZph12bPnCuMl2qvD4Vii9ioh09j/3+GX+LefPU91hZ4dbTXoNWoUgkA0nsAbiqaHHASjS2SNIG2oezfU80J/en++p6NOrhGn96hWDCqRQEzMucetNuup0Kk5MZq7GaWqQk9NZVohRaVQLNTnowv03Aguf4SBGTcuf5jaSgNtNZVEE0n6c3jZCWeAjbUmBuznHmZprHO8COsPkEkrBo2KDXUW6i1GaoxK4ql0S6jkCWfcS6OAXDBq1bTWmOU9fd+UU76/TM36c7x6PfO+IMMLlYtqo5pUMsmYw4tSIVBXacBWoU/z2xcYdE5vEGcgkrMevqOthjMLIh+7MzrgGqwmmqrSlNeBRSKWmbj7hr3cc+sVRT9nKZD27Gq1mtraWmpra+UhlD//+c/55S9/yc9//nNefPFFPvGJT+TUPHzxxRfp6upiw4YNANxxxx309PTcCmQa+zbgnoX//zvgZwv/XwR0gIY0/0sNlNYpxBp59sUCFosZdIlEgpMnT6LT6eS56rmgUCiIxBP864/+QN/UPK01lQxMu/L2RedCbaUBs0GbtUcOZ/SDp/eo51R9qir0dNRbiCdS9E85aaqqwBcME4xE2dFag2FBJTYUjTPv9uMJJ/JqxwPs3VDP6QmHLKNk94awLyTz1EqBzY0WKo0Gxh1eOXse51wSaEeTJW+2X2oMMWrVqJQK4olUVrvtmN2HTq1ixunDGYzRUm2mucqEWa9hYMZVkGST/i50VOi1S+i87lBsSeXj3Ny49HdpM+nY0VaDSiHg9IVRKvw4AxGm3UGm3YV7Dc59d+n581J0Mu1OL2KhWCJrscpqqskQ5vzYG/fxkZty89hXgmQyuYTrIQ2hfMc73sHRo0d561vfSiKRyLu/n5qaoqXlXPdqc3MzLBWg6AbeTFqH7k1AhSAIVaIoPi8Iwu+AGdLGfp8oimdKvf91E6+QEm3Lob2O2T38w+N9TLjDMnusxmygvc5CLJGgb9KZU4tMwuamKhy+cyIUEs5OOfNKFGca7sGuBtw+P7UVOuaD8SUqK8VwxaZGXujPH7bHkyJnpz2wUG5rrTFTbzHhC0fZ31mPQhCY8wTZ2VaDWimQiCcIhCMEo0nc4TgOXxhHnjq0xail2qyX++elZhqJWisRepQKBUNzniURUXN1BfFEMmcr73wgllPTPRMNVhNHB2cJReMoFQKd9VY66tPJwVLagw901vPS8FwWo8/hXVoChKXTauotRt60p54PvHZ30essB8lksqD6TDAYpLq6mn379p3vpT4N3Lcwkvn3pPvYk4IgdAFbOTf19deCIFwriuIfSjnpuhi7FNrPzs4yPDzMzp07i4ryPX1ylI//v/+Vy2f1NhPj8z4cvhAOX9oz6tQqdnfUoVIIDE7P4w2f89i726o4PeXOG/rmylRnItdD1VZbSb3FhCcYYWDamdczpuu7tQUNXYLFqKPGrKdCr0GlUpJMpdBpVGhUSqLxBNWVBoZn3Xmz57nQZDOREkUGZzzy30btvqzhjP5wjO7R9EIgCT1KCj0qhYLZHOoumTg75aTSqM25bdrTUcepCYf83S/uoa+q0NNeW0kylWIwh17e3vYaWTxjMY4PzxakDSsVAh+/+SAb9OELMq652HPd1NTExMQ5ZanJyUlIG7MMURSnSXt2BEEwAW8RRdEjCMIHgD+JohhYeO2XwFXAhTP2xWF8KpUiEokwNTW1pElmMURR5OtPHOXff/anrFX91LgDo1YlizBAugTXnbH37mqwUmU2QDLOC4P5xR4gXTeut5qYXbRvNWpVdNTlpo+O2b1yp1ilQUtdhRqj0cTgjEtO8qVnm+nlEqBKIVCzUKMXUglSySQqrS7NnvME8QQjeDJaRzMz7vs663lpaFaWy7KadMx7QzmVVSS0WPV4QhH8kez+hHlfKG+kIQs94mFXWw1z3hBd9RaSKZGh2dzilYFIXOb+ZyIfxz4TmdGTSiHQUWOi0qjDE4phNWg4Ppo/gkqkRAw5JLkhrcf3lfe8mhv3d3Hs2LELMq65WOnt4MGDDAwMMDIyQlNTEz/60Y8gLVYhQxCEasAlimKKtDjFfy28NA58QBCEL5IO468H/k+p97/mnj0Wi9Hd3V0S7TUQjvGp//o1//vS0JLXgpE4hzY1LlFzlSCKadkptVLJmcl5qkxaOhuriEQT9OVoaU2mRFprzFnGXmvWo1ErlzC2csEbii5EHWnBhau3NqNc6AH3hWNsa6le6PQKZdXo08jtlTIZdZAOW9Psu2y5rBqzgbbaSmKJJAPTTlnxZU9HHacnHHnbd09PzMvz03Lh4MYGjg2mDVXq4c8i9HhCWR51sZctpY9do1JQbzFhMWrRalQggscfxO4N0VhVuTC5tp5EKsXgtCtrcZfQN+Vk74b6LHq0Vq3k/37gdbxqZxtwYWazlyIjrVKpuO+++3j9619PMpnkve99L93d3acEQfhn4Kgoio+Tlp76oiAIIukw/sMLb38EeDVwknSy7rAoij8v9f5Xvc4OEI/HZVH93t5eNm3axMDAAFdddVVeYx+adfOX9/1iyf46E201lYzl0TLrqrfij8RyCk3oNSo2NVWhUSkYmvHgWsj4GnVpxZNAJE5nnRm7N7TEI0J6Ca2pNGQlw6LxJHanB39cxBeKsqu9lqEZN3UWI9FEkiabiUAePblc2N9VT8/o0ox7oSYUOCecUVNh4NTYHPY8hBYJmdJPmSiWX5BQV2nAolOg1uoYnHHT1WDj1LiD/Rle3qzXUFNpoEKvlev0wYWBFfP+0BKKrFIBW5ts9E64aLSZcHhDxJMpVMo0hyC90ATlfAOAWaciKaY5/watim/e/Qau2nIuz3XkyBEOHjxY9PMsB729vXR0dOSlxF533XUcPXq04IKQB5emeIWEiYkJJicnZcno4eHhgnueMxMO2mstTDr9WeqpmRhzeHOKPezvbODUuCOv0EQ4lh3ub2y0pUUYPEFqLUZSyRS9Y3OY9Wo2NNjQLmiSh6MLyrCeYFYWfTF2NFXSO2YnJaZr/ceH55ha4Gfr1Qq2NFnRaDRp3nkOtZdChjbnzc+th3R9XK9R8bve9MijlmozDTYT/lCM/innklbTM1NuLHo1nnA6NFcKsGdDfUmGnr6fEHNeqK5IpFtNTVr+bHc7wUiMznoLDm8IXziWN3pYDL1GSZPFQO9EOnKZdgXkCGHxsApJGz+eTDIw7aKrzszQnI+//rNOalRh3G43lZWVq+7RJRTbs4uiuGbXXg2sibEPDg4SCARk2iucq5nn+7JuOriJmw5uIhJL8IfT4/zmxDBPdY8y78s2sHj0nLEIwKESPRKAzaSnulJPhU6DKKbr51qlgNPnp8lqZHQ+wJyv8F4/EypFerpJpmjCiZE56i16Zhc0ysPxlPwgK4QFZVeTjjlPkGmnn53thRN5U05/3sy3VqVkc3NVFtEnc3ijSaehq9GKUiEwNJNeaGKJFLva6zg6OINOraDZmnv0lFopUGcxLajKLFBQgyHcgTDecIJ5f4RwLIlIeos16fSlVXoabQQXpsQm8nScSTDrNdRbTUuqIj1jjnQf/qLFNVMbX61SYNDrefCvX82mBhtut1se96XX64nFYkSj0RVpH+ZDoed3uTMJLwTWJIwPhUIoFIqsD3/8+HE2b95cUiuhfHFR5KHDv2fIC091jzBs96FSCNgq9IRjCdprLZwcO2ecKoVAvdVEhV6NUkxhNBqIxJN4Q1HmPAHCGfs/s15Do81E35SL3W3VdI/No1cr2dxSjUqhZHjWlZepB2A2aGiqMi+ZAwewtcHMmZn8IhdwTn1FpVTIo6Ly7bUbbSbm3MEsL21ZaO4otO3JhEIQ6FwQk3T6Q2hVSgRBIBxNUGHQoBIgGosTCEfxReJ4QvGCdXibSYfFpGN4gUOfqZMHaeZgV31aK2/c4VvSbVhl0lNh0OQtw+1ur81qk12M2koDD3z8Jroalg4CDYVCnDhxAoPBQCKRkKmt5+v1jx49yr59+3KeQxRFrrvuOk6cOLGSU1+6YbxGo1nCmFuJgIUgCGxttLC9WcEbNlVgbWzndyfHGZp1MuHwE42nRwwHI2nZI7snJLc3puHJed7WajOxeEJWdg0sJLjC8aQ8REEhCGxqsmE16Zf0skuKr7kMHaB/zp/TM0lospkQIasRR69Rsau9Go1KycicJ4ukM+0KZBlTk81EKiWWZOhVJj3VZgMmgxqVIt09ZjPp0alVhGMJVAoFA1MuQnm2TrnQaDXBQjusBEl+Wsrch6IJucYP0F5rpsZswBtK05MFQShYb+8etS8ZfiGhyWbiu5+4mdaayiWvSdNZNBoNe/fuJZFILPH6VVVVVFVVLdvrF/Leq6gZv2ZYs9LbkguV0Pm2GPF4HLfbTXV1NXv37kUQBDrq0yu5NxTl6ZOj/ObECE+fHC25kWVnWy1DM05CGZrlQzNuNjVY6M+oS6dEkf4MLne9xUhLjRmVQuDsVGGvn0yJtNVW5jT2DTUm5rxhgotGF4VjCXpGM8Ql6i1Umw3pQQ6zbibn/agUAhvqrcx5AngXeOgSz72qQi93fEUTSbwLuYbFFNSuBitOfxh3IMKmRiv9027UKgXbW6sxajVMOH2LKgfZaLDoCEVjeBbx4N2BSMFqyajdx6jdR0t1BRqVEotRR12lgWG7N2etHsAXji5h6nXUVvLdT9xMvTV/iSvT8FQqFTU1NdTU1Mhe3+l0cvr0aRKJhNzGajabS/L6+Yw9EonknOt2MWHdJsIsV3RSYtpVVFTQ2Ni45EuuNGi59YrN3HrFZuKJJC/2T/PrE8P8pnuEyTyEi0Mb00SZXItwrrbHTMx6grTUmDkyMItWrWRPR12afTbrykl4OTnqWDK2eF9nPSdL4LiLYtprSp6z0Waiq8HKrrZanL4ACrOG5ioz7mAEuydUMs99R1sNQzNumS4sPdzxhb5vCS3VZuotRhzutIqstPVOD5zwEcwzYeXE8ByNVhPTee5l48IoJm8oKmfWFYJAi1VPjaUCTziWFS1MzPuzyE2bm2z898duotpceCuYzwNnjmFubW2Vvf7s7Cz9/f0YDAbZ+Jfr9S92SSpYZ2Mv1bPPzc0xNDTEzp07mZmZKfo+tUrJNdtauGZbC/e+83q6B6d46LdHOWWP0DM6h0qhYGdb/oGGAGemXNRV6pjzLvXYCiEtNiG9PxFNyY006YdVR43VzLwvItecI/EEuzpqZU9XrLRVadBQqVOhVYJWrUZQKokmWVA/DdBUVcHpMTuBaJzOBhsGjZp5fzivIORi7O+qp3vYnnV836STnW21WXkPWJzkU9PVYKPSoKFvypnX0CFNW7Ua1Uzn2F1sa6lm1O5dMoopJYpMuMNMuNPRh82ko622ktQCoef0xDxVFToabRXc/9E3YjEW14srtca+2OtLnWuSroK01y/F61/sKjWwjmF8KcYuiiKDg4PyXHW1Wo3dbl92+L+1tYY37W3kCwcPYvcG+dWLZ3j0j6dy6qiduzZUm7RLjN2oVbOh3pJ3oUg/rBEm3On3NdpMNFeZ0zX2SSdWk5aOOisv9k9TtxBuG7RqlAphoZ8+it0TxBuKyaH5YuxqsckZ98U03raayrRufDDKwExuCm+h8DoQiS6ZsZb9ehyVUsHxkfSi2dVgpapCz7wvlFMt5tSkm1abnnHXua3D7vZazkw6iSWK/46Zk12qzXo2NVjpqKvkk7degdlQujDkchNxUkOLyWSira0tp9ePx+N5M/xlz555IZWqYBifSCTo6enBaDRmzVVf7thmyG68iXidbDQn+dFn7yCFYqGsN8Jve0aWZIgH5vxYjFo5LK+3GNGqVXnFJBZDo1KgWujUM+k07GyrQa9TE46mBx/MeYJFp8tkQqtSsqnRSk9GiD0440StFGQhRUnfDdI1/s56a3rRnHETCMc4sDG/oQOMzHkLEncyI5KDGxs5MjAtJwarKvR01FnSSrzTTpntptHqgPR3u7XBzMkFDsJiSGSlarMeg1aDQkgvgNJIp3lfmC1NVdxzyyF06vQMAUEQUCgURSepnm+9e7HX9/l89Pb2curUKVKpFFarlerqasxmM4IglCxJlUO4Ivs7KSxc0Qp8h/TQRxG4URTF0VI/05qU3lKplDzcQcL09DTRaJSOjo4lxweDQbq7u+no6Fgy5XV8fBxBELLaAku5/p/+9CcsFguJRILt27cvqY+mUiLHh2f5TfcwvzkxItNRJULH5sYq7L4g7kWJOLNeg1mvwqRVYjaleezBSFpjLZMdJmXsg5E4kViceDLFxgYbZqOWWVcgbyOHhEKltcwxT/lg0Kq4ekszvlCMaZd/UZUiG/VWI/Pe0BKm36GNjbyYIbahUiiotxpznkutUrCxwYZJp2HK6aPOakKlUHB8eJY6i4EKnRqlAIl4nGgiSSCaxBtJFJzM87q9G/jPD7wOjSq94CeTSURRlBdyURRRKpXyAiAhFAoxPDzMjh07Cn5Hy0E8Hqe3t1fO8LtcLpxOJz6fD5/PxzPPPEMsFuOb3/xm3nMkk0k2bdqUJVzx4IMPsm3bNjkUFgThYeAXGcIV7xFF8c6F154GPi+K4q8XGmRSoiiWrJ19wcN4h8NBf39/3gEQmVNhSkUikSAYDNLQ0EB7e3vO+1EoBPZ3NbC/q4G/efPVDEw6+J8nn6XfHmJPexWxZLotMx1up/CH0+F2NjsstwFtaa5iJoPjLmmoZc51ywz3z05mM92KldbmfLGsDrbFMBs0NFhNPNUzSscCK7GpqiIvhXfWHWRrfQVnZtOfR6UQ2NVRl2XokB4YUWU25DT2eCKV1Wa6od5KMiWyuamKvsl5JkugDGfi1is28eW7XiMnTjO9eSqVQhRFksmkbPiJREI+Zq158SqVKku8YmRkhEcffZQ//vGPXHXVVfzjP/4jN9xww5Jz5BKueOyxx9i2bVvmYTmFKwRB2AaoRFH8NYDU+bYcrGsYn2ns0pfkdDo5ePBg3mZ/aSpMqZCy+BqNJmcUsRjSQ9NRb+Xe992S9grj0/zmxCDHxjz0TPhyNmPkQ66M+5jdmxV6Q/a8MZNOw8ZG24JyS5wJh78g3XTeH17idSXYjGmCjCQBZVwQ25hy+mUKr5R0UygEBqed+MJxZv1xjNr0+OUN9VZZVmsxukfm2NxUxdk8Aykhvbj94fS5Nk6jVs32RiuxSIQZTwRPkTLpzfva+fw7r81bIZEMWTI+yetL/ys9L8lkconXXynysecEQWDDhg1cf/317Nq1i7vvvjvvjIRcwhU5hkDmFK4ANgEeQRB+AnQAvwE+I4piyQmtC1J6kwY0arVa9u/fX/DHWE4W3+FwMDAwwM6dOzl58mTR4yVDl0o1Umlm59aN7Ny6kVgsxuycnd8cO8vzA3P0zoWw+/L3eOdLhNm9oSUMs0wEImlV3D0ddSgRseqVbG5qYMYdyBt+D866lyQcO2or8Ufi2DMoxifHHLTaDIy7zv0tkDG6SSHAxkYrVqMeg0bFlMtfVKQj39ZPKQjs3lC3RKcvGI3TnTGgorPeSrVZj9MfXhK93Hn9Nu66poNTp3pJpVLYbDZqamrkvXEuZHr9cDjMyMgInZ2dcgQgGfz5GH6pAyKkYabngZzCFaRt9VpgL+lW14eAu4D7Sz3xupfeQqEQ3d3dtLa25tToyve+QpD07Obn5zlw4EBRyV/pPdLik+9H1Gg0qFVKdjRWcMfrriIUCvHi6RGeOTVB93SAQbsfUUz3Ue9sry2YCJtxF25oObixkWOD04hius4tZdybqipoqqpY0tjiWuTdt7VUM+7w5uw91+u0pMd8L0VKhIFpN7WVUbRqFbFEkoMbGwpSePunXezdUJfFANSoFGxprs4bEUiQZMAkKTBbhZ6O2kpSosh121v52M1pGan29nbi8fQQh4mJCfx+PxUVFVRXV1NVVZVTEyESidDd3c2WLVuwWCyysadSKfm/RCKBIAgolcplGX6xJphQKCTLUOdDLuGKxTZQQLhiEjghiuLwwms/A67kQht7vj17KBTi+PHj7Nixo+QVsFg2PpVK0dvbi1KpLBolAFk/vuTN8x03NDREIBBg3759qFQqDAYDN15XzRuuPYDf76dvZIKnToww5g7zTH/hBpopp3+JgUA6I31gY0OW6myV2SAn8DLDb7NBI3PBB6Zd9E+7FoZdVNM7np+sc3banZd6CtBsMxJNiHJtXaoYZFJ4R+2erOrFnCcob02MWhWtNZUyAzAfzAYNtZVGzPr0IhpPpgiEYwzPefjLG/bxgdftzTperVZTV1dHXV2dnBGfn59nfHwchUJBdXU11dXVGI1G2dC3bt0qP1u5wn0p1Jf+k44r5vVLGRBRrPSWS7jihz/8YdYxBYQrjgAWQRBqRFF0kO5rP1rwgouwZp49U61GFEWmp6cJBAJce+21y2InKRSKvJ49Go1y4sQJGhoaaG1tLXquUg09mUxy6tQpdDodu3fvXnKcIAiYzWYO7d7Ood3bCYfDTM7M8tTxIV4cnqd7yr+ETgosaW/VqpRsWdS1Bul9cUu1eckMNV8oJntOlULBxiYbLVUVnJ6YL8rKy5fM66yrXJj7vjQvsYTCuzAfbn6BwntoYyP9006qKgxy27HNpKPGbMSgVREOB1FrtITjKezeIL5QbInctCDAP73jev78lYUz55Kcc2VlJZ2dnUSjUebn5+UFOR6Ps2HDhoLlLyncV6lUWV5fer4y9/iLDb8UlZqVCFds376dUoQrRFFMCoLwaeApIf1AHgO+XfCCi7AmpTdIK9RIe+JTp06hUCjw+/1cddVVyzpPKBTi7Nmz7N2bver7fD5OnjzJli1bqKqqWvK+5557Lksso1RDj0aj9PT00NDQICl/LguJRAKHY55ne4d4tm+G7ik/Exn75V3ttfSM2tOltQodg3nGGEmSVIWQGcZLfey+UJT+qdzkms5qA0Pz5+5lZ1sNA9PuvDoA+VBTaaCt2oylQodAmhfv9Ke5+IsZcoWgVAh86d2v5s1XbVnW9TMRDAbp6emhpaWFYDCI2+1Gq9XKXr9Uvnqx0t7c3ByJRCJvCfiTn/wk73nPe7jmmmtW8jEu3a43CeFwmO7ubpqammhubub5559f9jlyhfGScOWePXvyhk7SXl+lUuVMxOVCIBCgt7eXjRs35lxASoFKpaKhoZ63NdTzltek54f1DIzz25NjnJj0EQpHaLQaESGvoUPauzfZTDlls9VKgR1t2aWxTIprpUFLZ0OaXDMw7ZL38TGUMltuf2d61lwuuq1CYCHcVqMSUqSSKZIIhOIp5v3pCToapYIJp5/22kp6Rh1sbLSxvbWaMbs3b7dfJjQqBf/xvtdxw/7OosfmQyAQ4OTJk+zcuTPLo4dCIebn5zlz5gzxeBybzUZ1dXXBFtdipb1YLFawrHdZM+jcbjenTp3KOellOchM0En76Ew6bT5IP0wpiTiA+fl5BgcH2blz56r9aNL8sFdeYeP6Q7sJhUIMjU/xm5cG+NOwE29Qlbe1NJkSqbcuNXajVk1bbWXB8dTejIk0KqWCNpueaouZaU+QvRvqUCmVvDQ8Q731nDgFpMN2dyCC3RNk1hMk31rUYNERiKTpvZ5gBItRy6mMDH5LVQUNVRVpCm8OFV6dWsXX776BVy7oxa0Efr+f3t5edu3ateT3MhgMtLa2ys0uLpeLmZkZ+vr6MBqNstfPl8hdvNf3+XzMzs6ybdu2vHv9UpRlLzTWzNidTmfeSS/LgbRnl8p1Op2uqHCl9D6JXlnIm0NaQmtubo59+/aVlMlfCaSyXo2lgms31fKBW67D5fbydPcgzw/a6Z70L9GP6x6x02A1yR1t1QuS06fzJNoWQ60UsBpUVJiMKJUKmm0VmA0aookkGxtsnJ1yFmTW5cLmpiqmnD45WojGk1SbNFnefMLpZ2LhvIspvKmUyP/7yBu5cnPxSkw++Hw+Tp8+ze7du4uKoSwmwAQCAebn5+nu7gagqqqK6upqKioqcj4jwWCQU6dOsWvXLkwmU15Cj8vluqiHOsIa7tkl0clMPPfcc1x99dXLPtezzz6LSqWipaWlpHKdKIqcOHECi8VCQ0ND3ggglUrR399PIpFg27Zta6ofJpGIvF4vu3btykr2JJNJXC4XL5wa5vdnpuieCjDkSJf1JLnmlmozsURyCbfeoFFRZzFSadShUStJpdL0XbsniCsQkX9ApSCwZ0M9x4Zm5HNajDo6FySjM7nt+bCjrYbBHHt8Aair0DDrL0yWsZl0/NfHbmZne23J39tieDwe+vr62L1793n3j8diMZxOJ/Pz8wQCASorK6mursZms6FSqeR8wOJtQiZSqRTPPfcc7373u+np6aGmpmYlt7Iue/Y1M/ZEIrEki744aVYK3G63rBRaynZASrAEg0FmZmZwOp2o1Wq5qUGKNKQRVJWVlXR0dKypflgqlaKvrw9BENi8eXPBRUUqMfUNT/BUzyjdkz5iIiRFBRq1GrVKQTKZIrAw5qmQiIYESatOyqzXVhpwByJZs9DVKgWbGqswaFVMOHzyKCoJezfU0TvmyDmTHWBbcxWnJ/Oz6qoq9Hz3E7ewtaW66P3mg9vt5uzZs+zZs+e8I8bFkNSQ5+fncblcKBQKQqEQ27ZtK2jAR48e5WMf+xiPPfYYbW0r3pa8/Iz9hRde4MCBAyVL7U5NTTE+Pk4ymeQVr3hF0ePzJeLC4TAOhwOHw0EymcRisTA/P5+z8Wa1IS0qVquVtra2ZS8q4XCY8akZfts9xJFhFyemfHlbYXMhPSq5Ygm9NR/dVkJ6+o0RdyBCpUHDsaG5goMfIO35c43Sthk1/N0Nm9ixoZGamhoqKyuX/T24XC4GBgbYvXv3qhv6Ykh8kLq6Ovx+P9FoVO5ys1qt8mJ94sQJPvShD/GTn/yEzs6VJxp5ORr70aNH2blzZ9E6uyiKnD17lnA4zM6dO3nxxReLhv/LScSdPn1aFiO02WzU1tau6AEshmg0Snd3Ny0tLauyqMhlvZPpst6JaT+TrvyZ7xqzHoNOI0+xyYTVpCMaSxTVnrtiUyMD0y466i0kkyL9U86872mrrWTC4c1KyLVWm/mfe26l0WqUQ2av10tFRQU1NTXYbLaCiVZI538GBwfZs2fPqqrF5kIoFKKnp4ft27fLCbdkMonb7WZ+fh63283x48ex2+38/Oc/57HHHmPTpk3ne9mXn7GXojCbSCTo7u7GbDbT1dWFIAhF9/pSsqRYIm52dpaxsTF27dqFXq+X98p2ux2fz4fZbKa2thabzbYSof8sSGW8TZs2FaVRrgSpVLqsd6J/jN+dHOPYmJsBR0g2tOaqChLJ1JJwPBOFRC3yva5RKWmrMqJWCswH40vKbJkjoTrrrfzPJ29ZoheXyYZzOp0olUqqq6upqalZ8mw4HA5GRkbYs2fPmiVPJUil4m3btuXswJTu/fDhw3z+859Ho9EgCAI/+MEP5E62FeLSNnYpg56Jnp4eOjo68pYoJAngDRs2UF9fL//9+eef54orrljisZdDfZWSYzt37kSlWlqEEEURr9eL3W7H5XKh1+upqakpWKLJB5fLRX9/Pzt27FjzDG0sFpNDTm8oylMnhumd8tFnDzLnzT3hVUKFXoMASzrsBMia8FIIHXWV1FYaF2bQu6g2GwiEo2yot/LAJ26hqqJ4Ei0SiTA/P4/D4SAajcrNL7FYjPHxcfbs2VPU+58vJEPPpNvmQn9/P+9617v4wQ9+wM6dO/F6vRgMhvO9v5efsZ8+fZqGhoaciTan00lfX19O3vyLL77I3r17s77QUg09lUpx+vRpVCoVmzZtKinjLumRORwO5ufnEQSB2tpaampqimaAZ2ZmmJiYYPfu3esWci4mAcViMaZn53j6xCB/GnLkLOtJWKyNlybs1C7h8JeCqgodW5qraak287dvuRqzYfmfX4q2xsfH8Xg8VFdXU1tbS3V19ZoZfCQS4cSJE0UNfWRkhHe84x1897vfXcLoPE+8/Iz97Nmzcl0zE+Pj48zMzOTdkx07dowdO3bIr5XKiIvFYvT09FBbW1sSdz4fIpGInOCLx+NyyJlZm5U67zweT97oYTUhkUq2b9+eN+SEpWW9E1N+hh0BWVFHr1Gh16hwBSLo1Eq6Gm05k2wSBKDWYlyQklIjCAKRWAJPIMKsJ8i+znr+34dvLDoSuxBmZmaYmppi9+7d8nc/Pz8vN79I4f5q5FgkQ5c65fJhfHyc22+/ne985zurPkOOS93Yc0lTDQ4OUlFRQV1dnXxMX19fXukoCZl7/UxDL9ZffPLkSbq6upYsLueDRCIhh5yBQEDO0trt6bLWli1b1nzel7RN2LVr17In7Ph8Ps4MT/BU9wjHJ3z0TvvY1V5D/7SbRpuJs1OudN3cYqSqwoBeoyQUCpISlARjKeY8wbzCka/c0cY3PnSDPCtvJZienpYX/sXPg9T8Mj8/Tzgcxmq1UlNTg8ViWdF3XqqhT01N8ba3vY1vfOMbK+KJlICXn7GPjo6i0WhobGyURzlXV1fnlY6SIO31DQZDSYk4p9PJwMDAmu+ZU6mUvP2QhBZqa2upqqpaM88+NzfH2NjYqmwTpLLe77qHGJrzM+KOYvfHsPtCRbvoFuMN+zr5j/e/FrVq5YnNyclJ7HY7u3fvLpogXZwhL4UGm4loNCo7kUL8jdnZWd761rfyH//xH1x//fXL/kwl4uVn7BMTE4iiiM1mo6enh66uLmpri7OpTp06RUNDgxw2FzL0yclJZmZm2LVr15rvmTNLa/X19fj9fux2O06nE41GIxN5Vus+JiYmsNvt7Nq1a9X3r/F4HIdjnj+cHOKPZ2c4MeVnyl04wSfhzVdt4UvvfhXK84hoxsfHcTqdS9iFpWBxjgWQw32j0bjkeSnV0O12O295y1v48pe/zGte85rlf6jScWkbuyiKxGLZWd7p6WmcTid+v1/mGpdynv7+fkKhEM3NzVit1pzGLooiAwMDRCKRgluC1YK0TchXWguFQvI+XxRF2fBX0mQjiiLDw8MEAgF27NixLp/t+PHj+JMqnuuf48SUj7Mz/iXjnwH+/Pod/NM7rzuv/fPo6KhcKVmNLVAsFpPD/WAwKIf7VquVRCLB8ePH2bhxY8GS6Pz8PG95y1v4l3/5l5zikauMl5exi6JIb28vLpeLq666qmTpKKnLyOPxYLfb8Xq9S+rhyWSSkydPYjKZ6OzsXPPRuRJts9RtgvTw2e12IpEIVVVVJTPJRFHkzJkzCILAli1b1vyzSYuYlJmWvObQ+BS/OTHMsTEPJ6d8hGIJPvC6vXzmree3hx0ZGcHv97Njx441yXWkUik53He5XEQiEZqammhvb8/7DLrdbt785jfzuc99jptvvnnV7ykHXj7GLklHxeNxDAYDW7duLen9uTLumfVwp9OJTqcjGAzS1ta2LG35lWJ2dpbx8XF27dq1ItpmMpnE6XTicDjw+XxUVlbKTLLFHjuZTNLb24vJZGLDhg1rbuhShr9Q40c0GmVmzs7RM6M0V6rlmrjFYlnW/UnRSjgcXvMmJEgvuC+99BItLS3E43Hm5+dJpVJyuG8ymRAEAa/Xy1ve8hY+/elP8+Y3v3lN7ykDl7axQ/rByJSOMpvNTE1NsX379oLvK5UR5/V6OXnyJBaLhWAwKO+Ta2trV51tJYoiY2NjuFwudu3atSoJOKn5QiLyGAwGuaYM6cRkXV3dihRzlgupm2w5GX6prOdwOLIosMUSlNKYr1gsxrZt29Z8EYvFYjJZK7MyIxn9/Pw8o6OjPPTQQ4yPj/Pxj3+cO++8c03vaREufWN3OBxZ0lF+v5+RkRF27dqV+2IlEmUgnTwZHh7OejhDoRB2ux2Hw4EgCLLhn28rpMTVTyaTbN26dU28kBQuS/cfCoWora2ls7NzzRs/JO75+TSZSGU9h8ORt9NQOq6/v59UKrUu25J4PM7x48eXGPpiuN1u7rrrLtRqNbOzs/zFX/wF99xzT97jVxmXtrFLI5i2bNkiJ6XC4TB9fX052UfLob5KHnbnzp15s9LRaBSHw4HdbieRSMhMrFzZ2UKQ8gFms3nNW2Hh3J5ZklJ2OBwkEgmqqqqora2Vw83VglTKW23u+eJOQ4lMNT09Lbf6rpehd3R0FGxTDYfD3H777bzjHe/gfe97H5Cuwa/1IpuBS9vY4ZzoZOa/u7u7lzCQSmXEpVKprGRVqR5WCtfsdjvhcLjkBJl0v01NTTQ2NpZ0rfOB1+vl9OnT7NixI6t/QLp/h8NBMBjM2iefT5QxPT3N9PQ0u3fvXlPuuXT/Q0NDJBIJ6urq5DzFWu3V4/E4J06coL29vaChRyIR3vnOd3Lbbbfxl3/5l+e9AL33ve/lF7/4BbW1tfT29i55XRRFPv7xj/Pkk09iMBh44IEH2Ldv38vP2JPJJEeOHOHKK688d4ESDT0ej9PT00N1dTWtra0r/lEWJ8gsFgu1tbVZfcpwzsOej/jkciARgYopsKRSKXmf7PF4qKiokIk8yynJnU9de7kQRZHTp0+j1Wrp6OjA6/XicDjkPMVKG47yQSqvtbW1FeRxxGIx7rzzTl772tfy0Y9+dFUijd///veYTCbe9a535TT2J598kq997Ws8+eSTvPDCC3z84x/nhRdeuPTVZTO14yF7lDKck++VXssHqeFjw4YNJZFwCkGpVMqaZFKbqN1up7+/XzYcpVIpd62th4ig1DxTigZe5nAEaZ8s5S+0Wq3csJPvPFIHYCAQYPfu3WueBU+lUpw6dQqj0Si3gdpsNmw2WxYZpru7G0EQssgwK4Fk6K2trQWflXg8znvf+16uv/76VTN0gOuuu47R0dG8rz/22GO8613vQhAErrzySjweD4IgNIiiWLzF8DyxbuOfgKzyWamJOLfbTV9fX9GGj5VAUn+VHjyfzycPm7Rarfj9fnQ63ZqGuOPj48zPz8tTZ5aDzMEJGzduzDIcQE5QSglMiXiUSCTYuXPnmu+ZpZJrRUVFziGbgiBgMpkwmUx0dHTI3HeJHLUcPgKkDf3EiRO0trbK/Rf5jvvABz7A/v37+dSnPrXm30Mmcg137O/vbwJeXsYuoVRDn56eZnJykr17965LssTj8ZBKpbj++uuJRCLY7XaOHz+eFQ2sFvVVKj9FIhH27NmzKh7WaDRiNBppb2+XDefs2bNEo1GqqqoIBoPodDq2bt26Lobe09Mjy3GVAq1WS1NTE01NTXJZb3p6mjNnzmA2m+V9fq5FUTL05ubmgoaeTCb5q7/6K7Zu3crf/d3frauhX2iseRifCVEUEUWRmZkZampq8noySR8+GAyyf//+ddlTSqU1yfAkj7NhwwY5s3zy5EmZ+prpMZcLKdGoUqnYsWPHmjxwmYYj1ZmTyaRcEVnLBFkymZTzKyslOimVSrl0l1nWGxkZkfkU1dXV6HQ6ksmknEjNFD3JdV8f+9jHaG5u5t57770ghp5ruCPpSa1rjnXz7FIibvv27bI8lF6vXyJMII2L0uv17Nq1a81/kEyWWr5ykF6vl4cOxGIxHA4HZ8+eJRaLUVVVRV1dXcklMckQLBZL0W6/1YD0fdbX19Pa2irnKaTx1kajUU7wrcZ2RTK82traVSMDZW5Xurq65Ikvp06dIpFIEI/HaWxsLGjoqVSKT33qU1itVr7whS9cMI9+yy23cN9993HHHXfwwgsvSJTkNQ/hYY2z8ZIOXS4xSCk5Mzc3x/z8PBqNBpvNxuzsLM3NzSXpw58vpNJaY2Pjiq4n9bbb7Xa5JFZbW5uXOhqPx+nu7qahoWFdPp8U2jY2NuYsHUpDEyTqsUqlykmEWc71pM+3HqXKZDLJ8ePHMRqNJBIJWV9AanrJHOf0mc98BoD/+3//75omJd/xjnfw9NNPMz8/T11dHf/0T/8kd3/efffdiKLIRz7yEQ4fPozBYOC///u/OXDgwKVfekskEsRisZL25w6Hg9OnT6NWq+Ws8mrukRdjtUtrUklMataprKyUm3UUCoU8UnjDhg0rHSSwLEihe1tbW8E9bCYWE2GWQ0SSFpampqY1l+eGcxFEXV2dvHBmRi0ulwulUsnRo0cZGxsjGo3yrW99a82rDyvEpW/s3/3ud9mwYUNO1ZFMOBwOhoaG5Dlr0kMnqb+sFu1Vgsfj4cyZM2tWWhNFUS7puVwutFotgUCA7du3r0vNXlpYurq6Vny9xUSkQg0vEoGlWBZ8tZDL0BdDGhN+zz33cOzYMTZt2sTdd9/NHXfcseb3twJc+sb+05/+lB/+8IecPXuWV7/61dx6660cPHgwK5SfmJjA4XCwc+fOnLXhaDSK3W7HbreTTCapqamhrq5uxckxu93OyMjIugwbgHMDLm02G36/H41GU7QWfj6QOAnFpJaWg8UNL1JmvKqqimQyKTPVzpcDUQpSqRTd3d3U1NQUzAmIosiXv/xlBgcH+e53v8v8/Dyzs7Ps2bNnze9xBbj0jV1COBzm8OHDPPLII3R3d3P99dfzxje+kV/84hfccccd7Nu3r6TwSkqO2e12YrEY1dXV1NXVlcx3Hx8fx+FwrInSSy44HA6Gh4ezFpbFzTqlqtaWAmmE8VqSgaQWY0kVJhKJ0NjYSEdHx5rrukuGXizLL4oi//mf/8nx48f54Q9/uCq/9eHDh/n4xz9OMpnk/e9/v5wDkDA+Ps673/1uPB4PyWSSL33pS9x4442lnv7lY+yZiEaj/OxnP+PTn/40tbW17N27lze/+c1cc801y/pRcvHd6+rqck7jlDqt4vH4uvROQ5ojICmk5jOCTNXazGadlejmeb1ezpw5s6ojpwtBal1ubm6WxTlWs9NwMaS6fVVVVVFD/+Y3v8mzzz7Lj3/841VZgJLJJJs2beLXv/41zc3NHDx4kAcffJBt27bJx3zwgx9k7969fOhDH+L06dPceOONBZl0i3Dp02VzQavVcurUKf7P//k/3Hzzzfzud7/j0Ucf5a//+q85dOgQt912G9dff33RH0mtVtPQ0EBDQwPJZJL5+XnGxsYIBALYbDbq6uqorKyUWVxGo5FNmzatS8lldHQUt9vNvn37CuYqdDodLS0tsqCClLuQFq/a2lrMZnPRe5bUZldjsmkpkFRZMyW5JAacw+Ggr69PLkvW1tbmHYdcKlKpFCdPnsRmsxU19Pvvv5+nn36aRx99dNUijRdffJGuri6Z7nvHHXfw2GOPZRm7IAj4fD4gvfCuRzViuVh3z54PiUSCP/zhDzz88MM888wz7N27l9tuu41Xv/rVy9pbZ4508nq9cpfVxo0b19yjr1YEITXr2O12/H4/VqtVLuktPqdENFmPoRRwbnJKsZxAIpGQP0O+klgpkAzdYrEUZeI98MADPPbYYzz22GOrmo955JFHOHz4MN/5zncA+N73vscLL7zAfffdJx8zMzPD6173OtxuN8FgkN/85jfs37+/1Eu8PD17PqhUKl71qlfxqle9imQyyXPPPccjjzzCP/3TP7Ft2zZuu+02Xvva1xZNzEnMK6PRSHd3N83NzUQiEZnAkFkOW01IDR86nY7t27eflydb3KzjdruZm5vj7NmzWV1udrtdphOvRw5CSv4Vm5wC6d+zrq6Ouro6+TM4HA76+/sxmUwlSW5LUVllZWVRQ//BD37Ao48+ys9//vP17EOX8eCDD3LXXXfxqU99iueff54777yT3t7ei6rUd9EYeyaUSiXXXnst1157LalUiiNHjvDwww/zpS99ia6uLm655RZuuOGGvEkoqS88s3lGKofNzc0xMDCw4tbQXEgkEvJ+8jxmdOeEQqGgqqqKqqqqrORYX18foiie76jgkhEMBuXppsttSFr8Gfx+Pw6Hg9HRUdRqtZykzIxMMpto2tvbC57/4Ycf5vvf/z5PPPHEiqs0hZCL4rq45Hf//fdz+PBhAK666ip5ft16VChKxUUTxpeCVCrFiRMneOSRR/jlL39JS0sLt9xyCzfeeKMcUs7NzTE6OipPas2FzNZQp9OZpf223M4ziYXX3Ny8LmQSOJcT2LBhgzwGWaVS5TSa1cBaZvnD4bBcnRBFUW5xHRkZkbvhCuFnP/sZ3/zmN/nFL35RNNpYKRKJBJs2beKpp56iqamJgwcP8sMf/jBLS/ENb3gDt99+O3fddRdnzpzhNa95DVNTU6VGeC/PbPxqQZKmfuSRR3jiiSfkh0Sj0fDv//7vJYe1EmVUou3qdDrZaIqdQ9q/rvaIqUL3Ojg4SDQaXZITWGw059usI0Ey9EKKs6sFqbQ6NDSEKIo0NjYWbHF94okn+I//+A+eeOKJgsMeVgNPPvkkn/jEJ0gmk7z3ve/l7//+7/mHf/gHDhw4wC233MLp06f5wAc+QCAQQBAEvvzlL/O6172u1NOXjb1UJJNJPvjBD/LSSy+h1WqpqKjglltu4eabb6ampmZZ++dMvr7kLXOp1Uqyy9u2bVszj5IJURTp6+srSb8tFovJhi/xEVaiX+fz+Th16hS7du1al3KepGij0+lob2+XE625JLd/9atf8cUvfpEnn3xyXViJa4yysZcKv9/P/fffz8c+9jEEQWBoaIhHH32Uxx57DI1Gwy233MKtt95KfX39sh72TAKMQqGQvWU4HObs2bPrVtOWxk7rdLplD8HIbNYJhUJys04xQQipbr/c4ZErRaZ01eLPKOVbHA4HDz30EM888wyzs7P88pe/ZPPmzWt+b+uAsrGfL0RRZHx8nEcffZSf/exnpFIpbr75Zm677Taam5uXZTSSmMXU1BThcJjW1laamprWvK4tqdtKLbHne67F3jJXdULSkF+vur009UatVtPV1VXwd3nmmWe49957efWrX83vf/97PvvZz3LTTTet+T2uMcrGvpqQRDMeffRRfvrTnxIOh3njG9/IrbfeWvK0lcnJSWZnZ9m6dStut1uWqZY8/mp7ealldC0GRSzuEJPKYUqlksHBQfbs2bMuJSxpe6JSqYoa+vPPP8+nP/1pfvGLX8jZcEmodCUoRoEF+PGPfywLXezevZsf/vCHK7pWEZSNfS1ht9v56U9/yk9+8hNcLhc33ngjt912W06WnSTS6PP52LlzZ1apTmK+zc3NEY1GZcM/X313qZNsPbL8UjlsbGwMu92O1Wqlvr5+VRVf8123r68PpVLJxo0bC35fR48e5WMf+xiPP/44ra2t533tUiiwAwMDvP3tb+e3v/0tVqsVu92+VqW0srGvF5xOJ4899hiPPvoos7OzvP71r+dNb3oTW7dulRlcWq22qFb94v3xciivmZB45+vV+w7Iuu579+4lHo9jt9uZn5/PylWspqeXpMAEQShKYz5x4gQf+tCH+OlPfypTVs8Xzz//PPfeey//+7//C8AXv/hFAD772c/Kx/zN3/wNmzZt4v3vf/+qXLMALi8G3YVEVVUV733ve3nve9+Lx+Ph5z//OV/4whcYGhpCpVLxmte8hn/4h38oyoZSqVTU19dTX18vU14nJibw+/1FVWwkSOW8fKOg1wIS5Xbv3r1oNBo0Gg0dHR10dHTIzTqnTp2SW4zPd8si0YqBoobe29vL3XffzSOPPLJqhg65VV5feOGFrGOke7zmmmtIJpPce++96zG+ec1QNvZFsFgs3Hnnndx2223ceuuttLe3Mzo6yjXXXMNrXvMabr31Vg4cOFDU8BdTXiWl1L6+vqKDKUqho64W5ubmGB8fz0u5zdWsI0k9S9yG5UQukpy1KIpFS4hnzpzh/e9/Pz/60Y/YtGnTij/jSpFIJBgYGODpp59mcnKS6667Tk6WXoooG3seKBQKPvOZz8jECKkn/9vf/jYf/ehHuf7667n11lu58sori9JtMwc7SImxubk5+vv75VnzarVablFda/KKhNnZWSYmJtizZ09JJCS1Wi3r2UmdhuPj43KjS75mHQkSKSiZTBYd6tjf38973vMefvCDH2Tto1cLpVBgm5ubueKKK1Cr1XR0dLBp0yYGBgaWjC+7VLAqe/ZiWc1oNMq73vUujh07RlVVFQ899NB5l5EuJCKRCL/+9a955JFHOHbsGFdffTVvetObuOaaa5ZFt5W47pKohs1mo7Gxkerq6jWXz5bmvO3Zs+e8x09LjS52ux2PxyMvYJkz5yVDj8fjRXXrR0ZGeOc738kDDzyQcwjoaqAUCuzhw4d58MEHZaWbvXv3cuLEibUg8VwaCbpSsprf+MY36Onp4Vvf+hY/+tGP+OlPf8pDDz20Crd/4RGLxfjd737HI488wvPPPy/35F933XUlZbKlxNiuXbtk5pvT6ZRltgvp668UU1NTskTTai8q0gKWOXO+pqYGv98vt/4WMvTx8XFuv/12vvOd76y5By1GgRVFkU996lMcPnwYpVLJ3//936+Vht2lYeylZDVf//rXc++993LVVVeRSCSor6+XZZleTkgkEvz+97/n4Ycf5g9/+IPck/+qV70qZyY737jkXDLb0uTT821llTT/du/evS7DNwKBAGfPniUQCMgeP1+zztTUFG9/+9v5+te/ztVXX72m93aR4dLIxpeS1cw8RqVSUVlZidPpXJfmkfWESqXi1a9+Na9+9atJJpP88Y9/5NFHH+Xee+9l+/bt3HbbbfzZn/0ZBoOBU6dOEYlEcs54y5yB1tnZSTAYlEdRnU932/j4OC6Xa9XGTRWDIAg4HA50Oh379++XVYN7enoA5ASmXq9ndnaW22+/nf/8z/+83Ax93VBO0K0RlEol1113Hddddx2pVIoXX3yRRx55hC9+8YuYTCaUSiUPPfRQSSG60WiUS2FSd1tPT48sWFlKDXx0dBSv18uuXbvWTVBhZGSEYDAoj7gyGAy0tbXR1tYmS1j98Y9/5G/+5m9IpVL89V//Nddee+263NvliPM29lKymtIxzc3NJBIJvF7vy6FTqWQoFAquvPJKrrzySv75n/+ZZ599lr1793LjjTfS2toq9+SXUm7T6/VZBmO32+UauGT4ixtXhoeHCQQC7Ny5c10N3e/3551lp9VqaW5uRqfTYbFYuP766/nVr37FsWPH+Na3vrUu93i54bz37KVkNb/+9a9z8uRJOUH3k5/8hB//+MercPuXHg4fPsxrX/talEql3JP/8MMP8+STT1JTU8Ott97KTTfdtGxCjdQLPjc3Rzwel8c4zc3NEYlEzlsqazkYHR3F5/OxY8eOgouL2+3mzW9+M5/73Oe4+eabgfPjukNpfHeARx99lLe+9a0cOXKEAwcOrPh6q4RLI0EHxbOakUiEO++8k+PHj2Oz2fjRj35UlA1V7Ef793//d77zne/I88n+67/+a9UlodYTEk/8kUcekVVXbrnlFm666aZl9+RL5JeRkRHi8ThNTU15ZbZXG2NjY3g8nqJRhNfr5S1veQuf/vSnefOb37wq1y6lMgTplug3vvGNxGIx7rvvvrKxL+CCcONL+dF+97vfccUVV2AwGPjmN7/J008//bIp50kjq6WefK1Wy80331xyT75ER02lUnR1deFyuZibmyMYDMp8/WL97CvB+Pg4bre7qKH7/X7e+ta38pGPfITbb7991a5fSmUI4BOf+ASvfe1r+cpXvsJXv/rVy8bYLx7pywxk6nRrNBpZpzsTr3rVq+S96ZVXXinNuX5ZQBAEurq6+Nu//Vv++Mc/8sADDwBw1113ccMNN/C1r32NiYkJci3UUoQAsGXLFtRqNXV1dezatYtDhw5htVqZnJzkT3/6E319fbhcrpznWS6kTH8xQw8Gg9xxxx188IMfXFVDh9yVoamp7NHnL730EhMTE7zxjW9c1WtfCrgos/GllPMycf/99/OGN7xhPW5t3SEIAm1tbdxzzz188pOflHvy7777biKRCDfddBO33norHR0d8tQUk8mUszdcktmuqamRWW+zs7OcPXuWyspK6urqlq3rDunavdPpZPfu3QXfGw6HueOOO7jzzju58847V/R9nA9SqRT33HOPvHhebrgojX05+P73v8/Ro0d55plnLvStrDkEQaCxsZGPfvSjfOQjH5F78u+55x7cbjcqlYpXvvKV/P3f/33REH2xvLNEd+3v76eiooK6urosums+TE5OMj8/X9TQI5EIf/7nf87b3vY23vOe96zo8xdDscqQpBv4yle+Ekj3Btxyyy08/vjjF0Mov+a4KPfspe69fvOb3/DRj36UZ5555qLS515vJBIJbr/9dlKpFLFYjLm5uaye/OXszTPprk6nE6PRSF1dXU6+/uTkJHa7vSgbLxaL8Rd/8Re87nWv46Mf/eiaJQlLqQxl4pWvfOVltWe/KD37wYMHGRgYYGRkhKamJn70ox8tkQM6fvw4f/mXf8nhw4cva0OHdHh866238q53vQtIa8g9/vjj/Ou//itjY2O89rWv5bbbbiuJUCMIAhaLBYvFIivYSGOudTqdbPjSGO1ihh6Px3nPe97DK1/5yjU1dEgzGO+77z5e//rXy5Wh7du3Z1WGLmdclJ4dipfz/uzP/oyTJ0/Kkk2tra08/vjjBc95idZgzwt+v58nnniCRx99lP7+frknf//+/cvemwcCAex2O9PT0yQSCTo7O6mrq8vb8JNIJHjf+97H3r17+exnP/uy64VYRVy+pbe1wCVcg101hEIhfvnLX/Loo4/S29sr9+RfccUVJTfFzMzMMD09zaZNm+TBjZlCHRJfP5lMcvfdd9PV1SULNpaRF2VjX01cwjXYNUFmT/5LL70k9+RfffXVefn6kqEvbo2VZLbtdjtut5vnnnuOiYkJ2tra+MIXvlA29OK4fOvsa4FyDTYbOp2Om2++me9+97scO3aMN73pTTz66KNcffXVfPSjH+Wpp54iFovJx8/OzjI1NZVzj67T6WhtbeXAgQPs3r2bvr4+/vSnP/HMM8/w3//93+v90crIg4syQXchcDnXYDUaDTfccAM33HBDVk/+Zz/7Wfbt20ddXR1+v58vf/nLRUcsf/WrX6WlpYWf/vSneDweRkdHz/v+Ljfq9JpBFMVC/71s8Nxzz4mve93r5H9/4QtfEL/whS/I//Z4PGJVVZXY1tYmtrW1iVqtVmxoaBCPHDlyIW73okAikRD/5V/+RWxubhb37Nkj3nHHHeKDDz4oOhwOMRgMZv3n9/vFT37yk+IHPvABMZlMruo9bNiwQRwaGhKj0ai4a9cu8dSpU1nH/Pa3vxWDwaAoiqL4jW98Q3z729++atdfJxSzw1X577Ix9ng8LnZ0dIjDw8PyQ9Pb25v3+Ouvv/6yNnRRTH9n73nPe0SPxyMmk0nx+eefF++55x5x165d4lve8hbxe9/7njg3NycGAgHxb//2b8V3v/vdYiKRWNV7KLZIL8ZLL70kXn311at6D+uAdTH2y2bPnlmD3bp1K29/+9vlGmyxkl0hHD58mM2bN9PV1cWXvvSlnMf8+Mc/Ztu2bWzfvp13vvOdK77WekOlUvFf//VfVFZWyj35//Zv/8bx48f57Gc/S29vL69//es5dOgQZ8+e5f777191qatSci2ZeDlTp88bRVaDMgqglBCzv79f3LNnj+hyuURRFMW5ubkLcatrhmQyKT722GOi3+9fk/M//PDD4vve9z753//zP/8jfvjDH8557Pe+9z3xiiuuECORyJrcyxqi7NkvdpTSnfftb3+bD3/4w1itVoCXHdtPoVBwyy23rJnWfSlKSJCmTn/+85/n8ccfX7Y23+WCsrGfB0oJMfv7++nv7+eaa67hyiuv5PDhw+t9m5c0MqnTsViMH/3oR0torxJ1+vHHH3/ZLaariXLpbY3xchshtN4ohe/+13/91wQCAd72trcBpVGnL0eUjf08cDmOELoQuPHGG7nxxhuz/vbP//zP8v//zW9+s963dEmiHMafB0oJMW+77TaefvppID39pb+/f1WnkZZRRqkoG/t5oJRy3utf/3qqqqrYtm0br3rVq/jKV75Skox2sZLe+Pg4r3rVq9i7dy+7du3iySefXPXPV8bLC5dNI8ylhFI69D74wQ+yd+9ePvShD3H69GluvPHGVaGmlnFBUG6EuVxRSklPEAR8Ph+QlmVubGy8ELdaxiWEsrFfhCilpHfvvffy/e9/n+bmZm688Ua+9rWvrfdtrhjFtijRaJTbb7+drq4urrjiinLEskooG/sligcffJC77rqLyclJnnzySe68805SqdSFvq2iSCaTfPjDH+aXv/wlp0+f5sEHH+T06dNZx9x///1YrVYGBwf55Cc/yd/+7d9eoLt9eaFs7BchSinp3X///bz97W8H4KqrriISiTA/P7+u97kSlLJFeeyxx3j3u98NwFvf+laeeuqpVdG2v9xRNvaLEKWU9FpbW3nqqacAOHPmDJFIhJqamgtxu8tCKVuUfCO+yzg/lI39IkQpJb1/+7d/49vf/ja7d+/mHe94Bw888EBZ/qmMgihWeivjZQRBEP4LuAmwi6K4I8frAvCfwI1ACLhLFMWXVvkergLuFUXx9Qv//iyAKIpfzDjmfxeOeV4QBBUwC9SI5Yf1vFD27JcXHgBuKPD6G4CNC/99EPjmGtzDEWCjIAgdgiBogDuAxUT2x4F3L/z/twK/LRv6+aNs7JcRRFH8PeAqcMitwP8s9Fj/CbAIgtCwyveQAD4C/C9wBvixKIqnBEH4Z0EQpMTE/UCVIAiDwD1AboH/MpaFciNMGZloAiYy/j258LeZ1byIKIpPAk8u+ts/ZPz/CPC21bxmGWXPXkYZlw3Kxl5GJqaAlox/Ny/8rYyXAcrGXkYmHgfeJaRxJeAVRXFVQ/gyLhzKe/bLCIIgPAi8EqgWBGES+EdADSCK4rdI76NvBAZJl97WZpB6GRcE5Tp7GWVcJiiH8WWUcZmgbOxllHGZoGzsZZRxmaBs7GWUcZmgbOxllHGZoGzsZZRxmaBs7GWUcZmgbOxllHGZ4P8HbELK0z1V54oAAAAASUVORK5CYII=", + "text/plain": [ + "
" + ] + }, + "metadata": { + "needs_background": "light" + }, + "output_type": "display_data" + } + ], + "source": [ + "n = 10\n", + "coords, elems3, dirichlet, neumann = maillage_carre(n)\n", + "show(coords, f(coords[:, 0], coords[:, 1]))\n", + "# show(coords, u_ex(coords[:, 0], coords[:, 1]))\n", + "# print(maillage_carre(3))" + ] + }, + { + "cell_type": "code", + "execution_count": 6, + "metadata": {}, + "outputs": [], + "source": [ + "def raideur(triangle):\n", + " M = np.zeros((3, 3))\n", + " x = triangle[:, 0]\n", + " y = triangle[:, 1]\n", + "\n", + " # calcul de alpha\n", + " mat_alpha = np.array(\n", + " [\n", + " [x[1] - x[0], x[2] - x[0]],\n", + " [y[1] - y[0], y[2] - y[0]]\n", + " ]\n", + " )\n", + " alpha = np.linalg.det(mat_alpha)\n", + "\n", + "\n", + " for i in range(3):\n", + " grad_eta_i = np.array(\n", + " [\n", + " y[(i+1)%3] - y[(i+2)%3],\n", + " x[(i+2)%3] - x[(i+1)%3]\n", + " ]\n", + " )\n", + " for j in range(3):\n", + " grad_eta_j = np.array(\n", + " [\n", + " y[(j+1)%3] - y[(j+2)%3],\n", + " x[(j+2)%3] - x[(j+1)%3]\n", + " ]\n", + " )\n", + "\n", + " M[i, j] = np.dot(grad_eta_i, grad_eta_j)\n", + "\n", + " return M / alpha / 2" + ] + }, + { + "cell_type": "code", + "execution_count": 7, + "metadata": {}, + "outputs": [ + { + "ename": "IndexError", + "evalue": "tuple index out of range", + "output_type": "error", + "traceback": [ + "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m", + "\u001b[0;31mIndexError\u001b[0m Traceback (most recent call last)", + "\u001b[1;32m/home/laurent/Documents/Cours/ENSEEIHT/S8 - Équations aux Dérivées Partielles/TP-EDP.ipynb Cell 8'\u001b[0m in \u001b[0;36m\u001b[0;34m\u001b[0m\n\u001b[1;32m 60\u001b[0m coords, elems3, dirichlet, neumann \u001b[39m=\u001b[39m maillage_carre(n)\n\u001b[1;32m 62\u001b[0m A \u001b[39m=\u001b[39m assemblage(coords, elems3)\n\u001b[0;32m---> 63\u001b[0m b \u001b[39m=\u001b[39m second_membre(coords, elems3, f)\n\u001b[1;32m 64\u001b[0m U_d \u001b[39m=\u001b[39m calcul_Ud(coords, dirichlet)\n\u001b[1;32m 65\u001b[0m b \u001b[39m-\u001b[39m\u001b[39m=\u001b[39m np\u001b[39m.\u001b[39mdot(A, U_d)\n", + "\u001b[1;32m/home/laurent/Documents/Cours/ENSEEIHT/S8 - Équations aux Dérivées Partielles/TP-EDP.ipynb Cell 8'\u001b[0m in \u001b[0;36msecond_membre\u001b[0;34m(coords, elem3, f)\u001b[0m\n\u001b[1;32m 22\u001b[0m mat_alpha \u001b[39m=\u001b[39m np\u001b[39m.\u001b[39marray([[x[\u001b[39m1\u001b[39m] \u001b[39m-\u001b[39m x[\u001b[39m0\u001b[39m], x[\u001b[39m2\u001b[39m] \u001b[39m-\u001b[39m x[\u001b[39m0\u001b[39m]], [y[\u001b[39m1\u001b[39m] \u001b[39m-\u001b[39m y[\u001b[39m0\u001b[39m], y[\u001b[39m2\u001b[39m] \u001b[39m-\u001b[39m y[\u001b[39m0\u001b[39m]]])\n\u001b[1;32m 23\u001b[0m alpha \u001b[39m=\u001b[39m np\u001b[39m.\u001b[39mlinalg\u001b[39m.\u001b[39mdet(mat_alpha)\n\u001b[0;32m---> 25\u001b[0m b[triangle] \u001b[39m+\u001b[39m\u001b[39m=\u001b[39m alpha \u001b[39m/\u001b[39m \u001b[39m6\u001b[39m \u001b[39m*\u001b[39m f(centre[\u001b[39m0\u001b[39;49m], centre[\u001b[39m1\u001b[39;49m])\n\u001b[1;32m 27\u001b[0m \u001b[39mreturn\u001b[39;00m b\n", + "\u001b[1;32m/home/laurent/Documents/Cours/ENSEEIHT/S8 - Équations aux Dérivées Partielles/TP-EDP.ipynb Cell 5'\u001b[0m in \u001b[0;36mf\u001b[0;34m(x, y)\u001b[0m\n\u001b[1;32m 1\u001b[0m \u001b[39mdef\u001b[39;00m \u001b[39mf\u001b[39m(x, y):\n\u001b[1;32m 2\u001b[0m \u001b[39m# return 2 * np.pi ** 2 * np.sin(np.pi * x) * np.sin(np.pi * y)\u001b[39;00m\n\u001b[0;32m----> 3\u001b[0m \u001b[39mreturn\u001b[39;00m np\u001b[39m.\u001b[39mones(x\u001b[39m.\u001b[39;49mshape[\u001b[39m0\u001b[39;49m])\n", + "\u001b[0;31mIndexError\u001b[0m: tuple index out of range" + ] + } + ], + "source": [ + "def assemblage(coords, elems3):\n", + " Ns = len(coords)\n", + " A = np.zeros((Ns, Ns))\n", + " for triangle in elems3:\n", + " M = raideur(coords[triangle])\n", + " for i, a in enumerate(triangle):\n", + " for j, b in enumerate(triangle):\n", + " A[a, b] += M[i, j]\n", + " return A\n", + "\n", + "\n", + "def second_membre(coords, elem3, f):\n", + " Ns = len(coords)\n", + " b = np.zeros(Ns)\n", + " for triangle in elem3:\n", + " coords_triangle = coords[triangle]\n", + " centre = np.mean(coords_triangle, 0)\n", + "\n", + " # calcul de alpha\n", + " x = coords_triangle[:, 0]\n", + " y = coords_triangle[:, 1]\n", + " mat_alpha = np.array([[x[1] - x[0], x[2] - x[0]], [y[1] - y[0], y[2] - y[0]]])\n", + " alpha = np.linalg.det(mat_alpha)\n", + "\n", + " b[triangle] += alpha / 6 * f(centre[0], centre[1])\n", + "\n", + " return b\n", + "\n", + "\n", + "def calcul_Ud(coords, dirichlet):\n", + " Ns = len(coords)\n", + " U = np.zeros(Ns)\n", + " # for d in dirichlet:\n", + " # x, y = coords[d].flatten()\n", + " # U[d] = u_d(x, y)\n", + " U[dirichlet.T] = u_d(coords[dirichlet, 0], coords[dirichlet, 1])\n", + "\n", + " return U\n", + "\n", + "\n", + "def tildage(A, b, coords, dirichlet):\n", + " A_tild = np.delete(A, dirichlet, 0)\n", + " A_tild = np.delete(A_tild, dirichlet, 1)\n", + " b_tild = np.delete(b, dirichlet, 0)\n", + " coords_tild = np.delete(coords, dirichlet, 0)\n", + "\n", + " return A_tild, b_tild, coords_tild\n", + "\n", + "\n", + "def untildage(x, dirichlet, U_d):\n", + " x_untild = np.zeros(U_d.shape[0])\n", + " not_dirichlet = np.setdiff1d(range(n*n), dirichlet)\n", + "\n", + " x_untild[dirichlet] = U_d[dirichlet]\n", + " x_untild[not_dirichlet] = x\n", + "\n", + " return x_untild\n", + "\n", + "n = 50\n", + "coords, elems3, dirichlet, neumann = maillage_carre(n)\n", + "\n", + "A = assemblage(coords, elems3)\n", + "b = second_membre(coords, elems3, f)\n", + "U_d = calcul_Ud(coords, dirichlet)\n", + "b -= np.dot(A, U_d)\n", + "\n", + "A_tild, b_tild, coords_tild = tildage(A, b, coords, dirichlet)\n", + "\n", + "x = np.linalg.solve(A_tild, b_tild)\n", + "x_untild = untildage(x, dirichlet, U_d)\n", + "\n", + "# show(coords_tild, x)\n", + "# print(coords.shape, x_untild.shape)\n", + "# print(coords, x_untild)\n", + "show(coords, x_untild)\n", + "show(coords, u_ex(coords[:, 0], coords[:, 1]))" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "**Partie II : maillage mixte et ajoût des conditions de Neumann**\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "e3 = np.array(\n", + " [[1, 2, 12], [2, 3, 12], [3, 4, 14], [4, 5, 14], [2, 15, 3], [3, 15, 4]]\n", + ").astype(int)\n", + "\n", + "e4 = np.array(\n", + " [\n", + " [0, 1, 12, 11],\n", + " [11, 12, 13, 10],\n", + " [12, 3, 14, 13],\n", + " [10, 13, 8, 9],\n", + " [13, 14, 7, 8],\n", + " [14, 5, 6, 7],\n", + " ]\n", + ").astype(int)\n", + "\n", + "dds = np.array([2, 15, 4, 6, 7, 8, 9, 10, 11, 0]).astype(int)\n", + "\n", + "nns = np.array([[4, 5], [5, 6], [0, 1], [1, 2]]).astype(int)\n", + "\n", + "ccs = np.array(\n", + " [\n", + " [0.0, 0.0],\n", + " [1 / 3, 0],\n", + " [0.53333333333333, 0.0],\n", + " [2 / 3, 1 / 3],\n", + " [1.0, 0.47],\n", + " [1, 2 / 3],\n", + " [1.0, 1.0],\n", + " [2 / 3, 1.0],\n", + " [1 / 3, 1.0],\n", + " [0.0, 1.0],\n", + " [0.0, 2 / 3],\n", + " [0.0, 1 / 3],\n", + " [1 / 3, 1 / 3],\n", + " [1 / 3, 2 / 3],\n", + " [2 / 3, 2 / 3],\n", + " [1.0, 0.0],\n", + " ]\n", + ")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "**Compléments d'analyse du système**\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3", + "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.2" + } + }, + "nbformat": 4, + "nbformat_minor": 4 +} diff --git a/test.py b/test.py new file mode 100644 index 0000000..4efa0f9 --- /dev/null +++ b/test.py @@ -0,0 +1,17 @@ +import numpy as np + +############################# Maillage mixte ################ +e3 = np.array([[1, 2, 12], [2, 3, 12], [3, 4, 14], [ + 4, 5, 14], [2, 15, 3], [3, 15, 4]]).astype(int) +e4 = np.array([[0, 1, 12, 11], [11, 12, 13, 10], [12, 3, 14, 13], [ + 10, 13, 8, 9], [13, 14, 7, 8], [14, 5, 6, 7]]).astype(int) +dds = np.array([2, 15, 4, 6, 7, 8, 9, 10, 11, 0]).astype(int) +nns = np.array([[4, 5], [5, 6], [0, 1], [1, 2]]).astype(int) +ccs = np.array([[0., 0.], [0.33333333333333, 0], [0.53333333333333, 0.], + [0.66666666666667, 0.33333333333333], [ + 1., 0.47], [1, 0.66666666666667], + [1., 1.], [0.66666666666667, 1.], [ + 0.33333333333333, 1.], [0., 1.], + [0., 0.66666666666667], [0., 0.33333333333333], [ + 0.33333333333333, 0.33333333333333], + [0.33333333333333, 0.66666666666667], [0.66666666666667, 0.66666666666667], [1., 0.]]) diff --git a/tp_edp_2021_2022.pdf b/tp_edp_2021_2022.pdf new file mode 100644 index 0000000..910b4af Binary files /dev/null and b/tp_edp_2021_2022.pdf differ