{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "%config InlineBackend.figure_format='retina'"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Verification and Validation (Fish and Belytschko, Section 8.2)\n",
    "\n",
    "Two of the most important things from any numerical simulation are:\n",
    "\n",
    "1. Verification: Are we solving the equations correctly.\n",
    "2. Validation: Are we solving the right equations for the problem that we are trying to model.\n",
    "\n",
    "__Validation__ is a significantly more involved issue, and cuts right to the core of how we try to model a physical phenomena. For example, if we are modeling a material then is the material elastic, viscous or plastic? Is the constitutive relation that we are writing correct? Are the material properties correct? Do the boundary conditions correctly represent the real problem we are modeling? This requires significant experience, connection with the actual experiments and finding parallels with other similar systems. We are not going to discuss this matter further. \n",
    "\n",
    "__Verification__ is a more well defined problem. We believe that certain equations that we have are **correct**. The more relatively modest goal is to ensure that numerical solution to the differential equation is the *correct* solution, i.e., the __strong form__ of the differential equation is solved correctly. However, since the finite element solution almost never gets the exact solution in the __strong__ sense, the problem of verification becomes a bit more involved. The most basic test that any correct FEA code should satisfy is called as the [__patch test__](https://en.wikipedia.org/wiki/Patch_test_(finite_elements)). \n",
    "\n",
    "The simple idea behind this test is as follows. If an analytical solution $\\theta$ could be exactly represent by the finite element solution $\\theta^h$, then even with a few number of arbitrary shaped and sized elements, one should __exactly__ get $\\theta = \\theta^h$. Since any _good_ finite element interpolation should at least be linear complete (linear elements, e.g., 3 node triangle), one can perform this test with respect to a __linear__ field. In the case of heat conduction problem that we are trying to solve, if the temperature field is of the form  \n",
    "$$T(x, y) = \\alpha_0 + \\alpha_1 x + \\alpha_2 y $$,\n",
    "where the coefficients $\\alpha_i$ are some arbitrary constants. Now this field $T(x,y)$ satisfies the heat conduction equation (which is a Poisson equation)\n",
    "$$\\nabla^2 T + s = 0$$,\n",
    "where $s$, the source term is zero. Hence, if we provide __essential__ boundary conditions $T(x, y) = \\bar{T}$ on $\\Gamma_T$, corresponding to this function, we have a boundary value problem whose solution is the same as that above. Since, the solution to a well-posed, linear, boundary value problem (BVP) is unique, this solution is __THE__ solution to the problem. Also, since, as mentioned earlier, any FEA formulation is linear complete, if we create a mesh like:\n",
    "<img src=\"Patch.png\" width=\"200\"/>\n",
    "On this particular mesh with arbitrary size elements, if we now now solve the BVP using FEA we should get an answer that should be extremely close to the exact solution (error $< 10^{-8}$). \n",
    "\n",
    "In FEniCS, since the quadrilateral element support is not so prominent as of now, we use triangular elements for the patch test. Using the following steps.  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "1. We create a mesh of our own as below. Please note the steps. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "# import the libraries\n",
    "\n",
    "import dolfin as df\n",
    "import matplotlib.pyplot as plt\n",
    "import mshr \n",
    "import numpy as np # numpy library for arrays etc.\n",
    "\n",
    "# define an editor to write the mesh\n",
    "%matplotlib inline\n",
    "mesh = df.Mesh()\n",
    "editor = df.MeshEditor()\n",
    "editor.open(mesh, 'triangle', 2, 2)\n",
    "editor.init_vertices(5)\n",
    "editor.init_cells(4)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "# add vertices\n",
    "editor.add_vertex(0, np.array([0.0, 0.0]))\n",
    "editor.add_vertex(1, np.array([1.0, 0.0]))\n",
    "editor.add_vertex(2, np.array([0.75, 0.25]))\n",
    "editor.add_vertex(3, np.array([1.0, 1.0]))\n",
    "editor.add_vertex(4, np.array([0.0, 1.0]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "# add cells\n",
    "editor.add_cell(0, np.array([0, 1, 2], dtype=np.uintp))\n",
    "editor.add_cell(1, np.array([2, 1, 3], dtype=np.uintp))\n",
    "editor.add_cell(2, np.array([2, 3, 4], dtype=np.uintp))\n",
    "editor.add_cell(3, np.array([2, 4, 0], dtype=np.uintp))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x113a5c9e8>,\n",
       " <matplotlib.lines.Line2D at 0x113a5ce10>]"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAf0AAAHwCAYAAACsUrZWAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAAWJQAAFiUBSVIk8AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAgAElEQVR4nO3de7CkdXng8e/DzSDMDMI4xhoMAxMuA25lFeQSNgokjkazCSa6f8QQllo1rtSCSkqJxgBxjZrEyGU1WZMoGK1UYqygJq5OBBYNl1gMm1QiyBDgIEYMBSgDw0Qm8Owf73vgTM/pc/p0v93v7fupOvU6/fb5za/bmfmefvpCZCaSJKn79qp7A5IkaTaMviRJPWH0JUnqCaMvSVJPGH1JknrC6EuS1BNGX5KknjD6kiT1hNGXJKknjL4kST1h9CVJ6gmjL0lST+xT9wZmKSLuAVYDczVvRZKkcW0Atmfm4Sv9xl5FH1i9//77H7xp06aD696IJEnjuP3229m5c+dY39u36M9t2rTp4K1bt9a9D0mSxnL88cdz6623zo3zvT6nL0lSTxh9SZJ6wuhLktQTRl+SpJ4w+pIk9YTRlySpJ4y+JEk9YfQlSeoJoy9JUk8YfUmSesLoS5LUE0ZfkqSeqCT6EfHaiLgiIr4WEdsjIiPiU2OudWhEfDwivhMRP4iIuYi4NCKeU8VeJUnqq6r+K3u/DvwY8BjwbeCYcRaJiI3AjcA64HPAN4ETgfOBV0bEqZn5UCU7liSpZ6oa778NOApYDfz3Cdb5KEXwz8vMMzPzwsw8A/gwcDTwvol3KklST1XySD8zr5v/3xEx1hoRcQSwGZgDPjJw+iLgTcBZEXFBZu4Yb6fVeeqpp7jnnns48MADWbt27di3W5LUbffeey/Pe97zeNaznsXee+9d616qGu9X4YzyuCUzn1p4IjMfjYgbKH4oOBm4ZqmFImLrkFNjPe2wmPe+971VLSVJ6olf+7VfY7/99qvt92/Sq/ePLo/bhpy/szweNYO9LOmJJ56oewuSpBbad999a/39m/RIf015fGTI+fnLD1puocw8frHLywnAi1e+td0t9VOaY35J0rzM3O3XdTeiSdFfzvw9lUtea0ZOOOEEbrnllt0uO/DAAzn77LNZu3ZtTbuSJDXFrl27+K3f+q2nf/1DP/RDNe6m0KTx/vwj+TVDzq8euF7jPPbYY1x11VU8+OCDdW9FklSzO++8c7dfP/e5z61pJ89oUvTvKI/DnrM/sjwOe86/Nscdd9zTz9MYfkkSwG233bbbr+se7UOzoj//tr/NEbHbviJiFXAqsBO4edYbW86P/MiP8PrXv97wS5KAYrS/bVvjHqPOPvoRsW9EHFN++t7TMvMuYAuwATh34NsuAQ4APtmE9+gv5rDDDjP8kiSgGO3v2rWr7m3soarP3j8zIq6MiCuBC8uLT5m/LCJ+d8HV1wO3s/h77d8CPABcHhFXR8T7I+Jaik/82wa8u4r9TovhlyTBnqP9pqjqkf5/BM4uv15RXnbEgsteO8oi5aP9E4ArgZOAC4CNwOXAKW343H3DL0n91tTRPlQU/cy8ODNjia8NC647N3jZwFr3ZeY5mfn8zNwvMw/LzPMz8+Eq9joLhl+S+qupo31o1gv5OsXwS1I/LRztb9iwob6NLMLoT5Hhl6R+GRztH3vssTXuZk9Gf8oMvyT1x8LR/tq1a1m3bl3NO9qd0Z8Bwy9J/bBwtN+0R/lg9GfG8EtStw2O9o877rgad7M4oz9Dhl+SumtwtN+Ez9ofZPRnzPBLUjcNjvab8Fn7g4x+DQy/JHVLG0b7YPRrY/glqTvaMNoHo18rwy9J3dCG0T4Y/doZfklqt7aM9sHoN4Lhl6T2astoH4x+Yxh+SWqntoz2weg3iuGXpHZp02gfjH7jGH5Jao82jfbB6DeS4ZekdmjTaB+MfmMZfklqtraN9sHoN5rhl6TmattoH4x+4xl+SWqmto32wei3guGXpGZp42gfjH5rGH5Jao42jvbB6LeK4ZekZmjjaB+MfusYfkmqV1tH+2D0W8nwS1J92jraB6PfWoZfkurR1tE+GP1WM/ySNFttHu2D0W89wy9Js9Pm0T4Y/U4w/JI0G20e7YPR7wzDL0nT1fbRPhj9TjH8kjQ9bR/tg9HvHMMvSdPR9tE+GP1OMvySVK0ujPbB6HeW4Zek6nRhtA9Gv9MMvyRVowujfTD6nWf4JWkyXRntg9HvBcMvSePrymgfjH5vGH5JGk9XRvtg9HvF8EvSynRptA9Gv3cMvySNrkujfTD6vWT4JWk0XRrtg9HvLcMvSUvr2mgfjH6vGX5JGq5ro30w+r1n+CVpcV0b7YPRF4ZfkgZ1cbQPRl8lwy9Jz+jiaB+MvhYw/JJU6OJoH4y+Bhh+SX3X1dE+GH0twvBL6rOujvbB6GsIwy+pr7o62gejryUYfkl90+XRPhh9LcPwS+qTLo/2wehrBIZfUl90ebQPRl8jMvySuq7ro30w+loBwy+py7o+2gejrxUy/JK6quujfTD6GoPhl9Q1fRjtg9HXmAy/pC7pw2gfjL4mYPgldUUfRvtg9DUhwy+p7foy2gejrwoYfklt1pfRPhh9VcTwS2qrvoz2weirQoZfUtv0abQPRl8VM/yS2qRPo30w+poCwy+pLfo02ocKox8Rh0bExyPiOxHxg4iYi4hLI+I5K1zn1RGxJSK+HRE7I+LuiPhMRJxS1V41fYZfUtP1bbQPFUU/IjYCW4FzgK8DHwbuBs4HboqIQ0Zc54PAXwEvBr4EXAbcCvwccENE/FIV+9VsGH5JTda30T5U90j/o8A64LzMPDMzL8zMMyjifzTwvuUWiIgfBn4V+Ffg2Mx8Q7nOa4FXAAH8ZkX71YwYfklN1bfRPlQQ/Yg4AtgMzAEfGTh9EbADOCsiDlhmqcPK/fxdZj6w8ERmXgc8CnT/x7AOMvySmqaPo32o5pH+GeVxS2Y+tfBEZj4K3AA8Gzh5mXXuBJ4AToyItQtPRMRLgVXAVyrYr2pg+CU1SR9H+1BN9I8uj9uGnL+zPB611CKZ+TDwTuB5wG0R8bGIeH9E/DmwBfgb4FdG2VBEbF3sCzhmlO/XdBh+SU3Rx9E+VBP9NeXxkSHn5y8/aLmFMvNS4OeBfYA3AhcCrwPuA64cHPurfQy/pLr1dbQPs3mf/vyPT7nsFSPeAfwFcCWwETgAOJ7inQCfjojfHuU3zMzjF/sCvjnODVC1DL+kOvV1tA/VRH/+kfyaIedXD1xvURFxGvBB4POZ+fbMvDszH8/MW4HXAP8CXFC+cFAtZ/gl1aWvo32oJvp3lMdhz9kfWR6HPec/72fK43WDJzLzcYr3/+8FvGilG1QzGX5Js9bn0T5UE/35SG+OiN3Wi4hVwKnATuDmZdZ5VnkcNmeZv/yJcTapZjL8kmapz6N9qCD6mXkXxavrNwDnDpy+hOJ5+U9m5g6AiNg3Io4pP8Vvoa+VxzdFxPqFJyLipyl+ePg34MZJ96xmMfySZqXPo32o7oV8bwEeAC6PiKvLt9pdC7yNYqz/7gXXXQ/cDlwzsMZfULwP/3nA7RFxVUR8MCI+D/w1xQsCL8zMhyrasxrE8Euatr6P9qGi6JeP9k+geNX9ScAFFK++vxw4ZZRQlx/s8yqKHxRuo3jx3gUUH+rzReAVmXlZFftVMxl+SdPU99E+VPiWvcy8LzPPycznZ+Z+mXlYZp5ffujOwuvNZWZk5oZF1tiVmZdm5smZuToz98nMdZn5M5m5paq9qrkMv6Rp6ftoH2bzPn1pRQy/pKo52i8YfTWS4ZdUJUf7BaOvxjL8kqriaL9g9NVohl/SpBztP8Poq/EMv6RJONp/htFXKxh+SeNytP8Mo6/WMPySVsrR/u6MvlrF8EtaCUf7uzP6ah3DL2lUjvZ3Z/TVSoZf0nIc7e/J6Ku1DL+kpTja35PRV6sZfknDONrfk9FX6xl+SYMc7S/O6KsTDL+khRztL87oqzMMv6R5jvYXZ/TVKYZfkqP94Yy+OsfwS/3maH84o69OMvxSfznaH87oq7MMv9Q/jvaXZvTVaYZf6hdH+0sz+uo8wy/1h6P9pRl99YLhl7rP0f7yjL56w/BL3eZof3lGX71i+KXucrS/PKOv3jH8Uvc42h+N0VcvGX6pWxztj8boq7cMv9QdjvZHY/TVa4Zfaj9H+6Mz+uo9wy+1m6P90Rl9CcMvtZmj/dEZfalk+KX2cbS/MkZfWsDwS+3iaH9ljL40wPBL7eFof2WMvrQIwy81n6P9lTP60hCGX2o2R/srZ/SlJRh+qbkc7a+c0ZeWYfil5nG0Px6jL43A8EvN4mh/PEZfGpHhl5rD0f54jL60AoZfqp+j/fEZfWmFDL9UL0f74zP60hgMv1QfR/vjM/rSmAy/NHuO9idj9KUJGH5pthztT8boSxMy/NLsONqfjNGXKmD4pelztD85oy9VxPBL0+Vof3JGX6qQ4Zemx9H+5Iy+VDHDL1XP0X41jL40BYZfqpaj/WoYfWlKDL9UHUf71TD60hQZfmlyjvarY/SlKTP80mQc7VfH6EszYPil8Tnar47Rl2bE8Esr52i/WkZfmiHDL62Mo/1qGX1pxgy/NDpH+9Uy+lINDL+0PEf71TP6Uk0Mv7Q0R/vVM/pSjQy/NJyj/eoZfalmhl/ak6P96TD6UgMYfml3jvanw+hLDWH4pWc42p8Ooy81iOGXHO1PU2XRj4hDI+LjEfGdiPhBRMxFxKUR8Zwx1vqJiPhsRNxfrnV/RGyJiFdVtV+pqQy/+s7R/vRUEv2I2AhsBc4Bvg58GLgbOB+4KSIOWcFavw58FXgp8CXgQ8AXgOcAp1WxX6npDL/6zNH+9FT1SP+jwDrgvMw8MzMvzMwzKOJ/NPC+URaJiNcB7wW+AhyRmedk5rsy802Z+RLg3RXtV2o8w68+crQ/XRNHPyKOADYDc8BHBk5fBOwAzoqIA5ZZZy/gg8DjwC9m5qOD18nMXZPuV2oTw6++cbQ/XVU80j+jPG7JzKcWnijDfQPwbODkZdb5ceBw4IvA9yLi1RHxzog4PyJOqWCfUisZfvWJo/3pqiL6R5fHbUPO31kej1pmnZeUx38FbgX+CvgAcClwY0RcHxEj/cgXEVsX+wKOGeX7paYx/OoDR/vTV0X015THR4acn7/8oGXWWVce3wzsD/wUsAp4IfBlihf2fWb8bUrtZvjVdY72p28W79Ofn83kMtfbe8H1X5uZ12TmY5n5DeA1wLeBl40y6s/M4xf7Ar457o2QmsDwq8sc7U9fFdGffyS/Zsj51QPXG+Z75fHuzPyHhScycyfFo32AE1e8Q6lDDL+6yNH+bFQR/TvK47Dn7I8sj8Oe8x9c5/tDzs//ULD/iPuSOsvwq2sc7c9GFdG/rjxuLt9297SIWAWcCuwEbl5mna8C/w4cGRH7LXL+heVxbvytSt1h+NUljvZnY+LoZ+ZdwBZgA3DuwOlLgAOAT2bmDoCI2Dcijik/xW/hOg8Cf0bxNMFvLDwXES8HXkHxFMGXJt2z1BWGX13gaH92qnoh31uAB4DLI+LqiHh/RFwLvI1irL/wk/TWA7cD1yyyztuBfwbeHRFfjYjfjYjPAP8HeBJ4Y2YOG/9LvWT41XaO9menkuiXj/ZPAK4ETgIuADYClwOnZOZDI67zQPn9HwZeAJxH8eE/fw38RGb6lj1pEYZfbeZof3Yqe8teZt5Xflb+8zNzv8w8LDPPz8yHB643l5mRmRuGrPNwZr49Mw8v1zkkM38uM5d7TYDUa4ZfbeRof7Zm8T59STNi+NU2jvZny+hLHWP41SaO9mfL6EsdZPjVBo72Z8/oSx1l+NV0jvZnz+hLHWb41WSO9mfP6EsdZ/jVRI7262H0pR4w/GoaR/v1MPpSTxh+NYmj/XoYfalHDL+awNF+fYy+1DOGX3VztF8foy/1kOFXnRzt18foSz1l+FUHR/v1MvpSjxl+zZqj/XoZfannDL9mydF+vYy+JMOvmXC0Xz+jLwkw/Jo+R/v1M/qSnmb4NU2O9utn9CXtxvBrGhztN4PRl7QHw6+qOdpvBqMvaVGGX1VytN8MRl/SUIZfVXC03xxGX9KSDL8m5Wi/OYy+pGUZfk3C0X5zGH1JIzH8Goej/WYx+pJGZvi1Uo72m8XoS1oRw6+VcLTfLEZf0ooZfo3C0X7zGH1JYzH8Wo6j/eYx+pLGZvi1FEf7zWP0JU3E8GsxjvabyehLmpjh1yBH+81k9CVVwvBrIUf7zWT0JVXG8Asc7TeZ0ZdUKcMvR/vNZfQlVc7w95uj/eYy+pKmwvD3k6P9ZjP6kqbG8PePo/1mM/qSpsrw94uj/WYz+pKmzvD3g6P95jP6kmbC8Hefo/3mM/qSZsbwd5uj/eYz+pJmyvB3k6P9djD6kmbO8HePo/12MPqSamH4u8XRfjsYfUm1Mfzd4Gi/PYy+pFoZ/vZztN8eRl9S7Qx/uznabw+jL6kRDH87OdpvF6MvqTEMf/s42m8Xoy+pUQx/uzjabxejL6lxDH87ONpvH6MvqZEMf/M52m8foy+psQx/sznabx+jL6nRDH8zOdpvJ6MvqfEMf/M42m8noy+pFQx/szjabyejL6k1DH8zONpvL6MvqVUMf/0c7beX0ZfUOoa/Xo7228voS2olw18PR/vtZvQltZbhnz1H++1m9CW1muGfLUf77Wb0JbWe4Z8NR/vtZ/QldYLhnz5H++1XWfQj4tCI+HhEfCcifhARcxFxaUQ8Z4I1z4qILL/eUNVeJXWT4Z8uR/vtV0n0I2IjsBU4B/g68GHgbuB84KaIOGSMNV8AXAE8VsUeJfWD4Z8OR/vdUNUj/Y8C64DzMvPMzLwwM8+giP/RwPtWslgUPz5+AngI+IOK9iipJwx/9Rztd8PE0Y+II4DNwBzwkYHTFwE7gLMi4oAVLHsecAbF5GDHpHuU1D+Gv1qO9ruhikf6Z5THLZn51MITmfkocAPwbODkURaLiE3AB4DLMvOrFexPUk8Z/mo42u+OKqJ/dHncNuT8neXxqOUWioh9gD8BvgW8a9wNRcTWxb6AY8ZdU1I7Gf7JOdrvjiqiv6Y8PjLk/PzlB42w1m8ALwL+a2bunHRjkgSGf1KO9rtjFu/Tn//TkUteKeJEikf3H8rMmyb5DTPz+MW+gG9Osq6k9jL843G03y1VRH/+kfyaIedXD1xvDwvG+tuA91SwJ0nag+FfOUf73VJF9O8oj8Oesz+yPA57zh/gwPL7NwH/tuADeZLiHQAAf1hedunEO5bUW4Z/ZRztd0sV0b+uPG6OiN3Wi4hVwKnATuDmJdb4AfDHQ77+X3mdvy1/PdHoX5IM/2gc7XfPxNHPzLuALcAG4NyB05cABwCfzMwdABGxb0QcU36K3/waOzPzDYt9AZ8vr3ZVedmfTbpnSTL8y3O03z1VvZDvLcADwOURcXVEvD8irgXeRjHWf/eC664Hbgeuqej3lqSxGP6lOdrvnkqiXz7aPwG4EjgJuADYCFwOnJKZD1Xx+0hS1Qz/4hztd1Nlb9nLzPsy85zMfH5m7peZh2Xm+Zn58MD15jIzMnPDiOteXF7/j6raqyQtZPj35Gi/m2bxPn1JajzDvztH+91k9CWpZPgLjva7y+hL0gKG39F+lxl9SRrQ9/A72u8uoy9Ji+hr+B3td5vRl6Qh+hh+R/vdZvQlaQl9C7+j/W4z+pK0jL6E39F+9xl9SRpBH8LvaL/7jL4kjajr4Xe0331GX5JWoKvhd7TfD0Zfklaoi+F3tN8PRl+SxtC18Dva7wejL0lj6kr4He33h9GXpAl0IfyO9vvD6EvShNoefkf7/WH0JakCbQ2/o/1+MfqSVJE2ht/Rfr8YfUmqUNvC72i/X4y+JFWsLeF3tN8/Rl+SpqAN4Xe03z9GX5KmpOnhd7TfP0ZfkqaoqeF3tN9PRl+SpqyJ4Xe0309GX5JmoGnhd7TfT0ZfkmakKeF3tN9fRl+SZqgJ4Xe0319GX5JmrO7wO9rvL6MvSTWoK/yO9vvN6EtSTeoIv6P9fjP6klSjWYff0X6/GX1Jqtmswu9oX0ZfkhpgFuF3tC+jL0kNMe3wO9qX0ZekBplW+B3tC4y+JDXONMLvaF9g9CWpkaoOv6N9gdGXpMaqKvyO9jXP6EtSg1URfkf7mmf0JanhJg2/o33NM/qS1ALjht/RvhYy+pLUEuOE39G+FjL6ktQiKw2/o30tZPQlqWVGDb+jfQ0y+pLUQqOE39G+Bhl9SWqp5cLvaF+DjL4ktdiw8N9///2O9rUHoy9JLbdY+D/2sY852tcejL4kdcBg+BdytK95Rl+SOuDJJ59k165drFq1ao9z++yzTw07UhP5J0GSWurJJ5/knnvu4Rvf+AZ33HEHO3fuXPR61157LZs2bWLt2rUz3qGaxuhLUouMGvpBV111FWeffbbh7zmjL0kNN2roV69ezbHHHstxxx3H+vXrmZub40//9E/ZtWvX06/qN/z9ZvQlqYHGDf3CF+wdfvjhvP71r+fTn/604Rdg9CWpMaoI/aD5V/UbfoHRl6RaTSP0gwy/5hl9SZqxWYR+kOEXGH1Jmok6Qj/I8MvoS9KUNCH0gwx/vxl9SapQE0M/yPD3l9GXpAm1IfSDDH8/GX1JGkMbQz/I8PeP0ZekEXUh9IMMf78YfUlaQhdDP8jw90dl0Y+IQ4HfBF4JHALcD1wNXJKZ3xvh+w8BXgO8GvgPwHrgCeAfgU8An8jMp6raryQN04fQDzL8/VBJ9CNiI3AjsA74HPBN4ETgfOCVEXFqZj60zDKvA36f4oeF64BvAc8Dfh74I+CnI+J1mZlV7FmSFupj6AcZ/u6r6pH+RymCf15mXjF/YUT8HvA24H3Am5dZYxvws8BfL3xEHxHvAr4O/ALFDwCfrWjPknrO0O/J8HfbxNGPiCOAzcAc8JGB0xcBbwLOiogLMnPHsHUy89ohl383Iv6A4geH0zD6kiZg6Jdn+Lurikf6Z5THLYPPuWfmoxFxA8UPBScD14z5e+wqj/8+5vdL6jFDv3KGv5uqiP7R5XHbkPN3UkT/KMaIfkTsA/xy+csvjfg9W4ecOmalv7+kdjL0kzP83VNF9NeUx0eGnJ+//KAx1/8A8ELgi5n55THXkNQDhr56hr9bZvE+/fm/TSt+1X1EnAdcQPFugLNG/b7MPH7IeluBF690H5Kay9BPn+HvjiqiP/9Ifs2Q86sHrjeSiDgXuAy4DfjJzHx4vO1J6hpDP3uGvxuqiP4d5fGoIeePLI/DnvPfQ0S8Ffgw8E8UwX9g/O1J6gJDXz/D335VRP+68rg5IvYaeI/9KuBUYCdw8yiLRcQ7KZ7H/3vg5Zn5YAV7lNRChr55DH+7TRz9zLwrIrZQvEL/XOCKBacvAQ4A/vf8e/QjYl9gI7ArM+9auFZEvIfio3y3Apsd6Uv9Y+ibz/C3V1Uv5HsLxcfwXh4RPwncDpwEnE4x1n/3guuuL8/fC2yYvzAizqYI/pPA14DzFvlLPJeZV1a0Z0kNYejbx/C3UyXRLx/tn8Az/8GdV1F8hv7lFP/BnVEesR9eHvcG3jrkOtcDV062W0lNYOjbz/C3T2Vv2cvM+4BzRrjeHM+8jW/h5RcDF1e1H0nNY+i7x/C3yyzepy+pxwx99xn+9jD6kipn6PvH8LeD0ZdUCUMvw998Rl/S2Ay9Bhn+ZjP6klbE0Gs5hr+5jL6kZRl6rZThbyajL2lRhl6TMvzNY/QlPc3Qq2qGv1mMvtRzhl7TZvibw+hLPWToNWuGvxmMvtQThl51M/z1M/pShxl6NY3hr5fRlzrG0KvpDH99jL7UAYZebWP462H0pZYy9Go7wz97Rl9qEUOvrjH8s2X0pYYz9Oo6wz87Rl9qIEOvvjH8s2H0pYYw9Oo7wz99Rl+qkaGXdmf4p8voSzNm6KWlGf7pMfrSDBh6aWUM/3QYfWlKDL00GcNfPaMvVcjQS9Uy/NUy+tKEDL00XYa/OkZfGoOhl2bL8FfD6EsjMvRSvQz/5Iy+tARDLzWL4Z+M0ZcGGHqp2Qz/+Iy+hKGX2sbwj8foq7cMvdRuhn/ljL56xdBL3WL4V8boq/MMvdRthn90Rl+dZOilfjH8ozH66gxDL/Wb4V+e0VerGXpJCxn+pRl9tY6hl7QUwz+c0VcrGHpJK2H4F2f01ViGXtIkDP+ejL4axdBLqpLh353RV+0MvaRpMvzPMPqqhaGXNEuGv2D0NTOGXlKdDL/R15QZeklN0vfwG31VztBLarI+h9/oqxKGXlKb9DX8Rl9jM/SS2qyP4Tf6WhFDL6lL+hZ+o69lGXpJXdan8Bt9LcrQS+qTvoTf6Otphl5Sn/Uh/Ea/5wy9JD2j6+E3+j1k6CVpuC6H3+j3hKGXpNF1NfxGv8MMvSSNr4vhN/odY+glqTpdC7/R7wBDL0nT06XwG/2WMvSSNDtdCb/RbxFDL0n16UL4jX7DGXpJao62h9/oN5Chl6TmanP4jX5DGHpJao+2ht/o18jQS1J7tTH8Rn/GDL0kdUfbwl9Z9CPiUOA3gVcChwD3A1cDl2Tm92a9TpMYeknqrqXC3zSVRD8iNgI3AuuAzwHfBE4EzgdeGRGnZuZDs1qnCQy9JPXHsPCffvrpdW9tN1U90v8oRajPy8wr5i+MiN8D3ga8D3jzDNephaGXpP5aLPxf+MIX6t7WbiaOfkQcAWwG5oCPDJy+CHgTcFZEXJCZO6a9Th3uuusu7r//fkMvST03GP6mqeKR/hnlcUtmPrXwRGY+GhE3UMT8ZOCaGawzE7fccsvT/3vbtm1LXnf16tWsX7+e7du3c9NNN017a5Kkmu233357RP9b3/pWTbt5RhXRP7o8DivfnRSxPoqlY13VOkTE1iGnjlnq+6Zl+/btbN++vY7fWpLUIN///vc56KCDavv996pgjTXl8ZEh5+cvX+5WVrWOJEmNtGrVqlp//1m8T3/+yeuc1TqZefyiCxQTgBdPuA8AXvayl3H99dezevVqTjjhBA4++OAqlpUkdcjjjz/O9ddfz44dOzj11FPZe++9a91PFdGffwS+Zsj51QPXm/Y6M3Haaadx2mmn1b0NSXl9y34AAAh3SURBVFLDveQlL6l7C0+rYrx/R3k8asj5I8vj0q92q24dSZK0iCqif1153BwRu60XEauAU4GdwM0zWkeSJC1i4uhn5l3AFmADcO7A6UuAA4BPzr+3PiL2jYhjyk/fG3sdSZK0MlW9kO8tFB+fe3lE/CRwO3AScDrFOP7dC667vjx/L0Xgx11HkiStQBXj/flH6ScAV1JE+gJgI3A5cMqon5df1TqSJGlPlb1lLzPvA84Z4XpzPPP2u7HXkSRJK1PJI31JktR8Rl+SpJ4w+pIk9YTRlySpJ4y+JEk9YfQlSeoJoy9JUk8YfUmSesLoS5LUE0ZfkqSeiMysew8zExEP7b///gdv2rSp7q1IkjSW22+/nZ07dz6cmYes9Hv7Fv17gNXAXEVLHlMev1nRen3ifTc+77vxed+Nz/tufFXfdxuA7Zl5+Eq/sVfRr1pEbAXIzOPr3kvbeN+Nz/tufN534/O+G1+T7juf05ckqSeMviRJPWH0JUnqCaMvSVJPGH1JknrCV+9LktQTPtKXJKknjL4kST1h9CVJ6gmjL0lSTxh9SZJ6wuhLktQTRl+SpJ4w+gtExKER8fGI+E5E/CAi5iLi0oh4Th3rtMmktzkiDomIN0TEX0bEP0fEzoh4JCL+NiL+W0R0+s/qNP7MRMRZEZHl1xuq3G+TVHnfRcRPRMRnI+L+cq37I2JLRLxqGnuvW4X/5r26vJ++Xf7dvTsiPhMRp0xr73WKiNdGxBUR8bWI2F7+HfvUmGvNtBd+OE8pIjYCNwLrgM9R/HePTwROB+4ATs3Mh2a1TptUcZsj4s3A7wP3A9cB3wKeB/w8sAb4LPC67OAf2Gn8mYmIFwD/COwNHAi8MTP/qMp9N0GV911E/DrwXuBB4K8o/iyuBV4EXJeZ76j8BtSown/zPgi8A3gIuJri/vtR4GeBfYBfzsyxgthUEfH3wI8BjwHfBo4BPp2Zv7TCdWbfi8z0q+jIl4EE/sfA5b9XXv4Hs1ynTV9V3GbgDOA/A3sNXP7DFD8AJPALdd/Wpt5/A98XwFeAu4DfKdd4Q923s8n3HfC68vp/A6xa5Py+dd/WJt535d/PJ4HvAusGzp1ernN33bd1Cvfd6cCR5d+108rb+ak6/j9Y8e9Z953XhC/giPIOvmeR6Kyi+GluB3DALNZp09csbjPwrvL3uKLu29uG+w84H3gKeClwcVejX+Hf272Au8vrPrfu29Wy++6kcp3PDTm/HXi07ts75ftyrOjX1YtOP0+6AmeUxy2Z+dTCE5n5KHAD8Gzg5Bmt0yazuM27yuO/T7BGU1V6/0XEJuADwGWZ+dUqN9pAVd13Pw4cDnwR+F75/PQ7I+L8rj4nTXX33Z3AE8CJEbF24YmIeClFvL5SyY67p5ZeGP3C0eVx25Dzd5bHo2a0TptM9TZHxD7AL5e//NI4azRcZfdfeV/9CcXTIe+afGuNV9V995Ly+K/ArRTP538AuBS4MSKuj4jnTrLRBqrkvsvMh4F3Urz+5raI+FhEvD8i/hzYQvF0ya9UsN8uqqUX+1S5WIutKY+PDDk/f/lBM1qnTaZ9mz8AvBD4YmZ+ecw1mqzK++83KF509p8yc+ekG2uBqu67deXxzRSj1p8C/g44DPgQ8ArgMxRj3K6o7M9dZl4aEXPAx4E3Ljj1z8CVmfnAuJvsuFp64SP90UR5nPSV41Wt0yZj3+aIOA+4gOIVrWdVuakWGen+i4gTKR7dfygzb5r6rtph1D97ey+4/msz85rMfCwzvwG8huLV2S/r8Kh/MSP/vY2IdwB/AVwJbAQOAI6neJ3EpyPit6e0x66bSi+MfmH+J6o1Q86vHrjetNdpk6nc5og4F7gMuA04vRwjdtHE99+Csf424D3Vba3xqvqz973yeHdm/sPCE+XEZH7CdOKKd9hcldx3EXEa8EHg85n59sy8OzMfz8xbKX5g+hfggog4ooI9d00tvTD6hTvK47DnTo4sj8Oee6l6nTap/DZHxFuB/wX8E0Xwvzv+9hqvivvvwPL7NwH/tuADeRK4qLzOH5aXXTrxjpuj6r+33x9yfv6Hgv1H3FcbVHXf/Ux5vG7wRGY+DnydojMvWukGe6CWXvicfmH+D+zmiNhr4SspI2IVcCqwE7h5Ruu0SaW3OSLeSfE8/t8DL8/MByveb9NUcf/9APjjIedeTPEP7t9S/CPTpdF/VX/2vkrxzpAjI2K/zHxi4PwLy+Pc5FtujKruu2eVx2EvdJy/fPA+VV29qPs9jk35YgUfkgDsS/EJTBsnWacrXxXed+8pr38LcHDdt6tt99+QtS+mo+/Tr/K+Az5VXv9/Dlz+corPPPg+cFDdt7dp9x3wX8rrfhdYP3Dup8v7bidwSN23d4r342ks8T79pvXCj+EtLfJxiLdTfPDE6RTjlR/P8uMQI2IDxat8783MDeOu0xVV3HcRcTbFC4GeBK5g8eex5jLzyuncivpU9WdvyNoXU4z4+/IxvOP+vV1H8b7oHwW+RjGWPozieekEfjEzPzP1GzRDFf293YsiXD8FPAr8JcUPAJsoRv8BvDUzL5vFbZqViDgTOLP85Q9TvMPjboo/OwAPZuavltfdQJN6UfdPSU36Al4AfILiM7efAO6leDHZwQPX20DxD8HcJOt06WvS+45nHpEu9fV/676dTb3/llh3/n7t5CP9Ku874GCKR1j3lOs8RPEP8cl138Ym33cUj2TfSjGG3k7xVMkDFJ93sLnu2zil+225f6/mFly3Ub3wkb4kST3hq/clSeoJoy9JUk8YfUmSesLoS5LUE0ZfkqSeMPqSJPWE0ZckqSeMviRJPWH0JUnqCaMvSVJPGH1JknrC6EuS1BNGX5KknjD6kiT1hNGXJKknjL4kST1h9CVJ6on/D1qMzmGalg8LAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "image/png": {
       "height": 248,
       "width": 254
      },
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# final formalities\n",
    "editor.close() # close the editor\n",
    "mesh.order # order the mesh\n",
    "df.plot(mesh) # plot"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Define Laplace equation. The exact solution is:\n",
    "\n",
    "$$\\theta = \\alpha_0 + \\alpha_1 x + \\alpha_2 y$$,\n",
    "\n",
    "we define $\\alpha_0 = 1, \\alpha_1 = 2, \\alpha_2 = 3$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "L2_error is =  3.624438086253896e-16\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAmIAAAH4CAYAAADpQ4FeAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAAWJQAAFiUBSVIk8AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAgAElEQVR4nOy9e7SlVXmn+7x1QbDkokTkBG0RgpStOagYlZgYgTRJTMeDF3qcPpEYEpPYOgIqjGi8tIU9HAdPYgQ0RhMvRHSMxJi200m8YBKGN7QdoxLbNoIgWERbCQKCgEUVVXueP+ac5WbXvq71rfV7916/Z4w9vqq11v6+t769a+9n/eY754xSCsYYY4wxZvpsUhdgjDHGGDOrWMSMMcYYY0RYxIwxxhhjRFjEjDHGGGNEWMSMMcYYY0RYxIwxxhhjRFjEjDHGGGNEWMSMMcYYY0RYxIwxxhhjRFjEjDHGGGNEWMSMMcYYY0RYxIwxxhhjRFjEjDHGGGNEWMSMMcYYY0QMImIR8YKIeFtEfCYivh8RJSI+MOK5HhkR742Ib0fEnojYFRGXRsRDh6jVGGOMMSYLWwY6z+uAU4B7gG8B20c5SUScCFwDHAP8FXAd8FTgAuDnI+IZpZTbB6nYGGOMMUbMUEOTrwAeCxwB/KcxzvMOqoSdX0o5u5Ty6lLKGcBbgZOBN41dqTHGGGNMEqKUMuwJI54FXA18sJTywjV83gnAjcAu4MRSyty85w4HvgMEcEwp5d4hazbGGGOMUZCpWf+MdrxqvoQBlFLuBj4HPBh4+rQLM8YYY4yZBJlE7OR2vH6J529ox8dOoRZjjDHGmIkzVLP+EBzZjnct8Xx//KiVThQRO5d46gnUCQW71lSZMcYYk4fjge+XUh6jLCIiPsiIk/NW4LpSyi9P4LwpySRiKxHtOE5T2+Y4ZOvDth738Iet5sWbYtj+uVHYtElbw6aYW/lFE2Zzgq/DZvF9UF8fYIu4hi2xX3r9WoP+67BZfB+2kuEeaH8mbBEPJn3thvvZfZ/+5yKw/bBD48knn7R1sBMm+rdNjUwi1hOvI5d4/ogFr1uSUsqpiz0eETu3HvfwJz/qkrVN7DzkkH1rev3QbHvQXun1AbYdoq3hyAfdJ70+wFGH7FaXwFFbfyC9/tFb9fNkHrblHun1H77lbun1AR62WXsPAI4R1/CwzdqfywA/sulQyXV/8ue+zZf+195dkosv4OSTtnLNJ350sPO1f9tg51sPZOoR+1o7LtUDdlI7LtVDNjH27tX66r17DuHePYdoa9irvf5dew7lrj2aH3qdO/cexp17D9PWcP+Dpde//f5t3H7/NmkNd+x7iPT63913ON/dd7i0hjv2P4Q79mvvw63i69+xfwt37Nf+bL5tTv8G0ax/MonY1e14VkQ8oK62fMUzgN3AF6ZdGFQZyyBk0uvvPSSFkKnJIGMZhEzJHfsekkLI1GSQsQxCpuS2ufssZGYspi5iEbE1Ira3VfQPUEq5EbiK2oT4sgWfdjGwDXi/eg2xDDKWQciUOB1rNSSQsQxCpsTpWCWDjGUQMmNGYZDv3Ig4Gzi7/fXYdjwtIq5of76tlHJR+/NxwLXAzVTpms9LqVscXR4RZ7bXPQ04nTok+doh6h2XLmPK3rF79xwi7R3rMqbsHbtrz6Hy3rE79x4m7R3rMqbsHbv9/m3S3rEuY8rese/uO1zeO3bH/odIe8e6jCl7x+7Yv0XaO9ZlTNU7ZtYnQ72FeCLwogWPndA+oErXRaxAKeXGiHgK8Ebg54FnU1fUvxy4uJRyx0D1DsLevVvkMgbaZv579x4ilzHQNvP3ZEwtZGoZA20z/x37HiKXMdA28/dkTC1kahkDbTP/bXP3WcbMqhlkaLKUsqOUEst8HD/vtbsWPrbgXN8spZxXSvk/SimHlFIeXUq5IJuEddw75t6xToahygzDlUrcO1bJMFSZYbhSiXvHzGrJ1Ky/rskgYxmETIl7x1oNCWQsg5Apce9YJYOMZRAyY5bDIjYgTsecjnUyyFgGIVPidKySQcYyCJkSp2NmOSxiEyCDjGUQMiVOx1oNCWQsg5ApcTpWySBjGYTMmIVYxCaE0zGnY50MMpZByJQ4HatkkLEMQqbE6ZhZiEVswmSQsQxCpsTpWKshgYxlEDIlTscqGWQsg5AZA7n2mtyweN0xrzvW8bpjXncMvO4YeN0x2Bjrju1jblCp3JdgU/lp40RsijgdczoGTsfA6Rg4Hes4HXM6NuvMnoiJZdu9Y+4d62SQsQxCpsS9Y5UMMpZByJS4d2x2mT0RA/bdpx+RzSBjGYRMidOxVkMCGcsgZEqcjlUyyFgGITOzxUyKGFQZUwuZ0zGnY50MMpZByJQ4HatkkLEMQqbE6dhsMbMi1lHLGDgdA6dj4HQMnI6B07FOBhnLIGRm4zPzIgZOxzoZZCyDkKnJIGMZhEyJ07FKBhnLIGRKnI6NRkTsioiyxMctqzzHry5zjv6xf9xa9XFQIvbdt4Uth+qmMkMVMvUyF4B8qQv1MheAdKmLLmPqpS7Uy1wA8qUu1MtcANKlLrqMqZe6UC9zAciXuljPy1yIuAu4dJHHV/vN9CXg4iWe+2ngDOBjI9T1ACxiC+jJmFLIvO6Y1x3reN0xrzsGXncMvO4YeKhyBO4spewY9ZNLKV+iythBRMTn2x//eNTzdzw0uQTqoUpw7xi4dwzcOwbuHQP3jnUyDFWqhytncdHTTETEE4CnA/8b+Ntxz2cRWwb3jlUyyFgGIVOTQcYyCJkS945VMshYBiEz6XlQRLwwIl4TERdExOkRsXmA8/5WO76nlOIesWng3jH3joF7xw7U4N4x947h3jHI0Tu2QdkeETsXe6KUcuoaznMscOWCx74REeeVUj41SmERcRjwQury8O8e5RwLcSK2SpyOVZyOOR0Dp2PgdKzjdMzpWFLeB5xJlbFtwI8D7wKOBz4WEaeMeN7/ABwFfKyU8s0B6nQitlacjjkdA6djB2pwOuZ0DKdjMLvp2P4Sg4ro/hIA160x+TqIUsrC2Y5fAV4SEfcAFwI7gOeOcOrfbMd3jV7dA3EiNgJOxypOx5yOgdMxcDrWcTrmdGwd8M52fOZaPzEi/i3wk8C3gI8OVZBFbAzUMgaeWQmeWQmeWQmeWQmeWdnJIGMWsrTc2o6j/LAYtEm/YxEbE6djlQwylkHI1GSQsQxCpsTpWCWDjGUQMpOO09rxprV8UkQcCpxLbdJ/z5AFWcQGQi1j4HQMnI6B0zFwOgZOxzoZZMxCNl0i4vER8bBFHn808Pb21w/Me3xrRGyPiBOXOe05wEOBjw7VpN+xiA2I07FKBhnLIGRqMshYBiFT4nSskkHGMgiZmRrnAN+OiI9FxDsi4s0R8WHgOuDHqP1dvz/v9ccB1wJ/v8w5e5P+2CvpL8TfGRPAMys9sxI8s/JADZ5Z6ZmVeGYlzO7MSgFXAycDT6IORW4D7gQ+S11X7MpSSlntySLiccBPMXCTfsciNiG8Z2XFe1Z6z0rwnpXgPSs73rNSv2flRqct1rrqBVtLKbuAWOb5a5d7flw8NDlh1EOV4N4xcO8YuHcM3DsG7h3rZBiq9HClAYvYVHDvWCWDjGUQMjUZZCyDkClx71glg4xlEDIz21jEpohaxsDpGDgdA6dj4HQMnI51MsiYhWx2mTkRK2Viw7yrwulYJYOMZRAyNRlkLIOQKXE6VskgYxmEzMweMydiAGXvZsrezdIa1DIGTsfA6Rg4HQOnY+B0rJNBxixks8VMf7XL3s3EIYPtUrBmPLOy4pmVnlkJnlkJnlnZ8czK9TOz8n42DSqv9/Pdwc61XpjJRGw+TscqTsecjoHTMXA6Bk7HOk7HzDSYeRHrZJAxtZC5d8y9Y50MMpZByJS4d6ySQcYyCJnZuFjE5uF0rJJBxjIImRKnY62GBDKWQciUOB2rZJAxC9nGxCK2CBlkTC1kTsecjnUyyFgGIVPidKySQcYyCJnZWFjElsDpWCWDjGUQMiVOx1oNCWQsg5ApcTpWUcvYfvEyTGZYLGIrkEHG1ELmdMzpWCeDjGUQMiVOxyoZZEwtZGZjYBFbBU7HKhlkLIOQKXE61mpIIGMZhEyJ07GKZcyMi0VsDWSQMbWQOR1zOtbJIGMZhEyJ07FKBhmzkJlRsYitEadjlQwylkHIlDgdazUkkLEMQqbE6VjFMmZGwSI2IhlkTC1kTsecjnUyyFgGIVPidKySQcYsZGYtWMTGwOlYJYOMZRAyJU7HWg0JZCyDkClxOlaxjJnVov8tvgHwnpXesxK8Z2XHe1Z6z0rwnpWQY8/KSbO/bB5UevcXbbihwInYQDgdqzgdczoGTsfA6Rg4Hes4HTPLYREbmAwyphYy9465d6yTQcYyCJkS945VMsiYhcwshkVsAjgdq2SQsQxCpsTpWKshgYxlEDIlTscqljGzEIvYBMkgY2ohczrmdKyTQcYyCJkSp2OVDDJmITMdi9iEcTpWySBjGYRMidOxVkMCGcsgZEqcjlUsYwYsYlMjg4yphczpmNOxTgYZyyBkSpyOVTLImIVstrGITRGnY5UMMpZByJQ4HWs1JJCxDEKmxOlYxTI2u8yciEWBTfdp/9kZZEwtZE7HnI51MshYBiFT4nSskkHGLGSzx8yJWCeDjGUQMjUZZCyDkClxOtZqSCBjGYRMidOximVstphZEYMqYxmETInTsUoGGcsgZGoyyFgGIVPidKySQcYsZLPBTItYJ4OMZRAyNRlkLIOQKXE61mpIIGMZhEyJ07GKZWzjYxFrOB1zOtbJIGMZhExNBhnLIGRKnI5VMsiYhWzjoo9BkrHpvk3MHTonu36XMfUm4soNxKEKmXoDcUC+ibh6A3FAuol4lzH1JuLqDcQB+Sbi6g3EAekm4l3G1JuIZ9tAfF/ZNKgs7yuzlw/N3r94FTgdczrWcTrmdAycjoHTsU6GdOx+/+reUPiruQwZZCyDkKnJIGMZhEyJe8daDQlkLIOQKXHvmNloWMRWwOmY07FOBhnLIGRqMshYBiFT4nSsYhkzQ2ARWyUZZCyDkKnJIGMZhEyJ07FWQwIZyyBkSpyOmY2ARWwNOB1zOtbJIGMZhExNBhnLIGRKnI5VLGNmVCxiI5BBxjIImZoMMpZByJQ4HWs1JJCxDEKmxOmYmU9E7IqIssTHLWs4zwsi4m0R8ZmI+H77/A8MWav+t+k6pcuYeqkL9TIXgHSpiy5j6qUu1MtcAPKlLpTLXEAVMvUyF4B8qQv1MheAfKkL5TIXUIVMucyFOcBdwKWLPL6WL87rgFPa53wL2D5AXQ/AIjYmXnfM646B1x0Drzt2oAavO+Z1x8ix7pjhzlLKjjHP8QqqgH0d+Bng6nGLWoiHJgfAvWPuHetkGKrMMFypJsNQZYbhSiXuHat4qHJ9U0q5upRyQymlTOoaTsQGxOmY0zFwOgZOxw7U4HTM6RhOx4Q8KCJeCPwb4F7gy8CnSym6X5KLYBEbGPeOuXes494x946Be8fAvWMd946tmu0RsXOxJ0opp67hPMcCVy547BsRcV4p5VMjVzcwHpqcEBmGKjMMV6rJMFSZYbhSiWdWthoSDFVmGK5U4pmVM8X7gDOpMrYN+HHgXcDxwMci4hRdaQ9E/5tyA+N0zOlYx+mY0zFwOgZOxzobJR3bVzYPKtj7ymaA69aYfB1EKeXiBQ99BXhJRNwDXAjsAJ47zjWGwonYFHA65nQMnI6B07EDNTgdczqG0zER72zHZ0qrmIdFbEp4ZqVnVnYyyFgGIVOTQcYyCJkSz6ysWMamyq3tqP3mn8dgZhARj4yI90bEtyNiT1vV9tKIeOgaz/OLEXFVRHwrInZHxE0R8RcRcdoghZYY5DSjkkHGMgiZmgwylkHIlDgdazUkkLEMQqbE6dhM0V3iJmkV8xjECiLiRGAncB7wReCt1H/kBcDnI+LoVZ7nzcDfAE8GPg5cBvwj8H8Bn2vTUMdm0169jGUQMiVOxyoZZCyDkKnJIGMZhEyJ07GKZWx8IuLxEfGwRR5/NPD29tcPzHt8a0Rsby4zdYb6LfQO4Bjg/FLK2/qDEfEH1FVp3wS8ZLkTRMSxwEXAvwL/Zynl1nnPnQ78A/BG5t28cegyNnfIxNZoW7kGrzvmdcfwumPgdccO1OB1x7zuGF53bADOAV4dEVcD3wDuBk4EfhE4FPgo8PvzXn8ccC1wM3VW5QEi4mzg7PbXY9vxtIi4ov35tlLKReMUO7aIRcQJwFnALuAPFzz9BuA3gXMj4sJSynL/ux9NTej+x3wJg7qybUTcDTx83HoXsmlvyGUMPLMSPLPSMys9sxI8sxI8s7KzUWZWCrgaOBl4EnUochtwJ/BZ6rpiV65hpfwnAi9a8NgJ7QOqvI0lYkOMj53RjleVUh5gE6WUu4HPAQ8Gnr7CeW4A9gJPjYgfmf9ERDwTOBz4uwHqPYhNeyPFcKUS945VMgxVZhiuVOLesVZDgqHKDMOVStw7tj4ppXyqlPIfSynbSylHlVK2llIeXkr5d6WU9y+UsFLKrlJKlFKOX+RcO9pzS30c9DlrZYjf/ie34/VLPH9DOz52uZOUUu4AXgU8AvhqRPxxRPy/EfEh4Crgk8BvDVDvkmSQsQxCpsS9Y5UMMpZByNRkkLEMQqbEvWMVy9jGZYjfNke2411LPN8fP2qlE5VSLo2IXcB7gd+Y99TXgSsWDlkuxVJbIwDbV/pc9465d6zj3jH3joF7x8C9Y5Crd8xsLKYRv/SYaUWziYjfAT4MXEFtrNsGnEqdgfnBiPj/JlTjQTgdczoGTsfA6VjH6ZjTMciRju0v2p/NZliG+A3TE68jl3j+iAWvW5SIeBbwZuAjpZRXznvqHyPiudShzwsj4p2llGXX/1hqa4SWlD15uc+dj9Mxp2Mdp2NOx8DpGDgdgxzpmNk4DBG5fK0dl+oBO6kdl+oh6/z7drx64ROllB9Q1yfbRJ0FMVWcjjkdA6dj4HSs43TM6RjkSMfM+meI3ypdnM6KiE3zZ05GxOHAM4DdwBdWOM+D2nGpJSr645K35E7HnI51nI45HQOnY+B0DJyO7SubBpXyfWX2dl4c+19cSrmROqvxeOBlC56+mNrn9f6+htgyK9h+ph1/MyKOm/9ERPwCVejuA64Zt+ZxcDrmdAycjoHTsY7TMadj4HTMjM5Qv0leShWkyyPiTOoKtU8DTqcOSb523muXWsH2w9R1wn4WuDYiPgLcAjyOOmwZwKtLKbcPVPPIOB1zOtZxOuZ0DJyOgdMxcDpmRmOQaKWlYk+hznZ8GnAhddbj5cBpq5GnNqT5bOqWSF8FntvO83TqdgQ/V0q5bIh6h8LpmNMxcDoGTsc6TsecjoHTMbM2BvvtUUr5JnXT75Vet4sfLmmx8Ln7gUvbx7rA6ZjTsY7TMadj4HQMnI6B0zGzemavK25COB1zOgZOx8DpWMfpmNMxcDpmVsYiNiDes9J7VnYyyFgGIVPiPStbDQlkLIOQKcmwZ6XJi0VsAmSQsQxCpsTpWCWDjGUQMjUZZCyDkClxOmayYhGbEE7HnI51MshYBiFT4nSs1ZBAxjIImRKnY2YhFrEJk0HGMgiZEqdjlQwylkHI1GSQsQxCpsTpmMmERWwKOB1zOtbJIGMZhEyJ07FWQwIZyyBkSpyOGbCITZUMMpZByJQ4HatkkLEMQqYmg4xlEDIlTseMmpkTsdAttwU4HQOnY50MMpZByJQ4HWs1JJCxDEKmxOnY7KL/bSRg833B/kN1C7BCFTL1IrCAfCFY9SKwgHQh2C5j6oVg1YvAAvKFYJWLwEIVMvUisIB8IVj1IrCAfCHY9bQI7P6yaVCR3+9Nv2eHzfcFm+8TJ1NOx5yONZyOOR0Dp2PgdAycjs0aMytiHbWMgXvHwL1j4N4xcO9YJ4OMZRAyJe4dM9Ni5kUMnI4dqCGBjGUQMjUZZCyDkClxOtZqSCBjGYRMidOxjY9FbB5qGQOnY+B0DJyOgdOxTgYZyyBkSpyOmUliEVuA07FWQwIZyyBkajLIWAYhU+J0rNWQQMYyCJkSp2MbE4vYEqhlDJyOgdMxcDoGTsc6GWQsg5ApyZCO7ZvBmYUbGX81l8HpWKshgYxlEDI1GWQsg5ApcTrWakggYxmEzJghsIitArWMgdMxcDoGTsfA6Vgng4xlEDIlGdIxs/6xiK0Sp2OthgQylkHI1GSQsQxCpsTpWKshgYxlEDJjRsUitkbUMgZOx8DpGDgdA6djnQwylkHIlDgdM6NiERsBp2OthgQylkHI1GSQsQxCpsTpWKshgYxlEDJj1oL+t8g6xntWes9K8J6VHe9Z6T0rwXtWQo49K6fF/rJp0DcA3mvSrBmnY60Gp2NOx3A6Bk7HDtTgdMzpmFkVFrGBUMsYuHcM3DsG7h0D9451MshYBiFT4t4xsxIWsQFxOtZqSCBjGYRMTQYZyyBkSpyOtRoSyFgGITNmMSxiE0AtY+B0DJyOgdMxcDrWySBjGYRMidMxsxgWsQnhdKzVkEDGMgiZmgwylkHIlDgdazUkkLEMQmZMxyI2YdQyBk7HwOkYOB0Dp2OdDDKWQciUOB2bLhFxbkSU9vHiNXzernmft/DjliFq079VnwG6jCmXuugypl7qQr3MBSBf6kK5zAVUIVMvcwHIl7pQL3MBSJe66DKmXupCvcwFIF/qYhaWuVASEY8C3gbcA4xiv3cBly7y+CBfOIvYFPG6Y153DLzuWMfrjnndMfC6YzBb645Nm4gI4H3A7cB/BS4a4TR3llJ2DFnXfGZvaFLrQe4d6zUkGKrMMFypJsNQZYbhSiXuHWs1JBiqzDBcaQbnfOAM4DxAZ9vLMHsiBmzaWz+UqGUM3DsG7h0D946Be8c6GWQsg5Apce/YcETE44BLgMtKKZ8e41QPiogXRsRrIuKCiDg9Igb75aF/Sy5k016YE/7sde9Yq8G9Y+4dw71j4N6xAzW4d8y9Y2MSEVuAK4F/AV4z5umObeeazzci4rxSyqfGPPdsJmLzcTpWcTrmdAycjoHTsY7TMadjQrZHxM7FPtZwjv8MPAn41VLKOO8q3gecSZWxbcCPA+8Cjgc+FhGnjHFuYMYTsfk4HXM6Bk7HOk7HnI6B0zFwOrYS+0sM+sZhfxk/FIiIp1JTsLeUUj4/zrlKKRcveOgrwEsi4h7gQmAH8NxxrmERm0dPxtRC5pmVnlnpmZUVz6z0zErwzEqYuZmV15VSTh3lE+cNSV4PvH7Qqh7IO6ki9sxxTzTzQ5OLkWGoUj1c6ZmVnlnZyTBUmWG4UolnVrYaEgxVZhiuNMvyEOCxwOOA++YvwAq8ob3mT9pji60NtlpubcexvyH0P+WT4nSs4nTM6Rg4HQOnYx2nY/p0bF/RvkFMzh7gPUs892Rq39hnga8B4wxbntaON41xDsAitiLuHXPvGLh3rOPeMfeOgXvHIEfvmDmY1pi/6BZGEbGDKmJ/Wkp597zHtwInAveXUm6c9/jjge+UUu5YcJ5HA29vf/3AuDV7aHIVeGZlJcNQZYbhSiWeWVnJMFSZYbhSTYahygzDlWbdcxxwLfD3Cx4/B/h2RHwsIt4REW+OiA8D1wE/BnwU+P1xL+5EbA04HXM6Bk7HOk7HnI6B0zFwOraBuRo4mZqinUbtB7uTOrR5JXBlKWXsX4YWsTXi3rGKe8fcOwbuHQP3jnXcO6bvHTPL0/aL3LHI47uAg4Z82mKtYy/YuhIemhyRDEOV6uFKz6z0zMpOhqHKDMOVSjyzstWQYKjSw5VmLVjExsC9Y5UMMpZByJS4d6ySQcYyCJmaDDKWQciMWQ0WsQHIIGNqIXM65nSsk0HGMgiZEqdjrYYEMmYhMythERsIp2OVDDKWQciUOB2rZJCxDEKmJoOMZRAyY5bCIjYwGWRMLWROx5yOdTLIWAYhU+J0rNWQQMYsZGYx9D+pNyCeWVnxzErPrATPrATPrOx4ZuXGm1k5VzYN+mZjrsxePjR7/+Ip4nTM6Rg4Hes4HXM6Bk7HwOmYeSAWsQnj3rFKBhnLIGRK3DtWySBjGYRMTQYZyyBkxljEpkQGGVMLmdMxp2OdDDKWQciUOB1rNSSQMQvZbGMRmyJOxyoZZCyDkClxOlbJIGMZhExNBhnLIGRmNpk5EYs52KztV00hY2ohczrmdKyTQcYyCJkSp2OthgQyZiGbPWZOxDoZZCyDkKnJIGMZhEyJ07FKBhnLIGRqMshYBiEzs8PMihhUGcsgZEqcjrUaEshYBiFTk0HGMgiZEqdjrYYEMmYhmw1mWsQ6GWQsg5CpySBjGYRMidOxSgYZyyBkajLIWAYhMxsbi1jD6ZjTsQM1JJCxDEKmJoOMZRAyJU7HWg0JZMxCtnGxiC0gg4xlEDI1GWQsg5ApcTpWySBjGYRMTQYZyyBkZuNhEVsEp2NOxw7UkEDGMgiZmgwylkHIlDgdazUkkLF9M7gN0EbGX81lyCBjGYRMTQYZyyBkSpyOVTLIWAYhU5NBxtRCZjYO+re6yekytl/4s2fTXv0G4oB0E/EuY+pNxNUbiAPyTcSVG4hDFTL1BuKAfBNx9QbigHQT8S5j6k3ElRuIZ2BuLgZ9gzI3p3/zP22ciK0Sp2NOx8DpGDgd6zgdczoGTsfM+FjE1oB7x9w7dqCGBDKWQcjUZJCxDEKmxL1jrQbLmBkRi9gIZJCxDEKmJoOMZRAyJU7HKhlkLIOQqckgYxYys1YsYiPidMzp2IEaEshYBiFTk0HGMgiZEqdjrQbLmFkDFrExySBjGYRMTQYZyyBkSpyOVTLIWAYhU5NBxixkZjVYxAbA6ZjTsQM1JJCxDEKmJoOMZRAyJU7HWg2WMbMCFrEBySBjGYRMTQYZyyBkSpyOVTLIWAYhU5NBxixkZiksYgPjdMzp2IEaEshYBiFTk0HGMgiZEqdjrQbLmFkEi9iEyCBjGYRMTQYZyyBkSpyOVTLIWAYhU5NBxixkZj4WsQnidMzp2IEaEshYBiFTk0HGMgiZEqdjrQbLmGlYxKZABhnLIGRqMshYBiFT4nSskkHGMgiZmgwyZiEz+reoM4L3rPSelQdq8J6V3rMS71kJ3rPyQA3reM/KuRKDvrGZKwZR784AACAASURBVPo37dNmsLfnEfHIiHhvRHw7IvZExK6IuDQiHjrCuX46Iv4yIr7TzvWdiLgqIp49dqG637+A0zFwOgZOx8DpWMfpmNMxcDo2ywzymyAiTgR2AucBXwTeCtwEXAB8PiKOXsO5Xgd8Gngm8HHgLcBfAw8FnjVEvZvlIpJDyJS4d6zVkEDGMgiZmgwylkHIlLh3rNVgGZs5hvrp8w7gGOD8Usrb+oMR8QfAK4A3AS9Z6SQRcQ7wX4C/A55XSrl7wfNbB6r3gIztlw7V6YcqQT9cqRyqhCpk6qFKQD5cqR6qBKTDlV3G1MOV6qFKQD5cqRyqhCpk6qFKYN0OV5q1Mfbb8Yg4ATgL2AX84YKn3wDcC5wbEdtWOM8m4M3AD4D/Z6GEAZRS7h+33oU4HXM6Bk7HwOlYx+mY0zFwOmamxxA/+c9ox6tKKQ94S99k6nPAg4Gnr3CenwQeA3wU+F5E/GJEvCoiLoiI0waoc0k2780hZErcO1bJIGMZhEyJe8cqGWQsg5CpySBjFrKNzRA/8U9ux+uXeP6GdnzsCuf5iXb8V+Afgb8BLgEuBa6JiE9FxMPHKXQlMshYBiFT4nSs1ZBAxjIImZoMMpZByJQ4HWs1WMbWRES8OSL+PiK+GRG7I+KOiPiniHjDGvvWIyJ+LSK+EBF3R8QP2nnOj4hBfkgO8dP+yHa8a4nn++NHrXCeY9rxJcBhwM8ChwNPAD5Bbd7/i9UUFBE7F/sAtq/0uU7HnI51MshYBiFT4nSskkHGMgiZmgwyZiFbNa8AtgGfBC4DPgjsA3YAX46IR63yPH8KvIc6YvfnwJ8Ah7Rz/nlEjP2LYho/XXqRK3VD95/4AbyglPI/29//OSKeS03cfiYiTiulfH4CdT6wmL36Rn7QN/OrG/nB64553TGvOwZedwy87lhnf/Fa7KvgiFLKQd8oEfEm4DXA7wIvXe4EEXE2cC7wDeCppZTb2uNbgQ8BzwdeBFwxTqFDfDV74nXkEs8fseB1S/G9drxpnoQBUErZTU3FAJ66UkGllFMX+wCuW+lz5+N0zOlYx+mY0zFwOgZOxzrqdMwsz2IS1vhQO560itM8rx3f0iWsnft+4PXtr789WoU/ZIif7F9rx6V6wPo/dqkesoXnuXOJ57uoTf27P4OMZRAyJe4dazUkkLEMQqYmg4xlEDIl7h0zI/JL7fjlVbz22Ha8aZHn+mNPjoiVWq+WZYifJle341kRsWn+zMmIOBx4BrAb+MIK5/k0dfz2pIg4pJSy8Ff/E9px1/glrx2vO+Z1xzped8zrjoHXHQOvO9ZRrztmliYiLgIeQh21ewrwU1QJu2QVn95TsMcs8twJ8/68nZUdZ0nGfntdSrkRuAo4HnjZgqcvpjbLvb+Uci/UsdWI2N5W459/ntuojXBHAv95/nMR8e+An6MOb3583JrHwemY0zFwOgZOxzpOx5yOgdOxgdm+zKS7tXIRdU3Tl1Ml7OPAWaWU767ic/+mHV8ZEQ/rD0bEFqrfdNa8leN8hvoJ8lLgGuDyiDgTuBZ4GnA6dUjytfNee1x7/maqvM3nle3zXhsRz6Rul/Ro4LnAfuA3SilLDV1ODadjTsc6TsecjoHTMXA61pm5dGxu4DdEA/8oK6UcCxARj6CuV3oJ8E8R8e9LKf+4wqf/GfBC4BeAr0bEf6cuOv+zwInU5blOovrJyAzylrqlYk+hzhx4GnAhtcjLgdNKKbev8jy3ts9/K/Ao4HzqgrF/C/x0KWVVy1dMC6djTsfA6Rg4Hes4HXM6Bk7HBuC6ZSbdjUQp5V9LKR+h7gR0NPD+VXzOHPAcaqp2C3UG5a8B36Kma91tbh21Lhhw+YpSyjepm36v9Lpd/HBJi8Wev4OajL1yqNomidMxp2Mdp2NOx8DpGDgd68xcOrYOKKXcHBFfBZ4YET8yfzbkEq/fB7ylfRwgIg4Dnkjtgf/ncWryYiQD4XTM6Rg4HQOnYx2nY07HwOlYUn60Hcd513gucCjwoXH3wbaIDYjXHfO6Y50MMpZByJR43bFKBhnLIGRqLGPTo00IPHaRxze1BV2PAa4ppXyvPb7oJML23BGLPPYT1F6ze4A3jluv/m3jBsSr8ntVfvCq/OBV+Tteld+r8kOOVflnhJ8Hfi8iPg3cSO3legTwM9RlJ24BfmPe65ebRPjJiNgNfAW4G3g88GxgD/C8Uspia4ytCYvYhHDvmHvHOu4dc+8YuHcM3DvWce/YxPk74I+p65ieQt3r+l7qKg5XApe3fvTV8GHg/6bOnjwM+DbwbuCS1vM+NhaxCeN0zOkYOB0Dp2Mdp2NOx8Dp2CQppXyFg9c1Xe71u1hiEmEp5feA3xumssVxj9gUcO+Ye8c67h1z7xi4dwzcO9Zx75ixiE2RDDKWQciUeGZlqyGBjGUQMjUZZCyDkCnxzEqjZuZELHQjM4DTMXA61skgYxmETInTsUoGGcsgZGosY7PJzIkY5EiGMsiY+h5kkDG1kDkdczrWySBjGYRMidMxo2AmRayjFhGnY07HOhlkLIOQKXE6VskgYxmETI1lbHbQvw0Uk2FWoWdWemYleGYleGZlxzMrPbMS1sfMylJi0DdRpejfGE+bmU7E5qNOhpyOOR3rOB1zOgZOx8DpWMfp2MbGIjaPDH1TGWRMfQ8yyJhayNw75t6xTgYZyyBkStw7ZiaJRWwR1CLidMzpWCeDjGUQMiVOxyoZZCyDkKmxjG08LGJLkCEZyiBj6nuQQcbUQuZ0zOlYJ4OMZRAyJRnSsf3Fv7o3Ev5qroBaRJyOOR3rZJCxDEKmxOlYJYOMZRAyY4bAIrYKMiRDGWRMfQ8yyJhayJyOOR3rZJCxDEKmJEM6ZtY/FrE1oBYRp2NOxzoZZCyDkClxOlbJIGMZhMyYUbGIrZEMyVAGGVPfgwwyphYyp2NOxzoZZCyDkClxOmZGxSI2ImoRcTrmdKyTQcYyCJkSp2OVDDKWQciMWQsWsTHIkAxlkDH1PcggY2ohczrmdKyTQcYyCJkSp2NmLVjEBkAtIk7HnI51MshYBiFT4nSskkHGMgiZMSthERuIDMlQBhlT34MMMqYWMqdjTsc6GWQsg5ApcTpmVkL/k2KDsfk+/QbioN9EXL2BOOg3EVduIA5VyNQbiAPyTcTVG4gD0k3Eu4ypNxFXbyAOyDcRV24gPimiDPvGK7Q/NiU4EZsAGZIhp2NOx8DpGDgd6zgdczpmcmIRmyBqEXHvmHvHOhlkLIOQKXHvWCWDjGUQMmM6FrEJkyEZyiBj6nuQQcbUQuZ0zOlYJ4OMZRAyJU7HTMciNiXUIuJ0zOlYJ4OMZRAyJU7HKhlkLIOQmdnGIjZFMiRDGWRMfQ8yyJhayJyOOR3rZJCxDEKmxOnYbDNzIpZhRoZaRJyOOR3rZJCxDEKmxOlYJYOMZRAyM3vMnIgBbN5T2LxHa2QZkqEMMqa+BxlkTC1kTsecjnUyyFgGIVPidGz2mEkR66hlDPQi4nTM6Vgng4xlEDIlTscqGWQsg5CZ2WCmRQycjh2oQS4i+nuQQcbUQuZ0zOlYJ4OMZRAyJU7HZoOZF7GOWsZALyJOx5yOdTLIWAYhU+J0rJJBxjIImdm4WMTm4XSs1SAXEf09yCBjaiFzOuZ0rJNBxjIImRKnYxsXi9giqGUM9CLidMzpWCeDjGUQMiVOxyoZZCyDkJmNhf6tVlK6jO1/kO6XUBcR9Sbi6g3EQb+JuHoDcUC6iXiXMfUm4uoNxAH5JuLKDcShCpl6A3FAvom4cgPx/UX/Bu0AZeDkPNO/bUo4EVsBp2NOx8DpWMfpmNMxcDoGOdIxszGwiK0C9461GuQior8HGWRMLWTuHXPvWCeDjGUQMmPGwSK2BtQyBnoRcTrmdKyTQcYyCJkSp2OVDDJmITOjYhFbI07HWg1yEdHfgwwyphYyp2NOxzoZZCyDkBmzVixiI6KWMdCLiNMxp2OdDDKWQciUOB2rZJAxC5lZCxaxMXA61mqQi4j+HmSQMbWQOR1zOtbJIGMZhMzoiIg3R8TfR8Q3I2J3RNwREf8UEW+IiKOnfZ7lsIgNgFrGQC8iTsecjnUyyFgGIVPidKySQcYsZDJeAWwDPglcBnwQ2AfsAL4cEY+a8nmWRP/WaYPgdcdaDV53zOuO4XXHwOuOdbzumH7dsRnliFLKQW/PI+JNwGuA3wVeOsXzLIkTsYFxOuZ0DJyOdZyOOR0Dp2PgdGzaLCZPjQ+140nTPM9yWMQmgHvHWg1yEdHfgwwyphYy9465d6yTQcYyCJmR8kvt+OUk5/HQ5CTZvKdIhyqhioh6qBL0w5XqoUrQD1cqhyqhCpl6qBKQD1eqhyoB6XBllzH1cKV6qBLwcOXSbI+InYs9UUo5dS0nioiLgIcARwJPAX6KKk+XKM6zGBaxCePesVaDe8fcO4Z7x8C9Yx33jm2M3rGYG7YNIob/r3kR8Ih5f/848KullO+KznMQHpqcEuqhStAP07l3zL1jnQxDlRmGK5W4d6ySYajSw5UHcV0p5dTFPtZ6olLKsaWUAI4FngecAPxTRDxZcZ7FsIhNEfeOtRrkIqK/BxlkTC1k7h1z71gng4xlEDIzOUop/1pK+QhwFnA08H7leeYzeyI2V9hyn1iGnI45HcPpWCeDjGUQMiVOxyoZZMxCNllKKTcDXwUeHxE/oj4PzKKINTLImFrIMiRDGWRMfQ8yyJhayJyOOR3rZJCxDEJmJsqPtuO4jZqDnGdmRQyqjGUQMjVqEXE65nSsk0HGMgiZEqdjlQwyZiEbjYjYHhHHLvL4prYQ6zHANaWU77XHt7bPOXGc84yK/u1PArbcV9h3qHBWo2dW1ho8s9IzK/HMSvDMyo5nVm6MmZUCfh74vYj4NHAjcDt1xuPPUJvsbwF+Y97rjwOuBW4Gjh/jPCNhEWv0ZEwtZF53rB7VQqaWMdALmdcd87pjXnes4nXH1h1/B/wx8AzgFOAo4F7geuBK4PJSyh1TPM+yWMQW4HQsRzLkdMzpGDgdA6djHadjTsdWSynlK8DL1vD6XcBBv3TXep5RmekesaVw71irwb1j8nvg3rGKe8fcOwbuHYMqY3PFv7o3Ev5qLkMGGVMLWYZZhRlkTH0PMsiYWsg8s9IzKzsZZEwtZGbjYBFbAadjrYYEMpZByJQ4HatkkLEMQqbE6VjFMmaGwCK2SjLImFrIMiRDGWRMfQ8yyJhayJyOOR3rZJAxC5kZB/3/onWEZ1a2GjyzUn4PPLOy4pmVnlkJnlkppQz85lA/ADR1nIiNgNOxHMmQ0zGnY+B0DJyOdZyOmfWIRWxE3DvWakggYxmETIl7xyoZZCyDkClx71jFMmbWgkVsTDLImFrIMiRDGWRMfQ8yyJhayJyOOR3rZJAxC5lZDRaxAXA61mpIIGMZhEyJ07FKBhnLIGRKnI5VLGNmJSxiA5JBxtRCliEZyiBj6nuQQcbUQuZ0zOlYJ4OMWcjMUljEBsbpWKshgYxlEDIlTscqGWQsg5ApcTpWsYyZxbCITYgMMqYWsgzJUAYZU9+DDDKmFjKnY07HOhlkzEJm5mMRmyBOx1oNCWQsg5ApcTpWySBjGYRMidOximXMdCxiUyCDjKmFLEMylEHG1Pcgg4yphczpmNOxTgYZs5AZi9iUcDrWakggYxmETInTsUoGGcsgZEqcjlUsY7PNYD8FIuKREfHeiPh2ROyJiF0RcWlEPHSMc54bEaV9vHiQOsUukkHG1EKWIRnKIGPqe5BBxtRC5nTM6Vgng4xZyGaTQX4CRMSJwE7gPOCLwFuBm4ALgM9HxNEjnPNRwNuAe4aocT5qGXI61mpIIGMZhEyJ07FKBhnLIGRKnI5VLGOzx1Dfce8AjgHOL6W8rT8YEX8AvAJ4E/CS1Z4sIgJ4H3A78F+Biwaq8wAZNvDecl+RbyAOSDcR7yKi3kRcvYE46DcRV28gDkg3Ee8ypt5EXL2BOCDfRFy5gThUIVNvIA6si03EY27YN5Sh+/aXMfZbsIg4ATgL2AX84YKn3wDcC5wbEdvWcNrzgTOoCdu949a4HOpkyulYq8HpmPweOB2rOB1zOgZOx8z0GOJ/+xnteFUp5QEuW0q5G/gc8GDg6as5WUQ8DrgEuKyU8ukB6luRDDKkvr57x1oNchHR34MMMqYWMveOuXesk0HGLGQbmyH+p5/cjtcv8fwN7fjYlU4UEVuAK4F/AV4zfmlrQy1DGYRQLWOgFxGnY07HOhlkLIOQKXE6VrGMbVyG+M46sh3vWuL5/vhRqzjXfwaeBPxUKWX3qAVFxM4lntq+0ue6d8y9YwdqcO+Ye8dw7xi4d6yTpXfMbCym8Var/zZf9qdYRDyVmoK9pZTy+YlXtQLqZMrpWKvB6Zj8HjgdqzgdczoGOdKxuTn9/wczHEN8N/XE68glnj9iwesOYt6Q5PXA68ctqJRy6hLX2Qk8ebXncTrmdOxADU7HnI7hdAycjnXU6ZjZOAzx9upr7bhUD9hJ7bhUDxnAQ9rnPw64b94iroU68xLgT9pjl45d8RpRJ1NOx1oNTsfk98DpWMXpmNMxyJGOmfXPEN9BV7fjWRGxaf7MyYg4HHgGsBv4wjLn2AO8Z4nnnkztG/ssVfokw5ZOx5yOHajB6ZjTMZyOgdOxjtMxMw5ji1gp5caIuIq6ltjLqKvhdy4GtgHvKqXcCxARW4ETgftLKTe2c+wGFt3CKCJ2UEXsT0sp7x633nFRy1AGIdy8p0hlDKqMqGUM9EKmljHQC5lSxqAKmVrGALmQqWUMkApZT8YsZGatDJWpvhS4Brg8Is4ErgWeBpxOHZJ87bzXHteevxk4fqDrT5UMMqQWQqdjrQanY07HcDoGTsc6TsfMWhmkyaAlW08BrqAK2IXU1Oty4LRSyu1DXCcb6r4t9461Gtw7Jr8H7h2ruHfMvWPg3jGzNgb7TimlfJO6JdFKr9vFD5e0WM15dwA7Rq1r0jgdczp2oAanY07HcDoGTsc6TsfMarCyD4RahjIIoXvH3DsG7h3ruHfMvWMwA71jZeARAf0gy9TRZtgbjAxDherre8/KVoN8mE5/DzIMVaqHK71npfes7Hio0iyFRWwCqGUogxCqZQz0IuLeMfeOdTLIWAYhU+LeMZMVi9iEyCBD6us7HWs1yEVEfw8yyJhayJyOOR3rWMbMfCxiE0YtQxmEUC1joBcRp2NOxzoZZCyDkClxOmYyYRGbAhlkSH19p2OtBrmI6O9BBhlTC5nTMadjHcuYsYhNEbUMZRBCtYyBXkScjjkd62SQsQxCpsTp2MYjIo6OiBdHxEci4usRsTsi7oqIz0bEr0fEqr7pI+JX5+99vcTH2NOCZ+4rH7qZ3ECOZSbUS2143bFWg9cd87pjeN0x8LpjHa87NhjnAH8EfIe6H/a/AI8Ange8G/iFiDinlLLSf7ovUbdqXIyfBs4APjZusTMnYgCbd9cfOPsP070TVMtQBiH0umNedwy87ljH64553TGYgXXHpsP1wHOAvy2lHPiGjojXAF8Enk+Vsr9c7iSllC9RZewgIuLz7Y9/PG6xMz002YVMRYahQvX13TvWapAP0+nvQYahSvVwpXvH3DvW8VDl6JRS/qGU8tfzJaw9fgvwzvbXZ416/oh4AvB04H8DfzvqeTozLWJQZSyDkKmvr65BLWOgFxH3jrl3rJNBxjIImRL3jm1Y7m/HcSLH32rH95RS3CM2FJt3z8mHKsG9Y+DeMfeOuXcM3DsG7h3rzGDv2PaI2LnYE6WUU0c9aURsAX6l/fXjI57jMOCFwBy132xsZj4Rm4/TMadjB2pwOia/B07HKk7HnI6B07GBuAR4AvDRUsonRjzHfwCOAj5WSvnmEEX5q7oITsecjkGOZMjpmNMxcDoGTsc62WQs5oZ909ZWNrhunORr0fNGnA9cCFwHnDvGqX6zHd81dlENJ2JL4HTM6diBGpyOye+B07GK0zGnYwBzRf+9uJ6IiJcBlwFfBU4vpdwx4nn+LfCTwLeAjw5Vn0VsBTLImFqG1Nf3zMpWg1xE9Pcgg4yphcwzKz2z0qyeiHg58HbgK1QJu2WM0w3apN+xiK0Cp2M5hFAtY6AXEadjTsc6GWQsg5ApyZCOmaWJiFcBb6WuBXZ6KeXWMc51KHVIcw54zzAVVixiayCDjKllSH19p2OtBrmI6O9BBhlTC5nTMadjZnEi4vXU5vydwJmllNuWee3WiNgeEScuc8pzgIdSG/0HadLv+LtnjXhV/hyTCbwqv1flB6/K3/Gq/F6V3/yQiHgR8EZgP/AZ4PyIg35f7CqlXNH+fBxwLXAzcPwSp+1N+mOvpL8Qi9iIeGalXgg9s7LV4JmVnlmJZ1aCZ1aaAzymHTcDL1/iNZ8CrljNySLiccBPMXCTfsdDk2Pg3rEcw6XqoUrQD9O5d8y9Y50MQ5UZhiuVuHdMSyllRyklVvh41rzX72qPHb/E+a5tzz9qyCb9jkVsADLImFqG1Nd371irQS4i+nuQQcbUQubeMfeOmfWDRWwgnI7lEEK1jIFeRJyOOR3rZJCxDEKmxOmYWQmL2MBkkDG1DKmv73Ss1SAXEf09yCBjaiFzOuZ0zOTGIjYBnI7lEEK1jIFeRJyOOR3rZJCxDEKmxOmYWQx/R0wQz6z0zErIMavQMys9sxI8sxI8s3Joogz7pjf075+njhOxCeN0zOnYgRqcjsnvgdOxitMxp2MmDxaxKZFBxtQypL6+e8daDXIR0d+DDDKmFjL3jrl3zOTAIjZFnI7lEEK1jIFeRJyOOR3rZJCxDEKmxOnYbDNzIhZzCX4JJ5AxtQypr+90rNUgFxH9PcggY2ohczrmdMzomDkRA9iyez9bdusaNcHpWL++uga1jIFeRJyOOR3rZJCxDEKmxOnY7DGTItZRyxg4Hes1KHE61mqQi4j+HmSQMbWQOR1zOmamy0yLGDgd66hlKIMQqmUM9CLidMzpWCeDjGUQMiVOx2aDmRexjlrGwOlYr0GJ07FWg1xE9Pcgg4yphczpmNMxM3ksYvNwOlZRy1AGIVTLGOhFxOmY07FOBhnLIGRKnI5tXCxii6CWMXA61mtQ4nSs1SAXEf09yCBjaiFzOuZ0zEwGi9gSOB2rqGUogxCqZQz0IuJ0zOlYJ4OMZRAyKdpfC2ZgLGIroJYxcDrWa1DidKzVIBcR/T3IIGNqIXM6liMdMxsDZ5yroMvYvsN0/+m6jKk3EVdvIA7aTcw37ynSDcShioh6A3HQbyKu3kAc9JuIKzcQhypk6g3EAfkm4soNxOXMDfxGPcGi69PGidgacDrmdAycjh2oQZ4M6e+B0zGnY+B0zIyHRWyNuHesopahDEKoljHQi4h7x9w71skgYxmEzJi1YhEbEbWMgdOxXoMSp2OtBrmI6O9BBhlTC5nTMadjZu1YxMbA6VhFLUMZhFAtY6AXEadjTsc6GWQsg5AZsxosYgOgljFwOtZrUOJ0rNUgFxH9PcggY2ohczrmdMysDovYQDgdq6hlKIMQqmUM9CLidMzpWCeDjGUQMmOWwiI2MGoZA6djvQYlTsdaDXIR0d+DDDKmFjKnY07HzNJYxCaA07GKWoYyCKFaxkAvIk7HnI51MshYBiEzZj4WsQmiljFwOtZrUOJ0rNUgFxH9PcggY2ohczrmdMw8EIvYhHE6VlHLUAYhVMsY6EXE6ZjTsU4GGcsgZMZYxKaEWsbA6VivQYnTsVaDXET09yCDjKmFzOmY0zFjEZsqTscqahnKIIRqGQO9iDgdczrWySBjGYTMzCazt+l3KWy+bx/7D9X907fs3i/dQByqkKk3EAftBt7qTcy7jCk3Ee8iot5EXL2BOOg3EVdvIA5INxHvMqbeRFy9gTiwrjYRjzLsm+vQv0edOjObiG2+b5/0+k7HKupkyulYq8HpmPweOB2rOB1zOjZrzKyIQZWxDEKmJoOMqWVIfX33jrUa5CKivwcZZEwtZO4dc+/YLDHTItbJIGNqIXM6lkMI1TIGehFxOuZ0rJNBxjIImdnYWMQaTscqGWRMLUPq6zsdazXIRUR/DzLImFrInI45HdvoWMQWkEHG1ELmdCyHEKplDPQi4nTM6Vgng4xlEDKz8bCILYLTsUoGGVPLkPr6TsdaDXIR0d+DDDKmFjKnYy0dK3oxzkxEHB0RL46Ij0TE1yNid0TcFRGfjYhfj4hVfxEj4gUR8baI+ExEfD8iSkR8YMh6Z2/5ijWQYZkLQLrURZcx9VIX6mUuQLvUxuY9RbrMBVQRUS9zAfqlLtTLXIB+qQvlMhdQhUy9zAUgXerCLMs5wB8B3wGuBv4FeATwPODdwC9ExDmllNV8E70OOAW4B/gWsH3oYp2IrYDTsYrTMadjkCMZcjrmdAycjplluR54DvDIUsovl1J+t5Tya1SJ+ibwfKqUrYZXAI8FjgD+0ySK9XfRKskgY2ohc+9YDiFUyxjoRcS9Y+4d62SQMQtZLkop/1BK+etSytyCx28B3tn++qxVnuvqUsoNq0zPRsLfPWvA6Vglg4ypZUh9fadjrQa5iOjvQQYZUwuZ0zGzBu5vR+0v83m4R2wE3Dvm3rF+fXDvmLpvyr1j7h3ruHdsw7A9InYu9kQp5dRRTxoRW4BfaX/9+KjnGRor/Ig4Has4HXM6BjmSIadjTsfA6ZhZlkuAJwAfLaV8Ql1Mx4nYmDgdczrWrw9Ox9TJkNMxp2Mdp2PTIeaGfUMe9VTXjZN8LXreiPOBC4HrgHOHPPe4WNsHwOlYxemY0zHIkQw5HXM6Bk7HTCUiXgZcBnwVOL2Ucoe4pAfg75ABySBjaiHzzMocQqiWMdCLiGdWemZlJ4OMWcg0RMTLgbcDX6FK2C3ikg7C3xkD43SskkHG1DKkvr7TsVaDXET09yCDjKmFzOnY7BERrwLeCnyJKmG3iktaFH9XTIgMMqYWMqdjOYRQLWOgFxGnY07H3HtJPAAAGxZJREFUOhlkzEI2eSLi9dTm/J3AmaWU25Z57daI2B4RJ06twHm4WX+CdBlTN/MrG/mhCpm6kR+0jfTqyQRdxpTN/F1E1M386kZ+0Dfzqxv5AWkzf5cxdTP/Rm/kVxERLwLeCOwHPgOcH3HQz75dpZQr2p+PA64FbgaOX3Cus4Gz21+PbcfTIqJ/7m2llIvGqdciNgU8s9IzK/v1wTMr1bMKPbPSMys7nlm5YXlMO24GXr7Eaz4FXLGKcz0ReNGCx05oH1DlbSwRcz46Jdw7VskwVKkeKlRf371jrQb5MJ3+HmQYqlQPV7p3bONRStlRSokVPp417/W72mPHj3Cugz5nrQz21Y+IR0bEeyPi2xGxJyJ2RcSlEfHQVX7+0RHx4oj4SER8PSJ2R8RdEfHZiPj1iBik1phT/wLSy5hayNw7lkMI1TIGehFx75h7xzoZZMxCNpsMIze1wW0ncB7wReoshZuAC4DPR8TRqzjNOcCfAE8D/gdwKfCX1FVw3w18KBYZ5B2FTT/Q/tRxOlbJIGNqGVJf3+lYq0EuIvp7kEHG1ELmdMwoGOor/g7gGOD8UsrZpZRXl1LOoArZycCbVnGO64HnAI8spfxyKeV3Sym/BmwHvgk8H3jeQPWy6Qd7UwiZEqdjFbUMZRBCtYyBXkScjjkd62SQMQvZ7DD2VzoiTgDOAnYBf7jg6TcA9wLnRsS25c5TSvmHUspfl1LmFjx+C/DO9tdnjVvvQjLIWAYhU5NBxtQypL6+07FWg1xE9Pcgg4yphczpmJkWQ3yVz2jHqxaRqLuBzwEPBp4+xjXub8eJGIvTMadjHbUMZRBCtYyBXkScjjkd62SQMQvZxmaINRVObsfrl3j+Bmpi9ljg79d68ojYAvxK++vH11zdGtj0g73MPVg3n9vrjlW87ph+qQ2vO9Zq8LpjXncMrzu2HDFXBn0Tr55Qp2CI3/hHtuNdSzzfHz9qxPNfQm3Y/2gp5ROr+YSI2LnEU9tX+tyejKmFTC1j4HXH1DKUQQi97pjXHQOvO9bJsu6Y2VhM46vaf4qv+bs3Is4HLgSuA84dsqiVyDBUmWG4Uk2GoUr1UKH6+u4dazXIh+n09yDDUKV6uDJD71jMXmi0oRkidumJ15FLPH/Egtetioh4GXAZ8FXqPlF3rPZzSymnLnHOncCTV3sep2NOxzpOx5yOgdMxcDrWUadjZuMwxG+2r7XjY5d4/qR2XKqH7CAi4uXA24GvUHdMv2X08sbH6ZjTMXA6Bk7HDtQgT4b098DpWI50zKx/hhCxq9vxrIWr30fE4cAzgN3AF1Zzsoh4FXX9sS9RJezWAWocG8+s9MzKjlqGMgihWsZALyKeWemZlR3LmBmHsUWslHIjcBV1x/KXLXj6YmAb8P5Syr0AEbE1Ira31fgfQES8ntqcv5M6HHnbuPUNTQYZyyBkajLImFqG1Nd3OtZqkIuI/h5kkDG1kDkdM6MyVPPRS4FrgMsj4kzgWupWRadThyRfO++1x7Xnb6bKGwAR8SLgjcB+4DPA+YvsaLSrlHLFQDWPjHvH3DvWce+Ye8fAvWPg3rGOe8fMWhnkN3kp5caIeApVpH4eeDbwHeBy4OJVNto/ph03Ay9f4jWfAq4Yr9rh8LpjXncMcsiQWgi97lirweuOed0xcqw7ZtYPg/0GL6V8k7rp90qv28UPl7SY//gOYMdQ9UwLp2NOxzpqGcoghE7HnI6B07GO0zGzGrw63EC4d8y9Y+DeMXDv2IEa5H1T+nvg3jH3jpmVsYgNiGdWemZlRy1DGYRQLWOgFxHPrPTMyo5lzCyFbjxrA+PeMfeOQY6hQvVwqXvHWg3uHXPvGBu0d6yUYQOAsoHuzSpxIjYhnI45HeuokymnY60Gp2Pye+B0rOJ0zMzHIjZhMshYBiFTk0HG1DKkvr57x1oNchHR34MMMqYWMveOmY5FbAo4HXM61lHLUAYhVMsY6EXE6ZjTsY5lzFjEpkgGGcsgZGoyyJhahtTXdzrWapCLiP4eZJAxtZA5HZttZk/E5sTryjgdczrWUMtQBiFUyxjoRcTpmNOxjmVsNpk9EQNi915it/Z/fQYZyyBkajLImFqG1Nd3OtZqkIuI/h5kkDG1kDkdmz1mUsQ6GWQsg5ApcTpWUctQBiFUyxjoRcTpmNOxjmVsdphpEQOnY+B0rJNBxtQypL6+07FWg1xE9Pcgg4yphczp2Gww8yLWySBjGYRMidOxilqGMgihWsZALyJOx5yOdSxjGxuL2Dycjjkd62SQMbUMqa/vdKzVIBcR/T3IIGNqIXM6tnGxiC1CBhnLIGRKnI5V1DKUQQjVMgZ6EXE65nSsYxnbeFjElsDpmNOxTgYZU8uQ+vpOx1oNchHR34MMMiYXsmIZ20h40+8ViN17KYfpdqrtMqbeRFy9gTgg3US8y5h6E3H1BuKg3cR8854i3UAcqoioNxAH/Sbi6g3EQb+JuHID8SzEXBk0NAjxWp8KnIitAqdjTsc6TsecjkGOZMjpmNMxszGwiK2BDDKWQciUuHesopahDEKoljHQi4h7x9w7ZtY/FrE14nTM6Vgng4ypZUh9fadjrQa5iOjvQQYZs5CZUbCIjUgGGcsgZEqcjlXUMpRBCNUyBnoRcTrmdMysTyxiY+B0zOlYJ4OMqWVIfX2nY60GuYjo70EGGbOQmdViERuADDKWQciUOB2rqGUogxCqZQz0IuJ0zOnYLBMRL4iIt0XEZyLi+xFRIuIDI5xnV/vcxT5uGapeL18xEF3G1EtdqJe5AORLXSiXuYAqZOplLkC7zIR6qY0uY8qlLrqIqJe6UC9zAfqlLtTLXABe6mK6vA44BbgH+BawfYxz3QVcusjj94xxzgdgERsYrzvmdcfA647164PXHVOvueV1x7zu2AzyCqqAfR34GeDqMc51ZyllxxBFLYWHJieAe8fcO9bJMFSpHipUX9+9Y60G+TCd/h5kGKr0cOXkKaVcXUq5oZSyLszXidgEcTrmdAycjvXrg9MxdTLkdMzpmFkzD4qIFwL/BrgX+DLw6VLKYO/0LWITxr1j7h3ruHdML4TuHWs1uHfMvWO52R4ROxd7opRy6pRrORa4csFj34iI80opnxriAh6anBIZhiozDFcq8czKinqoMMNwqXqoEvTDdJ5Z6ZmVZkXeB5xJlbFtwI8D7wKOBz4WEacMcREnYlPE6ZjTsY7TMadjkCMZcjrmdGws5sqwQUPd9Ps6QfJ1EKWUixc89BXgJRFxD3AhsAN47rjXmb1EbE6bRoDTMXA6Bk7H+vXVNTgdczoGTsfMmnhnOz5ziJPNnogB/GB3/RDimZWeWdnJIGNqGVJf3zMrWw1yEdHfgwwyZiFLz63tuG2Ik82miHXEMgZOx8DpGDgd69dX16CWMdCLiNMxp2NmRU5rx5uGONlsixg4HWtkkLEMQqYmg4ypZUh9fadjrQa5iOjvQQYZs5BNlojYGhHbI+LEBY8/PiIetsjrHw28vf11zdsmLYab9Ts/2A0PPkxagtcd87pj4HXH+vXB646p19zyumNed2w9EhFnA2e3vx7bjqdFxBXtz7eVUi5qfz4OuBa4mTobsnMO8OqIuBr4BnA3cCLwi8ChwEeB3x+iXovYfHoyJhQyz6z0zMqOZ1bqhdAzK1sNnlnpmZXriycCL1rw2AntA6p0XcTyXA2cDDyJOhS5DbgT+Cx1XbErh1q53yK2GE7HnI7hdKyjlqEMQuh0zOkY5EnHzPK0vSF3rPK1u4CDbmpbrHWQBVtXwj1iS+HeMcC9Y+DeMXDvGLh37EAN8r4p/T1Q946FfhUmMyAWsZXwzErPrMQzKztqGcoghGoZA72IeGZljpmVZmNgEVsNTscAp2PgdAxyyJD6+k7HWg0JZEx9DyxjZlwsYmvB6ZjTMZyOddQylEEI1TIGehFxOuZ0zIyHRWytOB0DnI6B0zHIIUPq6zsdazUkkDH1PbCMmVHwrMlR8cxKz6zEMys7nlnpmZXgmZWQY2blVJmbGzacSLAf9LRxIjYOTscAp2PgdAycjoHTsQM1OB1zOmZWjUVsCNw75t4x3DvWUctQBiFUyxjoRcS9Y+4dM6vDIjYUTscAp2PgdAxyyJD6+k7HWg0JZEx9DyxjZjksYkPjdMzpGE7HOmoZyiCEahkDvYg4HXM6ZpbGIjYJnI4BTsfA6RjkkCH19Z2OtRoSyJj6HljGzEIsYpPE6ZjTMZyOddQylEEI1TIGehFxOuZ0zDwQi9ikcToGOB0Dp2OQQ4bU13c61mpIIGPqe2AZM2ARmx5Ox5yO4XSso5ahDEKoljHQi4jTMadjxiI2XZyOAU7HwOkY5JAh9fWdjrUaEsiY+h5YxmaX2ROxuTmKOp1SXx+nY+B0DJyO9eura1DLGOhFxOmY07FZZfZErJFCxpyOpZCxDEKmJoOMqWVIfX2nY62GBDKmvgeWsdliZkUMqoylEDIxGWQsg5ApcTpWUctQBiFUyxjoRcTpmNOxWcKbflOFLJQbeHcZE9bQZUy9ibh6A3FAvom4cgNxqEKm3kActBt4qzcx7zKm3ES8i4h6E3H1BuKg30Q89QbiQ7f7eNPv2cXpWMXpmNMxcDrWr6+uwemY0zFwOrbRsYgtIIWMuXcshYxlEDI1GWRMLUPq67t3rNWQQMbU98AytjGxiC2C07FKBhnLIGRKnI5V1DKUQQjVMgZ6EXE61mRM/61gBsQitgwpZMzpWAoZyyBkajLImFqG1Nd3OtZqSCBj6ntgNg4WsRVwOlbJIGMZhEyJ07GKWoYyCKFaxkAvIk7HzEbBIrZKUsiY07EUMpZByNRkkDG1DKmv73Ss1ZBAxtT3wKxvLGJrwOlYJYOMZRAyJU7HKmoZyiCEahkDvYg4HTPrGYvYCKSQMadjKWQsg5CpySBjahlSX9/pWKshgYyp74FZf1jERsTpWCWDjGUQMiVOxypqGcoghGoZA72IOB0z6w2L2JikkDGnYylkLIOQqckgY2oZUl/f6VirIYGMqe+BWR9YxAbA6Vglg4xlEDIlTscqahnKIIRqGQO9iDgdM+sB7zU5IN6z0ntWgves7HjPSu9ZCTn2a/SelZOj7J9j7p57Bz3frOFEbGCcjlWcjjkdA6dj/frqGpyOOR0zebGITYgUMubesRQylkHI1GSQMbUMqa/v3rFWQwIZU98DkwuL2ARxOlbJIGMZhEyJ07GKWoYyCKFaxkAvIk7HTCYsYlMghYw5HUshYxmETE0GGVPLkPr6TsdaDQlkTH0PjB6L2JRwOlbJIGMZhEyJ07GKWoYyCKFaxkAvIk7HjJrBRCwiHhkR742Ib0fEnojYFRGXRsRDFedZCvWMjBQy5nQshYxlEDI1GWRMLUPq6zsdazUkkDH1PdhoDOETUfm1iPhCRNwdET+IiH+KiPMjYpCp6YOIWEScCOwEzgO+CLwVuAm4APh8RBw9zfOsxJBTbUfB6Vglg4xlEDIlTscqahnKIIRqGQO9iDgd2zgM6BN/CrwHeAzw58CfAIcAlwF/HhFjrw0zVCL2DuAY4PxSytmllFeXUs6g/sNPBt405fOsyNw996YQMilOxwCnY+B0DHLIkPr6TsdaDQlkTH0PNgBj+0REnA2cC3wDeHwp5cWllAuAJwL/DXg+8KJxCx1bxCLiBOAsYBfwhwuefgNwL3BuRGybxnnWSgYZSyFkYjLIWAYhU+J0rKKWoQxCqJYx0IuI07H1y4A+8bx2fEsp5bb+YCnlfuD17a+/PW69QyRiZ7TjVaWUB/wELaXcDXwOeDDw9CmdZ804HcPpWCODjGUQMjUZZEwtQ+rrOx1rNSSQMfU9WIcM5RPHtuNNizzXH3tyRBw1aqEwjIid3I7XL/H8De342CmdZ2QyyFgKIROTQcYyCJkSp2MVtQxlEEK1jIFeRJyOrTuG8omegj1mkedOmPfn7ausa1GG2AzvyHa8a4nn++MrGeNQ5yEidi7x1Cn38n2+cP/Hl/7k77VzbBau7PF9YFOClUVS1KDbIw+giK8PwPi9oGOR4R6UBN+KRX0bMnwd9CXovxcS3IO9t/8rwPHiMgBY8XfqCOcDti/1e7yUcuoqTjOUT/wN8B+BV0bEn5VS7gCIiC3AxfNeN9aqDtPYlbh/2477tmqI82yaY//+u/ne/1zxlerRGfX1105/R3CdtIqNie/t5PC9nRy+t5PjFOAh6iKA6+bYz909wRiO44c+4QJW6xN/BrwQ+AXgqxHx34EfAD8LnEhN1k5izN/YQ4hYN8sjl3j+iAWvm/R5ljTmbtirNGqzBnxvJ4fv7eTwvZ0cvreTY5lRn6lSSvlldQ1LMIhPlFLmIuI51CUvzm0f9wPXUGdLvp0qYreOU+wQIva1dlxqrPWkdlxqrHbo8xhjjDFmdhnMJ0op+4C3tI8DRMRh1GUsdgP/PFqZlSFG269ux7Mi4gHni4jDgWdQC/3ClM5jjDHGmNllGj5xLnAo8KG2nMXIjC1ipZQbgauoY7ovW/D0xcA24P2llHsBImJrRGxvq96OfB5jjDHGmIUM5SXtuSMWeewngEuAe4A3jlvvUM36L6WOmV4eEWcC1wJPA06nRn+vnffa49rzN3NwQ95azmOMMcYYsxhDecknI2I38BXgbuDxwLOBPcDzSimLrTG2JgaZCNzs8ynAFdR/6IXUGQWXA6eVUm6f5nmMMcYYM7sM6BMfBg6nzp58JfDjwLupWx59YohaoxT9Yn3GGGOMMbOIemk8Y4wxxpiZxSJmjDHGGCPCImaMMcYYI8IiZowxxhgjwiJmjDHGGCPCImaMMcYYI8IiZowxxhgjYt2LWEQ8MiLeGxHfjog9EbErIi6NiIcqzrORGPeeRMTREfHiiPhIRHw9InZHxF0R8dmI+PWFe4DNEpP4fouIcyOitI8XD1nvemLIexsRPx0RfxkR32nn+k5EXBURz55E7dkZ8OftL7b7+K32c+GmiPiLiDhtUrVnJSJeEBFvi4jPRMT32//fD4x4Lv8eW4es6wVd275Q1wDHAH8FXAc8lbqFwdeAZ6xm9dyhzrORGOKeRMRLgD8CvkPdhPVfgEcAzwOOBP4SOKes52/CEZjE91tEPAr4X8Bm4CHAb5RS3j1k3euBIe9tRLwO+C/AbcDfUL+PfwR4EnB1KeV3Bv8HJGbAn7dvBn4HuB34b9T7+2PAc6jb7v1KKWUkEVmPRMSXgFOo+xZ+C9gOfLCU8sI1nse/x9YrpZR1+wF8AijAby94/A/a4++c5nk20scQ9wQ4A/glYNOCx4+lSlkBnq/+t67He7vg8wL4O+BG4PfaOV6s/neu53sLnNNe/0ng8EWe36r+t67He9v+7+8HbgGOWfDc6e08N6n/rVO+r6cDJ7X/x89q9+ADiq+PP0TfA+oCRi4cTmjfXN9Y5Bf94dR3F/cC26Zxno30MY17ArymXeNt6n/ver+3wAXAHPBMYMesitiAPxM2ATe11z5c/e/K8DHgvX1aO89fLfH894G71f9e4X0eScT8e2x9f6znHp0z2vGqUsrc/CdKKXcDnwMeDDx9SufZSEzjntzfjvvGOMd6ZNB7GxGPAy4BLiulfHrIQtchQ93bnwQeA3wU+F7rZ3pVRFwwiz1MjaHu7Q3AXuCp/3979xdiVRXFcfy7zIrojzKVBBEOmYXgS0USFaWgRtCDQfQQVE9RIIRkJASWD0EWCCP1VFRC9ZJB0ENUFP2lf4QUlJGFjQ+BiZX9QUuI1cPaly4X7zgzZ9/Zs8/5feCwuXef2eyzOHP2Yu9zzjWz8/orzOx6Iml4O0uPu0XjWMVqTsQuS+W+IfXfp/LSOWqnTUYaEzNbCNyZPr4xmzYqli22KY4vEMu8DzXvWvVyxfaqVP4M7CHuD9sOTAAfm9n7ZnZ+k45WKEts3f1XYAtxr+heM3vazB4zs5eBt4il4Hsy9LdrNI5VbGHpDjSwKJW/D6nvfb94jtppk1HHZDuwEnjd3d+cZRu1yhnbh4kbx69z92NNO9YCuWK7JJX3Eks9a4HPgKXADuBGYDexjNQV2c5bd58ws0ngOeDuvqofgF3ufmi2newwjWMVq3lG7GQslU2fyMvVTpvMOiZmdh+wmXii546cnWqJacXWzFYRs2A73P2TkfeqHaZ73p7St/+t7v6Ou//l7t8AtxBPtt3Q4WXKE5n2NcHMHgReAXYBy4AzgSuJ+/JeMrMnRtTHLtM4No/VnIj1MvxFQ+rPGdhv1O20yUhiYmYbgZ3AXmBNWqbomsax7VuS3Adszde16uU6b39L5X53/6q/Is089mZxV824h/XKElszWw08Drzm7ve7+353P+rue4gk9ydgs5ldnKHPXaJxrGI1J2LfpXLYmvfyVA5bM8/dTptkj4mZbQKeAr4mkrCDs+9e1XLE9qz09yuAv/te4urAI2mfZ9J3E417XI/c14QjQ+p7idoZ0+xXG+SK7c2pfHewwt2PAp8T49LlM+1gx2kcq1jN94j1/pHXm9mC/idFzOxs4FrgGPDpHLXTJlljYmZbiPvCvgTWufvhzP2tSY7Y/gM8O6TuCmIQ+4i4OHdp2TLXefsB8TTvcjM7zd2PD9SvTOVk8y5XI1dsT0/lsIcdet8PxlympnGsZqXfn9FkYwYvsANOJd5YvKxJO13ZMsZ2a9r/C2Cs9HHNhy1XbIe0vY2OvkcsZ2yBF9P+jw58v454Z9sRYHHp460ttsBtad+DwIUDdTel2B4Dzi19vIVivJop3iOmcaydW9t+4uhb4oWBa4gp2Gs8/aSDmY0TT0AdcPfx2bbTFTlia2Z3ETfk/gs8yYnvT5h0912jOYr5Kdd5O6TtbcTypH7iqNk1YQnx7qVLgA+JJbOlxH1MDtzu7rtHfkDzSKZrwgIiYVgL/Am8SiRlK4hlSwM2ufvOuTim+cDMNgAb0scLiKdy9xPnHcBhd38g7TuOxrH2KZ0JNt2Ai4Dnid+BOw4cIG4IHxvYb5y4gE42aadLW9PY8v/szFTbe6WPs8bYTtFuL+adnBHLGVtgjJhN+DG18wsxwF1d+hhrji0xq7OJWCb7g1gGPkS8r2196WMsENOTXScn+/bVONbCreoZMREREZGa1fzUpIiIiEjVlIiJiIiIFKJETERERKQQJWIiIiIihSgRExERESlEiZiIiIhIIUrERERERApRIiYiIiJSiBIxERERkUKUiImIiIgUokRMREREpBAlYiIiIiKFKBETERERKUSJmIiIiEghSsREREREClEiJiIiIlKIEjERERGRQv4DCQfNK3VAZe8AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 2 Axes>"
      ]
     },
     "metadata": {
      "image/png": {
       "height": 252,
       "width": 305
      },
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "V = df.FunctionSpace(mesh, 'P', 1)\n",
    "\n",
    "# The exact solution\n",
    "u_D = df.Expression('a0 + a1*x[0] + a2*x[1]', a0 = 1, \\\n",
    "                 a1 = 2, a2 = 3, degree=1)\n",
    "\n",
    "# Define boundary condition\n",
    "def boundary(x, on_boundary):\n",
    "    return on_boundary\n",
    "\n",
    "bc = df.DirichletBC(V, u_D, boundary)\n",
    "\n",
    "# Define variational problem\n",
    "u = df.TrialFunction(V)\n",
    "v = df.TestFunction(V)\n",
    "f = df.Constant(0.0)\n",
    "a = df.dot(df.grad(u), df.grad(v))*df.dx\n",
    "L = f*v*df.dx\n",
    "\n",
    "# Compute solution\n",
    "u = df.Function(V)\n",
    "df.solve(a == L, u, bc)\n",
    "\n",
    "# Plot solution and mesh\n",
    "c = df.plot(u)\n",
    "plt.colorbar(c)\n",
    "#plot(mesh)\n",
    "\n",
    "print(\"L2_error is = \", df.errornorm(u_D, u, 'L2'))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The $L_2$ error is really small as expected. Also, since FEniCS is a well tested program, of course, it passes the patch test. \n",
    "\n",
    "This concludes our discussion on Poisson equation which is of [__elliptic__](https://en.wikipedia.org/wiki/Poisson%27s_equation). The next topic is solving the dynamic heat (or diffusion) equation using FEniCS which is extension of the Poisson equation to the PDEs of [__parabolic__](https://en.wikipedia.org/wiki/Heat_equation) type. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Below, I am attaching the entire patch code, which could be save as a .py (python) file. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "L2_error is =  3.624438086253896e-16\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAmIAAAH4CAYAAADpQ4FeAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAAWJQAAFiUBSVIk8AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAgAElEQVR4nOy9eZQkZ3mv+by1dbda3WpJaENbS7qC5oIHG3FZDGYRRixeBmO459iDjLleYOCMWGfsC8YIz+FcGC8IsDHYgGWMz5nraw/32kYCscjsDGdkMx4uEgKJ1i61pJa6pe6qrq7Kb/6IiKrs7IjKLSLetzJ/zzl1oiozK+Lr6KzMJ3/vG99nKSWEEEIIIUT7zHgPQAghhBBiWpGICSGEEEI4IRETQgghhHBCIiaEEEII4YRETAghhBDCCYmYEEIIIYQTEjEhhBBCCCckYkIIIYQQTkjEhBBCCCGckIgJIYQQQjghERNCCCGEcEIiJoQQQgjhhERMCCGEEMIJiZgQQgghhBO1iJiZvcLMPmRmXzWzg2aWzOxTI+7rHDP7hJndbWZHzGyvmV1lZifXMVYhhBBCiCjM1bSf3wGeDDwK3AnsGWUnZnYR8A3gdOC/ATcBTwPeCLzYzJ6VUnqwlhELIYQQQjhTV2nyzcDjgJ3A/zzGfj5MJmFXpJRellL67ZTSpcD7gccD7xl7pEIIIYQQQbCUUr07NHsecD3w1ymlVw3xexcCtwB7gYtSSp2u+3YA9wAGnJ5SOlTnmIUQQgghPIjUrH9pvr2uW8IAUkqPAF8HTgCe0fbAhBBCCCGaIJKIPT7f3lxx/w/y7eNaGIsQQgghROPU1axfByfl2wMV9xe37+q3IzO7oeKuJ5FdULB3qJEJIYQQcdgNHEwpXeA5CDP7a0a8OK8PN6WU/qcG9huSSCLWD8u34zS1zc7NzZ1y2mmnnVLHgIQQQoi2uf/++1lZWfEeBsCebVvtKY+/eL62HX7/B0dZXKq3dz06kUSsSLxOqrh/Z8/jKkkpXVJ2u5ndcNpppz3lta997dptXz/LOLjFyh6+xvzcar9DNsrWhaOuxwfYNuc7hu0Ly67HB9gxd8R7CJw4t+R6/JOcjw+wc/aw6/FPnvM9PsDOmUXvIXDyjO952Dnr+7oMsMu2tHq8/Qd2cs0//RQf/ehHueeee/a2evAKHn/xPN/43GNr299PvuhuvvP/+b/et0kkEft+vq3qAbs431b1kI3ET+xLfP78fjI2x8KC56eP7NPG9i1+T86HmXOVoYP5U/WkLX4icJhZAHYt+L0JHmGGXfN+b4D78+2p834XLj8CnDL3qNvxiyOfNveI2xiKj0WnzPqdhxXgdMfjF5wy65sMPWZma2vHuu2us1o7lmiPSM361+fby8zsmHHl01c8C1gEvjXugQ7NwdHcu05YhRfe1mHnkY2j0OXlOZaXfb310JEF3+MvL3Bo2XcMB46096JXxcPL23yPf/QEHj56gusYHjy63fX4+1dOZP/Kia5juH9lh+vxAfav+p6Dfasnss95DPtXfV+XH+gs8UCnnQ+I3//R7laOI9qldREzs3kz25PPor9GSukW4DqyJsQ39Pzau4HtwCfrmENs3wnGl86bGVrGgBAyFkHIPDlwZKu7kD28vC2EkHny4NHtIYTMk/tXdrgL2f7VE0MImSf7V+dCCFmT7H94J/seVHvzJFLLM9fMXga8LP/xzHz7TDO7Ov/+gZTS2/LvzwZuBG4jk65uXk+2xNEHzewF+eOeDjyfrCT5jjrGC+syduntHebTuox9/vyZvj1jhYx5lisPHVlwLVUWMuZZrjxwZKtrqRIyIfMsVRYy5lmufPDodtdSZSFjnuXK+1d2uJYqIRMyz1JlIWOe5cr9q3OupcpCxpooV37/1t2171PEoK5E7MeBV+dfL8pvu7DrtlcMspM8FXsqcDWZgL0VuAj4IPDMuteZHCcZA6VjoHQMlI6B0jFQOlagdKyZdExlycmlFhFLKV2ZUrINvnZ3PXZv7209+7ojpfSalNJZKaWFlNL5KaU3ppT2lz1+XOqQsQhC5np89Y4B6h0D9Y6BesdAvWNQb++YypKTTaRmfTfGlTFQOgZKx0DpGCgdA6VjBRFkLIKQjYvKkpONRCynLhmLIGSux1c6BigdA6VjoHQMlI7B+OlYd1nytFMaKQ4JRyRiXdQhY6B0DJSOgdIxUDoGSscKIshYBCEblu6y5OzsKhedd0fdwxLOSMR6qFPGIgiZ6/GVjgFKx0DpGCgdA6VjMHw61l2WvPDcO50nFxdNIBEroS4ZA6VjoHQMlI6B0jFQOlYQQcYiCNkgdJclH3/h3mYGI1yRiFVQt4xFEDLX4ysdA5SOgdIxUDoGSsegfzq2GcqSK3TW/h11fK3Q8f4ntY5EbAPqlDFQOgZKx0DpGCgdA6VjBRFkLIKQldFbltyisuREIhHrQxMyFkHIXI+vdAxQOgZKx0DpGCgdg/J0TGXJ6WAqRayzuvESRr3ULWOgdAyUjoHSMVA6BkrHCiLIWAQhg81RlhT1MJUiBrCyNMfK0uB/cE3JWAQhcz2+0jFA6RgoHQOlY6B0DDIZ++dbzln7WWXJyWZqRazAW8ZA6RgoHQOlY6B0DJSOFUSQMU8h2/ujC9e+V1lyspl6EYPh0rEmZSyCkLkeX+kYoHQMlI6B0jGY3nTswIGT2L//MQDMzq5w8jk/bH0Mmx0z22tmqeLr3gH38asb7KP4Wh13rL7v/MFYWZpjbmv/+LeQsUtv7zCf1mXs8+fPcHDLcP1nvSwvz7lO2FfI2PYty35jWF5g+4Lf8QsZO2lLPQv2jkIhY7sWFv3GcPQEds0fdjt+IWOnzh9yG8P+lRM5Ze5Rt+MXMnba3CNuYyhk7JRZv/Owb/VETnc8fiFjp8y289r8o1svWvv+7HPuYH7hKA90jvKYGf8PipuMA8BVJbcP+mT6DvDuivt+CrgUuHaEcR2DRKyHIhnrJ2RNyxjgLmTeMga4C5mnjEEmZN4yBrgLmbeMAe5C5iljkAmZt4wB7kLWhox1lyUvuOCWte8f6CxxKB1t/PgTxMMppStH/eWU0nfIZOw4zOyb+bd/Nur+C1SarGCQUmVTZcqCCKXKCOVKT9Q7lo8hQKkyQrnSE/WOZUQoVTZZruwtS5533m2NHUuMhpk9CXgGcBfwmXH3JxHbgEF6x9qQsQhC5np89Y4B6h0D9Y6BesdgsnvHysqSYmS2mNmrzOztZvZGM3u+mc3WsN/X5tuPp5TUI9YG/XrHmixTFqh3TL1joN4xUO8YqHesYBJ7x6rKklPGHjO7oeyOlNIlQ+znTOCvem77kZm9JqX05VEGZmbbgFcBHeBjo+yjFyViA9IvHWs6GQOlY6B0rEDpmNIxUDoGk5WOqSxZK38BvIBMxrYDPwZ8FNgNXGtmTx5xv/8e2AVcm1KqZZZdJWJDslE61kYyBkrHQOkYKB0DpWOgdKxgEtKxzViWXE1Wa5l2NRnATUMmX8eRUuq92vG7wOvM7FHgrcCVwC+MsOvfzLcfHX10x6JEbAQ2SsfaSMZA6RgoHStQOqZ0DJSOweZPx1SWbIWP5NvnDPuLZvZvgZ8E7gSuqWtAErEx8JYx0JWVoCsrQVdWgq6sBF1ZWRBBxoYVMpUlW2Nfvh3lxaLWJv0CidiYVKVjbctYBCFzPb7SMUDpGCgdA6VjsPnSsc1YltykPDPf3jrML5nZVuBysib9j9c5IIlYTXjLGCgdA6VjoHQMlI6B0rGCCDI2iJCpLFkfZvZEMzul5PbzgT/Of/xU1+3zZrbHzC7q/Z0uXgmcDFxTV5N+gUSsRsrSMQ8ZiyBkrsdXOgYoHQOlY6B0DOKnYypL1s4rgbvN7Foz+7CZvc/M/ha4Cfg3ZP1df9D1+LOBG4EvbrDPokl/7Jn0e5GINYC3jIHSMVA6BkrHQOkYKB0riCBjZUKmsmTtXA98GrgA+GXgLcBzga8BrwZ+NqU08GX3ZvYE4NnU3KRfIBFriN50zEvGIgiZ6/GVjgFKx0DpGCgdg5jpmMqS9ZJS+nJK6ZdSSntSSrtSSvMppdNSSi9MKX0ypZR6Hr83pWQppd0V+7sxv//cOpv0CyRiDeMtY6B0DJSOgdIxUDoGSscKIsjY/tU5lSWFRKwNutMxTxmLIGSux1c6BigdA6VjoHQMYqRj37vl4rXvVZacTiRiLeItY6B0DJSOgdIxUDoGSscKPGXsrr0XrH2vsuR0IhFrmSId85axCELmenylY4DSMVA6BkrHwCcde/TATg7uPxWAmdkVTjznrlaPL2IwnSLWMdLyrOsQvGUMlI6B0jFQOgZKx0DpWEGbMnbX3vUm/TPOvpP5+aMjzcovNjdT/b+dlmexhdovgBiYlaU57p6BL513tPGFwqsoZMx7EXHvBcQB90XEPRcQh0zIvBcQB9wXEfdeQBxwX0TccwFxyITMewFxoPFFxLvLko/d/aNj7tu/OjfWAuJtcZSZWuX1KPfXtq/NwnQmYl2k5Vn3dOzumXnXZAyUjoHSMVA6BkrHQOlYQZPpWG9Z8szzbj/uMUrHpoOpF7GCCDL2+TPm3WUsgpC5Hl+9Y4B6x0C9Y6DeMWiud6ysLFmFZGyykYh14Z2O3bdtxl3GQOkYKB0DpWOgdAyUjhXULWMblSXLUDo2uUjESpCMKR0DpWMFEWQsgpB5onQsI4KM1SFkg5Qlq1js6G170tD/aAWe6VgUGQOlY6B0DJSOgdIxUDpWMK6MDVOWFJOPRKwPkjGlY6B0rCCCjEUQMk+UjmVEkLFRhWzYsqSYbCRiA+CVjkWSMVA6BkrHQOkYKB0DpWMFw8rYOGVJMZlIxIZAMqZ0DJSOFUSQsQhC5onSsYwIMjaokKksKXqRiA2JRzoWTcZA6RgoHQOlY6B0DJSOFQwiYypLil4kYiMiGVM6BkrHCiLIWAQh80TpWEYEGasSMpUlRRkSsTFoOx2LKGOgdAyUjoHSMVA6BkrHCspkTGVJUYZmh6uBNtesLGTshfcdPWZtymvPnOfwzk4rYyhDa1ZqzcoCrVmpNStBa1bC8WtWTmJZcjXN1iq9q8l3lRsPlIjVRJvpWFky9pJ7j3LCQf//TqVjSsdA6RgoHQOlYwX7Vk9UWVJU4v/OPWFEkLGVJV8ZUu+YescKIshYBCHzRL1jGd4ydvOte9a+V1lSdCMRa4C20rEqGdu53HGXMVA6BkrHQOkYKB0DpWMP3Hbu2veTUpYU9SARa5AIMuYtZErHlI4VRJCxCELmidKxjLZlbPHgDg4/dDIANrPK3NkP1r6IuNi8SMQapo10bCMZA9xlDJSOgdIxUDoGSsdg+tKxB/eup2Enn30Ps/PZRU2SMQESsdaIIGPeQqZ0TOlYQQQZiyBknigdy2hDxrrLkqeef8cx942zZqWYDCRiLdJ0OtZPxkDpGCgdA6Vja2MIIGMRhMyTSU/HesuSp5x7V+njJGPTy1SKmK2a6/EjyJi3kCkdUzpWEEHGIgiZJ0rHMpqQsaqyZBlKx6aTqRQxgJkl3396k+nYIDIGSsdA6RgoHVsbQwAZiyBknkxiOrZRWbIKydh0MbUiBpmMRRCyJhhGxryFTOmY0rGCCDIWQcg8UTqWUYeMDVqWLEPp2PQw1SJWEEHGmhCyQWUMlI6B0jFQOrY2hgAyFkHIPJmEdGyYsmQVkrHJRyKWM6np2LAy5i1kSseUjhVEkLEIQuaJ0rGMUWVslLJkGUrHJhv/GCQYM0szdLb6LZ5dyFidi4iXLRT+knuPcu2Z8xxcOF4+V5bmmNvqt3g3ZELmvYA44L6IuPcC4oDrIuKFjHkvIu69gDjgvoi49wLigOsi4oWMDbqI+DhlySoiythKmqlVllfS9OVD0/cvHoBJTMeGScZA6ViB0jGlY6B0DJSOFQyajtVRlizj0c6WWvYj4iAR24AIMlankA0rY6DeMVDvGKh3bG0MAWQsgpB5sll6x+oqS4rJRyLWh0lLx0aVMW8hUzqmdKwggoxFEDJPlI5lVMlYE2VJMblIxAYkgozVJWSjyBgoHQOlY6B0bG0MAWQsgpB5EjUda6osKSYTidgQTFI6No6MeQuZ0jGlYwURZCyCkHmidCyjW8ZUlhTDIBEbgQgyVoeQjSpjoHQMlI6B0rG1MQSQsQhC5kmUdOzeh85QWTIAZrbXzFLF171D7OcVZvYhM/uqmR3Mf/9TdY7V/910k1LImPdUF+NOczHs1BbdFDLmOdVFIWPeU114T3MBuE914TnNBWRC5j3NBeA+1YX3NBeA+1QXntNcHLz9rLXvVZZ05wBwVcntwzxBfwd4cv47dwJ7ahjXMUjExmQS5h0bR8ZA846B5h0DzTu2NgbNOzbV84490iViW87d1/rxxTE8nFK6csx9vJlMwH4IPBe4ftxB9aLSZA1MQu/YOGVKUO9YQYRSZYRypTcRSpURypWeTGPv2PLB7Rx56CQgK0ueeM59tS4gLtonpXR9SukHKaXU1DEkYjUSQcbGEbJxZQzUOwbqHQP1jq2NIYCMRRAyT9rsHesuS25/7P3MzmeVinHXrBQjs8XMXmVmbzezN5rZ882s/rUEx8T/XXPC2Oy9Y+OWKUG9YwXqHVPvGKh3DKand6y7LLnjvLuPu3//6okDL5E05ewxsxvK7kgpXTLEfs4E/qrnth+Z2WtSSl8eeXQ1o0SsITZzOlZHMgZKx0DpGCgdWxuD0rGJTsfKypJlKB1rjb8AXkAmY9uBHwM+CuwGrjWzJ/sN7Vj83yknmM2cjtWRjIHSsQKlY0rHQOkYTG46VlWWrGJS0rGVNFurYK+kWYCbhky+jiOl9O6em74LvM7MHgXeClwJ/MI4x6gLJWItsFnTsbqSMVA6BkrHQOnY2hiUjk1cOtavLFmG0jEXPpJvn+M6ii4kYi2xWa+srFvGvIVMV1bqysqCCDIWQcg8mZQrKwctS1YhGWuVYk4R3yd/F7WZgZmdY2afMLO7zexIPqvtVWZ28pD7+Rkzu87M7jSzRTO71cz+i5k9s66xehJBxoYVsjplDJSOgdIxUDq2NoYAMhZByDwZNx0btixZhtKx1ihc4lbXUXRRixWY2UXADcBrgG8D7yf7R74R+KaZnTrgft4H/CPwFOCzwAeAfwb+R+DrZvaqWsbr17IFbM50rAkZ8xYypWNKxwoiyFgEIfNkM6djo5Qlq5CMjY+ZPdHMTim5/Xzgj/MfP9V1+7yZ7cldpnXqehf6MHA6cEVK6UPFjWb2R2Sz0r4HeN1GOzCzM4G3AfcB/0NKaV/Xfc8HvgT8Hl0nbxxmljOj6Cw0Nkdb/zFssln562rg70az8mtWftCs/Gtj0Kz8m25W/nHLkmUUMjYJzfxOvBL4bTO7HvgR8AhwEfAzwFbgGuAPuh5/NnAjcBvZVZVrmNnLgJflP56Zb59pZlfn3z+QUnrbOIMdO5YxswuBy4C9wJ/03P0u4BBwuZn1+7h1fj6e/7tbwiCb2ZbsRJ427nh7KYTMi82WjtWdjIHSsQKlY0rHQOkYbK50rI6yZBVKx0bmeuDTwAXALwNvIVue6GvAq4GfTSkN+snzx/PfeTXwovy2C7tue8W4g63DAC7Nt9ellI55N04pPQJ8HTgBeEaf/fwAWAaeZmaP6b7DzJ4D7AC+UMN4j2Nm2UIImSfD9I41IWOg3jFQ7xiod2xtDAFkLIKQeTJI71idZcky1Ds2PCmlL6eUfimltCeltCulNJ9SOi2l9MKU0id7lytKKe1NKVlKaXfJvq7M76v6Ou53hqWOd//H59ubK+7/Qb593EY7SSntB34LOAP4npn9mZn9JzP7G+A64PPAa2sYbyURZCyCkA1CkzLmLWRKx5SOFUSQsQhC5knkdKyJsmQVkrHJpY53m5Py7YGK+4vbd/XbUUrpKjPbC3wC+I2uu34IXN1bsqyiamkEYE+/31Xv2OC9Y030jBWod0y9Y6DesbUxqHcsZO9Yk2XJMvavnshi8v2AJOqnjfiliJn6mo2Z/W/A3wJXkzXWbQcuIbsC86/N7P9oaIzHoXRssHSsqWQMlI4VKB1TOgZKxyBeOtZ0WVJMB3W80xeJ10kV9+/seVwpZvY84H3A36eU3pJSujWldDil9M9kyxDcBbw1vzhgQ1JKl5R9ATcN8g8qUO/YYL1jTcoYqHcM1DsG6h1bG0MAGYsgZJ7cv7KDex86vbWypJhs6niX/36+reoBuzjfVvWQFfxsvr2+946U0mGy+clmgJ8YdoDjEkHGIgjZRrQhY95CpnRM6VhBBBmLIGSeeKdjR+5Yv4i/jbKkmFzqeFcpxOkyM5vpvnLSzHYAzwIWgW/12c+WfFs1RUVxu0vDinrH+veONdkzVqDeMfWOgXrH1sag3jG33rHlLhHjnIdaP34UVtJMrVK+kqZv5cWx/8UppVvIrmrcDbyh5+53k/V5fTKldAg2nMH2q/n2N83s7O47zOwlZEK3BHxj3DGPg9KxjdOxppMxUDpWoHRM6RgoHYP207HVR7ax+nB+vJkOC499sNYFxMV0Udc7+uvJFtL8oJn913zaiS+Rzap/M/COrscWM9h+sWcff0s2T9gZwI1m9pdm9j4z+3vgM2RN/7+dUnqwpjGPjHrHNu4da0PGQL1joN4xUO/Y2hgCyFgEIWuD7rLk/Fn7sbwsOe6alWI6qeXdPE/Fnkp2tePTgbeSXfX4QeCZg8hTXtJ8KZm8fY+sQf+tZBPBXgO8KKX0gTrGWxcRZCyCkJXRpox5C5nSMaVjBRFkLIKQedJGOtZdltxy7v3H3S8ZE8NQ27tHSukOskW/+z1uL+tTWvTedxS4Kv/aFKh3rLp3rI2esQL1jql3DNQ7tjYG9Y411jtWVpYsY9g1K8X0Mn1dcQ2hdKw8HWsrGQOlYwVKx5SOgdIxaCYdqypLVqF0TPRDIlYj6h0r7x1rU8ZAvWOg3jFQ79jaGALIWAQhq4t+Zcky1DsmNkIi1gARZCyCkHXjIWPeQqZ0TOlYQQQZiyBkntSRjg1alqxCMibKkIg1hNKx49OxtmUMlI6B0jFQOrY2hgAyFkHIRmXYsmQZSsdELxKxhokgYxGErMBLxryFTOmY0rGCCDIWQcg8GTUdG6UsWYVkTBRIxFpA6dix6ZiHjIHSMVA6BkrH1sYQQMYiCNmgjFuWLEPpmACJWKtEkLEIQga+MuYtZErHlI4VRJCxCELmyaDpWB1lySokY9ONRKxllI6tp2NeMgZKx0DpGCgdWxtDABmLIGQbUWdZsgylY9OL/7uRA1bfB5mRmVk290lgAfeJYO/bRmuTvvZSyJjnRLCFjHlPBOs9CSzgPhGs5ySwkAmZ9ySwgPtEsN6TwALHTQTbRFmyivtXdmyqSWBX00ytIr+qRb+nh9klY3bJOZlSOkZanuXe2Xm3ZAyUjoHSMVA6tjYGpWPHpWNNliXLUDo2XUytiBV4yxiodwwIIWPeQqbeMfWOFUSQsQhC5kl371jTZckqJGPTwdSLGCgdWxtDABm77tQtbjIGSsdA6RgoHVsbQwAZ8xayBx86tbWyZBlKxyYfiVgX3jIGSsfu2zIbQsa8hUzpmNKxgggyFkHI3Ljr5LVv2yhLViEZm1wkYj0oHcvHMOUyBkrHQOkYKB1bG0MAGXMRsjvXRezoYw+2f/wu7l/ZwaHOFtcxiPqRiFXgLWMw3elYJBnzFjKlY0rHCiLIWAQha4v0yBY4kP97Zzpw1oFa1qwUohuJ2AYoHcvHMOUyBkrHQOkYKB1bG0MAGWtFyLrKkpx5AJtff+2RjIm6kIgNgLeMwfSmY9FkzFvIlI4pHSuIIGMRhKxRusqSnP3QcXcrHRN1IBEbEKVj+RimXMZA6RgoHQOlY2tjCCBjTQhZWVmyCsmYGAeJ2JB4yxhMZzoWUca8hUzpmNKxgggyFkHIamWDsmQZSsfEqEjERkDpWD6GCDJ2j5+MgdIxUDoGSsfWxhBAxmoTsj5lySokY2JY/N9FNjGzS8bqVr/1ImH61qwsZOyyB49ka1N2Mhm79qzm16asQmtWZmjNSq1ZCZOxZuUwZckyqtasnERW00ytHwC01qQYGqVj+RhaTMeOS8ZyGdtxyPccKB1TOgZKx9bGsJnTsSHLklUoHRODIBGrCW8Zg+nqHSuTsZc+sBRCxryFTL1j6h0riCBjEYRsaEYsS5ah3jHRD4lYjSgdy8cQQMbS8mwrY6jCW8ZA6RgoHQOlYzBcOjZuWbIKyZioQiLWAN4yBtOTjlXJ2M6jnRAy5i1kSseUjhVEkLEIQtaXmsqSZSgdE2VIxBpC6Vg+hgAyFkHIvIkgYxGEzBOlY/kYAsjYhkJWY1myCsmY6EYi1jDeMgbTkY5tJGNACBnzFjKlY0rHCiLIWAQh66WpsmQZSsfaxcwuN7OUf/36EL+3t+v3er/urWNs/h/Vp4BCxjynuihkzHuqiyanuSib2uKlDyxxzWO2cnB+Zk3GbGG1sTH0Y2VpznWaC8iEzHuaC8B9qgvvaS4A16kuChnznurCe5oLYH2qiwbLklXsXzlxKqa58MTMzgU+BDwKjGK/B4CrSm6v5T9OItYimnes+XnH+skYZOmYt4yB5h3TvGOadwyCzTvWQlmyjGmad6xtzMyAvwAeBP4v4G0j7ObhlNKVdY6rG5UmW0a9Y/kYGixV9itTAuody4lQqoxQrvREvWP5GJxLlQ89tKu1smQVKlU2whXApcBrgNFn+W2Q6RSxDsz4fQgG1DsGzfaODSJjoN4xUO8YqHesIIKMeQnZwt07139oqSxZhnrH6sPMngC8F/hASukrY+xqi5m9yszebmZvNLPnm1ltbx7TKWI5EWTMW8gmOR0bRsYiCJk3EWQsgpB5onQsH4ODjC3cvWPt+0fPXKx/EfEhkYyNh5nNAX8F3A68fczdnZnv6z1kvWJfAn5gZs8dc7/AlIsYZDIWQci8iSBjTQjZoDIGSsdA6RgoHSuIIGNtCdnMo/PMHczOeZrpsHxm1qsVQcamVMj2mNkNZV9D7ON3gZ8AfjWlNE4T5F8ALyCTse3AjwEfBXYD15rZk8fYN6Bm/TVmlqHj+NqrKyvzMTRwZeUgDfwFurIyQ1dW6spKmJ4rK7vLkkdPPwRz68G2choAACAASURBVK9Bx11Z6UDkKytXk9X6wWE1jR8KmNnTyFKwP0wpfXOcfaWU3t1z03eB15nZo8BbgSuBXxjnGFOfiHWjdCxjEtOxYZIxUDoGSsdA6VjBpKdj3WXJ5cc+UvoYpWOtclNK6ZKyr36/2FWSvBl4Z4Nj/Ei+fc64O5KIlRBBxryFbBJ7x0aRsQhC5k0EGYsgZJ6odywfQwMyVlWWLGOYNSub4rBn6WZzcCLwOOAJwFL3BKzAu/LH/Hl+W9ncYIOyL9+O/YTwf5UPSiFj3uVKzTtW77xjw5QpCzTvmOYdA807VjBp845tVJasYm3eMRGRI8DHK+57Clnf2NeA7wPjlC2fmW9vHWMfgESsL+odm7zesVFlDNQ7pt4x9Y7BZPWODVKWLCNC75g4nrwxv3QJIzO7kkzE/jKl9LGu2+eBi4CjKaVbum5/InBPSml/z37OB/44//FT445ZpckBUO9YRoRSZV3lymHLlAURSpXe5Ur1jql3rCBCqXKccuUwZckqvEuVohbOBm4Evthz+yuBu83sWjP7sJm9z8z+FrgJ+DfANcAfjHtwidgQRJAxbyGbpN6xcWQsgpB5E0HGIgiZJ+ody8cwooyNUpYsI0LvmGiE64FPAxcAvwy8BXguWWnz1cDPppTGNgOJ2JAoHcuIIGN1CNmoMgZKx0DpGCgdK4ggY8MK2ahlySokY7FJKV2ZUrLusmR++9789t09t385pfRLKaU9KaVdKaX5lNJpKaUXppQ+mVKqpV9HIjYiEWTMW8gmJR0bV8YiCJk3EWQsgpB5onQsH8OAMlZHWbIMpWNiWCRiY6B0LCOCjI0rZOPIGCgdA6VjoHSsIIKM9ROyusqSVUjGxKBIxGoggox5C9kkpGN1yFgEIfMmgoxFEDJPlI7lY9hAxuouS5ahdEwMgkSsJpSOZUSQsXGEbFwZA6VjoHQMlI4VRJCxXiFrqixZhWRMbIRErGYiyJi3kG32dKwuGYsgZN5EkLEIQuaJ0rF8DF0y1nRZsgylY6IK/1fqCUSz8mds5ln5R5n0tQzNyq9Z+UGz8hdEmZX/vBbKklVM2qz8nTRT64eNTpq+fGj6/sUtonRsc6djdSRjoHSsQOmY0jHwT8fmD821WpYsQ+mY6EYi1jDqHcuIIGOjCFldMgbqHQP1joF6xwq8ZGz7PeslysOnLbVSlqxCMiZAItYaEWTMW8g2azpWt4xFEDJvIshYBCHzZFrTse33rovYo2ceHmuJpDpQOiYkYi2idCwjgowNK2R1yhgoHQOlY6B0rKAtGZs/NMeWR7Lz3ZlJHD4961cbd83KOpCMTS9TKWLm1zsNxJAxbyHbjOlYEzIWQci8iSBjEYTMk2lJx7rLkounLZLmjr2YKIKMScimj6kUMYBZ34uHlI7lRJCxYYSsbhkDpWOgdAyUjhU0KWO9ZcnS4ysdEy0ztSIGmYxFEDJPlI7lYwggYxGEzJsIMhZByDyZ1HSsqixZOYYAMiYhmw6mWsQKIshYBCHzJoKMDSpkTcgYKB0DpWOgdKygThnrV5YsPb7SMdECErEcpWNKx9bGEEDGIgiZNxFkLIKQeTJJ6dggZcnKMQSQMQnZ5CIR6yGCjEUQMm8iyNggQtaUjIHSMVA6BkrHCsaRsWHLkqXHVzomGkIiVoLSMaVja2MIIGMRhMybCDIWQcg82czp2ChlycoxOMvY4VXf54GoH4nYBkSQsQhC5k0EGesnZE3KGCgdA6VjoHSsYFgZG6csWXr8AOmYmBz8P+oGp5CxVcfXnpll/wXEAddFxAsZ815EfKMFxOtaKLyKQsa8FxH3XEAcMiHzXkAccF9E3HsBccB1EfFCxvotIl5HWbJyDEdPYNf8+GK3mel0rNYPKJ2O/4f/tlEiNiBKx5SOQf90rOlkDJSOgdIxUDpW0C8dq7MsWXp8pWNiTCRiQ6DeMfWOrY0hgIxFEDJvIshYBCHzJHrvWN1lycoxSMbEiEjERiCCjEUQMm8iyFiVkLUhY6B0DJSOgdKxgl4Za7IsWXp8pWNiBCRiI6J0TOnY2hgCyFgEIfMmgoxFEDJPoqVjTZclK8cgGRNDIBEbkwgyFkHIvIkgY2VC1paMgdIxUDoGSscKHl7e1lpZsvT4SsfEgEjEakDpmNKxtTEEkLEIQuZNBBmLIGSeeKdjWw7PtlqWrEIyJvohEauRCDIWQci8iSBjvULWpoyB0jFQOgbTnY6dfN+Wte8PnrrcWlmyDKVjYiMkYjWjdEzp2NoYAshYBCHzJoKMRRAyTzzSsV371kXsodOP1LZm5ThIxkQZErGGiCBjEYTMmwgy1i1kbcsYKB0DpWMwXenYlsOznPBo9v/dmUkceMz6i2EEGZOQiW4kYg2idEzp2NoYAshYBCHzJoKMRRAyT9pIx3rLkp2esqTSMREJiVgLRJCxCELmTQQZK4TMQ8ZA6RgoHYPJT8d6y5JVRJAxCZmQiLWE0jGlY2tjCCBjEYTMmwgyFkHIPGkiHduoLFmG0rHx6CRb+3BTx1cn+X9ob5vaRMzMzjGzT5jZ3WZ2xMz2mtlVZnbyCPv6KTP7OzO7J9/XPWZ2nZm9tK7xehFBxiIImTcRZGxmacZNxkDpGCgdg8lLx/qVJauIIGObWcjE6NQiYmZ2EXAD8Brg28D7gVuBNwLfNLNTh9jX7wBfAZ4DfBb4Q+AfgJOB59Uy3ubf4zZE6ZjSsbUxBJCxCELmTQQZiyBkntSVjg1alixD6ZjwoK5Xnw8DpwNXpJQ+VNxoZn8EvBl4D/C6fjsxs1cC/zvwBeDlKaVHeu6fr2m8zOYisur42jO7BKuOE1AXMtZxPQfG6la/+X0gE7LOgt8YZpZmuJ8ZrjsVLnvwCPNpXcauecxWDs4330GQlmexhdXGj1NFIWNzW1fcxlDI2MKC3xgOHVlg+xa/T0mFjG1f8BvDgSNbOWnLaJ9Uhy1LVvHw8jZ2LfhMAAvrMrZrvt3VAIQPY7/Cm9mFwGXAXuBPeu5+F3AIuNzMtvfZzwzwPuAw8Mu9EgaQUjo67nh7mXVPhpSOKR3LuD/NuyVjoHSsQOnY5k3HRi1LlqF0TLRFHR+1L82316WUjnnHyGXq68AJwDP67OcngQuAa4CHzOxnzOy3zOyNZvbMGsZZyexyDCHzRL1jGRFk7As7trnJGKh3DNQ7Bpuzd2ycsmQVEWRMQjbZ1CFij8+3N1fc/4N8+7g++/l3+fY+4J+BfwTeC1wFfMPMvmxmp40z0H5EkLEIQuaJ0jG4b342hIxFEDJvIshYBCHzZNB0rK6yZBlKxzYfZvY+M/uimd1hZotmtt/M/sXM3jVk37qZ2X8ws2+Z2SNmdjjfzxVmVsuLZB0idlK+PVBxf3H7rj77OT3fvg7YBvw0sAN4EvA5sub9/zLIgMzshrIvYE+/31U6pnSsIISM5T97yBgoHQOlY7A50rE6y5JVRJCxRc/G5s3Fm4HtwOeBDwB/DawAVwL/ambnDrifvwQ+Tlax+8/AnwML+T7/s5mN/UbRxqtLMch+fxXFK74Br0gp/b/5z//dzH6BLHF7rpk9M6X0zQbGeexglv0b+cG/md+7kR9wbeYvZMyjmf+++Vm+sHMbP31wkXnab+AvKGTMu5nfs5EfMiHzbuQH3Jv5vRv5gdJm/ibKkmUUMubZzC8GYmdK6bgnipm9B3g78B+B12+0AzN7GXA58CPgaSmlB/Lb54G/AX4ReDVw9TgDrePVvEi8Tqq4f2fP46p4KN/e2iVhAKSUFslSMYCn9RtQSumSsi/gpn6/243SMaVjBV7pWCFj3skYKB0DpWMQMx1rsixZhXc6JjamTMJy/ibfXjzAbl6eb/+wkLB830eBd+Y//i+jjXCdOkTs+/m2qges+MdW9ZD17ufhivsLUWv92R9BxiIImSfT3DsWTcYiCJk3EWQsgpB50t071kZZsowIvWNiaH4u3/7rAI89M9/eWnJfcdtTzKxf69WG1CFi1+fby/IpKNYwsx3As4BF4Ft99vMVsvrtxWZW9hf+pHy7d/Shjo7SMaVjBdMuY6B0DJSOQZx0rK2yZBWSsbiY2dvM7Eoze7+ZfZVsrtJ/JbsYsB9FCnZByX0Xdn3ft/98I8YWsZTSLcB1wG7gDT13v5usWe6TKaVDkNVWzWxPPht/934eIGuEOwn43e77zOyFwIvIypufHXfM4xBBxiIImSfTmo5FlLEIQuZNBBmLIGRebFu01suSZSgdq5U9G1x0NyxvI5vT9E3As8kc4rKU0v0D/O4/5tu3mNkpxY1mNkfmNwVDL+XYTV2vIK8HvgF80MxeANwIPB14PllJ8h1djz07v/82Mnnr5i35773DzJ5DtlzS+cAvAKvAb6SUqkqXraFZ+TUrf0Hbs/KXNvDff4RrTtvSagN/N5qVX7Pyg9+s/Kc/sL7gyoO7VlorS1bhPSt/63Rq/kBU8+fKlNKZAGZ2Btl8pe8F/sXMfjal9M99fv3/BF4FvAT4npn9Pdmk8z8NXEQ2PdfFZH4yMrW8cuep2FPJrhx4OvBWskF+EHhmSunBAfezL//99wPnAleQTRj7GeCnUkoDTV/RFkrHlI5B++nYcclYSrz0/iNuyRgoHStQOtZ+OnbaA+vnfN9jVmpbs3IclI6NzU0bXHQ3Eiml+1JKnyZbCehU4JMD/E4H+HmyVO1esiso/wNwJ1m6VrjNvlHHBTVOX5FSuoNs0e9+j9vL+pQWZffvJ0vG3lLX2JpE6ZjSsYI207HjkrFcxj570jYe3nAxsWZROqZ0DNpLx7YtGjsOZR8AVi3x4Cnr53ycNSvrYurSsU1ASuk2M/se8ONm9pjuqyErHr8C/GH+tYaZbQN+nKwH/r+PMyafWsYEonRM6Ri0m46VJWMvPrDIrkOtHL4SpWMZSseaT8e6y5L7T15hteeUKx0TFTw2347zqfFyYCvwN+Ougy0RqxFdWakrKwsiyNjMku+fdwQZ8xYyXVnZ7JWVvWXJKrxlDHRlZZvkFwSeWXL7TD6h6+nAN1JKD+W3l15EmN+3s+S2f0fWa/Yo8Hvjjtf/Y+MEoln5NSs/tDcrf1mZ8sUHFvnsSds4uDRDZ6tv7xhoVn7Nyl//rPwblSXL2GhW/rbQrPyt8WLg983sK8AtZL1cZwDPJZt24l7gN7oev9FFhJ83s0Xgu8AjwBOBlwJHgJenlMrmGBsKJWINoXRM6VhBG+lYVTK2c7XDzNKM0jGlY8BkpWP9ypJVKB2bCr4A/BlZU/7Lgf+VbDmi/WTTTjwxpfS9Aff1t2TrXr+KrHf9x4CP5fv43Ea/OChKxBpG6ZjSMWgnHdswGZvNZEzpmNKxSUnHBi1LlqF0bLJJKX2X4+c13ejxe6m4iDCl9PvA79czsnKUiLWA0jGlYwVNp2MbJWOA0jGUjhVs5nRs2LJkFUrHRAQkYi0SQcYiCJkn03BlZT8ZgxiN/BGEzJsIMhZByIZl1LJkGbqyUngjEWsZpWNKxwoiyFgEIfNE6VhGBBkbRsjGKUtW4S1joHRsWplKEbNVfxmJIGPe5yCCjHkLWZPp2CAyBkrHQOkYbJ50rK6yZBlKx4QHUyliBd4ionRM6VhBBBmLIGSeKB3LiCBjGwlZnWXJKrxlDJSOTRNTLWIQIxmKIGPe5yCCjHkLWVPp2KAyBkrHQOkYxE7HmihLlqF0bDBSsrW/21q+kv8H47aZehEr8BYRpWNKxwoiyFgEIfNE6VhGBBnrFrImy5JVeMsYKB2bdCRiXURIhiLImPc5iCBj3kLWRDo2jIyB0jFQOgax0rE2ypJlKB0TTSIRK8FbRJSOKR0riCBjEYTME6VjGRFk7NT710WsybJkFd4yBrDUln2K1pCIVRAhGYogY97nIIKMeQtZ3enYsDIGSsdA6Rj4pmPbF42TDmfPw7bKkmVESMfEZCER64O3iCgdUzpWEEHGIgiZJ0rHMjxk7KyH1p979+/qcNBz3TRipGNiMpCIDUCEZCiCjHmfgwgy5i1kdaZjo8gYKB0DpWPQfjp21v71593dJ2fP0boWEB8VpWOiDiRiQ+AtIkrHlI4VRJCxCELmidKxjDZkrLcsed/J68/PcdasrAvJmBgHidiQREiGIsiY9zmIIGPeQlZXOjaqjIHSMVA6Bs2nY71lydWS//IIMiYhE6MgERsRbxFROqZ0rCCCjEUQMk+UjmU0JWNlZcnS4ysdE5sQidgYREiGIsiY9zmIIGPeQlZHOjaOjIHSMVA6BvWnYxuVJSvHEEDGJGRiUCRiNeAtIkrHlI4VRJCxCELmidKxjLpkbJCyZOnxlY6JTYJErCYiJEMRZMz7HESQMW8hGzcdG1fGQOkYKB2DetKxQcuSlWMIIGMSMrERErGa8RYRpWNKxwoiyFgEIfNE6VjGqDI2Slmy9PhKxxrD0vrfeh1flrz/Re0jEWuACMlQBBnzPgcRZMxbyMZJx+qQMVA6BkrHYLR0bNSyZOUYAsjYpAqZGB2JWIN4i4jSMaVjBRFkLIKQeaJ0LGMYGRu3LFl6fKVjIhgSsYaJkAxFkDHvcxBBxryFbNR0rC4ZA6VjoHQMBkvH6ipLVo4hgIxJyARIxFrDW0SUjikdK4ggYxGEzBOlYxkbyVjdZcnS4ysdEwGQiLVIhGQogox5n4MIMuYtZKOkY3XKGCgdA6VjUJ2ONVGWrBxDABmTkE0vEjEHvEVE6ZjSsYIIMhZByDxROpbRLWNNlyVLj690TDgxlSJmHZg94nuNbIRkKIKMeZ+DCDLmLWTDpmN1yxgoHQOlY7CejrVRlqwcQwAZk5BNF1MpYgXeMgb+IqJ0TOlYQQQZiyBknigdyzjjgfX/h6bLkmUoHRNtMtUiBpmMeQtZhGQogox5n4MIMuYtZMOkY03IGCgdg+lOx048AifnfwdtlSWriCBjErLJZ+pFrMBbxsBfRJSOKR0riCBjEYTMk2lNx859eP37e3bAwRVfGVI6JppGItaF0rF8DO4i4n8OIsiYt5ANmo41JWOgdAymLx3rFrE7dmXbOtasHJcIMiYhm0wkYiV4yxj4i4jSMaVjBRFkLIKQeTIt6VhvWfLuncfeH0HGvIVsaWXe9fiifiRiFSgdy8fgLiL+5yCCjHkL2SDpWJMyBkrHYPLTsd6y5ErJ6VY6Foxka68PdXyR/D98to1ErA/eMgb+IqJ0TOlYQQQZiyBknkxyOlZWlqwigoxJyEQdSMQGQOlYPgZ3EfE/BxFkzFvI+qVjTcsYKB2DyUvH+pUly1A6JiYBidgQeMsY+IuI0jGlYwURZCyCkHkySenYIGXJKiLImIRMjIpEbEiUjuVjcBcR/3MQQca8hWyjdKwNGQOlYzAZ6dgwZckylI6JzYpEbES8ZQz8RUTpmNKxgggyFkHIPNnM6dgoZckqIsiYhEwMg0RsDJSO5WNwFxH/cxBBxryFrCoda0vGQOkYbM50bJyyZBlKx4SZvc/Mvmhmd5jZopntN7N/MbN3mdmpbe9nIyRiNeAtY+AvIkrHlI4VRJCxCELmyWZLx8YtS1YRQcYkZG68GdgOfB74APDXwApwJfCvZnZuy/upxP+j04RQyNjqFr83wkJEVh0nX55dhlXH150I52BmGTqu5yAv8Wz1+4BQyFhnYX0MhYz99MFF5lmXsc+etI2Ds/WL08zSDJ2tfusUFjJmC6tuY1hZmmNu64rb8SETsoWF6jHUWZYso5Cx7Vv8PiUdWl5g+4Lzp7TpY2dK6biP52b2HuDtwH8EXt/ifipRIlYzSseUjoHSsYLedKzNZAyUjkH8dKzusmQVSsemizJ5yvmbfHtxm/vZCIlYA6h3LB+Du4j4n4MIMuYtZL29Y23LGKh3DOL2jjVVlixDvWMC+Ll8+69B9qPSZJPMHkmupUrIRMS7VAn+5UrvUiX4lys9S5WQCVlRqmy7TAnrMuZdrvQuVQKu5cpCxhYWVhovS1Zx6MiCe6kSULmymj1mdkPZHSmlS4bZkZm9DTgROAl4KvBsMnl6r8d+ypCINYx6x/IxqHdMvWMc2zvmIWOg3jGI0zt27sPrY2iyLFmGesfqwTr1tkFY/X+abwPO6Pr5s8CvppTud9rPcag02RLepUrwL9Opd0y9YwWFkHmUKUG9YxCjd+ych9afi02XJauIUKpUufI4bkopXVL2NeyOUkpnppQMOBN4OXAh8C9m9hSP/ZQhEWsR9Y7lY3AXEf9zEEHGvIWs6B3zkjFQ7xj49Y7tOJI45Uj2/arRWlmyDPWOTT4ppftSSp8GLgNOBT7puZ9uplLEZvwqAoDSMVA6BkrHCiLIWAQh88QjHTv/kfXXwbu2w+HVetasHIcIMiYha5aU0m3A94AnmtljvPcDUypiAHNLzsmU0rFsDO4i4n8OIsiYt5DNLBv3pzk3GQOlY9BuOnb+wfXXv9t2rj//IshYBCETjfLYfDtuLFPLfqZWxCCTsQhC5o23iCgdUzpWEEHGIgiZJ22kY71lyTtPPPa5N+qalXUSQcYkZKNhZnvM7MyS22fyiVhPB76RUnoov30+/52LxtnPqOiqSTIhW9nqeFWjrqzMxqArK3VlJZmMfXHbNl6w2O7VlN3oyspmr6zsLUuuzJa/9vWblb9pdGXlpuXFwO+b2VeAW4AHya54fC5Zk/29wG90Pf5s4EbgNmD3GPsZCYlYTpGMeQuZ5h3Ltt5C5i1j4C9knjK2by6GjIHmHYP65x2rKkuW0T3vmBead2zT8QXgz4BnAU8GdgGHgJuBvwI+mFLa3+J+NkQi1oPSsRjJkNIxpWMRZAyUjkG96Vi/smQVSsfUOzYoKaXvAm8Y4vF7geOeiMPuZ1SmukesCvWO5WNQ75j7OZj23rFCxrx6xgrUO1Zf79igZcky1DsGy6u+zwNRPxKxDYggY95CFuGqwggy5n0OIsiYl5BFkTHQlZUw/pWVw5Qlq4ggY95CJiYHiVgflI7lYwggYxGEzJNpTseiyVgEIfNk1HRs1LJkGUrHxKQgERuQCDLmLWQRkqEIMuZ9DiLImIeQRZIxUDoGw6dj45Qlq4ggYxIyMQ4SsSFQOpaPIYCMRRAyT6Y1HYsoYxGEzJNh0rE6ypJlKB1zJK2/HtXxhf9bXOtIxEYggox5C1mEZCiCjHmfgwgy1raQRZMxUDoG/dOxOsuSVUSQsakVMjEyErERUTqWjyGAjEUQMk+mMR2LKmMRhMyTjdKxJsqSZSgdE5sNidiYRJAxbyGLkAxFkDHvcxBBxtoUsioZ2+X8NxlBxiIIWS9NlSWriCBjEjIxCBKxGlA6lo8hgIxFEDJPpi0dK5Oxyw4fDiFjEYTMk+50rI2yZBlKx8RmQCJWIxFkzFvIIiRDEWTM+xxEkLG2hGwjGZtZ9l0yLIKMRRCytsqSVUSQMQmZqEIiVjNKx/IxBJCxCELmyTSlY1UytrPTCSFjEYTMk/MOrH/fRlmyDKVjIioSsYaIIGPeQhYhGYogY97nIIKMtSFk/WQsgpB54pWO7Tza4dTl7LVopcWyZBURZExCJrqRiDWI0rF8DAFkLIKQeTIt6dhGMgaEkLEIQtYmuw+tX8l617YZlo7O17Jm5TgoHRORkIi1QAQZ8xayCMlQBBnzPgcRZKxpIRtExiIImSdtpmPdIrZ3+/q/21vGQOmYiIFErCWUjuVjCCBjEYTMk2lIx/rJGCgdg+bTsd6y5O0nHPvvHXXNyjpROia8qe1VwMzOMbNPmNndZnbEzPaa2VVmdvIY+7zczFL+9et1jdWTCDLmLWQRkqEIMuZ9DiLIWJNCNqiMRRAyT5pMx3rLkisz5efaW8ZA6Zjwo5ZXADO7CLgBeA3wbeD9wK3AG4FvmtmpI+zzXOBDwKN1jPGYfa/6iojSsXwMAWQsgpB5Munp2CAyBkrHoJl0rKosWYbSsQzJ2PRR11/+h4HTgStSSi9LKf12SulSMiF7PPCeYXZmZgb8BfAg8JGaxngMEWTI+/hKx/IxuIuI/zmIIGNNCdkwMhZByDypMx3rV5aswlvGQOnYMFhn/TWsji/zW6XMjbH/6s3sQuAyYC/wJz13vws4BFxuZtuH2O0VwKVkCduhcce4Ed4yFEEIvWUM/EVE6dhkp2ODyhgoHYN60rFBy5JlKB3L2CwyJsajjr/2S/PtdSmlY17VUkqPAF8HTgCeMcjOzOwJwHuBD6SUvlLD+PoSQYa8j690LB+Du4j4n4MIMtaEkA0rYxGEzJNx07FhypJVeMsYKB0TzVPHX/rj8+3NFff/IN8+rt+OzGwO+CvgduDt4w9tOLxlKIIQessY+IuI0rHJTceGkTFQOgajpWOjliXLUDqWIRmbXOp4Zp2Ubw9U3F/cvmuAff0u8BPAs1NKi6MOyMxuqLhrT7/fLURoZavfC/DcUnI9fiFjq1scx5CLyOpWtyEwuwyrjq99Ec7BzDJ0XM9B9hxc3VrfB4RCxl6wuMg86zJ23QkncHDmeGEoZKyz4PchZWZphs5Wv+aZQsZsYXWgx49TlqxiZWmOua0rY+9nHJaX51hY8BvDoSMLLPueAtEAbXzUKv4CN3wVM7OnkaVgf5hS+mbjo+qDdzKldCwfg9Ix93MwienYsMkYKB2DwdOxOsqSZSgdE5NIHX8hReJ1UsX9O3sedxxdJcmbgXeOO6CU0iVlX8BNw+wnggx5H1+9Y/kY3EXE/xxEkLE6hWxUGYsgZJ706x2rsyxZhbeMgX/vmJgc6vgL+X6+reoBuzjfVvWQAZyY//4TgKWuSVwT2ZWXAH+e33bV2CMeEm8ZiiCE3jIG/iKidGzy0rFRZAyUjkF1OtZEWbIMpWNiUqjjGXR9vr3MzGa6r5w0sx3As4BF4Fsb7OMI8PGK+55C1jf2NTLpcylbqndMvWNrY1Dv2ET1jg3bM1ag3rHy3rGmypJVqHdMbHbGFrGU0i1mdh3ZXGJvIJsNv+DdwHbgoymlQwBmNg9cBBxNKd2S72MRKF3CyMyu5BSYpwAAIABJREFUJBOxv0wpfWzc8Y6LtwxFEMLZI8lVxiCTEW8ZA38h85Yx8BcyTxmDTMi8ZQxwFzJbWG2lLFlGkYx5ClmRjEnIxLDU9VfyemAf8EEz+69m9p/M7EvAm8lKku/oeuzZwI3AF2s6dutEKBV6H1+9Y/kY3Mt0/ucgQqmyjnLlqGVKUO8YZDK2++D6z02WJavwLlWCesfE8NTyl5snW08FrgaeDryVLPX6IPDMlNKDdRwnGt4yFEEIvWUM/EVEvWOT0zs2joyBesd2L66nQW2UJctQ75jYbNT2l5JSuiOl9JqU0lkppYWU0vkppTemlPb3PG5vSslSSrsH3O+V+ePdy5JlRJAh7+MrHcvH4C4i/ucggoyNK2R1yFgEIWubnUc7nHo0L0sCt83Ntz6GbrxlDJSOicHwzbInCG8ZiiCE3jIG/iKidGwy0rFxZQymLx27YGm9Yf+urbOszFhtC4iPitKxFkjrr3t1fG084+hkIhGrkQgy5H18pWP5GNxFxP8cRJCxcYSsLhmLIGRt0F2W/NG2dQEbd83KOvCWMVA6JqqRiDWAtwxFEEJvGQN/EVE6tvnTsTpkDCY/HestS96+9XjxiiBj3kI28emYGAmJWENEkCHv4ysdy8fgLiL+5yCCjI0qZHXKWAQha4KysmQZSscyJGOiG4lYw3jLUAQh9JYx8BcRpWObOx2rS8ZgMtOxqrJkFRFkzFvIlI6JAolYC0SQIe/jKx3Lx+AuIv7nIIKMjSJkdctYBCGrg0HKkmUoHcuQjAmJWIt4y1AEIfSWMfAXEaVjmzcdq1PGYDLSsUHLklVEkDFvIVM6Vi9mdqqZ/bqZfdrMfmhmi2Z2wMy+Zma/ZmYDPenN7Fe7176u+Frtv6eN0f98y0RYosh7mSatWZmPQWtWbso1K8dZDqmMzb5m5bBlyTLK1qxsG61ZOVG8EvhT4B6y9bBvB84AXg58DHiJmb0ypdTvj+47ZEs1lvFTwKXAteMOdipFzFZhdrHD6ja/QNBbhiIIodas1JqVsDnXrKxbxmBzrlk5almyimLNSi+0ZuXEcDPw88BnUkprT2gzezvwbeAXyaTs7zbaSUrpO2Qydhxm9s382z8bd7BTXZqcXfRbJBdilAq9j6/esXwM7mU6/3MQoVQ5TLmy7jIlbL7esXHLkmWodyxDpcrRSSl9KaX0D90Slt9+L/CR/Mfnjbp/M3sS8AzgLuAzo+6nYKpFDDIZiyBk3sf3HoO3jIG/iKh3bPP1jjUhY7B5esfqKEtWEUHGvIVMvWONUPy5jhM5vjbffjylNHaEO/UiVhBBxrxlyPv4SsfyMbiLiP85iCBjgwpZkzIWQciqqLssWYbSsYwplLE9ZnZD2dc4OzWzOeBX8h8/O+I+tgGvAjpk/WZjIxHrQulYDCH0ljHwFxGlY5srHWtKxiBuOtZEWbKKCDLmLWRKx2rhvcCTgGtSSp8bcR//HtgFXJtSuqOOQel/tYQIjfygKytBV1bqysrNc2VlEw38BRGvrGyyLFmGrqzMWF2NlZ9Yp94PbZY9xW5KKV1S317BzK4A3grcBFw+xq5+M99+dOxB5cT6Hw2E0jGlY2tjUDrmfg42SzrWZDIGcdKxNsqSVSgdE8NiZm8APgB8D3h+Smn/iPv5t8BPAncC19Q1PolYHyLImLcMeR9fvWP5GNxFxP8cRJCxfkLWhox5C9mFj6z/W5ouS5ah3jExKGb2JuCPge+SSdi9Y+yu1ib9AonYACgdiyGE3jIG/iKidGxzpGNNyxj4pmO7l9stS1YRQcYkZHExs98C3k82F9jzU0r7xtjXVrKSZgf4eD0jzJCIDUEEGfOWIe/jKx3Lx+AuIv7nIIKMbSRkbclY20K2c7XDKavZv2EVuNPmWz1+L0rHRBlm9k6y5vwbgBeklB7Y4LHzZrbHzC7aYJevBE4ma/SvpUm/QM+eISlkzLuZ37uRHzQrv/eM9JqVP/6s/E028HfT5qz8u4+sp2F3zs+yYsZMLqSjLpNUB5qVXxSY2auB3yP7rPBV4Aqz494v9qaUrs6/Pxu4EbgN2F2x26JJf+yZ9HuRiI2Irqz0F0JdWZmPQVdWhr6ysk0Zg+avrOwuS+7dcuxbyDhrVtaBrqwUORfk21ngTRWP+TJw9SA7M7MnAM+m5ib9ApUmx0C9YzHKpd6lSvAv06l3LHbvWBtlyoImS5W9Zck7Fo7/LD/orPxNEqFUqXKlHymlK1NK1ufreV2P35vftrtifzfm959bZ5N+gUSsBiLImLcMeR9fvWP5GNxFxP8cRJCxMiFrW8aaELKysmTlGALIWAQhE6IfErGaUDoWQwi9ZQz8RUTpWNx0rE0Zg/rTsY3KkqXHVzqmdEz0RSJWMxFkzFuGvI+vdCwfg7uI+J+DCDLWK2QeMlaHkA1SlqwcQwAZiyBkQpQhEWsApWMxhNBbxsBfRJSOxUzH2pYxGD8dG6YsWXp8pWNKx0QpErEGiSBj3jLkfXylY/kY3EXE/xxEkLFuIfOSsVGFbNiyZOUYAshYBCGbFCytv87W8WX+n59bRyLWMErHYgiht4yBv4goHYuXjnnIGAyfjo1Tliw9vtIxpWNiDYlYS0SQMW8Z8j6+0rF8DO4i4n8OIshYIWSeMjaokI1blqwcQwAZiyBkYrqRiLWI0rEYQugtY+AvIkrHYqVjXjIGg6VjdZUlS4+vdEzp2JQjEXMggox5y5D38ZWO5WNwFxH/cxBBxmaXzF3GqoSs7rJk5RgCyFgEIRPTx1SK2MxqYm7RbwkMUDpWHN97DN4yBv4ionQsTjrmKWNQno41VZYsPb7SMaVjU8hUiliBt4yB0rFiDJ4oHcvH4C4i/ucggow9uDLvLmPdQtZkWbJyDAFkLIKQielgqkUMMhnzFjKlYzGE0FvGwF9ElI7FSMe8ZQwyIWurLFl6fKVjSsemhKkXsQJvGQOlY8UYPFE6lo8hgIx5n4MIMvZPsye4ytgFh9srS1YRQcYiCJmYXCRiXSgdy/CWoQhC6C1j4C8iSsf807F9M3OuMnb+SvtlyTKUjikdm2QkYiV4yxgoHSvG4InSsXwMAWTM+xxMo4zt6HQ4pbNelryL+doXER+WCDLmLWSdFb1tTxr6H61A6ViGtwxFEEJvGQN/EVE65puOecjY7qNH176/a25urSwZQcYiCJkQdSER64O3jIHSsWIMnigdy8cQQMa8z8G0yFh3WfK2uWNLYuOsWVkXEWRMQibqQAXnAShkbGWb3x9dIWOr2/xefOaWEitb/V58CxnzHMPskcTqFt83oNklWN3qePxcRFYXHMfgfA4KGeu0fA72zczxT5zA81YPM8+6jF13wgkcnKnvtaG3LHnHXPlbxcyy0Vnw+4BSyFhnq9+H1bQ8iy34f2B3o1PzB/WOfwWibZSIDYHSMaVjoHRsbQxKx1zSsTaSsaqyZBlKx5SOifGQiA2JescyvGUoghB6yxj4i4h6x3x6x5qWsY3KklVEkLEIQibEsEjERsRbxkDpWDEGT5SO5WMIIGPe5yCCjL3o0PgyNmhZsgylY0rHxPBIxMZA6ViGtwxFEEJvGQN/EVE61n461itj2xhfxoYpS1YRQcYiCJkQgyARqwFvGQOlY8UYPFE6lo8hgIx5n4MIMrZrxNeEUcqSZSgdUzomBkMiVhNKxzK8ZSiCEHrLGPiLiNKxdtOxMhl74crwMjZOWbKKCDIWQciEqEIiVjPeMgZKx4oxeKJ0LB9DABnzPgcRZGx2aTAZqqMsWYbSMaVjohqJWAMoHcvwlqEIQugtY+AvIkrH2kvHqmRsR1odSMbqKktWEUHGIgiZEN1IxBrEW8ZA6VgxBk+UjuVjCCBj3ucggoxVCVkTZckylI4pHRPHIhFrGKVjGd4yFEEIvWUM/EVE6Vg76dhGMgaUylhTZckqIshYBCETQiLWEt4yBkrHijF4onQsH0MAGfM+BxFkrFvImi5LlqF0TOmYkIi1itKxDG8ZiiCE3jIG/iKidKz5dKyfjEEmZG2VJauIIGMRhExMJ1O56LetOovI4qrrAuKQCZn3AuLgu4C39yLmhYx5LiJeiIj3IuLeC4iD/yLiTS0g3rtQeCFjn587gUcsex26YHE9DWujLFlGIWPei4h7LyAObKpFxC3V++Ha/D+jts7UJmKzSyv9H9QgSscyvJMppWP5GJSOuZ+DJtOxfsnYeZ31/rC2ypJVKB1TOjZtTK2IQSZjEYTMmwgy5i1D3sdX71g+hgAy5n0O2paxsztHOYX1suTdK/MDzzvWFOodU+/YNDHVIlYQQca8hUzpWAwh9JYx8BcRpWPNpWNlMva81cW1+++29bKkt4yB0jFQOjYNSMRylI5lRJAxbxnyPr7SsXwMAWTM+xy0IWPd3D4zf8zPG8071hZKx5SOTToSsR4iyJi3kCkdiyGE3jIG/iKidKz+dGxb6nByWmW+5L6qZ5y3jIHSMVA6NqlIxEpQOpYRQca8Zcj7+ErH8jEEkDHvczCOjG1LHR6/eoQXrhzi5SuP8tTOkdLHPXt18ZipLbpROpaPwVnGWNXbdj/M7FQz+3Uz+7SZ/dDMFs3sgJl9zcx+zcwGPolm9goz+5CZfdXMDppZMrNP1TneqZy+YlBml1ZY3ep3igoZ85zqopAx76kuvKe5AN+pNmaPJNdpLiATEe9pLsB/qgvvaS5gsKkutqUO53WOcl5a4fQKuepw/Kfx3qktepldMla3+n44mFk292kuANepLsSGvBL4U+Ae4HrgduAM4OXAx4CXmNkrU0qDPIl+B3gy8ChwJ7Cn7sFKxPpQJGPeQqZ5x/xlyFsINe9YPgbNO1Y579ig8nWfzXL7zDx32BxHbIbTOysbzjPWS5GMeQqZ5h0TG3Az8PPAZ1JKa/9BZvZ24NvAL5JJ2d8NsK83kwnYD4HnkoldrUjEBkTpmNKx4vigdMw7GVI6tp6ObZkfTb66GWTS1zKUjikdi0hK6UsVt99rZh8B3gM8jwFELKW0Jl7W0ETHErEhUDqWoXTMXwiVjuVjmNJ0bBsdzuUo57HCaSujyVcv48gYKB1TOrZpKC4Y9m0E70IiNgJKx5SOFccHpWPeydC0pGPHyJfVI1+9jCpjoHQMlI7VyB4zu6HsjpTSJaPu1MzmgF/Jf/zsqPupG4nYiCgdy1A65i+ESsfyMUxgOjaQfCXYxyy3M8+dzLE4P97f47gyBkrHlI6F5b3Ak4BrUkqf8x5MgURsTJSOKR0rjg9Kx5SOjX8ORpGvI13XPg5zZWUV48gYKB2D6UnHrFPvVEeW7eqmcZKv0v2aXQG8FbgJuLzOfY+LRKwGlI5lKB3zF0KlY/kYNlk6Nq58lVF1ZeWg1CFjoHRM6Zg/ZvYG4APA94AXpJT2Ow/pGCRiNaJ0TOlYcXxQOqZ0bONz0IR89TJuOjaujIHSMZiedCwiZvYm4P3Ad8kkbJ/zkI5DIlYzSscylI75C6HSsXwMgdKxNuSrjHHSsbpkDJSOKR1rFzP7LbK+sO8AL0wpPeA8pFIkYg2hdEzpWHF8UDo2zenYNjqcM3uUczornDbbnnz1Mk46VoeMgdIxUDrWFmb2TuD3gBuAyzYqR5rZPHARcDSldEtLQ1xDItYgSscylI75C6HSsXwMLaVja/I16ytfZYyajtUpY6B0TOlYc5jZq8kkbBX4KnBFyWSse1NKV+ffnw3cCNwG7O7Z18uAl+U/nplvn2lmxe8+kFJ62zjjlYi1gNIxpWPF8UHp2KSmYwPLV2eWO1fnuX2+PfnqZdR0rC4ZA6VjoHSsQS7It7PAmyoe82Xg6gH29ePAq3tuuzD/gkzeJGKbAaVjGUrH/IVQ6Vg+hhrSsWHl667VLvlazd4lvJdJ8pYxUDqmdKxeUkpXAlcO8fi9QOkL4rD7GoXa3hHN7Bwz+4SZ3W1mR8xsr5ldZWYnD/j7p5rZr5vZp83sh2a2aGYHzOxrZvZrZkNOER2UQsi8mFtcXUvIvJhd7NQ678woFELmeXzvMRRC5jqGJefjL68nZIOyjQ4Xzx7h+QuH+Lltj/ITC0eOk7BOgntXZ/l/lrfyD0sn8pXl7dy6ulCagHmfg5nl9YRsUPbNzPFPsyesrRVTyNiOirUu+1EImSeFkLkdf2lmLSET00Ut8YyZXQR8Azgd+G9kE6Y9DXgj8GIze1ZK6cE+u3kl8KfAPWSrm98OnEG2QvrHgJeY2StTSmO/e9hqhxoddGiUjmUoHVM6BpsjHRsr+Rrk+AHOwbDpWJ3JGCgdWxuD0rGpoy4T+DCZhF2RUvpQcaOZ/RHwZrKVzl/XZx83Az8PfCaltPYsNLO3A98GfpFMyvqulj4IM4ezj4CdE/yua1fvmHrHiuODesei9Y41LV+lY3A+B8P2jtUtY6DeMVDv2LQx9rufmV0IXAbsBf6k5+53AYeAy81s+0b7SSl9KaX0D90Slt9+L/CR/MfnjTveXgoh82J2aSVEudKbCKVK71Kh9/FnjyT3cuXskm+pbqt1eDz1lR1HwfscwHClyrrLlJDJmHe5cmbZQpQrxeRTRxxzab69rkSiHjGzr5OJ2jOAL454jOJvvBFjUTqmdKxA6dj0pWNbrcM580c5e36Fx8y1k3wNwmZKx5pIxkDpGCgdmwbqeOd/fL69ueL+H5CJ2OMYQcTMbA74lfzHzw49uiGYObzsLmOg3jH1jvkL4aT3jg0qX/evzHJHak++etlMvWNNyhiodyxq75h1Uq1VFev4X0TUNnW845+Ubw9U3F/cvmvE/b8XeBJwTUrpc4P8gpndUHHXnn6/q3RM6ViBtwxFEMJJSseGka+7js5z18ocy6nr+Rd0zco2GDQda0rGQOkYZDJmK/5XmIp6aePdvnjWDP3sNbMrgLeSXYV5eZ2D6ofSMaVjEEOGvIVwM6djY8tX9xgCrVnpxSDpWNMyBkrHxGRRx7t8kXidVHH/zp7HDYSZvQH4APA9shXTK9eJ6iWldEnFPm8AnjLofpSOKR0r8JahCEK4WdKxOuXruOM7rlm5NoZNkI41KWOgdExMFnW8w38/3z6u4v6L821VD9lxmNmbgPcD3yWTsH2jD298lI4pHYMYMuQthFHTsSblq3QMSsf6pmNtyBgoHRObnzre2a/Pt5eZ2UzPHGA7gGcBi8C3BtmZmf0WWV/Yd4AXppQeqGGMY6N0TOlYgbcMRRDCCOnY9uUOZ53Ynnz1onSsfzrWtIyB0jGx+Rn7XT2ldIuZXUd2ZeQbgA913f1uYDvw0ZTSIQAzmwcuAo6mlG7p3peZvZNsxfQbgMuGKUe2hdIxpWMQQ4a8hdAjHds60+GxW45y9tYVTl1oX77KUDq2cTrWloyB0jGxOanr3fz1ZEscfdDMXgDcCDwdeD5ZSfIdXY89O7//NmB3caOZvZpMwlaBrwJXmB33Ar83pXR1TWMeGaVjSscKvGUoghA2nY5FlK9elI5tnI61IWOgdExsTmp5J89TsaeSidSLgZeSrRn5QeDdAyZbF+TbWeBNFY/5MnD1eKOtD6VjSscghgx5C2Hd6dig8vXA8ix3HZnnnqV1+fJeJslbxsBfyLxlDJSOic1Dbe/gKaU7gNcM8Li9rE9p0X37lcCVdY2nLZSOKR0r8JahCEI4Tjo2jnwdM4Zga1a6jCFoOtaWjIHSMbF58Hv3njCUjikdgxgy5C2Ew6RjdcnXcWMIkAwpHStPx9qWMVA6JmIjEasRpWNKxwq8ZSiCEFalY03JV+kYlI65n4OydKxNGQOlYyI2ErEGUDqmdAxiyJC3EBbp2Py21Jp8HTeGAMmQ0rHj0zEPGQOlY7WT0tp7Tl37mzYkYg2hdEzpWIG3DHkJ4dbZDmedsMJjt61wytZ25asM72RI6djx6VjbMgZKx0Q8JGINo3RM6RhMTzoWTb56iZAMKR07Nh3zkjFQOiZiIBFrAaVjSscKJjEdG1S+Hlya5e7Dc9y7OMfivLOYKx1zPwfd6ZiHjIHSMREDiViLKB1TOgaTkY6NIl/LnfVzHnXNytbHoHRsLR3zlDFQOib8kIi1jNIxpWMFmy0dG1e+yoiwZqV3MqR0bD0d27fgI2OgdEz4MZ0ittrBFpdJ2/xe+ZSOKR2D+OlYE/LVi9KxfAxKx5hZ9pcxUDom2mU6RSwngoyB0jFQOhYpHWtDvspQOqZ0DDIZe4A5/mnOR8ZA6Zhol6kWMchkDHAXMm8ZA6Vj05yObZnrcOaJK5y5Y4VTtnVKH9OEfPWidCwfg9IxHljxlzFQOiaaZ+pFrEDpmNIxmK50LIp8laF0TOkYZDL2ZU7gueYjY6B0TDSPRKwLpWNKxwomNR0bVL72H57l3kfnuOtoe/LVi9KxfAxTno7dzxxfTv4yBjHSMTF5SMRKUDqmdAwmJx0bVr7ue3SO5dX1Y/7/7Z17sJxlfcc/v3MLSciFBEIqAiExGgY7oFQQqSheaNXWwQudaQsybbFlar3iaMdqxV6m2BYniG2lpUorzlgvY3UcqyhDFbW0o1ZblJuEYJFATq4kOYfkXH794302OezZPWd33+fd5333/X5mdjbZd3n22Sfn7H74vLvvO0L60ySpjmXXqYUsuYwxwailkTEoRx2ziGcTEuVAItYG1THVsQZVrGN55avVHFLLGKiO1bmOjRN2U3p6GYO0dUwMFhKxRVAdUx2DatSx2PLV6vEh7aE2VMfqXcfKImNQjjpWBmzWj75PxRqvbkjEOkB1THWsQdnqWNHy1W4OqWUMVMfqWsfKJmOgOibyIRHrAtUx1TFIX8eWjMxy6vIZ1q3pn3w1ozoW5qA6lmQNyiRjoDom8iER6xLVMdWxBv2sY0tGZlm3aoaTV09zwvI08tUK1THVMUizBmWUMVAdE90jEesR1THVMSi2jnUqX3sPDrFjYrRv8tWM6liYg+pY39egpYxNTfC10TQyBqpjonskYjlQHVMdaxCrjnUjX4/vH2Hn/hGm5spXQhFQHVMdg/6vwTwZs3LIGKiOic6QiEVAdUx1DHqvY7nlaw6pZUh1LMxBdayva9BOxm5nGfvHEr4mqI4lwcxeD7wIOAc4G1gBfNLdL+9ynO3A6W02P+7u6/PMs4FELBKqY6pjDTqpYzHlq5kyyFBqIVQdC3OoUR1rJWMv9QluP5JexkB1rM+8l0zADgKPAFtyjLUf2Nri9oM5xnwKErHIqI6pjkHrOlakfLUitQyVQQhVx+pVxxaSsQMMM5t0DVTH+sjbyQTsJ2Rl7I4cY+1z92tjTKodErECUB1THWuwdHqGE0+e7Zt8NVMGGUothKpjYQ41qWNtZYxlHDiSXsZAdaxo3P2oeJmV/xydErECUR2rZx0bG3XWrZlh3doZVq/ov3y1IrUMlUEIVcfqU8cWkzEguZBJxirDEjO7HDgNOAT8D/BNd5+J9QASsYJRHatHHSujfDVTBhlKLYSqY2EONahjC8oYwwwdSS9joDrWhi1m9r1WG9z93D7PZT3wiabbHjKz33L3b8R4AIlYn1AdG7w61ql87XtiiJ27hxnfO8zUdHjxXVpfGSqDEKqO1aOOdSJjkF7IJGOl5ePAncCPgAPARuAPgN8F/s3MLnD3H+Z9EIlYH1Edq34dyyNfcynbOStTzSG1jIHq2KDXscVkDFAdy8OsH31vizUecG+C8jUPd/9A0013A1eb2UHgGuBa4DV5H0cilgDVsWrVsVjy1Uzqc1ZCehkqgxCqjg1+HetUxiC9kFVSxurHR8lE7KIYg9VTxKZnYGISli1NNgXVsXLXsaLkqxWqY+mFUHUszGGA61gnMgaqY6Ijdobr5TEGq6eINUgsY6A6BuWpY0Mrh/omX82ojpVDCFXHBruOdSNjkF7IJGOl5YJwvS3GYPUWMchkDFTHalrHxsacE0+a5aR1s6xa1fpFryj5aoXqWHohVB0LcxjQOtapjIHqWB0ws1FgEzDl7g/Ouf0sYIe772m6/+nAR8Jfb40xB4lYA9Wx2tSxsslXM6pj5RBC1bHBrWPdyhikFzLJWOeY2aXApeGvjfNBXmBmt4Q/73L3d4Y/nwLcAzwMbJgzzGXAH5rZHcBDZN+a3AS8CjgO+DLw1zHmKxGbi+oYMJh1rBP5cod9e43x8SF27xpiciTtr4fqWHohVB0LcxjAOtaNjEF56hjlP1B8GTgHuLLpto3hApl0vZOFuQN4FvAcsl2Ry4F9wLfIjiv2CXePYscSsVaojg1EHetFvqamjr3KjUyV85yV/Sa1DJVBCFXHBrOO9SJjkFbIrI91vqqEc0Ne2+F9t9NCb8PBWqMcsHUxJGLtUB0DqlfH8spXK8pwzkrVsfRCqDoW5jBgdaxbGYP0dUwMFhKxxVAdK30dK0K+mklxzspmVMfKIYSqY4NXx3qVMZCQifxIxDpBdQwoVx3rh3y1QnWsHDKUWghVx8IcBqiO9SJjoDom8iMR6wbVsaR1bGyJs3Y9nLj+CCtPaH2fIuSrGdWxjNQyVAYhVB0brDqWR8ZAQiZ6QyLWLapjQP/q2DH58qTy1QrVsXLIUGohVB0LcxiQOtarjIHqmOgNiVivqI4VVsc6kq9Z2L8Hdj1mjO8f6Zt8NaM6lpFahsoghKpjg1PH8soY1EjIZmePBYpY49UMiVgeVMeAOHWsW/na/ThMH5WvGYYp5zkr+4nqWHohVB0LcxiAOpZHxkB1THSORCwGqmM91bF88jWfspyzUnVMdUx1bDDqWAwZAwmZWBiJWCxUx4DF61hs+Wom1Tkr56I6Vg4ZSi2EqmNhDhWvY3llDFTHxMJIxGKjOjavjhUtX61QHVMdazw+qI6pjuVbg1gyBhIyMR+JWBGojjF2HKz5uSnWnjbcN/lqRnUsQ3UsvRCqjoU5VLiOxZAxUB0T85GIFUnN6tjYcbDmFGPtKUOsPLH1G07R8tUK1THVscbjg+qY6ljvaxBTxkBCJjIkYkUz4HWsM/ly9o87u8aH+yZfzaiOZagrOhVqAAAPo0lEQVSOpRdC1bEwh4rWsVgyBqpjIkMi1i8GqI51I1+7f+bsedSZPgIQjg9T0nNW9gPVsYzUMlQGIVQdq24diy1jICGrMxKxflLhOta7fM2nTOesTIXqWDlkKLUQqo6FOVSwjsWUMVAdqzP1FLGZmbSPX5E6FlO+mkl5zsoGqmOqY43HB9Ux1bHu16AIGQMJWd2op4gBPjGJpZShktaxIuWrFapjqmNQDhlKLYSqY2EOFatjsWUMVMfqRm1FDDIZA9ILWeI6NsYUazaN9U2+mlEdUx1rkFqGyiCEqmPVq2NFyRhIyOpArUWsQR3r2OhSY+1pw6w5bZiV61q/UBQpX61QHVMdg3LIUGohVB0Lc6hQHStCxqACdWx29mjUiDVe3ZCIBepQx8ooX82ojqmONUgtQ2UQQtWxatWxImUMSi5komckYk0MWh3rWL4en2XPT6fZ838zTB8Ot5f4nJVFozqWoTqWXghVx8IcKlLHipIxCEKW9iVBFIBErAVVr2N55GsuZTtnZQpUx1THGo8PqmOqY52tQZEyZom/9C/iIxFbgCrVsVjy1Uzqc1aC6hiojkE5ZCi1EKqOhTlUoI4VKWNisJCILUKZ61hR8tUK1THVMVAdazw+qI6pji2+BpIx0QkSsQ4pSx0bXbusb/LVjOqY6lgD1bH0Qqg6FuZQ8jomGROLIRHrglR1bHTZEGvPGGPNxjFWrh9teZ+i5KsVqmOqY6A61nh8UB1THVt4DSRjYiEkYj3QjzpWNvlqRnVMdayB6lh6IVQdC3MocR2TjIl2SMR6pIg61rF8PTrFnoeOsGf7EaaHEr7qoToGqmOgOtZ4fFAdUx1rvwaSMdEKiVhO8taxnuQrvOBnlPOclf1GdUx1DMohQ6mFUHUszKGkdUwyJpqRiEWg2zqWX75aUIJzVqqOqY6B6ljj8UF1THWs9RpIxsRcJGIRWaiOFSJfzSQ4Z2UzqmOqYw1Ux9ILoepYmEMJ69igyJjPzDJ78FDU8eqGRCwyc+tYX+SrFapjqmOojjVILUNlEELVsXLWsUGRMZEPiVhkxpYPsWbzMtZuXsrKpy1peZ9C5KsZ1TFAdQxUx6AcMpRaCFXHwhxKVsckY0IiFoHSyFcrVMdUx1Ada5BahsoghKpj5atjkrF6IxHrkY7l65HD7H5gkr2Pev/kqxnVMUB1DFTHoBwylFoIVcfCHEpUxyRj9UUi1gXdyteeByeZnnzqBw+TnyZJdQxQHQPVsdQyVAYhVB0rVx2TjNWTaK+CZvZ0M/uYmT1qZofNbLuZbTWzE1KME4ux5UOsP+d4zrrsJM696mmc8aLV8yTMZ519P32SB2/fy3dv3sE9n9/FzrsPzZMwOPZh/mRMTB4rZImwySNHC1kqGkKWiuEnp48WslQ0hCwlwy1+R/rJyJN+VIhSziElw4f9aCFLNocnj9WhZHNI+5JwdA0aMjYV/kmWmvNSJlhB+t/XqhHDJyzjt83sLjM7YGYTZvbfZvYWM4tix1H+t9zMNgHfAdYBXwDuBc4D3gr8spld6O67+zXOYvjUwm+AMcrXgo+f6JyVT0F1THUM1bEGqmOqY1CeOjZ+XPsyJjojok/8E3AFsBP4F+AQ8DLgBuAiM7vM3XP9n0ysd4C/JXuyb3H3Gxs3mtmHgLcDfw5c3cdxFqVx3JOh45cDxctXy/H6cM7KBdFnxwB9dgz02TEohwylFkJ9dizMoQSfHdvDCN84br6M7dAnijolt0+Y2aVkEvYQcJ677wq3jwKfBl4HXAnckmeiuV/1zGwjcAmwHfibps3vJ7PHK8xseT/G6YaxFSOcfOZotN2OveATk+XYXZmYMuyqLMPuypSMTM4k3105PDlbit2VqR8/9RxS76qEcuyqTL27cs+T83dTbrSptJOqABF94rXh+vqGhAG4+xTwvvDXN+edbwy1fkm4vs3dn/IK6u4HzOzbZAvyfOD2PoyzKCc9eyVLVo2w8tTWmTd2+eoE1THVMVAda6A6pjoGqmOQydidQ8t44VhWxkRHxPKJ9eF6W4ttjduea2ar3X1fr5ON8Wr/rHB9f5vtD5A94Wey8BOONU5HtJMwAHdYtmaUZeePcur5K/M+VPdY4t+2MvyyJ12DmRLM4Qhein+H1BMg+e9DKf4dElOKNaj5HCRhXRHLJxoV7IwW2zbO+fMW4K5uJjiXGCK2Klzvb7O9cfvqPo2DmX2vzaazx8fHuemmmxYbQgghhCgl4+PjABsSTwOAQzzBXVNfiToesKXd+7i7n9vBMLF84kvArwPvMLNPufseADMbAT4w5365jurQj/0fDY/P+6GDGOMMTU9Pz+zYseOHOeci5rMlXN+bdBaDida2OLS2xaG1LY6zgeNTTwK4d5YZDrA39rgbYg/YRKc+8SngcuAVwI/N7IvABNm3JjeRlbXNkO/YIjFErGGWq9psX9l0v6LHaWvMDcPu0KhFF2hti0NrWxxa2+LQ2hbHAnt9+oq7/2bqObQhik+4+6yZvZrskBdXhMsU2WExrgQ+QiZiO/NMNoaI3Reun9lm++Zw3W5fbexxhBBCCFFfovmEu08D14fLUcxsKXAOMAn8qLdpZsT4WtId4foSM3vKeGa2AriQbKKLfZAt1jhCCCGEqC/98IkrgOOAT4fDWfRMbhFz9weB28j26b6pafMHgOXAP7v7IcgOhGZmW8JRb3seRwghhBCimVheErbNO3SCmT0PuA44CPxJ3vnG+rD+75PtM/2wmb0UuAc4H7iYLP390Zz7nhK2P8z8D+R1M44QQgghRCtiecnXzGwSuBs4AJwFvBI4DLzW3VsdY6wrohwxMdjnL5Ad5v984BqybxR8GLig0/NDxhpHCCGEEPUlok98FlhB9u3JdwA/D9wMnOXuX40xV8t5rkohhBBCCNEj6c4hIoQQQghRcyRiQgghhBCJkIgJIYQQQiRCIiaEEEIIkQiJmBBCCCFEIiRiQgghhBCJkIgJIYQQQiSi8iJmZk83s4+Z2aNmdtjMtpvZVjM7IcU4g0TeNTGztWZ2lZl93sx+YmaTZrbfzL5lZr/TfA6wOlHEz5uZXWFmHi5XxZxvlYi5tmb2QjP7nJntCGPtMLPbzOyVRcy97ER8vX1VWMdHwuvCNjP7jJldUNTcy4qZvd7MbjSzO83sifD7e2uPY+l9rIJU+oCu4bxQ3wHWAV8A7gXOIzuFwX3AhZ0cPTfWOINEjDUxs6uBvwN2kJ2E9afAycBrgVXA54DLvMo/hD1QxM+bmZ0K/C8wDBwPvNHdb4457yoQc23N7L3AnwK7gC+R/RyfCDwHuMPd3xX9CZSYiK+3HwTeBewG/pVsfZ8BvJrstHtvcPeeRKSKmNkPgLPJzlv4CLAF+KS7X97lOHofqyruXtkL8FXAgTc33f6hcPtH+znOIF1irAnwEuBXgaGm29eTSZkDr0v9XKu4tk3/nQFfBx4E/iqMcVXq51nltQUuC/f/GrCixfbR1M+1imsbfvdngMeAdU3bLg7jbEv9XPu8rhcDm8Pv8YvDGtya4t9Hl0Q/A6kn0PPEYWP44XqoxRv9CrL/uzgELO/HOIN06ceaAO8Jj3Fj6udb9bUF3grMAhcB19ZVxCK+JgwB28J9T0r9vMpwibi254dxvtBm+xPAgdTPN+E69yRieh+r9qXKn9F5Sbi+zd1n525w9wPAt4FlwPP7NM4g0Y81mQrX0znGqCJR19bMzgSuA25w92/GnGgFibW2LwDOAL4M7A2fZ3q3mb21jp9hCsRa2weAI8B5Znbi3A1mdhGZNHw9yozrhd7HKkyVRexZ4fr+NtsfCNfP7NM4g0Sha2JmI8Abwl+/0ssYFSba2oZ1/ATZbt735J9a5Ym1ts8L148D3yf7fNh1wFbgO2b2DTM7Kc9EK0iUtXX3PcC7yT4r+mMz+3sz+wsz+zRwG9mu4N+LMN+6ofexCjOSegI5WBWu97fZ3rh9dZ/GGSSKXpPrgGcDX3b3r/Y4RlWJubZ/TPbB8V9098m8ExsAYq3tunB9NdmunpcB/wmcDlwP/BLwGbLdSHUh2s+tu281s+3Ax4A3ztn0E+AWd9/Z6yRrjN7HKkyVi9hiWLjO+428WOMMEj2viZm9BbiG7Bs9V8Sc1IDQ0dqa2XlkFex6d/+Pwmc1GHT6czs85/6vd/fb3f2gu/8IeA3ZN9teVOPdlK3o+DXBzN4FfBa4BdgELAfOJftc3ifN7C8LmmOd0ftYiamyiDUMf1Wb7Sub7lf0OINEIWtiZm8CbgB+DFwcdlPUjdxrO2eX5P3A++JNrfLE+rndG663ufsP524I5bFRcc/reobVJcramtmLgQ8CX3T3d7j7NnefcPfvk0nuz4BrzGxjhDnXCb2PVZgqi9h94brdPu/N4brdPvPY4wwS0dfEzN4GfAS4m0zCHut9epUmxtoeH/77M4En5xzE1YH3h/v8Q7hta+4ZV4fYrwn72mxviNrSDuc1CMRa218J13c0b3D3CeC/yN6XntPtBGuO3scqTJU/I9b4Rb7EzIbmflPEzFYAFwKTwF19GmeQiLomZvZuss+F/QB4ubvvijzfKhFjbQ8D/9hm23PJ3sS+RfbiXKfdlrF+br9J9m3ezWY25u5HmrY/O1xvzz/lyhBrbZeE63Zfdmjc3rzmYmH0PlZlUh8/I8+FLg5gB4ySHbF4U55x6nKJuLbvC/f/LrAm9fMqwyXW2rYZ+1pqehyxmGsL3Bru/2dNt7+c7Jht+4DVqZ9v1dYW+LVw38eAU5q2vSKs7SSwNvXzTbTGL2aB44jpfWwwL4N2iqN7yA4YeDFZgn2Bh1M6mNkGsm9APezuG3odpy7EWFszu5LsA7kzwI20/nzCdne/pZhnUU5i/dy2Gftast2TOsVRvteEdWTHXnoGcCfZLrPTyT7H5MBvuPtnCn9CJSLSa8IQmTC8DDgAfJ5Mys4k221pwNvc/YZ+PKcyYGaXApeGv64n+1buNrKfO4Bd7v7OcN8N6H1s8EhtgnkvwKnAx8nOA3cEeJjsA+Frmu63gewFdHuecep0ybu2HKszC13+PfXzrOLaLjBuY81rWcRiri2whqwmPBTG2U32Bvf81M+xymtLVnXeRrab7Amy3cA7yY7Xdknq55hgTRd7ndw+5756HxvAS6WLmBBCCCFElanytyaFEEIIISqNREwIIYQQIhESMSGEEEKIREjEhBBCCCESIRETQgghhEiEREwIIYQQIhESMSGEEEKIREjEhBBCCCESIRETQgghhEiEREwIIYQQIhESMSGEEEKIREjEhBBCCCESIRETQgghhEiEREwIIYQQIhESMSGEEEKIREjEhBBCCCESIRETQgghhEjE/wN7TslyIHHijQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 2 Axes>"
      ]
     },
     "metadata": {
      "image/png": {
       "height": 252,
       "width": 305
      },
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# import the libraries\n",
    "\n",
    "import dolfin as df\n",
    "import matplotlib.pyplot as plt\n",
    "import mshr \n",
    "import numpy as np # numpy library for arrays etc.\n",
    "\n",
    "# define an editor to write the mesh\n",
    "%matplotlib inline\n",
    "mesh = df.Mesh()\n",
    "editor = df.MeshEditor()\n",
    "editor.open(mesh, 'triangle', 2, 2)\n",
    "editor.init_vertices(5)\n",
    "editor.init_cells(4)\n",
    "\n",
    "# add vertices\n",
    "editor.add_vertex(0, np.array([0.0, 0.0]))\n",
    "editor.add_vertex(1, np.array([1.0, 0.0]))\n",
    "editor.add_vertex(2, np.array([0.75, 0.25]))\n",
    "editor.add_vertex(3, np.array([1.0, 1.0]))\n",
    "editor.add_vertex(4, np.array([0.0, 1.0]))\n",
    "\n",
    "# add cells\n",
    "editor.add_cell(0, np.array([0, 1, 2], dtype=np.uintp))\n",
    "editor.add_cell(1, np.array([2, 1, 3], dtype=np.uintp))\n",
    "editor.add_cell(2, np.array([2, 3, 4], dtype=np.uintp))\n",
    "editor.add_cell(3, np.array([2, 4, 0], dtype=np.uintp))\n",
    "\n",
    "# final formalities\n",
    "editor.close() # close the editor\n",
    "mesh.order # order the mesh\n",
    "df.plot(mesh) # plot\n",
    "\n",
    "# create a BVP on FEniCS and solve\n",
    "\n",
    "V = df.FunctionSpace(mesh, 'P', 1)\n",
    "\n",
    "# The exact solution\n",
    "u_D = df.Expression('a0 + a1*x[0] + a2*x[1]', a0 = 1, \\\n",
    "                 a1 = 2, a2 = 3, degree=1)\n",
    "\n",
    "# Define boundary condition\n",
    "def boundary(x, on_boundary):\n",
    "    return on_boundary\n",
    "\n",
    "bc = df.DirichletBC(V, u_D, boundary)\n",
    "\n",
    "# Define variational problem\n",
    "u = df.TrialFunction(V)\n",
    "v = df.TestFunction(V)\n",
    "f = df.Constant(0.0)\n",
    "a = df.dot(df.grad(u), df.grad(v))*df.dx\n",
    "L = f*v*df.dx\n",
    "\n",
    "# Compute solution\n",
    "u = df.Function(V)\n",
    "df.solve(a == L, u, bc)\n",
    "\n",
    "# Plot solution and mesh\n",
    "c = df.plot(u)\n",
    "plt.colorbar(c)\n",
    "#plot(mesh)\n",
    "\n",
    "print(\"L2_error is = \", df.errornorm(u_D, u, 'L2'))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Satisfying patch test is a sufficient condition that the code is correct. However, it is also shown that satisfaction of patch test also implies that the element is convergent. The other approach to check convergence is comparing with known solutions. However, for many geometries there is no knowledge of closed form solution. In that case, one uses the __method of manufactured solutions__ that was discussed in one of the previous classes and also uploaded on the course website. "
   ]
  },
  {
   "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.7.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}