{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## The Steady State Equation\n",
    "-----------\n",
    "\n",
    "\\\\[\n",
    "\\nabla \\cdot q = h\n",
    "\\\\]\n",
    "\n",
    "\\\\[\n",
    "q = - k \\nabla T\n",
    "\\\\]\n",
    "\n",
    "where:\n",
    " * $T$ is a scalar quantity \n",
    " * $k$ is diffusion (or conductivity) coefficient\n",
    " * $q$ is the heat flux vector \n",
    " * $h$ is the source/sink term\n",
    "\n",
    "-----\n",
    "\n",
    "Here we consider 2D models in the region, $0 \\leqslant x \\leqslant 1 $ and $ y_{0}\\leqslant y \\leqslant y_{1}$\n",
    "\n",
    "with no variation in the x direction, i.e. $ \\frac{\\partial T}{\\partial x} = 0 $ and a constant value $ k $ across the domain. \n",
    "\n",
    "Leading the 1D general solution:\n",
    "\n",
    "$ T = -\\frac{h}{2 k}y^{2} + c_{0}y + c_{1} $\n",
    "\n",
    "where $c_{0}, c_{1}$ are arbitrary constants found by applying each model's boundary conditions\n",
    "\n",
    "Three models are presented below, each with an analytic solution that the numerical results are tested against."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "# analytic solution definitions\n",
    "def analyticTemperature(y, h, k, c0, c1):\n",
    "     return -h/(2.*k)*y**2 + c0*y + c1\n",
    "\n",
    "def exactDeriv(y, h, k, c0):\n",
    "    return -h/k*y + c0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "import underworld as uw\n",
    "import underworld.visualisation as vis\n",
    "import numpy as np\n",
    "uw.utils.matplotlib_inline()\n",
    "import matplotlib.pyplot as pyplot\n",
    "import matplotlib.pylab as pylab\n",
    "pyplot.ion() # needed to ensure pure python jobs do now hang on show()\n",
    "\n",
    "rank = uw.mpi.rank\n",
    "\n",
    "# for machines without matplotlib #\n",
    "make_graphs = True\n",
    "try:\n",
    "    import matplotlib\n",
    "except ImportError:\n",
    "    make_graphs=False\n",
    "\n",
    "# depth range\n",
    "y0 = -.60\n",
    "y1 = 1.3\n",
    "\n",
    "# build mesh and fields\n",
    "mesh = uw.mesh.FeMesh_Cartesian( elementType = (\"Q1\"), \n",
    "                                 elementRes  = (10, 20), \n",
    "                                 minCoord    = (0., y0), \n",
    "                                 maxCoord    = (1., y1))\n",
    "\n",
    "tField    = mesh.add_variable( nodeDofCount=1, dataType='double')\n",
    "topWall    = mesh.specialSets['MaxJ_VertexSet']\n",
    "bottomWall = mesh.specialSets['MinJ_VertexSet']"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "## Model 1)\n",
    "\n",
    " * a fixed temperature condition or `DirichletCondition` - topWall \n",
    " \n",
    " $ T(x,y_{1}) = T_{1} $\n",
    "\n",
    " * a heat flow condition or `NeumannCondition` - bottomWall.\n",
    "\n",
    " $ q \\cdot n_{b} = (\\,0.0\\,,\\, f\\,) \\cdot (\\,0.0\\,,\\,-1.0\\,) = - f$\n",
    " \n",
    "     **Note** The heat flow is calculated using the heat flux vector $q$ multiplied by the outward surface normal, $n$. \n",
    "     The bottom surface outward normal $n_{b}$ point along the negative j-axis \n",
    "\n",
    "\n",
    "When the `NeumannCondition` object is associated with the `SteadyStateHeat` system it defines a flux along a boundary such that:\n",
    "   \n",
    "$ q \\cdot n = \\phi $ at $ \\Gamma_{\\phi} $\n",
    "\n",
    "where:\n",
    "* $ \\Gamma_{\\phi} $ is the set of vertices along the surface of the domain, \n",
    "* $\\phi $ is the scalar flow associated with the vector flux $q$ along $\\Gamma_{\\phi}$\n",
    "\n",
    "\n",
    "---------------\n",
    "\n",
    "An example: Defining a scalar field's flux at the bottom wall in a 2D rectangular domain.\n",
    "\n",
    "The outward facing normal vector at the bottom wall $\\mathbf{n}\\mid_{(x,y_{0})}=\\left[0,-1\\right] $) and the imposed flux vector $k \\nabla T = \\left[0, \\phi\\right]$\n",
    "\n",
    "The `NeumannCondition` object definition of this condition would be: \n",
    "\n",
    "```\n",
    "nbc = uw.conditions.NeumannCondition( fn_flux= -1.0 * phi, variable=tField,\n",
    "                                      indexSetsPerDof=mesh.specialSets[\"MinJ_VertexSet\"] )\n",
    "```\n",
    "\n",
    "Applies a 'fn_flux' to the scalar 'variable' `MeshVariable` over the boundary vertices in the set 'indexSetsPerDof'. The factor -1 is from the vector multiplication with the outward facing normal vector.\n",
    "\n",
    "Here `phi` can be a `underworld.Function` or `underworld.MeshVariable` type.\n",
    "\n",
    "------\n",
    "\n",
    "Arbitrary constants are:\n",
    "\n",
    "$c_{0} = \\frac{1}{k} (\\,f + hy_{0}\\,) $\n",
    "\n",
    "$c_{1} = T_{1} + \\frac{h}{2 k}y_{1}^2 - c_{0}y_{1} $\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [],
   "source": [
    "T1 = 8.0       # surface temperature\n",
    "k = 6.7        # diffusivity\n",
    "h = 8.0        # heat production, source term\n",
    "f = 2. \n",
    "\n",
    "# analytic solution definitions\n",
    "# 1 dirichlet conditions (top) + 1 neumann (bottom)\n",
    "c0 = 1./k*(f+h*y0)\n",
    "c1 = T1 + h/(2.0*k)*y1**2 - c0*y1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [],
   "source": [
    "for ii in topWall:\n",
    "    tField.data[ii] = T1\n",
    "\n",
    "# flag the dirichlet conditions on the topWall only\n",
    "bc = uw.conditions.DirichletCondition(tField, indexSetsPerDof=(topWall) )\n",
    "\n",
    "# define neumann condition\n",
    "nbc = uw.conditions.NeumannCondition( variable=tField,\n",
    "                                      fn_flux=-f,\n",
    "                                      indexSetsPerDof=(bottomWall))\n",
    "\n",
    "# flag the dirichlet conditions on the topWall only\n",
    "bc = uw.conditions.DirichletCondition(tField, indexSetsPerDof=(topWall) )\n",
    "\n",
    "# define heat eq. system\n",
    "ss = uw.systems.SteadyStateHeat( temperatureField  = tField,\n",
    "                                 fn_diffusivity = k,\n",
    "                                 fn_heating     = h,\n",
    "                                 conditions     = [bc, nbc] ) "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [],
   "source": [
    "solver = uw.systems.Solver(ss)\n",
    "solver.solve()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAtQAAAF3CAYAAACFTdwtAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvNQv5yAAAIABJREFUeJzs3Xd8leXdx/HPLyAjKIKARdEk2Gpx1BkXirbi3rVatbFqVaLVarXDUfSxtU2r1vr4tGptFGfTuheuuq0TjRvrKConljqwilIDsq7njzsIRkYgnNwZn/frdb9ycu47J98cwdc3F9d9XZFSQpIkSdLSKck7gCRJktSRWaglSZKkVrBQS5IkSa1goZYkSZJawUItSZIktYKFWpIkSWoFC7UkSZLUChZqSZIkqRUs1JIkSVIrWKglSZKkVuied4AlNXDgwFRRUZF3DEmSJHVyTz/99PsppUGLu67DFeqKigrq6+vzjiFJkqROLiIKLbnOKR+SJElSK1ioJUmSpFawUEuSJEmtYKGWJEmSWsFCLUmSJLWChVqSJElqBQu1JEmS1AoWakmSJKkVLNSSJElSK1ioJUmSpFawUGuB6urqqKiooKSkhIqKCurq6jr195UkSVpa3fMOoPanrq6O6upqGhsbASgUChw9ahTdp09n/29+E2bMgE8/XfDHGTNg1qx5LxbxxcfNn4uA5Zbjrgcf5IIzz6T39OlUAJ8WCpwyahTdp01j/4MPhuWW+/zXLuOfefTo0TQ0NFBWVkZNTQ1VVVVF+V6SJKlziZRS3hmWSGVlZaqvr887RptodclLCRob4f335x0ffAAffwxTp37xaHr+H+PG0XPmTFYAVgB60o7+KaNnz+zo0wdWWAH69v3isbDnBw6EQYNgwADoPu93yea/QACUlpZSW1trqZYkqQuLiKdTSpWLu66oI9QR8UNgFBDAxSml8xZwzdeB84DlgPdTStsWM1NHsaBR4urqakpmzuTAbbeFf/8bJk2C9977fGH+z38+//n06Yv+RssvP6+ArrACrLAC/5w5k6nw2TENmNF0/Pb//g969MhKbY8en38892O3btlI8vy/rM19vKDn5syBmTPZbeRIepAV+PmPXsBvf/WrbBT800+zn6mxcd4vBh9/DK+/nn2ce8yeveifu3//rFwPHMigZ5/lvGnTeB+YDNnHxkbqfvpTqnbcMSviRRoZlyRJHV/RRqgjYj3gamAzsi52F3BUSmnCfNf0Ax4Ddk4pNUTEyiml9xb1up12hDolmDw5K8mTJnHKIYfQ64MPGAKsCp99HLSwr19ppaz4zT0GDPj853OPlVaaV5779IGSL449V1RUUCgUvvB8eXk5EydOXGY/ctG+b0owbdq8sv3xxzBlyrxfMiZP/tzHF+6/n4Fk7+1yC3q9Hj1gyJB5x2qrffHxKqtk10mSpE6jPYxQrw2MSyk1NgV6CNgHOHu+a74D3JhSagBYXJnu8P77X3jzTXjjjS9+nDgxG3Vt8htgDvAu8G+gQPabx9vAGWPGZCVu1VVh8OBstLX7svtPWVNTs8ApEDU1NcvsexT1+0ZAaWl2fOlLi718z/mK/IrwWblef9Ag/nTaafCvf332iw5PPw233poV9uYGD4Y11ljwscoqC/zlBZy/LUlSR1fMEeq1gVuALclmDdwH1KeUjp3vmrlTPdYlm677fymlKxf1unmMUC9R4fnwQ3jllez45z/nleY338xGROe3/PJZ2Ro6dN7RNOK55be+Rf2kScxq9vLFHiWeK6+Sl8f3XeI51CllI97zF+1//QsaGub9gvTWW9lUlrl69sz++85fsr/8Zca++ioHnX46H89X0J2/LUlS+9DSEeqi3pQYEYcDRwOfAC8Bn6aUjp/v/PlAJTAS6A08DuyWUnqt2etUA9UAZWVlmyxoWkCxLKhs9endmz//+tfsPWwYvPzyvAL9yivZnOa5uneH8vJ5RWpuaZ77eMCAhc7N9Ua5trXMi/yMGVnBfuON7Hj99c8/njr1s0tnAv8EXm46XgGmDB7MbRMmZNNyJElSLtpFof7cN4r4NfCvlNKF8z13MtA7pXR60+djgLtSStct7HXaeoS6oqKC/xYK/AAY1nR8laz9f2allWDttWHYsHkfhw3LynQrpmI4FaCTSilbbeX11zlk8835Ktn8qLWBr9BsHlZZWfZnau6x7rqw/vrZHHhJklRU7aJQz73JMCLKgLuBLVJKU+Y7vzZwPrAT0AN4EjggpTR+Ya/Z1oW6pKSEgSnxDvAm2ejhy8CrwMUPP5wV54ED2yyPOpfmN2IuR1aqtxk4kIt++MPP/wvIfP9awZe/DBtuCBtsMO/j6qu7GokkSctQe7gpEeCGiBhA9q/ax6SUpkTEUQAppYtSSi9HxF3AC2T34F2yqDKdh7KyMgqFAqXAp/M9X15eDltvnVcsdRLNb8ScCRRKSxlx3nkw/79GzJmTzct+8UV4/nl47rns4w03zLumf/95BXtuyV5nHVcfkSSpyNzYZTGcy6xia9XUnqlTv1iyX3hh3iokyy2XFevNNsuOzTeHtdaCkhKnFEmStBjtYspHMbT7VT6kvM2ena0w8/zz8Mwz8NRT2fHf/2bn+/bl7dVX56pXX+XRWbN4EngHf1GUJKk5C7WkeWbPzuZhP/kkPPkk4y+9lGEzZnw256uB7AaGf/brxym33gqbbgq9euUYWJKk/FmoJS1USUkJPVNiI2Bzsu1MNwO+PPeCHj2yKSIjRmTH8OGw4op5xZUkKRctLdQL3rpNUqdWVlbGdLKF388j27L0K8DGq60GN98Mxx6braV99tmw667Z0pAbbQTHHQfXXw/vvptrfkmS2hNHqKUuqMU3237yCTzxBDz8cHY8/vi8Gx7XXDMbvd5mG9huu2zZPkmSOhGnfEhapKW62XbGjOxGx7kF+5FH4MMPs3NrrQXbb58dX/96toyfJEkdmIVaUvHNmQPjx8N998G998JDD2Wj2iUlUFmZleuRI7M52PPd5OjKOZKkjsBCLantzZgB48Zl5fq++7LpIrNnZ2V6xAjYfnvumDGDb//613wyd+oILtknSWqfLNSS8vfxx/D3v2cF+9574aWXAHgPuBO4Hbgb+Ihs99GJEyfmFlWSpObay9bjkrqyvn1h992zA+Dtt/nuqquyM7A7cAgwC3gUuKNQyAr3OutARG6RJUlaUi6bJ6ntrLIKD5eXcxCwMjAcOAtYsekj660HQ4fC0UfD7bfDfKuQSJLUXlmoJbWpmpoaSktLmUO2DvapwFalpdz4+99DbS1suCFceWU2qj1gQLYO9oUXwqRJOSeXJGnBLNSS2lRVVRW1tbWUl5cTEZSXl1NbW8s+xx4Lo0ZlG8v85z9w991w5JEwYQIccwysthpsvjmcdRa89lreP4YkSZ/xpkRJ7d8rr8BNN8GNN8Lcv//rrgvf/GZ2bLSR864lScucq3xI6pzeeisbxb7xxmwFkTlzoLx8Xrneaivo1i3vlJKkTsBCLanze/99GDs2K9f33AOffgqDBsFee8F++2Vbond3MSNJ0tJpaaF2DrWkjmvgQPje97JSPXkyXHtttjPjNdfATjvBqqtm868ffhjmzKGuro6KigpKSkqoqKigrq4u759AktQJOEItqfOZPh3uuguuvhpuvRWmTaOxf38umTqVq2bNYu7/QdyhUZK0KE75kCSA//4Xxo7lnsMPZ9tp0+gBvA5c3XRMLStjYqGQb0ZJUrvklA9JAlh+eTjwQHaaPp0vAd8DJgAnAS8Ctzc0wC9/mS3PJ0nSUrBQS+oSysrKmAJcDuwMrAp8H/hvz57wP/8Da66ZrRBSWwtTpuQZVZLUwVioJXUJc3donGsycGVpKRPGjMmW4jvrrKxIH3kkDB4M++8Pd9wBs2blF1qS1CFYqCV1CQvbobGqqirbhfHEE2H8eHjqqWzHxvvug912g9VXh5/8BF58Me8fQZLUTnlToiQtyIwZcPvtcMUV2cdZs2DDDeGQQ+A734GVV847oSSpyLwpUZJao0ePbOfFm2+Gt9+G3/8+24HxhBNgyBDYc89sST6nhEhSl2ehlqTFGTgQjj0W6uuzaSEnnJBNDdlrr2zb89NOg4kT804pScqJhVqSlsS668LZZ0NDA9x0UzYNpKYG1lgDdt452wZ95sy8U0qS2pCFWpKWxnLLwd57Z/Or33wzG6V+6SX41reyGxlPOQVef/2zy932XJI6r6LelBgRPwRGAQFcnFI6byHXbQo8DhyQUrp+Ua/pTYmS2q1Zs7Itzy++GG67DebMgZEjeXjttdljzBg+mjbts0vd9lyS2r/ctx6PiPXIdvbdDJgB3AUclVKa0Oy6bsA9wHTgUgu1pE5h0iS47DK45BIoFHgPuAz4IzB3o/Py8nImOvdaktqt9rDKx9rAuJRSY0ppFvAQsM8CrjsWuAF4r4hZJKltDRkCp54Kr7/OLsAjwI+B14GbgJFAQ6GwqFeQJHUQxSzU44ERETEgIkqBXYHV578gIoYA3yQbtJGkzqdbN14uL+dbwFDgTGAr4F7gte7d4YILYOrUXCNKklqnaIU6pfQycBZwN9l0j+eA2c0uOw84KaU0Z1GvFRHVEVEfEfWTJ08uSl5JKpa5257/CziVbGRhVI8e9F99dfjBD7LR7OOOg1dfzTmpJGlpFHWVj5TSmJTSJimlbYAPgdeaXVIJXB0RE4F9gQsjYu8FvE5tSqkypVQ5aNCgYkaWpGWu+bbng8vL+fqllzLgjTdg3LhsPes//QmGDYMdd4SxY2F28/EHSVJ7VexVPlZOKb0XEWVkI9VbpJSmLOTay4HbvClRUpf07rvZDYx//GN2Q+PQoXD00XDEEdCvX97pJKlLag83JQLcEBH/AMYCx6SUpkTEURFxVJG/ryR1LF/6Eowena1pfd11UFYGP/1ptqb18ce7E6MktWNFHaEuBkeoJXUZzz4L554LV1+drWm9777w4x/DZpvlnUySuoT2MkItSVpaG20EV12VjVr/5Cfwt7/B5pvDiBFw883Os5akdsJCLUnt3WqrwVlnwVtvwXnnZR+/+c3sJsYLL4TGxrwTSlKXZqGWpI5ihRXghz+ECRPgmmtgpZXgmGOyedanngrvvJN3QknqkizUktTRdO8O3/42PPEEPPwwbLMN/PrXUF4Oo0ZlhVuS1GYs1JLUUUXA1lvDTTdlm8Icfng25/qrX4UDD4QXXqCuro6KigpKSkqoqKigrq4u79SS1OlYqCWpM1hzzWw+9cSJ2Q2Mt98OG2zASgcfzKqFAiklCoUC1dXVlmpJWsYs1JLUmQwenN3AWChwzoorsumcOTwGPADsADQ2NjJ69OicQ0pS52KhlqTOqH9/Tvz4Y8qB44GvkG1X+xSwSaGQrWstSVomLNSS1EmVlZXRCPwf8GXgCKAfcAPAeuvBlVfCzJk5JpSkzsFCLUmdVE1NDaWlpQDMAMYAG/fuzSPHHAPLLQeHHJLNvb74You1JLWChVqSOqmqqipqa2spLy8nIigvL+ePF1/M1uefD889B2PHwsorQ3U1rLUWjBljsZakpRAppbwzLJHKyspUX1+fdwxJ6hxSgjvvhNNPh/p6GDo02yTmu9/NRrElqQuLiKdTSpWLu84RaknqyiJg113hySfhttuy3RcPPzzb1vyyy2DWrLwTSlK7Z6GWJGXFerfd4Kmn4NZboV8/OOywrFhffrnFWpIWwUItSZonAvbYI5v+ccst0LcvfO97sPba2aogFmtJ+gILtSTpiyJgzz3h6afh5pth+eWzVUHWWQf+8hfXsZak+VioJUkLFwF77ZUV65tugt69oaoKNt442968g93YLknFYKGWJC1eSQnsvTc8+yzU1cHUqbD77rDNNvDII3mnk6RcWaglSS1XUgLf+Q68/DJceCFMmAAjRmQ3ND7/fN7pJCkXFmpJ0pLr0QO+/314/XU480x47DHYcMOsbE+YQF1dHRUVFZSUlFBRUUFdXV3eiSWpaCzUkqSlV1oKJ50Eb7wBp5wCt9zCnGHDaDz0UGYUCqSUKBQKVFdXW6oldVoWaklS6/XvD7/+Nbz+OleVlnLIrFlMAH4D9AcaGxsZPXp0ziElqTgs1JKkZWfwYL733/8yDLgBOBF4HTgBeKdQyDWaJBWLhVqStEyVlZXxJnAwsAHwBHAu8Fr37nDttS61J6nTsVBLkpapmpoaSktLARgP7Ars2bMnK6y6Kuy/PwwfDo8+mmtGSVqWLNSSpGWqqqqK2tpaysvLiQjKy8vZf8wY+r/xBlx6KRQKsPXWsO++2bJ7ktTBRepg//RWWVmZ6uvr844hSVpan3wCv/sdnH02zJgBRx8Np50GAwbknUySPicink4pVS7uuqKOUEfEDyNifES8FBHHL+B8VUS8EBEvRsRjEbFBMfNIktqBPn3gf/4H/vlPOPRQ+MMf4MtfhnPOgU8/zTudJC2xohXqiFgPGAVsRnZfyu4R8ZVml70JbJtS+hrwS6C2WHkkSe3MKqtAbW22w+Lw4fDTn8KwYXDNNd64KKlDKeYI9drAuJRSY0ppFvAQsM/8F6SUHkspfdj06RPAakXMI0lqj9ZbD+64A+6+G/r2hQMOyLYzf+aZvJNJUosUs1CPB0ZExICIKCW70Xv1RVx/OHBnEfNIktqzHXbISvQll8Brr0FlJYwaBe+9l3cySVqkohXqlNLLwFnA3cBdwHPA7AVdGxHfICvUJy3kfHVE1EdE/eTJk4uUWJKUu27d4PDDs/nVJ5wAl18Oa64J556b3cAoSe1QUW9KTCmNSSltklLaBvgQeK35NRGxPnAJsFdK6T8LeZ3alFJlSqly0KBBxYwsSWoPVlwxWwnkxRez+dU//jGsvz7cdVfeySTpC4q9ysfKTR/LyOZP/6XZ+TLgRuC7KaUvlG1JUhc3bFg2v3rsWJg9G3bZBfbYIxvBlqR2otgbu9wQEf8AxgLHpJSmRMRREXFU0/n/AQYAF0bEcxHhAtOSpM+LgN13h/Hjs7WrH3oI1l0XTjoJpk7NO50kubGLJKmDeecdOOWUbH714MFw5pnw3e9CiZv/Slq22sXGLpIkLXODB8Nll8G4cVBenm0Os9VW8OyzeSeT1EVZqCVJHdNmm8Fjj8EVV8Abb2TL7B13HNdefDEVFRWUlJRQUVFBXV1d3kkldXLd8w4gSdJSKymBgw+GPfeEU08lnX8+26bEVkABKBQKVFdXA1BVVZVrVEmdlyPUkqSOr18/OP989vjSlygAdcD9wDCgsbGR0aNH55tPUqdmoZYkdRp3vPsuWwJHAhsALwC/Ad4vFHLNJalzs1BLkjqNsrIy5gC1wFeBq4CTgVe7dYObboIOtrKVpI7BQi1J6jRqamooLS0F4H3gcGD7nj0pXXVV2GefbD3rN97INaOkzsdCLUnqNKqqqqitraW8vJyIoLy8nO+NGUP/N96Ac8+Fv/8d1lkHzjgDPv0077iSOgk3dpEkdR2TJsGPfwzXXJNta15bCyNG5J1KUjvlxi6SJDU3ZAhcfTXccQdMmwbbbANHHglTpuSdTFIHZqGWJHU9u+wCL72UjVZfcgmsvTZcf703LUpaKhZqSVLX1KcPnHMOPPkkrLIK7Lcf7L03vPVW3skkdTAWaklS17bJJlmpPuccuPfe7KbFP/wBZs/OO5mkDsJCLUlS9+7Z9I/x42GrreC442D4cHjhhbyTSeoALNSSJM01dCjceSfU1cGbb2aj1z/7WXYDoyQthIVakqT5RcB3vgMvvwwHHQS/+Q187WvwwAN5J5PUTlmoJUlakAED4LLL4L77ss+32w6OPhqmTs03l6R2x0ItSdKibLddNpf6hBPgoouy0ep77sk7laR2xEItSdLilJZmW5c/8gj06gU77gijRsFHH+WdTFI7YKGWJKmlhg+HZ5+FE0+ESy+F9dbLbmKU1KVZqCVJWhK9e8NZZ8Hjj0PfvrDrrnDooVxXW0tFRQUlJSVUVFRQV1eXd1JJbcRCLUnS0thsM3jmGRg9mjlXXcXWRx3F1woFUkoUCgWqq6st1VIXYaGWJGlp9ewJv/oVe6y8Mu+lxFjgKmAloLGxkdGjR+ccUFJbsFBLktRKd777LpsCpwP7Ay8BewMNDQ255pLUNizUkiS1UllZGTOBM4BNgEnATcD1paUwZUqu2SQVn4VakqRWqqmpobS0FIAXgS2AXy23HHtPm5atW33vvbnmk1RcFmpJklqpqqqK2tpaysvLiQiGlJcz9LLLKBk3DpZfHnbYAY49Fhob844qqQgipZR3hiVSWVmZ6uvr844hSVLLTJsGP/sZnHcerLUWXHklbL553qkktUBEPJ1SqlzcdUUdoY6IH0bE+Ih4KSKOX8D5iIjfR8SEiHghIjYuZh5Jktpc797wv/8L998P06dnm8OceirMmJF3MknLSNEKdUSsB4wCNgM2AHaPiK80u2wXYM2moxr4Y7HySJKUq298A154AQ4+GGpqslHq8ePzTiVpGSjmCPXawLiUUmNKaRbwELBPs2v2Aq5MmSeAfhGxShEzSZKUnxVXhMsug5tvhkmTYJNN4Le/hdmz804mqRWKWajHAyMiYkBElAK7Aqs3u2YI8NZ8n/+r6TlJkjqvvfbKRqd33RVOPBG+/nV44428U0laSkUr1Cmll4GzgLuBu4DngKX6FTwiqiOiPiLqJ0+evAxTSpKUk5VXhhtvhCuuyKaCrL8+XHopdLDFAiQV+abElNKYlNImKaVtgA+B15pdMonPj1qv1vRc89epTSlVppQqBw0aVLzAkiS1pYhsTvWLL8Kmm8Lhh8O3vw0ffJB3MklLoEWFOiJWi4hvND3uGRF9Wvh1Kzd9LCObP/2XZpfcChzctNrHFsBHKaW3W5xekqTOoKws2/zlzDOz+dUbbAAPPph3KkkttNhCHRGHkRXfS5qeKgduaeHr3xAR/wDGAseklKZExFERcVTT+TuAN4AJwMXA0UsSXpKkTqNbNzjpJHjiCSgthe22g5NPdnk9qQNY7MYuEfEc2dJ341JKGzU990JKaf02yPcFbuwiSer0PvkEjj8eLrkkWwnkL3/JNoWR1KaW5cYu01NKn/16HBHdgGhNOEmStAh9+sDFF8MNN8Cbb8JGG2Xl2hsWpXapJYX60Yg4EejVNI/6GuC24saSJEnss0+2AsgWW8CoUbDfft6wKLVDLSnUJwJTgVeAHwL3AaOLGUqSJDUZMgTuuQfOPhtuvTVbXu/++/NOJWk+iyzUTdM7Lksp/TGl9M2U0t5Nj+e0UT5JklRSAj/9aXbD4vLLw/bbZzcwesOi1C4sslCnlGYDa0TEcm2UR5IkLczGG8PTT0N1dTZiPXw4vP46dXV1VFRUUFJSQkVFBXV1dXknlbqU7i245nXg4Yi4Bfhk7pMppd8XLZUkSVqwPn3gootg553hsMOYud563D1nDoWm0epCoUB1dTUAVVVVeSaVuoyWzKFuAO4BSoFB8x2SJCkve+8Nzz3H8ylxxYwZ/Ano3XSqsbGR0aO93UlqK4sdoU4pndYWQSRJ0hIqK2P4p5/yc+BnwHDg28DLQENDQ57JpC5lsYU6Iu4BvrDwZUppx6IkkiRJLbZqeTmjCwUeBK4C6oEfAPevvnquuaSupCVzqE+d73Ev4FvAp8WJI0mSlkRNTQ3V1dXc09jIhsCfgUuBN1dbDaZOhRVWyDmh1Pktdg51SmncfMdDKaXjgG3aIJskSVqMqqoqamtrKS8v590IRpWV8fy++zL0iSeyVUGeeSbviFKnt9hCHRF95zv6RcRIoH8bZJMkSS1QVVXFxIkTmTNnDm8UCmxw3XXwwAMwbRpsuSX84Q9uWy4VUUumfLxENoc6gFnAm8CoYoaSJEmttM028NxzcOihcNxx2e6Kl14K/R0Tk5a1liybt0ZKqSyltHpKaWhKaTvg0WIHkyRJrTRwIIwdC7/7Hdx+O2y4IYwbl3cqqdNpSaFe0N+8J5d1EEmSVAQR8KMfwSOPZFuYjxgB55/vFBBpGVpooY6IlSNiA6B3RHwtItZvOrYm2+RFkiR1FJttlm1bvtNOcOyxcOCB2SogklptUXOodwMOA1YDLpzv+amAm71IktTRrLQS3HILnH02jB6dzbG+4QZYd928k0kd2kJHqFNKl6WURgCHp5RGzHfsmlK6rg0zSpKkZaWkBE4+Ge67D6ZMyUau//znvFNJHVpLth6/NiJ2AtYl29hl7vO/LmYwSZJURF//Ojz7LBxwAHz3u9kc6/POg169Fvulkj6vJetQXwgcAvwI6A0cBHylyLkkSVKxrbJKNlJ90knwpz/B1lvDm2/mnUrqcFqyysfWKaXvAP9JKZ0GbI6FWpKkzqF7dzjzzGxu9YQJ2e6KY8fmnUrqUFpSqKfP/RgRg5s+X7V4kSRJUpvbc89sm/I11sgen3IKzJqVdyqpQ2hJob4jIvoB5wDPAROBa4sZSpIk5WCNNeDRR6G6Ohu13n57eOedvFNJ7d4iC3VElAB3ppSmNK3sMRT4WkrpZ22STpIkta1evbL51FdcAU8+mU0BefzxvFNJ7doiC3VKaQ7wp/k+n5ZS+qDoqSRJUr4OPhieeAJ694Ztt4WLLnJ3RWkhWjLl44GI2KvoSSRJUvuy/vpQX59N/fj+9+GII2D69MV/ndTFtKRQHwrcFBHTIuKDiPgwIhylliSpK+jfP1v149RT4dJLYcQIaGjIO5XUrrSkUA8ElgOWBwY1fT6oJS8eESdExEsRMT4i/hoRvZqdL4uIByLi2Yh4ISJ2XdIfQJIkFVm3bvDLX8JNN8Grr8Imm8ADD+SdSmo3FluoU0qzgf2Ak5oerwJsuLivi4ghwHFAZUppPaAbcECzy04Frk0pbdR07sIliy9JktrM3ntnNyoOHAg77MDTVVVUlJdTUlJCRUUFdXV1eSeUctGSnRLPB74BfLfpqUbgoha+fnegd0R0B0qBfzc7n4C+TY9XXMB5SZLUngwbBk8+ScPGG7PJX/7Cbxoa6J0ShUKB6upqS7W6pJZM+RieUjqSpg1emlb56LG4L0opTSJbu7oBeBv4KKV0d7PLfg4cFBH/Au4Ajm15dEmSlIsVVmCbd9/lFGB/4HFgDaCxsZHRo0fgDRRuAAAfzElEQVTnm03KQUsK9cym9agTQEQMAOYs7osioj+wF9na1asCfSLioGaXHQhcnlJaDdgVuKrpezV/reqIqI+I+smTJ7cgsiRJKqaGt97iTGBnYDWgvulxgzcsqgtqSaG+ALgBGBQRvwAeAc5qwddtD7yZUpqcUpoJ3AgMb3bN4TTtuphSehzoRXbT4+eklGpTSpUppcpBg1p0P6QkSSqisrIyAO4BNgEKwO3A2X37ul61upyW3JR4JdnNg+cAHwD7pZSubsFrNwBbRERpRAQwEnh5AdeMBIiItckKtUPQkiS1czU1NZSWlgIwkWzE7Npu3fjJRx/Bt78Nn3ySZzypTbVkhBqyFTpmAjNa+jUppXHA9cAzwItNX1cbEWdExJ5Nl/0YGBURzwN/BQ5NyV9rJUlq76qqqqitraW8vJyIYOXycmZffjmccw7ceCNstRUUCnnHlNpELK6/RsRo4DvATUCQzYuuSyn9pvjxvqiysjLV19fn8a0lSVJL3HknHHgg9OgBN9yQbQYjdUAR8XRKqXJx17VktPlgYNOU0qkppdHAZmS7J0qSJH3RLrvAuHHZLovbbQcXX5x3IqmoWlKo3yZbT3qu7k3PSZIkLdhXv5qV6pEjoboafvADmDkz71RSUbSkUH8AvBQRl0TExWTzod+PiHMj4tzixpMkSR1Wv35w++3w4x/DBRfATjvB++/nnUpa5rov/hJubzrmeqJIWSRJUmfTrVt2o+L662cj1ZttBrfcAl/7Wt7JpGVmsYU6pTSmLYJIkqRO7OCDs2kg3/wmbLkl/PnPsPfeeaeSlonFTvmIiJ0j4qmIeC8iPoiIDyPig7YIJ0mSOpHNN4ennoJ11smK9a9+5SYw6hRaMof6fOBIYAgwiGwnQ7crlCRJS27IEHjoITjoIDjtNNh/fzeBUYfXkkL9L+C5lNLMlNLsuUexg0mSpE6qd2+48ko4+2y4/nrYZhuYNCnvVNJSa8lNiScCYyPiQeDTuU+mlH5frFCSJKmTi4Cf/hTWXjvbBGazzWDsWNh447yTSUusJSPUvwBmA/3IpnrMPSRJklpn993h0Uez1UBGjICbb847kbTEWjJCvXpKab2iJ5EkSV3T+uvDk0/CXnvBPvvAmWdmo9cReSeTWqQlI9R/i4jtip5EkiR1XYMHw4MPwn77wUknwRFHwIwZeaeSWqQlI9SHASdERCMwAwggpZRWKmoySZLUtfTuDX/9a7Ze9S9/CW+8ATfcACtZOdS+tWSEeiCwHLAiLpsnSZKKqaQEzjgj2/jlscdgiy3gtdfyTiUt0mILddMSefsBJzU9XgXYsNjBJElSF1ZVBfffD1OmZKX6gQfyTiQtVEt2Sjwf+Abw3aanGoGLihlKkiSJrbaCceNglVVgxx3hkkvyTiQtUEumfAxPKR0JTAdIKX0A9ChqKkmSJIChQ7OpHyNHwqhR2eofs91fTu1LSwr1zIgoARJARAwA5hQ1lSRJ0lwrrgi33QbHHAPnnAPf+hZXjxlDRUUFJSUlVFRUUFdXl3dKdWELXeUjIrqnlGYBFwA3AIMi4hfAt8k2e5EkSWob3bvD+efDV7/KnB/+kK+MHcu0OXNIQKFQoLq6GoCqqqp8c6pLipTSgk9EPJNS2rjp8brA9mRL5t2bUhrfdhE/r7KyMtXX1+f17SVJUs6OWHllfj95Mu8CuwKvND1fXl7OxIkT8wumTicink4pVS7uukWtQ/3Z9kQppZeAl5ZFMEmSpNa49P33eR64DXgM2Bv4O9DQ0JBrLnVdiyrUgyLiRws7mVI6twh5JEmSFqmsrIz6QoEtgDuAu4HvAY+VleUbTF3Wom5K7AYsD6ywkEOSJKnN1dTUUFpaykRgOPA48Bfg5s03h4VMZZWKaVEj1G+nlM5osySSJEktMPfGw9GjR9PQ0ED16qvzt9VWY8Nrr81WBLnwwuwmRqmNLGqEOhZxTpIkKTdVVVVMnDiROXPm8FqhwNCHH4af/Qwuvhj22AOmTs07orqQRRXqkW2WQpIkqTVKSqCmBmpr4Z57YJttYNKkvFOpi1hooW7aEVGSJKnjGDUq2wRmwgTYYgt48cW8E6kLaMlOiZIkSR3HzjvDww/DnDmw9dZw7715J1InV9RCHREnRMRLETE+Iv4aEb0WcM23I+IfTdf9pZh5JElSF7HhhvDEE1BWBrvsAldemXcidWJFK9QRMQQ4DqhMKa1HtgzfAc2uWRM4BdgqpbQucHyx8kiSpC5m9dXhkUey+dSHHAJnneWyeiqKYk/56A70jojuQCnw72bnRwEXpJQ+BEgpvVfkPJIkqStZcUW44w444AA4+WQ4/vhsKoi0DBVtkcaU0qSIOAdoAKYBd6eU7m522VoAEfEo2Qj2z1NKdxUrkyRJ6oJ69oS6Ohg8GM47D955J5sC0rNn3snUSRRzykd/YC9gKLAq0CciDmp2WXdgTeDrwIHAxRHRbwGvVR0R9RFRP3ny5GJFliRJnVVJCZx7Lvz2t3DttdmNix99lHcqdRLFnPKxPfBmSmlySmkmcCPZDqHz+xdwa0ppZkrpTeA1soL9OSml2pRSZUqpctCgQUWMLEmSOq0I+MlP4M9/nje3+t/NZ6NKS66YhboB2CIiSiMiyDaKebnZNTeTjU4TEQPJpoC8UcRMkiSpq6uqgttvhzfegOHD4dVX806kDq5ohTqlNA64HngGeLHpe9VGxBkRsWfTZX8D/hMR/wAeAH6aUvpPsTJJkiQBsOOO8OCDMG1aVqoffzzvROrAInWw5WMqKytTfX193jEkSVJn8PrrsNNO2dSPa6+F3XfPO5HakYh4OqVUubjr3ClRkiR1XV/+Mjz2GKyzDuy9N4wZk3cidUAWakmS1LWtvHI2/WP77eGII+BXv3IDGC0RC7UkSdLyy8PYsfDd78Jpp8Exx8Ds2XmnUgdRtI1dJEmSOpTlloMrrsg2gPntb+GDD7INYHr0yDuZ2jkLtSRJ0lwRcPbZMHAgnHQSTJkCN9wAffrknUztmFM+JEmSmjvxRLjkErjnHthhh2y0WloIC7UkSdKCHH44XHcdPP00bLutuypqoSzUkiRJC7PPPnDHHTBxImy9Nbf87ndUVFRQUlJCRUUFdXV1eSdUO2ChliRJWpSRI+H++5k+eTJb/OQnrFgokFKiUChQXV1tqZaFWpIkabE23ZTd+vZlBvAQsFXT042NjYwePTrHYGoPLNSSJEkt8MDbb7MV8C5wN7BL0/MNDQ35hVK7YKGWJElqgbKyMt4CtgZeBm4BvtP0vLo2C7UkSVIL1NTUUFpayvvAN4BHgDrgum23zTeYcmehliRJaoGqqipqa2spLy/nvxEcVVbGW5tswqZXXgm/+AWklHdE5cSdEiVJklqoqqqKqqqqeU/MmgWjRsHPfw4ffQS/+12226K6FAu1JEnS0ureHcaMgb594X//F6ZOhYsugm7d8k6mNmShliRJao2SEjjvPFhhBaipgU8+gSuugOWWyzuZ2oiFWpIkqbUi4Fe/guWXh1NOgcZGuPpq6NUr72RqA96UKEmStKycfDKcfz7ccgvsuWc2Wq1Oz0ItSZK0LB1zDFx2Gdx3H+y0U3azojo1C7UkSdKyduih2ZSPceNg5Eh4//28E6mILNSSJEnFsN9+cPPNMH48bLstvP123olUJBZqSZKkYtltN7jzTigUYMSI7KM6HQu1JElSMX3jG3DPPdm0jxEj4LXX8k6kZcxCLUmSVGxbbgkPPgjTp8M228CLL+adSMuQhVqSJKktbLgh/P3v2S6K224LTz2VdyItIxZqSZKktjJsGDz8MPTrB9tvD48/nnciLQNFLdQRcUJEvBQR4yPirxGxwO2CIuJbEZEiorKYeSRJknK3xhrw0EOw8sqw445ZwVaHVrRCHRFDgOOAypTSekA34IAFXLcC8ENgXLGySJIktSurr56V6iFDYOed4YEH8k6kVij2lI/uQO+I6A6UAv9ewDW/BM4Cphc5iyRJUvux6qpZqR46FHbdFe6+O+9EWkpFK9QppUnAOUAD8DbwUUrpc39SImJjYPWU0u3FyiFJktRufelL2ej0WmvBnnvCHXfknUhLoZhTPvoDewFDgVWBPhFx0HznS4BzgR+34LWqI6I+IuonT55crMiSJEltb9AguP9+WHdd2HtvuOWWvBNpCRVzysf2wJsppckppZnAjcDw+c6vAKwHPBgRE4EtgFsXdGNiSqk2pVSZUqocNGhQESNLkiTlYMAAuO8+2Ggj2HdfuOGGvBNpCRSzUDcAW0REaUQEMBJ4ee7JlNJHKaWBKaWKlFIF8ASwZ0qpvoiZJEmS2qd+/bIdFTfbDPbfH66+Ou9EaqFizqEeB1wPPAO82PS9aiPijIjYs1jfV5IkqcPq2xf+9jfYaiuoqoKrrso7kVogUkp5Z1gilZWVqb7eQWxJktSJffIJ7LVXNrf6kkvgsMPyTtQlRcTTKaXF7pPiTomSJEntTZ8+MHZstvHL4YfDn/6UdyItgoVakiSpPerdG26+GXbfHY46itNXWomSkhIqKiqoq6vLO53mY6GWJElqr3r14q/77sst3brxiw8/5AcpUSgUqK6utlS3IxZqSZKkduyU009n39mzuRH4PXAs0NjYyOjRo3NOprm65x1AkiRJC9fQ0EAC9geuISvVCbigoSHXXJrHQi1JktSOlZWVUSgUmAUcQFaq/wD079cv32D6jFM+JEmS2rGamhpKS0sBmEk2Uj22WzfO+PBDuOCCXLMpY6GWJElqx6qqqqitraW8vJyIYNXycv47Zky2TvUPfmCpbgfc2EWSJKkjmjED9tsPbr01K9VHH513ok7HjV0kSZI6sx494LrrYI894Jhj4I9/zDtRl2WhliRJ6qjmlurdd89GqC+6KO9EXZKFWpIkqSPr2ROuvx522w2+/323Kc+BhVqSJKmj69kTbrghK9VHHWWpbmMWakmSpM6geamurc07UZdhoZYkSeos5pbqXXeFI4+ESy7JO1GXYKGWJEnqTOaW6p13hupquPLKvBN1ehZqSZKkzqZXL7jxRthuO/je9+Caa/JO1KlZqCVJkjqj3r3hlltgq62gqgpuuinvRJ2WhVqSJKmz6tMHbr8dNt0U9t8fbrst70SdkoVakiSpM1thBbjzTlh/ffjWt+Duu/NO1OlYqCVJkjq7fv2yIj1sGOy1Fzz4YN6JOhULtSRJUlew0kpw772wxhrZVuWPPpp3ok7DQi1JktRVDBoE990HQ4bALrvAk0/mnahTsFBLkiR1JYMHZ6V60CDYaSd45pm8E3V4FmpJkqSuZrXV4P77oW9f2HFHePHFvBN1aBZqSZKkrqi8PCvVPXvC9tvDK6/knajDslBLkiR1VV/+claqI7JdFSdMyDtRh2ShliRJ6sq++tVsTvXMmVmpLhTyTtThFLVQR8QJEfFSRIyPiL9GRK9m538UEf+IiBci4r6IKC9mHkmSJC3AuuvCPffAxx9n0z/eeSfvRB1K0Qp1RAwBjgMqU0rrAd2AA5pd9mzT+fWB64Gzi5VHkiRJi7DhhtmOim+/DTvsAB98kHeiDqPYUz66A70jojtQCvx7/pMppQdSSo1Nnz4BrFbkPJIkSVqYLbeEW26Bf/4Tdt45G7HWYhWtUKeUJgHnAA3A28BHKaVFbR5/OHDngk5ERHVE1EdE/eTJk5d9WEmSJGVGjoTrrsvWp95jD2hsXPzXdHHFnPLRH9gLGAqsCvSJiIMWcu1BQCXw2wWdTynVppQqU0qVgwYNKlZkSZIkQVakr7oKHn6YScOHs2Z5OSUlJVRUVFBXV5d3unanmFM+tgfeTClNTinNBG4Ehje/KCK2B0YDe6aUPi1iHkmSJLXUgQfyxGGHMeT55/lNQwMlKVEoFKiurrZUN1PMQt0AbBERpRERwEjg5fkviIiNgD+Rlen3iphFkiRJS+iAe+/lBGBf4BIggMbGRkaPHp1vsHame7FeOKU0LiKuB54BZpGt6FEbEWcA9SmlW8mmeCwPXJd1bhpSSnsWK5MkSZJarqGhgfOAFYAzgKlkS7g1NDTkmqu9KVqhBkgpnQ6c3uzp/5nv/PbF/P6SJElaemVlZRQKBX4J9AV+AnwMXFxWlm+wdsadEiVJkrRANTU1lJaWAvBTsnm6o4GbNt88z1jtjoVakiRJC1RVVUVtbS3l5eVEBGeVlfHm8OFsdO21cMEFecdrNyKllHeGJVJZWZnq6+vzjiFJktQ1zZwJ++2XbQBz+eVwyCF5JyqaiHg6pVS5uOscoZYkSVLLLbccXH01bL89HHYY3Hhj3olyZ6GWJEnSkunVC26+GTbfHA48EO6/P+9EubJQS5Ikacn16QO33QZrrQV77QVPPZV3otxYqCVJkrR0VloJ/vY3GDgQdtkFXnkl70S5sFBLkiRp6a26KtxzD3TvDjvsAF1w0xcLtSRJklrnK1/JRqqnToUdd4TJk/NO1KYs1JIkSWq9DTaAsWOhUMimf3z8cd6J2oyFWpIkScvGiBFw/fXw3HOw994wfXreidqEhVqSJEnLzm67wRVXwAMPZEvqzZqVd6Kis1BLkiRp2aqqgv/7v2yt6upq6GA7cy+p7nkHkCRJUid03HHwn//AGWfAgAFw9tkQkXeqorBQS5IkqTh+/vOsVJ9zTlaqTz4570RFYaGWJElScUTA738PH3wAp5ySlepRo/JOtcxZqCVJklQ8JSVw+eUwZQoceST07w/77pt3qmXKmxIlSZJUXD16ZMvpbblldsPiAw/knWiZslBLkiSp+EpLs41fvvKVbI3q55/PO9EyY6GWJElS21hpJbjrLujbF3beGSZOzDvRMmGhliRJUttZffWsVE+fDjvtBJMn552o1SzUkiRJalvrrgu33QYNDbD77vDJJ3knahULtSRJktreVlvBNddAfT3stx/MnJl3oqVmoZYkSVI+9twTLroI7rwTjjiiw25R7jrUkiRJys+oUfD223D66TB4MJx1Vt6JlpiFWpIkSfk67bSsVJ99NqyyChx/fN6JloiFWpIkSfmKgPPPh/fegxNOgC99CQ48MO9ULVbUOdQRcUJEvBQR4yPirxHRq9n5nhFxTURMiIhxEVFRzDySJElqp7p1g7o62GYbOOQQuOeevBO1WNEKdUQMAY4DKlNK6wHdgAOaXXY48GFK6SvA/wIdb9KMJEmSlo1eveCWW2DYMGbuuSe7rbIKJSUlVFRUUFdXl3e6hSr2Kh/dgd4R0R0oBf7d7PxewBVNj68HRkZEFDmTJEmS2qt+/bixupq3P/2US995h6EpUSgUqK6ubrelumiFOqU0CTgHaADeBj5KKd3d7LIhwFtN188CPgIGFCuTJEmS2r8fnXMOO6ZEN+BvwMpAY2Mjo0ePzjnZghVzykd/shHoocCqQJ+IOGgpX6s6Iuojon5yJ9ieUpIkSQvX0NDAq8DuwDtAzPd8e1TMKR/bA2+mlCanlGYCNwLDm10zCVgdoGlayIrAf5q/UEqpNqVUmVKqHDRoUBEjS5IkKW9lZWUAjANGAO82e769KWahbgC2iIjSpnnRI4GXm11zK3BI0+N9gftT6qBb5EiSJGmZqKmpobS09HPPlZaWUlNTk1OiRSvmHOpxZDcaPgO82PS9aiPijIjYs+myMcCAiJgA/Ag4uVh5JEmS1DFUVVVRW1tLeXk5EUF5eTm1tbVUVVXlHW2BoqMNCFdWVqb6+vq8Y0iSJKmTi4inU0qVi7uu2MvmSZIkSZ2ahVqSJElqBQu1JEmS1AoWakmSJKkVLNSSJElSK1ioJUmSpFawUEuSJEmtYKGWJEmSWsFCLUmSJLWChVqSJElqhQ639XhETAYKeefIyUDg/bxDdGC+f63j+9c6vn+t4/vXOr5/reP71zod+f0rTykNWtxFHa5Qd2URUd+S/eS1YL5/reP71zq+f63j+9c6vn+t4/vXOl3h/XPKhyRJktQKFmpJkiSpFSzUHUtt3gE6ON+/1vH9ax3fv9bx/Wsd37/W8f1rnU7//jmHWpIkSWoFR6glSZKkVrBQtzMRsVJE3BMR/2z62H8h182OiOeajlvne35oRIyLiAkRcU1E9Gi79PlryfsXERtGxOMR8VJEvBAR+8937vKIeHO+93bDtv0J2l5E7BwRrzb9mTl5Aed7Nv1ZmtD0Z6tivnOnND3/akTs1Ja524sWvH8/ioh/NP1Zuy8iyuc7t8C/x11NC97DQyNi8nzv1RHznTuk6e/7PyPikLZN3j604P373/neu9ciYsp857r0n8GIuDQi3ouI8Qs5HxHx+6b39oWI2Hi+c/7ZW/z7V9X0vr0YEY9FxAbznZvY9PxzEVHfdqmLJKXk0Y4O4Gzg5KbHJwNnLeS6/y7k+WuBA5oeXwR8P++fqb29f8BawJpNj1cF3gb6NX1+ObBv3j9HG75f3YDXgTWAHsDzwDrNrjkauKjp8QHANU2P12m6vicwtOl1uuX9M7XD9+8bQGnT4+/Pff+aPl/g3+P/b+9+Y+QqqziOf3/QFAyKtFSxKWi6hgQNqUAa0mBTsCC2vKAkYFIVrUJSmxjfGBNC1heVxD/vSPwTY6xtQBMwLBBLiErrtqlGq4J/2mqxLK3R1soqrRgpWVSOL54z5DLu7M46Ozt3dn6f5GbvPXOfmeeePHf27N3nzgzS0mYOPwp8ZZK2i4Gj+XNRri/q9THVLX9N+38S2F7ZHugxCKwBrgIOtXj8JuB7gIBVwM8yPvBjr838XdPIC7C+kb/c/gOwpNfHMFuLr1DXzwbgvly/D7il3YaSBKwFRv6f9vPEtPmLiCMR8Uyu/xkYB6b90PZ56mpgLCKORsTLwIOUHFZVczoCXJ9jbQPwYERMRMQxYCyfb5BMm7+I2BMRZ3JzP3DxHPex7toZg628D9gVEaci4jSwC1jXpX7W1Uzz9wHggTnpWR+IiH3AqSl22QDcH8V+4AJJS/HYA6bPX0T8JPMD8/z9zwV1/VwUESdz/S/ARS32O1fSk5L2S2oUjRcCf4+If+f2cWBZF/taR+3mDwBJV1Ou6jxbCX8u/0V1r6RzutTPulgG/KmyPdmYeXWfHFsvUMZaO23nu5nm4E7K1a6Gyc7jQdNuDm/N83JE0iUzbDuftZ2DnG60HBithD0Gp9Yqvx57M9f8/hfAE5KekrS5R32aNQt63YFBJGk38JZJHhqubkRESGr1MSxvi4gTkoaAUUkHKYXOvDdL+SOvMnwL2BQRr2T4bkohvpDyMT93AffMRr9tsEm6HVgJXFsJ/895HBHPTv4MA+0x4IGImJD0ccp/TNb2uE/9aCMwEhH/qcQ8Bq3rJL2HUlCvroRX59h7M7BL0tN5xbsvuaDugYi4odVjkp6TtDQiTmbBN97iOU7kz6OS9gJXAg9T/h21IK8kXgycmPUD6LHZyJ+k84HHgeH8N17juRtXtyck7QA+PYtdr6MTwCWV7cnGTGOf45IWAG8Enm+z7XzXVg4k3UD5g+/aiJhoxFucx4NWzEybw4h4vrK5jXKvRKPtdU1t9856D+ttJufhRuAT1YDH4LRa5ddjr02SVlDO2/XVc7ky9sYlPUqZvtS3BbWnfNTPTqBxt/Am4LvNO0ha1JiKIGkJ8G7gd1Fm+e8Bbpuq/TzXTv4WAo9S5sWNND22NH+KMv960juX55FfAJeqfDrMQsov3OY7/as5vQ0YzbG2E9io8ikgy4FLgZ/PUb/rYtr8SboS+Dpwc0SMV+KTnsdz1vP6aCeHSyubNwOHc/0HwI2Zy0XAjRkbJO2cw0i6jHLz3E8rMY/B6e0EPpKf9rEKeCEvvHjstUHSW4FHgA9HxJFK/DxJb2isU/LX379ve31XpJfXLpS5qT8EngF2A4szvhLYluvXAAcpd3MfBO6stB+iFDVjwEPAOb0+phrm73bgX8CvK8sV+dho5vQQ8G3g9b0+pjnI2U3AEcpVqeGM3UMpAAHOzbE0lmNrqNJ2ONv9nnL1oefHU8P87Qaeq4y1nRlveR4P2tJGDr8A/DZztQe4rNL2jhybY8DHen0sdcxfbm8FvtjUbuDHIOUGzZP5O+E4ZVrCFmBLPi7gq5nbg8DKSluPvenztw04XXn/ezLjQznufpPn9nCvj6XTxd+UaGZmZmbWAU/5MDMzMzPrgAtqMzMzM7MOuKA2MzMzM+uAC2ozMzMzsw64oDYzMzMz64ALajOzmsnPvP2xpPWV2Pslfb8HfdkraeVcv66ZWT/xNyWamdVMRISkLcBDkvZQ3qs/D6zr5utWvmXVzMxmwAW1mVkNRcQhSY8BdwHnUb7Z8zVfCS1pHaXQPhv4W0RcL2kxsJ3yxQlngM0RcWCK+Fbg7Rn/o6Q7gB3Au4Cngdfla50NfJPyJUkBbI+Ie7uZAzOzfuGC2sysvj4L/BJ4mVLIvkrSm4BvAGsi4lgWzI02v4qIWyStBe4HrpgiDvBOYHVEvCTpU8CZiHiHpBX5+uS+yyLi8nz9C7p0zGZmfccFtZlZTUXEi5K+A/wzIiaaHl4F7IuIY7nvqYyvBm7N2KikCyWdP0Ucytehv5Tra4Av5X4HJB3I+FFgSNKXgceBJ2b7eM3M+pVvSjQzq7dXcummF6fbISJOU6aB7AW2ANu63Cczs77hgtrMrD/tB9ZIWg5QmfLxI+BDGbuOMrf6H1PEm+0DPpj7XQ6syPUlwFkR8TDwGeCqrhyVmVkf8pQPM7M+FBF/lbQZeETSWcA48F5gK7A9p2qcATZlk1bxZl8Ddkg6DBwGnsr4sow3LsTcPbtHZGbWvxQRve6DmZmZmVnf8pQPMzMzM7MOuKA2MzMzM+uAC2ozMzMzsw64oDYzMzMz64ALajMzMzOzDrigNjMzMzPrgAtqMzMzM7MOuKA2MzMzM+vAfwHakQ0H2+/9LgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7fa65fb5b310>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Numerical flux at y =  -0.6 is 0.241786540077\n",
      "Exact     flux at y =  -0.6 is 0.298507462687\n",
      "\n",
      "Abs. error = 2.581e-04\n",
      "Rel. error = 6.156e-06\n"
     ]
    }
   ],
   "source": [
    "# create numpy arrays for analytics, parallel check first\n",
    "\n",
    "yvals       = np.zeros(len(mesh.specialSets['MinI_VertexSet']))\n",
    "ycoord      = np.zeros_like(yvals)\n",
    "analytic    = np.zeros_like(yvals)\n",
    "\n",
    "ids = mesh.specialSets['MinI_VertexSet']\n",
    "yvals[:] = tField.evaluate(ids).T\n",
    "\n",
    "ycoord = tField.mesh.data[ids.data,[1]]\n",
    "analytic = analyticTemperature(ycoord, h, k, c0, c1)\n",
    "\n",
    "abserr = uw.utils._nps_2norm(analytic - yvals)\n",
    "mag = uw.utils._nps_2norm(analytic)\n",
    "relerr = abserr / mag\n",
    "\n",
    "# measure border flux, analytic is easy, parallel check needed for numeric result\n",
    "yspot = y0\n",
    "ana_flux = exactDeriv(yspot,h,k,c0)\n",
    "\n",
    "tmp = tField.fn_gradient.evaluate_global([0.2,yspot])\n",
    "if tmp is not None: num_flux = tmp[0][1]\n",
    "else: num_flux = 0.\n",
    "\n",
    "\n",
    "from mpi4py import MPI\n",
    "comm = MPI.COMM_WORLD\n",
    "# assuming order in the allgather is the same\n",
    "coords = comm.allgather(ycoord)\n",
    "numerical = comm.allgather(yvals)\n",
    "\n",
    "if rank == 0:\n",
    "    # build matplot lib graph of result only on proc 0\n",
    "\n",
    "    # 1st build exact solution hiRes\n",
    "    big = np.linspace(y0,y1)\n",
    "    cool = analyticTemperature(big, h, k, c0, c1)\n",
    "\n",
    "    pylab.rcParams[ 'figure.figsize'] = 12, 6\n",
    "    pyplot.plot(coords, numerical, 'o', color = 'black', label='numerical') \n",
    "    pyplot.plot(big, cool, color = 'red', label=\"exact\") \n",
    "    pyplot.xlabel('Y coords')\n",
    "    pyplot.ylabel('Temperature')\n",
    "    pyplot.show()\n",
    "\n",
    "\n",
    "if rank == 0:\n",
    "    threshold = 1.0e-4\n",
    "    yspot = y0\n",
    "    print(\"Numerical flux at y = \" ,yspot,\"is\", num_flux)\n",
    "    print(\"Exact     flux at y = \" ,yspot,\"is\", ana_flux)\n",
    "    print(\"\\nAbs. error = {0:.3e}\".format(abserr))\n",
    "    print(\"Rel. error = {0:.3e}\".format(relerr))\n",
    "    if relerr > threshold:\n",
    "        raise RuntimeError(\"The numerical solution is outside the error threshold of the analytic solution.\" \\\n",
    "                           \"The Relative error was \", relerr,\" the threshold is \", threshold)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Model 2)\n",
    "\n",
    " * a fixed temperature condition or Dirichlet BC - bottomWall \n",
    " \n",
    " $ T(x,y_{0}) = T_{0} $\n",
    "\n",
    " * a heat flow condition or Neumann BC - topWall.\n",
    "\n",
    " $ q \\cdot n_{u} = (\\,0.0\\,,\\, f\\,) \\cdot (\\,0.0\\,,\\,1.0\\,) = f$\n",
    " \n",
    "     **Note** The top surface outward normal $n_{u}$ point along the j-axis \n",
    "\n",
    "------\n",
    "\n",
    "Arbitrary constants are:\n",
    "\n",
    "$c_{0} = \\frac{1}{k} (\\,f + hy_{1}\\,) $\n",
    "\n",
    "$c_{1} = T_{0} + \\frac{h}{2 k}y_{0}^2 - c_{0}y_{0} $"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [],
   "source": [
    "T0 = 8.0       # surface temperature\n",
    "k = 2.2        # diffusivity\n",
    "h = -7.4       # heat production, source term\n",
    "f = 4.0        # temperature flow, implies negative gradient\n",
    "\n",
    "# analytic solution definitions\n",
    "# 1 dirichlet conditions (top) + 1 neumann (bottom)\n",
    "c0 = 1.0*(f+h*y1)/k\n",
    "c1 = T0 + h/(2.0*k)*y0**2 - c0*y0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [],
   "source": [
    "for ii in bottomWall:\n",
    "    tField.data[ii] = T0\n",
    "\n",
    "# define neumann condition\n",
    "nbc = uw.conditions.NeumannCondition( fn_flux=f, \n",
    "                                      variable=tField, \n",
    "                                      indexSetsPerDof=(topWall) )\n",
    "\n",
    "# flag the dirichlet conditions on the topWall only\n",
    "bc = uw.conditions.DirichletCondition(tField, indexSetsPerDof=(bottomWall) )\n",
    "\n",
    "\n",
    "# define heat eq. system\n",
    "ss = uw.systems.SteadyStateHeat( temperatureField  = tField,\n",
    "                                 fn_diffusivity = k,\n",
    "                                 fn_heating     = h,\n",
    "                                 conditions     = [bc, nbc] ) "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [],
   "source": [
    "solver = uw.systems.Solver(ss)\n",
    "solver.solve()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAtQAAAF3CAYAAACFTdwtAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvNQv5yAAAIABJREFUeJzs3XmcneP9//HXZ7JghAiCWGaGUPsWU1TtSe2Vqi01ilJjKV2U0m/aqraD9tcN/aJTpfQ7lgotal9qX2OLrQgyWWjFXh2yXr8/rhPGyDLJ5Mw9y+v5eJzHOec+95y853Ymebvnuq8rUkpIkiRJWjQVRQeQJEmSujMLtSRJktQBFmpJkiSpAyzUkiRJUgdYqCVJkqQOsFBLkiRJHWChliRJkjrAQi1JkiR1gIVakiRJ6gALtSRJktQBfYsOsLBWXHHFVFNTU3QMSZIk9XCPPvroGymlwQvar9sV6pqaGsaOHVt0DEmSJPVwEdHcnv0c8iFJkiR1gIVakiRJ6gALtSRJktQBFmpJkiSpAyzUkiRJUgdYqCVJkqQOsFBLkiRJHWChliRJkjrAQi1JkiR1QFkLdUR8JyKeiYinI+KyiFiyzetLRMQVETE+Ih6KiJpy5llUTU1N1NTUUFFRQU1NDU1NTUVHkiRJUhdRtkIdEasB3wRqU0obAX2AUW12OwJ4O6W0NvAb4OflyrOompqaqK+vp7m5mZQSzc3N1NfXW6olSZIElH/IR19gqYjoC1QCr7Z5fSRwcenxGGB4RESZMy2U0aNH09LSwpeBOafXW1paGD16dJGxJEmS1EWUrVCnlKYAvwQmAq8B76aUbmmz22rApNL+M4F3gRXavldE1EfE2IgYO3Xq1HJFnquJEyeyMXAV8Ns22yVJkqRyDvkYRD4DvSawKrB0RBy8KO+VUmpMKdWmlGoHDx68OGMuUFVVFU8BZwJHAQe22i5JkiSVc8jHCOCVlNLUlNIM4Gpgmzb7TAHWACgNCxkIvFnGTAutoaGByspKfgjcB/wB2HjJJWloaCg4mSRJkrqCchbqicDWEVFZGhc9HHiuzT7XAoeWHu8H3JFSSmXMtNDq6upobGxktepqDgJmVVRw58orU7fvvkVHkyRJUhdQzjHUD5EvNHwMeKr0ZzVGxE8iYu/Sbn8EVoiI8cAJwCnlytMRdXV1TJgwgeaUWO5vf2P55mY48cSiY0mSJKkLiC52QniBamtr09ixY4sN8d3vwq9/DWPGgGeqJUmSeqSIeDSlVLug/VwpcVGccQZsuSUccQS8/HLRaSRJklQgC/Wi6N8fLr88Px41CqZPLzaPJEmSCmOhXlRrrgkXXQSPPAInn1x0GkmSJBXEQt0R++wDxx8Pv/0tXHtt0WkkSZJUAAt1R/2//wfDhsFhh0Fzc9FpJEmS1Mks1B21xBJwxRUwc2YeTz1jRtGJJEmS1Iks1IvD2mvDBRfAgw/CD35QdBpJkiR1Igv14nLAAXD00fCLX8ANNxSdRpIkSZ3EQr04/eY3sMkmcMghMHly0WkkSZLUCSzUi9OSS8Jf/gIffggHHZTHVUuSJKlHs1AvbuuuC7//PdxzD5x6atFpJEmSVGYW6nKoq4Ovfx1OPx1uvLHoNJIkSSojC3W5nH02bLopHHyw81NLkiT1YBbqcllqKRgzJo+j3n9/mDat6ESSJEkqAwt1Oa29Nlx8MTzyCJxwQtFpJEmSVAYW6nL70pfgxBPh3HPh0kuLTiNJkqTFzELdGU4/HbbbDo48Ep59tug0kiRJWows1J2hXz+4/HJYZhnYd194//2iE0mSJGkxsVB3llVXhcsugxdeyGeqUyo6kSRJkhYDC3Vn2mkn+NnP8tnqc88tOo0kSZIWAwt1Zzv5ZNhrL/jOd+Chh4pOI0mSpA6yUHe2ioo8ld5qq+X5qd98s+hEkiRJ6gALdRGWXx6uvBL+/e+8kuLs2UUnkiRJ0iKyUBeltjYvT37TTdDQUHQaSZIkLSILdZHq6/MZ6lNPhVtvLTqNJEmSFoGFukgRcP75sMEGcNBBMHly0YkkSZK0kCzURVt6abjqKvjwQzjgAJg+vehEkiRJWggW6q5g3XXhj3+EBx6A73636DSSJElaCBbqruKAA+CEE+B3v4NLLik6jSRJktrJQt2V/PznsOOOcNRR8NhjRaeRJElSO1iou5K+feGKK2DFFeHLX3bRF0mSpG7AQt3VrLQSXH01vPYar+24I2tVV1NRUUFNTQ1NTU1Fp5MkSVIbFuqu6LOf5cFDDmHI009z5MSJpJRobm6mvr7eUi1JktTFWKi7qFG33srvge8DXy5ta2lpYfTo0QWmkiRJUlsW6i5q4sSJfBN4EPgTsF6r7ZIkSeo6LNRdVFVVFdOB/YAW4K/AMqXtkiRJ6jos1F1UQ0MDlZWVTAH2B9YGmvr0oeGnPy04mSRJklqzUHdRdXV1NDY2Ul1dzb0RNAwaxBdnzaLOIR+SJEldStkKdUSsGxFPtLq9FxHfbrPPjhHxbqt9flSuPN1RXV0dEyZMYPbs2Zz65ptQVwc//CHcdFPR0SRJklTSt1xvnFJ6HtgMICL6AFPIQ4HbuieltFe5cvQYEdDYCE8/DQcdBGPHwlprFZ1KkiSp1+usIR/DgZdSSs2d9Of1TJWVedEXyCsptrQUm0eSJEmdVqhHAZfN47XPRcSTEXFjRGzYSXm6r7XWgksvhXHj4MgjIaWiE0mSJPVqZS/UEdEf2Bu4ci4vPwZUp5Q2Bc4B/jaP96iPiLERMXbq1KnlC9td7LYb/PSnuVifdVbRaSRJknq1zjhDvTvwWErp321fSCm9l1J6v/T4BqBfRKw4l/0aU0q1KaXawYMHlz9xd/D978M++8CJJ8LttxedRpIkqdfqjEL9FeYx3CMiVomIKD3espTnzU7I1P1VVMDFF8N668H++8NLLxWdSJIkqVcqa6GOiKWBLwBXt9p2dEQcXXq6H/B0RDwJnA2MSslBwe22zDJw7bV5BpCRI+E//yk6kSRJUq8T3a2/1tbWprFjxxYdo2u54w7YZRfYa688C0iF6/VIkiR1VEQ8mlKqXdB+Nq+eYOed4Te/gWuugVNPLTqNJElSr1K2hV3UyY47Dp58En72M9hkkzyuWpIkSWXnGeqeIgL+939hm23gsMPgiSeKTiRJktQrWKh7kiWWyGOol18+X6T4+utFJ5IkSerxLNQ9zcorw9/+lsv0fvvB9OlFJ5IkSerRLNQ90RZbwIUXwj33wDe/WXQaSZKkHs2LEnuqr3wlX6T485/DppvCMccUnUiSJKlH8gx1T9bQAHvumc9S33VX0WkkSZJ6JAt1T9anDzQ1wdpr5/HUEyYUnUiSJKnHsVD3dAMH5gVfZszIM3/8979FJ5IkSepRLNS9wWc+A1dcAU8/DYceCrNnF51IkiSpx7BQ9xa77gq/+AVcdRX8+MdFp5EkSeoxnOWjNznhBHj2WfjpT2HddaGuruhEkiRJ3Z5nqHuTCDjvPNhhBzj8cLj//qITSZIkdXsW6t6mf/887KOqCr70JWf+kCRJ6iALdW+0wgpw3XV55o+99oL33is6kSRJUrdloe6t1lsPxoyBf/4TRo2CmTOLTiRJktQtWah7s+HD4dxz4cYb4cQTi04jSZLULTnLR29XXw/PPQe//W0+a3300UUnkiRJ6lYs1IJf/hJefBGOOy4vUz5iRNGJJEmSug2HfAj69IHLLoMNNoD99svjqiVJktQuFmplyyyTZ/5YYok888ebbxadSJIkqVuwUOtj1dVwzTUweTL/3nZb1qmupqKigpqaGpqamopOJ0mS1CVZqPVJW2/NvUccwcr//CffnziRlBLNzc3U19dbqiVJkubCQq1POfj66/kxcDhwUmlbS0sLo0ePLi6UJElSF+UsH/qUiRMnchqwHnAm8CLwt9J2SZIkfZJnqPUpVVVVAHwNeAi4FNiy1XZJkiR9zEKtT2loaKCyspIPgZHAq8DfgbO++c1ig0mSJHVBFmp9Sl1dHY2NjVRXV/NGBF9fdVWWHTCAkb//vdPpSZIktWGh1lzV1dUxYcIEZs+ezT+mTGGJG2+E5mb40pfgww+LjidJktRlWKjVPttuCxdfDPfeC4cdBrNnF51IkiSpS3CWD7XfgQfms9Qnnww1NXDmmUUnkiRJKpyFWgvnpJPglVfg5z+HNdeEo44qOpEkSVKhLNRaOBFwzjkwaRIceyyssQbssUfRqSRJkgrjGGotvL594fLLYbPN4IAD4LHHik4kSZJUGAu1Fs2AAfD3v8MKK8Cee4KrKEqSpF7KQq1FN2QI3HADfPBBHvbxzjtFJ5IkSep0Fmp1zIYbwtVXwwsvwL77wvTpRSeSJEnqVBZqddzOO8MFF8Add8CRR0JKRSeSJEnqNGUr1BGxbkQ80er2XkR8u80+ERFnR8T4iBgXEcPKlUdldsghcNppcMkl+V6SJKmXKNu0eSml54HNACKiDzAF+Gub3XYH1indtgLOK92rO/rhD2HChFyoV1stn62WJEnq4TprHurhwEsppeY220cCl6SUEvBgRCwXEUNSSq91Ui4tThHw+9/Dv/4FRx8NK60EI0cWnUqSJKmsOmsM9SjgsrlsXw2Y1Or55NI2dVf9+sGVV0JtLYwaBffdV3QiSZKksip7oY6I/sDewJUdeI/6iBgbEWOnTp26+MKpPJZeGq6/HqqqYK+94Jlnik4kSZJUNp1xhnp34LGU0r/n8toUYI1Wz1cvbfuElFJjSqk2pVQ7ePDgMsXUYrXiinDzzbDUUrDbbnmpckmSpB6oMwr1V5j7cA+Aa4FDSrN9bA286/jpHqSmBm68Ed57D3bdFd56q+hEkiRJi11ZC3VELA18Abi61bajI+Lo0tMbgJeB8cAfgGPLmUcF2HRTuOYaeOkl2HvvvKqiJElSD1LWWT5SSv8FVmiz7fxWjxPwjXJmUBew447Q1AQHHJAvVLzqKujbWRPMSJIklZcrJapz7LcfnHMOXHstHHOMqylKkqQew9OE6jzf+Aa89ho0NMCQIfCTnxSdSJIkqcMs1OpcP/1pLtU//Wku1cccU3QiSZKkDrFQq3PNWU3x9dfzGeuVV4Yvf7noVJIkSYvMMdTqfH37whVXwNZbw0EHwV13FZ1IkiRpkVmoVYzKSrjuOlhzTRg5EsaNKzqRJEnSIrFQqzgrrJBXUxwwAHbZBV58sehEkiRJC81CrWJVVcGtt8KsWfCFL8DkyUUnkiRJWigWahVv/fXhppvgrbd4d6ut2HyNNaioqKCmpoampqai00mSJM2XhVpdwxZbcOu3vkX/V1/lD5MnMyAlmpubqa+vt1RLkqQuzUKtLuPIP/+Z/YBNgWuBJYGWlhZGjx5dbDBJkqT5sFCry5g4cSI3AIcA2wN/IU+UPnHixEJzSZIkzY+FWl1GVVUVAJcDxwJfBP4EVK+xRnGhJEmSFsBCrS6joaGByspKAH4PnALUATevsw6kVGQ0SZKkebJQq8uoq6ujsbGR6upqIoLLq6t5Zq+9+Mztt8MPflB0PEmSpLnqW3QAqbW6ujrq6uo+3pASHH00nH46DBoEJ55YXDhJkqS5aFehjojVgXVSSv+IiCWAviml/5Y3mgREwLnnwrvvwkknwXLLwde/XnQqSZKkjyywUEfE4cBxwEBgKFANnAuMKG80qaRPH7jkEnjvPaivh4EDYf/9i04lSZIEtG8M9TeBrYH3AFJKLwArlTOU9Cn9+8OYMbDttlBXl1dWlCRJ6gLaU6g/TClNn/MkIvoAUb5I0jxUVsJ118FGG8E++8CddxadSJIkqV2F+r6I+B6wZETsBFwB/L28saR5GDgQbrkFhg6FvfaC++4rOpEkSerl2lOovwf8B/gn8C3gdsC1oFWcFVeE22+H1VeH3XeHhx8uOpEkSerF5luoS8M7LkopnZdS2iel9KXS49mdlE+au5VXzqV68GDYdVd4/PGiE0mSpF5qvoU6pTQLWCsi+nVSHqn9VlsN7rgDll0WvvAFeOqpohNJkqReqD3zUL8E3BMR1wAfzT2dUjq7bKmk9qquzqV6++1hxAi46y5Yb72iU0mSpF6kPWOoJwK3ApXA4FY3qWsYOjSX6gjYeWcYP77oRJIkqRdZ4BnqlNIPOyOI1CHrrpvHVO+4IwwfDnffnc9eS5IklVl7Vkq8FUhtt6eUdilLImlRbbgh3Hor7LRTvt19d54JRJIkqYzaM4b6B60eLwnsC0wrTxypgzbbLM9TPXx4Hv5x110wZEjRqSRJUg/WniEfD7XZdFdEtN0mdR2f/SzceGOeTm/EiLyi4mCH/UuSpPJY4EWJEbFsq9tyETEcGNQJ2aRF9/nPw9//Dq+8kqfUe+utohNJkqQeqj1DPp4hj6EOYCbwCnBkOUNJi8WOO8I118AXv5jPVt96Kyy3XNGpJElSD9OeQr1WSmlG6w0R0Z6vk4r3hS/AVVfBPvvkx7fcAoP8BYskSVp82jMP9dzGSz+8uINIZbPnnnD11TBuXC7Vb79ddCJJktSDzLNQR8RKEbEpsFREbBwRm5Ru25IXeZG6j732yqX6qafyhYqOqZYkSYvJ/IZu7AkcDqwOnNtq+38AF3tR97PnnvC3v+XhH8OHw223wQorFJ1KkiR1c/Ms1Cmli4CLIuKAlNJfOjGTVD67755L9Ze+9HGpXnHFolNJkqRurD3zUP8lInYFNiQv7DJn++nlDCaVzW675dk/Ro7Mpfr22y3VkiRpkbVnHupzgUOBE4ClgIOBtcucSyqvXXeFa6+FF17IKypOnVp0IkmS1E21Z5aPbVNKBwFvppR+CGxFOwt1aSGYMRHxz4h4LiI+1+b1HSPi3Yh4onT70cJ/C9Ii2mUXuO46ePFF3hk2jC3WWIOKigpqampoamoqOp0kSeom2lOoP5xzHxGrlJ6v2s73Pwu4KaW0HrAp8Nxc9rknpbRZ6faTdr6vtHiMGMFt3/42/SdP5pLJkxmcEs3NzdTX11uqJUlSu7SnUN8QEcsBvwSeACYAC7xIMSIGAtsDfwRIKU1PKb2z6FGl8vj6ZZexJ7Am8A9gZaClpYXRo0cXG0ySJHUL8y3UEVEB3JhSeieldCW5c2ycUvqfdrz3msBU8kwhj0fEBRGx9Fz2+1xEPBkRN0bEhgv9HUgdNHHiRO4E9gCqyaV6ldJ2SZKkBZlvoU4pzQZ+3+r5Byml9q6I0RcYBpyXUtoc+C9wSpt9HgOqU0qbAucAf5vbG0VEfUSMjYixU714TItZVVUVAHcBuwNrkEv1Fqu2d2STJEnqzdoz5OMfETFyEd57MjA5pTRn6fIx5IL9kZTSeyml90uPbwD6RcSn5i9LKTWmlGpTSrWDBw9ehCjSvDU0NFBZmRf/vIdcqlcHbp85EzxLLUmSFqA9hfow4K8R8UFEvBURb0fEAs9Sp5T+BUyKiHVLm4YDz7beJyJWiYgoPd6ylOfNhfkGpI6qq6ujsbGR6upqIoJJ1dXc96MfseyHH8J228H48UVHlCRJXViklOa/Q0SfuW1PKc1a4JtHbAZcAPQHXga+BhxY+vrzI+I44BhgJvABcEJK6f75vWdtbW0aO3bsgv5oqeMeeyxPrdevX15RcUOH+EuS1JtExKMppdoF7regQl16s1HAWiml0yNidWDllNKjiyHnQrNQq1M98wyMGAEzZsAtt8CwYQv+GkmS1CO0t1C3Z6XE3wE7AV8tbWoBzu9YPKmb2HBDuOceWHrpvKLi/fP9BYokSeqF2jOGepuU0lGUFngpzfLRv6yppK5k7bVzqR48OA8BueOOohNJkqQupD2FekZpPuoEEBErALPLmkrqaqqq4O67oaYG9tgDrr++6ESSJKmLaE+h/l/gKmBwRJwG3Av8vKyppK5oyBC46y7YaCPYZx8YM6boRJIkqQtYYKFOKV0C/IC89PhbwP4ppcvLHUzqklZYAW6/HT77WTjwQLjkkqITSZKkgrXnDDVAH2AGMH0hvkbqmQYOzDN+7LQTHHoonHde0YkkSVKB2jPLx2jgMmBV8gJyl0bE98sdTOrSll4a/v532GsvOPZY+OUvi04kSZIK0rcd+xwCbJ5SagGIiAbgceCMcgaTurwll4Srr4aDD4aTToL334dTT4W8+KckSeol2lOoX2uzX9/SNkn9+sGll+Yz1qedBm++CWedBRWOjJIkqbdoT6F+C3gmIm4mT523C/BIRPwaIKV0QhnzSV1fnz5wwQUwaBD8+te5VP/pT9Df6dolSeoN2lOory/d5niwTFmk7quiIo+jXnllOPnkXKqvugoGDCg6mSRJKrMFFuqU0h87I4jU7UXA974HK64IRx4Jw4fnBWBWXLHoZJIkqYzaM8vHbhHxSES8HhFvRcTbEfFWZ4STuqXDD88XKz75JGy3HUyaVHQiSZJURu25cup3wFHAasBgYMXSvaR5GTkSbr4ZXn0VttkGnnuu6ESSJKlM2lOoJwNPpJRmpJRmzbmVO5jU7e2wQ16qfMYM2HZbeOihohNJkqQyaE+h/h5wXUScFBHfnHMrdzCpR9hsM7jvPlhuOdh553zWWpIk9SjtKdSnAbOA5chDPebcJLXH0KG5VK+zDnzxi3DZZUUnkiRJi1F7ps1bI6W0UdmTSD3ZKqvk4R977w11dXlaveOOKzqVJElaDNpzhvrmiNi57Emknm7gwDzkY++94fjj8zLlKRWdSpIkdVB7CvXhwG0R8b7T5kkdtOSSMGZMnlrvJz+B+nouveQSampqqKiooKamhqampqJTSpKkhdCeIR+uSiEtTn375qXKhwyBhgZWuPBC3pw9mwQ0NzdTX18PQF1dXbE5JUlSuyzwDHVpirz9gZNLj4cAm5U7mNSjRcDPfsYpyy/P8NmzuQtYpfRSS0sLo0ePLjKdJElaCO1ZKfF3wE7AV0ubWoDzyxlK6i1+8fbbfBH4DPAgsH5p+8SJE4sLJUmSFkp7xlBvk1I6CvgQIKX0FtC/rKmkXqKqqoqbgO3JP1T3AzuUtkuSpO6hPYV6RkRUAAkgIlYAZpc1ldRLNDQ0UFlZyePA1sCrwC1A0557FhtMkiS12zwLdUTMuWDxf4GrgMERcRpwL/DzTsgm9Xh1dXU0NjZSXV3NpAhGrb46b62/Pp8/91w44wyn1ZMkqRuINI9/sCPisZTSsNLjDYERQAC3pZSe7ryIn1RbW5vGjh1b1B8vld+0aXlavUsvhaOOgt/9Ls8MIkmSOlVEPJpSql3QfvP7VzrmPEgpPQM8sziCSVqAJZaAP/8ZqqvzWepJk+CKK2DAgKKTSZKkuZhfoR4cESfM68WU0q/LkEcSQEUFnH56LtXHHgs77ADXX5+XMJckSV3K/C5K7AMMAJaZx01SuR11FFx7LTz/PGy9NTz3XNGJJElSG/M7Q/1aSuknnZZE0tztuSfcdVe+/9zn8tLlI0YUnUqSJJXM7wx1zOc1SZ1piy3goYdgjTVgt93gvPOKTiRJkkrmV6iHd1oKSQtWXQ333w+7757HVR9/PMycWXQqSZJ6vXkW6tKKiJK6kmWWgb/9DU48MU+nt+ee8M47RaeSJKlXa89KiZK6kj594P/9P7jgArjjjjyuevz4olNJktRrWail7uqII+C22+D112GrrfKFi5IkqdNZqKXubIcd4OGHYaWV8swff/xj0YkkSep1LNRSdzd0KDz4IOy8M3z96/Dd78KsWUWnkiSp1yhroY6I5SJiTET8MyKei4jPtXk9IuLsiBgfEeMiYlg580g91sCBeSXF44+HX/8aRo6E994rOpUkSb1Cuc9QnwXclFJaD9gUaLvM2+7AOqVbPeDkutKi6tsXzj4bzj0XbroJPv95mDCh6FSSJPV4ZSvUETEQ2B74I0BKaXpKqe38XiOBS1L2ILBcRAwpVyapVzjmmFyoJ0+GLbeEu+8uOpEkST1aOc9QrwlMBS6KiMcj4oKIWLrNPqsBk1o9n1zaJqkjRozI46oHDYLhw/Oc1SkVnUqSpB6pnIW6LzAMOC+ltDnwX+CURXmjiKiPiLERMXbq1KmLM6PUc627bp4BZLfd8tjqww+HDz8sOpUkST1OOQv1ZGBySumh0vMx5ILd2hRgjVbPVy9t+4SUUmNKqTalVDt48OCyhJV6pIED4Zpr4NRT4U9/gu22g0mTFvhlkiSp/cpWqFNK/wImRcS6pU3DgWfb7HYtcEhpto+tgXdTSq+VK5PUK1VUwI9/nJcsf/552GILF4GRJGkxKvcsH8cDTRExDtgMOD0ijo6Io0uv3wC8DIwH/gAcW+Y8Uu81cmQeArL88jB8OGMPOYSa6moqKiqoqamhqamp6ISSJHVLkbrZhUq1tbVp7NixRceQuq/33mPSTjuxxmOPcTFwNPAhUFlZSWNjI3V1dQUHlCSpa4iIR1NKtQvaz5USpd5m2WXZ/o03+BFwKHAv+UKGlpYWRo8eXWw2SZK6IQu11As1T5rET4EvAmsDjwI7ABMnTiw0lyRJ3ZGFWuqFqqqqAPg7sCXwBnAbcOpyyzlftSRJC8lCLfVCDQ0NVFZWAvACsBVwQ58+nPr223DwwfD++4XmkySpO7FQS71QXV0djY2NVFdXExEsX13Nfy66CH72M7j88rxk+bNtZ7mUJElz4ywfkj7pjjvgK1/JZ6nPPx+++tWiE0mSVAhn+ZC0aHbeGZ54Ampr4ZBDoL4ePvig6FSSJHVZFmpJnzZkCNx+O3z/+/CHP8A228D48UWnkiSpS7JQS5q7vn3h9NPh73+H5mYYNgyuuqroVJIkdTkWaknzt+ee8PjjsMEGsN9+8O1vw/TpRaeSJKnLsFBLWrDqarj7bvjWt+Css2D77cFFYCRJAizUktqrf3/47W/hyivzlHqbbw7XX190KkmSCmehlrRw9tsPHnsM1lgD9torX7g4Y0bRqSRJKoyFWtLCW3tteOABOPJIOPPMPATklVeKTiVJUiEs1JIWzVJLQWMjXHZZHgKy2WZw6aVFp5IkqdNZqCV1zKhR8OSTsNFGUFcHhx4K//lP0akkSeo0FmpJHVdTA3fdBT/6Efzf/+ULFh9+uOhUkiR1Cgu1pMWjb1847TS48848T/XnP5/HV8+eXXQySZLKykItafHabrs8BGSfffIMIF95EvRzAAAekUlEQVT4AkyZUnQqSZLKxkItafEbNAiuuAIuuAAefBA23RSuuaboVJIklYWFWlJ5RMARR+Q5q6uq4Etfgm98Az74oOhkkiQtVhZqSeW17rp5zurvfhfOPRdqa7n+jDOoqamhoqKCmpoampqaik4pSdIis1BLKr8lloBf/hJuvpkPpkzhC//zP4xqbqYiJZqbm6mvr7dUS5K6LQu1pM6zyy5ss8wyXAOcCdwDrAO0tLQwevToYrNJkrSILNSSOtWTU6ZwAPAVYF3gCeA4YFJzc6G5JElaVBZqSZ2qqqoKgMuBjYA7gXOAu5dcEizVkqRuyEItqVM1NDRQWVkJwGvAnsCx/fuzJcDGG8NFF0FKBSaUJGnhWKgldaq6ujoaGxuprq4mIqiurubzF15Iv2efhWHD4PDDYeRI+Ne/io4qSVK7ROpmZ4Jqa2vT2LFji44hqRxmz4azzsorLC69NJx3HhxwQNGpJEm9VEQ8mlKqXdB+nqGW1HVUVMB3vgOPPw5Dh8KBB8JXvgJvvVV0MkmS5slCLanrWX99uP9++OlPYcwY2HBDuPrqolNJkjRXFmpJXVPfvvCDH8DDD8Mqq8C+++bba68VnUySpE+wUEvq2jbfPJfqM8+EG27IZ68vuMCZQCRJXYaFWlLX168fnHwyjBuXC/aRR8Lw4TB+fNHJJEmyUEvqRtZZB26/HRob4bHH8rzVv/gFzJxZdDJJUi9moZbUvVRU5DPUzz4Lu++ez1xvuWWeGUSSpAJYqCV1T6uummf+GDMmX6j42c/CKafABx8UnUyS1MtYqCV1b/vum89WH3YY/PznsMkmcOedRaeSJPUiFmpJ3d+gQXnmj9tvz6st7rRTXsJ86tSik0mSeoGyFuqImBART0XEExHxqfXCI2LHiHi39PoTEfGjcuaR1MPtvDM89VQeV/3nP8NnPpOXL581q+hkkqQerDPOUO+UUtpsPuug31N6fbOU0k86IY+knqyyMs9ZPW4cDBsGxx4LW28NjzxCU1MTNTU1VFRUUFNTQ1NTU9FpJUk9QN+iA0hSWay/Ptx2G1xxBZxwAmmrrfiwTx/+M3MmCWhubqa+vh6Aurq6YrNKkrq1cp+hTsAtEfFoRNTPY5/PRcSTEXFjRGxY5jySepMIGDUK/vlPLhgwgENnzuR54AgggJaWFkaPHl1wSElSd1fuQr1tSmkYsDvwjYjYvs3rjwHVKaVNgXOAv83tTSKiPiLGRsTYqV5kJGlhLbssR73/PpsDzwIXAPcDmwMTJ04sNJokqfsra6FOKU0p3b8O/BXYss3r76WU3i89vgHoFxErzuV9GlNKtSml2sGDB5czsqQeqqqqiqeBHYCvAmsCjwAXDRgAb79daDZJUvdWtkIdEUtHxDJzHgO7AE+32WeViIjS4y1Led4sVyZJvVdDQwOVlZUA/B+wLtDYty9fff99WHdduOiiPOWeJEkLqZxnqFcG7o2IJ4GHgetTSjdFxNERcXRpn/2Ap0v7nA2MSimlMmaS1EvV1dXR2NhIdXU1EcFy1dUs+6c/UTF2LAwdmuet/uxn4e67i44qSepmorv119ra2jR27KemtJakRTd7Nlx2WV66fPJk+PKX4Re/yEVbktRrRcSj85n6+SOulChJFRVQVwfPPw8/+QncdBNssAGcdBK8+27R6SRJXZyFWpLmqKyEH/4QXnwxF+xf/QrWXjuvtjhzZtHpJEldlIVaktpadVW48EIYOxY23DCvtrjppvnMtSRJbVioJWlehg2Df/wDrr4apk2D3XfPt2efLTqZJPV4TU1N1NTUUFFRQU1NDU1NTUVHmicLtSTNTwTss08u0b/6FTzwAGyyST5r/e9/F51OknqkpqYm6uvraW5uJqVEc3Mz9fX1XbZUW6glqT3694cTToDx4+GYY6CxMc8C8oMfwDvvFJ1OknqU0aNH09LSwl7AtUB/oKWlhdGjRxecbO4s1JK0MFZcEc45J5+x3msvaGiAtdbK0+y1tBSdTpJ6hKHNzTwAXAesB1SXtk+cOLG4UPNhoZakRfGZz8Dll8Njj8HnPgcnn/zxjCAzZhSdTpK6p4ceghEjuB1YDfg6sAHwYunlqqqqwqLNj4Vakjpi883h+uvzCotDh+ax1eutB01NLmUuSe01bhyMHAlbbw3jxjH24IPZdKml+CMwZ9LSyspKGhoaikw5TxZqSVocttsul+rrr4dlloGDD4bNNoNrr4WUutXV6pLUaV58EQ46KP99eddd8LOfwcsvU/vnP3POH/5AdXU1EUF1dTWNjY3U1dUVnXiuXHpckha32bPhyis/WiRm6jrr8NWJE7l52rSPdqmsrOzS/zhIUllNmpRXpr3oIlhiCfjWt+DEE2H55YtO9gkuPS5JRamogAMPhGeegcZGZr78MjdNm8atwPalXbry1eqSVDavvw7f+U6+5uSSS+Ab34CXXoLTT+9yZXphWKglqVz69YMjj2TorFmcAGwE3FW6jQAmNjcXGk+SOs2UKblI19TA2WfnYXEvvABnnQWrrFJ0ug6zUEtSma1UXc1vgDWB44G1gFuBR/v3z2Ouu9nQO0lqt1degaOPztOLnnMO7L9/nnb0j3+E6uoFf303YaGWpDJraGigsrKSD4HfAUOB4/v1Y52BA/Nc1rW18Ne/OiuIpJ7jn/+EQw+FddbJ46S/9rV8AeLFF8O66xadbrGzUEtSmdXV1dHY2PjR1epDqqvZ+qKLGDBlClx4Ibz3Hnz5y7DppnDFFTBrVtGRJWnRPPlkvoZkgw3yxdnHHw8vvwznnw9rrll0urJxlg9JKtrMmblINzTAc8/lszf/8z95Kqm+fYtOJ0kL9vDDecq7667LU4cedxx8+9uw0kpFJ+sQZ/mQpO6ib1+oq4Onn85ndJZcMv+qdN114Xe/g//+t+iEkvRpKeW5o3fZBbbaCu67D047DZqb86wd3bxMLwwLtSR1FRUVsN9+8PjjcM01+R+j44+HNdbIZ6xffbXohJJ6obYLU1168cVw+eW5RO+4Yx7m8YtfwIQJ8KMfwaBBRUfudA75kKSu7IEH4Fe/yhct9umTh4GccAJssknRyST1Ak1NTdTX19PS0sJA4EjgmxGskVK+4PDb384XHC61VNFRy6K9Qz4s1JLUHbz0Up6v9cIL8xCQL3whF+tdd4WIotNJ6qFqamqoaG7mW8ARwADgDqBp8GD++K9/5d+s9WCOoZaknmTo0LwYwqRJcMYZebz17rvDxhvnkt1qWXNJ6rCU4N57+VVzM+OBY4Crgc2B4cBFb7zR48v0wvBISFJ3MmgQnHJKHqt48cV5GMgRR+QFEn72M5g6teiEkrqzGTM+Hh+93XbsXFHBGUANcCjwRGm3qqqqwiJ2RRZqSeqO+veHQw6BJ56AW2+FzTeHH/4QVl89L+l7772Q0qcuJmpqaio6uaSu6PXX4cwz82/DvvIVePddOPdcbm5s5PTKSl5rtWtlZSUNDQ2FRe2KnOBUkrqzCBgxIt+efRbOOw8uuQSamnh7jTV45F//4q0ZM0hAc3Mz9fX1QF5sRlIvN2fau/PPh6uvzmend9oJzj0X9tgDKioYBcxacklGjx7NxIkTqaqqoqGhwb9D2vCiREnqad5/Hy67jKePO46Npk/nP0ATcB4wDqiurmbChAmFRpRUoLfeykPGfv97eP75PJTssMOgvh7WW6/odF2KFyVKUm81YAAceSSbTJ/OVsBV5LGPTwL3Ads3N8OHHxYaUVInSylPw3noobDaanmWoOWXz8V6yhT49a8t0x1goZakHqqqupqHga8BqwHfAVYALoE81vqkk2D8+AITSiq7997LQzg22wy22SbPaX/44Xkxlvvvz9di9NA5pDuThVqSeqiGhgYqKysBeBv4LTBsqaW47ZRT8upmv/lNXphhu+3gggvyRUiSur+U8jLgRx4Jq64K3/hGnhGosTGvuPq//+viUIuZhVqSeqi6ujoaGxuprq4mIqiurqbxD39gxBlnwJgxMHFintP6jTfyP7yrrJJXYrz5Zpg1q+j4khbWSy/BaafB2mvDttvCpZfCgQfCww/Do4/mn/MBA4pO2SN5UaIk9XYp5X9wL744zz/79tv5rNbBB+fxlhtsUHRCqddramqa+0wbb78NV16ZZ/e5774888/OO+ehHF/+sgW6g1x6XJK08KZNg+uuy+X6xhvzmerPfjYX61GjYIUVik4o9TpNTU3U19fT0tIC5DmPR/bvz68224zqJ5/MP7frr59/Tuvq8jUSWiws1JKkjvn3v/OvjC++OF/A1K8ffPGLedGHPfaA0vhsSeVVU1NDc3Mzw4BDgK8AKwFvVlSwwnHH5bPRw4bls9NarCzUkqTF58knc7G+9NJctCsrYc89Yf/9c7leeumiE0o9T0rw+OOcscUW7At8BpgGXAP8GbgZmN7Nelx34zzUkqTFZ9NN8zy1U6bAHXfkXy3fdRcccACstFIu1n/5C/z3vwAueS4tqpTgoYfytJZDh8IWW3ASMAGoB1YBDgT+DqxaXV1gULXm0uOSpPbr0ycvTbzTTnDOOXDPPfmCqKuuyjOHLLUUzRtvzE1PPMHU6dNd8lxqj9mz86IrY8bkn6VJk/IQqxEj4Ac/4K/Tp3PYd7/70RhqgMrKShoaGgoMrdYc8iFJ6rhZs+Dee+Evf2Hq+eczePZsWoAbyCs13gws65Ln0sfm/MyMGQNXX53nh+7fH3bdFfbbL1+vMGjQR7vPc5YPlZVjqCVJhegbweeB/YF9gSHALOB+YLszzshjrzfayAuo1Pu89RbccgvccAPcdBNMnQpLLpmvQ9hvv/yzseyyRadUKxZqSVIh5sxIAPlCnVpgT+DL/fuz0fTpeaeqqlwe9twzDx9xxhD1RCnlC3pvuCHfHnggD+9YYQXYbTfYe+9cpp0rusvqEhclRsSEiHgqIp6IiE+14MjOjojxETEuIoaVM48kqfxaL3k+G3gY+HllJU9eeGG+qPEPf8hTfF1yCey1Vy4Xe+4J554LpSIudVULvOD2vffyEI6vfx1WWw023xxGj4YPP8z3DzyQZ8r5v//LF/VapnuEsp6hjogJQG1K6Y15vL4HcDywB7AVcFZKaav5vadnqCWp62vXeM9p0/JMIddfn28vvZS3b7hhvhhrp51g++0/MY5UKlLbBVYABiy1FFeMHs0e/fvnxZDuuQdmzoSBA2GXXfIZ6N12g1VWKTC5FlWXGPLRjkL9e+DOlNJlpefPAzumlF6b13taqCWpB0oJXnghF+sbb8xLKH/wQR5nvfnmeSnlnXaC7baDZZYpOq16qZqaGiY1N7MJsGPptj3w0f/ybbxxLtB77AGf+1yeqUPdWpcY8gEk4JaIeDQi6ufy+mrApFbPJ5e2SZJ6kwhYd1044QS49VZ4++189vrUU3OBPvvsPCxk0CDYemv4/vfzxV2lea9bcw5sLVazZsHjj8NvfsNZzc28ATwO/AbYABgDHAwweTKMGwdnnpl/s2KZ7lXKfYZ6tZTSlIhYCbgVOD6ldHer1/8OnJlSurf0/Hbg5JTS2DbvU0+ez5yqqqotmh1jJ0m9ywcfwP33wz/+kW8PP5x/rd6vH2y1VT5zvdVWXDV5Mod873ufmq+3sbHRKcbUPtOn52J8zz1w551w993wzjsAvNK3L7fNnMmdwF3AlNKXVDslZI/VJYZ8fOIPivgx8H5K6ZettjnkQ5K08N5/Pw8LueOOXLAffzwXbPKKcg+1uj0GrGzh0dykBOPH55UJH3443x5/PJdqgLXXhh13zLcddqDprrs+NYba/2Hr2dpbqMu2UmJELA1UpJT+U3q8C/CTNrtdCxwXEZeTL0p8d35lWpIkIM+MsOuu+Qb5DPZjj/HdbbdlS2Br8vLMADOAJ5ub4RvfyGezt94a1lnHebB7o3/9Cx555OPy/PDDH519prISamvhm9+ELbfMY6BXX/0TXz6nNLvAitoq2xnqiFgL+GvpaV/g0pRSQ0QcDZBSOj8iAvgdsBvQAnyt7XCPtjxDLUmal9ZzYK9MPlOzFbDDkkvy+b5985ltyDMwbLJJvohs443z4402clGNbmS+M8nMmpVnjXn66XwbNy4X6YkT8+t9+uT/7ltu+fFt/fWhb9nOM6qb6nJDPhYXC7UkaV7mNq3ZR7+SHzUKnnsu/3r/kUfgqafy7T//+fgNqqs/Ltpz7j/zmXYVLZeG7jyt/zuvDmwEDOvXjyO22oq1/vvf/N/5ww/zzhEwdGg++zynPG++uYsJqV0s1JKkXmmhim1K+azluHG5XM+5f/75fJYToH//fPZyvfVyMRs6NI+tHToUhgyBior5F3lLdce9/z688ko+6/zyy1x66qlUv/8+GwEDW+32Wp8+DBkxIv+2Yc5t/fVh6aWLSq5uzkItSdKimjYtn+WccxZ73Lh88dqECR8XbYClloK11uKWl1/mqQ8+4CVgPPASMBFYrYdfDLnYzsrPng2vvgovv/zxrVSeeflleP31T+z+FvAU8HSb27sRzJ49u8PflzRH4RclSpLUbS2xBGy2Wb61NmNGPqP90ku5YL/0Erz0Eqs+8wzbAq0HEcwGXm9uzsMLVl01n80eMuTjx3PuV1mlw3MWFzHcpO1Z+ebmZurr85ITH/3Zs2fDG2/Aa6/lCwLn3Fo/f+21vOT8tGkfv3lFBVRVwVprwd575/u11sq/FVhrLYZtvjnNc8ZDt1JdVVXW71maF89QS5LUQTU1NUxsbmYVYG1gKFADrDNgAAftuGM++/raa/Dvf+eS2dbgwblcDxoEyy2XL5psz/3SS3PplVdy5DHH0PLBBx+93WIdbjJtGrz3Hrz7br4vPf7OEUcw/c03WZY87GLZ0m3VpZZi5/XXz9/v669/8oz+HAMGfPw/E6uskseul8oya62Vy3T//vOM5BAbdRaHfEiS1EnaXfBmzcolc07BnnM/5/b227m4vvPOxwW2nf9OTyNPETi9dD+rTx9WnVNM+/XL93375gwzZuR5u2fOnPvj1tvmVojbmA68C7wHvANssccenyzMq6zy8fOVV86FuoO8CFSdwUItSVInKkvBmz07z0Iyp2C3vn/nHWhp4YennEI/oB/Qv3TrBywBfK2uLpfj6dPz/YwZuVT37ZtLduv7eW1beuk8neDAgfm+9HiHL36R5159lffIZX4OVw1UT2KhliSpF2g993Zr5S62DrtQb9DeQl3RGWEkSVJ5NDQ0UNlmTuXKykoaGhrK+ufW1dXR2NhIdXU1EUF1dbVlWr2WZ6glSermHE8slYdDPiRJkqQOcMiHJEmS1Aks1JIkSVIHWKglSZKkDrBQS5IkSR1goZYkSZI6wEItSZIkdYCFWpIkSeoAC7UkSZLUARZqSZIkqQMs1JIkSVIHdLulxyNiKtBcdI6CrAi8UXSIbszj1zEev47x+HWMx69jPH4d4/HrmO58/KpTSoMXtFO3K9S9WUSMbc968po7j1/HePw6xuPXMR6/jvH4dYzHr2N6w/FzyIckSZLUARZqSZIkqQMs1N1LY9EBujmPX8d4/DrG49cxHr+O8fh1jMevY3r88XMMtSRJktQBnqGWJEmSOsBC3cVExPIRcWtEvFi6HzSP/WZFxBOl27Wttq8ZEQ9FxPiIuCIi+nde+uK15/hFxGYR8UBEPBMR4yLiwFav/SkiXml1bDfr3O+g80XEbhHxfOkzc8pcXl+i9FkaX/ps1bR67ful7c9HxK6dmburaMfxOyEini191m6PiOpWr83157i3accxPCwiprY6Vl9v9dqhpZ/3FyPi0M5N3jW04/j9ptWxeyEi3mn1Wq/+DEbEhRHxekQ8PY/XIyLOLh3bcRExrNVrfvYWfPzqSsftqYi4PyI2bfXahNL2JyJibOelLpOUkrcudAN+AZxSenwK8PN57Pf+PLb/BRhVenw+cEzR31NXO37AZ4B1So9XBV4Dlis9/xOwX9HfRycerz7AS8BaQH/gSWCDNvscC5xfejwKuKL0eIPS/ksAa5bep0/R31MXPH47AZWlx8fMOX6l53P9Oe5Nt3Yew8OA383la5cHXi7dDyo9HlT099TVjl+b/Y8HLmz1vFd/BoHtgWHA0/N4fQ/gRiCArYGHStt7/WevncdvmznHBdh9zvErPZ8ArFj097C4bp6h7npGAheXHl8MfKm9XxgRAewMjFmUr+8hFnj8UkovpJReLD1+FXgdWOCk7T3UlsD4lNLLKaXpwOXkY9ha62M6Bhhe+qyNBC5PKU1LKb0CjC+9X2+ywOOXUvpHSqml9PRBYPVOztjVteczOC+7AremlN5KKb0N3ArsVqacXdXCHr+vAJd1SrJuIKV0N/DWfHYZCVySsgeB5SJiCH72gAUfv5TS/aXjAz387z8LddezckrptdLjfwErz2O/JSNibEQ8GBFzSuMKwDsppZml55OB1cqYtStq7/EDICK2JJ/VeanV5obSr6h+ExFLlClnV7EaMKnV87l9Zj7ap/TZepf8WWvP1/Z0C3sMjiCf7Zpjbj/HvU17j+G+pZ/LMRGxxkJ+bU/W7mNQGm60JnBHq81+BudvXsfXz97Ca/v3XwJuiYhHI6K+oEyLTd+iA/RGEXEbsMpcXhrd+klKKUXEvKZhqU4pTYmItYA7IuIpctHp8RbT8aN0luHPwKEppdmlzd8nF/H+5Gl+TgZ+sjhyq3eLiIOBWmCHVps/9XOcUnpp7u/Qq10HXJZSmhYRR5F/Y7JzwZm6o1HAmJTSrFbb/Ayq7CJiJ3Kh3rbV5m1Ln72VgFsj4p+lM97dkoW6ACmlEfN6LSL+HRFDUkqvlQrf6/N4jyml+5cj4k5gc+Aq8q+j+pbOJK4OTFns30DBFsfxi4hlgeuB0aVf48157zlnt6dFxEXAiYsxelc0BVij1fO5fWbm7DM5IvoCA4E32/m1PV27jkFEjCD/D98OKaVpc7bP4+e4t5WZBR7DlNKbrZ5eQL5WYs7X7tjma+9c7Am7toX5ORwFfKP1Bj+DCzSv4+tnr50iYhPyz+3urX+WW332Xo+Iv5KHL3XbQu2Qj67nWmDO1cKHAte03SEiBs0ZihARKwKfB55NeZT/P4D95vf1PVx7jl9/4K/kcXFj2rw2pHQf5PHXc71yuQd5BFgn8uww/cn/4La90r/1Md0PuKP0WbsWGBV5FpA1gXWAhzspd1exwOMXEZsDvwf2Tim93mr7XH+OOy1519GeYzik1dO9gedKj28Gdikdy0HALqVtvUl7foaJiPXIF8890Gqbn8EFuxY4pDTbx9bAu6UTL3722iEiqoCrga+mlF5otX3piFhmzmPy8eve/94WfVWkt0/eyGNTbwdeBG77/+3dT4hXVRTA8e/RRUkRYbly1bSpEJWYRYQMUQTWKog2GYQGMus2EdNialHuAiVaNI0QbUJqUQQRoaItWlTajKLGkBBB9IeEyJESOi7eMV4DjsXz5+89/H7gwZvz7uW+e7m/3+/wuG8usL7ik8BcnT8ILNK8zb0IPNeqP0GT1CwBB4Cbxt2nHo7fM8BF4Hjr2FrXDtaYngDeBW4dd5+uw5g9DnxL81RqpmKv0CSAADfXXFqquTXRqjtT9c7QPH0Ye396OH6fAT+15tqHFb/i5/hGO/7DGL4GnKyxOgTc06q7q+bmErBz3H3p4/jV37PAnhX1bvg5SPOC5o/1m/ADzbKEaWC6rgfwRo3tIjDZquvcu/r4zQHnWt9/X1Z8oubdN/XZnhl3X7oe7pQoSZIkdeCSD0mSJKkDE2pJkiSpAxNqSZIkqQMTakmSJKkDE2pJkiSpAxNqSeqZ+p+3n0fEY63YUxHxyRju5XBETF7vdiVpSNwpUZJ6JjMzIqaBAxFxiOa7+lVg+yjbbe2yKkn6H0yoJamHMvNERHwEvADcQrOz57+2hI6I7TSJ9lrg18x8JCLWA/M0GycsA7szc2GV+Cxwd8W/j4hdwH5gC3AaWFdtrQXeptkkKYH5zHx9lGMgSUNhQi1J/fUy8DXwF00i+4+I2AC8BUxl5tlKmC/XOZaZT0TEw8A7wNZV4gD3Adsy80JEPA8sZ+a9EbG52qfKbszMTdX+7SPqsyQNjgm1JPVUZp6PiPeAPzLzzxWXHwCOZObZKvtbxbcBT1bsYETcERG3rRKHZjv0C3U+BeytcgsRsVDx74CJiNgHfAx8eq37K0lD5UuJktRvf9cxSuevViAzz9EsAzkMTANzI74nSRoME2pJGqYvgKmIuAugteTjKLCjYg/RrK3+fZX4SkeAp6vcJmBznd8JrMnM94GXgPtH0itJGiCXfEjSAGXmLxGxG/ggItYAPwOPArPAfC3VWAaerSpXiq/0JrA/Ik4Bp4CvKr6x4pcfxLx4bXskScMVmTnue5AkSZIGyyUfkiRJUgcm1JIkSVIHJtSSJElSBybUkiRJUgcm1JIkSVIHJtSSJElSBybUkiRJUgcm1JIkSVIHlwCkHClrFeopvwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7fa65f970890>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Numerical flux at y =  1.3 is -4.41274249682\n",
      "Exact     flux at y =  1.3 is -4.57272727273\n",
      "\n",
      "Abs. error = 1.252e-03\n",
      "Rel. error = 4.705e-05\n"
     ]
    }
   ],
   "source": [
    "# create numpy arrays for analytics\n",
    "yvals       = np.zeros(len(mesh.specialSets['MinI_VertexSet']))\n",
    "ycoord      = np.zeros_like(yvals)\n",
    "analytic    = np.zeros_like(yvals)\n",
    "\n",
    "ids = mesh.specialSets['MinI_VertexSet']\n",
    "yvals[:] = tField.evaluate(ids).T\n",
    "\n",
    "ycoord = tField.mesh.data[ids.data,[1]]\n",
    "analytic = analyticTemperature(ycoord, h, k, c0, c1)\n",
    "\n",
    "abserr = uw.utils._nps_2norm(analytic - yvals)\n",
    "mag = uw.utils._nps_2norm(analytic)\n",
    "relerr = abserr / mag\n",
    "\n",
    "# measure border flux, analytic is easy, parallel check needed for numeric result\n",
    "yspot = y0\n",
    "ana_flux = exactDeriv(yspot,h,k,c0)\n",
    "\n",
    "tmp = tField.fn_gradient.evaluate_global([0.2,yspot])\n",
    "if tmp is not None: num_flux = tmp[0][1]\n",
    "else: num_flux = 0.\n",
    "\n",
    "from mpi4py import MPI\n",
    "comm = MPI.COMM_WORLD\n",
    "# assuming order in the allgather is the same\n",
    "coords = comm.allgather(ycoord)\n",
    "numerical = comm.allgather(yvals)\n",
    "\n",
    "if rank == 0:\n",
    "    # build matplot lib graph of result only on proc 0\n",
    "\n",
    "    # 1st build exact solution hiRes\n",
    "    big = np.linspace(y0,y1)\n",
    "    cool = analyticTemperature(big, h, k, c0, c1)\n",
    "\n",
    "    pylab.rcParams[ 'figure.figsize'] = 12, 6\n",
    "    pyplot.plot(coords, numerical, 'o', color = 'black', label='numerical') \n",
    "    pyplot.plot(big, cool, color = 'red', label=\"exact\") \n",
    "    pyplot.xlabel('Y coords')\n",
    "    pyplot.ylabel('Temperature')\n",
    "    pyplot.show()\n",
    "\n",
    "\n",
    "if rank == 0:\n",
    "    threshold = 1.0e-4\n",
    "    yspot = y1\n",
    "    print(\"Numerical flux at y = \" ,yspot,\"is\", num_flux)\n",
    "    print(\"Exact     flux at y = \" ,yspot,\"is\", ana_flux)\n",
    "    print(\"\\nAbs. error = {0:.3e}\".format(abserr))\n",
    "    print(\"Rel. error = {0:.3e}\".format(relerr))\n",
    "    if relerr > threshold:\n",
    "        raise RuntimeError(\"The numerical solution is outside the error threshold of the analytic solution.\" \\\n",
    "                           \"The Relative error was \", relerr,\" the threshold is \", threshold)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "## Model 3) \n",
    "\n",
    "2D, Steady State Heat Equation with Dirichlet BC at the top and bottom surfaces.\n",
    "\n",
    "$T(x,y_{1}) = T_{1}$\n",
    "\n",
    "$ T(x,y_{0}) = T_{0} $\n",
    "\n",
    "------\n",
    "\n",
    "arbitrary constants are:\n",
    "\n",
    "$ c_{0} = \\frac{1}{y_{1}-y_{0}} \\left[ T_{1}-T_{0}+\\frac{h} {2k}(y_{1}^2-y_{0}^2) \\right] $\n",
    "\n",
    "$c_{1} = T_{1} + \\frac{h}{2k}y_{1}^2 - c_{0}y_{1}$\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Model parameters\n",
    "T1 = 8.0   # top surface temperature\n",
    "T0 = 4.0   # bottom surface temperature\n",
    "k = 0.50   # diffusivity\n",
    "h = 10     # heat production, source term\n",
    "\n",
    "# arbitrary constant given the 2 dirichlet conditions\n",
    "c0 = (T1-T0+h/(2*k)*(y1**2-y0**2)) / (y1-y0)\n",
    "c1 = T1 + h/(2*k)*y1**2 - c0*y1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "# set boundary conditions\n",
    "for ii in topWall:\n",
    "    tField.data[ii] = T1\n",
    "for ii in bottomWall:\n",
    "    tField.data[ii] = T0\n",
    "\n",
    "# flag boundary conditions\n",
    "bc = uw.conditions.DirichletCondition(tField, indexSetsPerDof=(topWall+bottomWall) )\n",
    "\n",
    "# define heat eq. system\n",
    "ss = uw.systems.SteadyStateHeat( temperatureField  = tField,\n",
    "                                 fn_diffusivity = k,\n",
    "                                 fn_heating     = h,\n",
    "                                 conditions     = [bc] )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "solver = uw.systems.Solver(ss)\n",
    "solver.solve()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [],
   "source": [
    "# create numpy arrays for analytics\n",
    "yvals       = np.zeros(len(mesh.specialSets['MinI_VertexSet']))\n",
    "ycoord      = np.zeros_like(yvals)\n",
    "analytic    = np.zeros_like(yvals)\n",
    "\n",
    "ids = mesh.specialSets['MinI_VertexSet']\n",
    "yvals[:] = tField.evaluate(ids).T\n",
    "\n",
    "ycoord   = tField.mesh.data[ids.data,[1]]\n",
    "analytic = analyticTemperature(ycoord, h, k, c0, c1)\n",
    "\n",
    "abserr = uw.utils._nps_2norm(analytic - yvals)\n",
    "mag = uw.utils._nps_2norm(analytic)\n",
    "relerr = abserr / mag\n",
    "\n",
    "# measure border flux, analytic is easy, parallel check needed for numeric result\n",
    "yspot = y0\n",
    "ana_flux = exactDeriv(yspot,h,k,c0)\n",
    "\n",
    "tmp = tField.fn_gradient.evaluate_global([0.2,yspot])\n",
    "if tmp is not None: num_flux = tmp[0][1]\n",
    "else: num_flux = 0.\n",
    "\n",
    "from mpi4py import MPI\n",
    "comm = MPI.COMM_WORLD\n",
    "# assuming order in the allgather is the same\n",
    "coords = comm.allgather(ycoord)\n",
    "numerical = comm.allgather(yvals)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Numerical flux at y =  -0.6 is 20.1552435024\n",
      "Exact flux at y= -0.6 is 21.10526315789474\n",
      "\n",
      "Abs. error = 1.988e-04\n",
      "Rel. error = 3.573e-06\n"
     ]
    }
   ],
   "source": [
    "if rank == 0:\n",
    "    threshold = 1.0e-4\n",
    "    yspot = y0\n",
    "    print(\"Numerical flux at y = \" ,yspot,\"is\", num_flux)\n",
    "    print(\"Exact flux at y=\" ,yspot,\"is\", ana_flux)\n",
    "    print(\"\\nAbs. error = {0:.3e}\".format(abserr))\n",
    "    print(\"Rel. error = {0:.3e}\".format(relerr))\n",
    "    if relerr > threshold:\n",
    "        raise RuntimeError(\"The numerical solution is outside the error threshold of the analytic solution.\" \\\n",
    "                           \"The Relative error was \", relerr,\" the threshold is \", threshold)  \n",
    "        "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAtEAAAF3CAYAAABjZBdpAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvNQv5yAAAIABJREFUeJzs3Xd4leX9x/H3HZbGPdC6kmAdtXVVUdFaF4q4ilp3HLjitq3+nKm70apttW5SUbQet9a9cO8qCI6KoiIB1Cq4MQpI7t8f96FGBEkgyXPOyft1XefKyXNOkk+eK4Rv7vN9vneIMSJJkiSp5cqyDiBJkiQVG4toSZIkqZUsoiVJkqRWsoiWJEmSWskiWpIkSWoli2hJkiSplSyiJUmSpFayiJYkSZJaySJakiRJaiWLaEmSJKmVumYdoCWWXHLJWFVVlXUMSZIklbjhw4dPijH2nNPziqKIrqqqYtiwYVnHkCRJUokLITS05Hm2c0iSJEmtZBEtSZIktZJFtCRJktRKFtGSJElSK1lES5IkSa1kES1JkiS1kkW0JEmS1EoW0ZIkSVIrWURLkiRJrWQRLUlqE7lcjqqqKsrKyqiqqiKXy2UdSZLaTVFs+y1JKmy5XI6amhoaGxsBaGhooKamBoDq6uoso0lSu3AlWpI0d2KEr76CDz9k0PHHs1JjIxsCWwH9gA0aG7nj2GPh+edh2DB4+WV4/XV46y14912YMAE+/BA+/hi+/nqeorgKLqmjuRItSSUml8tRW1vLuHHjqKiooK6ubs6rwTHCp5+mwnbCBBg//rv7H3wAkyf/8NbYmD4OeHJ2n/fDD2HDDVsWvLwcevaEpZZKb2e+NT++1FKwwAL/+35dBZfU0ULM/wIsZL17947Dhg3LOoYkFbyZC0qA8vJy6gcNorpvX3j1VRgz5oeF8oQJqShurqwMll0WllkGFl4YFlwwFa4LLvj92wILcMyppzL244+ZDHwFRKAbsNxSS3H9NdfAtGnw7bfp7Yxb8/cbG2HixO/fPvoovf3mm1l/sz17wiqrcPPIkYz46itGA6OBt4FvgMrKSsaOHdv2J1lSSQshDI8x9p7T81yJlqQSUltbS1ljI32ANWbcGhtZc7/9oKnpuyd26ZIK5OWXh7XXhu23T/dXWCG9XX55+MlPoGvL/ptYd5FFGDSL4r3mb3+D/v3n/hua0TLSvKieODGtcL/zDowezcZffcVuM31YA/BWQwMcfjisskq6/fznUFkJIcx9HknKs4iWpGLV1ASjRsErr6QV5ldf5bGGBno1e8qXwGvAbU1NHHzRRbDGGrDyyqlA7tKlzaLMaJtodRvJnITw3Yp3r16zfMpGVVV83NDAysAq8L+3a3TvDjfcAJ999t2Te/aE9deHDTZIb9dbDxZffN4ySuqUbOeQpGIxdSoMHw5PPZVuTz/9XYHYtSusuip3vfsuzzc28iqpeG4gtVaUcmvDbFtY6uup3muvdOHim2+mPzZeeCHdRo36Xz83K62UCuoZxfXaa8N882X03UjKmu0cklTsvvoKnnvuu6L5+ee/m2Kx6qqwyy6w8cbwy1/Cz34G3bvzZS7H32dRUNbV1WX0TbS/Oa6CL7lkuv3qV3DYYenYF1+kiSEziuonnoDrr0+Pde0Ka62VCuottoC+fWHRRTP4ziQVMleiJalQfPZZKuZmFM0vvZQuvisrS6ujv/41bLJJKpyXWmq2n2aupnMI3nvvu6J6xm3y5NT2ssEGsPXW0K9fagHJt8J4rqXS09KVaItoScrSmDFw991w113w5JOpaO7RI7UWzCiaN9wwTcdQx5o2La3+P/hgug0fnlpAFlsMttyS5xdZhH2vu463mk0P+V8biYW0VLQsoiUpQ7NdoWxqSiucd92ViufXXksf8POfw29+A9tum1Y67cktPJMmwcMPp4L6oYfg/fcBeB14MH97Ali6hPvPpc7AIlqSMjLzhW7lwHbdu3P2hhuy0htvpPFsXbqkVebf/AZ22AF++tNsQ6t1YmSNsjL6kXZn3ASYnzQN5U5g77vuSq0fPXpkmVLSXLCIlqSMVFVV8VlDA7sAA4AtSQXWFyGw8G67pcJ5m21SW4CKVlVVFQ0NDQDMB2wK7AzsWlbGYk1NsMgisNNOsPvu6eLEbt2yjCuphVpaRJd1RBhJ6hSmTYN77+W8hgb+C1wJrA7UA32BnjHCjTfCXntZQJeAuro6ysvLgbRD4oPAH8rLuf+qq+C++2DHHeH229MfTMssA4ccAo8+CtOnZ5pbUtuwiJakeTVyJBxzTNrlb/vt2bKsjHpgPWBF4PfAo8AylZWZxlTbqq6upr6+nsrKSkIIVFZWUl9fz1777ZcK5yFDUuvOHXek1o5cLq1IL7ccHHlkmvPdfBdJSUXFdg5JmhsffJDmCl97bdrEo1u31Nu8777c8NlnHHT44bPe/MOpDZ1XY2Naob7xRrj3Xvjmm7TN+oEHwkEHpeJaUuZs55Cktvb116kA2mabtOr8f/8H888Pl12WiurbboMBA9hzv/1muUJpAd3JlZenDXJuvRU++iitTK+2Gpx+OlRWpvaPBx5wdVoqEq5ES9KcjB4NF1+cVp2/+AIqKmCffdJt1VWzTqdiN2YM/OMfMHgwTJwIvXrBwQfDAQfA0ktnnU7qdJzOIUnzoqkJhg6Fv/8d7r8fundPUxb23x823TTtIii1palT4V//gkGD4LHHUovQTjvBoYfCZptBCFknlDoF2zkkaW5MnpzaM37xC+jfH0aMgDPOgHHj0kr05ptbQKt9zPhD7dFHYdSodPHh0KGwxRbws5/B3/4GH3+cdUpJef5PIEkA776bepyXXx6OOAIWWgiuuw4aGuDUU31ZXR1rRtH83nvpj7cll4Rjj00XHx50EHf95S9UVVVRVlZGVVUVuVwu68RSp9M16wCSlJkY4YknUsvGXXell8t32QV+9zvo08eXz5W9+ef/rv/+1Vfhssv4dvBgtp82jSnAn4GXGhqoqakB8OJVqQPZEy2p85k2LU1GuOCCNJ5uiSWgpgYOPzytREsFrPcKK7DThAkcASwKDAXOAcZUVDA2v4OipLlnT7QkAblc7n8ve69UWcnzBx+cJmrsv3+6ePAf/4Dx4+Hssy2gVRReeu89/ghUAMeTdsV8FLh53Li0Q6Ij8qQOYREtqWTlcjlqamp4v6GBA2Nk6Lhx9LnySj4GuPvutAp90EHpJXOpSFRUVADwJXA+0AuoAZbq2hV++1v4+c/h6qvTtA9J7cYiWlLJOv3kk6lubGQ08A9gIrAtsO706bD99vY8qyjV1dVRXl7+v/enALnycp696iq46aa0qcsBB8CKK6aWpcmTswsrlTCLaEmlZ8oUuOIKHh43jnrgQ2AbYAPgfmDc+PGZxpPmRXV19Sx3xNxrn31gt91g+PC08+HKK8Mxx0BVVSqmp0zJOrpUUrywUFLpmDIFrroKzjkHxo/npe7dOXnqVB6c6WmVlZWMHTs2i4RSx3r++TSicejQtNPmmWfC3ntDly5ZJ5MKlhcWSuo8pk5NG6SstNJ3EzYefJBRgwfzVLOXvQHKy8upq6vLKKjUwfr0gYcegocfhqWWgoEDYa210kjHIlhEkwpZuxXRIYSrQggfhRBem8Vjx4YQYghhyfb6+pI6gRjTNsm/+EXaIKWyMhUMzzwD/fpRvffes3zZ21m66nT69oUXXoBbbkkjHgcMgF//Gp5+OutkUtFqz5XoIUD/mQ+GEFYA+gHj2vFrSyp1w4bBppvCzjtDjx5w333w1FOw1Vbfu2CwurqasWPH0tTUxNixYy2g1XnN2Ezotddg0CAYMyYV0jvskDZykdQq7VZExxifBD6ZxUMXkEZb+jqSpNYbPx723RfWWw/eeAOuuAJGjoRttnHahtQS3bqlzYXefhv+/Oe0Gr3WWunfldcKSC3WoT3RIYQBwHsxxpc78utKKgGTJ8Mpp8Aqq8DNN8NJJ6Ui4JBDoGvXrNNJxae8HE44Ia1IH3dcavVYZZW07f3HH2edTip4HVZEhxDKgZOBU1v4/JoQwrAQwrCJEye2bzhJhWv6dBg8OI3r+tOfYKed4M030w6DCy+cdTqp+C22GJx7bvqjdOBAuPTSVEwPGpT+/UmapY5cif4paWOll0MIY4HlgZdCCD+Z1ZNjjPUxxt4xxt49e/bswJiSCsbDD8M666RdBXv1gueeg+uvTxcQSmpbyy0H9fWpPWqNNeDQQ9N0jxdeyDqZVJA6rIiOMb4aY1wqxlgVY6wCJgDrxBj/21EZJBWJN95IOwputRV88UXahe2ZZ9J/6JLa1+qrw2OPQS4H772X/t0dfDBMmpR1MqmgtOeIuxuA54BVQwgTQggHttfXklQivvkmbQyx5ppp0sZ558GoUWkXNi8alDpOCLDXXukP2mOOgauvTi0eV1zB9f/8J1VVVZSVlVFVVUUul8s6rZQJdyyUVBgefTS9fPzWW2lHtb/+NW0OISl7//kPHHkkPP44I0LgsBj5d/6h8vJy56+rpLhjoaTiMHEi7Ldf2gyiqSltT/zPf1pAS4XkF7+ARx/lqCWXZKkYeR64ElgSaGxspLa2NuOAUseziJaUjRhhyBBYbTW44QaorU0bPmy5ZdbJJM1KCFz68cf8DDgP2BcYDRwGTGhoyDSalAWLaEntLpfLfa+H8q7zz4fNN4f994ef/QxGjEjj6+afP+uokn5ERUUFk4ETgLWAl4DLgJe6d3fXQ3U6FtGS2lUul6OmpoaGhga6xcjAhga2Pv54prz4Yhqn9eST6aViSQWvrq6O8vJyAEYBWwL7du/OKvPNB+uum/4YnjYt04xSR7GIltSuamtraWxsZFPgFeB04FbgV4svnsZmlflrSCoW1dXV1NfXU1lZSQiByspKtr7qKuZ7+2347W/TrqIbbAAvuzGxSp/TOSS1q8VD4K/A/sAYUv/kQ0AIgaampkyzSWpjt98Ohx0Gn3wCf/wjnHQSdO+edSqpVZzOISl7Dz/Mf7p0YR/gHGB1UgENqbdSUonZeWd4/fU02/3002H99dM1D1IJsoiW1PYaG+Hoo2GrrShfemk269GDk4Gv8w+Xl5dTV1eXZUJJ7WWJJdJuh3fcAR9+mArpU0+FqVOzTia1KYtoSW3rxRdhnXXg4ovh6KNZ5O23OWzw4O/1ULoxg9QJDBiQNmnZc0846yzo3RuGD886ldRm7ImW1DamTYO6unR1/jLLpBnQfftmnUpSIbjnHjjkkLQyfcIJaWW6R4+sU0mzZE+0pI7zxhuw0UZwxhlp1enVVy2gJX1n++3htddgn33g7LPTODwneKjIWURLmntNTXDRRfDLX8KYMXDLLWnL7kUXzTqZpEKz2GJw9dVw331pesf666e2ryJ4RVyaFYtoSXNn/HjYemv43e/S7oOvvQa77JJ1KkmFbptt0ir0llumC5AHDIBJk7JOJbWaRbSk1okxXXm/xhrw7LNwxRVw772pD1qSWqJnz9QnfeGF8OCDsNZa8NhjWaeSWsUiWlLLffklVFfD3nunrbpffjldLBRC1skkFZsQ0itZ//43LLRQuo7ij39023AVDYtoSS3zyitpRNVNN8GZZ8KTT8JKK2WdSlKxW3vtNPpu//3ThJ9NN4WxY7NOJc2RRbSkHxcjDB4MG2wAX3wBjzwCp5wCXbpknUxSqVhggfR75oYb0mzptdeGm2/OOpX0oyyiJc3eV1/BfvvBQQfBr34FI0fCZptlnUpSqdpjj/R7ZrXVYPfd4eCD0+8hqQBZREuatddfTyOorrsOTj89Xfyz9NJZp5JU6nr1Su1iJ5+cVqd793amtAqSRbSkH7r2WlhvvTR26qGH4LTTbN+Q1HG6dUv90UOHwuefw/rr88L++1NVWUlZWRlVVVXkcrmsU6qTs4iW9J2vv06tG/vtl4roESPSLFdJykLfvvDyy7y32mqsP2QIp40bR48YaWhooKamxkJambKIlpS8+Wa6eHDw4PQy6sMPw7LLZp1KUmfXsycbf/oppwP7A08BFUBjYyO1tbWZRlPnZhEtKV0R37s3vP8+3H9/ehm1a9esU0kSAA3jx3MGsAOwMjAc2AIYN25cprnUuVlES53Z1Klw+OGw115px7CRI6F//6xTSdL3VFRUAHAPsB7wIfAQ8KdFFkljOKUMWERLnUgul6OqqoqysjLWXWEFPlpjDbj8cjjuuLTl7vLLZx1Rkn6grq6O8vJyAN4C+gB3dunCyZ99lkbhTZ6caT51ThbRUieRy+WoqamhoaGBNWPk9gkTWHD0aJ4+4gg477x0NbwkFaDq6mrq6+uprKwkhMASlZV8PWQInHsu3HYb9OkDb72VdUx1MiEWwcsgvXv3jsOGDcs6hlTUqqqqaGhoYBdgCPAJsCPwcWUlY91iV1KxGjo0bdIyfXqaa7/99lknUpELIQyPMfae0/NciZY6ifENDZwO3AKMJPUVvoQX5kgqclttBcOHw4orwg47pM2hmpqyTqVOwCJa6gwmT+be+efnNOAq0lXtH+YfmnHBjiQVraoqeOYZ2GcfOOMMGDAAPvss61QqcRbRUql7913YaCP6ffMNx3XrxoHA1PxD5eXl1NXVZZlOktrG/PPDNdfAxRfDAw/A+uvD6NFZp1IJs4iWStnjj6edB8ePp+yBB1j76qv/d2FOZWUl9fX1VFdXZ51SktpGCHDkkfDoo/Dpp+mCwyeeyDqVSpQXFkql6oor4KijYKWV4K67YOWVs04kSR3nnXdgu+1gzBi48krYd9+sE6lIeGGh1FlNm5Y2UDnsMOjXD55/3gJaUufz05/Cc8/BxhvDfvvBqae6MYvalEW0VEomTUpXql9+OZxwQlqBXmSRrFNJUjYWWyz1R++/P5x1FlRXwzffZJ1KJaJr1gEktZExY9KW3ePGpVmp9jpLEnTvDoMHp1fkTj45/Y684w5Ycsmsk6nIuRItlYIXX4QNN4SPP04X1FhAS9J3QoCTToIbb4Rhw9IFh2++mXUqFTmLaKnY3XsvbLYZlJfDs8/CRhtlnUiSCtPuu8Njj8EXX6SFh8cfzzqRiphFtFTM/vGPtKnAaqulC2hWXTXrRJJU2DbcMF1w/ZOfpIuvr70260QqUhbRUjGKMV1pXlOTLiR8/PH0H4Ikac5WXDG9crfJJmlyxymnOLlDrWYRLRWbadO+u9L8wAPTBI4FF8w6lSQVl0UXhfvvT79H//Qn2GsvJ3eoVZzOIRWTL7+EXXaBhx6C009Pq9EhZJ1KkopTt26pLW7lleHEE+GDD+DOOx0NqhZxJVoqFh98kF56fOSRNK7ptNMsoCVpXoWQ5upffz088wxsvjl89FHWqVQELKKlYjBqVBrJ9NZbcM89cMABWSeSpNKy555w993wxhtpl8OxY7NOpAJnES0Vuqeegl/9CqZMgSeeSBuqSJLaXv/+6dW+SZPS793XXss6kQpYuxXRIYSrQggfhRBea3bs/BDCGyGEV0II/wohLNpeX18qCbfemqZvLLVUGmG37rpZJ5Kk0rbhhvDkk2laxyabpN+90iy050r0EGDmJbOhwOoxxjWB0cBJ7fj1peI2eDDstlsqnJ95Bnr1yjqRJHUOq6+efu8usQRsuSWPHn88VVVVlJWVUVVVRS6XyzqhCkC7FdExxieBT2Y69lCM8dv8u88Dy7fX15eK2t//DgcdlDYCGDo0/SKXJHWcXr3g6af5pGdPfn3++WzY0ECMkYaGBmpqaiyklWlP9AHA/Rl+fanwxJjmlf7+97DzzmnUUnl51qkkqXNaemk2bWriWSAHHJ4/3NjYSG1tbYbBVAgyKaJDCLXAt6Sfydk9pyaEMCyEMGzixIkdF07KSoxpTukpp8A++8BNN0GPHlmnkqRO7T8TJtAfuBu4FDg1f3zcuHHZhVJB6PAiOoQwENgeqI5x9ntsxhjrY4y9Y4y9e/bs2WH5pEw0NcERR8B558Ghh8KQIdDVvZAkKWsVFRV8A/yWdLHXGcBFQOUKK2QZSwWgQ4voEEJ/4HjgNzHGxo782lLB+vZbGDgQLr8cjjsOLrsMypw+KUmFoK6ujvLycqaT+lD/AhwFPLbccjB1arbhlKn2HHF3A/AcsGoIYUII4UDgEmAhYGgIYWQI4Yr2+vpSIcvlclRVVTFfCNy/yCLwz3/CWWfBuee6C6EkFZDq6mrq6+uprKyEELikooIRu+9O1XPPwY47wtdfZx1RGWm314tjjHvO4vDg9vp6UrHI5XLU1NQQGxu5A+jf2Mjx3bqxVq9eVFtAS1LBqa6uprq6+vsH+/aFQw6BAQPgjju8CLwTCj/SllwwevfuHYcNG5Z1DKlNVFVV8UlDA/cAGwMHA1cBlZWVjHWbWUkqHkOGwAEHwOabpy3DLaRLQghheIyx95yeZ+Ol1MEmNzTwCLAhsCepgAav9JakojNwIFx7LTz+OGy3HUyenHUidSCLaKkj/fe/PN2tG2sAOwM3N3uooqIio1CSpLm2995w3XVpq/Btt4Uvv8w6kTqIRbTUUcaPh1//mpXKyvhtjx7c0+yh8vJy6urqMosmSZoHe+4JN9wAzz4L/fvDF19knUgdwCJa6ggTJqSeuY8+ouujj7LX4MFUVlYSQqCyspL6+vofXrQiSSoeu+2WNsl64QXYemv4/POsE6mdeWGh1N7efx822wz++1946CHo0yfrRJKk9nLHHamgXnvt9Dt/0UWzTqRW8sJCqRB88EFagf7gA3jgAQtoSSp1O+4It90GL78MW24Jn3ySdSK1E4toqb18+CFssQW89x7cfz9stFHWiSRJHWGHHeBf/4JXX03zpD/+OOtEagcW0VJ7+OijVECPGwf33Qcbb5x1IklSR9p2W7jzThg1Kv1/MHFi1onUxiyipbY2cWJaeXj3XbjnHthkk6wTSZKy0L9/2oRl9OhUSH/0UdaJ1IYsoqW29PHHqQfu7bfTL87NN886kSQpS1ttBffeC++8k/5P+PDDrBOpjVhES23lk09SAf3mm+klvL59s04kSSoEW2yRro0ZOzYV1V5sWBIsoqW28Omn6Rfj66+n8Ub9+mWdSJJUSDbdNC2wvPmmG7KUCItoaV599lkqml99FW6/Pf1ylCRpZltuCbfeCiNGwPbbQ2Nj1ok0DyyipXnx+edpZ6qXX05zQbfbLutEkqRCtsMOcN118PTTsNNOMGVK1ok0lyyipbn15ZewzTbw0ktwyy3pF6MkSXOy++5w5ZVpR8Pdd4dp07JOpLlgES3NjcbGtOr8wgtw000wYEDWiSRJxeSAA+Cii1Kf9MCBMH161onUSl2zDiAVnWnTYLfd0ktx118PO++cdSJJUjE66ij46is46SQoL4f6eggh61RqIYtoqTWamtKKwb33wuWXwx57ZJ1IklTMTjwRJk+GujpYcEH4298spIuERbTUUjHC0Uen1eezz4ZDD806kSSpFJx1ViqkL7wQFloIzjwz60RqAYtoqaVOOw0uvRSOPTatHEiS1BZCgAsuSK0dZ50FCywAJ5yQdSrNgUW01BIXXph+sR1wAJx/vi+1SZLaVghwxRWpkD7xxNTaccQRWafSj7CIlubk2mvhD39IFxAOGmQBLUlqH126wDXXpAlQRx6ZVqQHDsw6lWbDEXfSj7nzzrT63Ldv6oXu6t+dkqR21K1bGp3arx9NBxzAET17UlZWRlVVFblcLut0asYiWpqdxx5LQ/DXXRf+9S/o0SPrRJKkzqBHD27cfXeeC4ELJk1i8xhpaGigpqbGQrqAhBhj1hnmqHfv3nHYsGFZx1BnMmwYbL45VFTAk0/CEktknUiS1IlUVVXxWUMDTwGVwCbAy0BlZSVjx47NNFupCyEMjzH2ntPzXImWZjZqFPTvD0sumbZktYCWJHWwcePG8TmwDfAZcD9QlT+uwmARLTXX0AD9+qXe56FDYbnlsk4kSeqEKioqAHgP6A/0AB4E1vb/pYLRoiI6hLB8CGHz/P0eIYQF2jeWlIGPPoKttoIvv4QHH4SVVso6kSSpk6qrq6O8vByAUcAOwArAQz16pDF4ytwci+gQwgHAXcCV+UOVwJ3tGUrqcF98kVo4JkxIW3qvtVbWiSRJnVh1dTX19fVUVlYSQuC9ykpe+P3vWfLdd2G33WDatKwjdnotWYk+GugDfAEQYxwNLNWeoaSOkMvlqKqqonsIPLn00jS9/DLcdhv86ldZR5MkierqasaOHUtTUxNjx45l0wsugMsug/vug0MOgSIYDlHKWjL09psY49SQ32AihNAFcLcJFbVcLkdNTQ2NjY0MBjb55hsO696djT/5hOqsw0mSNDuHHALvvw9nngnLLgt/+lPWiTqtlqxEPxNCOB6YL98XfRNwT/vGktpXbW0tjY2NnAIcAJwBXDF1KrW1tRknkyRpDk4/HQ46COrq0sq0MjHHOdH5lecaoB9pBfpBYFCMsan94yXOiVZbKysrY58YuQa4BhiYPx5CoKmpw360JUmaO99+CzvvDPfcA7femu6rTbTJnOh8AX11jPHyGONOMcYd8/etMlTU9uzZkyuBh4GDmx2fMVJIkqSC1rUr3Hgj9OkDe+2VNgZTh/rRIjrGOB1YMYTQrYPySO3vlVcY8sUXvBkCvwVmXN9cXl5OXV1dlskkSWq58nK4+27o1Qt+8xt47bWsE3UqLemJfgd4KoRwUgjh6Bm39g4mtYsJE2Dbbem2xBKMvvBCFsuPDqqsrKS+vp7qai8rlCQVkSWWgAceSAV1//4wfnzWiTqNlkznGJe/ledvUnH6/HPYdts0E/qpp9h5rbXY+Wj/HpQkFbnKylRI//rXqZB++mlYbLGsU5W8ORbRMcZTOiKI1K6mTYNddoFRo9xMRZJUetZcE+68E/r1S//fPfAAdLMbtz3NsYgOIQwFfjDCI8bYr10SSW0tRqipgYcfhquuSr9gJEkqNZttBv/4BwwcCEccAYMGQXBrj/bSknaOPza7Px/wW2BK+8SR2sGZZ8KQIXDaabD//lmnkSSp/ey3H4weDWefDauuCscem3WiktWSdo5/z3ToiRDCzMekwnT11Wko/cCBqYiWJKnUnXVWKqSPOw5WWgkGDMg6UUma43SOEMLCzW6LhhD6Anarq/A99FBq49h4qDF0AAAgAElEQVRqK6iv9yUtSVLnUFYG11wDvXunGdIvvZR1opLUknaO/5B6ogPwLfAu39+fQio8r7ySLqz4+c/TTk5eXCFJ6kzKy+Guu2D99WGHHeCFF2C55bJOVVJaMid6xRhjRYxxhRhjrxjjFsAzc/qgEMJVIYSPQgivNTu2eAhhaAjhrfxbV7TV9j78ELbfHhZeOE3iWHjhrBNJktTxfvKTtC34F1+kQvqrr7JOVFJaUkTPqv/5hRZ83BCg/0zHTgQeiTGuDDySf19qO1OmwE47waRJadTP8stnnUiSpOysuWbaHvzll2HvvaGpKetEJWO2RXQIYakQwlrA/CGENUIIa+ZvG9OCTVdijE8Cn8x0eABwTf7+NcCOc5lb+qEZo+yeey71gq27btaJJEnK3nbbwQUXwB13wImuX7aVH+uJ3g44AFgeuKzZ8S+Bud2AZekY4wf5+/8Flp7LzyP90F//Ctdem6Zx7Lpr1mkkSSocRx0Fb74J558Pq6wCBx2UdaKiN9siOsZ4NXB1CGG3GOPNbf2FY4wxhPCDTVxmCCHUADUAFRUVbf3lVWruvReOPz4Vz6e4yaYkSd8TAvz97/DOO3DYYdCrF/Ttm3WqohZinG0d+92TQtga+AVpsxUAYoxnt+DjqoB7Yoyr599/E9gsxvhBCGEZ4PEY46pz+jy9e/eOw4YNm2NOdVL/+Q9suCGsvDI89VS6IlmSJP3Q55/DRhvB+++n9sef/SzrRAUnhDA8xth7Ts9ryZzoy4D9gGOA+YG9gZXmMtdd+c9F/u2dc/l5pGTSpHTF8QILpAsJLaAlSZq9RRZJEzu6dUuTrCZNyjpR0WrJdI6NY4x7AR/HGE8BNqAFRXQI4QbgOWDVEMKEEMKBwJ+BrUIIbwFb5t+X5s7UqWkW9Pvvp4slnMQhSdKc9eqVFp4mTICdd06TrdRqLdls5ZsZb0MIPwE+Bpad0wfFGPeczUM24GjexQhHHglPPAHXXQcbbJB1IkmSiseGG8KQIbDnnnDooXDVVe7s20otKaLvCyEsCvwFGAlM57sxdVI2LrkE/vEPOOkkqK7OOo0kScVnjz1g1Cg488w0FvbII7NOVFR+9MLCEEIZsF6M8d/59+cH5o8xzjz/uV15YaG+56GHYJttUi/07bdDWUu6kiRJ0g80NcGOO8J998HDD8Nmm2WdKHNtcmFhjLEJGNTs/a87uoCWvufNN2G33WD11VMbhwW0JElzr6ws/X+68sppTOy4cVknKhotqUAeCyEMaPck0px8+mlafe7eHe66CxZcMOtEkiQVv4UXThfoT52aVqUbG7NOVBRaUkQPBP4VQvg6hPBJCOHTEIKr0epY336bVqDHjk0tHJWVWSeSJKl0rLoqXH89jBwJNTXpAn79qJYU0UsC3YAFgZ7593u2ZyjpB445JvVqDRoEG2+cdRpJkkrPdtvBWWdBLsfwvfemqqqKsrIyqqqqyOVyWacrOHMsomOM04FdgRPy95cB1m7vYNL/DBkCF1+cCun99886jSRJpevkkxm33nqsff31rNTQQIyRhoYGampqLKRn0pIdCy8BNgf2yR9qBK5oz1DS/7z0UppfufnmcO65WaeRJKm0hUD///6X14GbgF75w42NjdTW1mYYrPC0pJ1joxjjIeQ3XclP5+jerqkkgI8/TjspLbUU3HQTdG3JWHNJkjQv3pgwgR1JReIdQHn++Dgnd3xPS4roafl50REghLAE0NSuqaTp09MuSh98ALfdBj1tw5ckqSNUVFQwBtgD+AVwdbPj+k5LiuhLgduAniGEM4CnAV9XV/s69VQYOhQuvRTWWy/rNJIkdRp1dXWUl5fzEHAisBtwSrdu1NXVZZyssMzx9fEY47UhhOHAlvlDu8YYX2vfWOrU7rgDzj4bDjoo3SRJUoeprq4GoLa2lr82NPDr8nLO+PprwuKLZ5yssPzott//e1IIawIbk1o6nokxvtLewZpz2+9O5M0308rzz34GTz4J882XdSJJkjq3xkbYaKO0V8OLL6bdDUtYm2z7nf9EtcANwLLA8sD1IYST5j2ilORyOaqqqlg4BEavsUa6gvW22yygJUkqBOXl6VXirl1hwAD48susExWElvRE7wusF2P8Y4yxFliftIuhNM9yuRw1NTU0NDQwGPjptGnsPHUquSefzDqaJEmaoaoqTcoaPRr23ReanDHRkiL6A77fO901f0yaZ7W1tTQ2NnIsaUefE4H7p0xxFqUkSYWmb184//y0Kv2Xv2SdJnNz7IkOIdwOrAc8SOqJ7ge8CIwDiDEe084Z7YkuYWVlZWwWI0OB20lXAAOEEGjyr1xJkgpLjLD77nD77fDII7DpplknanMt7Yluye4V9+ZvMzw/16mkmWyw7LLc9N57vAkc0Oy4syglSSpAIcCVV8LIkbDHHjBiBPzkJ1mnykRLRtwN7ogg6oS++Ya7e/SgO7ATMDl/uLy83FmUkiQVqoUXTgMANtgA9toLHnqoU+4q3JLpHP1DCC+GED4KIXwSQvg0hPBJR4RTiTv6aJYcM4YRv/89UyorCSFQWVlJfX39/2ZUSpKkArTGGnD55fDYY3DaaVmnyURLeqLfJrWqvkqz7b5jjNPbN9p37IkuQVdeCQcfDCefDK46S5JUnA4+OP2ffs89sN12WadpE202JxqYAIyMMU6LMU6fcZv3iOq0hg2DI46Afv3gzDOzTiNJkubWRRfB2mvDPvukzVg6kZY0sBwP3B1CeByYMuNgjPGi9gqlEvbZZ7DbbukihOuvhy5dsk4kSZLm1vzzwy23wLrrwq67wtNPQ48eWafqEC1ZiT4DmA4sCvRsdpNaJ0Y48EAYPz4NbF9iiawTSZKkebXSSjBkSHql+dhjs07TYVqyEr1CjHH1dk+i0nfJJWmu5F/+An36ZJ1GkiS1lZ12SgX0X/8Kv/oV7Lln1onaXUtWoh8MIWzR7klU2mb8dbrDDnBMu+/PI0mSOto556QC+uCDYdSorNO0u5YU0QcAD4cQJjviTnNlRh/0Msukl3tCyDqRJElqa926pXbN8nLYZRf46qusE7WrlhTRSwLdgEVIvdBLYk+0WipGOOig1Ad9442w+OJZJ5IkSe1lueXS4IBRo+CQQ1IdUKLmWETnx9ntCpyQv78MsHZ7B1OJuPTStKvROefAhhtmnUaSJLW3LbeEM86AXA7q67NO025asmPhJcDmwD75Q43AFe0ZSiVi+PDUB7399vZBS5LUmdTWwtZbw9FHp3qgBLWknWOjGOMhwDcAMcZPgO7tmkrF7/PPUx/00kunPuiylvyoSZKkklBWBtddl+qAXXaBTz/NOlGba0llMy2EUAZEgBDCEjTb/lv6gRl90A0NqQ/aedCSJHU+Sy6ZNmJ57720T0SJ9UfPtogOIcyYIX0pcBvQM4RwBvA0cG4HZFOxuuwyuPXW1Ae90UZZp5EkSVnZYAM4+2z4179g0KCs07SpEGfzV0EI4aUY4zr5+78AtgQC8HCM8bWOiwi9e/eOw4YN68gvqbn10kvpAsItt4S777aNQ5Kkzq6pCbbdFp54Al58EVYv7D38QgjDY4y95/i8HymiR8QYf9nmyeaCRXSR+PxzWHddmDIFRoxIL+NIkiR9+CGstVZq8XzxxTRLukC1tIj+sW2/e4YQZjtSIcb4t7lKptIUY9qhaOzY9JemBbQkSZph6aXh2mvTxI5jjoErin/Q24+91t4FWBBYaDY36TuXX54uHqirS1t+SpIkNdevHxx/fOqNvu22rNPMsxb1RGfNdo4CN2IE9OkDffvCPffYBy1JkmZt6lTYeGN46y0YORIqK7NO9AMtbef4sWontGEelarJk9M86J4908s0FtCSJGl2undP42+nT4fqavj226wTzbUfq3j6dlgKFa+jj4Z33klbe9oHLUmS5mTFFVNLxzPPwJlnZp1mrs22iM7vTCjN3i23wNVXw0knwaabZp1GkiQViz33hIED4U9/gscfzzrNXJltT3QhsSe6AI0bl0bVrLIKPP00dOuWdSJJklRMJk9Oo3EnT4aXXy6YV7TboidamrXp02GffVIfUy5nAS1JklpvwQVTf/SkSXDAAUW3LbhFtFrvvPPgySfh4othpZWyTiNJkorVL3+Z6oq774ZLLsk6TatYRKt1XnwRTj01TeTYb7+s00iSpGJ39NGw3Xbwf/+X2jqKRCZFdAjhDyGE/4QQXgsh3BBCmC+LHGqlyZNhr71gmWXSTkPBKYiSJGkehZAGFSyxBJ9vsw2rVVRQVlZGVVUVuVwu63Sz1eFFdAhhOeBooHeMcXXSzoh7dHQOzYWjj4YxY+C662CxxbJOI0mSSkXPnjy8//4s9MEHHDt+PDFGGhoaqKmpKdhCOqt2jq7A/CGErkA58H5GOdRSzcfZbbJJ1mkkSVKJOSiX42zgIGC3/LHGxkZqa2szTDV7mYy4CyH8DqgDvgYeijFWz+I5NUANQEVFxboNDQ0dG1LfmTHObtVV4amnnMYhSZLaXFlZGV1i5AmgElgRmAqEEGhqauqwHAU74i6EsBgwAOgFLAssEELYe+bnxRjrY4y9Y4y9e/bs2dExNYPj7CRJUgeoqKjgW2BPYAtSAT3jeCHKop1jS+DdGOPEGOM04HZgowxyqCXOPTeNs7vkEvjpT7NOI0mSSlRdXR3l5eWMA0bnj5WXl1NXV5dlrNnKoogeB/QJIZSHEALQFxiVQQ7NyQsvwGmnwe67w777Zp1GkiSVsOrqaurr66msrCSEQGVlJfX19VRX/6DrtyBk1RN9BrA78C0wAjgoxjhlds932+8MfPllGoA+bVqa2bjoolknkiRJanct7Ynu2hFhZhZjPA04LYuvrRY6+mh49114/HELaEmSpJm4Y6F+6OabYcgQOPlk+PWvs04jSZJUcCyiBUAul6OqqoqKEPh8zz2Z9NOfpu29JUmS9AMW0SKXy1FTU8O4hgYGA12amtjsvffI3Xxz1tEkSZIKkkW0qK2tpbGxkcOArYBjgf98803B7hAkSZKUNYtoMW7cOFYCzgfuB+qbHZckSdIPWUSLqhVW4BpgCmm/+hkKdYcgSZKkrFlEi9v69GEj4Ajg/fyxQt4hSJIkKWsW0Z3dK6/wy3/9i4b11+fZioqi2CFIkiQpa5lstqICMXVq2s57scWovOcexvbsmXUiSZKkomAR3ZmdeWba0vvOO8ECWpIkqcVs5+is/v1vOOccGDgQfvObrNNIkiQVFYvozqixMbVxLL88XHhh1mkkSZKKju0cndFJJ8Ho0fDII7DIIlmnkSRJKjquRHc2jz4KF10ERx0FW2yRdRpJkqSiZBHdmXz+Oey/P6y8Mvz5z1mnkSRJKlq2c3Qmf/gDTJgAzzwD5eVZp5EkSSparkR3FnffDVdfDSecAH36ZJ1GkiSpqFlEdwaTJsHBB8Oaa8Jpp2WdRpIkqejZzlHqYoTDD4dPPoGHHoIePbJOJEmSVPQsokvdjTfCLbfA2WenlWhJkiTNM9s5StmHH8KRR8IGG8Bxx2WdRpIkqWRYRJeyo46CyZPTBYVdfdFBkiSprVhZlarbb09tHHV1sNpqWaeRJEkqKa5El6JPPkkXE669tm0ckiRJ7cCV6FJ0zDHw8cfwwAPQrVvWaSRJkkqOK9Gl5v774Zpr0qYqa6+ddRpJkqSSZBFdSr74Ag45JPVAn3JK1mkkSZJKlu0cpeSEE2DCBHj2WTdVkSRJakeuRJeKxx+HK66AP/wB+vTJOo0kSVJJs4guBY2NcNBB8NOfwllnZZ1GkiSp5NnOUQpOOQXeeQceewzKy7NOI0mSVPJciS52zz8PF1wAhx4Km22WdRpJkqROwSK6mE2ZAgccAMsvD+eem3UaSZKkTsN2jmJ21lkwalSaDb3wwlmnkSRJ6jRciS5WI0bAn/8M++0H/ftnnUaSJKlTsYguRtOmpTaOnj3hb3/LOo0kSVKnYztHMTr/fBg5Em6/HRZfPOs0kiRJnY4r0cXm9dfhjDNgt91gp52yTiNJktQpWUQXk+nT4cADYaGF4OKLs04jSZLUadnOUUwuvzzNhb7uOlhqqazTSJIkdVquRBeL996Dk0+Gfv1gr72yTiNJktSpWUQXi6OPTlM5Lr8cQsg6jSRJUqdmO0cxuOuuNInjnHNgxRWzTiNJktTpuRJd6L78Eo44AtZYA449Nus0kiRJIqMiOoSwaAjh1hDCGyGEUSGEDbPIURROPTX1Qw8aBN26ZZ1GkiRJZNfO8XfggRjjLiGE7kB5RjkK2/DhcNFFcOihsKF/Z0iSJBWKDi+iQwiLAJsAAwFijFOBqR2do+B9+y0cfHAaZXfOOVmnkSRJUjNZtHP0AiYCV4cQRoQQrgwhLJBBjoKUy+WoqqrimG7dYMQIntp1V1hkkaxjSZIkqZksiuiuwDrA5THGXwJfASfO/KQQQk0IYVgIYdjEiRM7OmMmcrkcNTU1NDU0cCZwD9D/yivJ5XJZR5MkSVIzIcbYsV8whJ8Az8cYq/Lv/xo4Mca43ew+pnfv3nHYsGEdlDA7VVVVNDQ0cBewBfBzYBxQWVnJ2LFjM80mSZLUGYQQhscYe8/peR2+Eh1j/C8wPoSwav5QX+D1js5RiMaNG8dvgR2AU0kF9IzjkiRJKhxZTec4CsjlJ3OMAfbPKEdB+cXyy3PR+PG8RBpfMkNFRUVWkSRJkjQLmRTRMcaRwByXyTub21ZdlaXHj+c3wPT8sfLycurq6rKMJUmSpJm4Y2GheO45VnnkEd7q149JlZWEEKisrKS+vp7q6uqs00mSJKmZrNo51Ny0aXDIIbDccvzs1lsZu9BCWSeSJEnSj7CILgR/+xu8+irccQdYQEuSJBU82zmyNmYMnHEG7LgjDBiQdRpJkiS1gEV0lmKEww+Hrl3h4ouzTiNJkqQWsp0jSzfdBA8+CBddBMsvn3UaSZIktZAr0Vn54gv4wx9g3XXTarQkSZKKhivRWTntNPjwQ7jrLujSJes0kiRJagVXorPwyiupB7qmBtZbL+s0kiRJaiWL6I4WIxxxBCy6KLgToSRJUlGynaOj/fOf8PTTcOWVsMQSWaeRJEnSXHAluiN99hkcdxz06QP77591GkmSJM0lV6I70imnwKRJcP/9UObfL5IkScXKSq6jjBgBl10Ghx0G66yTdRpJkiTNA4vojtDUlGZBL7EE/OlPWaeRJEnSPLKdoyMMGQLPP5/eLrpo1mkkSZI0j1yJbm+ffAInnAAbbwz77pt1GkmSJLUBi+j2VlsLn34Kl14KIWSdRpIkSW3AIro9vfgiDBoERx0Fa66ZdRpJkiS1EYvo9jJ9erqYcOml4fTTs04jSZKkNuSFhe1l8GAYNgxyOVhkkazTSJIkqQ25Et0eJk2Ck06CzTaDPffMOo0kSZLamEV0ezjpJPjiC7jkEi8mlCRJKkEW0W3t+efhyivh97+HX/wi6zSSJElqBxbRbWnGxYTLLgunnpp1GkmSJLUTLyxsS1dcASNGwE03wUILZZ1GkiRJ7cSV6LYycSL88Y/Qty/sumvWaSRJktSOLKLbyimnwJdfwsUXezGhJElSibOIbgsjR0J9PRx5JKy2WtZpJEmS1M4soudVjHD00bDEEu5MKEmS1El4YeG8uuUWeOopGDQIFl006zSSJEnqAK5Ez4vGRvi//4O114YDD8w6jSRJkjqIK9Hz4rzzYPx4uO466NIl6zSSJEnqIK5Ez61x4+Dcc2G33WCTTbJOI0mSpA5kET23jjsujbI7//ysk0iSJKmDWUTPjSeegJtvhhNOgIqKrNNIkiSpg1lEt9b06fC738EKK6TVaEmSJHU6XljYWldeCS+/DDfdBOXlWaeRJElSBlyJbo1PP4Xa2nQh4a67Zp1GkiRJGbGIbo0zzkiF9N//ni4qlCRJUqdkEd1Sr78Ol1wCBx+cNleRJElSp2UR3RIxwu9/DwsuCGedlXUaSZIkZcwLC1vi7rth6FC48ELo2TPrNJIkScqYK9FzMmUKHHMMrLYaHH541mkkSZJUADJbiQ4hdAGGAe/FGLfPKsccXXABvPMOPPggdOuWdRpJkiQVgCxXon8HjMrw689WLpejqqqK5ULgq5NPZvw660C/flnHkiRJUoHIpIgOISwPbAdcmcXX/zG5XI6amhoaGho4G+gaI9u+/jq5XC7raJIkSSoQWa1EXwgcDzRl9PVnq7a2lsbGRtYH9gMuAF775htqa2szTiZJkqRC0eFFdAhhe+CjGOPwOTyvJoQwLIQwbOLEiR2UDsaNGwfAgsCzQN1MxyVJkqQsVqJ/BfwmhDAWuBHYIoRw3cxPijHWxxh7xxh79+zAsXIVFRUAPJoPOnmm45IkSVKHF9ExxpNijMvHGKuAPYBHY4x7d3SO2amrq6O8vPx7x8rLy6mrq5vNR0iSJKmzcU70TKqrq6mvr6eyspIQApWVldTX11NdXZ11NEmSJBWIEGPMOsMc9e7dOw4bNizrGJIkSSpxIYThMcbec3qeK9GSJElSK1lES5IkSa1kES1JkiS1kkW0JEmS1EoW0ZIkSVIrWURLkiRJrWQRLUmSJLWSRbQkSZLUShbRkiRJUitZREuSJEmtVBTbfocQJgINWefIwJLApKxDFDnP4bzx/M0bz9+88fzNG8/fvPH8zZtiPn+VMcaec3pSURTRnVUIYVhL9m7X7HkO543nb954/uaN52/eeP7mjedv3nSG82c7hyRJktRKFtGSJElSK1lEF7b6rAOUAM/hvPH8zRvP37zx/M0bz9+88fzNm5I/f/ZES5IkSa3kSrQkSZLUShbRGQshLB5CGBpCeCv/drHZPG96CGFk/nZXs+O9Qgj/DiG8HUK4KYTQvePSZ68l5y+EsHYI4bkQwn9CCK+EEHZv9tiQEMK7zc7t2h37HWQjhNA/hPBm/ufmxFk83iP/8/R2/uerqtljJ+WPvxlC2LojcxeKFpy/Y0IIr+d/3h4JIVQ2e2yW/5Y7kxacv4EhhInNztNBzR7bL//v/a0Qwn4dm7wwtOD8XdDs3I0OIXzW7DF//kK4KoTwUQjhtdk8HkIIF+XP7yshhHWaPebP35zPX3X+vL0aQng2hLBWs8fG5o+PDCEM67jU7STG6C3DG3AecGL+/onAubN53uTZHL8Z2CN//wrgsKy/p0I7f8AqwMr5+8sCHwCL5t8fAuyS9ffRweesC/AOsCLQHXgZ+PlMzzkcuCJ/fw/gpvz9n+ef3wPolf88XbL+ngrw/G0OlOfvHzbj/OXfn+W/5c5ya+H5GwhcMouPXRwYk3+7WP7+Yll/T4V2/mZ6/lHAVc3e79Q/f/lzsAmwDvDabB7fFrgfCEAf4N/5453+56+F52+jGecF2GbG+cu/PxZYMuvvoa1urkRnbwBwTf7+NcCOLf3AEEIAtgBunZuPLxFzPH8xxtExxrfy998HPgLmOES9hK0PvB1jHBNjnArcSDqPzTU/r7cCffM/bwOAG2OMU2KM7wJv5z9fZzLH8xdjfCzG2Jh/93lg+Q7OWMha8vM3O1sDQ2OMn8QYPwWGAv3bKWehau352xO4oUOSFYkY45PAJz/ylAHAtTF5Hlg0hLAM/vwBcz5/McZn8+cHSvz3n0V09paOMX6Qv/9fYOnZPG++EMKwEMLzIYQZheISwGcxxm/z708AlmvHrIWopecPgBDC+qTVm3eaHa7Lv/R0QQihRzvlLCTLAeObvT+rn5v/PSf/8/U56eetJR9b6lp7Dg4krWrNMKt/y51JS8/fb/P/Lm8NIazQyo8tZS0+B/k2ol7Ao80Od/afv5aY3Tn256/1Zv79F4GHQgjDQwg1GWVqM12zDtAZhBAeBn4yi4dqm78TY4whhNmNS6mMMb4XQlgReDSE8CqpsCl5bXT+yK8k/BPYL8bYlD98Eqn47k4ax3MCcGZb5JZCCHsDvYFNmx3+wb/lGOM7s/4MndbdwA0xxikhhENIr4pskXGmYrQHcGuMcXqzY/78qUOEEDYnFdEbNzu8cf7nbylgaAjhjfzKdlGyiO4AMcYtZ/dYCOHDEMIyMcYP8kXeR7P5HO/l3475//buLsSKMgzg+P9JscKyNLsQK3BDsBDT8EJKTPoQ7UKCDKQvSUGE7iKI2C6si+ouKCIiU4ggwlTaECpMxYKi79SyTBIiL7JMilRM6Oli3hPTwV3Pod095+j/B8POec68O/M+vLP77Ow7ZyJiJzAH2ET1b6ax5WrhFcChYe9Ahw1H/iJiArAV6C//nmt878ZV7JMRsQF4eBgPvVsdAq6svT7duGls81NEjAUuAY602PZs11IOIuJWqj/0bsrMk434IOfyuVTEnDF/mXmk9nId1b0PjbYLm9ruHPYj7G7tnIPLgQfrAcdfSwbLseOvRRExi+rcXVI/n2vj73BEbKGantSzRbTTOTpvAGjc4bsCeLN5g4iY2JhmEBGTgRuBb7Kapb8DWDZU+7NcK/kbB2yhmuP2RtN7U8rXoJpPfdq7jc8ynwDTo/pkl3FUv2ib79Kv53UZsL2MtwFgeVSf3jENmA58PErH3S3OmL+ImAO8CCzNzMO1+GnP5VE78u7QSv6m1F4uBfaV9XeARSWPE4FFJXYuaeX8JSJmUN389mEt5vhrzQBwf/mUjnnA7+WCi+OvBRFxFbAZuC8z99fi4yPi4sY6Vf56+3dup+9sPNcXqnmm7wHfA9uASSU+F1hX1m8A9lDdhb0HWFVr30dVxBwANgLnd7pPXZi/e4FTwJe1ZXZ5b3vJ6V7gVeCiTvdplPJ2O7Cf6gpUf4k9QVX0AVxQxtOBMr76am37S7vvqK4ydLw/XZi/bcDPtfE2UOKDnsvn0tJC/p4Cvi552gHMqLVdWcblAeCBTvelG/NXXq8Fnm5q5/ir8vAa1ac0naKa17wKWAOsKe8H8HzJ7x5gbq2t4+/M+VsHHK39/Pu0xPvK2PuqnN/9ne7L/118YiQ78U4AAAImSURBVKEkSZLUJqdzSJIkSW2yiJYkSZLaZBEtSZIktckiWpIkSWqTRbQkSZLUJotoSeqw8nm0H0TEklrsroh4uwPHsjMi5o72fiWp1/jEQknqsMzMiFgDbIyIHVQ/m58EFo/kfmtPO5UktckiWpK6QGbujYi3gEeA8VRP2PzP45gjYjFVcT0G+DUzb4mIScB6qgcZHAdWZ+buIeJrgatL/MeIWAlsAK4DvgUuLPsaA7xM9eCiBNZn5jMjmQNJ6iUW0ZLUPR4HPgf+oipe/xURlwMvAQsy82ApkhttvsjMOyLiZuAVYPYQcYBrgfmZeSIiHgKOZ+Y1ETGr7J+y7dTMnFn2f+kI9VmSepJFtCR1icw8FhGvA39m5smmt+cBuzLzYNn2txKfD9xZYtsj4rKImDBEHKrHkJ8o6wuAZ8t2uyNid4n/APRFxHPAVuDd4e6vJPUybyyUpO7yd1lG0rEzbZCZR6mmeOwE1gDrRviYJKmnWERLUm/4CFgQEdMAatM53gfuKbGFVHOl/xgi3mwXcHfZbiYwq6xPBs7LzE3AY8D1I9IrSepRTueQpB6Qmb9ExGpgc0ScBxwGbgPWAuvLNIzjwIrSZLB4sxeADRGxD9gHfFbiU0u8cbHl0eHtkST1tsjMTh+DJEmS1FOcziFJkiS1ySJakiRJapNFtCRJktQmi2hJkiSpTRbRkiRJUpssoiVJkqQ2WURLkiRJbbKIliRJktr0D6+gpdmirKsiAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7fa65fb55910>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "\n",
    "if rank == 0:\n",
    "    # build matplot lib graph of result only on proc 0\n",
    "\n",
    "    # 1st build exact solution hiRes\n",
    "    big = np.linspace(y0,y1)\n",
    "    cool = analyticTemperature(big, h, k, c0, c1)\n",
    "\n",
    "    pylab.rcParams[ 'figure.figsize'] = 12, 6\n",
    "    pyplot.plot(coords, numerical, 'o', color = 'black', label='numerical') \n",
    "    pyplot.plot(big, cool, color = 'red', label=\"exact\") \n",
    "    pyplot.xlabel('Y coords')\n",
    "    pyplot.ylabel('Temperature')\n",
    "    pyplot.show()\n"
   ]
  },
  {
   "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.5.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}