{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "### Example of a 1D Channel driven flow\n",
    "\n",
    "Model Setup\n",
    "-----------\n",
    "\n",
    "2D, Stokes Equation with noslip BC at top and bottom boundary and a lateral pressure gradient driving the flow, a.k.a. Poiseuille Flow.\n",
    "\n",
    "\\\\[\n",
    "\\frac{\\partial \\tau}{\\partial y} = \\mu \\frac{\\partial^{2} \\mathbf{u}}{\\partial{y}^{2}} = \\frac{\\partial p}{\\partial x}\n",
    "\\\\]\n",
    "\n",
    "\\\\[\n",
    "\\nabla \\cdot \\mathbf{u} = 0\n",
    "\\\\]\n",
    "\n",
    "with $ x_{a} \\leqslant x \\leqslant x_{b} $ and $ 0.0 \\leqslant y \\leqslant h $\n",
    "\n",
    "Boundary conditions:\n",
    "\n",
    " * $\\mathbf{u}(x,y=h) = \\mathbf{u}(x,y=0) = \\left[0,0 \\right]$\n",
    " * $P(x_a) = P_a$\n",
    " * $P(x_b) = P_b $\n",
    "\n",
    "------\n",
    "\n",
    "A 1D solution in $y$-axis, described by\n",
    "\n",
    "$ \\mathbf{u}(x,y) =  \\left[ \\frac{1}{2 \\mu} \\frac{\\partial p }{\\partial x} ( y^{2} - h y ), 0.0 \\right]$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We implement the above boundary conditions using:\n",
    " * a `DirichletCondition` object for $\\mathbf{u}(x,y=1) = \\mathbf{u}(x,y=0) = \\left[0,0 \\right]$\n",
    " * a `NeumannCondition` object for $P(x_a) = P_a$ & $P(x_b) = P_b $\n",
    "\n",
    "The `NeumannCondition` object, used with the `Stokes` object, defines a stress along a boundary such that:\n",
    "  * $ \\sigma_{ij} n_{j} = \\phi_{i} $ on $ \\Gamma_{\\phi} $\n",
    "\n",
    "  where  \n",
    "    * $n$ is the surface normal pointing outwards,\n",
    "    * $ \\sigma_{ij} = \\tau_{ij} - \\delta_{ij} P$ is the prescribed stress tensor, which is multiplied my $ n $ at $ \\Gamma_{\\phi} $ to produce $\\phi_{i}$, a surface traction on the given boundary."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import underworld as uw\n",
    "from underworld import function as fn\n",
    "import glucifer\n",
    "import math\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"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Setup parameters\n",
    "-----"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "# domain height\n",
    "h = 1.0\n",
    "# Set a constant viscosity.\n",
    "viscosity = 1.4\n",
    "\n",
    "# position of walls and associated pressure on walls\n",
    "xa = -1.0\n",
    "pa = 4.0\n",
    "\n",
    "xb = 1.0\n",
    "pb = 3.0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "mesh = uw.mesh.FeMesh_Cartesian( elementType = (\"Q1/dQ0\"), \n",
    "                                 elementRes  = (128, 128), \n",
    "                                 minCoord    = (xa, 0.), \n",
    "                                 maxCoord    = (xb, h))\n",
    "\n",
    "velocityField   = mesh.add_variable(         nodeDofCount=2 )\n",
    "pressureField   = mesh.subMesh.add_variable( nodeDofCount=1 )\n",
    "appliedTractionField = mesh.add_variable(    nodeDofCount=2 )\n",
    "\n",
    "# initialise velocity, pressure field\n",
    "velocityField.data[:] = [0.,0.]\n",
    "pressureField.data[:] = 0."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "jWalls = mesh.specialSets[\"MinJ_VertexSet\"] + mesh.specialSets[\"MaxJ_VertexSet\"]\n",
    "iWalls = mesh.specialSets[\"MinI_VertexSet\"] + mesh.specialSets[\"MaxI_VertexSet\"]\n",
    "allWalls = iWalls + jWalls"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "-0.5\n"
     ]
    }
   ],
   "source": [
    "vBC = uw.conditions.DirichletCondition( variable = velocityField, \n",
    "                                        indexSetsPerDof = (jWalls, allWalls) )\n",
    "\n",
    "dp_dx = (pb-pa)/(xb-xa)\n",
    "\n",
    "# This stress is multiplied by the wall normal to produce a traction force.\n",
    "#############\n",
    "# Remember total stress = deviatoric - isotropic. \n",
    "# Thus +pressure is a negative stress.\n",
    "#############\n",
    "\n",
    "# The left wall normal unit vector is (-1,0)\n",
    "# The right wall normal unit vector is (1,0)\n",
    "\n",
    "# (-press) * normal_j = surface_force\n",
    "appliedTractionField.data[mesh.specialSets[\"MinI_VertexSet\"].data] = (pa,0.0)\n",
    "\n",
    "appliedTractionField.data[mesh.specialSets[\"MaxI_VertexSet\"].data] = (-1*pb,0.0)\n",
    "\n",
    "nbc = uw.conditions.NeumannCondition( fn_flux=appliedTractionField, \n",
    "                                      variable=velocityField,\n",
    "                                      indexSetsPerDof=(iWalls, None) )\n",
    "print dp_dx"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "stokes = uw.systems.Stokes( velocityField = velocityField, \n",
    "                            pressureField = pressureField,\n",
    "                            conditions    = [vBC, nbc],\n",
    "                            fn_viscosity   = viscosity, \n",
    "                            fn_bodyforce   = 0.0 )\n",
    "\n",
    "solver = uw.systems.Solver( stokes )\n",
    "solver.solve()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<img src='data:image/png;base64,iVBORw0KGgoAAAANSUhEUgAAAoAAAAHgCAIAAAC6s0uzAAA+SUlEQVR4AezVMQHAIAAEsda/6GfCArckJvJv+wCAtwQMAAEBA0BAwAAQEDAABAQMAAEBA0BAwAAQEDAABAQMAAEBA0BAwAAQEDAABAQMAAEBA0BAwAAQEDAABAQMAAEBA0BAwAAQEDAABAQMAAEBA0BAwAAQEDAABAQMAAEBA0BAwAAQEDAABAQMAAEBA0BAwAAQEDAABAQMAAEBA0BAwAAQEDAABAQMAAEBA0BAwAAQEDAABAQMAAEBA0BAwAAQEDAABAQMAAEBA0BAwAAQEDAABAQMAAEBA0BAwAAQEDAABAQMAAEBA0BAwAAQEDAABAQMAAEBA0BAwAAQuAEf9spf55Kc6MNdPfOFm30RARIiIUFIyxKQACJBBCsQEQloCQEJRIK0IoKMGwAugZSIW0CQgQRXMy5c/+yyu895z7sz455dfo96evy6XWW33aceAAAAACwEAgYAAAAuAAIGAAAALgACBgAAAC4AAgYAAAAuAAIGAAAALgACBgAAAC4AAgYAAAAuAAIGAAAALgACBgAAAC4AAgYAAAAuAAIGAAAALgACBgAAAC4AAgYAAAAuAAIGAAAALgACBgAAAC4AAgYAAAAuAAIGAAAALgACBgAAAC4AAgYAAAAuAAIGAAAALgACBgAAAC4AAgYAAAAuAAIGAAAALgACBgAAAC4AAgYAAAAuAAIGAAAALgACBgAAAC4AAgYAAAAuAAIGAAAALgACBgAAAC4AAgYAAAAuAAIGAAAALgACBgAAAC4AAgYAAAAuAAIGAAAALgACBgAAAC5ghYCJ6OrXBAAAAB5lgRkrL9e8TPn4i9XDftUXs3eTxuade73veq8D5Jk+0qf1P29It0S1VHs0Ys+83To9intnvraYqOGPtmFMrIFtJRN5xtPZn9U4SeKP+9THkOOk85iWIfLoeH7zq30gZEs7fHyRe7HPeuWzJPUQ74S8kRe8lYSeu0tPvKB9tsypPfbYz6veC1P78yRE+qn9HOMuIYX7r0fG9Aw95HQlRX/idYzOfj6mpc3FIC91WttxJXyIemRt48j+sm3BJ2uLOjHNkvfquLZhrtN9OyYZ13bjlWO1EWK7lBeZX99DeGjE2objmALvfBh55ZtNHYlaktKjhtmLvs6xjkZI/2a2tIcpSa9h+plR0S+8hfijeZHe334OJQKtJFD68dX+f/7i/W0JiwRMpUi5F8WmsuJ7oz/WumEk/8KsexftJOCtHSz38S3qWLysnX9Yg253+X5yt6XtXdqWG+lIOvt4UubTpwAA8A7A6f7amaTYSfXkcHCaYyf1qJVSv8u/WqMLs6lxqpUaIr21KKtcPVkZrV8hnUO10BouVPd06MV63AlSvussnIt7jXpVyBxsa1pZwRcJ2PZDt8hkWU2sjUmor1gkKhtQuqQpNlw3r4vQd1oH7Mmpbb+7wpODzzV8/ETD9xQhMq38Y//wtqUHBQAArw29ln1paEkxtEpYLeeOPA+JfvElk1Rx5mo9K+E0FePqUeJa0WvNLRpLYW+ryCUsaYGiBHtoQpVRNNRx8XpLr6LdeNdqXjRh/bP65pXJqU21hEUC5qpVOZy6tfUtd93YMLGpUbZMT6OYeqnfyYRqbXWw7T83Xda0pudttm8/VRocPGnYToksaeRvM7Y8skKyw2YPa1PkYQAA8Omm19fzp1oiow6Lf6U2ai1P8fKn9ujdK6QV1130LTa1R5SmtWpa7WgaVnOYL7XB0ig12AdbuRYzmNA1pfSKT7zQi+xLrC3sIiG1ntdUdZTIm13DyzZ5kYB1b0Kx3PZHTbyHHdtT97SdcjzqEiUPManbYfao0b7WqCPtY7HGLQ1nj1LScD4Pa7J/a2x/ZtMDAMC7Dd026xE+iZ4fUS+m1oxavnX7eq3NxZJ5J9VkraU6wOUnfVrvNVOJsruLfDW36Fjqb2lW37wY7zoPjxqmSCIajpWYk+p/NUTsqyW9ji6eaQWrBGwONlOy3TfdE/Ocyni3/U4mdvtuveGireauCUfz+aaO4+1eSvd3e9Rhn7R9mT6Am29DxjmIfdJ2BwCAzwpjTcx4YWx+bQHqOimnXlBVsWTF0xR78Ln4lGr1l+dVgWaGMKMNp67wKPb75jPpMkzDlFfrhmkm2bYXZllbyiYh5m7y4s01JxGV7bDEt8kqAesO6c7tviXtAP1R6Sa2zeTRvl3Aau46ptCs1VNnOzzsKmWhWuZI3r43T6hfQ/uqbPAt28LEAIDPEKN9Xb3dwfHHlmRM3b7DcBWqKpZz+Y2QbdvVwVKqTZwxRf2zeKDqWuepIqkG2MUc1Z3ytLCV7KZe17lMp8nq4xfq9uLV2jSsD8mL+762fu+L5iklLj2jqmGTcRNnffRKn9aG7CHtv/vP/tt/yxi7bLxc2jDb5aulzZltvHndTrPK26+wLNH+q78NjueWfIsDTR/b6Qm1NTzFi4/+umbLX3z/z/VaM9fLb/7p5Tf+uGiur/x+zUSV//v8x4smeu+n9Voz1/+/9+N6rZnry5/74ZqJKt/90g/WTPTRBx/+5GsfrpnrN9/+Tr3eYMJcyyb+8L1v3Yo6q3nVZOz97steCy2kX9TuIkki/suPvi7FXq9WO6VgywAdowljjNu3atj+rNcrNzfvelmj9uhgGWAj//6zD/7x86/Kn0V73LL8os4lITqRP5UBr0zOq6hr5bc+B1H59RfccM2L7aDa3R81X4YRPcrOVI+63XPOPCw38kexjVF95PhV9iTx4dhRt/Wkdxva+cpPH29Mb92//bTOY8gUeDKmZYg89js4Br7uah8I2drP8exF7sU+65XPklgBeMaMz3/BW0noubv0xAt6JeXUHnvYyqLVr/jzJET6tYdzFFn56z9Q2buWoYecrsSrqpSXSHKcOtK2AqAz9qVOazuuhA9Rj6xtHNlfti34ZG3JMXmWvFfHtQ1zne7bMcm4thuvHKuNENslHkXY5vWQmOXOcUyBx82fDiimox67hZpTzvYttZwl9Wfo8NH2kWn2XhS1KrcBxfeE4vPT9431U/rpUKQy6dq2RI2R0f/65fsLzFjZF8zxBL4xcTn+5YxufHxH+EaeNtPdIQAAAF4DHq/c2fV5o2pvNwbMljCZhu/jGZ9EDlOIpGf/0DW79PKaaZ8Hj9v5XBmfctF+AwDAZx1++CmNPZOwjwPOK7g7oZmYo813Sj13Xaeca82wUMCc3rSUjUiu9t58GHNq4RZV7zZ4DqEYGv0ynuZtrVEtQ0s4pNGQfiDU19DONIcQjA4AeNehp+x4Cw5jneTqBfxoU2pOjbFWe0WTzPK0ts8t2Ost75qneJhUbk4alurL205c5Dnt4l8u+nCXSaj2l6bmCNw3LrRpZgmsy5LwtaV8lYD33bd58qV4Th/ZkdjmFD1ov/QQmvO20cFNwCfmHk1sQs2jJgfbn9m+NiWNJ8Jj4mxiOBgA8D+B1DqOKmvFklO57U6N8iyjkmijVqpQt0HbW6qpU73dqcrBC7XljFn8fxlgo9XHrn8NYXVwYfJ6H8uvT2q/aJjFx1XXZeE+LhMwuUcHDdu+kjtVtthMlseYGofhLsX9EQ1v6aO4ocnJwT6AZvVaz33gYADAOw8dhPes6PNY0sqdfLnlqpzs2+q0Fl1ulf9YvylnCwerKb3sV19yr+3uVG2EbnTkHhapqygbFW65Rw3XbLXB6w5ilYCz29rhm2WzUyeB9ZD0qJ1SiajJfHeiZsseVmgn0qO2Ll24FQAAGqOoKBTb7S4tJqm68phHy3Iuq6FYTp5uw8yy9sgLeZjVwjyKvYetRzq9obWfq18LVcuKv2vnpGEbU5PsC0v9QgFzs+zhUYaS/05VeoxqYwbN3wixYdnZfZTsvnxH/nXEUrl/HzAxAAAM5GJsHjQ/8mBWMge78ah771BTefR0sq9V55B3F65nkyjq9i3RaDHSUgfXhRTVc+ndvhQajfRWWSXg4psQr3iQWBYqhz6zSoeQOIVpp472PZ3IOrOtp69AFtDOnwb7DtMt2jwAAPjUcCYwmmsrZyVwsmPYtUl3cvA8iZpeha61OsaTjWyNvA4N2Xbi6uDNJd2mXcq+aJ5ZnzzsCN06Oj4P6THP3y46Zr43DgAAwCeFDu5kK98UtZfPii8dKzF9klJPT8lngm8r4W3wct1UNL0sP+U7jt7cOB7IYyaeZ8+Z+cYhw8EAAPBGiHI6FP5ee/ngBjpVBd2Q9p1pWco/HyfgucqvVK+xL52NHh/WtiIa/NjmPDLsVNnn5oWDAQDgteFJfIKWVybtuFNq6fCUz0bTcB0M8k7W8lUCJoqtYLmz3qdNsw3ikOhp407ayb53TMyH+y1J37I1AACAT8BUWo/1OFSQqm4zNE8651zpc62WJEyaiojnkPjT3UNDgeeFm7FKwPflJxvQTJx26LTRd46Gfeek3nYmx8atxbRUdPNDAQAA8MYYa2wv4dRrfJaxlmceSjYPf9PNi5suJIQH3+f8NnbZBswC/i975RNbxXHH8TfP9qkRVZtDm1MRlWhVkFrXcOFkqlSWqOw04VIiWpdTCT0gJRfHXBroofEpUUEKSXvAuUAPOE2hQrRISFQER4kSqVYR58Kh6S1/TsY73bc7OzszO/v+mJk3qPl8hNbzZuc339/vN8N+46KM1nBE8+mskbZNtpp39afZNGdzrwf3MWPHerFhAIDIyOo7XH/jhfsRFvaXWnlwZaump+qAwoMLaxV2iKzN3pWLz3gNuGObsWqeYZNm6c1+1G+9bvnIxNkVAADacO1Vup9h2RIjzEFjmSz82PVmW0vaErnDjNF8e4zLgB2XFS1lyqoh5QKP9XbsDsk+ZxTMkgEAICai+YGXykT7f4eFHVmbrg6s3Lf3LCzWMRbXhsdY9eQYtQyX9Tprr3xRPrun7+rZ7Dd7Gl0RAZvUffFOrfXagXhHMXHsmh5vXTgUcuum1k8v1lp/PhJVa3L2vB4/vHUirtb3X621/vlyPKGpby3r8ea/fxe1qKkdddM2P38jqtaTTyzq8X8/W42qtfepn+nxxw8uxhOa+85zevzXu5ejFvXLffN6/Mf3r0TVWj44p8e/vXE9qta5Z2b1+IV3bsYTWjtafWM7nWfevj1sWG2h9qCk8A3RmF4/sU+PZ859UJpP7SGysmRZebb07B0VIWV0OSFEtvRt87f7tAbaX4tnb94YWAs6nh30c4CKL6r+WUqU/+o1cuAmXvWRBp5N9H15hJKdrlbrZfhsh2l11dX+p+CtYoSSfZuIjhz1YoQ6UzFqlwYUqL5IUjhfJz0jy8+K7GRS6J+eEFl+u+r1xbMXkpXfI/VVEsVKNaNDvJnkgeWaQt2/Rm+rBsXmZqpObs1MZCNqmNzslXWxOmFPbioTJ1aYvWrmZml5+9bcxM6tpeQq2yqk7JKZpFm+CpHWoMrNOg4n0Gl+Vmcr7IvRqS9G1W2zXnUZdLaNnrh9UwP7lNsuhlW1MM7UF2WX7Nx884A2Tk6PwRlzumPQGIzlvho5xAAAAL6sDG0Fom+oaF8WlfEasDC8Vkr36e+A3Ga/AQDgy0O7OTT9VT4ejjJGAxaOE3fsDjwmDQEAgP9DRMN/krvO2AzYqLTZBtFYAwAAEBPThmUK7xmbAQMAADwuiNQJdDBgAACAJGDAAAAACcCAAQAAEoABAwAAJAADBgAASAAGDAAAkAAMGAAAIAEYMAAAQAIwYAAAgARgwAAAAAnAgAEAABIwmToBD90z95yZ7JU9sbRevONqvXYghtDEsWvOzNaFQ5GKmnj2kqv17pFIWpOz552Zh7dOxNL6wYqrtbEcQ2hq5ylnZvP+q5GKmtrhtmvz8zciaT35xKIz88lnq5G09jzlXrmPH1yMIfTj3Yedmat31yIVtTgz78z84f0rkbSWD845M2duXI+kdXZh1pk5/s7NGEJrR92v68Lq7UhF3XlhvzOz79wHkbS2jZBSRtcQIlvalf/Rv4cd5M/eX1EMGq/aQlSgGLy4bYEpKuo1cnC2PvWRS3YG6nU12FbJQu9Q7VOsl+GzHabVVVeHPFOzihFK9m0iOnLUixHqTMWoXRpQYHkbpTTG9kz+P7t8ZlLon56Q3nwxI82oXkhWLCq/EL019Q51iDeTPLBcU6j71+ht1aDY3EzVya2ZiWxEDZObvbIuVifsyU1l4sQKs1fN3Cwtb9+am9i5tZRcZVuFlF0ykzTLVyHSGlS5WcfhBDrNz+pshX0xOvXFqLpt1qsug8620RO3b2pgn3L/i6Emqx6qEOdo7NyMvjVL3jg5PQZnzOmOQQMAAAAcMGAAAIAEYMAAAAAJwIABAAASgAEDAAAkAAMGAABIAAYMAACQAAwYAAAgARgwAABAAjBgAACABGDAAAAACcCAAQAAEoABAwAAJAADBgAASAAGDAAAkAAMGAAAIAGTqRNw6Z65553PTu8Nr/XSul/r9QNhhSaOXfPOb104FLyoiWcv+bXefT641uTBN73zD2+dCK81veLX2lgOKzS185R3fvPBSvCipnb4G7X5xfngWl9/YtE7/8mnq8G1vvfNI975j+5fCiv09O7D3vmrd9eCF/WLmXnv/FvrV4NrvXxwzjt/+sb14FpnF2a988fXboYVunzU/11dWL0dvKj3ju/3zs+c/TC41qMgpJTRNYTIlnblf/TvYQf5s/dXFIPGq7YQFSgGL25bYIqKeo0cnK1PfeSSnYF6XQ22VbLQO1T7FOtl+GyHaXXV1SHP1KxihJJ9m4iOHPVihDpTMWqXBhRY3kYpjbE9k//PLp+ZFPqnJ6Q3X8xIM6oXkhWLyi9Eb029Qx3izSQPLNcU6v41els1KDY3U3Vya2YiG1HD5GavrIvVCXtyU5k4scLsVTM3S8vbt+Ymdm4tJVfZViFll8wkzfJViLQGVW7WcTiBTvOzOlthX4xOfTGqbpv1qsugs230xO2bGtin3P9iqMmqhyrEORo7N6NvzZI3Tk6PwRlzumPQAAAAAAcMGAAAIAEYMAAAQAIwYAAAgARgwAAAAAnAgAEAABKAAQMAACQAAwYAAEgABgwAAJAADBgAACABGDAAAEACMGAAAIAEYMAAAAAJwIABAAASgAEDAAAkAAMGAABIwGTqBCy6Z+71eZud3htS66X1VqHXD4Sta+LYtbZXW6s/CSn03J/6vN36y/MBtSZ/9Faftw//8euQWtMrrUIbywGFcqZ2nmp7tflgJaTQV/u1aPOL8wG1vvaVxT5v//Pp2wG1vvuNI22vPrp/KaBQztO7D7e9uvKvtYBCP5+Z7/P2zfWrAbWWZuf6vH3l738LqPX7hdm2V8fXbgYUyrl8tPW7Or/6XkCh27/a3+ftD89+GLauR0FIKaNrCJEt7cr/6N/DDvJn768oBo1XbSEqUAxe3LbAFBX1Gjk4W5/6yCU7A/W6GmyrZKF3qPYp1svw2Q7T6qqrQ56pWcUIJfs2ER056sUIdaZi1C4NKLC8jVIaY3sm/59dPjMp9E9PSG++mJFmVC8kKxaVX4jemnqHOsSbSR5YrinU/Wv0tmpQbG6m6uTWzEQ2oobJzV5ZF6sT9uSmMnFihdmrZm6WlrdvzU3s3FpKrrKtQsoumUma5asQaQ2q3KzjcAKd5md1tsK+GJ36YlTdNutVl0Fn2+iJ2zc1sE+5/8VQk1UPVYhzNHZuRt+aJW+cnB6DM+Z0x6ABAAAADm0G/D/26idEqiOB4/h73dOX/HNP8eBeFESEOcTemAWHEP9APExaYrzEnc0YJMsEZkEPEgUN2YjkMidFQnRzUWHFrDDgPxC9CEazEvHiePAsgldhvaivdvrV+1NV3a/e634zVTPL90PoVE/3q19V9fP9AADAIqKAAQDwgAIGAMADChgAAA8oYAAAPKCAAQDwgAIGAMADChgAAA8oYAAAPKCAAQDwgAIGAMADChgAAA8oYAAAPKCAAQDwgALWNPbfcZbV3HPNXdaO886y/i+1/njQXdZbX/ve7vI2vn6ns6yvPvjE93aXt0uTm5xl/T79vu/tmkIhxKJnhGF0aM38/7L3RYPGsceWeaKjo/mXiybJXo23RtCBe/Y1R8fHkknCfBJhX78aquQ29163Z70+O26dJMlIB7YtNz/7xRZ0+S/5PPH3Re9ZVR6MbP2nJevV7eni1aan2ncjPYOR9oz9AF/NHS48FmMQBsJ6Y7RWH7FnvXw6M9Ap5T9laB51a8W0Lei/p+y3sbwbhVDG+l/m/2XL10iEf3hzjyXr2fNzgfL9+LV7b0TxjPIJIbpnl0XEb/OP8tx1K3fbD/DBkwvZJHLaZBBPHsW5yRrij7UgfSVb1+6yZ115NNv3lNTEbLNyU0ZQtsEv2h1L0On/XDXWZmb1Pbc0MV9D/MdDm7dbsr6/eSNZbXqJPKV0EOjHm04utIH86ERni/0Ap2ZvZVNF+WpD/cYI8hsjPW11v90LRXBxYsye1TlzV9mOdpL5ZotujOSPyYH8OrXREtQ+eT/Qz03fSHfah/s2OGjGYKkVcL8nl+yeqg/rigVcMpsaWruAh3lYD1vAZTsKFrCAh9zy4AWs7WKALQ9TwIv4m4aDnlKwgAWcve1zifIsq1/AciyfufJ5HRV8p34Bi56rqqxtuAJWE/t358IVcKD9WL1bDheqgI2fw7jQOPyaBZyvtudMzHOrXcDpJcZPo69taRRww0EGAAAwUMAAAHhAAQMA4AEFDACABxQwAAAeUMAAAHhAAQMA4AEFDACABxQwAAAeUMAAAHhAAQMA4AEFDACABxQwAAAeUMAAAHhAAQMA4MHyKODGd3PusvbfcRPU3HNt/j9HWTvOuwmaN/Lhj+6yRn9wltVa9Y2joLe+drapd9+edJb13qrPnWWNr9/pJuirDz5xtqlvt33sLOunTzc7y7o0uclN0O/T7zvbVEWhEGLRM8IwOrRm/n/Z+76DxrHH9nmio6P5VUWzZa/GW/3LjQP3SrKOjyWThPkkwrp+LVTJbe69bs96fXbcOkmSkQ5sW25+9ktJ1uWJZJ74+6L3rKoNRrb9bA96dXu6eLXpqfbdSM9gpD1TkjV3uPBYjEEYCOuN0Vp9xJ718ulM9VPSfsrQPOrWiumSrBenLbexvBuFUMb6X+b/ZcvXd9740h707Pm5QPl+/Nq9N6J4RvmEEN2zyyLit/lHee66lbvtWQ+eXMgmkdMmg3jyKM5N1hB/rAXpK9m6dpc968qj2b6npCZmm5WbMoKyDX7R7tizTv12VV2bmdX33NLEfA1BcPCj7fag72/eSFabXiJPKR0E+vGmkwttID860dliz5qavZVNFeWrDfUbI8hvjPS01f12LxTBxYkxe1bnzF1lO9pJ5pstujGSPyYH8uvURntW++R99dz0jXSnfbhvg4NmDJZUARc8uWT3VH1YVyzgktnU0NoFPMzDetgCLtuRMk+9Ah5+y4MXsLaLAbY8TAEv4m8aDnpKJRusXsBqq/W5RHmW1S9gOZbPXPm8jgq+U7+ARc9VVdY2XAGrif27s2dtwxWwuraCLYcLVcDGz2FcaBx+zQLOV9tzJua51S7g9BLjp9HXtjQKuOEgAwAAGChgAAA8oIABAPCAAgYAwAMKGAAADyhgAAA8oIABAPCAAgYAwAMKGAAADyhgAAA8oIABAPCAAgYAwAMKGAAADyhgAAA8oIABAPBgeRRw47s5d1n77zjLau655ihox/nmjn+5yRr58Ec3Qd2s0R+cZbVWfeMu680pN0Hvvj3pbFPvrfrcWdb4+p3Osv7253E3Qd9u+9jZpn76dLOzrEuTm5xl3f/7n5xlVREKIRY9IwyjQ2vm/5e97ztoHHtcOlV0dDS5qmi27NV4q3+5ceBeSdDxsWSSMJ9EWNevhSq5zb3X7Vmvz45bJ0ky0kHhlpu7/l16gK8vTyTzxPOL3rOqNhjZ9nNp1qvb0wWrTU+170Z6BiPtmZKgucOFx2IMwkBYb4zW6iP2rJdPZ6qfkvZTqlsOw9aK6dIDfPnitOU2lnejEMpY/8v8v2z5+s4bX5ZmPXt+Lvt+/Nq9N6J4RvmEEN2zyyLit/lHee66lbvtQQ+eXMgmkdMmg3jyKM5N1hB/rAXpK9m6dpc968qj2b6npCZmm5WbMoLk4K/tTukBnvrtqro2M6vvuaWJ+RqC4OBH20uz/nHzRveq9BJ5Sukg0I83nVxoA/nRic4We9DU7K1sqihfbajfGEF+Y6Snre63e6EILk6M2bM6Z+4q29FOMt9s0Y2R/LH7ndtTG0sPsH3yvnpu+ka60z7ct8FBMwZLqoALnlyye6o+rCsWcMlsamjtAh7mYT1UAdu2HGYzhAtSwMNvefAC1nYxwJaHKeBF/E3DQU+pZIPVC1httT6XKM+y+gUsx/KZK5/XUcF36hew6LmqytqGKGC12wq7s2dtwxWwuraCLYdaJ9UoYOPnMC40Dr9mAeer7TkT89zqFbByifHT6GtbGgXccJABAAAMFDAAAB5QwAAAeEABAwDgAQUMAIAHFDAAAB5QwAAAeFBWwP9jr25Cozr3OI7PMy8rLfYuiov2ggiCgoUaXxaRQjBQoTHa6qZaNWAvJDQLe3HRULuwuYWC2dxgaIlZGReiIMNVI1gtVOpNGjC4ibncrlWwy4DdTHKenjnv58yZ58wkz5nHwvezGE9m/s/z/59njvMDAAA5IIABADCAAAYAwAACGAAAAwhgAAAMIIABADCAAAYAwAACGAAAA8qmBwgVv/0ts8Ya3amn17n5jEbj+3XdV+nMPXXB6nSfnkZHb2TWrN7+VEuv8oGpzJqVR8N6enWNZTR6el5LI1tlS8ZWtRdjehptyj6c2h+XtfR6c8NAZs3L5ataem3ffFxd8OTZdS2NbL3bjqkL7ixVtTQ61dWfWTP564yWXiM9BzNrLjy4r6XXpcM96oLB6kMtjWw3T3arC/qvzGlp9N/BvZk1XRMLuu5rnYSUMvceQlgjW+1/gr9bvbBf6/+K4oWlMHqDj5ot8RaK7OK0guIXs9a4/awIbxMR1sjsadO6N78oDdwN07fpJt7H/sVabrl05Jr9unr7RLiPUy/bmbbFWy6///3KL58rj9o/1ebfQpO1yVsuv/vdyuJX6ctTNxEF2e6D4V9U3vmy9vzi2k9JtHHUlY1DtVeT6u/UfRqljFzH37H/Z7uvlhTBn4klb71x+vflaek8E0G981pfYjlF7i9EvSbcoRAsSZ3EXujWON29mvfe/sRN32ATd1vvwtk8Omr97USj+CSyYVUwW9+Oj930TT2laMfgZoOBU2bzJkmsFe4O/9h36PL8TONssV6p5xbZxJshPlvilr/u/WD0px/Daf0l7ilFh4x+Nd7mMnbhzybis4ULf/ioZ6j6c3RaK5xWxB+MQvhg+KcdvV/vYQimbTiTW6e7+6fnCrEvveFbVj8Y3pv+GXpLwlVBu8fDe9z0jZ5b/Ebq2y6e3dWBZCz8VQK4jR/rdQewn3mdCODWNtEQwP4dFToQwC0ctbYAVt6y5gBe7ymJdk+p0IEATvyW5RrAydZ5BrBiNu0BHP+VzzGAk9PmGcCNh59fACfPLc8ADmd7PQK42IEeAAAggQAGAMAAAhgAAAMIYAAADCCAAQAwgAAGAMAAAhgAAAMIYAAADCCAAQAwgAAGAMAAAhgAAAMIYAAADCCAAQAwgAAGAMAAAhgAAAPKpgcIFb/9TV1gje7U1uvcfEav8f26epXO3FMXrE73aet19EZGr9ufamlUPjClLlh5NKzrpspdYxm9np7X1auyJWOr2osxbb02ZRxR7Y/LWhq9uWFAXfBy+aqum9q++bi64Mmz67p69W47pi64s1TV1evU7n51weSvM1oajfQcVBd88+C+rpu6dLhHXTBYfair182T3eqC/itzunrNDu5VF3RNLOjqtU5CSpl7DyGska32P8HfrV7Yr/V/RfHCUpi+wUfNlngLRXZxWkHxi1lr3H5WhLeJCGtk9rRp3ZtclAbuFqLp23QT72P/Yi23XDpybfXWCfcwvX2cetnytK3fcvn971d++Vx51P6pNv8WmqxN3nL53e9WFr9KX566iSjIdh8M/6Lyzpe15xfXfkqi1aOubByqvZrM/E7dp1HKyHX8Hft/tvtqSRH8mVjy1hunf1+els4zEdQ7r/UlllPk/kLUa8IdCsGS1EnshW6N092ree/tT9z0DTZxt/UunM2jo9bfTjSKTyIbVgWz9e342E3f1FOKdgxuNhg4ZTZvksTa+qef7Ts0NT8j02aL9Uo9N3+TcIb4bIlb/rr3g9Gffgyn9Ze4pxQdMvrVeJvL2IU/m4jPFi784aOeoerP0WmtcFoRfzAK4YPhn3b0fr2HIZi24Uxune7un54rxL70hm9Z/WB4b/pn6C0JV7k1j4f37J5YCEf1zy1+I/VtF8/u6kAyFv4qAdzGj/W6A9jPvNwDuOVNNARwJHpzD+AWjlpbACtvWXMAr/eURLun1IkATvyW5RrAydZ5BrBiNr0BHF2edwAnp80zgBsPP78ATp5bbgEcm+31COBiB3oAAIAEAhgAAAMIYAAADCCAAQAwgAAGAMAAAhgAAAMIYAAADCCAAQAwgAAGAMAAAhgAAAMIYAAADCCAAQAwgAAGAMAAAhgAAAMIYAAADCibHiCm+K//Kz61Rndqa3RuXl1gje/X1at05p66YHW6T1uvozdUjW6f0NXIVj4wpfh05dGwtkZdY+qClafndfWqbMnYqvZiTFuvTaojqr2a1NXI9rcNA4pPXy5f1dVo++bj6oInz67r6tW77Zi64M5SVVevU7v7FZ9enp/R1cg20nNQ8ek3D+7ranTpcI+6YLD6UFevmye71QX9V+Z09Zod3Kv4tGtiQVej9RNSytx7CGGNbLX/Cf5u9cJ+rf8rnIuGj5ot8RaK7OKGguI/56x/d8eairBGZk+b1r3JRWng7uqVD1vYxPvYv1jLLZeOXFu9dcK/L/99945am1bLLfuv/qm2+J2GBe18y6mbiIJs98FwLip/H6k9v7iuUxKtnlJl45CXvsobdJ9GKSPX8Xfs/9nuqyVF8GfKkvr7zjsyuqq+xHKK3F+Iek24Q7gkdRJ7oVvjdE+vCbb1LpzNo6MmZmucRDascms+3HF05n/VZrNFK/1NgqYi0SjcJNKxECn+bN+hqfkZmTZbrFfqufmbhDPEZ2tyy/60/hL3lKJDRr8ab4mMXfizxb6OxMLE4VvhtCL+YBTCB8M/7ej9eg9DMG38TP5zqvvw9Gzs3LyL+LesfjC8N/0z9JaEq9yax8N7dk8sJJ58mXbLi2d3dSAZCwRwk9/uTgRwy5toCOBI9BLAawlgDack2j0lAniNAayeTW8AR5cTwGsI4JRzyy2AY7O9HgFc7EAPAACQQAADAGAAAQwAgAEEMAAABhDAAAAYQAADAGAAAQwAgAEEMAAABhDAAAAYQAADAGAAAQwAgAEEMAAABhDAAAAY0GoA/8le2cRodZVx/D0z76xsMNqFdiVhgUZIdAQ2rAZTMwlmprZspEFHVlJckLQL6bCx4MLOqo2QlFYXTDfggqkVDEFJSDCUadq0iRORtbCw7vqxGuYe73s/zj1f9/2Ac+5Z9Pdb3Dnvuec5/+d5zpn7BwAAgIBgwAAAAAnop07AYOrMvSFvs9O7Q2q9tN4q9Pr+sHVNH73W9mpr9cchhZ7705C3W395PqBW/4dvDXn78B+/Cqk1u9IqtLEcUChnZvuptlebD1ZCCn11WIs2vzgfUOtrX1ka8va/n74dUOs73zjc9uqj+5cCCuU8vfNQ26sr/1oLKPSzPQtD3r65fjWg1sm5+SFvX/n73wJq/X5xru3VsbWbAYVyLh9p/a4urL4XUOj2L/cNefuDsx+GretxEFLK6BpCZCd35H/U73EH+XPwVxQD51VbSBUoRi9uW6CLimaNHJ2tT33ikq1B9boePFLJQu1Q71Osl+GzHafVdVfHPFO9iglK9m0ienLSixHqTMWkXRpRYHkbpdTG5kz+n10+MynUT0/IYL6YkXrUICQrFpVfiMGaZocmxJtJHliuKdT9a9S21aDYXE/Vys3NRDpR4+RmrmyKVQl7cqsysWKF3is3N0PL2zd3EzO3lpLrbOuQskt6knr5VYg0BnVuxnFYgVbzsyZbYV6MXnMx6m7r9VaXQWXr9MTuWzUwT3n4xagm6x5WIdbRmLlpfXNL3jgx24Ez5kx1oAEAAAAWGDAAAEACMGAAAIAEYMAAAAAJwIABAAASgAEDAAAkAAMGAABIAAYMAACQAAwYAAAgARgwAABAAjBgAACABGDAAAAACcCAAQAAEoABAwAAJAADBgAASAAGDAAAkIB+6gRsps7c885np3eH13pp3a/1+v6wQtNHr3nnty4cDF7U9LOX/FrvPh9cq3/gTe/8w1vHw2vNrvi1NpbDCs1sP+Wd33ywEryomW3+Rm1+cT641tefWPLOf/LpanCt737zsHf+o/uXwgo9vfOQd/7q3bXgRf18z4J3/q31q8G1Xj4w750/feN6cK2zi3Pe+WNrN8MKXT7i/64urt4OXtR7x/Z55/ec/TC41uMgpJTRNYTITu7I/6jf4w7y5+CvKAbOq7aQKlCMXty2QBcVzRo5Oluf+sQlW4PqdT14pJKF2qHep1gvw2c7Tqvrro55pnoVE5Ts20T05KQXI9SZikm7NKLA8jZKqY3Nmfw/u3xmUqifnpDBfDEj9ahBSFYsKr8QgzXNDk2IN5M8sFxTqPvXqG2rQbG5nqqVm5uJdKLGyc1c2RSrEvbkVmVixQq9V25uhpa3b+4mZm4tJdfZ1iFll/Qk9fKrEGkM6tyM47ACreZnTbbCvBi95mLU3dbrrS6Dytbpid23amCe8vCLUU3WPaxCrKMxc9P65pa8cWK2A2fMmepAAwAAACwwYAAAgARgwAAAAAnAgAEAABKAAQMAACQAAwYAAEgABgwAAJAADBgAACABGDAAAEACMGAAAIAEYMAAAAAJwIABAAASgAEDAAAkAAMGAABIAAYMAACQAAwYAAAgAf3UCXiYOnPPmsle2RVL68U7ttZr+2MITR+9Zs1sXTgYqajpZy/ZWu8ejqTVnztvzTy8dTyW1vdXbK2N5RhCM9tPWTOb91+NVNTMNrtdm5+/EUnrySeWrJlPPluNpLXrKfvKffzgYgyhH+08ZM1cvbsWqailPQvWzB/evxJJa/nAvDVz5sb1SFpnF+esmWPv3IwhtHbE/rourt6OVNSdF/ZZM3vPfRBJ65ERUsroGkJkJ3fkf9TvcQf5c/BXFAPnVVtIFShGL25boIuKZo0cna1PfeKSrUH1uh48UslC7VDvU6yX4bMdp9V1V8c8U72KCUr2bSJ6ctKLEepMxaRdGlFgeRul1MbmTP6fXT4zKdRPT8hgvpiRetQgJCsWlV+IwZpmhybEm0keWK4p1P1r1LbVoNhcT9XKzc1EOlHj5GaubIpVCXtyqzKxYoXeKzc3Q8vbN3cTM7eWkuts65CyS3qSevlViDQGdW7GcViBVvOzJlthXoxeczHqbuv1VpdBZev0xO5bNTBPefjFqCbrHlYh1tGYuWl9c0veODHbgTPmTHWgAQAAABYYMAAAQAIwYAAAgARgwAAAAAnAgAEAABKAAQMAACQAAwYAAEgABgwAAJAADBgAACABGDAAAEACMGAAAIAEdGjA0hkAAABERg51HpEusc4MWHg8WFrlJ+wDAAB86UhrP50Z8Ph9AAAACI90niWpbLhbAxaiKbCtEz1tDcYMAACPh+si0hkkMZtuDbhnmXHZANk8/dYrWrwZAAC+fIjJPFM6/iG8Niy6tuF+Rzq5v+aOWz5Fe4Hlgl5v6vS/1Vz2m11RezL14p1G67X98YSmj15T460LB+MJDbR+crHR+vPhqFr9ufNq/PDW8bha33u10frny/GEZr61rMab//ld1KJmtjVN2/z8jahaTz6xpMb/+2w1qtbup36qxh8/uBhPaP7bz6nxX+9ejlrUL/YuqPEf378SVWv5wLwa//bG9aha556ZU+MX3rkZT2jtSPONfebt2wF2FJrBSuPN+vG9arz33AeyJcK24Q4RUkYXFUJkv95R+a7+1Af6T/V0B82zZ+/QFuuV80YZm4tqE9GskSM38apPNPBsUt+ZxylZqB3qfYr1Mny247S67urwU/BWMUHJvk1ET056MUKdqZi0SyMKLG+jlNrYnMn/s8tnJoX66QkZzBczUo8ahGTFovILMVjT7NCEeDPJA8s1hbp/jdq2GhSb66laubmZSCdqnNzMlU2xKmFPblUmVqzQe+XmZmh5++ZuYubWUnKdbR1SdklPUi+/CpHGoM7NOA4r0Gp+1mQrzIvRay5G3W293uoyqGydnuit1pqm3QSpn5GZm3mNVRX6DXSjem19q3U3Tsx24Iw5Ux1oDFDfDv2CWwssM9aXtfZCBE0y1sYAAOBFOmNlBaLlO2z6vbGFqJ8DS2n5mOshwtm5C+Ot6cqAh+Nts+u+unnLaCYp8GAAgK6RPksurbTta9yYQ8s3u7RhYwMzpPwjEn3v+xOu/z975a6jS07E8a4+h3AzIgICIhKE0EJCAIgEEaxARCSgJUQQQIQQAaS8APAIpEQ8AhFE7OPgolwXu9yXmW9u7tk5/99p9fHndpXddk/9ngCN+yM7R/ryskO1TVsTP0CE9ACvQ7EAAHAVfN9z3o7iVqmp/+Co32yDOYVYiymbmSyAtdsG88FyqCeewUQB23vRTn3ZxDaojfHxqdFDkraXMa109v7zrdxrOE7o4IRzPwAAgKej9qNuya0OOZkyXEtNqF0azblpGIeb2azahNBM3OS99TBP24BZAnZZhsTo5ChY93UJg/o9NfZRLWQzYK/eQxnTSXuvYQAAAM9F1hwfaZgPZb0jOdVU42ZwfxOPE/UsKWSeckfWudORvugNQuOdFQ/1ybsGnSfv50C7xpiI7/xWAAAA3ASH7fiRgh61ui/JFFd+xG0EH65mCNisjCbuzmQBH8HpOrbj0Z6lgIfMlQ+TxhOEbAEA4Pkhr69dmqxVl8OedFSIl40ZeN91G7wzDG0tclXRnyVgcZ5dS94DOlXygXrpRJb8GAcfaPhwHZAxAAA8FdMsLepcrau5tr5UnQ31HtfyUSaXlPv3k+YhfdEuYDsN8l3g0Y7S30YO7cizsSbbodI23Bptom1IGtbu3fvss2x0T/0JAAB8uLTSTa46vweqNOJwW4yN2kvEEZUdTGNavXMyCJVdiI2UfDqePKmGrMxlnCgHrZKRlqKLqivk2dV9roDrW+72oApySds5nmFG+tflQKgcaenEwRu/Hjq4nbG5vMv+3MEwMQDgQ+ZAlt7hlT7urJJb0vDW7HrmIfEmbXQx1X8my6i/OplV7sHWycG2sMKuG1uNyX6tgvZH0ioTy/osAa/+vnW3uR1O0rNpuHXvxxhFHbwcafjQqZsGHe3s1sHUP4D9eBgXAAD2ZFkm+6opu5iXXNrdl9Wqo3/Pyy2JAdg9ygJZ+Q6VW80OXTcHS4g6V1xbNMQe2KIooiRboRP1vAyzBFw3ws25tDeUPck7LHvgnm5nSbE3zZ28WBDpBtq1rn3LyKag3ticOZ2YdSNv70wLPPgW5PjidQhmBgC8WUJXZkvKTuX9yC5jplpXa8Vms2yM5laPm7a1AG/0V2VpytTp7Kl5VGtyLcLFh1Ir5T7YGhEiP4vKuDbEJNzn11aNWO3VZjFLwEXfXfxqGnbV6bvXR8mmFNbs6k0NP1jenTmNAk4h6Rj74ENlZgdbLoqPIods4rKkYWIAwKvncY5prmrOyrVzM46SU6PHHazxauWcqJk4XNsyswqkmZiiisuPog5eVahe49X3xQarfdmT1v9WdUhhXtW3VcPNK+yrWCaW8FkCtvORXRENywata3cwtZ0bheomDps2+XE7lsOQEwFvVtKEeuzLdv7xKZ2ZlXZpAQDgtcLPl8oU2lTak/NxjQzzVi9qxaTNqig5eAkt+D3/TGnD3CJUsSmZ77N9xTdemy1QnSCDmbkQiY9FxjK4tHJPz7hD9zNLwHdw//s+yGp0awhcCQAAzwclO/LDq/ZhCN+YijUH3+OUoerzwxb5EqxXLyDo20zj9UCtdvYhnM6HH5kVAADACblk7+v4Bt7F0p3qzT1MvNglfTmympjO1VrHx2M1sN041jPXybMEzONr8clb2k7o0/VPn61//O/urGKniTxPjN9m7p0aQu3b2J/psv72X+Ph0+EnQ/wMh/Pu03++zBbvJvrx3+WaM9f77/7t/Xf+Ommur/95zkTCF778+0kTffRLuebM9cWPfi7XnLm+9qWfzplI+OFXfzJnok+/+ckvvvXJnLn+8P0fyDVnrr/86Ht3PN3okG6r9MvOx8Y/fvbtnJbH0kv7Gkz97uPZ3OsqMPHuijf9+1cf/+fXH++dTybqqvAQwxUOFqe8+GxEVH73FX8/2S251tUbcvl7x73129VCrL1aSGxSW/xZlDVkOttbTodAJma7mtHbMUd/b8RE9agj8zL0H6whP72lcbD+3We1D9kEHoxpGSKPjud94FNXe0OI/5HR8YvcFfugVz5KQuPZvcQLniWhh+7SPS9onypzao89rN+73IuVKj4Jqf3awzmqhpT8p1n3rmXoIYcrkUAbo7Mfj2lp25++ztiXulnbfiW8i7plbePI/rJtwQdri+KxmSXv1X5tw1yH+7ZPMq7t5JVjtRFiu5QXmV/fQ3hsHB3HNlCrQ087rLaXUpu69btWWk7ePa0ln8r4tEHJyiWHxweZP6qhKC5c/O2obOetPSUXfqsl6pv6SAdY2x599ptvTDCj8H7CHBURZz3sVhKKelF61u65ult258HEuV2iCtOqnyX108sNt2k0bLru0XTK8cV4tian5ci+kbOec73xcX0EAIBXDIV3HxvtWFGnlpHCzaH50KdU7ojSiWv5ruVca75XffKEOoJMGlpwrVJbtS1SdrXSqobN9H0Vq+YUqej4pcWqKUzq1bLkU5BbiCSR9OvC2TU8h2kCXt1eHKYUKVb7FtsA9fGqj9oZxshm2WbiFtIMvRyNz48MG58/I/s6soOt3+07fnB030cMBwMA3gTnktYqbU+5KVXrpT1hV2KLz+U811JTbFEH1+fiQZevWdNKOHsp1sZqPaTq1VZo2GfWsW5cFXH9veqKig4pKnt/B/WNyt6VXFzlk5glYNEtraNi7a4N6a97sHSh+tmaYpfuxbrDLBupnRHS7G7+a5bt2TjJex01HHB8O13Gix5e9NDOwG3SNjUAALwhRg03bVkZ74+ajIcqadYkLY5mTQ4vp6QiPzFice2KAk3V1bFsj7THrFnTmDW1INsSqoN9MTq5KkJ7Wv+ysaxpm4uvjtwPvM4t5dMELO9XQqhr7JudBqueqT9dw3aDUJuA/TR0my1baPgOBy/hS5srX5n+aYSDF/NrlbEeGZ0ekM0CEwMAXiuc7k/ORKZnU2PHaqfV3Wh7hSTXtoRslrSKJdTBq8q4lmrVpRnUfakaNp9r1bcs7nkZL8XdDB8Cdg1bHbeh5BquY2S8LV5T1tTNNHOYJuCsWNls8SW52+oe6lb9r1SJdk+HQZfYZNqFmNctVckhy1bD2cGZrYM5vk2OI2vnK2elDvYPbzyo5n4AAHib0NASebkBa03uhfMsmLyuVtdqxczF2DUshVwcrNYuUa19Fvemi9mWIjVahPpObRqVmFpbg1zhySWsVb+uWcKrhqvsQ830JgX8bq17WX3JdnRuvlX7s1aL7cGqxuXQdpjSo2x8+NLa686F1mgT9Uc1RueNGEoHmmVPMXWkJrfxRr8L7AsAeP3QXYoc4Dt9av2msexgTmVYa7P7LyzoVVRLOBfSp+caLlpmVRIUVZuqL9WalLQq/72TvJJTfxYfr2qRbKxRXrh7rdeWeLiuw8YsnnYSswRsorWdli0q4UXbQH+ahVrCnRwiNJuqmFfz8eLmdn9rLI32zWTf+1lz/ExHuTahJvGP1K/GZ45hcDAA4HPC7RpeDkbSVsPjr6Zbv2ujqpRzqNZjdXDYmlo6VlOIhmup1UecHGxZStjdIlzDNcpKf33SNGzWN22LaFuIu0WHWClXu8+r5Pp+E6aBnAAAAHxOmGTGOdMAAAAAIAMBAwAAABcAAQMAAAAXAAEDAAAAFwABAwAAABcAAQMAAAAXAAEDAAAAFwABAwAAABcAAQMAAAAXAAEDAAAAFwABAwAAABcAAQMAAAAXAAEDAAAAFwABAwAAABcAAQMAAAAXAAEDAAAAFwABAwAAABcAAQMAAAAXAAEDAAAAFwABAwAAABcAAQMAAAAXAAEDAAAAFwABAwAAABcAAQMAAAAXAAEDAAAAFwABAwAAABcAAQMAAAAXAAEDAAAAFwABAwAAABcAAQMAAAAXAAEDAAAAFwABAwAAABcAAQMAAAAXAAEDAAAAFwABAwAAABcAAQMAAAAX8FgB/7+9OspSkwegAGw8s5h/r/9Dt9YldCcpkqIg4DhDPRfq9z2giYBX4skFADZQwAAQoIABIEABA0CAAgaAAAUMAAEKGAACFDAABChgAAhQwAAQoIABIEABA0CAAgaAAAUMAAEKGAACFDAABChgAAhQwAAQoIABIEABA0CAAgaAAAUMAAEKGAACFDAABChgAAhQwAAQoIABIEABA0CAAgaAAAUMAAEKGAACjl3ApZR0BABiDl1hH+kAW/3342fXwl0Pl1Ntb05tWGo/OR22E6bD/pzRsD/jerfhnn/uNkxelvx6ebvb6Ta8XH4aXd7udrtkMdttOPot0+E82xNRu6k6ijoaDtkeRu1PmA6H7/161HbCLepk+ETU+5WaZ3sq6urC1U+i3mJMFm64pH4t6uQPM/yppn+S5aijbBdtA1o8Xl43HS+hHtz/mePn3/LCn1BeHv6J44YMZXuGT77lheG7Y3nxoyv//zod2TkdAADekQIGgAAFDAABChgAAhQwAAQoYAAIUMAAEKCAASBAAQNAQKm1pjNsSF9KOgIAMYeusHM6wFYVgLeU7p+tDl/AAHBEChgAAhQwAASUWms6AwC8nZ0WcCmlOy5mm390N9OGzXhyn78UgDtrO/Y39v/FRtiJczrAgvboO+MHt/bR4sl1MF8SAPZsbcd+Zv+/2/wXZ/ZjjwW8/UmV3l+8IQBH0Yr58cwefKQDrFp8Xt3MuFnXZvb8xAFY0/btu139X3VOB1j2oDtr7/EMAAfV2vfbHTyvjweFkrXHAt74sL69bABk1V57k87ych/pAMuuJdrWYFzJ7aPx2tzNdG/m5wBwRNf9f763H323LwfNDQCHpoABIEABA0CAAgaAAAUMAAG/AdNFA1s3xCNPAAAAAElFTkSuQmCC'>"
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig = glucifer.Figure()\n",
    "velmagfield = uw.function.math.sqrt( uw.function.math.dot( velocityField, velocityField ) )\n",
    "fig.append( glucifer.objects.VectorArrows(mesh, velocityField, arrowHead=0.2, scaling=0.9) )\n",
    "# fig.append( glucifer.objects.Mesh(mesh) )\n",
    "fig.append( glucifer.objects.Surface( mesh, pressureField ) )\n",
    "fig.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "ids = mesh.specialSets[\"MaxI_VertexSet\"]\n",
    "coords = mesh.data[ids.data] # xcoords\n",
    "V = velocityField.evaluate(ids)\n",
    "gradV = velocityField.fn_gradient.evaluate(ids)\n",
    "\n",
    "u = V[:,0] ; v = V[:,1]\n",
    "du_dx = gradV[:,0] ; du_dy = gradV[:,1]\n",
    "dv_dx = gradV[:,2] ; dv_dy = gradV[:,3]\n",
    "\n",
    "strainRate = fn.tensor.symmetric( velocityField.fn_gradient )\n",
    "devstress = 2.0 * viscosity * strainRate"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "def exact_vx(y):\n",
    "    ana_u = 1.0 / (2.0 * viscosity)* dp_dx * (y**2 - h*y)\n",
    "    return ana_u\n",
    "\n",
    "def exact_shearSR(y):\n",
    "    shearSR = dp_dx / (2.0*viscosity) * (y - h/2 )\n",
    "    return shearSR"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAtoAAAF3CAYAAACbhOyeAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvNQv5yAAAIABJREFUeJzt3Xt8nHWZ///3NW1aGimhDXVXCzNTpTxok0KxpVtXRaQCheXUQim7ZVsQqIqwVL/7XXGzK4cv0RVd6CooFgSLjRysjYJy8EcUcFEopViaiEqF3KEgW0i74ZBCC3P9/shMmExmkslhzq/n49FHM/fcc88nmcx9X/nM9bkuc3cBAAAAGF2hQg8AAAAAKEcE2gAAAEAOEGgDAAAAOUCgDQAAAOQAgTYAAACQAwTaAAAAQA4QaAMAAAA5QKANAAAA5ACBNgAAAJADBNoAAABADowt9ABGywEHHODRaLTQwwAAAECZe+KJJ15x9ymD7Vc2gXY0GtWmTZsKPQwAAACUOTMLstmP1BEAAAAgBwi0AQAAgBwg0AYAAAByoGxytAEAANDf3r17tX37dr355puFHkrJ2WeffXTggQeqqqpqWI8n0AYAAChj27dv18SJExWNRmVmhR5OyXB3dXZ2avv27Zo2bdqwjkHqCAAAQBl78803VVtbS5A9RGam2traEX0SQKANAABQ5giyh2ekPzcCbQAAACAHCLQBAADQq2lrk6KrowpdEVJ0dVRNW5sKPaSSRaANABWICymAdJq2Nmnl3SsVdAVyuYKuQCvvXjmic8Sll16q66+/vvf25Zdfrq9//etasGCB3F1/+ctfdMghh+ill14ajW+hqFB1BABKWNPWJjW0NKijq0PhmrAaFzRKUp9tJ04/Ufc8c0+f22u3rFX33m5J6r2QPtLxSJ/9sjlW44JGLZu1rGDfP4DR1dDS0HtuSOje262GloZhv9eXLl2qVatW6XOf+5wk6c4779T999+vLVu26Prrr9d9992nK664Qn/913894vEXG3P3Qo9hVMydO9c3bdpU6GEAwKjINoBODpglqSpUJTPTnnf29G6zmCkSRBREAnnIZTK5+p/7E9sT+78w7QXZmL7HSlVdVa0Vh68YNEAnIAcK5+mnn9aMGTOy2jd0RSjj+SF2WWzYY5gxY4ZaWlr08ssv68ILL9QjjzyiXbt2qb6+XvPnz9ePf/zjYR8719L9/MzsCXefO9hjmdEGgCKQHFhPnjBZr+15rTfADboCnfuTc2Vm2rt3b2/QfMOmG/pdEPfG9va5bTHTog2LVN9Wr9a6VjUvbpaH0k+wJILs1P0HSjLs3tutGzbdIMXUO67EWJPHn2nGnOAbKC7hmrCCriDt9pFYsmSJ1q9fr5deeklLly6V1FPfOxQK6X/+538Ui8UUCpVfRnP5fUcAUORS86Mv/PmFfXIiO3d39ptF3hvbq71792rRhkVafutyLdqwSMpicikSRFTfVq+Qh1TfVq9IEBnV/SVJMfUZ19tvv91v/ImAPDnv89yfnKsDrj6APHGgiDQuaFR1VXWfbdVV1b2fVA3X0qVLdfvtt2v9+vVasmSJ3n77bX3qU5/SbbfdphkzZuiaa64Z0fGLFYE2AOTQYEF1x64O3fuje7X7rd2DHms4QXAQCdRa16qYxdRa16og0jNTZepbGzZxO9P+ozGu5Nl3i5mm/nmqdr6xs8+Cqwt/fiGLNIECWjZrmdacvEaRmohMpkhNRGtOXjPiT5/q6ur02muvaerUqXrf+96nr3zlK/rYxz6mj370o7rmmmt000036emnnx6l76J4kKMNAKMkNa86XQ51cn70UNI6Mu2vUN8ANpsc7XR51cljHUqOtsl6Z7RH8n0k9k/NH8+UA07KCZC9oeRooz9ytAGgwBIlsZIreaTLoU6+nToTvHnOZrVPa097/KpQlWyMqXlxszbP2awgEmjC+AlZLUQ8cfqJumfSPVKXFKmJZAxUPxL+SO/jfLbrlgW3pD9WmgA9eVxjx47tH+wnBdADfd+pP69Eyklie2LmWxLBNoCix4w2AAxT8gx2yEJ6x9/pN3s8kOSZ3ba6Nm1YvKH3MVWhKu03fj/t3L2z6Ct5DLVCylBn8hOSf7aT3zNZ+47bt+h+FkAxYkZ7ZEYyo02gDQBZyCotJMsAMjV95JAXDtExpx6je/5c3ukRfSqrjJ+sSc9M0p8P/HPGtJFkg/1s0/1hUm4/P2C4CLRHhtQRAMihbNNC0qVEBNOCgXOOJ4X172f8e0UEhctmLevzfQ4lp32wNJu9sb3q3N0pifQSAMWDQBsAUvQGgLs6NK9znrYduK1fp7R0M6+Jih2JWdcdB+/QZ474DAv5MkgNvKW+eeLJwXfqz3awaijde7t1yc8v0be++y1trN2o8CR+9gDyj0AbAJIkZq93v7W7N1Vhat3UrPKIPeT66ek/1eY5m+WHudYcO/KSWJVmoOC7eXGzXvj4C9p24Db5W4Pnv3+86eN9Xj9muQHkG4E2gIqWmr7w+p7X1b23W9EgOmhFkHSl6Eaj3iz6SpdykpzKk06mVJNL7r2kKBeUAihPNKwBULESAVtyt8LePN9BGrdUV1XrM3M/M+pNHTC41IYatRNqNW7MuD77ZHr9Ond39nm9V969kqY4QKqmJikalUKhnv+beI8MF1VHAFSUdCX5MqGcXOlI98nEzjd2ZlVqsXZCLa8tytqQqo40NUkrV0rdSZ8YVVdLa9ZIy4b3vrj00kt10EEH6XOf+5wk6fLLL9fVV1+tH/7whzrttNMkScuWLdOZZ56pU089dVjPkUuU9xOBNoDBrduyTld96yr9aeqfsqrbnEBKSOnJJr0kVeIPqx0H79CaU3m9UT6GFGhHo1KQZrFxJCK1tw/r+Z988kmtWrVKDz30kCRp5syZ+s53vqNrr71WP/nJT9TV1aXZs2frmWee0dixxZfVPJJAm9QRAGWraWuToqujCl0RUvSaqO5ZcY+W3rxUizYsksUs4+NqJ9SSElLiUtNLIjUR1U6ozbh/ok738luXa+EdC7Xixyt6fm9WR0ktQWXp6Bja9iwcccQR2rFjh1588UVt2bJFkyZN0sc//nE988wzevnll3Xbbbfp9NNPL8oge6TK7zsCAPWf0bSnTNOfmj5ou/Pqqmr91wn/RWBdBoayiDLT4klqcqPihMPpZ7TD4REddsmSJVq/fr1eeuklLV26VJK0fPlyrVu3TrfffrtuueWWER2/WDGjDaAs9Jm9Xh3VJfde0iegGmhx4xgbw+x1BRholnug34/uvd265N5L+vx+McuNstXY2JOTnay6umf7CCxdulS333671q9fryVLlkiSzjnnHK1evVpSTzpJOWJGG0DJS9e5MZWHXM2Lm3u6NSYtjiP/urIMNMud7vcjoXN3J50nURkSCx4bGnrSRcLhniB7mAshE+rq6vTaa69p6tSpet/73idJ+qu/+ivNmDGjd0FkOWIxJICSk67CRCIIGgwVJpBqKJVokvG7hFIxpMWQedTd3a1Zs2Zp8+bNqqmpKfRwMhrJYkhmtAGUlOTOjdmUbktG/jXSSZ7lHkq1ks7dnX1KCDLLDWTvgQce0HnnnafPf/7zRR1kjxSBNoCS0tDS0Kc9emtda8b26Mw4YqgSvx/ZfGKSqFSS/HvY0NLA7xiQhU9+8pMK0i26LDME2gCKXvJH+y7Pqj06s9cYrmyrlaStVBJqV3R1lD/uAEii6giAIpfaJl1KXyGC2tfIlUzVStL9HpqMFu8AerEYEkBRi66Opq0iktwefcL4CQTWyKt0awUUUu8fg8kiNRG1r2rP/yCBuGJdDFkq6AwJoGyk1sNOF2RLPeX6gmmBwpPCBNnIu8Qsd3hSuPf3MF2QLfWUA6T+NlCZyNEGUDTS1cM2GbOEKEqpudwD/mGYlEqSeCyAkXvwwQc1btw4/e3f/m2hh5IWM9oACqbP7PU1UV193dXa/dbuPvu4XCbrs626qlqNC0bWpQwYbY0LGlVdVT3gPrvf2q2rr7ta0WvoMgmMhgcffFC/+c1vCj2MjHIaaJvZQjP7o5ltM7NL09w/3szuiN//mJlFU+4Pm9nrZvbPuRwngPxLXuSomDTnljk67cbTtGjDIlmsb2DtchY6ouilLppMlSgHeNqNp2nOLXOkmFgwiaLlMVf7g+3y2Ois5Vu3bp3mzZun2bNn69Of/rSCIND06dP1yiuvKBaL6WMf+5h+8YtfSJJOO+00zZkzR3V1dVqzZk3vMe677z596EMf0uGHH64FCxaovb1dN9xwg6699lrNnj1bv/71r/s9b1dXlyKRiGKxmCTpjTfe0EEHHaS9e/fqyCOP1IMPPihJ+tKXvqSGhoZR+V6T5Sx1xMzGSLpe0rGStkt63MzucvffJ+12nqRd7n6wmZ0l6WuSlibdf42ke3M1RgCF09DS0JsikrZMWlK5PtJEUCqS00lSU0ky/Z537+2m/jaKisdcG87eoLY721R3Zp0Wr1ssC/X/4zFbTz/9tO644w498sgjqqqq0oUXXqiHHnpIX/ziF/XZz35W8+bN08yZM3XcccdJkm6++WZNnjxZu3fv1pFHHqnTTz9dsVhMF1xwgR5++GFNmzZNO3fu1OTJk/WZz3xG++67r/75n9PPydbU1Gj27Nl66KGH9IlPfEI/+9nPdPzxx6uqqkrf//73dcYZZ+hb3/qW7rvvPj322GPD/h4zyeWM9jxJ29z9WXffI+l2Saem7HOqpLXxr9dLWmBmJklmdpqk5yS15XCMAPJkoEWO6cqkJZAmglKVmkoy0O950BX0eX8ww41CCh4O1HZnm/wdV9udbQoeHlljmZaWFj3xxBM68sgjNXv2bLW0tOjZZ5/V+eefr1dffVU33HCDvvGNb/Tu/81vflOHH3645s+fr+eff17PPPOMHn30UR111FGaNm2aJGny5MlZP//SpUt1xx13SJJuv/12LV3aM6dbV1enf/zHf9RJJ52km2++WePGjRvR95lOLhdDTpX0fNLt7ZL+JtM+7v62mXVJqjWzNyV9UT2z4aSNACVusEWOHnI1L27W5jmb9dqhrym8T5iGHyh5/bpMTgrroWUPafMfNiuIBH26mSbqb0ti0SQKLnJURHVn1vXOaEeOiozoeO6uFStW6Ktf/Wqf7d3d3dq+fbsk6fXXX9fEiRP14IMP6oEHHtBvf/tbVVdX6+ijj9abb745ouc/5ZRT9K//+q/auXOnnnjiCR1zzDG9923dulX777+/duzYMaLnyKRYq45cLulad389PsGdlpmtlLRSksLhcH5GBmDIktNEEhKLHJOD7R2H7NCavyP/GuUjbZfJ2Er53r5BdmplHdJJUEgWMi1et1hzVs5R5KjIiNJGJGnBggU69dRT9fnPf17vfe97tXPnTr322mv6xje+oWXLlikSieiCCy7Qz372M3V1dWnSpEmqrq7WH/7wBz366KOSpPnz5+vCCy/Uc8891yd1ZOLEiXr11VcHfP59991XRx55pC655BKddNJJGjNmjCRpw4YN2rlzpx5++GGddNJJ2rhxo/bff/8Rfa+pcpk68oKkg5JuHxjflnYfMxsrqUZSp3pmvq82s3ZJqyT9q5ldlPoE7r7G3ee6+9wpU6aM/ncAYFR0dHWk3c4iR1SadF0mM9XfzvS+AfLBQqbo0dERB9mSNHPmTF111VU67rjjdNhhh+nYY49Ve3u7Hn/8cX3xi1/UsmXLNG7cON1yyy1auHCh3n77bc2YMUOXXnqp5s+fL0maMmWK1qxZo8WLF+vwww/vTf84+eST1dzcnHExZMLSpUu1bt263se98soruvTSS3XTTTfpkEMO0UUXXaRLLrlkxN9rqpx1howHzn+StEA9AfXjkv7B3duS9vmcpFnu/pn4YsjF7n5mynEul/S6u39DA6AzJFBcmrY29X5kHrKQ3vF3+u3DIkcgc/3t2gm12nfcvqRRYcToDDkyRdkZ0t3flnSRpPslPS3pTndvM7MrzeyU+G7fU09O9jZJX5DUrwQggNLTtLVJK3+6UvY7k2JKG2SzyBHoka7+dlWoSq/tea23/KX9zrTyp5QBBEpNTnO03f0eSfekbPty0tdvSloyyDEuz8ngAORMw//XoIV3LFR9W71a61rVvLhZHnKNsTGKeYzZOSBJv0WTNWG9vud1de7u7K29nXgvNVSTtw1k0tjYqB/96Ed9ti1ZsiQn9bGzlbPUkXwjdQQonOQ0kXBNWPY70/JblyvkIcUspluX36r2ae0ymWKXxQo9XKDoha4IyeWKPhft917y2U46CYaE1JGRKcrUEQCVIbnDo8sVdAUZ6wWHa6gOBGQj8V5J915Kfq/RVRIobsVa3g9AiUhbui+pLnaiXjA52UD2Ghc09taeT30vJaMMILLl7hqoZDLSG2nmBzPaAIZkoA6PyTzk8tkuhUTpPmCIkssAKiT5bO8XZCcEXQEdJTGgffbZR52dnSMOGiuNu6uzs1P77LPPsI9BjjaArKV2eJTSN9uQKN0HjLaB/rBNqK6q5o9a9LN3715t3759xB0WK9E+++yjAw88UFVVVX22Z5ujTeoIgKxl0+FRonQfkAvJ6SSZkEqCdKqqqjRt2rRCD6MikToCIGt0eAQKJ7WrZCZ0lASKB6kjADLqLdu3q0PzOudp24Hb1PlWZ7/9SBMB8i9jR8nxtTp4+8HaWLtR4UmUAARygfJ+AEYkkY/dsatDizYs0vHfPl6f+OEnNM7G9dmPNBGgMNJ2lFSVjvnhMTr+28dr0YZF6tjVQQlAoIAItAGklcjHjgQR1bfVK+QhzWydqUNfOJQ0EaAIpKaSRGoimvHiDM1onaGQh1TfVq9IEOnN2waQfyyGBCCpf3fHxEfSiYYZiRbQW/96q2Kr6O4IFINls5b1+UM3dFlIB9cd3Pt+TTSLCroCRVdH6SgJ5Bk52gAGLdtnMVMkiCiIBApPCpOPDRSp6OqoOnZ19L5fE7W301UG4tMoYPjI0QaQtYHK9kk9zWfap7VrwvgJ5GMDRaxxQaMmjJ+g9mntGYNsSaSTAHlCoA2Asn1AmUiXt52uoZREGUAgH0gdASpUck52yEJ6x9/ptw9l+4DSl7EM4IRa7TtuX/K2gWEgdQRARomc7KArkMvTBtmU7QPKQ9oygKEqvbbntd5zQNAVUAYQyAECbaACpcvJlqQxNoY0EaDMpEsn2W/8ftrzzp4++5G3DYw+yvsBFahjV4eiQbRPVQJJinlMscso3QeUm35lAK/oO8+WXFkIwOgh0AYqQJ8a2RPDOuunZ2n6U9PVWteq5sXNvcF2uCZc4JECyIfkWvkWMy3asEj1bfX6ff3vFVKI1u3AKCF1BChzqfnY9pRp+lPT+3SOk8jJBipJct52avfXSBAhZxsYJQTaQJlLzcdOdHqMWUzPHfGcgkhATjZQYZLztoNIoN/X/14xi/XpJknONjBylPcDylzoilC/OrqJfMxnb35WFrICjQxAsQhdFurXTVLqaXbDug2gP8r7ARWqaWuToqujCl0RUnR1VJMnTO63j4dcPtsJsgFIksKTwn26SSZMnjC5z/mEVBJgaAi0gTKSmo8ddAV69a1XNW7MuD77kY8NIBm1toHcINAGyki6+th7Y3s1cdxEWqkDyIha20BuUN4PKCMdXR1pt+/cvVOv/MsreR4NgFIyWK3thEznGQD9MaMNlJFMdbCpjw1gqDifACNHoA2UsNSFjydOP7FfniX52ACGI13ednVVtU6cfiILJIEsEWgDJSqx8LFjV4ciz0XUsatDa7es1YrDV5CPDWDE0uVtrzh8hdZuWdvnvMMCSSAz6mgDJSq6OqqOXR29rZMT7dTDk8JqX9Ve6OEBKEOcd4Ae1NEGylxHV0ef1smJduosVAKQK5x3gKEh0AZKSHJOdshCfdqpJ1ons1AJQK6Ea8Jpzzs0tgHSI3UEKBGJnOzUOtmJdupBJNCE8RPIyQaQM4nz0O63dveed8aOHSsz61Nzu7qqmnMRyhqpI0CZSdeMRpJCY0IKpgUKTwpzYQOQU4kFkuFJ4d7zDo1tgMxoWAOUiEw5kDGPKXZZLM+jAVCpaGwDZI8ZbaBE0DwCQDHi3ARkRqANFCma0QAoBeka21SFqvT6ntdZHImKR6ANFKHEgqOgK5DLFXQFNKMBUJRSG9vUTqiVmalzd2fv+YumNqhUVB0BilB0dVRBV9Bve6QmQlMIAEWN8xcqAVVHgBKWaRERi4sAFDvOX8C7CLSBItAnH/uaqGa9OEsWs377sbgIQLFLd56ymGnWi7MUvYamNqgsBNpAgSXnYysmzblljk678TSd0XxGn2CbhY8ASkHq4kiLmU5vPl2n3Xia5twyR4qJvG1UDAJtoMCSG9FEgojq2+oV8pBmtM7QvM55LHwEUFJSF0fO65ynma0zFfKQ6tvqFQkikmhqg8pAwxqgwJLzFoNIoNa6VtW31autrk2//eZvZaH+KSQAUMySm9p4zLXkoSWqa6tTa12rgsi7CyXJ20a5I9AGCixcE+5doe8hV/PiZm2es1l+mBNkAyh5FjJtOneTnnjqCQWRQB56t9oZ605Q7kgdAfJssEY0HnLtOGSHGo8lHxtAeWg8tlE7DtnRJ8iurqrWidNP7HM+JGcb5YZAG8gjGtEAqESpeduRmohWHL5Ca7es7XM+ZIEkyg0Na4A8opEDAPTgfIhSRsMaoAjRyAEAenA+RCUg0AbyKNPCHxYEAag0nA9RCQi0gRxLXvz4+p7XNW7MuD7304gGQCVKbWwjsUAS5YfyfkAOJRY/JhrSdO7uVFWoSrUTarVz906Fa8JqXNDIwkcAFSdx3mtoaVBHV4fCNWGdOP1Erd2ytvecmVggmbw/UEpYDAnkUPSaqOwp61c7lsU+ANBf8gJJi5kiQURBJFB4UphzJopKtoshmdEGcsRjrrm3zO3thta8uLk32GaxDwD0lzg3Wsy0aMMi1bfV954/gVJEjjaQI8HDgera6hTykOrb6hUJIr33sdgHAPpLnBsjQUT1bfW95895nfMKPDJgeAi0gVGUvPDx6M1H662j3lLMYmqta1UQ6fk4lMWPAJBeYoFkEAnUWteqmMX0+/rfa9uB21gciZJEjjYwSlIXPkpS9ZhqnT/mfP10/E/V8VoHix8BYBBNW5t6Fkju6tCsl2bpD1P/oD2+p/f+6qpquuei4LLN0SbQBkYJXc4AYHRxXkWxKorOkGa20Mz+aGbbzOzSNPePN7M74vc/ZmbR+PZ5Zva7+L8tZrYol+MERgNdzgBgdHFeRanLWaBtZmMkXS/pBEkzJf29mc1M2e08Sbvc/WBJ10r6Wnx7q6S57j5b0kJJ3zUzKqSgqNHlDABGF+dVlLpczmjPk7TN3Z919z2Sbpd0aso+p0paG/96vaQFZmbu3u3ub8e37yOpPPJbUNYydTlj4SMADA/nVZS6XAbaUyU9n3R7e3xb2n3igXWXpFpJMrO/MbM2SVslfSYp8AaKQnKFkejqqCRpzclrFKmJyGSK1ERYsAMAI7Bs1rK051VJtGlHSSjadAx3f0xSnZnNkLTWzO519zeT9zGzlZJWSlI4zMdIyJ/UCiOJNsFrTl7DAh0AGEXLZi3rM2GR6fyb2BcoJrmc0X5B0kFJtw+Mb0u7TzwHu0ZSZ/IO7v60pNcl1ac+gbuvcfe57j53ypQpozh0YGANLQ19yvhJUvfebjW0NBRoRABQGTj/opTkMtB+XNJ0M5tmZuMknSXprpR97pK0Iv71GZJ+6e4ef8xYSTKziKRDJbXncKzAkLASHgAKg/MvSknOAu14TvVFku6X9LSkO929zcyuNLNT4rt9T1KtmW2T9AVJiRKAH5W0xcx+J6lZ0oXu/kquxgoMFSvhAaAwOP+ilOS0jra73+Puh7j7B929Mb7ty+5+V/zrN919ibsf7O7z3P3Z+PYfuHudu8929w+5+09yOU5gML0LHy8Laf5F83XiB09kJTwAFECmSiQnfvBEzb9ovkKXsUASxSOngTZQDhILbzp2dWjRhkU6/tvHq/PfOrVi1goqjABAnqWrRLJi1gp1/lunjv/28Vq0YZE6dnVo5d0rCbZRcLRgBwaRaAEcfS6q5bcuV8hDillM9194vx697tFCDw8AKt78i+br+G8f33t+vnX5rWqf1k6rduRMUbRgB8pBYoFNEAnUWteqmMXUWteqjbUbCzwyAIAkbazd2Of8HEQCSSyQROEVbR1toFiEa8IKugJ5yNW8uFmb52xWEAkUnsTCGwAoBuFJ4T7nZw/1fFrPAkkUGjPawCCSF954yNU+rV0Txk9g4SMAFInGBY2aMH6C2qe19wbZLFBHMSDQBtJIbq/e0NKgFYez8BEAihWt2lGsWAwJpEht7yv1zIwQXANA6eBcjlxiMSQwTLT3BYDSx7kcxYBAG0hBe18AKH2cy1EMCLSBFLT3BYDSx7kcxYBAGxUveeFjdHVUJ06nvToAlLp0rdqrQlV6fc/rLI5E3hBoo6IlFssEXYFcrqAr0Nota6kyAgAlLrUSSe2EWpmZOnd39p7vadOOXKPqCCpaor16Ktr2AkB54XyP0UTVESALLJYBgMrA+R6FQKCNihaeGFb0uagsZn23s1gGAMpKuvO6xUzzXp4nj5XHp/soPgTaqFgec138wMVafutyLdqwqDfYZuEjAJSf1MWRFjMt+ckSnfCdE7Th7A0E28gJAm1UlOQKIx/+pw/r9ftfV8hDqm+rVySIsPARAMpU6uLIeZ3zNLN1phSTnrr9KX3gUx+gEglG3dhCDwDIl9R2vBtrNypSH9HM1pk6bOlhevbmZ2UhG+QoAIBStWzWst6JFI+5rnn2Gr1636tqrWtVEAnkXa6Vd6/s3RcYKaqOoGKkW3FusZ5Zjd9+87cE2QBQYaLXRGVPWU+QHXo3HqISCQaTbdURZrRRMdKtLPeQa+OUjQTZAFCBOl7rkE/rP+FIJRKMFnK0UTFoxwsASMZ1AblGoI2Kka4dLxVGAKBycV1ArhFoo2wlVxiJro5KUp8V51QYAYDKllqJJHFdkNTn+kElEgwXiyFRllIrjEg9sxQE1gCAgXD9QDZowY6K1tDS0OckKUnde7vV0NJQoBEBAErz5wc4AAAgAElEQVQB1w+MJgJtlKVMK8ZZSQ4AGAjXD4wmAm2UJVaSAwCGg+sHRhOBNsoSK8kBAMPB9QOjadBA28z+08zq8jEYYCSSq4w0tDRoxeErqDACABgSKpFgNA1adcTMzpd0rnq6SN4i6TZ378rD2IaEqiOVjVXiAIBc4RqDVKNWdcTdb3L3j0haLikq6Skz+6GZfWLkwwRGR0NLg3a/tVvR56KyWE87dVaJAwBGQ6ISicWs9zrDNQbZyCpH28zGSDo0/u8VSVskfcHMbs/h2ICsdezq0KINi7T81uVatGFRb7DNKnEAwEh1dHXIYtbvOsM1BoPJJkf7Wkl/kHSipK+4+xx3/5q7nyzpiFwPEMjGvM55qm+rV8hDqm+rVySISGKVOABg5MI1YUWCSL/rDNcYDCabGe2nJM1290+7+8aU++blYEzAkF10wUX6w6w/KGYxtda1KogErBIHAIyKxgWN2nHwDrXWtfZeZ3YcvINrDAY1Not9znb3W5I3mFmLuy8oxkWRqExnH362dKt03Y3XaWPtRoUnhdW4oJFFKgCAEUtcSxqqG7T5qc3yw1xrjmUhJAaXseqIme0jqVrSryQdLcnid+0n6T53PzQfA8wWVUcqS9PWJjW0NKijq0PhGoJqAEB+cR2qbNlWHRloRvvTklZJer+kzUnbX5V03ciGBwxfapmloCvQyrtXShInOQBAznEdQrayqaN9sbt/K0/jGTZmtCtHdHVUQVfQb3ukJqL2Ve35HxAAoKJwHcKIZ7TN7Bh3/6WkF8xscer97r5hhGMEhiVTOSXKLAEA8oHrELI1UOrIxyX9UtLJae5zSQTaKIhwTTjtTAJllgAA+cB1CNnKWN7P3S+L/39umn+fyt8Qgb4aFzSquqq6zzZK+QEA8oXrELKVTcOar5jZ/km3J5nZVbkdFvCupq1Niq6OKnRFSNHVUUnSmpPXKFITkckUqYlozcmUWQIA5MeyWcvSXock9bleNW1tKvBIUWjZLIZ80t2PSNm22d0/lNORDRGLIctT6spuqWfWgMAaAFBMuF5VlmwXQ2bTGXKMmY1POvAESeMH2B8YNQ0tDX1OWpLUvbdbDS0NBRoRAAD9cb1COtl0hmyS1GJmie6Q50pam7shAe9iZTcAoBRwvUI6g85ou/vXJF0laUb83/9z96tzPTBAyryCm5XdAIBiwvUK6WSTOiJJT0p6SNKD8a+BvGBlNwCgFHC9QjrZVB05U9JGSWdIOlPSY2Z2Rq4HBnjM9ZHOj+i7f/ddKowAAIpapkok/1D3D2p/sF0eG7j4BMpTNlVHtkg61t13xG9PkfSAux+eh/Fljaoj5aFpa5MaWhrUsatDZ//sbH3wdx9U/dJ6LV63WBayQg8PAICsecx1zUnX6NX7XlVbXZs2nbtJjcc2MllUBkaz6kgoEWTHdWb5OGBIEqWRgq5AkSCiaU9Ok2JS6x2tCh7u34ELAIBidtP3b9Kr972qkIdU11Yne8q08u6V1NeuINkEzPeZ2f1mdo6ZnSPp55Luye2wUImSSyMFkUCtda2KWUx/nv1nRY6KFHh0AAAMTeP/NvZey1rrWhVEAkr+VZhBU0ckycxOl/SR+M1fu3tzTkc1DKSOlL7QFSG53v19tJgpEkQURALFrogVcGQAAAxd6IqQFFPvtcxDPdc4kyl2Gde1UpZt6kg2dbTl7j+W9OMRjwoYQLgmrKDr3RQRD7nap7UrUsNsNgCg9CSua+3T2vttR2XImDpiZq+Z2atp/r1mZq/mc5CoDJRGAgCUE65ryBhou/tEd98vzb+J7r5fPgeJ8tW0tUnR1VGFrgipoaVBKw5fQSk/AEBZyFTyT1LvtS+6OsriyDKWbY72RyVNd/dbzOwASRPd/bmcj24IyNEuPYkqI4kFkFLPX/oE1wCAcsW1rzyMWnk/M7tM0hclfSm+aZykdSMbHtC3ykgCq7EBAOWMa19lyaa83yJJp0h6Q5Lc/UVJE7M5uJktNLM/mtk2M7s0zf3jzeyO+P2PmVk0vv1YM3vCzLbG/z8m228IpaOjq2NI2wEAKHVc+ypLNoH2Hu/JL3FJMrP3ZHNgMxsj6XpJJ0iaKenvzWxmym7nSdrl7gdLulbS1+LbX5F0srvPkrRC0g+yeU6UlkyrrlmNDQAoV1z7Kks2gfadZvZdSfub2QWSHpB0YxaPmydpm7s/6+57JN0u6dSUfU6VtDb+9XpJC8zM3P3J+My5JLVJmmBm47N4TpQQVmMDACoN177Kkk2g7ZL+Wz11tA+R9GV3/1YWj5sq6fmk29vj29Lu4+5vS+qSVJuyz+mSNrv7W1k8J0pIptXYLAYBAJQrrn2VJZuGNftK+pSknZLukPRUTkeUxMzq1JNOclyG+1dKWilJ4TAfuRS7pq1NamhpUEdXh8I1YTUuaNSyWcs4uQAAKkrqtS9R6jb1+ojSN+iMtrtf4e51kj4n6X2SHjKzB7I49guSDkq6fWB8W9p9zGyspBpJnfHbB0pqlrTc3f+cYWxr3H2uu8+dMmVKFkNCoSTKGQVdgVyuoCvQyrtXUjsUAFDRuD6Wt2xSRxJ2SHpJPYHwe7PY/3FJ081smpmNk3SWpLtS9rlLPYsdJekMSb90dzez/SX9XNKl7v7IEMaIIkU5IwAA+uP6WN6yqaN9oZk9KKlFPfnTF7j7YYM9Lp5zfZGk+yU9LelOd28zsyvN7JT4bt+TVGtm2yR9QVKiBOBFkg6W9GUz+138XzbBPYoU5YwAAOiP62N5yyZH+yBJq9z9d0M9uLvfI+melG1fTvr6TUlL0jzuKklXDfX5ULzCE8Oyp0xBJJCH3u1GSjkjAEAlC9eEFXQFvbctZooEEflhg3fuRvEbNNB29y8Ntg8wEI+5Ln7gYr1636tqrWtV8+JmecgpZwQAqHiNCxp7W7JbzLRowyLVt9Vrv4X7yVe5LGSFHiJGYCg52sCwBA8HeuMXbyjkIdW31SsSRChnBACA+pb7iwQR1bfVK+QhvfGLNxQ8HAx+ABQ1Am3kRKJUUeiKkI7efLTec9x7ZGNMh511mJ69+Vm1r2onyAYAQD3Bdvuqdj1787M67KzDZGNM7znuPTp689EKXRFSdHWUKiQlynq6q5e+uXPn+qZNmwo9DOjdUkXJq6irx1Rr9ftX6/xzzudjMAAAMvCY66bv36RVL65S9ztJ19Gqaj4JLiJm9oS7zx1sP2a0MerSlip6p1uNrzYSZAMAMAALmRpfbewTZEuU/CtVBNoYdZQqAgBg+LiOlg8CbYy6TCX7KOUHAMDguI6WDwJtjLrGBY2qrqrus41SfgAAZIfraPkg0MaoSy5VZDJK+QEAMARcR8sHVUcwYk1bm9TQ0qCOrg6Fa8JqXNDIyQAAgFHG9bZ4ZFt1JJsW7EBGqaX8gq5AK+9eKUm8+QEAGCVcb0sTqSMYkbSl/ChBBADAqOJ6W5oItDEilCACACD3uN6WJgJtjAgliAAAyD2ut6WJQBsjQgkiAAByj+ttaSLQxohQgggAgNzjeluaKO8HAAAADEG25f2Y0caQNW1tUvSaqKadM03Ra6Jq2tpU6CEBAFCR1m1Zp/kXzVfospCiq7kmFxvqaGNImrY2aeVPV2rhHQtV31av1rpWreymjicAAPm2bss63b38bh2/9XhNrZuq5sXN1NYuMsxoY0gaWhr03m3vVX1bvUIeUn1bvd677b3U8QQAIM+uu/E6Hbr10N7rcSSIUFu7yBBoY0g6ujoURAK11rUqZjG11rUqiATU8QQAIM821m7sdz2WqK1dTEgdwZCEa8IKugI1L27W5jmbFUQCecgVqYkUemgAAFSU8KRwv+uxRG3tYsKMNoYkUcfTQ672ae3ykFPHEwCAAmhc0KgJ4yf0Xo8lamsXGwJtDAl1PAEAKA5ck4sfdbQxoKatTWpoaVBHV4fCNWE1LmjkDQwAQJHiup0f2dbRJkcbGTVtbdLKu1eqe2+3JCnoCigbBABAkeK6XXxIHUFGDS0NvW/WBMoGAQBQnLhuFx8CbWSUqTwQZYMAACg+XLeLD4E2MspUHoiyQQAAFB+u28WHQBsZJUr5JaNsEAAAxYnrdvEh0EZGlA0CAKB0cN0uPpT3AwAAAIYg2/J+zGijV9PWJkVXRxW6IqTo6qiatjYVekgAAGCEuL4XDnW0IYnamwAAlCOu74XFjDYkUXsTAIByxPW9sAi0IYnamwAAlCOu74VFoA1J1N4EAKAccX0vLAJtyGOuhv0aVD2G2psAAJSTtLW1x1SrYb8Geaw8Ks8VMwLtCucx14azN+gvK/+iK39zpSITqb0JAEC56Fdbe2JEV/7mSv1l5V+04ewNBNs5Rh3tCtf+YLtu/eSt8ndcNsa0/IHlih4dLfSwAABADnDdHx3U0caAEjU1P/CrD2jb4dukkFR3Zp0iR0UKPTQAAJAjkaMiqjuzTgpJ2w7fpg/86gPU1s4h6mhXoD41NUPSupPW6ZAjDtG/XfxvspAVengAACBHLGTa/cXduq36Nv1p6p/kIae2dg6ROlKBoqujCrqCftsjNRG1r2rP/4AAAEDeEAeMHKkjyIiamgAAVC7igPwh0K5A1NQEAKByEQfkD4F2BUpbU5Oa2QAAVATigPwh0K5A/WpqUjMbAICKQRyQP1QdqRBNW5vU0NKgjq4OhWvCalzQyIIHAAAq1LJZy/oF1uliBYLvkSHQrgB9yvlJlPEBAAB9ECvkBqkjFaChpaH3jZPQvbdbDS0NBRoRAAAoJsQKuUGgXQEo4wMAAAZCrJAbBNoVgDI+AABgIMQKuUGgXQEo4wMAAAZCrJAbBNoVgDI+AABgIMQKuWHuXugxjIq5c+f6pk2bCj0MAAAAlDkze8Ld5w62HzPaZahpa5Oiq6MKXRFSdHVUTVubCj0kAABQYognRo462mWGOpgAAGCkiCdGBzPaZYY6mAAAYKSIJ0ZHTgNtM1toZn80s21mdmma+8eb2R3x+x8zs2h8e62Z/crMXjez63I5xnLT0dUhi5miz0VlMeuzHQAAIBvJcUNyXEE8MTQ5Sx0xszGSrpd0rKTtkh43s7vc/fdJu50naZe7H2xmZ0n6mqSlkt6U9O+S6uP/kKXwxLDm3DJH9W31aq1rVfPiZnnIqYMJAACyFq4JK+gKZDHTog2LeuOKJ859otBDKym5nNGeJ2mbuz/r7nsk3S7p1JR9TpW0Nv71ekkLzMzc/Q13/2/1BNwYgob9G1TfVq+Qh1TfVq9IEKEOJgAAGJJEXe1IEOkTVzTsT+rIUOQy0J4q6fmk29vj29Lu4+5vS+qSVJvDMZW98885X/st3E8xi6mtrk1+mFMHEwAADEmirrYf5mqra1PMYtpv4X46/5zzCz20klLSVUfMbKWklZIUDpMaIUkWMn3hZ19Q8HCgyFERWcgGfxAAAECKZbOWadmsZfJVTlwxTLmc0X5B0kFJtw+Mb0u7j5mNlVQjqTPbJ3D3Ne4+193nTpkyZYTDLR8WMkWPjvJmAAAAI0ZcMXy5DLQflzTdzKaZ2ThJZ0m6K2WfuyStiH99hqRferm0qswTiskDAIB8Ie4Ympyljrj722Z2kaT7JY2RdLO7t5nZlZI2uftdkr4n6Qdmtk3STvUE45IkM2uXtJ+kcWZ2mqTjUiqWVDyKyQMAgHwh7hg6K5cJ5Llz5/qmTZsKPYy8iq6OKugK+m2P1ETUvqo9/wMCAABli7jjXWb2hLvPHWw/OkOWsExF4ykmDwAARhtxx9ARaJewTE1oaE4DAABGG3HH0BFol7BEMflkNKcBAAC5QNwxdATaJSxRTD5SE5HJFKmJ0JwGAADkBHHH0LEYEgAAABgCFkOWIWpXAgCAYkN8kllJt2CvJNSuBAAAxYb4ZGDMaJeIhpaG3l/ihO693WpoaSjQiAAAQKUjPhkYgXaJoHYlAAAoNsQnAyPQLhHUrgQAAMWG+GRgBNolgtqVAACg2BCfDIxAu0RQuxIAABQb4pOBUUcbAAAAGALqaJeB3rqUl4U0/6L5WrdlXaGHBAAAMKh1W9Zp/kXzFbqssmtrU0e7SCXqUu5+a7cWbVik+rZ63f3ru6VbpbMPP7vQwwMAAEhr3ZZ1unv53Tp+6/GaWjdVzYubK7a2NjPaRSpRlzISRFTfVq+Qh3To1kN13Y3XFXpoAAAAGV1343U6dOuhCnlI9W31igSRiq2tTaBdpBL1J4NIoNa6VsUspta6Vm2s3VjgkQEAAGS2sXZjn9gliASSKrO2NqkjRSpcE1bQFchDrubFzdo8Z7OCSKDwJOpSAgCA4hWeFO4Tu3iop/BGJdbWZka7SCXXpfSQq31auyaMn0BdSgAAUNQaFzRqwvgJap/W3htkV2ptbQLtIkVdSgAAUIqIYd5FHW0AAABgCKijDQAAABQQgXaR6G1Oc0VlF3YHAADlp1LjHKqOFIFEc5ruvd2SpKArqNjC7gAAoLxUcpzDjHYRSDSnSVaphd0BAEB5qeQ4h0C7CGQq4F6Jhd0BAEB5qeQ4h0C7CGQq4F6Jhd0BAEB5qeQ4h0C7CCQ3p0mo1MLuAACgvFRynEOgXQQo7A4AAMpVJcc5NKwBAAAAhoCGNQAAAEABEWgXQKUWbQcAAEiohHiIhjV5VslF2wEAAKTKiYeY0c6zSi7aDgAAIFVOPESgnWeVXLQdAABAqpx4iEA7z8ITw4o+F5XFrO/2CijaDgAAIPWPeyxmij4XVXhiecVDBNp55DHXxQ9crOW3LteiDYt6g+1KKdoOAAAg9W1iYzHTog2LtPzW5br4gYvlsfIoPS0RaOdV8HCgN37xhkIeUn1bvSJBpKKKtgMAAEh9m9hEgojq2+oV8pDe+MUbCh4OCj28UUPVkTyKHBVR3Zl1aruzTYedeZievflZWcgGfyAAAECZWTZrmZbNWiaPuTbs2aC2O9tUd2adIkdFCj20UUNnyDzzmCt4OFDkqAhBNgAAgEovPqIzZJFILcb+w7YfKnp0tCR+iQAAAPLBQqbo0fE4qYya2JA6kkOVUowdAABgpMoxbmJGO4cqpRg7AADASJVj3ESgnUOVUowdAABgpMoxbiLQzqFMTWhoTgMAANBXOcZNBNo5lFyMPYHmNAAAAP2VY9xEoJ1DycXYTUZzGgAAgAzKMW6ijjYAAAAwBNTRBgAAAAqIQHuUpTaoKfVC6wAAAIVUyrEVDWtGUTkWWgcAACiUUo+tmNEeReVYaB0AAKBQSj22ItAeReVYaB0AAKBQSj22ItAeReVYaB0AAKBQSj22ItAeReVYaB0AAKBQSj22ItAeReVYaB0AAKBQSj22omENAAAAMARF0bDGzBaa2R/NbJuZXZrm/vFmdkf8/sfMLJp035fi2/9oZsfncpwAAADAaMtZoG1mYyRdL+kESTMl/b2ZzUzZ7TxJu9z9YEnXSvpa/LEzJZ0lqU7SQknfjh+vqPz3f1yo7ZPHKmam7ZPH6r//40KpqUmKRqVQqOf/Cy/se7upqf8+2W7jWByrVI5VSmPlWByrVI5VSmPlWBwrh8f6w9JP9o+/ilTOUkfM7MOSLnf34+O3vyRJ7v7VpH3uj+/zWzMbK+klSVMkXZq8b/J+mZ4v36kj//0fF+qIL39H1XtNgSKKKNBbIVfVmLEas/ftzA+sqpLMpD17hr6NY3GsUjlWKY2VY3GsUjlWKY2VY3GsHB6rJ3J9N/7qrnI9eeVn9dFLv535+UZZtqkjuQy0z5C00N3Pj9/+R0l/4+4XJe3TGt9ne/z2nyX9jaTLJT3q7uvi278n6V53X5/p+fIdaG+fPFZTd8W0QYvUpnrVqVWL1SxTeeS8AwAAFCOX9Yu/XpgU0oE7B5joHGVFkaOda2a20sw2mdmml19+Oa/P/f5d7yhQRG2qlyukNtUrUCSvYwAAAKg06eKv9+96p9DDSiuXgfYLkg5Kun1gfFvafeKpIzWSOrN8rNx9jbvPdfe5U6ZMGcWhD+7FSWMUUaA6tcoUU51aFVGQ1zEAAABUmnTx14uTim4pn6TcBtqPS5puZtPMbJx6FjfelbLPXZJWxL8+Q9IvvSeX5S5JZ8WrkkyTNF3SxhyOdcja/2Wluqtci9Ws5bpVi9Wst0Kud6rGDvzAqipp3LjhbeNYHKtUjlVKY+VYHKtUjlVKY+VYHCunx+obf3VXudr/ZeXAz1co7p6zf5JOlPQnSX+W1BDfdqWkU+Jf7yPpR5K2qSeQ/kDSYxvij/ujpBMGe645c+Z4vv36q5/15yeN8Xckf37SGP/1Vz/rvm6deyTibtbz/2c/2/f2unX998l2G8fiWKVyrFIaK8fiWKVyrFIaK8fiWDk81tNnLugff+WZpE2eRSycs8WQ+UbDGgAAAORDRSyGBAAAAIoVgTYAAACQAwTaAAAAQA4QaAMAAAA5QKANAAAA5ACBNgAAAJADBNoAAABADhBoAwAAADlAoA0AAADkAIE2AAAAkANl04LdzF6WFBTo6Q+Q9EqBnhv5w+tc/niNKwOvc2XgdS5/hXyNI+4+ZbCdyibQLiQz25RNv3uUNl7n8sdrXBl4nSsDr3P5K4XXmNQRAAAAIAcItAEAAIAcINAeHWsKPQDkBa9z+eM1rgy8zpWB17n8Ff1rTI42AAAAkAPMaAMAAAA5QKA9BGa20Mz+aGbbzOzSNPePN7M74vc/ZmbR/I8SI5HFa/wFM/u9mT1lZi1mFinEODEyg73OSfudbmZuZkW9qh3pZfM6m9mZ8fd0m5n9MN9jxMhkcc4Om9mvzOzJ+Hn7xEKMEyNjZjeb2Q4za81wv5nZN+O/B0+Z2YfyPcZMCLSzZGZjJF0v6QRJMyX9vZnNTNntPEm73P1gSddK+lp+R4mRyPI1flLSXHc/TNJ6SVfnd5QYqSxfZ5nZREmXSHosvyPEaMjmdTaz6ZK+JOkj7l4naVXeB4phy/K9/G+S7nT3IySdJenb+R0lRsn3JS0c4P4TJE2P/1sp6Tt5GFNWCLSzN0/SNnd/1t33SLpd0qkp+5wqaW386/WSFpiZ5XGMGJlBX2N3/5W7d8dvPirpwDyPESOXzXtZkv6fev5YfjOfg8OoyeZ1vkDS9e6+S5LcfUeex4iRyeY1dkn7xb+ukfRiHseHUeLuD0vaOcAup0q61Xs8Kml/M3tffkY3MALt7E2V9HzS7e3xbWn3cfe3JXVJqs3L6DAasnmNk50n6d6cjgi5MOjrHP/Y8SB3/3k+B4ZRlc37+RBJh5jZI2b2qJkNNGOG4pPNa3y5pLPNbLukeyRdnJ+hIc+Gev3Om7GFHgBQiszsbElzJX280GPB6DKzkKRrJJ1T4KEg98aq56Pmo9Xz6dTDZjbL3f+3oKPCaPp7Sd939/80sw9L+oGZ1bt7rNADQ2VgRjt7L0g6KOn2gfFtafcxs7Hq+ZiqMy+jw2jI5jWWmX1SUoOkU9z9rTyNDaNnsNd5oqR6SQ+aWbuk+ZLuYkFkycnm/bxd0l3uvtfdn5P0J/UE3igN2bzG50m6U5Lc/beS9pF0QF5Gh3zK6vpdCATa2Xtc0nQzm2Zm49SzqOKulH3ukrQi/vUZkn7pFCovJYO+xmZ2hKTvqifIJp+zNA34Ort7l7sf4O5Rd4+qJxf/FHffVJjhYpiyOWf/RD2z2TKzA9STSvJsPgeJEcnmNe6QtECSzGyGegLtl/M6SuTDXZKWx6uPzJfU5e5/KfSgJFJHsubub5vZRZLulzRG0s3u3mZmV0ra5O53Sfqeej6W2qaepP2zCjdiDFWWr/HXJe0r6Ufxda4d7n5KwQaNIcvydUaJy/J1vl/ScWb2e0nvSPq/7s6nkCUiy9f4/0i60cw+r56FkecwAVZ6zOw29fxRfEA83/4ySVWS5O43qCf//kRJ2yR1Szq3MCPtj86QAAAAQA6QOgIAAADkAIE2AAAAkAME2gAAAEAOEGgDAAAAOUCgDQAAAOQAgTYAlDgza4/XgR7JMU4zs5mjNSYAAIE2AJSUeNfZXDhNEoE2AIwiAm0AGAYzu9LMViXdbjSzS9Lst9zMnjKzLWb2g/i2qJn9Mr69xczCg2z/vpndYGaPSbrazGrN7Bdm1mZmN0my+H7vMbOfx5+r1cyWphnPBWb2eHyfH5tZtZn9raRTJH3dzH5nZh9MecxPzWx5/OtPm1lTlj+jz5vZzfGvZ8XHVG1ml5vZWjP7tZkFZrbYzK42s61mdp+ZVcUf8+X4WFvNbE2869vY+Laj4/t81cwasxkPAOQbDWsAYBjMLCppg7t/yMxCkp6RNC+5s6CZ1UlqlvS37v6KmU12951mdrek9e6+1sw+pZ4W76cNsP37kg6QdKq7v2Nm35T0irtfaWZ/J+lnkqZI+rikhe5+Qfz5a9y9K2XctYkxmtlVkv7H3b8Vf46fufv6NN/rX0l6RD3d1r4naX78+/i/kpal+fE87O7/FP+5PCjpWkkNki5x90fM7HJJn5T0CfXMov9W0unufq+ZNUta6+4/Sfy84mP4gaQ73f3u+M91vaSL1dOt9W/cfc+gLxoA5Bkt2AFgGNy93cw6zewISX8l6ck07buPkfQjd38l/pid8e0flrQ4/vUPJF09yHbFj/NO/OujEvu5+8/NbFd8+1ZJ/2lmX1NP0PzrNEOvjwfY+0vaVz3tqwf7Xv/HzL4s6VeSFiW+D3f/unoC3UyPi5nZOZKekvRdd38k6e573X2vmW1VT/vs+5K+h2j860+Y2b9IqpY0WVKbpLvjbbZ/oJ4/MD5MkA2gWBFoA8Dw3STpHEl/LenmHD/XG4Pt4O5/MrMPSTpR0lVm1uLuV6bs9n1Jp7n7lngQfHSWzz9LUqek9yc2DDajHf96uqTXkx8X91Z8zDEz2+vvfrwakzTWzPaR9G1Jc939+fgs+D4p4/lfSe/NcvwAkHfkaAPA8DVLWijpSKWfGf6lpCVmVitJZjY5vv03ks6Kf71M0q8H2Z7qYUn/ED/mCZImxb9+v6Rud1+nnpnmD6V57HuJqpcAAAFdSURBVERJf4nnQScHya/F7+vHzOZJOkHSEZL+2cymST0z2u4+O82/f4o/rkbSN9UzA19rZmdk+H7SSQTVr5jZvpJ6H2tmi9Uzw32UpG+Z2f5DOC4A5A0z2gAwTO6+x8x+Jel/k9I6ku9viy/Ue8jM3pH0pHpmwC+WdEt8Rvhl9eQ+a4Dtqa6QdJuZtaknOO+Ib5+lngWNMUl7JX02zWP/XdJj8eM/pneD69sl3Whm/yTpDHf/sySZ2XhJN0o6191fNLP/I+lmMzsmaRY6k2slXR+faT9P0q/M7OFBHiNJcvf/NbMbJbVKeknS4/HxHCDpPyQtiM90XyfpvyStyOa4AJBPLIYEgGGKL/bbLGmJuz9T6PEAAIoLqSMAMAzW09xlm6QWgmwAQDrMaAMAAAA5wIw2AAAAkAME2gAAAEAOEGgDAAAAOUCgDQAAAOQAgTYAAACQAwTaAAAAQA78//3KhZ1trqXJAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7fe856a047d0>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "if make_graphs:\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",
    "    pylab.rcParams[ 'figure.figsize'] = 12, 6\n",
    "    pyplot.plot(coords[:,1], u, 'o', color = 'green', label='vx')\n",
    "    pyplot.plot(coords[:,1], v, 'o', color = 'red', label='vy')\n",
    "    big = np.linspace(0.0,h)\n",
    "    pyplot.plot(big, exact_vx(big), 'D', color = 'purple', label='exact_vx', markersize=2)\n",
    "    pyplot.legend()\n",
    "    pyplot.xlabel('y coords at x=xmax')\n",
    "    pyplot.ylabel('velocity')\n",
    "    pyplot.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Relative error against analytic solution: 1.40525963472e-05\n"
     ]
    }
   ],
   "source": [
    "ana_u = exact_vx(coords[:,1])\n",
    "\n",
    "abserr = uw.utils._nps_2norm(ana_u - u)\n",
    "mag = uw.utils._nps_2norm(ana_u)\n",
    "relerr = abserr / mag\n",
    "\n",
    "from mpi4py import MPI\n",
    "comm = MPI.COMM_WORLD\n",
    "\n",
    "\n",
    "if uw.rank() == 0:\n",
    "    threshold = 1.0e-4\n",
    "    print \"Relative error against analytic solution:\", 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": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Further analytics\n",
    "# pyplot.plot(coords[:,1], du_dy, 'o', color = 'purple', label='du_dy')\n",
    "# pyplot.plot(coords[:,1], du_dx, '+', color = 'black', label='du_dx')\n",
    "# pyplot.plot(coords[:,1], dv_dy, 'x', color = 'orange', label='dv_dy')\n",
    "# pyplot.plot(coords[:,1], dv_dx, '.', color = 'red', label='dv_dx')\n",
    "# pyplot.legend()\n",
    "# pyplot.xlabel('y coords at x=xmax')\n",
    "# pyplot.ylabel('velocity gradients')\n",
    "# pyplot.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Further analytics\n",
    "# pyplot.plot(coords[:,1], sr, '-', label='exact_shearSR')\n",
    "# pyplot.plot(coords[:,1], strainRate.evaluate(ids)[:,2], 'o', color = 'purple', label='sr_shear')\n",
    "# pyplot.plot(coords[:,1], devstress.evaluate(ids)[:,2], '+', label='tau_shear')\n",
    "# pyplot.legend()\n",
    "# pyplot.xlabel('y coords at x=xmax')\n",
    "# pyplot.ylabel('strain rate')\n",
    "# pyplot.show()"
   ]
  }
 ],
 "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
}