{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "Model Setup\n",
    "-----------\n",
    "\n",
    "2D, Heat Equation with Dirichlet BC at top and bottom boundary.  No variation in the x direction.\n",
    "\n",
    "\\\\[\n",
    "\\frac{\\partial T}{\\partial t} = -v \\cdot \\frac{\\partial T}{\\partial y} +k \\frac{\\partial^2 T}{\\partial y^2}+H\n",
    "\\\\]\n",
    "\n",
    "\\\\[\n",
    "\\frac{\\partial T}{\\partial x} = 0\n",
    "\\\\]\n",
    "\n",
    "with $0 \\leqslant x \\leqslant 1 $ and $ y_{0}\\leqslant y \\leqslant y_{1}$\n",
    "\n",
    "$T(y_{1}) = T_{1}$\n",
    "\n",
    "$ -k\\nabla{T}\\mid_{y_{0}} = \\left[\\,0.0,\\,f\\,\\right] $\n",
    "\n",
    "------\n",
    "\n",
    "Effectively a 1D problem in $y$-axis, described by the analytic function\n",
    "\n",
    "$c_{0} = ( -\\frac{f}{k}-\\frac{h}{v} ) \\cdot \\frac{k}{v} \\cdot \\exp \\left[-\\frac{v}{k}y_{0} \\right]$\n",
    "\n",
    "$c_{1} = T_{1}-c_{0}\\exp \\left[\\frac{v}{k}y_{1} \\right] - \\frac{h}{v}y_{1}$\n",
    "\n",
    "\n",
    "$T = c_{0} \\exp \\left[\\frac{v}{k}y\\right] + \\frac{h}{v}y + c_{1}$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We implement the above boundary conditions using: \n",
    " * a `DirichletCondition` for $T(y_{1})=T_{1}$\n",
    " * a `NeumannCondition` object for $ k\\nabla{T}\\mid_{y_{0}} = \\left[\\,0.0,\\,f\\,\\right] $\n",
    "\n",
    "When the `NeumannCondition` object is associated with the `AdvectionDiffusion` object it defines a $T$ flux along a boundary such that:\n",
    "  * $ -k\\nabla T \\cdot n = h $ on $ \\Gamma_{h} $\n",
    "  \n",
    "    where:  \n",
    "    * $\\Gamma_{h}$ is the set of degree of freedom vertices IndexSets along the surface of the domain,  \n",
    "    * $ n $ is the unit normal facing outwards from the surface (at $n\\mid_{y_{0}}=\\left[0,-1\\right] $)  \n",
    "    * $ h $ is the scalar flow associated with the flux vector $ -k \\nabla T $ along $\\Gamma_{\\phi}$.\n",
    "\n",
    "An example implementation  \n",
    "\n",
    "```\n",
    "nbc = uw.conditions.NeumannCondition( fn_flux=[h], variable=tField,\n",
    "                                      indexSetsPerDof=mesh.specialSets[\"MinJ_VertexSet\"] )\n",
    "```\n",
    "\n",
    "Applies a scalar flow $h$ to the `tField` over the boundary vertices in the set `indexSetsPerDof`. The outward facing normal of this boundary set is used to calculate the $h$.\n",
    "\n",
    "Here $h$ can be an `underworld.Function` or `underworld.MeshVariable` type."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import underworld as uw\n",
    "import glucifer\n",
    "import numpy as np\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 = -0.60\n",
    "y1 = 1.3\n",
    "\n",
    "T1 = 8.0       # surface temperature\n",
    "k = 6.70       # diffusivity\n",
    "h = 8.0        # heat production, source term\n",
    "f = 2.0        # heat flux vector\n",
    "v = 2.47       # j-axis velocity component"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "# analytic solution definitions with velocity\n",
    "c0 = (-f/k-h/v) * k/v * np.exp(-v/k*y0)\n",
    "c1 = T1 - c0*np.exp(v/k*y1) - h/v*y1\n",
    "\n",
    "def analyticT(y):\n",
    "    return c0*np.exp(v/k*y) + h/v*y + c1\n",
    "\n",
    "def analytic_dT_dy(y):\n",
    "    return v/k*c0*np.exp(v/k*y) + h/v"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "# # without velocity\n",
    "# v = 0;\n",
    "# c0 = (-f+h*y0)/k\n",
    "# c1 = T1 + h/(2.0*k)*y1**2 - c0*y1\n",
    "\n",
    "# # analytic solution definitions\n",
    "# def analyticT(y):\n",
    "#      return -h/(2*k)*y*y + c0*y + c1\n",
    "\n",
    "# def analytic_dT_dy(y):\n",
    "#     return -h/k*y + c0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "# build mesh and fields\n",
    "mesh = uw.mesh.FeMesh_Cartesian( elementType = (\"Q1\"), \n",
    "                                 elementRes  = (10, 10), \n",
    "                                 minCoord    = (0., y0), \n",
    "                                 maxCoord    = (1., y1))\n",
    "\n",
    "tField    = mesh.add_variable( nodeDofCount=1, dataType='double')\n",
    "tDotField = mesh.add_variable( nodeDofCount=1, dataType='double')\n",
    "vField    = mesh.add_variable( nodeDofCount=2, dataType='double')\n",
    "\n",
    "# set entire tField to T1\n",
    "tField.data[:] = T1\n",
    "\n",
    "# set constant velocity field\n",
    "vField.data[:] = (0.0,v)\n",
    "\n",
    "# define neumann condition - flux!\n",
    "nbc = uw.conditions.NeumannCondition( fn_flux=-f, \n",
    "                                      variable=tField, \n",
    "                                      indexSetsPerDof=(mesh.specialSets['MinJ_VertexSet']) )\n",
    "\n",
    "# flag top boundary nodes with dirichlet conditions\n",
    "bc = uw.conditions.DirichletCondition(tField, indexSetsPerDof=(mesh.specialSets['MaxJ_VertexSet']) )\n",
    "\n",
    "# define heat eq. system\n",
    "ss = uw.systems.AdvectionDiffusion( phiField       = tField,\n",
    "                                    phiDotField    = tDotField,\n",
    "                                    velocityField  = vField,\n",
    "                                    fn_diffusivity = k,\n",
    "                                    fn_sourceTerm  = h,\n",
    "                                    conditions     = [bc, nbc] ) "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Number of iterations 1634\n"
     ]
    }
   ],
   "source": [
    "# evolve to a <1e-5 relative variation. Assume that's steady state solution.\n",
    "er = 1.0\n",
    "its = 0              # iteration count\n",
    "tOld = tField.copy() # old temperature field\n",
    "\n",
    "while er > 1e-5 and its < 2000:\n",
    "    \n",
    "    tOld.data[:] = tField.data[:] # record old values\n",
    "    \n",
    "    dt = ss.get_max_dt()          # get time steps\n",
    "    ss.integrate(dt)              # integrate in time (solve)\n",
    "    \n",
    "    absErr = uw.utils._nps_2norm(tOld.data-tField.data)\n",
    "    magT   = uw.utils._nps_2norm(tOld.data)\n",
    "    er = absErr/magT              # calculate relative variation\n",
    "    \n",
    "    its += 1\n",
    "\n",
    "print \"Number of iterations\",its\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAtoAAAF3CAYAAACbhOyeAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvNQv5yAAAIABJREFUeJzs3XecXFX9//HXZxPaRkJLgFCyoXdBWEMPXSBIVRBZKYoEUEH8ibQISlks+BVRaQuCIKEJIr23gEBkQ02khZJAABNKaKEEcn5/nFmzCZtkw+7snZ19PR+P+5iZO3dmPnMfs/j25NzPiZQSkiRJkjpXTdEFSJIkSdXIoC1JkiSVgUFbkiRJKgODtiRJklQGBm1JkiSpDAzakiRJUhkYtCVJkqQyMGhLkiRJZWDQliRJksrAoC1JkiSVQe+iC+hM/fr1S4MGDSq6DEmSJFWx0aNHv5FS6j+346oqaA8aNIjm5uaiy5AkSVIVi4jx7TnOqSOSJElSGRi0JUmSpDIwaEuSJEllYNCWJEmSysCgLUmSJJWBQVuSJEkqA4O2JEmSVAYGbUmSJKkMDNqSJElSGRi0JUmSpDIwaEuSJEll0LvoArq9hx6Cjz6CJZeE/v1h8cWhV6+iq5IkSVLBDNoddeyxcM89Mx7X1MASS+TQ3RK+W25b7i+5JCy9NAwYAAsvXFjpkiRJKh+Ddkc1NcHLL8PkyTBp0udvn3gi3779dtuv79MnB+6W4D3r1rK/Xz+I6NrvJkmSpC/MoN0BI0aMYPjw4UyYMIGBAwfS2NhIw2GHtX3wtGnwxhs5gP/3v/D66/Daa3lruf/443DLLfDee59//fzzw3LLzbwtv/zMj5dcMo+oS5IkqXBlC9oRcQHwdWBSSmnt0r49gV8CawCDU0rNs3ntS8B7wGfApyml+nLV+UWNGDGCYcOGMXXqVADGjx/PsGHDAGhoaPj8C+abb8Yo9dx88MHMAfy112DiRHjllbw9+GB+/Mknn/+MZZedEcLr6mDQoBnbwIGw0EId+t6SJElqn0gpleeNI4YA7wMXtwraawDTgXOBI+cStOtTSm/My2fW19en5uY237LTDRo0iPHjx39uf11dHS+99FL5C5g+PY+Qt4Tvl1+ecf+VV2DChLxv2rSZX7fUUjOH71nDuEFckiRpjiJidHsGgss2op1SGhkRg2bZ9xRAVMFc4wkTJszT/k5XUzPjwsr112/7mM8+y6PhL70E48fn25btkUfgmms+Pyq+7LKw0kp5W3nlGfdXWgkWW6y830mSJKmKVOoc7QTcFhEJODel1FR0QbMaOHBgmyPaAwcOLKCa2ejVa8b87c02+/zz06fn6Skt4fuFF2DcOHj+ebj55vxca4svPnPwXmUVWHVVWG21/JwkSZL+p1KD9mYppYkRsSRwe0Q8nVIa2daBETEMGAZdG3IbGxtnmqMNUFtbS2NjY5fV0GE1NbDMMnnbZJPPP//BBzOH7+efz/dHjYIrr8xBvUW/fjlwz7qtuGK+kFOSJKmHqcignVKaWLqdFBHXAIOBNoN2abS7CfIc7a6qseWCx891HWnrQsjuqk8fWGedvM3qk0/gxRfh2WfhmWdmbDfcABdcMOO4Xr1y2G4J3muuOWPr27frvoskSVIXq7igHRF9gJqU0nul+18DTiq4rDY1NDRUV7BuhzZbGh555MwHTZmSQ/esIfyOO/Iqmi2WWw7WWiuH7rXWytsaa8Aii3Ttl5IkSSqDcnYduQzYEugH/Bf4BfAW8CegPzAFeCyltH1ELAOcn1IaGhErAteU3qY3cGlKqV3zMbqy60hPNGtLQ8jTZZqamtr3fzg++yzPBR87Fv7znxm3Tz0FH34447jllpsRvtdeG7785XzfjiiSJKkCtLfrSNmCdhEM2uVVtpaGLQG8dfgeO3bmAF5Tky++XHfdHLxbtoEDXTFTkiR1KYO2Ol1NTQ1t/V4igumtL4zsLNOn54sxn3gib48/nm9feGHGMYsskueQzxrAa2s7vx5JkiQqoI+2qk+XtzSsqcm9vFdeGfbYY8b+996DMWNmBO8nnoCLL56xdH1NTZ7rvf76M7b11vPiS0mS1KUM2mq3imlpuPDCsPHGeWuRUp5+8thj8OijeUGeO+6Av/1txjGrrDJz+P7KV2CJJbq2dkmS1GMYtNVuFd3SMAJWWCFvu+8+Y//rr88I3o88knuAX3HFjOfr6qC+HgYPztsGG+QgL0mS1EHO0VbP89ZbM8L36NHQ3JwX44Ec2NdYI4fur3413375yy66I0mS/seLIaV58cYbOXD/+98ztsmT83Pzz5+nmbQE78GD89LzdjuRJKlHMmhLHZESTJgwc/AePTovSw95bvdGG+Wl6zfZJIfwPn2KrVmSJHUJu45IHRGR52/X1THik08Y/ve/88oHH7D1gAGcNHQoG6UEDzwAN96Yj+/VK7cY3GSTfJHmJpvk1zvqLUlSj+WItjQHc10N86234KGHcuh+8MF8sWXLqPeAATNC92ab5U4n881X0DeRJEmdxakjUieY59UwP/0Unnwyh+4HHsjbiy/m52prc/AeMgQ23xw23NCFdSRJ6oYM2lIn6JTVMF97De6/H+67D0aOzAvspJRHt+vrZwTvTTeFRRft5G8gSZI6m0Fb6gTzPKLdHlOmwL/+lUP3fffBww/nkfCI3EpwyJC8bbkl9OvXofolSVLnM2hLnWCuc7Q7w9SpeW53S/B+4AH48MP83LrrwlZbwdZb5/C9yCKd85mSJOkLM2hLnWTEiBFduxrmJ5/knt533w133ZWD90cfQU1NXrmyJXhvtpktBSVJKoBBW6oWH32UO5u0BO9Ro2DaNOjdO19Q2RK8N9kEFlig6GolSap6Bm2pWn3wQZ7j3RK8m5th+vTcwWTIEPja1/K25pr28ZYkqQwM2lJP8c47eX737bfn7emn8/5lloHttsvbttvCUksVW6ckSVXCoC31VBMmzAjdt9+eF9WBfGFly2j3ZpvBggsWW6ckSd2UQVsSfPYZPPpoDty33ZannEyblkP2FlvA0KGw446wyipFVypJUrdh0Jb0ee+/n6eZ3HYb3HILPPNM3r/yyjlwDx2aA/hCCxVbpyRJFcygLWnuXngBbr4ZbropX1j50Uc5ZG+11YzR7hVXLLpKSZIqikFb0rz58EO4994cum+6CZ5/Pu9fbbUcuocOzV1N5p+/2DolSSqYQVvSFzZixAjOO+oo1n31VfZYcEE2++wzek2bBgsvnEe5d945B+/FFy+6VEmSulx7g3bvrihGUvfRetn5e4E/fvQRSyy0EFcddhhbvvceXH89XHkl9OqVu5fssksO3l5QKUnSTBzRljSTQYMGMX78+M/tr6ur46WXXsqL4zQ3w3XX5e3JJ/MBq6+eQ/cuu8BGG+UgLklSFXLqiKQvpKamhrb+uxARTJ8+/fMveOmlPMp93XVwzz3w6afQrx/stBPstlvu211bW/a6JUnqKu0N2jVdUYyk7mPgwIHztJ9Bg+Cww3Kv7jfegMsvz+H62mth992hf3/45jfh0kvzKpaSJPUQBm1JM2lsbKR2lhHo2tpaGhsb5/7iRRaBb30LRoyASZNy+N5//7xQTkNDDt1Dh8L558PkyWX6BpIkVQaDtqSZNDQ00NTURF1dHRFBXV0dTU1NNDQ0zNsbzTcfbLstnHUWTJyYw/bhh8PTT8NBB8HSS+d+3X/6E7z8cnm+jCRJBXKOtqSulRI8/jhccw384x8wZkzeP3gw7LEH7Lmni+RIkiqaF0NK6h6efTaH7quvhocfzvvq6/MUlD33hLq6YuuTJGkWXgwpqXtYdVU4+mj4979zB5PTToMI+NnP8oWWG20Ev/+900skSd2OQVtS5airgyOPzKH7+efh17+GadPgpz+FgQNh003hjDPynG9Jkipc2YJ2RFwQEZMiYkyrfXtGxNiImB4Rsx1uj4gdIuKZiBgXEceUq0ZJFWzFFfNI9+jR8Nxz0NgIH3wARxwByy0Hm2+eL6T873+LrlSSpDaVc0T7r8AOs+wbA+wBjJzdiyKiF3AmsCOwJvDtiFizTDVK6g5WXhmOOw4eeyx3LTnpJJgyJXcxWWYZ2H57uPhiePfdoiuVJOl/yha0U0ojgbdm2fdUSumZubx0MDAupfRCSukT4HJg1zKVKam7WW01OP74vPT7mDFw7LF5xHv//WGppfJFlNdeC598UnSlkqQerhLnaC8LtL7q6ZXSPkma2VprwSmn5Pnc//oXHHgg3HVXXvp96aVh2DC4915oa+l4SZLKrBKD9jyJiGER0RwRzZNdaU7qmSJgk03gz3+GV1+Fm27KK1BeeilsuWW+yPKoo3L/7ipqaSpJqmyVGLQnAsu3erxcaV+bUkpNKaX6lFJ9//79y16cpAo333yw445wySX5QslLL4X11oPTT8+366yTWwi+9lrRlUqSqlwlBu2HgVUiYoWImB/YG7iu4JokdUd9+sC3vw3XX5+D9VlnQd++eXR7ueVyIL/8cvjww6IrlSRVoXK297sMeBBYLSJeiYgDI2L3iHgF2Bi4MSJuLR27TETcBJBS+hT4EXAr8BRwZUppbLnqlNRD9OsHhx4KDzwAzzyTL6IcOzYH8QED8nzuf/3LqSWSpE7jEuySeq7p0+Gee+Cii/IS8B98kFsJ7rcf7LtvXplSkqRZuAS7JM1NTQ1svTVcdBFXnHEGP11iCe4aNw5OOAFWWAG22gr++ld4//2iK5UkdUMGbUk93ogRI/je4Yfz+zffZBugDjhxvvl496mn4LvfnTG1ZNQop5ZIktrNoC2pxxs+fDhTp0793+MJwC+nTePLCywA998Pe+4JI0bARhvlriV/+AO88UZxBUuSugWDtqQeb8KECW3vf/ll2HRTuOCC3LXk3HNzJ5Of/ASWXTavQnnbbS6II0lqk0FbUo83cODAue/v23fG9JEnnsgdTO64A7bfPs/nPvFEmE1glyT1TAZtST1eY2MjtbW1M+2rra2lsbGx7Re0TB959VW44gpYbbUctAcNysH76qth2rTyFy5JqmgGbUk9XkNDA01NTdTV1RER1NXV0dTURENDw5xfuMACsNdeefrICy/A8cfDU0/BN78JAwfmx45yS1KPZR9tSepMn30Gt9wC55wDN94IETB0KBx8cF6JslevoiuUJHWQfbQlqQi9esFOO+Vl3198EY47DpqbYeedYcUV4ZRT8oWVkqSqZ9CWpHKpq4OTT87TR666ClZdNU8nGTgwTy+54w47lkhSFTNoS1K5zTcffOMbcPvt8OyzcMQReen37bbLF1L+7nfw1ltFVylJ6mQGbUnqSqusAqedBq+8khfBGTAAfvaz3Jf7+9+Hxx8vukJJUicxaEtSERZcEPbZB0aOzOF6333h0kthvfVg883hyittEShJ3ZxBW5KK9uUvQ1MTTJyYp5FMnJhXnRw0KM/xfv31oiuUJH0BBm1JqhSLLQY//Sk891zuWrLOOnDCCfniye98Bx56CKqoJaskVTuDtiRVml694Otfz/24n346L/d+3XWw8cYweDBcdBF89FHRVUqS5sKgLUmVbLXV4Iwz8nSSM8+EDz6AAw7IrQNPPBEmTSq6QknSbBi0Jak7WHhh+MEPYOzYvOR7fT388pd5WsmBB8KYMUVXKEmahUFbkrqTiNx/+8Yb4amn4Lvfhcsuy/O5t9sObr7ZRXAkqUIYtCWpu1p9dTj7bHj5ZTj1VPjPf2DoUFhrLTj3XJg6tegKJalHM2hLUne3xBJw7LHw4otwySVQWwuHHJKnlQwfDq++WnSFktQjGbQlqVrMPz80NEBzM9x7b1745le/yv2499sPnnyy6AolqUcxaEtStYmAIUPgmmtyT+5DD4V//CMvjDN0KNx9t/24JakLGLQlqZqttFJuDzhhApxyCoweDVtvDRtuCFddBZ99VnSFklS1DNqS1BMsvnierz1+fL5QcsoU2HPP3Kf77LPhww+LrlCSqo5BW5J6kgUXhGHDuPT44zm4f39GPf88/OAHfLT00nDyyfDmm0VXKElVw6AtST3MiBEjOOiQQ2iaPJmNgCHA3e+/DyeckDuV/PjH8NJLBVcpSd2fQVuSepjhw4cztVWP7fuAodOns92AAbDXXnkqycorw7775t7ckqQvxKAtST3MhAkT2tx/5+uvw4UXwgsv5FHta67Ji9984xv5IkpJ0jwxaEtSDzNw4MA5719uOfi//8sXTh5/PNx1F9TXww47wMiRXVipJHVvBm1J6mEaGxupra2daV9tbS2NjY0zH7jEEnDSSTlw//rX8OijsMUWeSGcW26xF7ckzYVBW5J6mIaGBpqamqirqyMiqKuro6mpiYaGhrZf0LcvHH10XuL9j3/MwXvHHfMo99VXw/TpXfsFJKmbiFSmEYmIuAD4OjAppbR2ad/iwBXAIOAlYK+U0tttvPYzoGWt4AkppV3a85n19fWpubm548VLkmbvk0/gkkvyKPdzz8Eaa8Axx8C3vw3zzVd0dZJUdhExOqVUP7fjyjmi/Vdgh1n2HQPcmVJaBbiz9LgtH6aU1itt7QrZkqQuMv/88L3vwVNPwWWXQe/esP/+efGbv/wFpk0rukJJqghlC9oppZHAW7Ps3hW4qHT/ImC3cn2+JKnMevWCvfeGxx+H667Lc7q//31YdVU4/3wDt6Qer6vnaC+VUnqtdP91YKnZHLdgRDRHxEMRYRiXpEoWATvvDP/+N9xwA/TvDwcdlAP3eeflqSaS1AMVdjFkypPDZzdBvK4072Uf4A8RsdLs3icihpVCefPkyZPLUaokqT0iYKedYNQouPFGWHJJGDYsB+6mJgO3pB6nq4P2fyNiAEDpdlJbB6WUJpZuXwDuAb4yuzdMKTWllOpTSvX9+/fv/IolSfMmAoYOhYcegptugqWXhoMPhlVWgXPPNXBL6jG6OmhfB+xfur8/cO2sB0TEYhGxQOl+P2BTwDWAJam7ichtAB98MPfdXmYZOOSQHLjPOQc+/rjoCiWprMoWtCPiMuBBYLWIeCUiDgR+DWwXEc8B25YeExH1EXF+6aVrAM0R8ThwN/DrlJJBW5K6qwjYfnt44AG49VZYdlk49NAcuM87z4smJVWtsvXRLoJ9tCWpG0gJ7rgDTjghTy9ZaSX45S9zH+5evYquTpLmqhL6aEuS9HkRsN12eYT7hhtg4YVh331hnXXgqqtcaVJS1TBoS5KK0dKlZPRo+Pvf874998xLu994Yx75lqRuzKAtSSpWTQ1885vw5JNw8cXwzjvw9a/DppvCXXcVXZ0kfWEGbUlSZejVK08hefrp3Abw5Zdhm21g663zNBNJ6mYM2pKkyjLffHmhm+eegz/8AcaOzaPbO+0Ejz5adHWS1G4GbUlSZVpwQfjxj+GFF+BXv8r9uNdfH/bZJ++TpApn0JYkVbY+feCYY3K4PvZY+Oc/YfXV4bDDYFKbCwxLUkUwaEuSuodFF4VTT4Vx4+C734Wzz57Rg/u994quTpI+x6AtSepellkmXyw5dmxecfLEE3Pg/tOf4JNPiq5Okv7HoC1J6p5WWy0vcDNqFKy1Fhx+OKyxBlx6qYveSKoIBm1JUvc2eHDut33zzXmVyYYG2GADuPVWF72RVCiDtiSp+4uAHXaARx6BSy6BKVPy4223zfskqQAGbUlS9aipySPaTz8NZ5wBTzyRl3Tff/+8AI4kdSGDtiSp+iywQJ6zPW4cHHUUXHEFrLoq/PzndiiR1GUM2pKk6rXIIvDrX+cR7j32gMZGWHnl3LXk00+Lrk5SlTNoS5Kq36BBjBg6lF2WXpqRkybBIYcwZdAguOkmL5iUVDYGbUlS1RsxYgTDhg3j+tdfZwtgN2Dyq6/CTjvBdtvBY48VXaKkKmTQliRVveHDhzN16tT/Pb4WWCslfrHYYjlkr79+Xm1y4sTiipRUdQzakqSqN2HChM/tmwacPGVKvmDyyCPzQjerrJKXdG8VyiXpizJoS5Kq3sCBA2e/f9FF4be/zRdM7rxzXtJ9tdXgssucvy2pQwzakqSq19jYSG1t7Uz7amtraWxsnLFjhRVyG8CRI6F/f9hnH9h8c2hu7uJqJVULg7Ykqeo1NDTQ1NREXV0dEUFdXR1NTU00NDR8/uDNN4eHH4bzz4fnnstLvH/ve/D6611fuKRuLVI7/lksIpYDVkkp3R0RCwC9U0oflL26eVRfX5+aHXmQJHWWd9+FU06BP/whL4Lz85/DEUfk+5J6rIgYnVKqn9txcx3RjojvAdcB55d21ZEv2JYkqbr17Zvnb48dC1tvDcccA2utBf/8p/O3Jc1Ve6aOHA5sBLwLkFJ6FliynEVJklRRVlkFrr0Wbrstj2bvvnvuvz1mTNGVSapg7QnaH6WUPml5EBG9gChfSZIkVajttoPHH4c//QkeeQTWXRd++EN4++2iK5NUgdoTtP8VEUcBC0bEVsAVwA3lLUuSpArVuzf86Ef5QslDD4VzzoFVV4W//AWmTy+6OkkVpD1B+yjgPeBp4MfAncDwchYlSVLFW2IJ+POf88j26qvD978PG29sO0BJ/zPHoF2aJnJhSunslNLuKaXdSvf9v+ySJEGePjJyJFx8MYwfn9sBHnwwvPlm0ZVJKtgcg3ZK6TNgxYiYr4vqkSSp+4mAffeFZ56BH/84TyNZddU8reSzz4quTlJB2jN15Hngvog4NiIOb9nKXZgkSd3OIovA6afDY4/BOuvkOdwbbggPPVR0ZZIK0J6gPQG4HagF+rfaJElSW9ZeG+6+Gy69FF59Nc/dPvBAmDy56MokdaF2rQz5hd884gLg68CklNLapX2LkzuXDAJeAvZKKX2uL1JE7A/8vPTwlJTSRXP7PFeGlCRVnPfeg5NOyqtLfulLcPLJeaS7V6+iK5P0BXXmypC3R8Rts27trOOvwA6z7DsGuDOltAq5g8kxbXzm4sAvgA2BwcAvImKxdn6mJEmVY+GF4bTTcv/t9deHww7L00kcGJKqXnumjvwcOL60NZLb/D3enjdPKY0E3ppl965Ay+j0RcBubbx0e+D2lNJbpdHu2/l8YJckqftYc0244w64/HKYODF3JznsMHjnnaIrk1Qmcw3aKaVRrbZ7U0qHA0M68JlLpZReK91/HViqjWOWBV5u9fiV0j5JkrqvCPjWt+Dpp/OKkmeemXtwX3EFlHEqp6RitGfqSN9W26IRsQ3QKdM4Up4g3qH/skTEsIhojojmyV5kIknqDhZZJC/j/u9/wzLLwN57ww47wLhxRVcmqRO1Z+rIWGBM6fZR8qqQB3XgM/8bEQMASreT2jhmIrB8q8fLlfZ9TkqpKaVUn1Kq79/fZiiSpG6kvj6H7T/+ER58MHcrOflk+PjjoiuT1AnaE7RXTCkNTCktn1JaIaW0NfCvDnzmdcD+pfv7A9e2ccytwNciYrHSRZBfK+2TJKm69OqV52o//TTsuiuccEJebfLuu4uuTFIHtSdoj2pj37/b8+YRcRnwILBaRLwSEQcCvwa2i4jngG1Lj4mI+og4HyCl9BZwMvBwaTuptE+SpOq0zDJ5rvbNN8O0abD11rDffjCprX/4ldQdzLaPdkQsCQwALgf2AqL0VF/g/JTS6l1S4Tywj7YkqSp8+CE0NsJvfwt9+uTbAw+EmvaMj0kqt87oo70T8Gfy/OizgDNL23HkVn+SJKkcFloITjkl995ed10YNgy22gqeeaboyiTNg9kG7ZTShSmlzYEDU0qbt9qGppT+3oU1SpLUM62xRp6rff758MQTOXQ3NsInnxRdmaR2aNcS7BGxPbAWsGDLvpTSqWWs6wtx6ogkqWq9/jocfjj8/e+wzjpw3nl5hUlJXa4zl2A/i9wd5P8BCwHfAVbucIWSJKn9ll4arrwSrr0W3noLNt4YjjgC3n+/6MokzUZ7rqrYLKW0D/BmSul4YEMM2pIkFWOXXeA//4FDD839t9daC266qeiqJLWhPUH7o5bbiFi69HiZ8pUkSZLmqG/fvHz7/ffnriQ77QT77GMrQKnCtCdo3xQRiwK/Ax4DXgKuLGdRkiSpHTbZBB59FH75S7jqqnzx5EUXQTuuv5JUfnMM2hFRA9ycUppS6jSyArBOSum4LqlOkiTN2QILwC9+AY89BquvDgccAF/7Grz4YtGVST3eHIN2Smk6cG6rxx+6QqMkSRVozTXhvvvgrLNg1KjcmeTPf4bp04uuTOqx2jN15O6I2LXslUiSpI6pqckXSY4ZA5ttBocdlhe6GTeu6MqkHqk9QfsA4JqI+DAi3oqItyPCUW1JkirVwIFw8808eNBBvHvffUxdZRVOXnxxLv3b34quTOpR2hO0+wHzAV8C+pce9y9nUZIkqWNGXHop244YwRopcSdw/Ntvs+IBB3D9aacVXZrUY8w1aKeUPgP2BI4u3R8ArFfuwiRJ0hc3fPhwpk6dyqvALuTV5ladPp2vHXUU/Pa38OmnBVcoVb/2rAz5Z2ArYN/SrqnAOeUsSpIkdcyECRNmejwCWBO4EeDoo3NrwDFjCqhM6jnaM3Vkk5TSwZQWril1HZm/rFVJkqQOGThw4Of2/Rf4fwMHwuWX5/Z/668Pp5wC06Z1fYFSD9CeoD2t1E87AUTEEoC9giRJqmCNjY3U1tbOtK+2tpbGU0+Fb30rL+O+++5w/PGw4Ybw+OMFVSpVr/YE7TOBq4H+EXEicD/wm7JWJUmSOqShoYGmpibq6uqICOrq6mhqaqKhoSEf0L8/XHEFXH01vPoqfPWrcOqpzt2WOlGkdizTGhFrAduWHt6ZUqrISV319fWpubm56DIkSepe3ngDfvAD+Pvf8+j2RRfBaqsVXZVUsSJidEqpfm7HtWdEG6AXMA34ZB5eI0mSuoN+/eDKK/Pc7eeeg698Bf74R1eVlDqoPV1HhgOXAcsAywGXRsSx5S5MkiR1sW99K3ci2Xpr+PGPYdttYfz4oquSuq32jE7vB3w1pfTzlNJwYDB5tUhJklRtBgyA66+H88+Hhx+GddaBv/wF2jHVVNLM2hO0XwN6t3rcu7RPkiRVowg48EB48knYYAP4/vdh553hNf/nX5oX7QnabwFjI+L8iDgPeBJ4IyJ+HxG/L295kiSpMIMISCzrAAAdxklEQVQGwZ13whln5Nu11srzuCW1S++5H8KNpa3FQ2WqRZIkVZqaGjj8cNh+e9h/f/j2t+Gaa+DMM/NFlJJma65BO6X0l64oRJIkVbDVVoP774fTToNf/ALuvRcuvBB23LHoyqSK1Z6uIztExMMRMSki3oqItyPira4oTpIkVZDeveHYY/NFkv37w9Ch8KMfwdSpRVcmVaT2zNH+M3AwsCzQH+hXupUkST3RuuvmsP2Tn+QpJPX18MgjRVclVZz2BO1XgMdSStNSSp+1bOUuTJIkVbAFF4Tf/x5uvx3eeQc22gh+8xv4zIggtWjPxZBHAddHxD3Axy07U0p/LFdRkiSpm9h229wG8OCD4Zhj4Kab4OKLoa6u6MqkwrVnRPtE4DNgUfKUkZZNkiQJFl88L+F+0UXw6KPw5S/DJZe4yI16vPaMaC+fUlq77JVIkqTuKwL22w823xz23TdvN9wAZ58Niy1WdHVSIdozon1rRGxd9kokSVL3t8IKufVfYyNcfXUe3b7rrqKrkgrRnqD9PeCOiHjf9n6SJGmuevWC446DBx+E2lrYZhs48kj4+OO5v1aqIu0J2v2A+YBF6KT2fhHx44gYExFjI+KINp7fMiLeiYjHStsJHfk8SZJUgJa2f4ceCv/3f7kzydNPF12V1GXmGrRLrfz2BI4u3R8ArPdFPzAi1gYOAgYD6wJfj4iV2zj0vpTSeqXtpC/6eZIkqUB9+sBZZ8F118HLL8MGG8AFF3ihpHqE9qwM+WdgK2Df0q6pwDkd+Mw1gFEppakppU+Be4E9OvB+kiSp0u28Mzz+OGy4IRx4IHz72zBlStFVSWXVnqkjm6SUDgY+AkgpvQXM34HPHANsHhFLREQtMBRYvo3jNo6IxyPi5ohYqwOfJ0mSKsGyy+YFbk49Fa66CtZbDx54oOiqpLJpT9CeFhE1QAKIiCWA6V/0A1NKTwG/AW4DbgEeI/fpbu0RoC6ltC7wJ+Cfs3u/iBgWEc0R0Tx58uQvWpYkSeoKvXrBscfC/fdDTQ0MGZI7lLiipKrQbIN2RLT02D4TuBroHxEnAveTg/IXllL6S0ppg5TSEOBt4NlZnn83pfR+6f5NwHwR0W8279WUUqpPKdX37+86OpIkdQsbbZQXt9lrL/j5z/MKk6+8UnRVUqea04j2vwFSShcDPwd+Rw7Fe6aULu/Ih0bEkqXbgeT52ZfO8vzSERGl+4NLdb7Zkc+UJEkVZpFFYMQI+Otf4eGHYd114dpri65K6jRzWhkyWu6klMYCYzvxc68uTUGZBvwwpTQlIg4pfdY5wDeBQyPiU+BDYO+UvDxZkqSqEwH77w8bb5wvkNxtN/jBD+B3v4OFFiq6OqlDYnb5NSJeAX4/uxemlGb7XFHq6+tTc3Nz0WVIkqQv4uOPYfjw3HN77bXh8sthLfshqPJExOiUUv3cjpvT1JFewJeAhWezSZIkdZ4FFsgj2TffDJMmwVe/as9tdWtzmjrymgvFSJKkLrfDDrnn9ne+k3tu33NPXvTmS18qujJpnsxpRDvm8JwkSVL5LL003HornHgiXHJJHt1+8smiq5LmyZyC9jZdVoUkSdKsevWCE06AO+/Mq0gOHgznn+9UEnUbsw3apRUgJUmSirXVVvDYY7DZZnDQQXlKyXvvFV2VNFftWRlSkiSpWEstBbfcAiefnLuR1NfDE08UXZU0RwZtSZLUPfTqlVeRvOuuPKK94YbQ1ORUElUsg7YkSepettgiTyUZMgQOPhj22QfefbfoqqTPMWhLkqTuZ8klc7/txka48so8leSxx4quSpqJQVuSJHVPNTVw3HFw993wwQew0UZw3nlOJVHFMGhLkqTubciQPJq9xRYwbBhXLbwwfSIYNGgQI0aMKLo69WAGbUmS1P3178+l3/kOp/buzR4ffMADQO/x4xk2bJhhW4WJVEX/vFJfX5+am5uLLkOSJBVg0KBBjB8/nu2BEUBv4ADg0bo6XnrppSJLU5WJiNEppfq5HeeItiRJqgoTJkwA4FZgfeAZ4BrgR+PHw6efFliZeiqDtiRJqgoDBw783/0JwObAmcCRANtsA6+9Vkxh6rEM2pIkqSo0NjZSW1v7v8efAEfV1vKvQw+F5mb4ylfg3nuLK1A9jkFbkiRVhYaGBpqamqirqyMiqKuro6mpiU3POgtGjYJFFskj26edZgtAdQkvhpQkST3Du+/CgQfCVVfBbrvBhRfCoosWXZW6IS+GlCRJaq1v37yK5Omnww035NUkH3+86KpUxQzakiSp54iAI46Ae+6BDz+EjTeGSy8tuipVKYO2JEnqeTbdFB55JI9qNzTAT34C06YVXZWqjEFbkiT1TEstBXfeCYcfDn/4A2y3HUyaVHRVqiIGbUmS1HPNNx+ccQb87W+5M8kGG8DDDxddlaqEQVuSJOk734F//Qt69YLNN4cLLii6IlUBg7YkSRLA+uvnhW023zy3AfzBD+CTT4quSt2YQVuSJKlFv35w881w1FFw9tmw1Vbw6qtFV6VuyqAtSZLUWu/e8JvfwBVX5D7bG2yQp5VI88igLUmS1Ja99oKHHoI+ffLI9llnuXS75olBW5IkaXbWXjvP295uO/jhD/Pc7Y8+KroqdRMGbUmSpDlZdFG4/no44QS48ELYckt47bWiq1I3YNCWJEmam5oaOPFEuOoqGDMmryhpv23NhUFbkiSpvb7xDXjgAZh//twG8JJLiq5IFayQoB0RP46IMRExNiKOaOP5iIg/RsS4iHgiItYvok5JkqTP+fKX82j2xhvDvvvmVoCffVZ0VapAXR60I2Jt4CBgMLAu8PWIWHmWw3YEViltw4Czu7RISZKkOenXD267LS9qc9ppsPPOMGVK0VWpwhQxor0GMCqlNDWl9ClwL7DHLMfsClycsoeARSNiQFcXKkmSNFvzzQdnngnnngu33w4bbQTPPFN0VaogRQTtMcDmEbFERNQCQ4HlZzlmWeDlVo9fKe2TJEmqLMOGwV13wVtvwYYb5pUlJQoI2imlp4DfALcBtwCPAV94YlNEDIuI5ohonjx5cidVKUmSNA823zzP215hBfj61+F3v3NxGxVzMWRK6S8ppQ1SSkOAt4FnZzlkIjOPci9X2tfWezWllOpTSvX9+/cvT8GSJElzU1cH99+fO5P87Gew337w4YdFV6UCFdV1ZMnS7UDy/OxLZznkOmC/UveRjYB3Ukp2hpckSZWtTx+44go4+eTc+m+LLWBim2OF6gGK6qN9dUT8B7ge+GFKaUpEHBIRh5Sevwl4ARgHnAf8oKA6JUmS5k0E/Pzn8M9/wlNPweDBMHp00VWpAJGqaP5QfX19am5uLroMSZKk7Mknc+u/SZPyCPceszZaU3cUEaNTSvVzO86VISVJksplnXVg1ChYd908d/tXv/IiyR7EoC1JklROSy0Fd98N++wDxx0HBxwAH39cdFXqAr2LLkCSJKnqLbhgnjqyxhpw/PHw/PNwzTVgx7Sq5oi2JElSV2i5SPKKK/LFkYMHw9ixRVelMjJoS5IkdaW99oKRI+Gjj2DjjeGWW4quSGVi0JYkSepqX/0q/PvfsNJKsNNO8Kc/eZFkFTJoS5IkFWH55eG++2CXXeDww+GHP4Rp04quSp3IoC1JklSUL30Jrr4ajj4azj4bhg6FKVOKrkqdxKAtSZJUpJoa+PWv4cIL4d57YaON4IUXiq5KncCgLUmSVAkOOADuuAMmT85he9SooitSBxm0JUmSKsWQIfDAA7DwwrDllnlaibotg7YkSVIlWW01eOghWG892HNP+L//syNJN2XQliRJqjT9+8Ndd8E3vgFHHgk/+hF8+mnRVWkeGbQlSZIq0UIL5VUkjzoKzjoLdtsN3n+/6Ko0DwzakiRJlaqmBn7zGzjnnLyC5JAh8OqrRVeldjJoS5IkVbqDD4brr4fnnoMNN4Qnnyy6IrWDQVuSJKk72HFHuP/+fGHkppvCbbcVXZHmwqAtSZLUXay7bu5IssIKeRXJ888vuiLNgUFbkiSpO1luObjvPthuOzjoIMbssgsr1NVRU1PDoEGDGDFiRNEVqsSgLUmS1N307QvXX89zW2/N2tdfzykTJtA7JcaPH8+wYcMM2xXCoC1JktQd9e7NduPGcQzQANwE9AWmTp3K8OHDi61NgEFbkiSp25rw8sv8BtgX2AIYCQwAJkyYUGhdygzakiRJ3dTAgQMBuAQYCqwIPAhsNWBAgVWphUFbkiSpm2psbKS2thaAO4AhwALATe+8k1sBqlAGbUmSpG6qoaGBpqYm6urqiAjerqvjodNPZ4HlloNtt4V//KPoEnu0SCkVXUOnqa+vT83NzUWXIUmSVKw334Sdd849t//4R/jRj4quqKpExOiUUv3cjnNEW5IkqdossQTccQfssgscdhgcfTRMn150VT2OQVuSJKka1dbC1VfDoYfCb38L++0Hn3xSdFU9Su+iC5AkSVKZ9OoFZ56ZV5McPhxefz3P2+7bt+jKegRHtCVJkqpZBBx3HFx4Idx7LwwZAq++WnRVPYJBW5IkqSc44AC44QYYNw422QSefbboiqqeQVuSJKmn2H77PKo9dSpsuimMHl10RVXNoC1JktSTbLBBXsymTx/Ycku4666iK6pahQTtiPhJRIyNiDERcVlELDjL8wdExOSIeKy0fb+IOiVJkqrSqqvCAw/AoEGw445w1VVFV1SVujxoR8SywOFAfUppbaAXsHcbh16RUlqvtJ3fpUVKkiRVu2WWgZEj4atfhb32gnPOKbqiqlPU1JHewEIR0RuoBbz0VZIkqastthjcdhsMHZr7bZ98MlTRquFF6/KgnVKaCPwOmAC8BryTUrqtjUO/ERFPRMRVEbF8lxYpSZLUU9TWwjXX5AVtTjgBDj/cVSQ7SRFTRxYDdgVWAJYB+kTEd2Y57HpgUErpy8DtwEVzeL9hEdEcEc2TJ08uV9mSJEnVa775cp/tn/4U/vxnaGhwFclOUMTUkW2BF1NKk1NK04B/AJu0PiCl9GZK6ePSw/OBDWb3ZimlppRSfUqpvn///mUrWpIkqarV1MDvfpeXa7/8cth5Z3j//aKr6taKCNoTgI0iojYiAtgGeKr1ARExoNXDXWZ9XpIkSWXys5/BBRfAHXfANtvAG28UXVG3VcQc7VHAVcAjwJOlGpoi4qSI2KV02OGl9n+PkzuUHNDVdUqSJPVY3/1unrf9xBOw+eYwYULRFXVLkaroytL6+vrU3NxcdBmSJEnVYeRI2GUXWHhhuP12WH31oiuqCBExOqVUP7fjXBlSkiRJbRsyJC/ZPm1aHtl+5JGiK+pWDNqSJEmavXXXhfvuy20At9oqL9+udjFoS5Ikac5WWSUH7AED4Gtfg1tvLbqibsGgLUmSpLlbfvk8Z3u11XLrv6uvLrqiimfQliRJUvssuSTcfTcMHgx77ZUXudFsGbQlSZLUfosumqeObLstfO97cMYZRVdUsQzakiRJmjd9+sB118Eee8ARR8BJJ0EVtYzuLAZtSZIkzbsFFoArroADDoBf/AJ++lPD9ix6F12AJEmSuqneveEvf4G+feH00+Hdd+Hcc6FXr6IrqwgGbUmSJH1xNTXwhz/AIovAySfnsH3JJTD//EVXVjiDtiRJkjomIs/TXmQROPJIeO+93P6vtrboygrlHG1JkiR1jp/+FM47L3cl2XHHHLh7MIO2JEmSOs/3vw+XXgr/+ldeRXLKlKIrKoxBW5IkSZ1r773h73+H0aNhm23gzTeLrqgQBm1JkiR1vt13h3/+E8aOhS23hP/+t+iKupxBW5IkSeUxdCjceCO88AJssQVMnFh0RV3KoC1JkqTy2WYbuOUWePVVGDIExo8vuqIuY9CWJElSeW2+Odx+O7z1Vg7b48YVXVGXMGhLkiSp/DbcEO66Cz74IIftp54quqKyM2hLkiSpa3zlK3DvvTB9ep6z/cQTRVdUVgZtSZIkdZ211oKRI/MS7VttlVsAVimDtiRJkrrWqqvmsN23L2y9NTz4YNEVlYVBW5IkSV1vxRVz2F5ySdhuuzylpMoYtCVJklSM5ZfPYbuuDnbcEe68s+iKOpVBW5IkScUZMADuuQdWXhm+/vXcBrBKGLQlSZJUrP7982j2qqvCzjvDrbcWXVGnMGhLkiSpeC1he401YNdd82qS3ZxBW5IkSZWhXz+44w5Yc80ctm+6qeiKOsSgLUmSpMqxxBI5bK+9Nuy+O9xwQ9EVfWEGbUmSJFWWxRfPYXuddWCPPeD664uu6AsxaEuSJKnyLLZYDtvrrQff+AZce23RFc0zg7YkSZIq06KLwm23wfrrwze/CddcU3RF86SQoB0RP4mIsRExJiIui4gFZ3l+gYi4IiLGRcSoiBhURJ2SJEkq2KKL5nZ/9fWw115w9dVFV9RuXR60I2JZ4HCgPqW0NtAL2HuWww4E3k4prQycDvyma6uUJElSxVhkkRy2v/pV+Na3uO/wwxk0aBA1NTUMGjSIESNGFF1hm4qaOtIbWCgiegO1wKuzPL8rcFHp/lXANhERXVifJEmSKknfvnDrrUxaaSU2/tOfGDx+PCklxo8fz7BhwyoybHd50E4pTQR+B0wAXgPeSSndNsthywIvl47/FHgHWKIr65QkSVKFWXhhtvrwQx4ELgW+Vdo9depUhg8fXmBhbSti6shi5BHrFYBlgD4R8Z0OvN+wiGiOiObJkyd3VpmSJEmqQE+98go7AncD77faP2HChIIqmr0ipo5sC7yYUpqcUpoG/APYZJZjJgLLA5SmlywCvNnWm6WUmlJK9Sml+v79+5exbEmSJBVt4MCBfAB8Dbhxlv2VpoigPQHYKCJqS/OutwGemuWY64D9S/e/CdyVUkpdWKMkSZIqUGNjI7W1tTPtq62tpbGxsaCKZq+IOdqjyBc4PgI8WaqhKSJOiohdSof9BVgiIsYB/w84pqvrlCRJUuVpaGigqamJuro6IoK6ujqamppoaGgourTPiWoaKK6vr0/Nzc1FlyFJkqQqFhGjU0r1czvOlSElSZKkMjBoS5IkSWVg0JYkSZLKwKAtSZIklYFBW5IkSSoDg7YkSZJUBgZtSZIkqQwM2pIkSVIZGLQlSZKkMjBoS5IkSWVQVUuwR8RkYHzRdRSgH/BG0UV0Y56/jvH8dYznr2M8fx3nOewYz1/HdNfzV5dS6j+3g6oqaPdUEdGcUqovuo7uyvPXMZ6/jvH8dYznr+M8hx3j+euYaj9/Th2RJEmSysCgLUmSJJWBQbs6NBVdQDfn+esYz1/HeP46xvPXcZ7DjvH8dUxVnz/naEuSJEll4Ii2JEmSVAYG7W4iIhaPiNsj4rnS7WKzOe6ziHistF3Xav8KETEqIsZFxBURMX/XVV+89py/iFgvIh6MiLER8UREfKvVc3+NiBdbndv1uvYbFCMidoiIZ0q/m2PaeH6B0u9pXOn3NajVc8eW9j8TEdt3Zd2Voh3n7/9FxH9Kv7c7I6Ku1XNt/i33JO04fwdExORW5+n7rZ7bv/T3/lxE7N+1lVeGdpy/01udu2cjYkqr5/z9RVwQEZMiYsxsno+I+GPp/D4REeu3es7f39zPX0PpvD0ZEQ9ExLqtnnuptP+xiGjuuqrLIKXk1g024LfAMaX7xwC/mc1x789m/5XA3qX75wCHFv2dKu38AasCq5TuLwO8BixaevxX4JtFf48uPme9gOeBFYH5gceBNWc55gfAOaX7ewNXlO6vWTp+AWCF0vv0Kvo7VeD52wqoLd0/tOX8lR63+bfcU7Z2nr8DgD+38drFgRdKt4uV7i9W9HeqtPM3y/GHARe0etyjf3+lczAEWB8YM5vnhwI3AwFsBIwq7e/xv792nr9NWs4LsGPL+Ss9fgnoV/R36IzNEe3uY1fgotL9i4Dd2vvCiAhga+CqL/L6KjHX85dSejal9Fzp/qvAJGCuzeir2GBgXErphZTSJ8Dl5PPYWuvzehWwTen3titweUrp45TSi8C40vv1JHM9fymlu1NKU0sPHwKW6+IaK1l7fn+zsz1we0rprZTS28DtwA5lqrNSzev5+zZwWZdU1k2klEYCb83hkF2Bi1P2ELBoRAzA3x8w9/OXUnqgdH6giv/7Z9DuPpZKKb1Wuv86sNRsjlswIpoj4qGIaAmTSwBTUkqflh6/AixbxlorUXvPHwARMZg8CvR8q92NpX/mOj0iFihTnZVkWeDlVo/b+t3875jS7+sd8u+tPa+tdvN6Dg4kj461aOtvuSdp7/n7Runv8qqIWH4eX1vN2n0OSlOWVgDuarW7p//+2mN259jf37yb9b9/CbgtIkZHxLCCauoUvYsuQDNExB3A0m08Nbz1g5RSiojZtYupSylNjIgVgbsi4kly+Kl6nXT+KI1I/A3YP6U0vbT7WHJAn5/ciuho4KTOqFuKiO8A9cAWrXZ/7m85pfR82+/QY10PXJZS+jgiDib/68rWBdfUHe0NXJVS+qzVPn9/6hIRsRU5aG/Wavdmpd/fksDtEfF0aYS82zFoV5CU0razey4i/hsRA1JKr5WC4KTZvMfE0u0LEXEP8BXgavI/afUujTouB0zs9C9QsM44fxHRF7gRGF76p8CW924ZDf84Ii4EjuzE0ivVRGD5Vo/b+t20HPNKRPQGFgHebOdrq127zkFEbEv+P4NbpJQ+btk/m7/lnhR05nr+Ukpvtnp4PvlajJbXbjnLa+/p9Aor27z8De4N/LD1Dn9/7TK7c+zvr50i4svkv90dW/89t/r9TYqIa8hTobpl0HbqSPdxHdBy5fL+wLWzHhARi7VMaYiIfsCmwH9SvrLgbuCbc3p9lWvP+ZsfuIY85+6qWZ4bULoN8vzuNq+irjIPA6tE7lgzP/l/jGftPtD6vH4TuKv0e7sO2DtyV5IVgFWAf3dR3ZVirucvIr4CnAvsklKa1Gp/m3/LXVZ5ZWjP+RvQ6uEuwFOl+7cCXyudx8WAr5X29STt+fslIlYnX7D3YKt9/v7a5zpgv1L3kY2Ad0qDMv7+2iEiBgL/APZNKT3ban+fiFi45T75/HXf/80t+mpMt/Zt5HmvdwLPAXcAi5f21wPnl+5vAjxJvrr8SeDAVq9fkRx0xgF/BxYo+jtV4Pn7DjANeKzVtl7pubtK53QMcAnwpaK/Uxedt6HAs+SRrOGlfSeRgyHAgqXf07jS72vFVq8dXnrdM+TRisK/TwWevzuA/7b6vV1X2j/bv+WetLXj/P0KGFs6T3cDq7d67fdKv8txwHeL/i6VeP5Kj38J/HqW1/n7y+fhMnL3qWnkedYHAocAh5SeD+DM0vl9Eqhv9Vp/f3M/f+cDb7f6719zaf+Kpd/e46W/7+FFf5eObK4MKUmSJJWBU0ckSZKkMjBoS5IkSWVg0JYkSZLKwKAtSZIklYFBW5IkSSoDg7YkdROlfr33R8SOrfbtGRG3FFDLPRFR39WfK0ndiStDSlI3kVJKEXEI8PeIuJv83/BTgR3K+bmtVpWVJM0Dg7YkdSMppTERcT1wNNCHvJLpTEtjR8QO5ADeC3gjpbRNRCwOXEBeDGIqMCyl9MQc9v8SWKm0f0JEfA+4EFgXeBpYqPRZvYC/kBd/SsAFKaXTy3kOJKm7MGhLUvdzIvAI8Ak54P5PRPQHzgOGpJReLAXpltc8mlLaLSK2Bi4G1pvDfoA14f+3d4esVYZhGMf/1z7AELQtaRPGgskwhsVPYJvBJmabYPAzaFgZLi+4trIgYxaLZWVNu4KCuA2Lt+HcwjjgLHuQ9/D/wYGX69yH933axeE+HNar6jzJU+Csqm4nWev707MrVbXa97826MySNDkWbUmamKo6TbIL/Kiqn3Nv3wWOqupTz37tfB140NnbJNeTLF+Sw+wv4c/7egN42XPHSY47/wjcSvIK2AcOrvq8kjRV/hhSkqbpV79GOv3XQFV9Y7ZOcgg8AbYHP5MkTYZFW5IWy3tgI8lNgAurI++Ah53dY7a7/f2SfN4RsNlzq8BaX98AlqrqDfAcuDPkVJI0Qa6OSNICqaovSR4De0mWgM/AfeAF8LpXPs6AR/2Rv+XztoCdJCfACfCh85XO/3xx8+xqTyRJ05Wq+t/PIEmSJC0cV0ckSZKkASzakiRJ0gAWbUmSJGkAi7YkSZI0gEVbkiRJGsCiLUmSJA1g0ZYkSZIGsGhLkiRJA/wGpusz/2MsFioAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f56728f6090>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Numerical flux at y =  -0.505 is -0.42331216287\n",
      "Exact flux at y= -0.505 is -0.424589706231\n",
      "Abs. error  0.128516520089  magnitude (will change in parallel due to shadow space)  34.3695103704\n",
      "Rel. error  0.00373925955605\n"
     ]
    }
   ],
   "source": [
    "numeric     = np.zeros(mesh.specialSets['MinI_VertexSet'].count)\n",
    "ycoord      = np.zeros_like(numeric)\n",
    "analytic    = np.zeros_like(numeric)\n",
    "\n",
    "ids = mesh.specialSets['MinI_VertexSet']\n",
    "numeric[:] = tField.evaluate(ids).T\n",
    "\n",
    "ycoord = tField.mesh.data[ids.data,[1]]\n",
    "analytic = analyticT(ycoord)\n",
    "\n",
    "abserr = uw.utils._nps_2norm(analytic - numeric)\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",
    "offset = 0.5*(mesh.maxCoord[1]-mesh.minCoord[1])/mesh.elementRes[1]\n",
    "yspot = y0+offset\n",
    "ana_flux = analytic_dT_dy(yspot)\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",
    "\n",
    "# assuming order in the allgather is the same\n",
    "coords = comm.allgather(ycoord)\n",
    "numerical = comm.allgather(numeric)\n",
    "\n",
    "rank = uw.rank()\n",
    "\n",
    "if rank == 0 and make_graphs:\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 = analyticT(big)\n",
    "\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",
    "    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",
    "if rank == 0:\n",
    "    threshold = 1.0e-2\n",
    "    print \"Numerical flux at y = \" ,yspot,\"is\", num_flux\n",
    "    print \"Exact flux at y=\" ,yspot,\"is\", ana_flux\n",
    "    print \"Abs. error \", abserr, \" magnitude (will change in parallel due to shadow space) \", mag\n",
    "    print \"Rel. error \", 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)  "
   ]
  }
 ],
 "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.14"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}