{
 "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",
    "**Note:** depending on the choice of $k$ the other terms of the equation must be dimensionaly consistent.\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": {
    "collapsed": true
   },
   "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": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import underworld as uw\n",
    "import glucifer\n",
    "import numpy as np\n",
    "uw.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.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",
    "  * $ q \\cdot n = \\phi $ on $ \\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",
    "An example implementation  \n",
    "\n",
    "```nbc = uw.conditions.NeumannCondition( fn_flux=phi, variable=tField,\n",
    "                                      nodeIndexSet=mesh.specialSets[\"MinJ_VertexSet\"] )\n",
    "```\n",
    "\n",
    "Applies `phi` as a flow to the `tField` over the boundary vertices in the set `nodeIndexSet`.\n",
    "\n",
    "Here `phi` can be an `underworld.Function` or `underworld.MeshVariable` type.\n",
    "\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": 3,
   "metadata": {
    "collapsed": true
   },
   "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": 4,
   "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": 5,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "solver = uw.systems.Solver(ss)\n",
    "solver.solve()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAtoAAAF3CAYAAACbhOyeAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzs3Xl4VNXBx/HvCSA6qBUBqaKZKZXiLkrEfcN9w/1VGyta\nJWqta91qXKp909Zq9XVFI7i1cbdasS5VXKmCBgEF94WACxXFPSoI5/3jDooUJJBM7mTy/TzPfTJz\n703ml/tE++v13HNCjBFJkiRJLass7QCSJElSKbJoS5IkSQVg0ZYkSZIKwKItSZIkFYBFW5IkSSoA\ni7YkSZJUABZtSZIkqQAs2pIkSVIBWLQlSZKkArBoS5IkSQXQMe0ALaV79+4xl8ulHUOSJEklbuzY\nsR/EGHss6rySKdq5XI76+vq0Y0iSJKnEhRAamnKeQ0ckSZKkArBoS5IkSQVg0ZYkSZIKwKItSZIk\nFYBFW5IkSSoAi7YkSZJUABZtSZIkqQAs2pIkSVIBWLQlSZKkArBoN0NdXR25XI6ysjJyuRx1dXVp\nR5IkSVKRKJkl2FtbXV0dVVVVNDY2AtDQ0EBVVRUAlZWVaUaTJElSEQgxxrQztIiKiopYX1/fap+X\ny+VoaGhgIPAl8DbwHtArm2Xy5MmtlkOSJEmtK4QwNsZYsajzvKO9hKZMmQLAcCCX3zcHmNbQAAMG\nwKqrLnjr1Qs6d04ptSRJklqLRXsJlZeX09DQwO7AasCq+W2NLl04oGtXeOUVGDkSPv30v795pZUg\nl1vwls1CJvODn11XV0d1dTVTpkyhvLycmpoah6tIkiQVGYv2EqqpqaGqqopJjY1Myu/LZDLUXn01\nzFt6P/0U3nkH3n77u23KFGhogHHj4O67YebM7//wBRXxn/4UVl+dm0aNouroox0bLkmSVOQco90M\nLXJnec4cmDYNJk9e8NbQ8L0iPgt4C3h9vu3LVVbh0cmToVOn5v9ikiRJWqimjtG2aBe7OXPgvffg\njTfg9df54+GH81NgdaAPsNy853bokAw9WX116NMH1ljju61XLwghlV9BkiSplFi0S9Tc2U7m6kFS\nujft1o2/HH00vP56sr366vfHhy+7LPTt+/3yvcYaSSH34UxJkqQmc9aREjV3bPjcMdrTgS8yGY65\n5JLvjw2PMRmS8vLL39+efBLmXVinrAx6905K99prwzrrJNsaa8DSS//X5/sgpiRJUtNYtNuYuaV2\nkWU3BFh55WTbdtvvH/vii+SO97wF/MUX4cEHYdas5JyysuRu97rrflu+73nzTY4+5xw++/JLwAcx\nJUmSfohDR/SdWbPgtddg4sTvb6+/ntwhB74CXgJeyG/jgQ9XXZXnpk5NL7ckSVIrcoy2Wk5jI7z8\nMoP792dtYB1gXZL5w7+1yiqwwQbQr993W+/eyZ1xSZKkEuIYbbWcTAY23JDHs1lunOdBzK7A+sDA\nrl05a7vtYPx4eOABmD07OWHZZWH99b9fvtdd14cvJUlSu2DRVpPN/yDmR8AzmQxHXHbZdw9ifvVV\nMt57/PhkGzcObrwRrrgiOd6pE6y3HlRUwEYbJV/XWsv5vyVJUslx6IgWyxLNOjJnDrz1Fjz3HIwd\nC/X1yfbJJ8nxpZdO7nbPW7779k3mBW/O50qSJBWAY7RV3ObMSRbhmVu66+uTEv7FF8nxLl1gww1h\nwACemDWLX9bW8sZXX3377ZlMhtraWsu2JElqdRZttT2zZyfTDs4t3s8+m9wF//prAKYCo/Pb08CH\n5eW8Ms+YcUmSpNZg0VZpmDmTAZ07swl8u/WeewhYaqONYJNNYNNNk6+5nEvNS5Kkgmpq0XbuNRW3\npZbi/WyWy4BK4KdAT2BPYNjyyydDTIYPh5//PJlOcOWVYd994f/+L7kr/s03qcaXJEntl0VbRa+m\npoZMJvPt+/eBhzMZfnTllfDoo8lDlePHw9ChsOOOyUwnJ56YPFi5wgqw/fbwu9/Bww/D55+n9ntI\nkqT2xaEjahMWe9aRd96Bf/8bnnwSRo2CCROS1S07dEgW1tlii2TbfHP48Y9b7xeRJEltnmO0pXl9\n8gmMHv1d8R4zJpnzG5KpBLfZ5rttvuLt1IKSJGleFm3ph8ycmcxo8sQT8PjjSQH/7LPk2BprJIV7\n222584MPOOSUU75dpAecWlCSpPYu9aIdQrgW2B14P8a4Tn7fisCtQA6YDPxPjPGjBXzvbOCF/Nsp\nMcZBi/o8i7aa5ZtvkrHdjz2WjPt+8slvx3O/CDwGPAo8DkwHstkskydPTimsJElKUzHMOnI9sPN8\n+04HRsYY+wAj8+8X5MsYY7/8tsiSLTVbx47Jw5OnnAL33QcffQRjxnAa0AD8Arid5EHMF4ATGxrg\n3nu/uwsuSZI0n4IV7RjjE8CM+XbvCdyQf30DsFehPl9qlo4dYcAAbs1m2RVYEdiY5P8ZvgscGQLs\nsQesuCJsuSWcey489RTMmpVqbEmSVDxae3q/njHG9/Kvp5FMibwgS4cQ6kMIo0MIlnGlZu7Ugt8A\nzwDnA3tnMtw1fDiMHAknn5w8VHnuuckMJt26waBBcOml8NJLyUwnkiSpXeqY1gfHGGMIYWEtJBtj\nfCeE0Bt4JITwQozxjflPCiFUAVUA5eXlBUyr9mruA4/zzzpy0NwHIQcOhD/+EWbMgEceSebqfvhh\nGDEiOb7KKsk83jvtlMzx3b17Sr+JJElqbQWddSSEkAPunedhyFeAbWKM74UQVgYeizH2XcTPuD7/\nM+74ofN8GFJF5a23vivdDz+cFPEQknHgu+wCO++cvO7QIe2kkiRpMRXDw5ALcg8wOP96MPCP+U8I\nIXQNIXTOv+4ObE4y8YPUdvzkJzBkCNx6K7z/fjKH9znnQFkZnHcebLoprLQSHHQQ3HADTJuWdmJJ\nktTCCla0Qwg3A08DfUMIb4cQDgf+BOwQQngN2D7/nhBCRQhhWP5b1wTqQwgTSGZU+1OM0aKttqtD\nB9h446RoP/00TJ8ON9+cPEz56KNw6KGw8sqw4YZwxhnJ1IKzZlFXV0cul6OsrIxcLkddXV3av4kk\nSVoMLlgjpWnOHHj+ebj/fnjggWTZ+Nmz+TqT4Z6vv+bu2bO5H/gIF8qRJKlYpL5gTWuzaKskfPIJ\njBzJrYccwtZffMGPgdnAv4F7gbGrrMLIt99OxntLkqRUWLSlNqysrAxipIJkedXdgQ3nHuzdG3bf\nPdm23hqWWiq1nJIktUfF+jCkpCYoLy8nAs8C5wD9gVWBM1ZcEdZcE2prv5sucL/94Prr4YMP0ows\nSZLmY9GWitDchXLm9VEmw9qXXpos/f7hh8lc3QcdlDxgedhh0LMnbLstXHIJNDSklFySJM1l0ZaK\nUGVlJbW1tWSzWUIIZLPZ7z8ImckkQ0euvhrefhvq65MZS6ZPhxNOgFwumcXkvPOShy1LZIiYJElt\niWO0pVLz2mtw993J9vTTScnu3Rv22ivZNtvMhXIkSWoGH4aUlCyEc889SekeORJmzoQePWDQINhn\nn2R5eB+mlCRpsfgwpCT48Y+hqgruuy8ZVnLLLbDddnDbbbDbbsnqlIMHJ+O+v/76229zsRxJkprP\nO9pSe/T11/Dww3D77fCPf8DHH8Pyy8OgQTzeowd7X3UVH3355benu1iOJEnfceiIpKaZORMeeSQp\n3XffDTNm8BkwArgDuB/4Cshms0yePDnNpJIkFQWHjkhqmqWWgp13huHDYdo0dgJuBnYA/g5MB24B\nNmpogHnuckuSpB9m0Zb0nU6deCWb5UhgZWA74K/ANsDtkIzp/sUv4J//TO6ES5KkhbJoS/qeuYvl\nzAYeAX4FrL7MMjz8298mC+T885/JHN4rr5w8aPnIIzB7dsqpJUkqPhZtSd+zoMVyrrrmGrb/wx+S\npd+nTUtmKdl1V7j55mQWk1VXheOPh9GjXRxHkqQ8H4aUtOQaG5OpA2++ObnT/fXXyaqUBxyQ3P1e\nbz0IIe2UkiS1KB+GlFR4mQzstx/ceSe8/z7ccAOsuSZceCH06wfrrgt/+hNMnZp2UkmSWp1FW1LL\nWH55OOSQ5A73tGkwdCissAL89reQzcK228K118Inn6SdVJKkVmHRltTyuneHo46CUaPgjTfg3HPh\nnXfg8MOT1SoPOABGjIBZs9JOKklSwVi0JRVW795w1lnwyiswZgwccUQyU8mgQcnMJb/+dbI/Rpd+\nlySVFB+GlNT6Zs2CBx+Ev/0tWQL+q6/4tGdPLpkxg2GzZjElf5pLv0uSipEPQ0oqXp06JXNx33JL\nMp772muZ+PHHnDVrFg3Aw0AlEBsbqa6uTjmsJElLxqItKV0/+hEcdhhbzJxJFjgLyAF/A6YBZzY0\nwL//7fzckqQ2x6ItqSiUl5czBfhfoA+wFXAncFAIsMUW0Lcv/OEPThUoSWozLNqSisLcpd8BIvAk\n8OtMhnuvuQauuw5WWQWqq5OpAnfcMVkk58svU80sSdIPsWhLKgoLWvq9traWAw4/HA49FB57DF5/\n/bsZTH7+82TWkmOOgXHj0o4vSdJ/cdYRSW3PnDlJ8b72WrjjjmTp9/79k6kDf/7zZPEcSZIKxFlH\nJJWusjIYODCZHvC99+DSS2HmTDj66OQu92GHwVNP+QClJClVFm1JbVvXrnDssTBhQrLwTWVlcpd7\n881hnXXg4ovhgw/STilJaocs2pJKQwgwYADU1iZ3uYcNg+WWg5NOgl694MADYeTIZNiJJEmtwKIt\nqfQsuywcfjiMHg3PPw9HHQX/+hdsvz387GdwwQXe5ZYkFZxFW1JpW3dduOQSePddqKtL7m6femry\n9eCD4d//pu5vfyOXy1FWVkYul6Ouri7t1JKkEuCsI5Lan0mT4Kqr4MYb4dNPmRgCV8bI34DPgEwm\nQ21tLZWVlWknlSQVoabOOmLRltR+ffEFp5WX8z8zZtAf+ByoA4YCH2ezTJ48OdV4kqTi5PR+krQo\nXbpwwUcfUQFsBNwK/AIYD9zS0AA33ODqk5KkJWbRltSulZeXA1APHAH0Ao4DunfsmKxI2asXnHIK\nvPVWeiElSW2SRVtSu1ZTU0Mmk/n2/cfA8EyGMdddl6w+uf32yVzcP/0p7LknPPSQC+FIkprEoi2p\nXausrKS2tpZsNksIgWw2mzwIefDBsPXWcNttMHkyVFcn0wXuuCOstRZcfjl89lna8SVJRcyHISWp\nqb7+Gm6/HS67DJ55JlkQ59BD4ZhjoG/ftNNJklqJD0NKUkvr3DmZe3vMmGTbay+4+mpYYw3YaScY\nMQJmz047pSSpSFi0JWlJDBiQzMM9dSr87/8mc3MPGgR9+sBFF8Enn6SdUJKUMou2JDXHSisl47ff\neisZVrLqqvCb3yRfjz8e3ngj7YSSpJRYtCWpJXTqBPvtB088AWPHwt57w9ChyR3uvfeGxx93thJJ\namcs2pLU0jbcMBlWMnkynHEGPPkkbLMN9O+f7J85M+2EkqRWYNGWpEJZZZVk/PbUqVBbm8xaMngw\nZLPJ/unTqaurI5fLUVZWRi6Xo66uLu3UkqQWYtGWpEJbZhkYMgQmToQHH4R+/eCss/imVy9mDh5M\nl4YGYow0NDRQVVVl2ZakElGwoh1CuDaE8H4IYeI8+1YMITwUQngt/7XrQr53cP6c10IIgwuVUZJa\nVQjJgjf33w8vvsitnTtzwOzZTAL+CWwDNDY2Ul1dnW5OSVKLKOQd7euBnefbdzowMsbYBxiZf/89\nIYQVgXOAjYEBwDkLK+SS1GatuSa/+OILVgPOBPoDjwLPAps2NMA336QaT5LUfAUr2jHGJ4AZ8+3e\nE7gh//oGYK8FfOtOwEMxxhkxxo+Ah/jvwi5JbV55eTkzgBogCwwBlgVuBlh9dbjkEvj88xQTSpKa\no7XHaPeMMb6Xfz0N6LmAc3oBU+d5/3Z+nySVlJqaGjKZDABfA8OAimWW4bETT4TVVoMTToDy8mSe\n7mnTUs0qSVp8qT0MGWOMQLMmlQ0hVIUQ6kMI9dOnT2+hZJLUOiorK6mtrSWbzRJCIJvNcvU117DN\nRRclUwI+/TQMHAh//GMyU8kRR8BLL6UdW5LURK1dtP8TQlgZIP/1/QWc8w6w2jzvV83v+y8xxtoY\nY0WMsaJHjx4tHlaSCq2yspLJkyczZ84cJk+eTGVl5XcHN9kE7rgDXn01Kdk33QRrrQV77AGjRqUX\nWpLUJK1dtO8B5s4iMhj4xwLOeRDYMYTQNf8Q5I75fZLUPq2+OlxxBUyZAueeC6NHw5ZbJtt997ni\npCQVqUJO73cz8DTQN4TwdgjhcOBPwA4hhNeA7fPvCSFUhBCGAcQYZwC/J3n4/lngvPw+SWrfuneH\ns8+Ghga49NKkeO+2WzIv9803O1OJJBWZEEvkTkhFRUWsr69PO4YktZ5Zs5KC/ac/JWO3e/eGU06B\nQw+FpZdOO50klawQwtgYY8WiznNlSElqqzp1gkMOSVacvOuu5I730UdDLgfnnw+ffpp2Qklq1yza\nktTWlZXBXnslY7cfeQTWWw9OP/27qQHfX9Bz55KkQrNoS1KpCAG23Rb+9S+or4cddvhuasBjj4Wp\nUxf9MyRJLcaiLUmlqH9/uP12ePllqKyEq66Cn/4UqqrgzTepq6sjl8tRVlZGLpejrq4u7cSSVHIs\n2pJUyn72Mxg2DN54A4YMgRtvZE6fPjB4MJ0bGogx0tDQQFVVlWVbklqYRVuS2oPy8mQu7jff5Nou\nXdh79mxeAm4C1gYaGxuprq5OOaQklRaLtiS1J6usQtXnn5MDzgd2ByYCdwLdGhrSTCZJJceiLUnt\nTHl5OdOBM4AscC4wEBgLsPvuyewlkqRms2hLUjtTU1NDJpMB4CPgd8CayyzD+P33T0r2ppsmM5aM\nGpVmTElq8yzaktTOVFZWUltbSzabJYRANpvlwmuuod9tt8HkyXDhhfD887DllknhfuqptCNLUpvk\nEuySpP/W2AhDh8Kf/5wseLPjjvC73yV3uyWpnXMJdknSkstk4De/gTffhAsugHHjYLPNYOedHcMt\nSU1k0ZYkLVyXLnDyyfDWW3D++TB2bHJXe5dd4Jln0k4nSUXNoi1JWrQuXeDUU5PC/ac/wbPPwsYb\nw267Ja8lSf/Foi1Jarpll4XTTksK9x/+kAwjGTAgmRZw7Ni000lSUbFoS5IW33LLwW9/m8xSUlMD\nTz8NFRWw774waVLa6SSpKFi0JUlLbrnl4Iwzkjvcv/sdPPQQrLsuHHwwvP562ukkKVUWbUlS8y2/\nPJxzTlK4TzkF/v53WGMNGDIEpk5NO50kpcKiLUlqOd26JbOTvPkm/OpXcOONsPrqcPzx8J//AFBX\nV0cul6OsrIxcLkddXV3KoSWpMFywRpJUOFOmwHnnwfXXQ+fOTBo4kJ1GjuSdL7/89pRMJkNtbS2V\nlZXp5ZSkxeCCNZKk9JWXw7Bh8NJLsNderHnvvUz68kvOApbLn9LY2Eh1dXWaKSWpICzakqTC69MH\n6uroB4wEzgPeBE4AOgNTpkxJM50kFYRFW5LUaj7NZtkX2Ah4DrgYeAU4acUVYfbsVLNJUkuzaEuS\nWk1NTQ2ZTIZ6YCdgO+CDsjIu/PBDWH99uOceKJFnhyTJoi1JajWVlZXU1taSzWYJIfBGNsvLN9wA\nt90GM2fCnnvCllvCqFFpR5WkZnPWEUlScZg1C667Lln45r33kmXd//CHZAEcSSoizjoiSWpbOnWC\nqqpkRck//hGefDIZTnLIIclS75LUxli0JUnFJZOB009PFr055RS4/Xbo2xdOOAGmT087nSQ1mUVb\nklScVlwxWWXytddg8GC4/PJklck//hHmWfBGkoqVRVuSVNxWXRVqa2HiRNhmGzjjDPjZz5LVJp0S\nUFIRs2hLktqGNdaAf/wDHn8cVl4ZDjsM+veHf/0r7WSStEAWbUlS27LVVjB6NNx8M3z6Key0U7JN\nmJB2Mkn6Hou2JKntKSuDAw+El16Ciy6CZ5+FDTaAQw+Ft99OO50kAU0s2iGEVUMI2+Zfdw4hdCls\nLEmSmqBzZzjxRHjjDfjNb5K73H36JOO4P/kk7XSS2rlFFu0Qwi+Be4Bh+V1Z4B+FDCVJ0mLp2hUu\nuABeeQX22SeZmWT11ZOZSmbNSjudpHaqKXe0jwM2AT4FiDG+CqxUyFCSJC2RXA7q6qC+PllR8thj\nYb314J//pO5vfyOXy1FWVkYul6Ouri7ttJJKXMcmnPNVjHFmCAGAEEIHIBQ0lSRJzdG/P4wcCSNG\nwMknw+678+OyMpadM4cINDQ0UFVVBUBlZWW6WSWVrKbc0f53COFUYOn8OO1bgXsLG0uSpGYKAQYN\ngokTOa9rVzacM4cJwFCgB9DY2Eh1dXXKISWVsqYU7VOBz4CXgeOBkYD/ZpIktQ1LLcXvPv6Y1YEr\ngCOA14BTgGkNDalGk1TafrBo54eJXBdjHBpj3DvGuFf+9ZxWyidJUrOVl5czg+Ru0TrAE8CfgVc7\ndoQ77oAYU80nqTT9YNGOMc4GeocQOrVSHkmSWlxNTQ2ZTAaAV4BBwO6dO7PcyivD/vsni+DU16ea\nUVLpacrQkTeAJ0MIvw0hHDd3K3QwSZJaSmVlJbW1tWSzWUIIZLNZDho+nK5vvQVXXw2vvgobbQSH\nHALvvJN2XEklIsRF/OeyEMLvF7Q/xnhWQRItoYqKiljv3QhJ0pL49NNk7u2LL4aOHZMFb046CZZe\nOu1kkopQCGFsjLFikectqmi3FRZtSVKzvflmssLk3XdD797wl7/AnnsmM5hIUl5Ti3ZTVoZ8KITw\nr/m3lokpSVIR6d0b7roLHnoIllkG9t4bdtwRJk1KO5mkNqgpY7TPBM7KbzUk0/xNKGQoSZJStf32\nMH48XHpp8pDk+uvDccfBRx+lnUxSG7LIoh1jHDPP9niM8Thgq+Z8aAjh+BDCxBDCpBDCCQs4vk0I\n4ZMQwvj8dnZzPk+SpMXWsWOyhPtrr8GQIXDFFdCnD1x1FcyenXY6SW1AU4aOLD/PtkIIYTug65J+\nYAhhHWAIMABYH9g9hLD6Ak59MsbYL7+dt6SfJ0lSs3TvDkOHwnPPwdprw9FHJ0u8P/542skkFbmm\nDB2ZBEzMfx1HsirkkGZ85prAmBhjY4zxG+BxYJ9m/DxJkgpv/fXhscfgttuSISTbbAP/8z8wZUra\nySQVqaYU7d4xxvIY42oxxp/EGAcC/27GZ04EtgwhdAshZIBdgdUWcN6mIYQJIYT7QwhrN+PzJElq\nGSEkC9y89BL87ndw772wxhpQUwNff512OklFpilFe8wC9j2zpB8YY3wJOB/4F/AAMB6Yf7Dbc0A2\nxrg+cBlw94J+VgihKoRQH0Konz59+pJGkiRp8WQycM458PLLsOuucOaZsM46cP/9aSeTVEQWWrRD\nCCuFENYHlgkhrBtCWC+/bQFkmvOhMcbhMcb+McatgI+AV+c7/mmM8fP86/uATiGE7gv4ObUxxooY\nY0WPHj2aE0mSpMVXXg533AEPPggdOiSle6+94K230k4mqQj80B3t3YDLgVWBK4Er8tsZJFP9LbEQ\nwkr5r+Uk47Nvmu/4j0NIVgcIIQzI5/ywOZ8pSVLB7LgjPP88/OlP8PDDsNZacO658OWXaSeTlKKF\nFu0Y43Uxxi2Bw2OMW86z7RpjvL2Zn3tnCOFFYARwTIzx4xDCUSGEo/LH9wMmhhAmAJcCB8ZSWcJS\nklSalloKTjstGU4yaFAyhnvttWHECOrq6sjlcpSVlZHL5airq0s7raRW0KQl2EMIOwFrA0vP3Rdj\n/EMBcy02l2CXJBWVkSOTebhfeon7y8r49Zw5vJk/lMlkqK2tpbKyMtWIkpZMSy7BfiUwGDgJWAY4\nGFjQvNeSJGmu7baD8eOpWWEFtpgzh0nAuST/Q9rY2Eh1dXXKASUVWlNmHdkixvhz4MMY41nAxli0\nJUlatKWW4qxPPqEv8HfgbOBFYA9givNvSyWvKUX7q7lfQwg/zr9fpXCRJEkqHeXl5bwHVAJbA18A\n9wAPLr00NDSkmk1SYTWlaN8XQlgBuJBkzuvJwG2FDCVJUqmoqakhk0lmxX0C6Aec2akT286Zk8xO\ncv75MGtWqhklFcYPFu0QQhlwf4zx4/xMIz8B1o0xntEq6SRJauMqKyupra0lm80SQqBXNsua111H\nx1degR12gNNPh3794Ikn0o4qqYUtctaREML4GGO/VsqzxJx1RJLUJo0YkcxO0tAAgwfDBReAi7BJ\nRa3FZh0BHg0h7NkCmSRJ0vz22AMmTUrubNfVQd++UFsLc+aknUxSMzWlaB8K3BVC+DKEMCOE8FEI\nYUaBc0mS1H506QJ//CNMmADrrQdHHgmbbw7jx6edTFIzNKVodwc6AcsCPfLv/W9akiS1tLXWgkcf\nhRtugDfegP794cQT4bPP0k4maQkssmjHGGcD+wOn5V+vTPLQtCRJamkhwCGHwCuvQFUVXHIJrLkm\n3H132skkLaamrAx5ObAt8Iv8rkbgqkKGkiSp3evaFYYOhaefhm7dYO+9k+3tt9NOJqmJmjJ0ZLMY\n45HkF66JMc4AlipoKkmSlNh4Y6ivT+bbfvDBZHjJZZfB7NlpJ5O0CE0p2rPy82lHgBBCN8BHoSVJ\nai2dOsGpp8LEibDZZnDcccnXCRPSTibpBzSlaF8B3An0CCGcC4wCzi9oKkmS9N9694b774ebboLJ\nk5OHJU89Fb74Iu1kkhagKQ9D3gicSbIE+wxg/xjjLYUOJkmSFiAEOOggeOklOOywZIGbddaBBx5I\nO5mk+TTljjZAB2AWMHMxvkeSJBXKiivCNdckS7cvvTTssktSwP/zn7STScpryqwj1cDNwCrAqsBN\nIYTfFjqYJElqgi23TBa2Oe88+PvfYY01kgLuypJS6ppyd/oQYKMY45kxxmpgAMlqkZIkqRh07gxn\nnQUvvAAbbJDMvz1wIPf85S/kcjnKysrI5XLU1dWlnVRqV5pStN8DOs7zvmN+nyRJKiY/+xmMHAnD\nhjHz2WfZ4eSTOaChgQ4x0tDQQFVVlWVbakVNKdozgEkhhGEhhGuAF4APQggXhRAuKmw8SZK0WEKA\nww9n864HE6jFAAAePUlEQVRduY9kmrAxwAZAY2Mj1dXV6eaT2pGOiz6Ff+a3uUYXKIskSWohY999\nl/2AvUnm6X2GZPqw8xoaUs0ltSeLLNoxxuGtEUSSJLWc8vJyGhoauAt4FLgAOB04oGNHeOwx2Gab\nNONJ7UJTZh3ZOYTwbAjh/RDCjBDCRyGEGa0RTpIkLZmamhoymQwAHwNDgF07d6b7iivCttsmD0x+\n/HGqGaVS15Qx2pcDRwK9gB5A9/xXSZJUpCorK6mtrSWbzRJCIJvNUjl8OMu99RacfDIMHw5rrQV3\n3512VKlkhRjjD58QwmPAwBhjUU/IWVFREevr69OOIUlS21BfD4cfDs8/D/vtB5dfDj17pp1KahNC\nCGNjjBWLOq8pd7RPBUaEEE4JIRw3d2t+REmSlJqKiqRs19TAiBHJ3e26OljEDThJTdeUon0uMBtY\ngWTIyNxNkiS1ZZ06wRlnwLhxyRzcBx8Me+4J776bdjKpJDRler/VYozrFDyJJElKx5prwqhRcMkl\nUF0Na68NF18Mgwcn83JLWiJNuaP9YAhhYMGTSJKk9HToACedlIzZXmcdOOww2G03mDo17WRSm9WU\nov1L4OEQwudO7ydJUonr0wcefxwuvTT5uvbacM01jt2WlkBTinZ3oBPwI5zeT5Kk0ldWBsceCy+8\nkDw0WVUFO+4IkyennUxqUxZZtGOMs4H9gdPyr1cG+hU6mCRJSlnv3vDwwzB0KIweDeuuC1deCXOK\nesZfqWg0ZWXIy4FtgV/kdzUCVxUylCRJKhJlZXDUUTBxImy6KRxzDGy3HbzxRtrJpKLXlKEjm8UY\njwS+AogxzgCWKmgqSZJUXLJZePDBZLz2c8/BeuvBFVd4d1v6AU0p2rNCCGVABAghdAP8p0qSpPYm\nBDjiCJg0CbbYAn79a9hhB2hoSDuZVJQWWrRDCHPn2L4CuBPoEUI4FxgFnN8K2SRJUjFadVV44AG4\n+mp45plk7PawYc5MIs3nh+5oPwMQY7wROBO4EPgI2D/GeEsrZJMkScUqhGQ2khdegP79YciQZN7t\nd95JO5lUNH6oaH+7FFSMcVKM8ZIY4//FGCe2Qi5JktQW5HIwcmQy7/Zjj8E66/DU0UeTy2YpKysj\nl8tRV1eXdkopFSEu5D/zhBDeBi5a2DfGGBd6LA0VFRWxvr4+7RiSJLVfr73G9N12o8drr3E3cCTw\nPpDJZKitraWysjLlgFLLCCGMjTFWLOq8H7qj3QFYFlhuIZskSdJ3+vRh46+/5mRgZ2ASsB/Q2NhI\ndXV1utmkFHT8gWPvxRjPa7UkkiSpzZs8dSp/Ae4DrgduB24BjnVmErVDTRqjLUmS1BTl5eUAvARs\nBlQD+wAvlpXBiBEpJpNa3w8V7e1aLYUkSSoJNTU1ZDIZAGYDfwC2WnppOq66KgwalMxO8tlnqWaU\nWstCi3Z+BUhJkqQmq6yspLa2lmw2SwiBbDbLscOG0fXVV+G002D4cFh/fRg1Ku2oUsEtdNaRtsZZ\nRyRJagNGjYJDDoHJk+HUU+Hcc6Fz57RTSYulJWYdKZgQwvEhhIkhhEkhhBMWcDyEEC4NIbweQng+\nhLBhGjklSVIL22ILmDABDj8czj8fBgxIFr2RSlCrF+0QwjrAEGAAsD6wewhh9flO2wXok9+qgKGt\nGlKSJBXOcsvBNdfAPffAtGlQUQEXXgizZ6edTGpRadzRXhMYE2NsjDF+AzxO8kDyvPYEboyJ0cAK\nIYSVWzuoJEkqoD32gIkTk6XbTzkFBg5MhpRIJSKNoj0R2DKE0C2EkAF2BVab75xewNR53r+d3ydJ\nkkpJjx5w551w/fUwbhystx5cdx2UyDNkat9avWjHGF8Czgf+BTwAjCeZAWixhRCqQgj1IYT66dOn\nt2BKSZLUakKAwYPh+edhww3hl7+EvfeG999PO5nULKk8DBljHB5j7B9j3Ar4CHh1vlPe4ft3uVfN\n75v/59TGGCtijBU9evQoXGBJklR4uRw88kgyXvv++2HddeG++9JOJS2xtGYdWSn/tZxkfPZN851y\nD3BIfvaRTYBPYozvtXJMSZLU2srK4De/gbFjoWfPZPz2McdAY2PayaTFlkrRBu4MIbwIjACOiTF+\nHEI4KoRwVP74fcCbwOvANcCvUsopSZLSsM468MwzcOKJcOWVycwk48alnUpaLC5YI0mSittDDyVj\nuD/4AGpqkjveZWndK5SKfMEaSZKkJtthh2RRmz32SFaT3H57mDp10d8npcyiLUmSil+3bnDHHTB8\neDKkZL314Lbb0k4l/SCLtiRJahtCSKb+GzcOfvYzOOCAZEjJp5+mnUxaIIu2JElqW/r0gVGj4Kyz\n4G9/g3794Kmn0k4l/ReLtiRJans6dYLzzoMnnkhWkdxySzj7bG668UZyuRxlZWXkcjnq6urSTqp2\nrGPaASRJkpbY5pvDhAlw7LHw+9/Tu6wM5swhAg0NDVRVVQFQWVmZbk61S97RliRJbdvyy8MNN3Bc\n9+6sOWcOE4AD8ocaGxuprq5OM53aMYu2JEkqCZd/+CH9gBeBW4BrgS7AlClTUs2l9suiLUmSSkJ5\neTmTga2A/wUGA88Bu/TsmWYstWMWbUmSVBJqamrIZDJ8A5wFbAd0CYF7pk+Hv/wF5sxJOaHaG4u2\nJEkqCZWVldTW1pLNZgkh8FY2y1NXXkmHPfaAk0+GXXaBadPSjql2JMQY087QIioqKmJ9fX3aMSRJ\nUrGJEa6+Gk48MXlw8vrrk9ItLaEQwtgYY8WizvOOtiRJKm0hwFFHQX099OwJu+6alO6vv047mUqc\nRVuSJLUPa68NY8bAr38N//d/sMkm8MoraadSCbNoS5Kk9mOZZeCyy+Cee2DqVOjfH268Me1UKlEW\nbUmS1P7ssUeyomT//jB4cLJ9/nnaqVRiLNqSJKl96tULRo6Es8+Gv/4VKiqS8i21EIu2JElqvzp2\nhHPPTQr3p5/CxhvD0KHJTCVSM1m0JUmStt0Wxo9Pvv7qV7D//vDxx2mnUhtn0ZYkSQJYaSX45z/h\nz3+Gf/wDNtggmaVEWkIWbUmSpLnKyuCUU+DJJ5PhI1tsARde6PLtWiIWbUmSpPltsgmMGweDBiXF\ne/fdYfr0tFOpjbFoS5IkLUjXrnDHHXDFFfDII9CvHzz2WNqp1IZYtCVJkhYmhOThyNGjYdllYbvt\n4H//16EkahKLtiRJ0qL06wdjx8JBB8FZZ8EuuziURItk0ZYkSWqKZZdNFraprYXHH09mJRk1Ku1U\nKmIWbUmSpKYKAYYMSYaSLLMMbLMNnH++Q0m0QBZtSZKkxTV3KMm++8Lpp/POhhvSb7XVKCsrI5fL\nUVdXl3ZCFYGOaQeQJElqk5ZfHm65hWeWWYb1b7iBe4ADgNENDVRVVQFQWVmZakSlyzvakiRJSyoE\n/uexx9gMmAU8AZwINDY2Ul1dnW42pc6iLUmS1AxTpkzhOaA/MAK4CLgL+LShIdVcSp9FW5IkqRnK\ny8sB+ATYFzge2BWY0KEDPPtsismUNou2JElSM9TU1JDJZL59fymwfefOdO3aFTbfHC6/HGJML6BS\nY9GWJElqhsrKSmpra8lms4QQyGazHDl8OMu+8grstBMceyz8/Ofw+edpR1UrC7FE/h9WRUVFrK+v\nTzuGJEnSd+bMSebZPvNM6NsX7rwT1lwz7VRqphDC2BhjxaLO8462JElSoZSVwW9/Cw89BB98ABtt\nBLfemnYqtRKLtiRJUqENHAjjxiUL3Rx4IBx3HMycmXYqFZhFW5IkqTX06gWPPgonngiXXQZbbw1v\nv512KhWQRVuSJKm1dOoEF10Et98OkybBBhskw0pUkizakiRJrW2//ZI5tnv2TGYm+f3vkwcnVVIs\n2pIkSWno2xfGjIHKSjj7bNh9d/jww7RTqQVZtCVJktLSpQvceCMMHQojR0L//q4mWUIs2pIkSWkK\nAY46CkaNSlaQ3GILuOaatFOpBVi0JUmSisFGG8Fzz8E220BVFQwZAl99lXYqNYNFW5IkqVh06wb3\n3QdnnAHDhsGWW8KUKWmn0hKyaEuSJBWTDh2gpgbuugteeSUZtz1yZNqptARSKdohhBNDCJNCCBND\nCDeHEJae7/ihIYTpIYTx+e2INHJKkiSlZq+9kgcjV1oJdtwR/vznZAy32oxWL9ohhF7AcUBFjHEd\noANw4AJOvTXG2C+/DWvVkJIkScVg7hSA++4Lp50G++8Pn32Wdio1UVpDRzoCy4QQOgIZ4N2UckiS\nJBW3ZZeFW2+FCy+Eu++GAQPg5ZfTTqUmaPWiHWN8B7gQmAK8B3wSY/zXAk7dN4TwfAjhjhDCaq0a\nUpIkqZiEAL/5TbJc+4cfJjOU/P3vaafSIqQxdKQrsCfwE2AVoEsI4eD5ThsB5GKM6wEPATcs5GdV\nhRDqQwj106dPL2RsSZKk9G27bTIF4FprJcNJTj+dm268kVwuR1lZGblcjrq6urRTKi/EVh5UH0LY\nH9g5xnh4/v0hwCYxxl8t5PwOwIwY449+6OdWVFTE+vr6Fs8rSZJUdL7+Go4/Hq6+mkfLyth/zhzm\nLt6eyWSora2lsrIy1YilLIQwNsZYsajz0hijPQXYJISQCSEEYDvgpXlPCCGsPM/bQfMflyRJatc6\nd4arruKUbt3YdM4c6oH184caGxuprq5OM53y0hijPQa4A3gOeCGfoTaEcF4IYVD+tOPy0/9NIJmh\n5NDWzilJklTs/jJjBluSzDLxFN9N4zbFRW6KQqsPHSkUh45IkqT2JpfL0dDQwErA7cBWJDNODC0v\n542GhnTDlbBiHjoiSZKkFlBTU0Mmk+F9YHvgcuBkYNTyyyezkyhVFm1JkqQ2qrKyktraWrLZLN+E\nwIXZLE8PGcLKr76aTAE4YULaEds1h45IkiSVmjFjYJ994OOP4dpr4YAD0k5UUhw6IkmS1F5tvDGM\nHQsbbAAHHpgs3z57dtqp2h2LtiRJUin68Y/hkUfgqKPgz3+GXXeFGTPSTtWuWLQlSZJK1VJLwdCh\ncM018NhjybjtF15IO1W7YdGWJEkqdUccAY8/Dl9+CZtuCnfckXaidsGiLUmS1B5sskkybnu99WD/\n/eHss2HOnLRTlTSLtiRJUnux8srw6KPwy1/C738P++4Ln32WdqqSZdGWJElqTzp3hmHD4JJLYMQI\n2GwzePPNtFOVJIu2JElSexMCHHccPPAAvPNO8pDko4+mnarkWLQlSZLaq+23h2eeSaYC3GEHuOIK\nKJHFDIuBRVuSJKk9W311ePrpZJ7tX/8ajjwSZs5MO1VJsGhLkiS1d8svD3ffDWeckcy5vd128P77\naadq8yzakiRJgrIyqKmBm29OpgGsqIBx49JO1aZZtCVJkvSdAw+EUaOSsdqbbw633ZZ2ojbLoi1J\nkqTv23BDqK9Pvh5wAJx5povbLAGLtiRJkv5bz54wcmSyfHtNDey3H3zxRdqp2hSLtiRJkhasc2eo\nrYWLL4Z//AO22AKmTk07VZth0ZYkSdLChQAnnAD33pusILnRRjB6dNqp2gSLtiRJkhZtl12S+ba7\ndIFttuHfxxxDLpejrKyMXC5HXV1d2gmLjkVbkiRJTbPWWvDMM/znJz9h8yuvZEhDA8RIQ0MDVVVV\nlu35WLQlSZLUdN26sUVjI7VANXAHkAEaGxuprq5ON1uRsWhLkiRpsbwxdSpHAicAewKjgFWBKVOm\npJqr2Fi0JUmStFjKy8sBuATYHegNPAsM6tkzxVTFx6ItSZKkxVJTU0MmkwHgAWBToDEE7vzwQ7jp\nplSzFROLtiRJkhZLZWUltbW1ZLNZQgg0ZrM8d+WVdNhsM6ishOpqV5IEQowx7QwtoqKiItbX16cd\nQ5Ikqf2aOROOOQaGDYN99oG//hXyd75LSQhhbIyxYlHneUdbkiRJLWOppZKVJC+6CO66C7beGt57\nL+1UqbFoS5IkqeWEACeeCHffDS+9BAMGwIQJaadKhUVbkiRJLW/QIBg1CmKELbZIlnBvZyzakiRJ\nKox+/eCZZ6BvX9hzT7jkkqR4txMWbUmSJBXOKqvA448nd7hPOCF5WPKbb9JO1Sos2pIkSSqsLl3g\nzjvh1FNh6FDYbTf45JO0UxWcRVuSJEmFV1YG55+fTP33yCOw2Wbw1ltppyooi7YkSZJaz+GHw4MP\nwrvvwsYbw9NPp52oYCzakiRJal0DB8Lo0bD88rDttnDzzWknKgiLtiRJklpf375J2R4wAH7+czjv\nvJKbkcSiLUmSpHR07w4PPQSHHALnnAOHHpos414iOqYdQJIkSe1Y585w/fXw058mZXvq1GSGkq5d\n007WbN7RliRJUrpCgLPPhr/+NVlNskRmJLFoS5IkqTgcfHAylOQ//4FNNklWlWzDLNqSJEkqHltv\nDU89lSxys802cNddaSdaYhZtSZIkFZc11khmJFlvPdh3X7j44jY5I4lFW5IkScVnpZXg0Udhn33g\npJPg2GPhm2/STrVYUinaIYQTQwiTQggTQwg3hxCWnu945xDCrSGE10MIY0IIuTRySpIkKUXLLAO3\n3QYnnwxXXAF77w2ff552qiZr9aIdQugFHAdUxBjXAToAB8532uHARzHG1YGLgfNbN6UkSZKKQlkZ\nXHABXHkl3HcfbLUVf7/8cnK5HGVlZeRyOerq6tJOuUBpDR3pCCwTQugIZIB35zu+J3BD/vUdwHYh\nhNCK+SRJklRMjj4aRoxg1osvstFxx7FcQwMxRhoaGqiqqirKst3qRTvG+A5wITAFeA/4JMb4r/lO\n6wVMzZ//DfAJ0K01c0qSJKnI7Lore624ImUx8m9gYH53Y2Mj1dXVaSZboDSGjnQluWP9E2AVoEsI\n4eAl/FlVIYT6EEL99OnTWzKmJEmSitD906axMfAC8NE8+6dMmZJSooVLY+jI9sBbMcbpMcZZwN+B\nzeY75x1gNYD88JIfAR/O/4NijLUxxooYY0WPHj0KHFuSJElpKy8v5x1gC2DcfPuLTRpFewqwSQgh\nkx93vR3w0nzn3AMMzr/eD3gkxjY4eaIkSZJaVE1NDZlM5nv7MpkMNTU1KSVauDTGaI8hecDxOZK7\n/mVAbQjhvBDCoPxpw4FuIYTXgZOA01s7pyRJkopPZWUltbW1ZLNZQghks1lqa2uprKxMO9p/CaVy\no7iioiLW19enHUOSJEklLoQwNsZYsajzXBlSkiRJKgCLtiRJklQAFm1JkiSpACzakiRJUgFYtCVJ\nkqQCsGhLkiRJBWDRliRJkgrAoi1JkiQVgEVbkiRJKgCLtiRJklQAJbMEewhhOtCQdo6UdAc+SDtE\nG+b1ax6vX/N4/ZrH69c8Xr/m8fo1X1u9htkYY49FnVQyRbs9CyHUxxgr0s7RVnn9msfr1zxev+bx\n+jWP1695vH7NV+rX0KEjkiRJUgFYtCVJkqQCsGiXhtq0A7RxXr/m8fo1j9evebx+zeP1ax6vX/OV\n9DV0jLYkSZJUAN7RliRJkgrAot1GhBBWDCE8FEJ4Lf+160LOmx1CGJ/f7pln/09CCGNCCK+HEG4N\nISzVeunT15TrF0LoF0J4OoQwKYTwfAjhgHmOXR9CeGuea9uvdX+DdIQQdg4hvJL/uzl9Acc75/+e\nXs//feXmOfbb/P5XQgg7tWbuYtGE63dSCOHF/N/byBBCdp5jC/xnuT1pwvU7NIQwfZ7rdMQ8xwbn\n/3l/LYQwuHWTF4cmXL+L57l2r4YQPp7nmH9/IVwbQng/hDBxIcdDCOHS/PV9PoSw4TzH/Ptb9PWr\nzF+3F0IIT4UQ1p/n2OT8/vEhhPrWS10AMUa3NrABfwZOz78+HTh/Ied9vpD9twEH5l9fBRyd9u9U\nbNcP+BnQJ/96FeA9YIX8++uB/dL+PVr5mnUA3gB6A0sBE4C15jvnV8BV+dcHArfmX6+VP78z8JP8\nz+mQ9u9UhNdvWyCTf3303OuXf7/Af5bby9bE63cocPkCvndF4M381675113T/p2K7frNd/6xwLXz\nvG/Xf3/5a7AVsCEwcSHHdwXuBwKwCTAmv7/d//018fptNve6ALvMvX7595OB7mn/Di2xeUe77dgT\nuCH/+gZgr6Z+YwghAAOBO5bk+0vEIq9fjPHVGONr+dfvAu8Di5yMvoQNAF6PMb4ZY5wJ3EJyHec1\n73W9A9gu//e2J3BLjPHrGONbwOv5n9eeLPL6xRgfjTE25t+OBlZt5YzFrCl/fwuzE/BQjHFGjPEj\n4CFg5wLlLFaLe/0OAm5ulWRtRIzxCWDGD5yyJ3BjTIwGVgghrIx/f8Cir1+M8an89YES/vefRbvt\n6BljfC//ehrQcyHnLR1CqA8hjA4hzC2T3YCPY4zf5N+/DfQqYNZi1NTrB0AIYQDJXaA35tldk//P\nXBeHEDoXKGcx6QVMnef9gv5uvj0n//f1CcnfW1O+t9Qt7jU4nOTu2FwL+me5PWnq9ds3/8/lHSGE\n1Rbze0tZk69BfsjST4BH5tnd3v/+mmJh19i/v8U3/7//IvD/7d1fiBVlGMfx7y/FCsvSLBL7gxuC\nhZiGF1KiZibahQQZSFmSgi10F4HIdmFeVHdBERFtChVI+I9WhDJbxYKisj9qaWYKkYSWSZGKST5d\nvO+padtdj3Rmz876+8Cwc56Z95yZh3d2n53zzswWSTslLW3SNjXE4GZvgP1D0lbg2m4WtRVfRERI\n6ul2MTdGxGFJLUCnpN2k4mfAa1D+yGckXgcWRcTZHF5OKtCHkG5FtAxY2YjtNpO0EJgMTC+E/3Ms\nR8R33b/DBWsTsCYiTkt6lPTtyswmb1MVLQDWRcSfhZj7n/UJSXeSCu2phfDU3P+uAd6VtC+fIa8c\nF9r9SETM6mmZpCOSRkXEj7kQPNrDexzOPw9K2g5MAtaTvtIanM86XgccbvgONFkj8idpGLAZaMtf\nBdbeu3Y2/LSk1cATDdz0/uowcH3hdXf9prbOD5IGA1cAx+psO9DVlQNJs0j/DE6PiNO1eA/H8oVU\n6JwzfxFxrPCynXQtRq3tjC5ttzd8C/u38zkGFwCPFQPuf3XpKcfuf3WSNIF07M4tHs+F/ndU0kbS\nUKhKFtoeOlIdHUDtyuVFwFtdV5A0vDakQdJI4A7g60hXFmwD5vfWfoCrJ39DgI2kMXfruiwblX+K\nNL6726uoB5hPgLFKd6wZQvpj3PXuA8W8zgc6c3/rABYo3ZVkDDAW+LiPtru/OGf+JE0CXgbmRcTR\nQrzbY7nPtrx/qCd/owov5wF78/w7wOycx+HA7By7kNRz/CJpHOmCvQ8LMfe/+nQAD+e7j0wBfs0n\nZdz/6iDpBmAD8FBE7C/Eh0q6vDZPyl91/+Y2+2pMT/VNpHGv7wHfAluBETk+GWjP87cDu0lXl+8G\nlhTat5AKnQPAWuDiZu9TP8zfQuAM8EVhmpiXdeac7gHeAC5r9j71Ud7uAfaTzmS15dhKUmEIcEnu\nTwdy/2optG3L7b4hna1o+v70w/xtBY4U+ltHjvd4LF9IUx35ewb4KudpGzCu0HZx7pcHgEeavS/9\nMX/59Qrg2S7t3P9SHtaQ7j51hjTOegnQCrTm5QJezPndDUwutHX/O3f+2oHjhd9/n+Z4S+57X+bj\nu63Z+/J/Jj8Z0szMzMysBB46YmZmZmZWAhfaZmZmZmYlcKFtZmZmZlYCF9pmZmZmZiVwoW1mZmZm\nVgIX2mZmFZHv1/uBpLmF2P2S3m7CtmyXNLmvP9fMrEr8ZEgzs4qIiJDUCqyVtI30O/xpYE6Zn1t4\nqqyZmZ0HF9pmZhUSEXskbQKWAUNJTzL916OxJc0hFeCDgJ8j4i5JI4BVpIdBnASWRsSuXuIrgJty\n/HtJi4HVwK3APuDS/FmDgFdJD38KYFVEPFdmDszMqsKFtplZ9TwFfAb8QSpw/ybpauAVYFpEHMqF\ndK3N5xFxr6SZwGvAxF7iALcAUyPilKTHgZMRcbOkCfnzyeuOjojx+fOvLGmfzcwqx4W2mVnFRMQJ\nSW8Cv0fE6S6LpwA7IuJQXveXHJ8K3JdjnZKukjSslzikR8KfyvPTgOfzersk7crxg0CLpBeAzcCW\nRu+vmVlV+WJIM7NqOpunMp041woRcZw0nGQ70Aq0l7xNZmaV4ULbzGxg+QiYJmkMQGHoyPvAgzk2\ngzR2+7de4l3tAB7I640HJuT5kcBFEbEeeBK4rZS9MjOrIA8dMTMbQCLiJ0lLgQ2SLgKOAncDK4BV\necjHSWBRbtJTvKuXgNWS9gJ7gZ05PjrHayduljd2j8zMqksR0extMDMzMzMbcDx0xMzMzMysBC60\nzczMzMxK4ELbzMzMzKwELrTNzMzMzErgQtvMzMzMrAQutM3MzMzMSuBC28zMzMysBC60zczMzMxK\n8BeIPkuO/gB22AAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f448e4d9f50>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Numerical flux at y =  -0.6 is -0.355222812219\n",
      "Exact     flux at y =  -0.6 is -0.298507462687\n",
      "\n",
      "Abs. error = 3.009e-04\n",
      "Rel. error = 6.742e-06\n"
     ]
    }
   ],
   "source": [
    "# create numpy arrays for analytics, parallel check first\n",
    "\n",
    "yvals       = np.zeros(mesh.specialSets['MinI_VertexSet'].count)\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": 7,
   "metadata": {
    "collapsed": true
   },
   "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": 8,
   "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": 9,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "solver = uw.systems.Solver(ss)\n",
    "solver.solve()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAtMAAAF3CAYAAABnkcdUAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzs3XecnVPix/HPmUwirk6CIHNvsqzoEoMgutV2lf2tEi7L\ntllWX1abtcoabS2i78WKMkpC9FWi9zJRI4ggE10QFkOKOb8/nhuGTZlkynNn5vN+ve5r7n3uvZPv\nPK878XVynnNCjBFJkiRJ864s7QCSJElSR2WZliRJkuaTZVqSJEmaT5ZpSZIkaT5ZpiVJkqT5ZJmW\nJEmS5pNlWpIkSZpPlmlJkiRpPlmmJUmSpPlkmZYkSZLmU3naAeZFr169Yi6XSzuGJEmSOrkxY8Z8\nHGPsPbfXdagyncvlqKurSzuGJEmSOrkQQn1zXuc0D0mSJGk+WaYlSZKk+WSZliRJkuaTZVqSJEma\nT5ZpSZIkaT5ZpiVJkqT5ZJmWJEmS5pNlWpIkSZpPlmlJkiRpPqVapkMIh4UQXg4hjA0hXBtC6Jlm\nnh+rra0ll8tRVlZGLpejtrY27UiSJEkqIamV6RDC8sDBQGWMcXWgGzA0rTw/VltbS1VVFfX19cQY\nqa+vp6qqykItSZKk76Q9zaMcWDCEUA5kgPdSzvOd6upqGhoa6AHsUjzW0NBAdXV1mrEkSZJUQlIr\n0zHGd4EzgUnA+8DnMcZ7fvy6EEJVCKEuhFA3efLkdss3adIkAKqAkcAlwAJNjkuSJElpTvNYAtgJ\n6AcsBywUQtjrx6+LMRZijJUxxsrevXu3W76KigoALgROBn4PPAIMXm65dssgSZKk0pbmNI+tgLdi\njJNjjNOBUcCGKeb5gZqaGjKZDI3AccDOwMrA/f/9LzzwQLrhJEmSVBLSLNOTgMEhhEwIIQBbAq+k\nmOcH8vk8hUKBbDZLCIHns1keOuMMei6/PPzsZ3DWWRBj2jElSZKUohBTLIQhhBOB3YEZwHPA72OM\nU2f3+srKylhXV9de8Wbtiy9g331h1CgYOhQuvRQWWijdTJIkSWpVIYQxMcbKub0u1dU8YozHxxgH\nxBhXjzHuPaciXTIWWQRuuAFOPRVGjIDBg2HChLRTSZIkKQVpL43XMYUARx8Nd94J770HlZVwxx1p\np5IkSVI7s0y3xNZbQ10d9OsHO+wAJ50EjY1pp5IkSVI7sUy3VL9+8NhjsNdecPzxsPPO8PnnaaeS\nJElSO7BMt4ZMBq64As47L5n6se668PLLaaeSJElSG7NMt5YQ4MADkzWo//tfWG89uO66tFNJkiSp\nDVmmW9uQIfDsszBwIOyxBxxyCEyblnYqSZIktQHLdFtYbrlkhPrQQ+Hcc2HzzZNVPyRJktSpWKbb\nSvfucPbZyVSPF16AQYPgoYfSTiVJkqRWZJlua7vvDk89BYsvDltuCf/8p9uQS5IkdRKW6faw2mrw\n9NOw005wxBGw227JtuSSJEnq0CzT7WXRRZNtyP/xD7jppmT5vHHj0k4lSZKkFrBMt6cQkpHpe++F\nKVOS5fOuvz7tVJIkSZpPluk0bLZZsnzeWmvB0KFw2GEwfXraqSRJkjSPLNNpWX75ZPm8gw+Gc86B\nLbZw+TxJkqQOxjKdph49YNgwuOaa7zd6uf/+tFNJkiSpmSzTpWCPPeCZZ2DJJeFnP4NTTqH2qqvI\n5XKUlZWRy+Wora1NO6UkSZJ+pDztACpaddWkUFdVQXU1S5WV8d/GRiJQX19PVVUVAPl8Pt2ckiRJ\n+o4j06Vk4YWhtpa/Lrkkmzc28ixQWXyqoaGB6urqNNNJkiTpRyzTpSYETpkyhSFAAB4D9i8+NWnS\npPRySZIk6X9YpktQRUUFdcAgYDRwIVALDFhhhVRzSZIk6Ycs0yWopqaGTCbDp8AOwLHA7sAT337r\nromSJEklxDJdgvL5PIVCgWw2CyFwTTbLA8ccw2IzZiTbkF9zTdoRJUmSBIQYY9oZmq2ysjLW1dWl\nHSM9770Hu+8Ojz4K++8PZ58NCyyQdipJkqROJ4QwJsZYObfXOTLdkSy3XLKpy1/+AhddBBttBG++\nmXYqSZKkLssy3dF07w5nnAG33AJvvAGDBsGoUWmnkiRJ6pIs0x3VjjvCc8/ByivDr34FhxwC06al\nnUqSJKlLsUx3ZLkcPPJIUqTPPReGDIGJE9NOJUmS1GVYpju6Hj3gnHPgxhth/HgYODCZAiJJkqQ2\nZ5nuLP7v/+DZZ+EnP4Gdd4bDD4fp09NOJUmS1KlZpjuT/v3hscfgwAPhrLNgk03ALcglSZLajGW6\ns1lgATjvPBgxAl5+GdZeG26/Pe1UkiRJnZJlurPadddk2kc2CzvsAEce6bQPSZKkVmaZ7sxWXBGe\neCLZLfEf/4DNNnPahyRJUiuyTHd2PXvChRfCtdfCSy8l0z5uvTXtVJIkSZ2CZbqrGDo0mfbRrx/s\ntBMceihMnZp2KkmSpA7NMt2VrLgiPP54ssnLsGGw0UbJluSSJEmaL5bprmaBBZJNXm66Cd58M9nk\n5brr0k4lSZLUIVmmu6qdd4bnn4c11oA99oA//hG+/jrtVJIkSR2KZborq6iABx+EY46BQgHWW4/b\nzjiDXC5HWVkZuVyO2tratFNKkiSVLMt0V9e9O5xyCtx1F99MmsQWRx3F5vX1xBipr6+nqqrKQi1J\nkjQblmklttmGjRdZhCeBy4GrgIWBhoYGqqur080mSZJUoizT+s6Y995ja+A4YA9gDLA2MMmNXiRJ\nkmbJMq3vVFRU0AicDGwOZIAngRMWXxxiTDWbJElSKbJM6zs1NTVkMhkAHgHWAkZ368bfpkyBHXaA\nyZNTzSdJklRqLNP6Tj6fp1AokM1mCSGwSDbL58OHw/nnw733wlprwX33pR1TkiSpZITYgf75vrKy\nMtbV1aUdo2t68cVkS/JXX4Wjj4YTT0xWApEkSeqEQghjYoyVc3udI9NqnjXXhGeegd//Hk49FTbZ\nBN56K+1UkiRJqbJMq/kWWijZ3OX66+GVV2DttZP7kiRJXZRlWvNut92SrchXWy2Z+vG738FXX6Wd\nSpIkqd1ZpjV/cjl4+GGorobLL4d11kkKtiRJUhdimdb8Ky+Hk09OVvj44gtYf30YNsw1qSVJUpdh\nmVbLbb45vPACbLMNHHoo/Pzn8OGHaaeSJElqc5ZptY5eveCWW5I1qR94IFn94z//STuVJElSm0q1\nTIcQFg8h3BBCeDWE8EoIYYM086iFQoADDoC6OlhmmWSE+uCD4Ztv0k4mSZLUJtIemR4G3BVjHECy\ne/UrKedRa1htNXj6aTjkEDjvPFh3XRg7Nu1UkiRJrS61Mh1CWAzYBLgMIMY4Lcb4WVp51Mp69oRz\nzoE774TJk6GyMinWXpwoSZI6kTRHpvsBk4HLQwjPhRAuDSEslGIetYVtt022It9yy2TKhxcnSpKk\nTiTNMl0ODAIuijEOBL4Cjv7xi0IIVSGEuhBC3eTJk9s7o1rD0kvD7bcnI9P33+/FiZIkqdNIs0y/\nA7wTY3yq+PgGknL9AzHGQoyxMsZY2bt373YNqFYUAhx4YHJx4tJLJyPUhxzixYmSJKlDS61Mxxg/\nAN4OIaxcPLQlMC6tPGonq68OzzyTTPk491xYbz146aW0U0mSJM2XtFfzOAioDSG8CKwNnJJyHrWH\nnj2TnRLvuCOZP11ZCWefTe1VV5HL5SgrKyOXy1FbW5t2UkmSpDkqT/MPjzE+D1SmmUEp2n77ZFT6\nD3+AP/+Z5crKmNHYSATq6+upqqoCIJ/Pp5tTkiRpNtIemVZXt/TScPPNHL3kkqzX2MiLwK7Fpxoa\nGqiurk4znSRJ0hxZppW+EDhjyhTWBl4HRgBXAIsCkyZNSjWaJEnSnFimVRIqKiqYAAwBTgDywAvA\nr5ZeOs1YkiRJc2SZVkmoqakhk8kwAziRpFR/GwLXf/QRVFfDtGkpJ5QkSfpflmmVhHw+T6FQIJvN\nEkLg/WyWukKBst/+Fk45BTbYAF59Ne2YkiRJPxBijGlnaLbKyspYV1eXdgy1t5tvht//Hhoa4Mwz\nYf/9k01gJEmS2kgIYUyMca6rzjkyrdK3887JEnqbbgoHHJAsqffee2mnkiRJskyrg+jTB/7zHzj/\nfHjoIVhjDRgxIu1UkiSpi7NMq+MIIRmZfv55WHFF2H132HNPmDIl7WSSJKmLskyr4/npT+Gxx+Ck\nk2DkyGSU+p570k4lSZK6IMu0OqbycjjuOHjySVh0UdhmGzjwQPjqq7STSZKkLsQyrY5tnXVgzBg4\n7DC44AIYOBCeeirtVJIkqYuwTKvjW3BBOOssuP9+mDoVNtwwGbWePj3tZJIkqZOzTKvz2HxzePFF\n2HtvOPlkGDwYxo1LO5UkSerELNPqXBZbDIYPh1GjYNIkGDQoGbVubEw7mSRJ6oQs0+qcfvlLGDsW\ntt4aDj88GbV+8820U0mSpE7GMq3Oa5ll4JZb4PLLk7Wp11wTLr4YYkw7mSRJ6iQs0+rcQoB9901G\nqTfcEPbfP1lG7+23004mSZI6Acu0uoa+feHuu+Gii+Dxx2H11ZO51Y5SS5KkFrBMq+sIAfbbL1nx\nY+214Te/gR135MbzzyeXy1FWVkYul6O2tjbtpJIkqYMoTzuA1O7694cHHoBzz2XGkUey2e23Mxio\nB+rr66mqqgIgn8+nGlOSJJU+R6bVNZWVwaGHsnXv3rwOXAeMAHoBDQ0NVFdXp5tPkiR1CJZpdWkP\nvv8+Q4CjgZ2AscWvkyZNSjWXJEnqGCzT6tIqKir4FjgdWAd4F7gZuDGTgU8+STWbJEkqfZZpdWk1\nNTVkMhkgGZVeHzi5e3d2+uYbWG01uPnmVPNJkqTSZplWl5bP5ykUCmSzWUIILJ/N0u/yyymrq4M+\nfZKdFPfcEz7+OO2okiSpBIXYgdbZraysjHV1dWnHUFcxfTqcdhr8/e+wxBLJGtX/939pp5IkSe0g\nhDAmxlg5t9c5Mi3NTvfucNxxUFcHyy8Pv/oVDB0KkyennUySJJUIy7Q0N2uuCU89lYxQjxqVzKW+\n4Ya0U0mSpBJgmZaao3t3+OtfYcyYZGvyXXeF3XZzlFqSpC7OMi3NizXWgCefhJqaZKWPVVeFkSPT\nTiVJklJimZbmVffucOyx8OyzkMslI9S77AIffJB2MkmS1M4s09L8Wn11eOIJOPVUuP32ZJT6yiuh\nA62QI0mSWsYyLbVEeTkcfTQ8/zyssgrssw9svz24HbkkSV2CZVpqDQMGwMMPw7nnwiOPJCt+XHQR\nNDamnUySJLUhy7TUWrp1g4MOgrFjYfBg+NOfYPPN4fXX004mSZLaiGVaam25HNxzD1x2GbzwQrJO\n9ZlnwowZaSeTJEmtzDIttYUQ4Le/hXHjYJtt4C9/gQ03TEatJUlSp2GZltrScsvBTTfBddfBxIkw\naBCceCJMm5Z2MkmS1Aos01JbCwF23z0Zpd51VzjhBFhnHXj6aWpra8nlcpSVlZHL5aitrU07rSRJ\nmgflaQeQuoxevaC2FvbYA/bbjzh4MJ9168bHM2YQgfr6eqqqqgDI5/PpZpUkSc3iyLTU3n7xCxg3\njqsWWoj9Z8zgZWDb4lMNDQ1UV1enmU6SJM0Dy7SUhkUXZd+vvmJj4CvgTuBqoBcwyQ1fJEnqMCzT\nUkoqKip4HBgIHA/sCrwCHLrkkm5JLklSB2GZllJSU1NDJpNhGnASsDYwoayMsz75JFlO7623Uk4o\nSZLmxjItpSSfz1MoFMhms4QQaMhmeWP4cLjgAnjySVh9dfjnP93sRZKkEhZiB/rn5MrKylhXV5d2\nDKntvf02HHAA3HZbsozepZfC2munnUqSpC4jhDAmxlg5t9c5Mi2Vor594ZZbYMQIeOcdqKyEo46C\nhoa0k0mSpCYs01KpCiHZ5OWVV2CffeCMM2C11eCuu9JOJkmSippVpkMIK4QQNi/eXyCEsFDbxpL0\nnSWWgMsugwcfhJ49YbvtYOhQ+OCDtJNJktTlzbVMhxB+C9wKXFo8lAVuactQkmZh003h+efhpJPg\npptgwAD417+gsTHtZJIkdVnNGZk+GBgM/BcgxjgeWLotQ0majQUWgOOOg5degkGDYL/9YOONYezY\ntJNJktQlNadMfxNjnDbzQQihGxDaLpKkufrpT+G++2D4cHjtNRg4EI49Fr7+Ou1kkiR1Kc0p04+F\nEI4EehbnTV8P3N5aAUII3UIIz4UQWu17Sl1CCMmFia++CnvtBaeemqxNPXp02skkSeoymlOmjwS+\nAF4FDgHuA6pbMcMhJLsoS5ofvXrB5ZfD/fdDeTlsvTXk8/DRR2knkySp05tjmS5O6bg8xnhRjPGX\nMcadi/db5YqnEMIKwM/5/uJGSfNr883hhRfg+OPhhhuSCxQLBS9QlCSpDc2xTMcYvwX6hxC6t9Gf\nfw7JyLf/tZdaQ8+ecMIJSaleay344x9ho42SVUAkSVKra840jzeAR0IIx4QQDp55a+kfHEL4BfBR\njHHMXF5XFUKoCyHUTZ48uaV/rNQ1DBiQTPu48kp4441kS/LDDoMvvkg7mSRJnUpzyvQkYDSQAXo3\nubXURsCOIYSJwHXAFiGEq3/8ohhjIcZYGWOs7N27Nf5YqYsIAfbeO1nto6oKhg1LSvbIkdRefTW5\nXI6ysjJyuRy1tbVpp5UkqUMKMca0MxBC2Aw4Isb4izm9rrKyMtbV1bVPKKmzeeop2H9/eO45RpeV\nsV9jI28Wn8pkMhQKBfL5fKoRJUkqFSGEMTHGyrm9rjk7II4OIdzz41vrxJTUbtZfH55+mhOWWIL1\nGxt5GTgO6AE0NDRQXd2ai/RIktQ1zHVkOoSwfpOHPYFfAVNjjH9py2Cz4si01HJlZWUsGyNnAUOB\n14A/AQ+EQKMrf0iSBLTiyHSM8akmt4dijAcDm7RKSkntrqKigveBPYCtSf4SuA+4acEF4f33U80m\nSVJH05xpHos2uS0eQtgSWKIdsklqAzU1NWQyGSC5sngNoKa8nF9MmwYrrwznnAMzZqSaUZKkjqI5\nq3m8DIwtfn2OZPfDP7RlKEltJ5/PUygUyGazhBBYNpslN3w43V55BYYMSZbQGzQIHnkk7aiSJJW8\n5syZ7h5jnP6jY+UxxnYfunLOtNTGYoRbboFDDoFJk5Kl9c44A5ZdNu1kkiS1q1abMw08NYtjT897\nJEklLwTYeWd45RWorobrr0+mfgwb5tQPSZJmYbZlOoSwdAhhLWDBEMIaIYQ1i7chJBu4SOqsMhk4\n+WQYOxY22AAOPTTZRfHRR9NOJklSSZnTyPTPgfOBFYALgQuKt2NJlqeV1NmttBLceSfceCNMmQIb\nbwz77AMffph2MkmSSkJz5kzvFmMc0U555sg501KKvvoKamrgzDNhwQWTkev994fy8rSTSZLU6po7\nZ7pZ24mHELYBViPZtAWAGOMpLUo4HyzTUgl47TU46CAYPRrWXBPOOw82cel5SVLn0prbiV8I7AP8\nGVgQ2AtYscUJJXVMK68Md98NI0cmUz823RT23BPeeSftZJIktbvmrOYxJMa4J/BJjPE4YH0s01LX\nFgLssgu8+iocdxyMGgUDBsCpp8LUqWmnkySp3TSnTH8z82sIYdni4+XaLpKkDiOTgZNOgnHj4Gc/\ng2OPhdVWg9tvTzuZJEntojll+j8hhMWBM4HngYlASVyQKKlE9O8PN92UTP8oL4cddoCf/xxefz3t\nZJIktak5lukQQhlwZ4zxsxjjSKAfsEaM8dh2SSepY9l6a3jxxWTFj0ceSUapjz4avvwy7WSSJLWJ\nOZbpGGMj8K8mj7+OMX7a5qkkdVw9esDhh8P48cmFiaefnly0eM011F59NblcjrKyMnK5HLW1tWmn\nlSSpRZozzeOBEMJObZ5EUuey7LIwfDg8/jj06QP5PLl99mHJ+npijNTX11NVVWWhliR1aM0p0/sC\nN4UQvg4hfBpCmBJCcHRaUvNssAE8/TRHLbkkP21spA4oAEsDDQ0NVFdXpxxQkqT515wy3QvoDiwM\n9C4+7t2WoSR1MmVl/GPKFFYCzib5P/TXgcOB9+vr00wmSVKLzLVMxxi/BXYFjire7wOs3dbBJHUu\nFRUVfA4cAawOPEyyRNCr5eVw223QjN1YJUkqNc3ZAfF8YHNg7+KhBuDitgwlqfOpqakhk8kAMB7Y\nAdhpgQVYcumlYccdYZtt4OWXU80oSdK8as40jw1jjH+kuHlLcTWPHm2aSlKnk8/nKRQKZLNZQghk\ns1l2u+wyFps4EYYNg2eegbXWgoMOgk+9LEOS1DGEOJd/Wg0hPAVsANTFGAeFEJYC7o0xDmyPgE1V\nVlbGurq69v5jJbWHTz6B44+Hiy6CxRZLdlbcb79kExhJktpZCGFMjLFybq9rzsj0BcCNQO8QwonA\no8DpLcwnST+01FJw/vnwwgswaFAyQr3WWnDPPWknkyRptppzAeKVwF9JrhX6FNg1xnhdWweT1EWt\nvjqMHg033wxTpyZzqbffHsaNSzuZJEn/ozkj0wDdgOnAtHl4jyTNnxBgp52SCxL/+c9k45c114QD\nDoDJk9NOJ0nSd5qzmkc1cC2wHLACcE0I4Zi2DiZJLLAA/PnPMGEC7L8//OtfsOKK8I9/JKPWkiSl\nrDmjzL8G1o0x/jXGWA2sR7LngiS1j1694LzzYOxY2GQTOPJIWGUVGDnS9aklSalqTpl+H2h6OX15\n8Zgkta8BA5INXkaPhkUWgd12g403hqefTjuZJKmLak6Z/hR4OYRwaQjhEuAl4OMQwlkhhLPaNp4k\nzcJWW8Gzz8IllyRTQNZfH/baCyZNSjuZJKmLac4607+b0/MxxstaNdEcuM60pP/xxRdw+unJhYoA\nhx8ORx2VjFxLkjSfmrvO9FzLdCmxTEuarUmT4NhjobYWeveGE06AP/wBundPO5kkqQNqtU1bQgjb\nhhCeCSF8FEL4NIQwJYTgXr+SSktFBVx9dTJ/epVVkmX01lgjWa86Rmpra8nlcpSVlZHL5aitrU07\nsSSpE2jONI8JwG4kc6UbZx6PMX7bttH+lyPTkpolxuRCxaOOgldf5aOVV2bXiRN5uMlyeplMhkKh\nQD6fTzGoJKlUteZ24u8Az8cYp8cYv515a3lESWojIcCOO8JLL8HFFxNef52Hpk7lOqB/8SUNDQ1U\nV1enmVKS1Ak0Z2R6PeB44EHgu2GdGOO5bZpsFhyZljQ/Fg2Bw4EjgO7AhcDfgSkh0NjYOMf3SpK6\nptYcmT4R+BZYHOjd5CZJHcKS2SwnACsBVwAHAW8ANYstBt98k2Y0SVIH15wy3TfGuGOMsTrGeNzM\nW5snk6RWUlNTQyaT4X2gClgLeLKsjGM++wxWXhmuugq+dfaaJGneNadM3x1C2KLNk0hSG8nn8xQK\nBbLZLCEEvsxm+eTKK+G++5Ktyn/9axg0CO680+3JJUnzpDlzpqcAiwENwDQgADHGuGTbx/sh50xL\nanWNjTBiBFRXw5tvwqabJpvArL9+2skkSSlqzTnTvUiu2VmMZK50L5wzLamzKCuDoUPhlVfgvPOS\nr4MHwy67wGuvpZ1OklTi5lqmi8vg7QocVbzfB1i7rYNJUrvq0QMOPBAmTEh2T7z7blhtNdhvP3j/\n/bTTSZJKVHN2QDwf2BzYu3ioAbi4LUNJUmoWWQSOPx7eeAP+9Cf497/hJz9JpoF8/nna6SRJJaY5\n0zw2jDH+EfgGIMb4KdCjTVNJUtqWXhrOPRdefRV++Us45RTo3x/OOsvl9CRJ32lOmZ4eQigDIkAI\nYSmabCsuSZ1a//5QWwvPPguVlXD44clyev/+N8yYkXY6SVLKZlumQwjlxbsXADcCvUMIJwKPAqe3\nQzZJKh0DBybzqO+9F5ZZBn73O1h9dRg5MlkRRJLUJc1pZPppgBjjlcBfgTOBKcCuMcbr2iGbJJWe\nLbeEp56CUaOgWzfYbTdYd1246y7XqJakLmhOZTrMvBNjfDnGOCzGeE6McWw75JKk0hVCMo/6xRfh\niivg009hu+1gs83gscfSTidJakflc3iudwjhz7N7MsZ4VhvkkaSOo1u3ZPfEoUPh0kvh73+HIUNg\n++2hpgbWdhVRSers5jQy3Q1YGFhkNjdJEiRrVP/pT8ka1aedBk88kcyx3mMPbj3zTHK5HGVlZeRy\nOWpra9NOK0lqRbPdTjyE8GyMcVA755kjtxOX1CF89hmceSYzzjwTpk7lcuAk4B0gk8lQKBTI5/Mp\nh5QkzUlrbCce5vCcJGl2Fl8cTj6Zwb17cwHwa2ACcC6waEMD1dXV6eaTJLWaOZXpLdsthSR1Qs++\n+y6HAisBVwD7A28CB9XXw+TJqWaTJLWO2Zbp4k6HkqT5VFFRAcDbwB+BlYERwKEA/frBsccmK4FI\nkjqs5uyA2CZCCH1DCA+EEMaFEF4OIRySVhZJags1NTVkMpnvHr8J/CmT4T9nnAE77JBcrNivH5xw\nAnz+eWo5JUnzL7UyDcwADo8xrgoMBg4IIayaYh5JalX5fJ5CoUA2myWEQDabpVAosMNf/gLXXpus\nU73VVnDiiUmpPvVU+PLLtGNLkubBbFfzaG8hhFuA82OMo2f3GlfzkNQpPfcc/O1vcPvt0Ls3HHVU\nstTeggumnUySuqzWWM2j3YQQcsBA4Kl0k0hSCgYOhNtuS9anXnttOOII6N8fhg2Dr79OO50kaQ5S\nL9MhhIWBG4FDY4z/ncXzVSGEuhBC3WSvfpfUmQ0eDPfcAw8/DKuuCoceaqmWpBKXapkOIXQnKdK1\nMcZRs3pNjLEQY6yMMVb27t27fQNKUho23hjuuw8eeuiHpfqccyzVklRi0lzNIwCXAa/EGM9KK4ck\nlaxNNvlhqT7sMEu1JJWYNEemNwL2BrYIITxfvG2fYh5JKk2WakkqWamV6RjjozHGEGNcM8a4dvH2\nn7TySFLJs1RLUslJ/QJESdI8mlWp7tcP/vlP+OqrtNNJUpdimZakjqppqV5ttWRJvVwu2fzlv/+l\ntraWXC5HWVkZuVyO2tratBNLUqdjmZakjm5mqX7sMVh3XTj2WKb26cPEfffl8/p6YozU19dTVVVl\noZakVmalaWoIAAAb7ElEQVSZlqTOYsMN4T//gWee4cEYqZ4xg3rgFKAX0NDQQHV1dcohJalzsUxL\nUmdTWcl233zDmsCdwFHAROBMYFp9fZrJJKnTsUxLUidUUVHBS8BQYFWS3bEOAd4COOggePvtFNNJ\nUudhmZakTqimpoZMJgPAa8A+wNo9e/L2ppvCxRfDT34Cf/gDTJiQak5J6ugs05LUCeXzeQqFAtls\nlhAC2WyWYy69lBUffBDeeCMp0lddBSuvDEOHwgsvpB1ZkjqkEGNMO0OzVVZWxrq6urRjSFLn8MEH\nyYYvF14IX3wB228Pxx4LG22UdjJJSl0IYUyMsXJur3NkWpK6qmWXhdNOg0mT4OST4emnYciQZKm9\nu+6CDjTYIklpsUxLUle3+OJQXQ319TBsGEycCNttB+usAyNGwLffpp1QkkqWZVqSlMhk4OCDk4sS\n//1vaGiA3XeHVVaByy6DadPSTihJJccyLUn6oR494De/gZdfhhtugEUWgd//PlkB5KyzkvnVkiTA\nMi1Jmp1u3eBXv4K6Orj7blhxRTj8cOjbN7lQ8YMP0k4oSamzTEuS5iwE2HpreOCB5CLFn/0suXAx\nm4WqKhg/Pu2EkpQay7QkqfnWXRdGjkwK9G9/C1deCQMGJCPYTz0FQG1tLblcjrKyMnK5HLW1tSmH\nlqS24zrTkqT59+GHcN55cMEF8NlnfDhgAPu/9RY3TZ363UsymQyFQoF8Pp9iUEmaN64zLUlqe8ss\nk6xRPWkSnH02M15/nVFTp/IisDfQHWhoaKC6ujrloJLUNizTkqSWW2QROPRQ+n37LXsXD10JvAkc\nCXxeX59eNklqQ5ZpSVKrWS6b5WpgTWA74DXgdOCdEOCQQ+Ctt1LNJ0mtzTItSWo1NTU1ZDIZAO4C\ntgI26NmTjzbcEC68MFleb9dd4cknU80pSa3FMi1JajX5fJ5CoUA2myWEQDab5cBLL6Xfo48mo9JH\nHAGjR8MGG8BGG8GoUW5XLqlDczUPSVL7+vLLZLvyc85JCvZPfgKHHprsurjQQmmnkyTA1TwkSaVq\n4YXh4IOTtapHjoTeveGgg77fWfHdd9NOKEnNZpmWJKWjvBx22QWeeAIeeww23zzZWTGXgz33THZb\nlKQSZ5mWJKVvww3hxhthwoRklPqOO2D99ZPj118P06ennVCSZskyLUkqHf37w1lnwTvvwLnnwkcf\nwdChyfHTT4dPP007oST9gGVaklR6FlkkGaEePx5uuw1WXhmOPhpWWAH22w9eeSXthJIEWKYlSaWs\nrAx+8Qu491548UXI52H4cFh1Vdh2W7jzTmhsTDulpC7MMi1J6hjWWAMuuQTefhtOPjkp19tvD6us\nkkwJ+fxzamtryeVylJWVkcvlqK2tTTu1pE7OMi1J6lh694bqapg4EWprYaml4JBDmL7MMny1775k\n6uuJMVJfX09VVZWFWlKbskxLkjqmHj2SJfQefxzq6ri5vJxfz5jBOOBeYCdgakMD1dXVKQeV1JlZ\npiVJHd8667B7QwN9gaOBlYCbgTeAPerr4eOPU40nqfOyTEuSOoWKigo+Bk4H+gO/JCnTp0KyCshv\nfgNjxqSYUFJnZJmWJHUKNTU1ZDIZAL4lGZneIZPh9tNOS4r0iBFQWQkbbABXXQXffJNqXkmdg2Va\nktQp5PN5CoUC2WyWEALZbJZCocAvjjoKLroI3n0Xzj472fjl179ORqv/8hd44420o0vqwEKMMe0M\nzVZZWRnr6urSjiFJ6shihPvvTwr2zTfDt9/CNtvA/vvDz38O5eVpJ5RUAkIIY2KMlXN7nSPTkqSu\nJQTYcku44Qaor4cTToCXXoKdd4Z+/ZI1rD/4IO2UkjoIy7Qkqetafnk4/vikVI8alWwAc9xx0Lcv\n7LYbPPBAMpItSbNhmZYkqbwcfvlLuOceGD8eDj442cJ8iy2SrcvPOQc++STtlJJKkGVakqSmVloJ\n/vnP5ILF4cNhscXgsMOSUey99oKHH3a0WtJ3LNOSJM3KggvCPvvAk0/CCy/AH/4At98Om26aTAc5\n66zvNoOpra0ll8tRVlZGLpdzC3OpC3E1D0mSmquhAUaOhEIh2ca8Rw8mrrMO+z37LHdPnfrdyzKZ\nDIVCgXw+n2JYSS3R3NU8LNOSJM2PsWPhkkv4/PzzWayxkfFAAbgC+BjIZrNMnDgx1YiS5p9L40mS\n1JZWXx2GDWO5xkZ+DXwEnAm8C1wPDKivT9awltSpWaYlSWqB3tksVwEbA6sBFwJbAndBsm713/4G\nb72VYkJJbckyLUlSC9TU1JDJZAAYBxwGrLjggjxy4IHJhYonnwz9+8NWW8G118I336SaV1LrskxL\nktQC+XyeQqFANpslhEA2m+X8Sy5h4/POg7vvhokT4cQT4Y03YM89oU8fOPBAeO65tKNLagVegChJ\nUntobEx2VLzssmS3xalTYeBA+O1vIZ+HJZZIO6GkJrwAUZKkUlJWBltuCddcA++9B+efnxw/6KBk\ntHr33eHOO2HGjHRzSponlmlJktrbkkvCAQfAs88mtz/8Ae67D7bfHvr2hSOPhJdfTjulpGawTEuS\nlKaBA+G885LR6lGjYL314Oyzk6X31l0XLrgAPvkk7ZSSZsMyLUlSKejRA375S7jlFnj33aRQz5iR\nXKzYpw/ssgvcdhtMn552UklNWKYlSSo1Sy8Nhx6arPjx/PNJoX7kEdhxR1hhBfjzn+H556m9+mpy\nuRxlZWXkcjlqa2vTTi51Oamu5hFC2BYYBnQDLo0xnjan17uahySpy5o+He66C664Am69FaZPZ1wI\nXBEj1wDvAJlMhkKhQD6fTzut1OE1dzWP1Mp0CKEbMB74GcnfAc8Ae8QYx83uPZZpSZKATz6h+qc/\nZbtPP2UI0Ag8BFwNPN23Ly9NmpRuPqkT6AhL460HTIgxvhljnAZcB+yUYh5JkjqGpZbi1ClT2Bjo\nDxwPLAdcBjzz9tuw227J6PW0aanGlLqCNMv08sDbTR6/Uzz2AyGEqhBCXQihbvLkye0WTpKkUlZR\nUQHAW8DJwABgXeDaRRaBBx+EnXZKLlz805/giSegA23SJnUkJX8BYoyxEGOsjDFW9u7dO+04kiSV\nhJqaGjKZzA+Ojctk6HHRRclqIHfcAVtvDZdfDhtuCCuuCNXVrl8ttbI0y/S7QN8mj1coHpMkSXOR\nz+cpFApks1lCCGSz2e8vPuzePdkA5tpr4cMPk0L9k5/Aaacl61evuSaceiq89VbaP4bU4aV5AWI5\nyQWIW5KU6GeAPWOMs/1fZi9AlCSpBT78EEaOTEr2448nx9ZfH/bYI5ln3adPuvmkElLyFyDGGGcA\nBwJ3A68AI+ZUpCVJUgsts0yyZvVjjyWj0qedBt98k6xpvcIKsNVWcNllMGVK2kmlDiPVdabnlSPT\nkiS1gXHj4LrrkhHrCROSaSLbbQe77w477ACLLJJ2QqndlfzItCRJKhGrrgonnQTjx8Mzz8BBB0Fd\nHeTz0Lt3ss35tdfCF19895ba2lp3X5RwZFqSJM1KY2Myr3rkyOT2/vvQsydsvz2P9unD//3730z+\n+uvvXu7ui+psSn4HxPlhmZYkKQWNjck86xEj4IYb4IMP+Bq4AxhR/NoAZLNZJk6cmGZSqdU4zUOS\nJLWOsjLYeGM47zx45x02JdltcSOSMj25+HX9+nr48ss0k0rtzjItSZKar1s36rNZDiLZIGJT4HJg\nE+B6gF69kt0Xr7zSVUHUJVimJUnSPJm5+2Ij8DDJOrcrLrggo//6V/jjH+HZZ2GffWDppWGbbeBf\n/0rWuJY6Icu0JEmaJ7PaffHiSy7hZ3//OwwbBpMmwVNPweGHw5tvwn77JRvCbLpp8vzbb6f9I0it\nxgsQJUlS24kRXnoJRo2CG2+EsWOT4+uuC7/6Fey8M6y8croZpVlwNQ9JklR6xo9PivWoUcma1pCU\n6Z12Sm6DBycXPEops0xLkqTS9vbbcOutcMst8MADMGNGsuX5DjskxXrLLWHBBdNOqS7KpfEkSVJp\n69sXDjgA7rkHJk+Ga66BzTaD669PCnWvXslUkCuvhE8+Adx5UaWnPO0AkiRJLL447LFHcps6FR58\nMBmxvvXWZEpIt258uNJKvPjGG3SfPp0I1NfXU1VVBeDOi0qN0zwkSVLpihHGjIGbb+bVM85gwPTp\nALwG3F68vV1RwYT6+jRTqhNyzrQkSepUysrKqIiRnwO/ALYAFgA+BxbbbTf4xS9gu+2S6SFSCzln\nWpIkdSoVFRXUAxcC2wNLATsDdy68MDz8MPz618lGMRttBKeemizJ14EGDdUxWaYlSVKHMHPnxZm+\nAkZnMnx78cXw7rvJUnt/+1sy5/rYY2HNNSGXSzaNueUW+OKL1LKr83KahyRJ6jBqa2uprq5m0qRJ\nVFRUUFNTM+uLD997D/7zH7jjDrj3XvjyS+jeHTbeGLbdNpkOstpqEEL7/xDqEJwzLUmSBDBtGjz2\nGNx5Z3KbuQtj377fF+stt4RFF003p0qKZVqSJGlW3n4b7r47KdajRyfTP8rLk7nW220H22yTTBFx\nJ8YuzQsQJUmSZqVvX/j97+HGG5PNYB58EI44Aj77DI4+GgYOhD59IJ+HK65IpozghjGaNUemJUmS\nZnrvvWSO9T33JKPWH30EwGcrrMDVH3zAHTNm8DDQAGQyGQqFghvGdFJO85AkSWqJxsZkeb3Ro3n4\nuONY95tvWBCYCjwGjAZeWnZZbn/3XaeEdEKWaUmSpFZSVlbGAjEyBPgZsDWw9swnl1wSNt8cttgi\nuZDxpz91lZBOwDnTkiRJraSiooJvgHuBo4CBwDLAIUstBTvuCE8/DQccAAMGwAorwN57w/DhMGlS\nmrHVDizTkiRJc/HjDWMAvsxkWG/YMLj8cqivhwkT4F//Stayvvtu+M1vIJuFFVeEqiq4/vrv5mCr\n83CahyRJUjM0e8MYSLYxHzsW7r8f7rsPHnoI/vvf5LnVV4fNNoNNN4VNNkm2QFfJcc60JElSqZgx\nA559NinW998Pjz8ODQ3Jc6uskhTrmbc+fX7w1nkq8Wo1lmlJkqRSNW0ajBmTjFg/9BA8+miy5TnA\nSit9N3I96pNP2PuYY2iYWbxxSb72YpmWJEnqKGbMgOee+75cP/IIfP45AG8ADwOPFm/jgWw2y8SJ\nE1OL2xVYpiVJkjqqb7+FF1/kz4MGsQkwBOhVfGoySan+5T/+AUOGwKBB0KNHalE7K8u0JElSB5fL\n5aivrwdgZZJSvRGwWXk5/WbMSF7Usyesv35SrIcMgQ02gMUWSytyp+E605IkSR1c0yX5XgMuAw7M\nZHh8+HB4/3244QbYb7/kYsbTToPttoMlloC11kqODx8Or72WrC6iNuHItCRJUglr9moeX30FTz2V\nXMz42GPJ/eK8a5ZYAgYPTkatBw+G9dZz9HounOYhSZLUlTU2wquvwhNPwJNPJl/HjUtGqUOAVVf9\nvlxvsEGye2PZDyctdOVl+SzTkiRJ+qHPP0+2Pp9ZsJ98EqZMSZ5bdFFYZx1Yd11Ybz1ueucd9jrm\nGBq+/vq7t3elZfks05IkSZqzxkYYPz4p1888kxTtF1+E6dMB+BB4GnimyW3hLrIsX3PLdHl7hJEk\nSVIJKitLpncMGAC/+U1y7Jtv4MUXOWj99akE1gV+zverVrxZXw+7756MYK+zDgwcCIsvnk7+EmCZ\nliRJ0vd69oT11uO2bJbzi8vyLQKsQ1KsN8tk6P/kkzBixPfv6d8/We+66a137zTStzvLtCRJkv5H\nTU0NVVVVNDQ08AXwIPB0JsNahQLk8/DRR8mujc8++/3thhu+/wZ9+/5vwe7TJ7n4sROxTEuSJOl/\nzLzIcLareSy9NGyzTXKbacoUeP75HxbsW2/9fp3rXr2SNbCb3lZZ5X92cOxIq4h4AaIkSZLazpdf\nwgsvwJgxycWNL7wAY8cmc7MBysuTQl0s1/d9/DG/O/dc6lNeRcTVPCRJklSaZsyA119PinXT23vv\nffeS94HHgV2Kj7PtvIqIq3lIkiSpNM0cjV5lFRg69PvjH3/Mlr17syaw1o/eMmnSpPZM2GyWaUmS\nJJWGXr14I5vl/uIqIk1VVFSkEGjuyub+EkmSJKl91NTUkMlkfnAsk8lQU1OTUqI5s0xLkiSpZOTz\neQqFAtlslhAC2Wy2pLcw9wJESZIk6UeaewGiI9OSJEnSfLJMS5IkSfPJMi1JkiTNJ8u0JEmSNJ8s\n05IkSdJ8skxLkiRJ8ymVMh1C+EcI4dUQwoshhJtCCIunkUOSJElqibRGpkcDq8cY1wTGA8eklEOS\nJEmab6mU6RjjPTHGGcWHTwIrpJFDkiRJaolSmDP9W+DOtENIkiRJ86q8rb5xCOFeYNlZPFUdY7yl\n+JpqYAZQO4fvUwVUFR9+GUJ4rbWzdgC9gI/TDtHBeQ5bxvPXMp6/lvH8tYznr2U8fy3Tkc9ftjkv\nCjHGtg4y6z84hH2BPwJbxhgbUgnRQYQQ6pqzN7xmz3PYMp6/lvH8tYznr2U8fy3j+WuZrnD+2mxk\nek5CCNsCRwKbWqQlSZLUUaU1Z/p8YBFgdAjh+RDCxSnlkCRJkuZbKiPTMcYV0/hzO7BC2gE6Ac9h\ny3j+Wsbz1zKev5bx/LWM569lOv35S23OtCRJktTRlcLSeJIkSVKHZJkuESGEJUMIo0MIrxe/LjGb\n131bnGf+fAjh1ibH+4UQngohTAghXB9C6NF+6dPXnPMXQlg7hPBECOHl4lb2uzd5bngI4a0m53bt\n9v0J0hFC2DaE8Frxc3P0LJ5foPh5mlD8fOWaPHdM8fhrIYRt2jN3qWjG+ftzCGFc8fN2Xwgh2+S5\nWf4udyXNOH/7hhAmNzlPv2/y3D7F3/fXQwj7tG/y0tCM83d2k3M3PoTwWZPn/PyF8O8QwkchhLGz\neT6EEM4tnt8XQwiDmjzn52/u5y9fPG8vhRAeDyGs1eS5icXjz4cQ6tovdRuJMXorgRtwBnB08f7R\nwOmzed2Xszk+AhhavH8xsH/aP1OpnT/gp8BKxfvLAe8DixcfDwd2SfvnaOdz1g14A+gP9ABeAFb9\n0Wv+BFxcvD8UuL54f9Xi6xcA+hW/T7e0f6YSPH+bA5ni/f1nnr/i41n+LneVWzPP377A+bN475LA\nm8WvSxTvL5H2z1Rq5+9Hrz8I+HeTx13681c8B5sAg4Cxs3l+e5JN5QIwGHiqeLzLf/6aef42nHle\ngO1mnr/i44lAr7R/hta6OTJdOnYCrijevwLYublvDCEEYAvghvl5fycx1/MXYxwfY3y9eP894COg\nd7slLD3rARNijG/GGKcB15Gcx6aantcbgC2Ln7edgOtijFNjjG8BE4rfryuZ6/mLMT4Qv1/+80lg\nhXbOWMqa8/mbnW2A0THGT2OMU4DRwLZtlLNUzev52wO4tl2SdRAxxoeBT+fwkp2AK2PiSWDxEEIf\n/PwBcz9/McbHi+cHOvnff5bp0rFMjPH94v0PgGVm87qeIYS6EMKTIYSZhXEp4LMY44zi43eA5dsw\naylq7vkDIISwHslozhtNDtcU/0nq7BDCAm2Us5QsD7zd5PGsPjffvab4+fqc5PPWnPd2dvN6Dn5H\nMso106x+l7uS5p6/XxV/L28IIfSdx/d2Zs0+B8XpRf2A+5sc7uqfv+aY3Tn28zfvfvz3XwTuCSGM\nCclO1x1aKkvjdVVhDlusN30QY4whhNkts5KNMb4bQugP3B9CeImk4HR6rXT+KI4sXAXsE2NsLB4+\nhqSE9yBZxuco4KTWyC2FEPYCKoFNmxz+n9/lGOMbs/4OXdZtwLUxxqkhhD+S/CvJFiln6oiGAjfE\nGL9tcszPn9pFCGFzkjI9pMnhIcXP39Ike468Whzp7pAs0+0oxrjV7J4LIXwYQugTY3y/WPY+ms33\neLf49c0QwoPAQOBGkn9+Ki+OHq4AvNvqP0DKWuP8hRAWBe4Aqov/bDfze88c1Z4aQrgcOKIVo5eq\nd4G+TR7P6nMz8zXvhBDKgcWAT5r53s6uWecghLAVyf/wbRpjnDrz+Gx+l7tSmZnr+YsxftLk4aUk\n10bMfO9mP3rvg62esLTNy+/gUOCApgf8/DXL7M6xn79mCiGsSfK7u13T3+cmn7+PQgg3kUxb6rBl\n2mkepeNWYOYVwfsAt/z4BSGEJWZOPwgh9AI2AsbFZDb/A8Auc3p/J9ec89cDuIlkDtwNP3quT/Fr\nIJlvPcurkzuZZ4CVQrISTA+S/+D++Kr+pud1F+D+4uftVmBoSFb76AesBDzdTrlLxVzPXwhhIPAv\nYMcY40dNjs/yd7ndkpeG5py/Pk0e7gi8Urx/N7B18TwuAWxdPNaVNOf3lxDCAJKL5J5ocszPX/Pc\nCvy6uKrHYODz4sCLn79mCCFUAKOAvWOM45scXyiEsMjM+yTnr2P/NzftKyC9JTeSeaj3Aa8D9wJL\nFo9XApcW728IvERy1fZLwO+avL8/SZmZAIwEFkj7ZyrB87cXMB14vslt7eJz9xfP6VjgamDhtH+m\ndjpv2wPjSUakqovHTiIpfwA9i5+nCcXPV/8m760uvu81klGH1H+eEjx/9wIfNvm83Vo8Ptvf5a50\na8b5OxV4uXieHgAGNHnvb4ufywnAb9L+WUrx/BUfnwCc9qP3+flLzsO1JKs6TSeZ9/w7/r+9Owax\n4grDMPx+WgVBRGNlpVbKokEsUsgiBmHtBLHRQkixbJ0mhFhoYxtIkBTqCulCSJoQCClU1MJCDdkE\nTBUhpYpC0BUt/C3mKDcLrjh63TvwPjAw/HOGM2eYO3xczr0H5oC5djzA6XZ//wR2j5zr8/f6+3cW\neDDy/rve6lvas/dH+3x/udJjedvNFRAlSZKknpzmIUmSJPVkmJYkSZJ6MkxLkiRJPRmmJUmSpJ4M\n05IkSVJPhmlJmhDt/2yvJjkwUjuc5NcVuJZLSXa/734laWhcAVGSJkRVVZI54IckF+ne0aeAmXH2\nO7J6qiTpDRmmJWmCVNVfSX4GPgfW0K3Y+b9lnpPM0IXs1cC9qvokyXpgnm5BhEVgtqoWlqmfALa2\n+r9JPgXOAzuBv4EPWl+rgXN0CyAVMF9VX43zHkjSkBimJWnynARuAk/pQuxLSTYCZ4DpqrrdwvKL\nc36vqoNJ9gHfAR8tUwfYDuypqsdJPgMWq2pbkh2tf1rbTVU11fpfN6YxS9IgGaYlacJU1aMk3wMP\nq+rJksMfA5er6nZre7/V9wCHWu1Ckg1J1i5Th25588dtfxr4urVbSLLQ6v8AW5J8A/wC/PauxytJ\nQ+YPECVpMj1r2zg9el2DqnpAN/XjEjAHnB3zNUnSoBimJWlYrgHTSTYDjEzzuAIcbbW9dHOp/1um\nvtRl4EhrNwXsaPsfAquq6kfgOLBrLKOSpIFymockDUhV3U0yC/yUZBVwB9gPnADm2/SMReBYO+VV\n9aW+Bc4nuQXcAm60+qZWf/HlyxfvdkSSNGypqpW+BkmSJGmQnOYhSZIk9WSYliRJknoyTEuSJEk9\nGaYlSZKkngzTkiRJUk+GaUmSJKknw7QkSZLUk2FakiRJ6uk5rpTGEM+wlWkAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f448e1ee290>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Numerical flux at y =  1.3 is -8.04936746828\n",
      "Exact     flux at y =  1.3 is -8.20909090909\n",
      "\n",
      "Abs. error = 4.497e-04\n",
      "Rel. error = 2.647e-05\n"
     ]
    }
   ],
   "source": [
    "# create numpy arrays for analytics\n",
    "yvals       = np.zeros(mesh.specialSets['MinI_VertexSet'].count)\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": 11,
   "metadata": {
    "collapsed": true
   },
   "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": 12,
   "metadata": {
    "collapsed": true
   },
   "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": 13,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "solver = uw.systems.Solver(ss)\n",
    "solver.solve()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# create numpy arrays for analytics\n",
    "yvals       = np.zeros(mesh.specialSets['MinI_VertexSet'].count)\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": 15,
   "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.1052631579\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": 16,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAtEAAAF3CAYAAABjZBdpAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzs3Xd4leX9x/H3HZbGPdC6kmAdtXVVUdFaF4q4ilp3HLji\ntq3+nKm70apttW5SUbQet9a9cO8qCI6KoiIB1Cq4MQpI7t8f96FGBEkgyXPOyft1XefKyXNOkk+e\nK4Rv7vN9vneIMSJJkiSp5cqyDiBJkiQVG4toSZIkqZUsoiVJkqRWsoiWJEmSWskiWpIkSWoli2hJ\nkiSplSyiJUmSpFayiJYkSZJaySJakiRJaiWLaEmSJKmVumYdoCWWXHLJWFVVlXUMSZIklbjhw4dP\nijH2nNPziqKIrqqqYtiwYVnHkCRJUokLITS05Hm2c0iSJEmtZBEtSZIktZJFtCRJktRKFtGSJElS\nK1lES5IkSa1kES1JkiS1kkW0JEmS1EoW0ZIkSVIrWURLkiRJrWQRLUlqE7lcjqqqKsrKyqiqqiKX\ny2UdSZLaTVFs+y1JKmy5XI6amhoaGxsBaGhooKamBoDq6uoso0lSu3AlWpI0d2KEr76CDz9k0PHH\ns1JjIxsCWwH9gA0aG7nj2GPh+edh2DB4+WV4/XV46y14912YMAE+/BA+/hi+/nqeorgKLqmjuRIt\nSSUml8tRW1vLuHHjqKiooK6ubs6rwTHCp5+mwnbCBBg//rv7H3wAkyf/8NbYmD4OeHJ2n/fDD2HD\nDVsWvLwcevaEpZZKb2e+NT++1FKwwAL/+35dBZfU0ULM/wIsZL17947Dhg3LOoYkFbyZC0qA8vJy\n6gcNorpvX3j1VRgz5oeF8oQJqShurqwMll0WllkGFl4YFlwwFa4LLvj92wILcMyppzL244+ZDHwF\nRKAbsNxSS3H9NdfAtGnw7bfp7Yxb8/cbG2HixO/fPvoovf3mm1l/sz17wiqrcPPIkYz46itGA6OB\nt4FvgMrKSsaOHdv2J1lSSQshDI8x9p7T81yJlqQSUltbS1ljI32ANWbcGhtZc7/9oKnpuyd26ZIK\n5OWXh7XXhu23T/dXWCG9XX55+MlPoGvL/ptYd5FFGDSL4r3mb3+D/v3n/hua0TLSvKieODGtcL/z\nDowezcZffcVuM31YA/BWQwMcfjisskq6/fznUFkJIcx9HknKs4iWpGLV1ASjRsErr6QV5ldf5bGG\nBno1e8qXwGvAbU1NHHzRRbDGGrDyyqlA7tKlzaLMaJtodRvJnITw3Yp3r16zfMpGVVV83NDAysAq\n8L+3a3TvDjfcAJ999t2Te/aE9deHDTZIb9dbDxZffN4ySuqUbOeQpGIxdSoMHw5PPZVuTz/9XYHY\ntSusuip3vfsuzzc28iqpeG4gtVaUcmvDbFtY6uup3muvdOHim2+mPzZeeCHdRo36Xz83K62UCuoZ\nxfXaa8N882X03UjKmu0cklTsvvoKnnvuu6L5+ee/m2Kx6qqwyy6w8cbwy1/Cz34G3bvzZS7H32dR\nUNbV1WX0TbS/Oa6CL7lkuv3qV3DYYenYF1+kiSEziuonnoDrr0+Pde0Ka62VCuottoC+fWHRRTP4\nziQVMleiJalQfPZZKuZmFM0vvZQuvisrS6ujv/41bLJJKpyXWmq2n2aupnMI3nvvu6J6xm3y5NT2\nssEGsPXW0K9fagHJt8J4rqXS09KVaItoScrSmDFw991w113w5JOpaO7RI7UWzCiaN9wwTcdQx5o2\nLa3+P/hgug0fnlpAFlsMttyS5xdZhH2vu463mk0P+V8biYW0VLQsoiUpQ7NdoWxqSiucd92ViufX\nXksf8POfw29+A9tum1Y67cktPJMmwcMPp4L6oYfg/fcBeB14MH97Ali6hPvPpc7AIlqSMjLzhW7l\nwHbdu3P2hhuy0htvpPFsXbqkVebf/AZ22AF++tNsQ6t1YmSNsjL6kXZn3ASYnzQN5U5g77vuSq0f\nPXpkmVLSXLCIlqSMVFVV8VlDA7sAA4AtSQXWFyGw8G67pcJ5m21SW4CKVlVVFQ0NDQDMB2wK7Azs\nWlbGYk1NsMgisNNOsPvu6eLEbt2yjCuphVpaRJd1RBhJ6hSmTYN77+W8hgb+C1wJrA7UA32BnjHC\njTfCXntZQJeAuro6ysvLgbRD4oPAH8rLuf+qq+C++2DHHeH229MfTMssA4ccAo8+CtOnZ5pbUtuw\niJakeTVyJBxzTNrlb/vt2bKsjHpgPWBF4PfAo8AylZWZxlTbqq6upr6+nsrKSkIIVFZWUl9fz177\n7ZcK5yFDUuvOHXek1o5cLq1IL7ccHHlkmvPdfBdJSUXFdg5JmhsffJDmCl97bdrEo1u31Nu8777c\n8NlnHHT44bPe/MOpDZ1XY2Naob7xRrj3Xvjmm7TN+oEHwkEHpeJaUuZs55Cktvb116kA2mabtOr8\nf/8H888Pl12WiurbboMBA9hzv/1muUJpAd3JlZenDXJuvRU++iitTK+2Gpx+OlRWpvaPBx5wdVoq\nEq5ES9KcjB4NF1+cVp2/+AIqKmCffdJt1VWzTqdiN2YM/OMfMHgwTJwIvXrBwQfDAQfA0ktnnU7q\ndJzOIUnzoqkJhg6Fv/8d7r8fundPUxb23x823TTtIii1palT4V//gkGD4LHHUovQTjvBoYfCZptB\nCFknlDoF2zkkaW5MnpzaM37xC+jfH0aMgDPOgHHj0kr05ptbQKt9zPhD7dFHYdSodPHh0KGwxRbw\ns5/B3/4GH3+cdUpJef5PIEkA776bepyXXx6OOAIWWgiuuw4aGuDUU31ZXR1rRtH83nvpj7cll4Rj\nj00XHx50EHf95S9UVVVRVlZGVVUVuVwu68RSp9M16wCSlJkY4YknUsvGXXell8t32QV+9zvo08eX\nz5W9+ef/rv/+1Vfhssv4dvBgtp82jSnAn4GXGhqoqakB8OJVqQPZEy2p85k2LU1GuOCCNJ5uiSWg\npgYOPzytREsFrPcKK7DThAkcASwKDAXOAcZUVDA2v4OipLlnT7QkAblc7n8ve69UWcnzBx+cJmrs\nv3+6ePAf/4Dx4+Hssy2gVRReeu89/ghUAMeTdsV8FLh53Li0Q6Ij8qQOYREtqWTlcjlqamp4v6GB\nA2Nk6Lhx9LnySj4GuPvutAp90EHpJXOpSFRUVADwJXA+0AuoAZbq2hV++1v4+c/h6qvTtA9J7cYi\nWlLJOv3kk6lubGQ08A9gIrAtsO706bD99vY8qyjV1dVRXl7+v/enALnycp696iq46aa0qcsBB8CK\nK6aWpcmTswsrlTCLaEmlZ8oUuOIKHh43jnrgQ2AbYAPgfmDc+PGZxpPmRXV19Sx3xNxrn31gt91g\n+PC08+HKK8Mxx0BVVSqmp0zJOrpUUrywUFLpmDIFrroKzjkHxo/npe7dOXnqVB6c6WmVlZWMHTs2\ni4RSx3r++TSicejQtNPmmWfC3ntDly5ZJ5MKlhcWSuo8pk5NG6SstNJ3EzYefJBRgwfzVLOXvQHK\ny8upq6vLKKjUwfr0gYcegocfhqWWgoEDYa210kjHIlhEkwpZuxXRIYSrQggfhRBem8Vjx4YQYghh\nyfb6+pI6gRjTNsm/+EXaIKWyMhUMzzwD/fpRvffes3zZ21m66nT69oUXXoBbbkkjHgcMgF//Gp5+\nOutkUtFqz5XoIUD/mQ+GEFYA+gHj2vFrSyp1w4bBppvCzjtDjx5w333w1FOw1Vbfu2CwurqasWPH\n0tTUxNixYy2g1XnN2Ezotddg0CAYMyYV0jvskDZykdQq7VZExxifBD6ZxUMXkEZb+jqSpNYbPx72\n3RfWWw/eeAOuuAJGjoRttnHahtQS3bqlzYXefhv+/Oe0Gr3WWunfldcKSC3WoT3RIYQBwHsxxpc7\n8utKKgGTJ8Mpp8Aqq8DNN8NJJ6Ui4JBDoGvXrNNJxae8HE44Ia1IH3dcavVYZZW07f3HH2edTip4\nHVZEhxDKgZOBU1v4/JoQwrAQwrCJEye2bzhJhWv6dBg8OI3r+tOfYKed4M030w6DCy+cdTqp+C22\nGJx7bvqjdOBAuPTSVEwPGpT+/UmapY5cif4paWOll0MIY4HlgZdCCD+Z1ZNjjPUxxt4xxt49e/bs\nwJiSCsbDD8M666RdBXv1gueeg+uvTxcQSmpbyy0H9fWpPWqNNeDQQ9N0jxdeyDqZVJA6rIiOMb4a\nY1wqxlgVY6wCJgDrxBj/21EZJBWJN95IOwputRV88UXahe2ZZ9J/6JLa1+qrw2OPQS4H772X/t0d\nfDBMmpR1MqmgtOeIuxuA54BVQwgTQggHttfXklQivvkmbQyx5ppp0sZ558GoUWkXNi8alDpOCLDX\nXukP2mOOgauvTi0eV1zB9f/8J1VVVZSVlVFVVUUul8s6rZQJdyyUVBgefTS9fPzWW2lHtb/+NW0O\nISl7//kPHHkkPP44I0LgsBj5d/6h8vJy56+rpLhjoaTiMHEi7Ldf2gyiqSltT/zPf1pAS4XkF7+A\nRx/lqCWXZKkYeR64ElgSaGxspLa2NuOAUseziJaUjRhhyBBYbTW44QaorU0bPmy5ZdbJJM1KCFz6\n8cf8DDgP2BcYDRwGTGhoyDSalAWLaEntLpfLfa+H8q7zz4fNN4f994ef/QxGjEjj6+afP+uokn5E\nRUUFk4ETgLWAl4DLgJe6d3fXQ3U6FtGS2lUul6OmpoaGhga6xcjAhga2Pv54prz4Yhqn9eST6aVi\nSQWvrq6O8vJyAEYBWwL7du/OKvPNB+uum/4YnjYt04xSR7GIltSuamtraWxsZFPgFeB04FbgV4sv\nnsZmlflrSCoW1dXV1NfXU1lZSQiByspKtr7qKuZ7+2347W/TrqIbbAAvuzGxSp/TOSS1q8VD4K/A\n/sAYUv/kQ0AIgaampkyzSWpjt98Ohx0Gn3wCf/wjnHQSdO+edSqpVZzOISl7Dz/Mf7p0YR/gHGB1\nUgENqbdSUonZeWd4/fU02/3002H99dM1D1IJsoiW1PYaG+Hoo2GrrShfemk269GDk4Gv8w+Xl5dT\nV1eXZUJJ7WWJJdJuh3fcAR9+mArpU0+FqVOzTia1KYtoSW3rxRdhnXXg4ovh6KNZ5O23OWzw4O/1\nULoxg9QJDBiQNmnZc0846yzo3RuGD886ldRm7ImW1DamTYO6unR1/jLLpBnQfftmnUpSIbjnHjjk\nkLQyfcIJaWW6R4+sU0mzZE+0pI7zxhuw0UZwxhlp1enVVy2gJX1n++3htddgn33g7LPTODwneKjI\nWURLmntNTXDRRfDLX8KYMXDLLWnL7kUXzTqZpEKz2GJw9dVw331pesf666e2ryJ4RVyaFYtoSXNn\n/HjYemv43e/S7oOvvQa77JJ1KkmFbptt0ir0llumC5AHDIBJk7JOJbWaRbSk1okxXXm/xhrw7LNw\nxRVw772pD1qSWqJnz9QnfeGF8OCDsNZa8NhjWaeSWsUiWlLLffklVFfD3nunrbpffjldLBRC1skk\nFZsQ0itZ//43LLRQuo7ij39023AVDYtoSS3zyitpRNVNN8GZZ8KTT8JKK2WdSlKxW3vtNPpu//3T\nhJ9NN4WxY7NOJc2RRbSkHxcjDB4MG2wAX3wBjzwCp5wCXbpknUxSqVhggfR75oYb0mzptdeGm2/O\nOpX0oyyiJc3eV1/BfvvBQQfBr34FI0fCZptlnUpSqdpjj/R7ZrXVYPfd4eCD0+8hqQBZREuatddf\nTyOorrsOTj89Xfyz9NJZp5JU6nr1Su1iJ5+cVqd793amtAqSRbSkH7r2WlhvvTR26qGH4LTTbN+Q\n1HG6dUv90UOHwuefw/rr88L++1NVWUlZWRlVVVXkcrmsU6qTs4iW9J2vv06tG/vtl4roESPSLFdJ\nykLfvvDyy7y32mqsP2QIp40bR48YaWhooKamxkJambKIlpS8+Wa6eHDw4PQy6sMPw7LLZp1KUmfX\nsycbf/oppwP7A08BFUBjYyO1tbWZRlPnZhEtKV0R37s3vP8+3H9/ehm1a9esU0kSAA3jx3MGsAOw\nMjAc2AIYN25cprnUuVlES53Z1Klw+OGw115px7CRI6F//6xTSdL3VFRUAHAPsB7wIfAQ8KdFFklj\nOKUMWERLnUgul6OqqoqysjLWXWEFPlpjDbj8cjjuuLTl7vLLZx1Rkn6grq6O8vJyAN4C+gB3dunC\nyZ99lkbhTZ6caT51ThbRUieRy+WoqamhoaGBNWPk9gkTWHD0aJ4+4gg477x0NbwkFaDq6mrq6+up\nrKwkhMASlZV8PWQInHsu3HYb9OkDb72VdUx1MiEWwcsgvXv3jsOGDcs6hlTUqqqqaGhoYBdgCPAJ\nsCPwcWUlY91iV1KxGjo0bdIyfXqaa7/99lknUpELIQyPMfae0/NciZY6ifENDZwO3AKMJPUVvoQX\n5kgqclttBcOHw4orwg47pM2hmpqyTqVOwCJa6gwmT+be+efnNOAq0lXtH+YfmnHBjiQVraoqeOYZ\n2GcfOOMMGDAAPvss61QqcRbRUql7913YaCP6ffMNx3XrxoHA1PxD5eXl1NXVZZlOktrG/PPDNdfA\nxRfDAw/A+uvD6NFZp1IJs4iWStnjj6edB8ePp+yBB1j76qv/d2FOZWUl9fX1VFdXZ51SktpGCHDk\nkfDoo/Dpp+mCwyeeyDqVSpQXFkql6oor4KijYKWV4K67YOWVs04kSR3nnXdgu+1gzBi48krYd9+s\nE6lIeGGh1FlNm5Y2UDnsMOjXD55/3gJaUufz05/Cc8/BxhvDfvvBqae6MYvalEW0VEomTUpXql9+\nOZxwQlqBXmSRrFNJUjYWWyz1R++/P5x1FlRXwzffZJ1KJaJr1gEktZExY9KW3ePGpVmp9jpLEnTv\nDoMHp1fkTj45/Y684w5Ycsmsk6nIuRItlYIXX4QNN4SPP04X1FhAS9J3QoCTToIbb4Rhw9IFh2++\nmXUqFTmLaKnY3XsvbLYZlJfDs8/CRhtlnUiSCtPuu8Njj8EXX6SFh8cfzzqRiphFtFTM/vGPtKnA\naqulC2hWXTXrRJJU2DbcMF1w/ZOfpIuvr70260QqUhbRUjGKMV1pXlOTLiR8/PH0H4Ikac5WXDG9\ncrfJJmlyxymnOLlDrWYRLRWbadO+u9L8wAPTBI4FF8w6lSQVl0UXhfvvT79H//Qn2GsvJ3eoVZzO\nIRWTL7+EXXaBhx6C009Pq9EhZJ1KkopTt26pLW7lleHEE+GDD+DOOx0NqhZxJVoqFh98kF56fOSR\nNK7ptNMsoCVpXoWQ5upffz088wxsvjl89FHWqVQELKKlYjBqVBrJ9NZbcM89cMABWSeSpNKy555w\n993wxhtpl8OxY7NOpAJnES0Vuqeegl/9CqZMgSeeSBuqSJLaXv/+6dW+SZPS793XXss6kQpYuxXR\nIYSrQggfhRBea3bs/BDCGyGEV0II/wohLNpeX18qCbfemqZvLLVUGmG37rpZJ5Kk0rbhhvDkk2la\nxyabpN+90iy050r0EGDmJbOhwOoxxjWB0cBJ7fj1peI2eDDstlsqnJ95Bnr1yjqRJHUOq6+efu8u\nsQRsuSWPHn88VVVVlJWVUVVVRS6XyzqhCkC7FdExxieBT2Y69lCM8dv8u88Dy7fX15eK2t//Dgcd\nlDYCGDo0/SKXJHWcXr3g6af5pGdPfn3++WzY0ECMkYaGBmpqaiyklWlP9AHA/Rl+fanwxJjmlf7+\n97DzzmnUUnl51qkkqXNaemk2bWriWSAHHJ4/3NjYSG1tbYbBVAgyKaJDCLXAt6Sfydk9pyaEMCyE\nMGzixIkdF07KSoxpTukpp8A++8BNN0GPHlmnkqRO7T8TJtAfuBu4FDg1f3zcuHHZhVJB6PAiOoQw\nENgeqI5x9ntsxhjrY4y9Y4y9e/bs2WH5pEw0NcERR8B558Ghh8KQIdDVvZAkKWsVFRV8A/yWdLHX\nGcBFQOUKK2QZSwWgQ4voEEJ/4HjgNzHGxo782lLB+vZbGDgQLr8cjjsOLrsMypw+KUmFoK6ujvLy\ncqaT+lD/AhwFPLbccjB1arbhlKn2HHF3A/AcsGoIYUII4UDgEmAhYGgIYWQI4Yr2+vpSIcvlclRV\nVTFfCNy/yCLwz3/CWWfBuee6C6EkFZDq6mrq6+uprKyEELikooIRu+9O1XPPwY47wtdfZx1RGWm3\n14tjjHvO4vDg9vp6UrHI5XLU1NQQGxu5A+jf2Mjx3bqxVq9eVFtAS1LBqa6uprq6+vsH+/aFQw6B\nAQPgjju8CLwTCj/SllwwevfuHYcNG5Z1DKlNVFVV8UlDA/cAGwMHA1cBlZWVjHWbWUkqHkOGwAEH\nwOabpy3DLaRLQghheIyx95yeZ+Ol1MEmNzTwCLAhsCepgAav9JakojNwIFx7LTz+OGy3HUyenHUi\ndSCLaKkj/fe/PN2tG2sAOwM3N3uooqIio1CSpLm2995w3XVpq/Btt4Uvv8w6kTqIRbTUUcaPh1//\nmpXKyvhtjx7c0+yh8vJy6urqMosmSZoHe+4JN9wAzz4L/fvDF19knUgdwCJa6ggTJqSeuY8+ouuj\nj7LX4MFUVlYSQqCyspL6+vofXrQiSSoeu+2WNsl64QXYemv4/POsE6mdeWGh1N7efx822wz++194\n6CHo0yfrRJKk9nLHHamgXnvt9Dt/0UWzTqRW8sJCqRB88EFagf7gA3jgAQtoSSp1O+4It90GL78M\nW24Jn3ySdSK1E4toqb18+CFssQW89x7cfz9stFHWiSRJHWGHHeBf/4JXX03zpD/+OOtEagcW0VJ7\n+OijVECPGwf33Qcbb5x1IklSR9p2W7jzThg1Kv1/MHFi1onUxiyipbY2cWJaeXj3XbjnHthkk6wT\nSZKy0L9/2oRl9OhUSH/0UdaJ1IYsoqW29PHHqQfu7bfTL87NN886kSQpS1ttBffeC++8k/5P+PDD\nrBOpjVhES23lk09SAf3mm+klvL59s04kSSoEW2yRro0ZOzYV1V5sWBIsoqW28Omn6Rfj66+n8Ub9\n+mWdSJJUSDbdNC2wvPmmG7KUCItoaV599lkqml99FW6/Pf1ylCRpZltuCbfeCiNGwPbbQ2Nj1ok0\nDyyipXnx+edpZ6qXX05zQbfbLutEkqRCtsMOcN118PTTsNNOMGVK1ok0lyyipbn15ZewzTbw0ktw\nyy3pF6MkSXOy++5w5ZVpR8Pdd4dp07JOpLlgES3NjcbGtOr8wgtw000wYEDWiSRJxeSAA+Cii1Kf\n9MCBMH161onUSl2zDiAVnWnTYLfd0ktx118PO++cdSJJUjE66ij46is46SQoL4f6eggh61RqIYto\nqTWamtKKwb33wuWXwx57ZJ1IklTMTjwRJk+GujpYcEH4298spIuERbTUUjHC0Uen1eezz4ZDD806\nkSSpFJx1ViqkL7wQFloIzjwz60RqAYtoqaVOOw0uvRSOPTatHEiS1BZCgAsuSK0dZ50FCywAJ5yQ\ndSrNgUW01BIXXph+sR1wAJx/vi+1SZLaVghwxRWpkD7xxNTaccQRWafSj7CIlubk2mvhD39IFxAO\nGmQBLUlqH126wDXXpAlQRx6ZVqQHDsw6lWbDEXfSj7nzzrT63Ldv6oXu6t+dkqR21K1bGp3arx9N\nBxzAET17UlZWRlVVFblcLut0asYiWpqdxx5LQ/DXXRf+9S/o0SPrRJKkzqBHD27cfXeeC4ELJk1i\n8xhpaGigpqbGQrqAhBhj1hnmqHfv3nHYsGFZx1BnMmwYbL45VFTAk0/CEktknUiS1IlUVVXxWUMD\nTwGVwCbAy0BlZSVjx47NNFupCyEMjzH2ntPzXImWZjZqFPTvD0sumbZktYCWJHWwcePG8TmwDfAZ\ncD9QlT+uwmARLTXX0AD9+qXe56FDYbnlsk4kSeqEKioqAHgP6A/0AB4E1vb/pYLRoiI6hLB8CGHz\n/P0eIYQF2jeWlIGPPoKttoIvv4QHH4SVVso6kSSpk6qrq6O8vByAUcAOwArAQz16pDF4ytwci+gQ\nwgHAXcCV+UOVwJ3tGUrqcF98kVo4JkxIW3qvtVbWiSRJnVh1dTX19fVUVlYSQuC9ykpe+P3vWfLd\nd2G33WDatKwjdnotWYk+GugDfAEQYxwNLNWeoaSOkMvlqKqqonsIPLn00jS9/DLcdhv86ldZR5Mk\nierqasaOHUtTUxNjx45l0wsugMsug/vug0MOgSIYDlHKWjL09psY49SQ32AihNAFcLcJFbVcLkdN\nTQ2NjY0MBjb55hsO696djT/5hOqsw0mSNDuHHALvvw9nngnLLgt/+lPWiTqtlqxEPxNCOB6YL98X\nfRNwT/vGktpXbW0tjY2NnAIcAJwBXDF1KrW1tRknkyRpDk4/HQ46COrq0sq0MjHHOdH5lecaoB9p\nBfpBYFCMsan94yXOiVZbKysrY58YuQa4BhiYPx5CoKmpw360JUmaO99+CzvvDPfcA7femu6rTbTJ\nnOh8AX11jPHyGONOMcYd8/etMlTU9uzZkyuBh4GDmx2fMVJIkqSC1rUr3Hgj9OkDe+2VNgZTh/rR\nIjrGOB1YMYTQrYPySO3vlVcY8sUXvBkCvwVmXN9cXl5OXV1dlskkSWq58nK4+27o1Qt+8xt47bWs\nE3UqLemJfgd4KoRwUgjh6Bm39g4mtYsJE2Dbbem2xBKMvvBCFsuPDqqsrKS+vp7qai8rlCQVkSWW\ngAceSAV1//4wfnzWiTqNlkznGJe/ledvUnH6/HPYdts0E/qpp9h5rbXY+Wj/HpQkFbnKylRI//rX\nqZB++mlYbLGsU5W8ORbRMcZTOiKI1K6mTYNddoFRo9xMRZJUetZcE+68E/r1S//fPfAAdLMbtz3N\nsYgOIQwFfjDCI8bYr10SSW0tRqipgYcfhquuSr9gJEkqNZttBv/4BwwcCEccAYMGQXBrj/bSknaO\nPza7Px/wW2BK+8SR2sGZZ8KQIXDaabD//lmnkSSp/ey3H4weDWefDauuCscem3WiktWSdo5/z3To\niRDCzMekwnT11Wko/cCBqYiWJKnUnXVWKqSPOw5WWgkGDMg6UUma43SOEMLCzW6LhhD6Anarq/A9\n9FBq49h4qDF0AAAgAElEQVRqK6iv9yUtSVLnUFYG11wDvXunGdIvvZR1opLUknaO/5B6ogPwLfAu\n39+fQio8r7ySLqz4+c/TTk5eXCFJ6kzKy+Guu2D99WGHHeCFF2C55bJOVVJaMid6xRhjRYxxhRhj\nrxjjFsAzc/qgEMJVIYSPQgivNTu2eAhhaAjhrfxbV7TV9j78ELbfHhZeOE3iWHjhrBNJktTxfvKT\ntC34F1+kQvqrr7JOVFJaUkTPqv/5hRZ83BCg/0zHTgQeiTGuDDySf19qO1OmwE47waRJadTP8stn\nnUiSpOysuWbaHvzll2HvvaGpKetEJWO2RXQIYakQwlrA/CGENUIIa+ZvG9OCTVdijE8Cn8x0eABw\nTf7+NcCOc5lb+qEZo+yeey71gq27btaJJEnK3nbbwQUXwB13wImuX7aVH+uJ3g44AFgeuKzZ8S+B\nud2AZekY4wf5+/8Flp7LzyP90F//Ctdem6Zx7Lpr1mkkSSocRx0Fb74J558Pq6wCBx2UdaKiN9si\nOsZ4NXB1CGG3GOPNbf2FY4wxhPCDTVxmCCHUADUAFRUVbf3lVWruvReOPz4Vz6e4yaYkSd8TAvz9\n7/DOO3DYYdCrF/Ttm3WqohZinG0d+92TQtga+AVpsxUAYoxnt+DjqoB7Yoyr599/E9gsxvhBCGEZ\n4PEY46pz+jy9e/eOw4YNm2NOdVL/+Q9suCGsvDI89VS6IlmSJP3Q55/DRhvB+++n9sef/SzrRAUn\nhDA8xth7Ts9ryZzoy4D9gGOA+YG9gZXmMtdd+c9F/u2dc/l5pGTSpHTF8QILpAsJLaAlSZq9RRZJ\nEzu6dUuTrCZNyjpR0WrJdI6NY4x7AR/HGE8BNqAFRXQI4QbgOWDVEMKEEMKBwJ+BrUIIbwFb5t+X\n5s7UqWkW9Pvvp4slnMQhSdKc9eqVFp4mTICdd06TrdRqLdls5ZsZb0MIPwE+Bpad0wfFGPeczUM2\n4GjexQhHHglPPAHXXQcbbJB1IkmSiseGG8KQIbDnnnDooXDVVe7s20otKaLvCyEsCvwFGAlM57sx\ndVI2LrkE/vEPOOkkqK7OOo0kScVnjz1g1Cg488w0FvbII7NOVFR+9MLCEEIZsF6M8d/59+cH5o8x\nzjz/uV15YaG+56GHYJttUi/07bdDWUu6kiRJ0g80NcGOO8J998HDD8Nmm2WdKHNtcmFhjLEJGNTs\n/a87uoCWvufNN2G33WD11VMbhwW0JElzr6ws/X+68sppTOy4cVknKhotqUAeCyEMaPck0px8+mla\nfe7eHe66CxZcMOtEkiQVv4UXThfoT52aVqUbG7NOVBRaUkQPBP4VQvg6hPBJCOHTEIKr0epY336b\nVqDHjk0tHJWVWSeSJKl0rLoqXH89jBwJNTXpAn79qJYU0UsC3YAFgZ7593u2ZyjpB445JvVqDRoE\nG2+cdRpJkkrPdtvBWWdBLsfwvfemqqqKsrIyqqqqyOVyWacrOHMsomOM04FdgRPy95cB1m7vYNL/\nDBkCF1+cCun99886jSRJpevkkxm33nqsff31rNTQQIyRhoYGampqLKRn0pIdCy8BNgf2yR9qBK5o\nz1DS/7z0UppfufnmcO65WaeRJKm0hUD///6X14GbgF75w42NjdTW1mYYrPC0pJ1joxjjIeQ3XclP\n5+jerqkkgI8/TjspLbUU3HQTdG3JWHNJkjQv3pgwgR1JReIdQHn++Dgnd3xPS4roafl50REghLAE\n0NSuqaTp09MuSh98ALfdBj1tw5ckqSNUVFQwBtgD+AVwdbPj+k5LiuhLgduAniGEM4CnAV9XV/s6\n9VQYOhQuvRTWWy/rNJIkdRp1dXWUl5fzEHAisBtwSrdu1NXVZZyssMzx9fEY47UhhOHAlvlDu8YY\nX2vfWOrU7rgDzj4bDjoo3SRJUoeprq4GoLa2lr82NPDr8nLO+PprwuKLZ5yssPzott//e1IIawIb\nk1o6nokxvtLewZpz2+9O5M0308rzz34GTz4J882XdSJJkjq3xkbYaKO0V8OLL6bdDUtYm2z7nf9E\ntcANwLLA8sD1IYST5j2ilORyOaqqqlg4BEavsUa6gvW22yygJUkqBOXl6VXirl1hwAD48susExWE\nlvRE7wusF2P8Y4yxFliftIuhNM9yuRw1NTU0NDQwGPjptGnsPHUquSefzDqaJEmaoaoqTcoaPRr2\n3ReanDHRkiL6A77fO901f0yaZ7W1tTQ2NnIsaUefE4H7p0xxFqUkSYWmb184//y0Kv2Xv2SdJnNz\n7IkOIdwOrAc8SOqJ7ge8CIwDiDEe084Z7YkuYWVlZWwWI0OB20lXAAOEEGjyr1xJkgpLjLD77nD7\n7fDII7DpplknanMt7Yluye4V9+ZvMzw/16mkmWyw7LLc9N57vAkc0Oy4syglSSpAIcCVV8LIkbDH\nHjBiBPzkJ1mnykRLRtwN7ogg6oS++Ya7e/SgO7ATMDl/uLy83FmUkiQVqoUXTgMANtgA9toLHnqo\nU+4q3JLpHP1DCC+GED4KIXwSQvg0hPBJR4RTiTv6aJYcM4YRv/89UyorCSFQWVlJfX39/2ZUSpKk\nArTGGnD55fDYY3DaaVmnyURLeqLfJrWqvkqz7b5jjNPbN9p37IkuQVdeCQcfDCefDK46S5JUnA4+\nOP2ffs89sN12WadpE202JxqYAIyMMU6LMU6fcZv3iOq0hg2DI46Afv3gzDOzTiNJkubWRRfB2mvD\nPvukzVg6kZY0sBwP3B1CeByYMuNgjPGi9gqlEvbZZ7DbbukihOuvhy5dsk4kSZLm1vzzwy23wLrr\nwq67wtNPQ48eWafqEC1ZiT4DmA4sCvRsdpNaJ0Y48EAYPz4NbF9iiawTSZKkebXSSjBkSHql+dhj\ns07TYVqyEr1CjHH1dk+i0nfJJWmu5F/+An36ZJ1GkiS1lZ12SgX0X/8Kv/oV7Lln1onaXUtWoh8M\nIWzR7klU2mb8dbrDDnBMu+/PI0mSOto556QC+uCDYdSorNO0u5YU0QcAD4cQJjviTnNlRh/0Msuk\nl3tCyDqRJElqa926pXbN8nLYZRf46qusE7WrlhTRSwLdgEVIvdBLYk+0WipGOOig1Ad9442w+OJZ\nJ5IkSe1lueXS4IBRo+CQQ1IdUKLmWETnx9ntCpyQv78MsHZ7B1OJuPTStKvROefAhhtmnUaSJLW3\nLbeEM86AXA7q67NO025asmPhJcDmwD75Q43AFe0ZSiVi+PDUB7399vZBS5LUmdTWwtZbw9FHp3qg\nBLWknWOjGOMhwDcAMcZPgO7tmkrF7/PPUx/00kunPuiylvyoSZKkklBWBtddl+qAXXaBTz/NOlGb\na0llMy2EUAZEgBDCEjTb/lv6gRl90A0NqQ/aedCSJHU+Sy6ZNmJ57720T0SJ9UfPtogOIcyYIX0p\ncBvQM4RwBvA0cG4HZFOxuuwyuPXW1Ae90UZZp5EkSVnZYAM4+2z4179g0KCs07SpEGfzV0EI4aUY\n4zr5+78AtgQC8HCM8bWOiwi9e/eOw4YN68gvqbn10kvpAsItt4S777aNQ5Kkzq6pCbbdFp54Al58\nEVYv7D38QgjDY4y95/i8HymiR8QYf9nmyeaCRXSR+PxzWHddmDIFRoxIL+NIkiR9+CGstVZq8Xzx\nxTRLukC1tIj+sW2/e4YQZjtSIcb4t7lKptIUY9qhaOzY9JemBbQkSZph6aXh2mvTxI5jjoErin/Q\n24+91t4FWBBYaDY36TuXX54uHqirS1t+SpIkNdevHxx/fOqNvu22rNPMsxb1RGfNdo4CN2IE9OkD\nffvCPffYBy1JkmZt6lTYeGN46y0YORIqK7NO9AMtbef4sWontGEelarJk9M86J4908s0FtCSJGl2\nundP42+nT4fqavj226wTzbUfq3j6dlgKFa+jj4Z33klbe9oHLUmS5mTFFVNLxzPPwJlnZp1mrs22\niM7vTCjN3i23wNVXw0knwaabZp1GkiQViz33hIED4U9/gscfzzrNXJltT3QhsSe6AI0bl0bVrLIK\nPP00dOuWdSJJklRMJk9Oo3EnT4aXXy6YV7TboidamrXp02GffVIfUy5nAS1JklpvwQVTf/SkSXDA\nAUW3LbhFtFrvvPPgySfh4othpZWyTiNJkorVL3+Z6oq774ZLLsk6TatYRKt1XnwRTj01TeTYb7+s\n00iSpGJ39NGw3Xbwf/+X2jqKRCZFdAjhDyGE/4QQXgsh3BBCmC+LHGqlyZNhr71gmWXSTkPBKYiS\nJGkehZAGFSyxBJ9vsw2rVVRQVlZGVVUVuVwu63Sz1eFFdAhhOeBooHeMcXXSzoh7dHQOzYWjj4Yx\nY+C662CxxbJOI0mSSkXPnjy8//4s9MEHHDt+PDFGGhoaqKmpKdhCOqt2jq7A/CGErkA58H5GOdRS\nzcfZbbJJ1mkkSVKJOSiX42zgIGC3/LHGxkZqa2szTDV7mYy4CyH8DqgDvgYeijFWz+I5NUANQEVF\nxboNDQ0dG1LfmTHObtVV4amnnMYhSZLaXFlZGV1i5AmgElgRmAqEEGhqauqwHAU74i6EsBgwAOgF\nLAssEELYe+bnxRjrY4y9Y4y9e/bs2dExNYPj7CRJUgeoqKjgW2BPYAtSAT3jeCHKop1jS+DdGOPE\nGOM04HZgowxyqCXOPTeNs7vkEvjpT7NOI0mSSlRdXR3l5eWMA0bnj5WXl1NXV5dlrNnKoogeB/QJ\nIZSHEALQFxiVQQ7NyQsvwGmnwe67w777Zp1GkiSVsOrqaurr66msrCSEQGVlJfX19VRX/6DrtyBk\n1RN9BrA78C0wAjgoxjhlds932+8MfPllGoA+bVqa2bjoolknkiRJanct7Ynu2hFhZhZjPA04LYuv\nrRY6+mh49114/HELaEmSpJm4Y6F+6OabYcgQOPlk+PWvs04jSZJUcCyiBUAul6OqqoqKEPh8zz2Z\n9NOfpu29JUmS9AMW0SKXy1FTU8O4hgYGA12amtjsvffI3Xxz1tEkSZIKkkW0qK2tpbGxkcOArYBj\ngf98803B7hAkSZKUNYtoMW7cOFYCzgfuB+qbHZckSdIPWUSLqhVW4BpgCmm/+hkKdYcgSZKkrFlE\ni9v69GEj4Ajg/fyxQt4hSJIkKWsW0Z3dK6/wy3/9i4b11+fZioqi2CFIkiQpa5lstqICMXVq2s57\nscWovOcexvbsmXUiSZKkomAR3ZmdeWba0vvOO8ECWpIkqcVs5+is/v1vOOccGDgQfvObrNNIkiQV\nFYvozqixMbVxLL88XHhh1mkkSZKKju0cndFJJ8Ho0fDII7DIIlmnkSRJKjquRHc2jz4KF10ERx0F\nW2yRdRpJkqSiZBHdmXz+Oey/P6y8Mvz5z1mnkSRJKlq2c3Qmf/gDTJgAzzwD5eVZp5EkSSparkR3\nFnffDVdfDSecAH36ZJ1GkiSpqFlEdwaTJsHBB8Oaa8Jpp2WdRpIkqejZzlHqYoTDD4dPPoGHHoIe\nPbJOJEmSVPQsokvdjTfCLbfA2WenlWhJkiTNM9s5StmHH8KRR8IGG8Bxx2WdRpIkqWRYRJeyo46C\nyZPTBYVdfdFBkiSprVhZlarbb09tHHV1sNpqWaeRJEkqKa5El6JPPkkXE669tm0ckiRJ7cCV6FJ0\nzDHw8cfwwAPQrVvWaSRJkkqOK9Gl5v774Zpr0qYqa6+ddRpJkqSSZBFdSr74Ag45JPVAn3JK1mkk\nSZJKlu0cpeSEE2DCBHj2WTdVkSRJakeuRJeKxx+HK66AP/wB+vTJOo0kSVJJs4guBY2NcNBB8NOf\nwllnZZ1GkiSp5NnOUQpOOQXeeQceewzKy7NOI0mSVPJciS52zz8PF1wAhx4Km22WdRpJkqROwSK6\nmE2ZAgccAMsvD+eem3UaSZKkTsN2jmJ21lkwalSaDb3wwlmnkSRJ6jRciS5WI0bAn/8M++0H/ftn\nnUaSJKlTsYguRtOmpTaOnj3hb3/LOo0kSVKnYztHMTr/fBg5Em6/HRZfPOs0kiRJnY4r0cXm9dfh\njDNgt91gp52yTiNJktQpWUQXk+nT4cADYaGF4OKLs04jSZLUadnOUUwuvzzNhb7uOlhqqazTSJIk\ndVquRBeL996Dk0+Gfv1gr72yTiNJktSpWUQXi6OPTlM5Lr8cQsg6jSRJUqdmO0cxuOuuNInjnHNg\nxRWzTiNJktTpuRJd6L78Eo44AtZYA449Nus0kiRJIqMiOoSwaAjh1hDCGyGEUSGEDbPIURROPTX1\nQw8aBN26ZZ1GkiRJZNfO8XfggRjjLiGE7kB5RjkK2/DhcNFFcOihsKF/Z0iSJBWKDi+iQwiLAJsA\nAwFijFOBqR2do+B9+y0cfHAaZXfOOVmnkSRJUjNZtHP0AiYCV4cQRoQQrgwhLJBBjoKUy+Woqqri\nmG7dYMQIntp1V1hkkaxjSZIkqZksiuiuwDrA5THGXwJfASfO/KQQQk0IYVgIYdjEiRM7OmMmcrkc\nNTU1NDU0cCZwD9D/yivJ5XJZR5MkSVIzIcbYsV8whJ8Az8cYq/Lv/xo4Mca43ew+pnfv3nHYsGEd\nlDA7VVVVNDQ0cBewBfBzYBxQWVnJ2LFjM80mSZLUGYQQhscYe8/peR2+Eh1j/C8wPoSwav5QX+D1\njs5RiMaNG8dvgR2AU0kF9IzjkiRJKhxZTec4CsjlJ3OMAfbPKEdB+cXyy3PR+PG8RBpfMkNFRUVW\nkSRJkjQLmRTRMcaRwByXyTub21ZdlaXHj+c3wPT8sfLycurq6rKMJUmSpJm4Y2GheO45VnnkEd7q\n149JlZWEEKisrKS+vp7q6uqs00mSJKmZrNo51Ny0aXDIIbDccvzs1lsZu9BCWSeSJEnSj7CILgR/\n+xu8+irccQdYQEuSJBU82zmyNmYMnHEG7LgjDBiQdRpJkiS1gEV0lmKEww+Hrl3h4ouzTiNJkqQW\nsp0jSzfdBA8+CBddBMsvn3UaSZIktZAr0Vn54gv4wx9g3XXTarQkSZKKhivRWTntNPjwQ7jrLujS\nJes0kiRJagVXorPwyiupB7qmBtZbL+s0kiRJaiWL6I4WIxxxBCy6KLgToSRJUlGynaOj/fOf8PTT\ncOWVsMQSWaeRJEnSXHAluiN99hkcdxz06QP77591GkmSJM0lV6I70imnwKRJcP/9UObfL5IkScXK\nSq6jjBgBl10Ghx0G66yTdRpJkiTNA4vojtDUlGZBL7EE/OlPWaeRJEnSPLKdoyMMGQLPP5/eLrpo\n1mkkSZI0j1yJbm+ffAInnAAbbwz77pt1GkmSJLUBi+j2VlsLn34Kl14KIWSdRpIkSW3AIro9vfgi\nDBoERx0Fa66ZdRpJkiS1EYvo9jJ9erqYcOml4fTTs04jSZKkNuSFhe1l8GAYNgxyOVhkkazTSJIk\nqQ25Et0eJk2Ck06CzTaDPffMOo0kSZLamEV0ezjpJPjiC7jkEi8mlCRJKkEW0W3t+efhyivh97+H\nX/wi6zSSJElqBxbRbWnGxYTLLgunnpp1GkmSJLUTLyxsS1dcASNGwE03wUILZZ1GkiRJ7cSV6LYy\ncSL88Y/Qty/sumvWaSRJktSOLKLbyimnwJdfwsUXezGhJElSibOIbgsjR0J9PRx5JKy2WtZpJEmS\n1M4soudVjHD00bDEEu5MKEmS1El4YeG8uuUWeOopGDQIFl006zSSJEnqAK5Ez4vGRvi//4O114YD\nD8w6jSRJkjqIK9Hz4rzzYPx4uO466NIl6zSSJEnqIK5Ez61x4+Dcc2G33WCTTbJOI0mSpA5kET23\njjsujbI7//ysk0iSJKmDWUTPjSeegJtvhhNOgIqKrNNIkiSpg1lEt9b06fC738EKK6TVaEmSJHU6\nXljYWldeCS+/DDfdBOXlWaeRJElSBlyJbo1PP4Xa2nQh4a67Zp1GkiRJGbGIbo0zzkiF9N//ni4q\nlCRJUqdkEd1Sr78Ol1wCBx+cNleRJElSp2UR3RIxwu9/DwsuCGedlXUaSZIkZcwLC1vi7rth6FC4\n8ELo2TPrNJIkScqYK9FzMmUKHHMMrLYaHH541mkkSZJUADJbiQ4hdAGGAe/FGLfPKsccXXABvPMO\nPPggdOuWdRpJkiQVgCxXon8HjMrw689WLpejqqqK5ULgq5NPZvw660C/flnHkiRJUoHIpIgOISwP\nbAdcmcXX/zG5XI6amhoaGho4G+gaI9u+/jq5XC7raJIkSSoQWa1EXwgcDzRl9PVnq7a2lsbGRtYH\n9gMuAF775htqa2szTiZJkqRC0eFFdAhhe+CjGOPwOTyvJoQwLIQwbOLEiR2UDsaNGwfAgsCzQN1M\nxyVJkqQsVqJ/BfwmhDAWuBHYIoRw3cxPijHWxxh7xxh79+zAsXIVFRUAPJoPOnmm45IkSVKHF9Ex\nxpNijMvHGKuAPYBHY4x7d3SO2amrq6O8vPx7x8rLy6mrq5vNR0iSJKmzcU70TKqrq6mvr6eyspIQ\nApWVldTX11NdXZ11NEmSJBWIEGPMOsMc9e7dOw4bNizrGJIkSSpxIYThMcbec3qeK9GSJElSK1lE\nS5IkSa1kES1JkiS1kkW0JEmS1EoW0ZIkSVIrWURLkiRJrWQRLUmSJLWSRbQkSZLUShbRkiRJUitZ\nREuSJEmtVBTbfocQJgINWefIwJLApKxDFDnP4bzx/M0bz9+88fzNG8/fvPH8zZtiPn+VMcaec3pS\nURTRnVUIYVhL9m7X7HkO543nb954/uaN52/eeP7mjedv3nSG82c7hyRJktRKFtGSJElSK1lEF7b6\nrAOUAM/hvPH8zRvP37zx/M0bz9+88fzNm5I/f/ZES5IkSa3kSrQkSZLUShbRGQshLB5CGBpCeCv/\ndrHZPG96CGFk/nZXs+O9Qgj/DiG8HUK4KYTQvePSZ68l5y+EsHYI4bkQwn9CCK+EEHZv9tiQEMK7\nzc7t2h37HWQjhNA/hPBm/ufmxFk83iP/8/R2/uerqtljJ+WPvxlC2LojcxeKFpy/Y0IIr+d/3h4J\nIVQ2e2yW/5Y7kxacv4EhhInNztNBzR7bL//v/a0Qwn4dm7wwtOD8XdDs3I0OIXzW7DF//kK4KoTw\nUQjhtdk8HkIIF+XP7yshhHWaPebP35zPX3X+vL0aQng2hLBWs8fG5o+PDCEM67jU7STG6C3DG3Ae\ncGL+/onAubN53uTZHL8Z2CN//wrgsKy/p0I7f8AqwMr5+8sCHwCL5t8fAuyS9ffRweesC/AOsCLQ\nHXgZ+PlMzzkcuCJ/fw/gpvz9n+ef3wPolf88XbL+ngrw/G0OlOfvHzbj/OXfn+W/5c5ya+H5Gwhc\nMouPXRwYk3+7WP7+Yll/T4V2/mZ6/lHAVc3e79Q/f/lzsAmwDvDabB7fFrgfCEAf4N/5453+56+F\n52+jGecF2GbG+cu/PxZYMuvvoa1urkRnbwBwTf7+NcCOLf3AEEIAtgBunZuPLxFzPH8xxtExxrfy\n998HPgLmOES9hK0PvB1jHBNjnArcSDqPzTU/r7cCffM/bwOAG2OMU2KM7wJv5z9fZzLH8xdjfCzG\n2Jh/93lg+Q7OWMha8vM3O1sDQ2OMn8QYPwWGAv3bKWehau352xO4oUOSFYkY45PAJz/ylAHAtTF5\nHlg0hLAM/vwBcz5/McZn8+cHSvz3n0V09paOMX6Qv/9fYOnZPG++EMKwEMLzIYQZheISwGcxxm/z\n708AlmvHrIWopecPgBDC+qTVm3eaHa7Lv/R0QQihRzvlLCTLAeObvT+rn5v/PSf/8/U56eetJR9b\n6lp7Dg4krWrNMKt/y51JS8/fb/P/Lm8NIazQyo8tZS0+B/k2ol7Ao80Od/afv5aY3Tn256/1Zv79\nF4GHQgjDQwg1GWVqM12zDtAZhBAeBn4yi4dqm78TY4whhNmNS6mMMb4XQlgReDSE8CqpsCl5bXT+\nyK8k/BPYL8bYlD98Eqn47k4ax3MCcGZb5JZCCHsDvYFNmx3+wb/lGOM7s/4MndbdwA0xxikhhENI\nr4pskXGmYrQHcGuMcXqzY/78qUOEEDYnFdEbNzu8cf7nbylgaAjhjfzKdlGyiO4AMcYtZ/dYCOHD\nEMIyMcYP8kXeR7P5HO/l3475//buLsSKMgzg+P9JscKyNLsQK3BDsBDT8EJKTPoQ7UKCDKQvSUGE\n7iKI2C6si+ouKCIiU4ggwlTaECpMxYKi79SyTBIiL7JMilRM6Oli3hPTwV3Pod095+j/B8POec68\nO/M+vLP77Ow7ZyJiJzAH2ET1b6ax5WrhFcChYe9Ahw1H/iJiArAV6C//nmt878ZV7JMRsQF4eBgP\nvVsdAq6svT7duGls81NEjAUuAY602PZs11IOIuJWqj/0bsrMk434IOfyuVTEnDF/mXmk9nId1b0P\njbYLm9ruHPYj7G7tnIPLgQfrAcdfSwbLseOvRRExi+rcXVI/n2vj73BEbKGantSzRbTTOTpvAGjc\n4bsCeLN5g4iY2JhmEBGTgRuBb7Kapb8DWDZU+7NcK/kbB2yhmuP2RtN7U8rXoJpPfdq7jc8ynwDT\no/pkl3FUv2ib79Kv53UZsL2MtwFgeVSf3jENmA58PErH3S3OmL+ImAO8CCzNzMO1+GnP5VE78u7Q\nSv6m1F4uBfaV9XeARSWPE4FFJXYuaeX8JSJmUN389mEt5vhrzQBwf/mUjnnA7+WCi+OvBRFxFbAZ\nuC8z99fi4yPi4sY6Vf56+3dup+9sPNcXqnmm7wHfA9uASSU+F1hX1m8A9lDdhb0HWFVr30dVxBwA\nNgLnd7pPXZi/e4FTwJe1ZXZ5b3vJ6V7gVeCiTvdplPJ2O7Cf6gpUf4k9QVX0AVxQxtOBMr76am37\nS7vvqK4ydLw/XZi/bcDPtfE2UOKDnsvn0tJC/p4Cvi552gHMqLVdWcblAeCBTvelG/NXXq8Fnm5q\n5/ir8vAa1ac0naKa17wKWAOsKe8H8HzJ7x5gbq2t4+/M+VsHHK39/Pu0xPvK2PuqnN/9ne7L/118\nYiQ78U4AAAImSURBVKEkSZLUJqdzSJIkSW2yiJYkSZLaZBEtSZIktckiWpIkSWqTRbQkSZLUJoto\nSeqw8nm0H0TEklrsroh4uwPHsjMi5o72fiWp1/jEQknqsMzMiFgDbIyIHVQ/m58EFo/kfmtPO5Uk\ntckiWpK6QGbujYi3gEeA8VRP2PzP45gjYjFVcT0G+DUzb4mIScB6qgcZHAdWZ+buIeJrgatL/MeI\nWAlsAK4DvgUuLPsaA7xM9eCiBNZn5jMjmQNJ6iUW0ZLUPR4HPgf+oipe/xURlwMvAQsy82Apkhtt\nvsjMOyLiZuAVYPYQcYBrgfmZeSIiHgKOZ+Y1ETGr7J+y7dTMnFn2f+kI9VmSepJFtCR1icw8FhGv\nA39m5smmt+cBuzLzYNn2txKfD9xZYtsj4rKImDBEHKrHkJ8o6wuAZ8t2uyNid4n/APRFxHPAVuDd\n4e6vJPUybyyUpO7yd1lG0rEzbZCZR6mmeOwE1gDrRviYJKmnWERLUm/4CFgQEdMAatM53gfuKbGF\nVHOl/xgi3mwXcHfZbiYwq6xPBs7LzE3AY8D1I9IrSepRTueQpB6Qmb9ExGpgc0ScBxwGbgPWAuvL\nNIzjwIrSZLB4sxeADRGxD9gHfFbiU0u8cbHl0eHtkST1tsjMTh+DJEmS1FOcziFJkiS1ySJakiRJ\napNFtCRJktQmi2hJkiSpTRbRkiRJUpssoiVJkqQ2WURLkiRJbbKIliRJktr0D6+gpdmirKsiAAAA\nAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f448e255890>"
      ]
     },
     "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"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.13"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}