{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Joule heating model\n",
    "The thoretical information comes from [here](https://reference.wolfram.com/language/PDEModels/tutorial/Multiphysics/ModelCollection/JouleHeating.html). \n",
    "\n",
    "## Electroestatic model\n",
    "In this model the electric potential field $V(x,y,z)$ is assumed to be independent of the temperature $T$. That is, the electrical conductivity $G$ is kept at a constant value at all time. By Guass's law the stationary potential field satisfies Poisson's equation:\n",
    "\n",
    "$$\\nabla \\cdot (-G \\cdot \\nabla V)=0$$\n",
    "\n",
    "## Heat transfer model\n",
    "The heat equation is used to solve for the temperature field in a heat transfer model:\n",
    "\n",
    "$$\\rho C_p \\frac{dT}{dt}+\\rho C_p v\\cdot \\nabla T + \\nabla \\cdot(-k \\nabla T) = Q$$\n",
    " \n",
    "For a steady-state heat transfer model the transient term in the equation are set to zero. Since a solid is modeled, the internal velocity also vanishes and the heat equation simplifies to:\n",
    "\n",
    "$$ \\nabla \\cdot(-k \\nabla T) = Q $$\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Solving the PDE Model\n",
    "In order to solve the PDE model, first the electroestatic model will be solved first to simulate/obtain the potential field $V$ of the wire. The heat transfer model is then constructed to show the heating effect of the electric current.\n",
    "\n",
    "## The electroestatic model\n",
    "There are two types of the boundary conditions involved in the electrostatics model. At both ends of the wire an electric potential difference $V_0=0.2$ is applied.\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "from fenics import *"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x7f89569e9160>,\n",
       " <matplotlib.lines.Line2D at 0x7f89569e92e8>]"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYQAAABJCAYAAAA9vFENAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAHCpJREFUeJztnXtsXNd95z+Hj6E4FIfkDMnh25T1sh62Y0u2bLcpHHsLO01qB029dopFtoESx9lku2kehpwCeaFBWzQu0sJFF4Gd3U22eW0cLAw3hYu6zUOWZEl0HMuSLVu2LL7fFN/DGY1O/7j3jO8dDjmXc88lZ6TzAQjM48w9v3N57++c8zvf87tCSonBYDAYDGUbbYDBYDAYigPTIRgMBoMBMB2CwWAwGGxMh2AwGAwGwHQIBoPBYLAxHYLBYDAYANMhGAwGg8HGdAgGg8FgAEyHYDAYDAabio02YC00NjbK7u7ujTbDYDAYSoaenp5xKWWTl7Il1SF0d3dz8uTJNf/ua1/7mut9NBqloiL4po+Ojrreh8NhNm/evO71AjQ2NlJWFvyEMLvu2tpaqqur171eIQSNjY0IIda97vr6ekKh0LrXW1lZSX19/Ya0eaPuqU2bNhGJRAKvN5FIMDMz4/osFotRXl4eeN2qzV/5ylcK+r0Q4oLXsiXVIRTCpUuXln02OTmZeV1eXk5nZ2cgTiv74l1YWGBhYSHzPhaL0dTUpP0Gnpubc9UDMD4+nnkdCoXo6OigqqpKa72wvM2zs7PMzs5m3sfjcRoaGrS3ObteKSVjY2OZ9+FwmLa2NiorK7XWm6vuixcvut63t7cH4rSy602lUq42RyIRWlpatDstKeWyup33lBCCzs5OampqtNYLy9ucSCRIJBKZ99FolKamJu2Dn+np6WUdwsTEROZ1ZWUlHR0dbNq0SWu9kHuAFxRXfIdQUVFBZWUl+/fv54477qCvr4/e3l76+voYGhoinU7zzjvv0NjYSGdnJ11dXXR2dhKNRn07reeff54jR47wpS99iaGhoUy9fX19zM/PMzExwfz8PJ2dnZm/9vZ2305rfn6eb37zm7z//e9nx44drjaPjo6STCY5f/488Xg8U29XVxd1dXW+6gV4+umnGRwc5JFHHmFgYMDV5qWlJUZGRpibm8uc587OTlpbW307rYGBAZ588kk+8pGPEIvFXPWOj4+zsLDA22+/TVtbm6vNOpzWU089RSgU4oEHHqC/vz9T98DAAKlUioGBAebm5lz1Njc3+3Zar776Kk8//TSf+tSnkFK62nzx4kVmZmZYWFigvb09U29HR4eWwc/jjz/Ojh07uOuuuzJ19vX1MTg4SDqdpre3l2g0mvk/d3V1EYvFfN9Thw8f5vnnn+fQoUOMj49n6u3t7WVubo7JyUnm5+fp6OjI1Nve3u57xnbp0iW+8Y1vcPfdd3PDDTe47qnh4WFSqRTnz5+nubnZ5Ud0zNieffZZXn/9dV/H8MoV3yE42bx5M7t27WLXrl2ANaIaHBzM/GNfe+01fv3rXwNQU1PjctStra0FT4nLy8vp6Oigo6MDsEZYU1NT9Pb2Zup+8803ASgrK6O1tdXlPAoNMwkhqK+vp76+nuuvvx6ApaUll9N6+eWXOXHiBGCNKJ31xuPxgp1WZWUl3d3dqDUfNWJ3Oq3XXnsNsDpt5bTUX6FOSwhBLBYjFotx0003AdbMzOm0jh8/ztGjRwFrROlss58w06ZNm9i2bRvbtm0DIJ1OMzIykmnzhQsXePXVV4F3Z2lOp1XojE0IQXNzM/F4nFtuuQWwZmbOc/3CCy9w+PBhAJqamjL1dnV1+XJaNTU1XHfddVx33XWA5TgHBwcz9Z49e5aXX34ZgOrqate5bmtrK/ieKisro729nfb2dm677TaklFy8eNHlqH/+859nzk9LS4vLUfuZsUUiEfbs2cOePXsASCaT9Pf3Z9p86tQpenp6AMvnOAc/QczYdHJVdQjZVFZWcs0113DNNdcAltMaHx/PXFC9vb2Znrm8vHyZ0wqHwwXVK4QgGo0SjUZ5z3veA1hOy+moT548ybFjxwBoaGhw3Uh+wkxVVVVs3bqVrVu3AnD58mWX0+rr6+P06dOA22l1dnb6CjMpp9Xc3Mz+/fsBy2k5HfWRI0e4fPky8K7TUm32E2YKh8Ps3LmTnTt3ApbTcs7Y3nzzTX7zm98AllPPdlqFztjKy8tpa2ujra0t47Smp6dd5/oXv/hF5vyoGZty1H6cVm1t7TKnNTAw4Pofv/TSS4DltJxt9uO0KioqMvaDdU9NTEy4HPUbb7yROT9q8KOcZqEzNiEEDQ0NNDQ0cMMNNwBWOMl5fb300kscP34cgLq6Opej9jNjC4VCXHvttVx77bWAdU+Njo66Zi9nzpwBLJ+TPWMLIsxUKFd1h5CNEIKmpiaamprYt28fYMXjnRfz0aNHeeGFFwC0hpnC4TA7duxgx44dgDW6dDqtt956i1deeQVwOy2/YSY1I2ltbeXAgQMAGael6v7lL3+JlNLltHSEmWpra9m9eze7d+8GyIRXVL1nzpzJOC3njK2rq8tXmKmioiJzLLCc1uTkpMtR55qxqf+zjhmb02k5R5fOGVtdXZ22MFMoFGLLli1s2bIFsJzW2NiYy2nlmrH5DTOpxf3GxsbMjG1+fn7VGZvTUfudsW3fvp3t27cD1j01PDycqff8+fOcOnUKsAZKusJMZWVltLS00NLSkpmxzczMuPzI4cOHUc+iCSLMVCimQ8jDWsJM4XDYdTGn0+mC6/UTZorFYr7aXFdXx/XXX19QmEmN8AvBS5hJzdiyw0x+4vFewkwnTpzIzNiyw0x+HjKVHWa6fPkyw8PDnsJMyWSy4HrLysqIx+PE43HXjM1rmMlPm/2EmXKJRLyiZvmFhJlaW1sLrhdyh5mcg5/Vwkx+/MhaEaX0xLT9+/dLHbJTP7HxbFSYycuFWldXV3CYKReLi4vLFC25KC8v16pmUmEmL0SjUa1qptnZWebm5vKWq6qq0qpmunTpkkvBsxqNjY3a1EwqzLS4uJi3bE1NjVY1UzKZdClpVqO5uVlbbFyFmVKpVN6ykUhEq5opkUgwNTWVt5wKgeryIyrMtJo//vKXv1zQ9SyE6JFS7vdU9krvEJQ6YCXUSFyn05qbm2NwcDBvuVgspkXNpFBhpmzJaTahUEiLmsnJ9PS0p04iHo9rUTMpVJgp32g5HA5rUTMpVJjJi8Nsb2/X7rT6+/vzzsYikYhWp6XkptPT06uWE0JoUzMp1BpbPqLRqBY1k0KFmebn51ctp9YGdO4/mZmZYXh4OPPexz4Ezx2Cp5CREOJe4G+BcuBJKeVfZn1fBXwX2AdMAA9KKd+xv3sMOAikgT+RUj5nf/4d4IPAqJRyrxc7CsEpOz1w4IBrejgyMpKRyGXH8erq6rTJTg8dOuQKM/X19ZFIJJiYmCCRSGhTMymcstOtW7e62jw+Pk4ymeTChQva1ExOlOz04YcfdoWZ+vv7SaVSjIyMsLi46JoS65ixOWWn9fX1rnM9NTXFwsICFy5c0KZmcqJkpx/+8Idd51pJMAcGBrSqmRRKdvrII49krmPV5tnZWWZmZkgkEq4wU0dHhxanpWSnd955p6veoaEhpJT09fVpVTMplOz00UcfXbYGsri4yOTkJIuLi9rUTAo1sLzrrrvYu3fvMil3KpXiwoULWtVMiqKSnQohyoG/B34X6AdOCCGekVKecRQ7CExJKbcJIR4C/gp4UAixG3gI2AO0Af8qhNghpUwD/xt4AqsjWRfq6uqoq6tj716r/1laWnLF8V555ZXMTuja2tplcjE/EsyNUjOp2HixqZl6e3tzxsaNmmntqPUAo2a6stVM64GXbvNW4JyU8m0AIcQPgfsBZ4dwP/BV+/VPgCeEdVfdD/xQSrkEnBdCnLOPd1RK+UshRLeORhRKVVVVTrmY02kpCabaieh01H6cViFqJnVBX2lqJnUjGTXTladmUmEmo2bSr2YKAi8dQjvQ53jfDxxYqYyU8pIQYhqI2Z8fy/pte8HWBoxTLnbrrbcCltNyOupf/epXmYUf50jLb5jJj5qpFDfNgT81U9Cb5pxqJrW7WZ3zIDfNraZm0rlpTqmZVL1BbporVM3U2dnpa8a2UZvm/KiZdIWZCqXoZadCiIeBh4HM9HA98RNm8ivBLDTMFI1GC64316a5xcVF18W8WpjJj0jBz6Y5P7FpL2Em54wtO8zkBz9hJj8SzLKyssymuQMHDqw5zOQHP2EmP/eUnzCTH9npSmEm5+BntTCTnzav2dZ8N7AQ4nbgq1LKe+z3jwFIKf/CUeY5u8xRIUQFMAw0AYecZZ3l7PfdwLNeF5V1yE4rKytpbm7WKsEcHR31dHM2NTVpVzN5kZ1WV1cHombyQjwe165mcibLW4na2lrtaiavUtsg1ExeZKcNDQ3a1UzOpIgrUVFREYiayYvsNBaLaVUzzc/Pe5KdVlVVac2mq8JMRS87tR38G8DdwABwAvgjKeVpR5lPA9dLKR+xF5X/QEr5n4UQe4DvY60btAHPA9vtReV16RC8yE6DkGB6uZFisZgWNZNChZnydU6hUEiLmkmhwkxebqR4PK41BbgKM+UjHA5r338yNjbmqXNqa2vTLsH00iFHIhGtTkuFmfJ1Tkp2qjM2Pjs76ynrZzQa1brTV4WZ8nVOlZWVWtRMChVmcmaRLQrZqb0m8BngOSzZ6XeklKeFEF8HTkopnwGeAr5nLxpPYimLsMv9GGsB+hLwaUdn8APgTqBRCNEPfEVK+dQa25q/gQ7Z6f79+13Tw7GxMdLpNP39/cvieLW1tb7rVrLTRx991BVm6u/vz2z6SSaT2tRMCqfstLu727WgNzU1lUnGpUvN5ETJTg8ePOiaEquOamRkZFmbdWyaU7LThx56iEgk4gp9qKyfAwMD2tRMTpTs9EMf+pDr+hoeHuby5csMDg5qVTMplOz0k5/8JKlUytXmhYUFZmZmSKVSrvUAXSnAlez0ve99r6vekZERpJT09/drVTMplOz0i1/8IiMjI66w3tLSEpOTkySTSW1qJoVTdrpr1y7X/1ltouvv79emZnJSVLJTACnlz4CfZX32ZcfrBPDACr/9BrBsiC6l/MiaLNWAio3feOONgBUbdzqtnp4eXnzxRcB6yEm2XMyPBLNQNZOf5FfO2HixqZnOnTuXMzZu1Exrp7y8PDOggfVTMwEbpmaqrKy86nIzrQdFv6gcJNXV1TnlYrkkmFVVVcuclo7kV2tRMzmd1pWkZlrNaRk109q5GtVMfnIzlaqaKQiKx5IiwCkXu/322zNOy+moz507B7zr1HWFma7GTXN+1EzrmQI8e3QZ5Ka51dRMQWyac3bGQW2aK1TNpKTchXK1bZrTgekQVsHptAoJM/mVYBYaZmpoaPDVZj+b5vxQXV1d8KY5P3gJM6kZW64wkx/8hJn8ZMH0u2muUNa6aS4SiWTq9Ctr9rNprqWlxVebs8NMasamzvdqYab1zDd3xSe3A7fstKamJpDkV15kp62trVqVF9PT055kp2r2oVvNlA/lSHSrmbwoe2KxmHY1kzPR2EqEQqFA1ExeZKfNzc3a1UxeMryGw+FA1ExeZKctLS3as+k6lT0rUVtbqz2b7uDgYPHLTouJoGSnQTit7Ady5yIWi2lRMymWlpY8yRFDoZAWNZNChZnyZYQEy2npTgHuZX9AOBwOJAX40tJS3rKtra3anZaXLKuRSES701LPIV8NIYQ2NRO8mwLcy+AnGo1qTwHuZfBTWVmp9fGYKszkTPVeFLLTUkfJTvft28fNN9/smh5OTk5mwhLZC3o6nJaSnX7hC19Y9vD1S5cuMTExQTqd1qZmUijZ6b333ktXV5erzTMzMySTSYaGhrSpmZwo2enHPvYx15R4aGgoE/bS+aQ5hZKdPvjgg2zevNkV+pifn2dhYYGRkRFtaiYnSnZ63333uVJ+qPz2w8PDWtVMCiU7/cQnPkEymXTJmhOJBDMzM6TTaW1qJidKdnrHHXe4zvX4+DhSSoaGhrSqmRRKdvr5z3/etVitsumqe1qXmkmhBpbve9/72Llzp+ueunjxIqlUiqGhIW1qJidFJzu9EnDGxm+++WbgXbmYuqCPHTvGkSNHAGv07ryY/YSZNlLNpGLjRs10ZauZ2tra1l3NBGyYmqmqqqrgJ82VqpppPbhqOoRcZMvFlNPKJRcLh8MuR60r+dVa1EzqwjJqJu8YNdPVp2ZaawrwUlUzBcFV3SFks5LTck4Pz549C1hOS1eYya+aqVQ3zW1UCvC1qpmyR5dBbppbTc0UxKY5Ve/VngI8W81UiinAdWA6hFXwG2byg58wU319fcH1+t005wc/YSY/GSG9hJmcM7bsMJMf/ISZ/GbT9bNprlDWummuoaFBm+zUz6a5eDzuq825wkzOe2q1MJORna6ADtlpQ0NDIMmv8slO1U2lW83kRXnR1NSkXc00MDCQt9ymTZsCUTN5kZ22tLRoVzN5UXBt3rw5EDWTF9mpn9h4LmZnZz3JTuvr6wNRM+WTnapZum41kxfZqUpMqQsVZjKy0zUQlOw0Ho9rlYuNj4+TSCTylo3FYtqdlpcsq6FQKJAU4F50483NzdpTgHvJshoOhwNJAe7l/mlpadGeTdeLrDkSiWhPAe5lL4YQIpAU4AsLC3nLNjQ0aN1/kkgkPHWKKq2+7hTgyWQy85mRnWrAKTu98cYbXXG86elp0uk0Y2Njy+J4OiSYSnb6uc99zjU9HBwc5PLlyxktuS41k0LJTu+55x46OjpcbZ6fnyeZTDI2NqZNzeREyU4/+tGPutqssmCOjo5qfdKcQslOH3jgAWpqalyhj0QikXEmutRMTpTs9AMf+IBLdqo659HRUa1qJoWSnR48eJClpSWX7DSVSmU6DV1qJiePP/4427dv57bbbnOd66mpqcz/WaeaSaFkp5/97Gdd4Z6BgQHS6TRTU1MIIbSpmRRO2em2bdtcbZ6dnSWVSjE2NqZNzeTEyE4DQD2qrqWlhVtuuQWAmZkZ1z/28OHDmRFfc3Oz62L2E2baSDWTio0bNVOwaiZgmZppYWHBtR4QpJqpo6PDqJk8PmmuVNVM68FV0yHkIhKJsHfv3ozTUs8JUP/YU6dO0dPTA1ix4WynVeiUuFA1k7qwdKuZnDfweqqZnPnsr0Q1k5cU4OutZsoWBqx3CnCnBNOomfRsmtOJpw5BCHEv8LdYD8h5Ukr5l1nfVwHfBfYBE8CDUsp37O8eAw4CaeBPpJTPeTnmRhAKhXJKMJ1O68yZM4DltHSFmXKpmZTTUn9BbprL5bRUvWbTXHiZ09K5aW5ycnLFGduVsmmu2FKAq7pXUzPpTAGeTqddM7agNs3pIO+VLYQoB/4e+F2gHzghhHhGSnnGUewgMCWl3GY/QvOvgAeFELuxnp62B+sRmv8qhNhh/ybfMTccpwSzkDCTH1ZyWl7CTH5GeH7DTH7IFWZyjrRWCzP5lWBmz9iyH7binLFlh5kKxem0VgozZTst57n2m013tU1zzhlbrjCTnzavFGZaacbmdNR+8BNmam5uLrheNcsvNMy0nnh5pvLtwFellPfY7x8DkFL+haPMc3aZo/YzmIeBJuCQs6wqZ/9s1WPmwq/sVGn6daLCTKsl/Oru7taqMoF3w0z5ZKdBtFltmlvvepXTcib8Wq+65+bm8spOg6hXOa18qjXddSunlU9hE0SbVZhpNdVaY2OjrzTvuVAztnyy0yDarMJMGy079TL3bQf6HO/7gQMrlbGfwTwNxOzPj2X9tt1+ne+YAAghHgYeBnyPELxk4ywE56YVFWZyjliXlpZc8jFdhMNh13rC3NzcMjliUG1ua2vLvFZT4vWoNxKJuLJZTk9PL6sriLpVBk+FUpUEXS9Yi9UK5bSys6wGUbfKkaRYWlpalmU1qDY3NTVlXqsZm3OvTzqdDqTuTZs2udq8sLCwbNAVVJtbW1szr3PdU+uRA6noF5WllN8Gvg3WDKGQYxSq3zUYDIarCS8dwgDgDGR12J/lKtNvh4zqsBaXV/ttvmMuo6enZ1wIccGDzbloBPLv2ipOStX2UrUbjO0bhbFdP9d4LeilQzgBbBdCbMFy2g8Bf5RV5hngvwJHgT8E/k1KKYUQzwDfF0L8Ddai8nbgOCA8HHMZUsqmfGVWQghx0mscrdgoVdtL1W4wtm8UxvaNJW+HYK8JfAZ4Dksi+h0p5WkhxNeBk1LKZ4CngO8JIc4Bk1gOHrvcj4EzwCXg01LKNECuY+pvnsFgMBi8UlK5jPxQyr13qdpeqnaDsX2jMLZvLMWzRS54vr3RBvigVG0vVbvB2L5RGNs3kKtmhmAwGAyG1bmaZggGg8FgWIWS7BCEEPcKIc4KIc4JIQ7l+L5KCPEj+/sXhRDdju8esz8/K4S4x+sxi9V2IUSnEOLfhRBnhBCnhRD/o1Rsd3xXLoT4tRDi2VKyXQhRL4T4iRDidSHEa/au/lKx/U/t6+VVIcQPhBD+871rslsIEbOv6TkhxBNZv9knhDhl/+bvhAhmt5Zu24UQYSHEP9nXymkhxIbnbsuJlLKk/rBUSW8B1wIh4DfA7qwy/w34n/brh4Af2a932+WrgC32ccq9HLOIbW8FbrbL1AJvlIrtjt99Dvg+8GypXDP2d/8H+Lj9OgTUl4LtWNkCzgPVdrkfA39cRHbXAL8NPAI8kfWb48BtWNL1fwbeX2TnPKftQBh4n+Na+VUQtvv9K8UZwq3AOSnl21LKJPBD4P6sMvdj3awAPwHutkcS9wM/lFIuSSnPA+fs43k5ZlHaLqUcklK+BCClnAVe4930IEVtO4AQogP4APBkADYHZrsQog74HSzJNVLKpJQy//NMi8B2u1wFUC2sjaRhYLBY7JZSzkspDwOuBE5CiFYgIqU8Ji3P+l3gQ5rtDsR2KeWClPLf7ddJ4CWsDblFRSl2CLlyK2U7QFduJcCZWynXb70cUwdB2J7BnrbeBLyo0eZldq1UP4XZ/i3gUaDwdKX5CcL2LcAY8L/scNeTQoiaUrBdSjkAfBPoBYaAaSnlvxSR3asd05lVsRjv07wIIeqB3wee922pZkqxQzDkQAixGXga+KyUMv8Dd4sAIcQHgVEpZc9G21IAFcDNwD9IKW8C5rGz+xY7QogGrBHuFqwMAjVCiP+ysVZdHdgzsh8AfyelfHuj7cmmFDuEteRWUv+AfLmVvBxTB0HYjhCiEqsz+Ecp5U8DsNtlV3b9ucp4tP23gPuEEO9gTcvvEkL83xKxvR/ol1Kq2dhPsDoI3QRh+38Czkspx6SUKeCnwB1FZPdqx3SGWYrxPs3Ht4E3pZTf0mCnfjZ6EWOtf1gjs7exRjdqwWdPVplP417w+bH9eg/uRba3sRaQ8h6ziG0XWLHUb5Xaec/67Z0Et6gciO1YC4M77ddfBf66FGzHSjV/GmvtQGDFwv97sdjt+P6Pyb+o/HvFdM7z2P7nWAO3siCucy1t32gDCvyH/R6WmuYt4M/sz74O3Ge/3gT8P6xFtOPAtY7f/pn9u7M4VvlzHbMUbMdSNEjgFeBl+0/7TRLUeXd8fycBdQgBXjPvAU7a5/7/Aw0lZPvXgNeBV4HvAVVFZvc7WHnR5rBmY7vtz/fbNr8FPIG9ubbYbceaZUgs0Ye6Tz8e1PVe6J/ZqWwwGAwGoDTXEAwGg8EQAKZDMBgMBgNgOgSDwWAw2JgOwWAwGAyA6RAMBoPBYGM6BIPBYDAApkMwGAwGg43pEAwGg8EAwH8AaDlVDbUtrc4AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "#Creating mesh and define function space\n",
    "Nx = 0.013\n",
    "Ny = 0.0015\n",
    "eleX = 10\n",
    "eleY = 5\n",
    "mesh = RectangleMesh(Point(0,0),Point(Nx,Ny),eleX, eleY,'left')\n",
    "V = FunctionSpace(mesh, 'P', 1)\n",
    "plot(mesh)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Boundary conditions\n",
    "tol = 1E-14 # tolerance for coordinate comparisons\n",
    "#at y=1\n",
    "def Dirichlet_boundary1(x, on_boundary):\n",
    "    return on_boundary and abs(x[1] - Ny) < tol\n",
    "#at y=0\n",
    "def Dirichlet_boundary0(x, on_boundary):\n",
    "    return on_boundary and abs(x[1] - 0) < tol\n",
    "#at x=0\n",
    "def Dirichlet_boundarx0(x, on_boundary):\n",
    "    return on_boundary and abs(x[0] - 0) < tol\n",
    "#at x=20\n",
    "def Dirichlet_boundarx1(x, on_boundary):\n",
    "    return on_boundary and abs(x[0] - Nx) < tol\n",
    "\n",
    "#bc0 = DirichletBC(V, Constant(0), Dirichlet_boundary0)\n",
    "#bc1 = DirichletBC(V, Constant(0), Dirichlet_boundary1) \n",
    "bc2 = DirichletBC(V, Constant(0), Dirichlet_boundarx0)\n",
    "bc3 = DirichletBC(V, Constant(0.003), Dirichlet_boundarx1)\n",
    "bcs = [bc2,bc3]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x7f8956953eb8>,\n",
       " <matplotlib.lines.Line2D at 0x7f8955f23278>]"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYQAAABWCAYAAADPPOFDAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJzsvXmQI9l93/nJO3HfBaAKVQCqu7q7+j5nOAcpahg8JYqiZZlcybosWbLWWm/YG7thaSN2ZXnDthwbYVsryRYlSpYoLQ8dvCSKx5Aczgw5w+lz+p6+gDq67gMoAIUbuX+8bFQfVQVMN4bkzOIbgYgq4GXi5UNm/vL93vf7/UmWZdFHH3300Ucf8ve7A3300UcfffxgoB8Q+uijjz76APoBoY8++uijDxv9gNBHH3300QfQDwh99NFHH33Y6AeEPvroo48+gH5A6KOPPvrow0Y/IPTRxxaQJCkrSVJZkqTiXa8nJUmy7nuvKEnSR+xt/rv9+Yfu29d/st//+e/LwfTRRxdQv98d6KOPH3B80LKsZ+/8I0lSyv7Tb1lWY4ttrgE/C3ze3kYF/hFw843rZh99PDr6M4Q++ug9vgg8LUlSwP7/fcB5YO7716U++uiMfkDoo4/eo4KYHXzU/v9ngT/7/nWnjz66Qz8g9NHH9vicJEk5+/W5u95fuuv9nCRJ4/dt92fAz0qS5Ad+CPgcffTxA47+GkIffWyPH99iDSG8zRoClmW9KElSBPjfgb+1LKssSdIb2tE++nhU9ANCH328cfhz4P8Afvj73ZE++ugG/YDQRx9vHH4HeAF4/vvdkT766Ab9NYQ++ng45O7TIfyr+xtYlrViWdbXrX7RkT7eJJD652offfTRRx/QnyH00UcfffRhox8Q+uijjz76APoBoY8++uijDxv9gNBHH3300QfQDwh99NFHH33YeFPpECRJ2pISJfHGq0AttmdkPdiHHvXprt1YVqtj460UsT3lk3Xqh/Q9eNbo2AcJkHr2M2yJVod+yN+DsejYB3ssHhZbX3obTQCr2eEaUb4H12mnPsgS3YvGH/6qaTW3/1xRHnrX2+LuQ2uIPixZlhXpZts3VUAwcLBDOUCutUTOWmSdAgASMl4piF8KEzCH8KtRDNmx6T4k0+j8RY7Ntz09/1lKtRX2ht/FamOe1fI0+fXbNK266J/qxu9M4HcmCDiGcQWGkDvcGFtG55+gaW60yecnOXP6v7Jj5wdQVZN8boJ8Pku5vCyOT1bw+BJ4A0l8gRTeQBJNd9n76Xxjahqdr5SGIfHaV/4b9XKB4eMfpLiQobiQpbQ0SatRA0B3+XDF0vYrhSMQR7rvxtjUO34VrW3aFCevk/3L/0r0Pf8ASZYpT2dZn85Qz9ljoSiY8RHMZArHSBrHcArF6dr8e7Qu+mJsfnO4/Tu/A60Wgfe+l/JkhsqtDLWJSay6cLZQA36M0RTGaApzNI2ajDwwFg9A7xT4QdY37jjrZ64y/39/gvAvfxir3qRybYLqtUkai6sASKqKOTaEY88Ijj0jmLsTKC5xnpt6veN3ObQtXTracOlVXvmlT6C4DJIfOUbu4gz5SzMUrs7Tqou+OmMeggfihA7ECB6M40kGkOR7zzm3Wt32e/xaedvPr391gq/8+rf58X9/jEqhztTZFabOLbM2K7ZTDZmhAwHGjnlIH/GTPOTD4dn8BPAp6x2PO6QWN33/V5+5wsCwzgf+cZjsmRwXTld47UKFurhEGBzROHDcZP8xkwPHTY6MNbZ8mLsDv9w5QAXkjXvcJz5d4Ff+5fJEx41svKkCgozCkLyDIXkHADWrSs5aJG+JADHZusbE+lUAnLIXvxoloETxq1Fcsq/jYHcDVdYJO1KEHLsBaFktipUFVstT5Nanya1PM7d2BQBlQsfnSuB3DxNwD+NzJVCVLgJSF3C5BgiFdjM4eAKAslRibXWCtdUs+dUJbmdeZPrWtwBwuCJ4A0k8A2k8oRQOd6Q3Y6Gb+BJ78CX2ANBqNSmvzFBcyFBYylKcucHqjTMAyLqJayCJK5bCHU3jjCZB781YGJEYruROAkefBKBeXKM8nWkHiNXvPMfqi98AQI9EcYykMe0AoQXDPRkL2eHAOT6OeUiMhdVoUJu+TeVWluqtLOWr1ymdFGMhOUyM0STGzhTmjiR6egTZ6CI6dgF9JIYxmsD77sfFd5VWKF+donx1kvLVSVY++6KYTUgS+nAEx54RvPsGce4dRo/6ezIWmtck/OQOwk+K67RVa1C4vkD+4gyly9PMvzzJ1JdfE209BsEDMREgDsQJjA/07K4U3eMjssPL8Y+OApCfW2faDg7TZ1d49g8yWGIoiO1ykz7iF6+jfgKDZk/GwhdUefzdPt7/XjEdqFZbXLtQ5cKpChdOV/j2syW+9JfiwdYfkDl8TOPIcZ0jx3X2H9QxzO+t/9WbKiDcD10yGJASDJAAoGU1KZjr5Brz5BrzLNanmKldB0CTDPzKAIHGEH49hk8bQJEe/fBlScbriOF1xEgG7ZtzPU9ufZrV2gy54hS3Zp9HTD0lPI4ofvew/RrBMEKP3AcA3fAQju0nHNsPQLNZp5ifJr+aZW11guX5y8xPnwJA1V14Qyk8oRTeUAq3P4Gs9GAsZAVXeBhXeJiw8UNYlkWtuEpp9hal+SzFuQxzp74KWCDJOMKDOOMpXPE0zsE0utv/yH0A0NxetD2H8O45BEBdqlGZmaQ8maUymaFw6VXyp18GQHG5cYykMVJpHMkUxmACSX30sZBUFSOVxEgl4RkxFo2lZaq3slQmMlRvZsl//ivkAWQZfXgQY2dKvHakUAfcj9wHADXoxfPkPjxP7gOgVa5SuXHbDhBTFF68QP6r9nkRcOMaT+AcH8a1N4FjNIakPnpeQ9ZVfPsG8e0bxKsfwLIsStN5li/MsXJhlpULc8x/RzzESqpMaHeIyMFo++UIbT5jf73wxZz43u9k3/vF/cJZzTP5ap7M2RyZszlOf3GW73xqWrSNGqSP+Bk/5mTnUS/De1wo6qPfnA1D5sBxBweOi2OyLIvJW3UunKpw/UyRs6dqPPesmCFpOuw7IALE0RM6h4/p+CNvbPrxTR0Q7ocsKQTUKAE1CojBXm/lWW0skGvaQaI4BdhpJi1CQIsT0GP4tRiG4uxJPxyaD4fPR8xxGIBGs0q+NM1qcZJccYqZ5XNMLZ4EwNC9+D0j+LxJ/J4R3K4osvToF6GiaPiCaXzBNCDWHor1ZdaWsxSWs6ytZFmZvQSAJKu4Awm8oRSuaBpPJIVmbJ5aeT2QJAnDE8TwBAnuOg5Ao1pmfT5LcT5LcSHDyuVXWD7/IgCaJ4ArnsIZH8UVT2GG4vSC9yDrOs7UTpypnQBYrRa1xXnKkxnKUyJIFK9cEH1WVYyEnWZKpjFHtk4zvR5IkoQWCaNFwrjecRSAZmmd2q0JKjcmqN7MUnz+ZQpfF2OhRoIYO5MYY0mMsRTa0EDnNFMXkB0GzgOjOA+Ip2ar2YK5GdavTFO6MkXp8hT574hZtqSrOHcN4hofJnBwEPf4EKrbfOQ+SJKEe9iPe9hP8gNiRlXLV9oBIndxhtf++gpXPnkRAHfCw8BdAcI3ZjyQZnoYmC6VXU+G2PWkeChrNS1mrxXbASJ7Jse5L88DYDhl0gc9jB31sPOol9HDHpyeHjw4SBLJHTrJHTqhnxJpq5XlJudO1zhzssbZUzX+/E9K/MkflAAYHVU4cULnxHGdEyd0duxQejKTuYO3VEC4H5Ik4VL8uBQ/CXYBUNdb5GpzrNbnWK3NMrF+nuz6OQCcio+AHsffGCZgDOHSgr1JrSgGIe8OQl57Cm21KJbnyRUnWV2fJleYZH5ZnPyKrOPzJPB5RvB7k/g8w0g9+JkkScbpjeL0RomlRTqhVilSWMmytpyhsJxl5voLWNeeA8D0RPBE0njCKTyRFKanR2kmw4F3ZBzvyDhNHaxmk/LSDOuzGUqzGYq3b5K7dhYQaSZnPIlzKIVzMI0jnkTpQZpJkmWMaBwjGsd/QqSZapU1KhMZyhPiCT73wnPkviXSTFokKmYPO1KYo6OooVBPxkJxOXEcGMdxQJRSsBoNapMzVG9mqd7KUL50ndJLYiwkh4mxcwRjLIk5lkIfHUbWH/3BQVJkzHQURzpK6APHAKgvFyhdmRJB4vIUC3/9HRb+0gIJHMkInr1DuPcm8OxLYMR6k4rVfSbxp1PEn07hVqs0a01WXlti8fwCi+fnuf3SNLe+dAMAw6sTOxgmfihC/HCE6N4gqtmDGa4iMTTuYWjcw9M/NQxAc2GV62fWuHGmwI0za/ztf5tup5mGdonZw5ETOnuOuRgY0noyFsGQwjPvcfDMe8QsolqxuHRBBIiLp2t89asVPv1psR4SCEgcP74RIA4e1OARnmvf0gFhM+iygwEzzYApnpxbVpN8fYFcfY7V2hwLlSy3y+IJSZNN/MYgAXMQvzlEqxNtoEvIkozXGcfrjJMwngKgUs2RW5skV5gkV5ggM/0t7qSZ3O4YPl8Sny+JJPeOmqCbbkKD+wkNbqSZ1orTFJayFBYzrE5fZPHWKwCohksEh3AK51Aaq9l5kbEbSIqCMzqMMzpM+PA7sCyLemGV0myG9dkMxdkMC9/ZSDOZkUGcQylcg+me/R4AqseLe/8h3PtFmqlVq1GdttNMExmKF8+zduq7AChuN0Y6jZlKYabTWM3e9ENSVYzREYzREdCfFmmmxRWq17NUr09QvT5B/m++JtJMioyejGPuSmLuGqFZ2n6h9fVAC3nwP70X/9N7AWhVajRuTlK8dJvC5WmWn7vCwpfEQ5QWcOHel8Czd4iBQ1Fazc6L4d1A0RUiB6JEDkThp0WaqTC1xuL5efIXZ5h9dZGJF2cAkFWZyJ5AO0Csr1R60geAYNzg8R+J8PiPCJJOudggc77YDhIvf2GR5z4pfv9gVGX8mIvxYy72HHPRo9MCw5Q4esLg6AkDv2xhWRY3bzY5ebLGyVM1Tp6s8bWviTSTrsORAwZPnDB44jGTpeXX93u8qczt7qadmrg2pZp2nFZ3jOAW6821rvqjyga6uk047uJpwdqiSbNZpVYvddUP0wx0R23s5unlboorFtXicld9UHQH6t1ppi6+aqtjv78/zWqZZrm7sdD8wQcpr11+z3awLIvGSndjITudyK7tzosurrkt+tMqlWkVOzNfANSBYFfs323Y3He1ubsTLapz+a76oPkcqO6NWZ3cxXfJHaiesiRucpV8jeparat++BPOB9JMCl0wuTr0t9W0WJrenhV1B4GQgtOz9Q/SVX+2eH9lpcXa2pZ9PW1Z1vGOO+dNHBDuhoaOTwqjonVeEFS6uELuIwjXWmVylZk2vXQreMwoLj3UDlSW2sV3dZELtRRJrAEU51hfX9y2raa78PhHUNQHUytWFxxwq0OXG5UiheWJNr10KzjCCUx/hO3utJ2+a6s2VqtJZXGG2ur2Y6G4PDgGR7qiGnfTF+6j/DUKBarZLFZj+9mSnhxGC4fbQ2F1QR3sZulEki2sRpPaxCyN+e2DlRLwYO4cQjYfpFcqcucbkdKhz7XlIsXL01iN7fcV2BfFFfduG4BVaft96PLm492st1i8usLa7e0fHjxRB4P7/bjMzr+DJnWeCRv39Wd5rs6V06WOOoQDx02ig2p7LAyp85RC70IXoUsKlarF2fNVJqeb0OuAIEnS+4D/AijAH1mW9R/u+9xA1JA9BiwDH7EsK2t/9uvALwJN4F9YlvUV+/0/Bn4UWLAsa383nXVKHuu4+i5y1iI5a4lca5ECOe6IR9zYawA23dSU3Q/k9B5Vh1BtFHk8/j+Ql1YEk2h9mlx5mnpTTNk1xYHfMST0CP4UXtcQirx1kHpYHcLBQz+PyxUln8uSz0+QK0xSWpvlTprJ5Y3jCyTx2loE0xHouQ7BatbZ9Z5/Rml5iuK80CIUF7M0q+IJVjXduGIpXNEUrlgaZ2T4HjZTr3QIIz/9P6L5ApSnM6xPZSlPZ6guznEnzWTEBnGMCC2COZxG8z3IZnpUHYJsGER/4ReozE5SuZWheitLNTNBqyzOC9njxkwLLYKxJ4k+MrT9w8tD6hAG/+2vIrsdVF6bpHpNpJlqUwuikSJjpOM49gy39QhqwNNzHYIZ97H3N97H2pU58pdmyF+cIX95lmZJPEQYIWdbhxA6EMM3Fka+i83UKx3Cr/zNMyiazNS5ZabOCLrpckZoBmRVIrHXS/qIj/RRQTf1hB+8NzyqDmHskJNf/bcJZs8vc+F0hQunKlw6W6FcEudSJKYIHcIxk6cfk9i9V0Pdhs30kDqErgNCx7uRJEkK8HvAu4Fp4KQkSV+wLOvyXc1+EVi1LGunJEkfBX4b+IgkSXuBjwL7gEHgWUmSdlmW1QT+O/C7iEDSNUzJSUxKEiMJCjSsOnlrmZy1RN5aZKZ2k6maWAMwJKcIDuoAfiWKRwnRiwy8IqsEHMMEnMOksdlMtZV2cMitT7FYvAELIEkKXmfc1iKM4HMNY2iPzloBME0/Zuww0dhhmqZMo15hLTfJmk01nZs+zczESwDopg9v2KaaBlO4fPGerEfIqoYnOoonarNWrBaV/CJrK1lKcxlK81nyWbFgLikqzsgwrmgKdyyNOZxCdfSGXqn7Q+j+EL794rxvVsqUb2dZn86yPpMhf+YVct+1GTz+AI7htBCrjaTQB3rHZjLHdmCOCfKA1WpRn5sXVNNbWaqZLOvn7bHQVPTUMMaOlGAT7UihbJdm6hYSaNEQWjSE5x1H0IwGzWKZ8mtTVGwtQv6rp8j9raDdagMBXHuHcO0dxjk+jDkS6QmDR3XoBI+OEDw6Agg2Uym7TPXqJMs21XTmuVsAKKZKYHygLVjTDwfQPT3QqEgQTLoJJt0c+lASgPVcjelzy0yfW2H23CLf/uQ03/rTSQBCw452cEgf9TMw6qIXNwyHW+H4006OPy1+32bT4tZrNc6frHDxtNAjfPPvSvwO4HBKHDysceSE0CIcPqrj3ibN1Gt0s6j8GHDDsqxbAJIkfQr4EHB3QPgQ8Jv2338F/K4kHs0/BHzKsqwqkJEk6Ya9v5csy3r+roLlD38AkkZIihEiBoDkNCm0Vm0twgKrjTnm6xkAFFR85SgBLYZfj+PXomhyD1grkoTLCOEyQiQCYlGy1lhntTlHzqaaTi68wsS8uDk7jSB+9wh+9zC+QAqno0cMHs0kGNlFMCIYVVarSbEwZweILPnlLEvTrwIgqwae4AjeoK1FCI6gar2gFMo4/FH0aIzw+NsAqK8XKM1nKM0JLcLihedZePWbABj+AZyDthYhnsbw92YsFNOBe8c47h3jNA3BZqrOzdhU0wzr2RsULtiiOcPEGEniGElhJtOYw0lkozdsJn0wjj4Yx/P0EwA08mtUpoQOoXojy9pXvwVfFrMBLT7Q1iEY4yOoAz1iM7kduI/twn3MPi/qDSqZuXaAKJ7LkHtOBCrZZeDancC5N4FrfBjnrkFk89FFc5Ii494RYXDcS/rDIiFQXiqxcn5W0E3Pz3L9L85wrWnxkgT+0YBNMx0gcjCKe9DTk7Fw+nV2vTPOrnfGCaolGrUW05fXyJwRVNOrLyxx6vOzADi8aptmuvOoh/QBN7rZC0q4xNheg7G9Bj/xcz4A5mcaZM+ucfaUoJp+7P8p3tEPMrZHbWsR3vmYxtBQb6mmd6ObgDAETN31/zTw+FZtLMtqSJKUB0L2+y/ft+3QQ/e2C0iSjFcJ4VVCjBiCJVFplVi1xWq51iI3S2fAnrK51aCtRYjj12I4FE9PrG901cmAZzcDfqFobrYarK0LoVquOMVi/hozy+dgAjTVIWimnhF83hG87iEUuYscRgcIG4shPL4hhlJP0TRlquurrC1n23qEqavP0k4z+eK4B1J4wkKLYLgCj9wHAM3pwZ8+iD99EIBWo8b64jSFRUE1Xbt1kdXLgs2kmC5bqCaChDE0jNwLoZiiYA4NYw4NE3hCsJkauRURICZFmmnlG18FywJZpJnMpAgQjmQadZM008NA9XlxDRzAdewAAK1qjVp2iurNCao3sqyfvkDxBTEWstctaKY7hRZBTw32RjSnqTh2JXDsShD4sScxtBq1udU2zbR0ZZrCnwuVO4qMIx3Ft3/IZhMl0EO9mdU5wi6GntnJ0DNCH9JYr7N6ZZ7CpWkWz8+T/epNrn9WzPbNkKOtRRg96iO8O4iiPfqTs6rLpA77SR3288OI2f7SxDqZs0KwNnl2hfPPCfsPRZNI7nXZAUIECV+4Nwrz6KDK3mEHH/gxka4uFVucP1tvB4gv/k2ZT39CpK9iMfkeLcLeveq2aabXgx942qkkSb8M/DKAxsM9tZmyi7g+SlwfRTINGq06+fo8q/VZVmtzzFSuMVUWIi1DdhFwDOE3hgiYg3j0gY5+RN1AkVUC7hECbnsKbVmsV5dZqdwmXxB006XV1+xjVvC6BvF7R/B5RvBERtH13lyEhjNAxBkgMnwEgEa9QnFlkrXlDGvLWRYzp5i//h0AdKevHRw84RROf6/STDru+CiO5Eaaqbq62KaalmazrGU20kyO6HBbi+AcTKE6H30sJElCC4TQAiG8h47T0qBZLlOZylKxtQhrp14h/9JGmslIC5qpmU6jx2K9EYoZOubuHZi770ozzS5QzWbadNPyaVtAqKno6QTGWApjLIljbwLF/ehpJkmSMOJBjHiQwDMiaDeKZdav3mb98hSlK1PMf+kcc58TimYj5mvrEDx7h3Ake5RmcmpEjiVIPy4onq1mi/ytHIvn51k4P8/i+Xkmv5nlNKAaCtH9IeKHI8QPRYgdDGN4ejCTkSQiKReRlIvHPjyIT1mnuFrnxlmhQ7hxpsA3/mKWr/6JoLwOjJjsO+Zgj003Tew0kHswFi63zBNvN3ji7eKe12xaXLva4Oqpaptq+sUvCnqt0ylx5IjWDhDvOqHhfcg0U8dFZUmSngB+07Ks99r//zqAZVn//q42X7HbvCRJkgrMARHgX9/d9u529v8p4G+7XVQ2JIf1mPru7fu7xYJw+/NNPGMsLIqNZVZrc+Tqs+Tq8/d87lC9+I0hZkvCo+gdiV/C6mIabXUhlrH0jTa1eol8YYpcYYJ8YZJqrbDRb0nG50vSaFQpFmdIpd9FLHa0/XnT6HwCtLpo09BarOdmKSxmKSxlKC1P3fO54Q7hjqZZviluDgd+4jc23U9T73xRbLeQ2ygXKc1mKc3fYv12lsb6XWOhKDgH09QLq9Ryywy868fw7D6w9fd0s3i9yU9ltZrUZmcoT2SoTGSp3r53LLRwGCOdpnhSqM6Hf+M3aGmdF/2sLRan7935xqJyM1egemOiHSBahQ0WjaSpGLuT1CbnaK2VCP/TD2PuG93Yjd6ZuWJonReVdWqs35yncHmawqXbrN+89xoxE0ECB+LM/r0IXk/8xT/ZdD9urTNN1KVu3WZ9cZ3ipWlmX11k9twSlfzGArTmUIkdCjN7dpFGtcmH/+NxBvdvPcv1q521Gx75wTaNWouJy6V2gJi6ei+raXjMYM9RF1/79AoAf/DcHnxy58XpgNxZP+G9a1F5ZqbJyVN1Tp6scepUjXx+4zOPW+JtJ0xevVhlYbHV9aIylmVt+0LMIm4BaUAHXgX23dfmnwP/zf77o8Bn7L/32e0Ne/tbgHLXdingYqc+3NXe6r/6r/6r/+q/XtfrVLf32I6PsPaawK8BX0Gsuf+xZVmXJEn6LfuLvgB8HPiEvWi8gggK2O0+g1iAbgD/3GYYIUnSJ4F3AmFJkqaB/9OyrI936s9mEDqECH4pjGn6tm0raV1kybQHHykbVo0ry9/YdjOPGSXgSOAxo0iSjNWFtYDVhXlY665c6WtXP4s9hJtC011tqql2X5rJ0rp4au9ieGpUmfru57Zt4wgncMfSOEJx2MKbqZvvsrZqY1nc/sqntt1WcXlwJlIYyTSKa3tml9VFJsxSrQfea5XLLH/+89tupyeHMUfTaPEYkixtup8H0EUbSbVnES2LpY/9zbZtlYBH0E13DyPfl2bS1M6UUl3pQqtQLnL997+1bZvAviihAzE86eCWaab7ef33w6lsPYOwWhbf+K3vbru9J+pg+EiQXUfcODr4ETnkzjMa9yZP9oVckz/5dzPbbnfH+jq5Q0OSwS11nqk5u6CduiQx9W42LX71f+lOTHkHbyphmlPyWE9rH6RmVduW1zlriTVrhRbiJumQPW3L64AaxSXfa+nbCx3Ck0M/Q8vUKFQXya1PtemmlbpQOCuShs8h1gD87hF87gSasjmD51F0CKHQbmq1Emv5CVZLk6ytTlDIT2HZihiHM4w3kMQbTOELJNFDsY7shNerQ9jzgf8Jq9VkfXVW1EWwTevq60LJKmtG2/baFUvjGkii6GIseqlDcCV30igWWL8tbK/L0xnKs9PtKiVaeADHsF0XYSSNFrqXzdQLHUL8V36FpnKv7XX1VpbmmjgvJNPESI1g7BI0UyM9grzV+fiwOoT/61cxRhM0VteoXp+kdiNL5bVJKrdmwbaV0BPC9trcM4JjzzCekc4MnterQzj4b3+MVr1J4Yawvb7zqq3aGhW3Tmh/jOAd2+u9A6i2aK5nOoTPPkNkh5e1+TLT51aYOivopnOv5bGaFpIE0Z3ue7QIwYTjnrHohQ7hX/2nJO56gWuXhO21oJqWydm2Eh6fzP5jJo8fVzhyQmf/IR1zC9vr77sO4QcRumQQkYaI2ISlltVkzVoVtRGUVZYa08zUhRGWKunC9loVQcJvDaFIPWDwSDJeM4rXjDISFGNdqa/dpUWY5tbci4gZG7ht2+uAS9hem3qPTMF0F+HIXgLDYhmm1axTyN9uaxGWF64wf/s0AKruxBNMtq2vPYFhZKU3bCZXKIErlCA6/nYAyrUcxfmM0CLMZZk78zWwLJAkHMFBUTQnkcYVT6N7esNmUt0evLsP4t1ts5nqNSqz05TmMpQnMxSvXmTtrM1mcrkxh1PtIKEnh3vD4FEUjOQIRnIEfthmMy2vUM1k23qE/Bef5Q6bSU/EbaqpqI+gBnrEZgp4UR/bj//tdq2Kao3KjRnKVyaE7fVLl8g/a58XfhfOcZtmOp7AsSM6gUdFAAAgAElEQVSOrPWAPKAp+Mbj+Mbj8JPHREpiaZHl80KHsHxhjvmXxe8hKTK+XWFCB2IMHg4TORjFGe6N+7A36mDve4fY+15xv6itN7h9YZXl83NkzuY4+6V5XvrMbQA8YV0EBztAjO9TUHvAZtINif1HTfYfFQ9DlmUxnam3xWoXTlf4z9+wg6UGe/dv1EU4clwnHHmDyqvdhzdlQLgfsqTgl8L4CSO7nLbtdYFcY47VptAjLFWEz7lUlPFqYfy27XVAi/fM9trUvMR9e4n7BN21rlvkS9OCalqaYnb5PNOLNlND8wgdgi+F3zuC2xXrie21rGj4gil8wRQgTrxyaYm11Sy5tQkKy1lW58TiuCQruP2Jdl0EbyiFbHgeuQ8AuidA0BMguNO2eq5VKM1PUJrLUJzPsPLaSVqXvg2A5va36yK4BtOYoR6xmTQd58goxpjNZmq1qC0vUpnM2HTTDKWrNptJVTGGhu+1vXb1iM0UDqGFQ7hPCDfRZrNE9dYk1ZtZKjeyFF98hcI3xFgooQDGjiTmbmF9rSV6x2Zy7kvh3JcC7LGYXqR8dZLatQnWr0yx9pLNctNVnGNxURfBDhIEe/EQJeEa8uEa8jHyftv2eq3CyqX5th4h87lL3PyMmPm4Bz1tHULkYBT/6IPV1R4GulMl/XiEY0+J677VtJi7Ydten8mRPZvj/FeFult3yKQPuNl51MvYUS87jnhwentjez08qjM8qvOBn/QCoKyV2jTTs6dqfPLPSvzpH4oF6+GkwtETOk+d0DlxQmPnTrUnbKb78ZYICPdD2F57cSlehmzb61qrQq65QI4lcrU5ptYvMrEuRFpOxdsOEH4tjtsc7JHttU7IO0rIu0GvLJQX2mK1XGmK+VWh75NlDZ87YVNNk/g9w6hqb7znne4ITneEiCnkI/VqSegQVoQeYfbmi8xcF7lf0xO2La+F9bXpHeiNOEo38Q7vxjssdBlWq0lxbbZte70+myV/XThoypqBM5bEGU/hHE7jGNxIMz0KJFnGiEQxIlF8x4RorlEsUJ6y00wTGXLffp7c898ERJqpHSCSaZShHlVXczpw7N+NY789Fo0mtekZoUW4nqH62k3WXxFjITkMjB0j7boIxujw1mmm1wFJljFGohgjUcwfFXU76isFuy6CqI2w+NmXWfwrQUF2jIRx7xvCc8f2Ot6b6mq61yT2RJLYE0kAWvUmtZszNs10gdlXZsh8+SYAmlsncmCAkSNBYociRPeF0By9sb0e3O1hcLeHpz4qbK/zCxUyZ/LMnFvk+uk1/v4Pp/m7phCKDY452XlECNbGjnkJpqyejIU/IPPD7zb54XeLc71Wtbh8cUOL8MI3q3z+r0S6zO+TOGbTTE8c1zh0SMfh6ME965H38CaBLpsMyCNEzTHATjPVF1mtz5GrzbJUm2SmYpf1WzXxG3H85hABYxCfEeuNUEyS8TpjeJ0xRgYeA2CdErnCJHnb+jozvVFdze0cwO9J4gml8PlSmGaPKooZLkKD+wgN2lW0mnWKuWkhWFudYHXmCouZO9XVnLjDyXaAcIeGEWSzR4MkKzgHEjgHEoQPiTRT7Y7t9YzQIiyc/Bq8ItJMZnhwQ4swlEb39i7N5Bk/iOvgRpqpOj0tZhETGUqXL1I4bQvF3G7MZFJoEVIpjOEepZlUBSM1jJEahnfZttdry7bltaCa5j/39Y0003BcCNbGhzF3JVFD2xMpuoUW9OB7ahzfU6I+Q6tSZ/36DOtXpqhcnWLlhass/r14iFL9Tjs4iNoIrp2xXpwWyJpCeP8A4f0D8FNihlu8XbhHi/Dy74vZvqxIhO/YXtsvV6RH1dUGTA6/z+SHfkTMmCulJpnzhXZdhFf+bolvfVpQbwMRVegQjgstQnrcgdoFgaMTdEPi8DFRKe0XfkWMxcpEg1Mna23r669/3a6upsH+/RonTui88zEnT5wwiT5EmulNtagsSZK1Wz6GvI3vzGY6g3uwJcvIYr2xxmp9lvx9OoT7sdP/JLrD26G3YHWRh72fZdRs1sgXpx7QIWyGsV0/hmSL5lpdnIBWF2rGliaBZVEuLG2qQ7gfg0feh2Y+mFrphkHUqU2zVmF9YZL125l7dAibIfqef7BlaqUrNtNWP5VlUVta2FSHcD8C738/srdz+vFhWEat9TLVm5MP6BA2Q+gXPti2Q1fUzjoETemmTUucF5NLm+oQ7sfoLz2F5nlwVmcqXbiHdmijlIrMXVh6QIdwP3Snyrv+5b4tbapcXTCItmIZtZoWMzfWN9Uh3I+f/d/iRHydGUTdsIwc9y0qr65anDr1oA7hPvz/y/66jz766KOPLfHWZBnJKESkIXLWIlU2qGd36xCCzgQO2bNp8RwAOs0gAMzN89XfmvpDLFqkvMfETKJ8L8/4jg7B70zgcwwiOztPX1t6F7RTY+PxdXX1Jlcufwa3O47DESKfn6B210zibh2CL5DCcGykVrpRKndFO9Ulzn/m3wAwsPcdFBcyrC/d+/R8x9V0w/b6wZRbV7TTbTJ1hZsXmXn2r3Ak0qhON+vTGZrrGxTAOzoEM5XGHEmj+bZOM3U1o9mCdjr5b8RYeN/+dirZDLWp6Xs+v6NDMEZT6MMJcHaxSNwN7VTbeLIvvvgqq5/8Mo5Du5BUhcq1CVqFDcqk0CGMtK2v1dDGDLcbpbLZxUzDqVf59k/+IQBDHzpE/tIMxRv31qy4o0MIHozj3xVB3qRmiHMbpTKAT91e0Xvxr69z8g8vMv6eQRrVFtPnlinnN47xjg5h91E36aN+vJGt12Q2Uyrfj4C6+Qzhnzwh1gff99Mhrp8ucPPqvcd1R4dw4JjJ2F6diNG50I63C9qpT964sD72pwV++z93V8wI3mQBwcDBQVWUnGxaDdv22q6NYC2wZN2G4qvokkPQTG09glcJtf2IJKWLBTl18xt52Jmi2iiyJ/ROLIdh216vkivbdRHWp5hcPc3k6mkkZLyuOH6bZup3D2Nom6RWugoIG20cDlEQfHTHewmFxKJko1ElVxFU0/xqltzyDZbnhYWAbnjbwcEVS+P2DW7L4OkmIMiGhCc+htWsM/LYhwCR36yuLQotwkKWwmKWxYsvsHjxBbFecFeAcMXSaA438iPqEPSA8LyJ/ND7cSWFQVqzWqF8e0LURpjOULz1GoXXLgCgev2iJsKIoJoaA3EkuxjSo+gQjKRwRw3/+I/T0i3hRzS/cI/t9do3n4dvPg+qipFMbFBNd6RQPJuI5l6nDkEfEmMR+Ml3YYwmxFiUyjQnspSvTlK+OkXpzDWKL4nzQo3428FBOTAobK+3KR6ld6FDMHQN92gYM+5j9//8DCCEYqWJ5bYOYe3SbW586lX41KvIukJgrx0gDsQI7o+he02cHdZlXB3So8FRsabyjl/dQ2SHCHzlfI3pV1eYOrvC9Lllrn1zlsu2y2ww4RCW10d8pI76ie10txk8vi7SacEtAlg8qTN2yMk/+60EQblIs2mRvVa7h2r6yT/I8UnAdEgcOCyopkeP6xw6quP1Pfh7dKdD2Bi/0eTru8W/qQLC3VAklaAUJUgUsBefyLOm5YWzaXOB+XoWEDMLnxIhoEYJSAnb9ro3DB6XEcRlBBnyi0XJWmOdXPm2KJ5Tvc3U4ikmFoThq8MI4HeNEHAP43cP4zIjj9wHAFU1CIR3EgiLm6JltSgV5lhbyZJfnWBtNcvS3AW4ImipnuAInlBa6BGCI6jaoy/ESZKE6RvA9A0QHnucpiEJP6J5YXldmhcBYuH8cwAYvgjOwRTOuNAiGIEesZkME/fobtyjG2ym9dVZypMZKlOCalq4aBeu1w0ciSTmiG1aN5JENnrDZtLjMfR4DM9Tgs3UXCtQuaNFyGZYe/YF+MpzAKjRCMbOFObOlNAiDIR74riruByYR8ZwHRFECqvRpJqdswPEJOsXsxReuMACIDsNnLuHcI0ncO4dxrlrCMXRC/KAhDsdxp0OM/TBg3j1KpWlEssX51i5Y3v9/57DskVznlSA6CFBNR04GMWd6I3ttcOnM/aOGGPvEDb5zXqL9WtCh5A5m+Pad5Y5/cUN2+vkIR/pI372H3eQPujGcPTG9nrHuMGOcYMf/8ciaC3ONdrB4erpdT7++0U+ZrOZdu5S79EiJEbeeC3CmzYg3A9JkvDgx2cMMmwIlkSlVbJrIsyTa86TqZ7nVlWwJNxqEL8Wa2sRHIq3NzQ61cmAZ4wBzxiWQ6PVarK2PkuuJKimS2vXmV2xmRqKKWyvvcL62useQlF64dgo4/YO4vYOMph6EoBqOcdqaZKCbX09ffXr3GEzOX0xvEEhVnPF0xiuQG9otw43vtR+fClbNNeos740bYvVMqxlLrN6RRjDKaYTZ2xDi+AYGEZWeyOaMwcTmIMJeJtgM9Vzq4JqOiGCxMq3vgbPCTaTHhvEYdtem8k0mr83bCbF68F16ACuQwewjBatWp3axDTVG1mqN7OUz12i9G0xFrLHJYrm3KGapoa6s1zpAElVMHcOYe4cIvCjTwg200KOxo0MpStTrF+ZZv6Tz4vTQpZwpKNCi7A3gXwgjhHpTKToBmbYxdA7dzD0TuHw2qjUyV1ZaAvWJr+R4cbnBePPDDqIHLC1CIeiBHeH4NFPCxRNJnnIR/KQj3f+fBLLslieKt+jRfj7F27y94CiSoyMu9jZro3gxT/QG9vrSEzlmR9188yPugkpZUqlFhfObVBNv/SFMp/5C5H6iwzIPHZig2q6b5+G1gM20914ywSEzWDKLmJ6mpieBkR1tTU5T862vZ6r3GC6LPJ8huy0A0Qcv5zEqw/0RigmK/jdCfzuBESxba9XbB3CJLnSNDcnrwHiRu5xDeK3g4TPk8Tole21w08kECSSEJzzZqNKYWXSrouQYXHqDHMZUcBHc3jv0SK4AtunmbqFrGq4Y2ncMfv30CyqucUNLcJMlkJW/B6SrOAYGMaZ2KCa9sL2GkDzB9D8AbwHbNFcpUJ5bsJ2Nc2wduYk+ZeFUEz1+cUMIpkWNQniG2mmR4Gsa5hjacwxMRaiutpiu3BO9VaW8lm7BpWqYqSH2rbXxs7k5mmm1wlJktCiATzDbgI/LNxim8UK66/dpnRF2F6vfO0cy397kklAH/DaltcJ3PuGcCa3TzN1C9XUCB8ZInxEKIldcoV8ZpWF8wss2lTTqW9NAKAYCtF9IeKHwm3ba9PXmyJX4REn4REnJz40CEApV2fxwnybavrcp+b52p+KWUQ4YdhiNQ/HHlMZGTN7Y3vtknnbUwZve2rD9vrGtQZnT9U4c7LGq6dq/N3fiTUUh0Pi8GGtHSCOHdPxbZJmej14SweE+6FKGiFjiJAhTjzLsig2Vlitz5KrzbFan2O+egsK30aWVHxGjIBdF8FvDPakD5Ik4TJDuMwQQ+HDtAyVWn2dfGGKfEHURp6ef4XJWVsQZAbx+ZP4fCl8viR3rDAeFYpq4B8Ywz9gpxOsFqX8LPn8BIWlLIXFDCtT5wGRZnKHRjbEakOpnvRBkiTMwABmYIDgXiGau2N7facuwvLZF1g69Rwg1gycgxsBoldQTBPn2G6cY3aaqdmkOjdLxQ4QlcksxQu2UMwwMEdGMGwtgplM9qQPorpaFH0wiuftj4Peopm/Y3s9QeX6BGtfeRG+JASEaiyCuWcEc5fwRqJHbEHFbeI5tgPPMbs+Q6NJ+dY8tWsTFC9Ns/bqJMvfFIFKcRq49wy2BWvGwVBP+iDJEv4dQfw7guz6sFA0ry+ts2RrEVYvzHH2z65wuin6ERz1EbMDRPxwb9KwAC6/xuA7gxx6ZxAQtteTV0pct22vL38nx8tfWOQTgNMjs+eo0CHsOeZi16HeuB8oisTucY3d4xof/RkXftlidrbJyVM1oUc4Vef3fq9I004z7d6l8tRjJk+cMHnihPG6T4s3He1URccnBVG2imUdFqQ62QBUWqWOOgQQjCKH5oPtsr1dPD1ZmzxVWFaTQmmmow5B05x4vAkUWdt0Pw/stxtHz7v2Uy3nOuoQAJyhBLorcM9QWF08qHTTptVqUJmf7qhDUJxuzPjm1dV60ZdGLtdRhwCgjyRQAx3STF0sDG7WplVvUMtMd9QhKH43enoISVWR5S4Wp7voj3JnPxZU5/MddQgA3vEYRvjeWZ3aRX9Uafs2mtykUWmwcGllWx0CgGfAJLbHt2l1NU3uQn8hbdPGslicqnbUIQDsO2oQGtj+3qRv91132mzyXmm9xblz9b4O4Q40dAzshdFOvi/dTOvuqpDWtOqUG93Rtlx6EOnuNFMXefhublZIEuXyMq1WZ6aHprvQ9C38iLq5MW7T5VajSrW02nkniFrJHT14uvgpHuiPBbXVxbaj63ZQXB5Up2vbY3rYvljVKo3V7sZCiw7ce152059ugoZlUZ9dglYX1tR+N8oWojm5C3nPdqdys1ihtrR9wL4DT6qzH1Gn/shscrwWrNzq7jr1DJiYXh2lQ+ABumtzV39Ka02W5zrTeAFSY9o9Y6FudlwPfFcXaElcvX5PH96aOgQTJ3uUE+RtqmneWqZFkzo1FDT8UlikeNQobtnfVvHejUe1v640ChyMvJ+V5nzb1bRcFydiub6GzzFIwJnA70jgC6bQtqCw3sHD2l+P7/1HqKpJPj9BPjfBWmEaq9WgXiuhqGabauoNJHG6B5AkmabZIx2CbX/datRIvu0n2lTT4kKWmh0w6sUczoEk7rttr417x6JX9tdDH/45JEVhfTpDeTpDZXYKq9mkWSqgGGabZmqOpNHDm7OZHtX+WtI0Qh/8IOWpjG19naGZE+dFI5cXttejKczRNPruBPIWWpc2HtL+Ovq//gxWo0nl2iTV1yaoZmeg0aSZK6K4TNvyWugRtHhIpO30LhSyr8P+2hjwkP7Zt5G/OEPu4gz5SzPUlsVTdHmxRHBflOCBGKGDcQJ7o6iOewe/V/bXP/2xp6gUarb19QpzV3O0GhaFhQqGRyN1RDCJUkf8hEccm54Xj2p/nd7r4Md/KcLEmVUunK5w8VSFlSXxuy3NN9l3VOgQDhw3eepoC4dj+2v0Ie2vO25zB2+qgCAhE5EHiSDy+S2rRcG2vc5ZS6xY88yVxeKTio5fHRAvJYpPjaD2wvYaCY8ewe1IMIJwr6zUC3dpEabJLL2EhQVT4DYH8Ns0U797BIfeG1MwTXMSCu0mHBaMqrrWorh2m7xte726+BoLt88AoGoOvP4k7ohwNHUHR1B6YXstSThDQzhDQwyMPw1ArZSjuJChsJylOJdl7qxt9YyEGYzhtnUIrlgaJdgbNpPidOFK7sSza4PNVJmd3tAiXLvM2jmbweNwipoItu21MTiMrPWCzSRjJBJoo0N4f0j4ETVXV9t1ESq3suS//Cx52wJcS8Qxd6QEm2hHCjXUIzaTz40xmsB1QvhUKZSp3JyhYlNNi9+9wtrXxXmheF2Yu4fx7Bd0U8fOOHIv2EyKjHdPDO+eGMP/8CiWZVGZWyN/UfgirVyY5eofnwQLJEXCtzMsAsSBOMGDMdzxHtCIAHfEIP22COPvFmuG9XKD2xdXmT67wtS5Fc59eZ6X/9K2vQ7ppGwtQvpogKFxD6regwVzTWLXYRdvO2rxkV8S65YzkxtU0wunynz3Odv2WoU9+zZsr48e14lEvze21+3+fk+/rceQJRmfFMJHiCS2OMrRJNecF1TTxjw37theI+FRQgRqg21XU1N5dKYGgKl5iGnjxLzi5txo1ciXZwSbqTjF3OpFppeE97yuukVdBDtAuLQEci8YPIoqiuEExCKnZVlU1pftACGCxMrlq4BgM7n8ibbltSeUQjd7ZHvt8hNMH8G35y7b64VJQTWdz7By/TRLl8WCueryCpppPI0znsYR6R2byTmcxjmcJgQ0dIv68mLb8ro8laX0ml24XlEwBoftBWJRPEdx98b2Wg0GcQeDuI+LsWiVK1SzE1QmM1RvZCm+dIrCc2IslIDPFqsJPYKWHugNm8nQcO5N4txrnxetFrXbS5SvTlK5OkX56iRzJ+3zQlNw7Izj2jss6KZ7Eqi+R18clSQJR9yHI+7D+yPC+bdWqLJ6aZ7lC4JqOvHFK9z6KyEgdMXctuX1AAOHovhGA8g9YDNpDpXUiQipE2Lh2S8Xmb9ZatNMM2dzXHhW2F6rhszIAS97jrnZedTDjsMe3P7ePEQNJTWGkhrv+wlxza3lmlw6U+HGmSJnT9X4zJ+X+MTHxYwqMayI4HBC58gJnWO7lTfE9voOugoIkiS9D/gviBTWH1mW9R/u+9wA/gw4BiwDH7EsK2t/9uvALwJN4F9YlvWVbvb5MJAkCafixal4GdQFe6beqgrb64YIElPrl5lYF+wZh+JpU00DWhy3GuwN/17WCblSBB0bDJ5iZbFte71anGIhJ2oSyLKG1z1kU02T+DzDHdNM3UCSJByuMA5XmFjCLuAjlSmsTNhU0yyzt77DzI3nATBdITyhFO5oGk8khcM7sGnK7fVC0U28iV14E8KG3Gq1KK/MiroIi4JJlL8hdBmSquOMjWwEiVgSSe/NWOjhAfTwAL6jNpupVKQylW0Hidx3XoAXngNAC0UwbS2CI5lGi/RGNCc7TBzjuzEP2+dFs0ltelbYXt/IUr2eZf2kPRamLmyvd9pahB0jyI4e2V4PD2AMD8C7xXmhlFZZvzpN6fI061emWPr8d7H+WlCQjaEQzr0JAgeEq6k51JtZne4xiL5thOjbRgBoNZrkry+xfGGOtYu3mT8zS/artu21SyO8f6MugvOwG9356DdnWZaIj7mJj7l58iNC3b22WCVz1g4QZ3J85eO3+dIfiDTN4E5HW4ew86iHgRGzJ2Ph9Ss88YyLH323uN5qNYurl+qcsbUI33mxyhc/K9JkPp/EsaM6x09onDiuc+RIb2yv76BjQJDEKunvAe8GpoGTkiR9wbKsy3c1+0Vg1bKsnZIkfRT4beAjkiTtRdRX3gcMAs9KkrTL3qbTPnsCTTaIyMNENOFzbhkqa/UlcvU5VmuzrNRuM1u5DtjV1bQYAecwfnMQnxFH7ZHttccRxeOIMhw5AUClViBXmmS1fJv82iQTt18ke/t5QMLliNhitSQ+7wi6EelNmslwEYzvJRgXBXxarQal3G3WljKsLWfJzV9lcVLMZBTdgSecxBMWAcIVHEZRe6FclXGGh3CGhwjqdpqpmLMtr0VdhIVTd6WZIrE2zdQ5mEbz9kg053Lj3rMf9x6RZmpQpzoz3aaalq5epnBmI81kJlOYo0LRrA/3KM2kKMLGIpmAZ54SaaaVHNUbWSqZjLC9/sI32pXmtOEY5h0twlgKPd6bWZ0WcON7Yg++J+5UV6tTvjFL6fIUpSvTrL18jdWv2WJKnxPPXhEcPPuGcO2MIXdhv9IJsqoQGI8SGI/i/qk9WJZFabbIwqvzbS3C+T86AxZ8Q5EI7wq0tQjxwxHcAz2qrhYxOPSeKIfeIxwQzFqBzPkiN2yq6akvL/H8ZwTDyhvW2HnEw8HjBuPHXIzuc6D1IM2k6xIHj+gcPKLz8/9UzPanJpqcPVXj0mnhbPqN/yjWWVQV9u/TOG5rEU6c0AnEH/67u/klHwNuWJZ1C0CSpE8BHwLuvnl/CPhN+++/An5XElfth4BPWZZVBTKSJN2w90cX+3xDIEsKfj2KX4+Sch0SFcWaa6zaASJXn+N6ToiSxHrBAAFTaBGqjc70sm5h6h5i+j4GYoeAO7bX0+TWJskXJplfusjteVGTQNc9+HxJfL5kT57a70CWVTzBJJ5gkiHEiVeqL1NYzFK0tQi5mY00kzM4hCcslMSNcqEnRWsAdLcffdcR/LuOANCsVVmfn2B9NkNxLkPuymlWXrXTTG5vO0BYzc4Lnd1C1jQc9owAxFjUlxapTNiK5oks61ft01NRMBKJdl2EZqmEbPRGHKWGAqihAK63i/OiVa60ba+r1ycovniawtfF07sS9Akdwu4krVJnE7ZuIRsarn0juPaJp3erZcHcPMXL0xQuTVO4fJvVl8RDlKQpuMbiePYNET4YpZYvY8YfvT6DJEm4Bz24Bz2Mvl9YstQKVZYuLrJ2cZrZV5e4/LmbnP+UEHV64q52gCgu9m4sDIfCnsd97HlcHFOrtWF7fSdInPnaCiDqF4wdcoraCMdcFPOdmXDdQJIkRlIqIymVn/tHYsa8utri9BlRF+HUyTp//okSf/RHon1qROWJEwZvO2FwK/P6rpGOtFNJkv4h8D7Lsn7J/v9ngMcty/q1u9pctNtM2//fBB5HBImXLcv6c/v9jwN/b2+27T636Iu1pf5go9UjfQyCarodFEnrilL6KLCwaLU69ON12lxYD9HnVmN7xoesbnEj7OHwWFYLq7H9WMj61jfkXhGrrdr2YyF1HRS66NFW49dsYdW3v8gl8/WdFw9zKrfK2zuSKo7NZ1BdfVXHS1iMX6vRolnbno2lO7e+X0hdnRmd21TXt++D6ZK7+q5HYUfX6xa1rX+Stw7tVJKkXwZ+2f632qRx8ZF22IO7Q9Oq92I/YWDpkfrR7Fzg443GXQHjkY/nkfrR4Wb9OvFQx2JVe9qHh4ZVeeC82PZ43gglUrPcHRf/IfC6fpvaeu9mkQ+LSmnbgPG9uG66ltN3ExBuA8N3/Z+w39uszbQkSSrgQywub7dtp30CYFnWx4CPAUiSdKrbSPeDjrfSscBb63jeSscCb63jeSsdC/zgHU83CemTwJgkSWlJknTEIvEX7mvzBeDn7L//IfANS+SivgB8VJIkQ5KkNDAGvNLlPvvoo48++vgeouMMwbKshiRJvwZ8BUER/WPLsi5JkvRbwCnLsr4AfBz4hL1ovIK4wWO3+wxisbgB/HPLspoAm+2z94fXRx999NFHt3izeRn9sp1CetPjrXQs8NY6nrfSscBb63jeSscCP3jH86YKCH300Ucffbxx6B2pvY8++nZdw8cAAAQ9SURBVOijjzc1vm8BQZKk90mS9JokSTckSfrXm3xuSJL06f+vvXMJkauIwvD3m9GYiZoJWYUoONkIk42JLnwh8QFi1OjShYssshBFfCyEkI2KO10MYUCRgPgMxhhcRARBstCFCRpfGZMokwkyGjAbo4nIKPwuqnq405PpaWdudefC+aChuPXo+ruq77l1z61zc/4hSddW8nbk4yck3dNtmyWpW4+kayQdlPSDpHFJTzZVSyVvmaSvJR0or2LW95aYa0OS9kk6LumYpJsbrOXpPMeOStojqZ4dh12wWD2S1uT/xzlJY211bpD0fa6zSyq8aaiQFkmDkj7Kc2xc0pLD+yyI7Z5/SI7kCWA96Z0P3wIjbWUeA17N6YeB93J6JJdfDgzndpZ102bD9KwFNuUyVwI/9kJPCS2Ves8A7wIHmjzXct4bwPacvgwYaqIWYB0wCazI5fYC2xowNiuB24BHgbG2OoeBm0j7uD4G7m2iFmAQuKMyxz4rraVfK4SZcBi2p4FW6IoqD5L+dJDCYdyVLf1MOAzbk0ArHEY3bZaidj22T9s+AmD7T+AY6c/bOC0Akq4G7gN290BDldr1SFoF3E56ug7b07Z/b6KWXG4AWKG0h2gQ+LWwjhaL1mP7vO3Pgb+rhSWtBa6y/YXTmfRN4KGiKhK1a7H9l+2DOT0NHCHt2SpGvwzCOqD6PsIp5p7sZsrY/hc4C6zpULebNktRQs8MeWm5EThUY5/no5SWUeBZ6OK1UPVSQs8wcAZ4Pd8C2y2pnljqnaldi+1fgJeBn4HTwFnbnxTp/VyWoqdTm1MLtFmCElpmkDQEPAB8uuSediCcyhc5kq4APgCesv1Hv/uzGCTdD/xm+6t+96UmBoBNwCu2NwLngZ76rOpC0mrSleswKSLxSkmP9LdXQZW8ctsD7HIOCFqKfhmE/xMOo/WDLBQOo5s2S1FCD5IuJRmDd2zvL9LzuZTQciuwVdIp0lL6Tklvl+j8BSihZwqYst1ase0jGYjSlNByNzBp+4ztf4D9wC1Fej+Xpejp1Gb1tkqvzgMltLR4DfjJ9mgN/exMaWfLPA6YAeAk6aqk5YDZ0FbmcWY7YPbm9AZmO8dOkhw6C7bZMD0i3f8cbfrYtNXdTG+dykX0kBx81+X0c8BLTdRCiko8TvIdiHSP+4mLfWwq+dtY2Km8pcFaXiRdFF7SkzHpxZfM8wNuIT05MwHszMdeALbm9OXA+yTn12FgfaXuzlzvBBWv+4XabKoe0lMHBr4Dvsmf4hO71NhU8jfTQ4NQcK5dD3yZx+dDYHWDtTwPHAeOAm8ByxsyNqdIoXLOkVZtI/n4jVnLBDBG3oDbNC2kVYZJD5S0zgHbS2qIncpBEAQBEE7lIAiCIBMGIQiCIADCIARBEASZMAhBEAQBEAYhCIIgyIRBCIIgCIAwCEEQBEEmDEIQBEEAwH8eWTN0RPc6pwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "#4 Defining variational problem and its solution\n",
    "rho = 1.43E-7 #ohm*m\n",
    "G = 1/rho\n",
    "voltage = TrialFunction(V)\n",
    "v = TestFunction(V)\n",
    "f = Constant(0)\n",
    "a = dot(G*grad(voltage), grad(v))*dx\n",
    "L = f*v*dx\n",
    "\n",
    "# Compute solution\n",
    "voltage = Function(V)\n",
    "solve(a == L, voltage, bcs)\n",
    "\n",
    "# Plot solution and mesh\n",
    "plot(voltage,)\n",
    "plot(mesh, color='k', title=\"FEM\")\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "#saving the solution to a VTK file\n",
    "vtkfile = File('solution/sol.pvd')\n",
    "vtkfile << voltage"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ 0.003 ,  0.003 ,  0.0027,  0.003 ,  0.0027,  0.0024,  0.003 ,\n",
       "        0.0027,  0.0024,  0.0021,  0.003 ,  0.0027,  0.0024,  0.0021,\n",
       "        0.0018,  0.003 ,  0.0027,  0.0024,  0.0021,  0.0018,  0.0015,\n",
       "        0.0027,  0.0024,  0.0021,  0.0018,  0.0015,  0.0012,  0.0024,\n",
       "        0.0021,  0.0018,  0.0015,  0.0012,  0.0009,  0.0021,  0.0018,\n",
       "        0.0015,  0.0012,  0.0009,  0.0006,  0.0018,  0.0015,  0.0012,\n",
       "        0.0009,  0.0006,  0.0003,  0.0015,  0.0012,  0.0009,  0.0006,\n",
       "        0.0003,  0.    ,  0.0012,  0.0009,  0.0006,  0.0003,  0.    ,\n",
       "        0.0009,  0.0006,  0.0003,  0.    ,  0.0006,  0.0003,  0.    ,\n",
       "        0.0003,  0.    ,  0.    ])"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "volt_nodal = voltage.vector()\n",
    "volt_array = volt_nodal.get_local()\n",
    "volt_array"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "volt_array(       0,       0) = 0.003\n",
      "volt_array(  0.0013,       0) = 0.003\n",
      "volt_array(  0.0026,       0) = 0.0027\n",
      "volt_array(  0.0039,       0) = 0.003\n",
      "volt_array(  0.0052,       0) = 0.0027\n",
      "volt_array(  0.0065,       0) = 0.0024\n",
      "volt_array(  0.0078,       0) = 0.003\n",
      "volt_array(  0.0091,       0) = 0.0027\n",
      "volt_array(  0.0104,       0) = 0.0024\n",
      "volt_array(  0.0117,       0) = 0.0021\n",
      "volt_array(   0.013,       0) = 0.003\n",
      "volt_array(       0,  0.0003) = 0.0027\n",
      "volt_array(  0.0013,  0.0003) = 0.0024\n",
      "volt_array(  0.0026,  0.0003) = 0.0021\n",
      "volt_array(  0.0039,  0.0003) = 0.0018\n",
      "volt_array(  0.0052,  0.0003) = 0.003\n",
      "volt_array(  0.0065,  0.0003) = 0.0027\n",
      "volt_array(  0.0078,  0.0003) = 0.0024\n",
      "volt_array(  0.0091,  0.0003) = 0.0021\n",
      "volt_array(  0.0104,  0.0003) = 0.0018\n",
      "volt_array(  0.0117,  0.0003) = 0.0015\n",
      "volt_array(   0.013,  0.0003) = 0.0027\n",
      "volt_array(       0,  0.0006) = 0.0024\n",
      "volt_array(  0.0013,  0.0006) = 0.0021\n",
      "volt_array(  0.0026,  0.0006) = 0.0018\n",
      "volt_array(  0.0039,  0.0006) = 0.0015\n",
      "volt_array(  0.0052,  0.0006) = 0.0012\n",
      "volt_array(  0.0065,  0.0006) = 0.0024\n",
      "volt_array(  0.0078,  0.0006) = 0.0021\n",
      "volt_array(  0.0091,  0.0006) = 0.0018\n",
      "volt_array(  0.0104,  0.0006) = 0.0015\n",
      "volt_array(  0.0117,  0.0006) = 0.0012\n",
      "volt_array(   0.013,  0.0006) = 0.0009\n",
      "volt_array(       0,  0.0009) = 0.0021\n",
      "volt_array(  0.0013,  0.0009) = 0.0018\n",
      "volt_array(  0.0026,  0.0009) = 0.0015\n",
      "volt_array(  0.0039,  0.0009) = 0.0012\n",
      "volt_array(  0.0052,  0.0009) = 0.0009\n",
      "volt_array(  0.0065,  0.0009) = 0.0006\n",
      "volt_array(  0.0078,  0.0009) = 0.0018\n",
      "volt_array(  0.0091,  0.0009) = 0.0015\n",
      "volt_array(  0.0104,  0.0009) = 0.0012\n",
      "volt_array(  0.0117,  0.0009) = 0.0009\n",
      "volt_array(   0.013,  0.0009) = 0.0006\n",
      "volt_array(       0,  0.0012) = 0.0003\n",
      "volt_array(  0.0013,  0.0012) = 0.0015\n",
      "volt_array(  0.0026,  0.0012) = 0.0012\n",
      "volt_array(  0.0039,  0.0012) = 0.0009\n",
      "volt_array(  0.0052,  0.0012) = 0.0006\n",
      "volt_array(  0.0065,  0.0012) = 0.0003\n",
      "volt_array(  0.0078,  0.0012) = 0\n",
      "volt_array(  0.0091,  0.0012) = 0.0012\n",
      "volt_array(  0.0104,  0.0012) = 0.0009\n",
      "volt_array(  0.0117,  0.0012) = 0.0006\n",
      "volt_array(   0.013,  0.0012) = 0.0003\n",
      "volt_array(       0,  0.0015) = 0\n",
      "volt_array(  0.0013,  0.0015) = 0.0009\n",
      "volt_array(  0.0026,  0.0015) = 0.0006\n",
      "volt_array(  0.0039,  0.0015) = 0.0003\n",
      "volt_array(  0.0052,  0.0015) = 0\n",
      "volt_array(  0.0065,  0.0015) = 0.0006\n",
      "volt_array(  0.0078,  0.0015) = 0.0003\n",
      "volt_array(  0.0091,  0.0015) = 0\n",
      "volt_array(  0.0104,  0.0015) = 0.0003\n",
      "volt_array(  0.0117,  0.0015) = 0\n",
      "volt_array(   0.013,  0.0015) = 0\n"
     ]
    }
   ],
   "source": [
    "coor = mesh.coordinates()\n",
    "if mesh.num_vertices() == len(volt_array):\n",
    "    for i in range(mesh.num_vertices()):\n",
    "        print('volt_array(%8g,%8g) = %g' % (coor[i][0], coor[i][1], volt_array[i]))\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This is the projection of the voltage function; the soltion. This is:\n",
    "$$\\nabla u$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.tri.tricontour.TriContourSet at 0x7fc67e808278>"
      ]
     },
     "execution_count": 38,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYQAAABWCAYAAADPPOFDAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAHCtJREFUeJztnX20JEd5n5/fzOzM3bva1bKrHNv6QBIW5iDZCQiM4wTHGGwj4YD8h4lF4g8ZCZkYbBwSHBQSGzgQ44ScCBsIAYHxFxICbM4Gg+UQQTAQSSiQYK2EYn2srUUCa3e12tVd7czOnTd/VPfcmr7d090z3fNxt55z7rk93dU1Vd019ev3fauqZWYEAoFAINCYdwECgUAgsBgEQQgEAoEAEAQhEAgEAhFBEAKBQCAABEEIBAKBQEQQhEAgEAgAQRACgYVDjt+V9JikO+ZUBpN0kff5xZI+WfDcl0r6aH2lC9RFEIRAoGIkXSXpi1Nk8Xzgx4Bzzex5FRVrWt4OvKNIQjP7b8Alkv5uvUUKVE0QhEBg8TgfOGBma9NmFFkbU/3OJX0/cKaZ3VbitBuBa6f53sDsCYIQKISk75Z0RNKl0eezJT0q6QVjznmVpHskHZd0t3fuMyV9XtJRSfslvcw758OS3ivpM5KekPQlSd8p6frIhfINSc/20h+QdF2U/2ORq2UlUYb7orLvk3S2d8wkvVrSX0VleY8kecdfGZX/MUm3SDo/71xJzwTeB/xgVP6jGdfm7Kg8R6LyvSrafzVwg3f+W1LObUr6T5IOSXpQ0muj8rSi45+X9HZJXwJOAE+T9AvevXhA0i8m8nyDpEckPSzplYmvvBz4n17aC/zv877zGu+czwM/kVb3wAJjZuEv/BX6A14F3A2sArcA7xyT9uXAN4HvBwRchHvy3QbcB/wboA28EDgOPCM678PAIeA5wApwK/Ag8HNAE3gb8Dnvew4AdwHnAXuALwFvi469MMrrUqAD/A7wBe9cAz4F7AaeCjwKXBYduyIq5zOBFvBvgS8XPPcq4Is51/ILwHujOj4rOv+FRc4HXh3dh3OBpwCfjcrTio5/Hvgb4JKo7NtwnfN3R/fih3FCcWmU/jLg28D3AjuAj0T5XRQd/xjwBu/7L/C/z/vOa7zPe6I0u+bdbsNfid/4vAsQ/pbrD9gH/CXwdaAzJt0twOtS9v8Q8C2g4e27EXhztP1h4APesV8G7vE+fx9w1Pt8AHi19/klwP3R9geB/+AdOwM4BVwQfTbg+d7xm4E3RtufAa72jjWiTvT8AufmdejnAevATm/fbwIfLnj+rcAvep9/NEUQ3ppzHz8Z3x/gQ8A7vGPfkxCE/564xkUEYVuU5qnzbrPhr/hfcBkFyvIB3JPk75hZF0DSD0XujSck7Y/SnQfcn3L+2cBDZjbw9v01cI73+dve9pMpn89I5PlQIq/YLXR29BkAM3sCOJz4rm952ye8vM8H3hW5g44CR3BP10XOzeNs4IiZHU+U+5yM9Gnn+3V+KCXNyD5Jl0u6LXJRHcUJ51kZ+f01ozwG7CxYtpg4farLLLCYBEEIFEbSGcD1uCfvN0vaA2Bmf2FmZ0R/l0TJH8K5KJI8DJyXCHQ+FedempTzEnk97H2X7/ffAewt+F0P4Z7Cd3t/283sywXOzVtC+GFgjyS/ky1zDR7BuYtizktJMyyDpA7wCeCdwHeY2W7g0ziBi/NLXkOfr+Oshpg42L3q7fvOxDnPxAXGj2VXI7BoBEEIlOFdwJ1mdg3wp7jgaRY3AP9K0nOiYOtFUVD2dtzT9K9J2hYFpV8K3DRFuV4j6dxIoN4ExGPgbwR+QdKzok7x3wO3m9mBAnm+D7hO0iUAks6U9PKC5fk2cK6kdtpBM3sI+DLwm5JW5IZnXg38YcH8bwZeJ+kcSbuBf52Tvo2LoTwK9CVdDvx4Ir+rJF0saRX4jcT5n8bFHeLyP4oTr5+JAtyvZLP4/zDO7RZYIoIgBAoh6Qpc8PGfR7teD1wq6Z+lpTezj+HGrn8EFzT+JLDHzHo4AbgcF/B9L/BzZvaNKYr3EeDPgQdwbqq3RWX4LPDvcE/Hj+A6rSuLZGhmfwL8FnCTpGO4wPXlBctzK7Af+JakQxlpXoHzxT8M/AnwG1F5i/ABXH2/DnwN12H3cXGJTUSuqV/BdfyPAf8UFwuKj38GZ/ndiguk35o4/6vA45J+wNv9KuANOBfcJTiBS9bvvxasT2BBkFl4QU5geZF0ABfMLNqZbjmiJ/73mdn5uYkn/44fB37JzH6yQNqXAj9rZv+krvIE6iEIQmCpOR0FQdJ24EdwVsJ34Cyg28zsV+dasMDSE1xGgcDyIeAtOPfP14B7gF+fa4kCW4JgIQQCgUAACBZCIBAIBCKCIAQCgUAAcOucLA3bOjuss7pn3sVAZdxsg/wkE+edg22s05ZPwUeDrDwLlbvAtcjLR4Pi18caJeofn1Pmmvkkrl9mPinX2VKSWtb9SEvr72t4++TlFaeRQQMkdx0loxFtNxpGU+4mNeVvD4bbLdZpeftb0U1tCZophVvH6Ee3rB8Vbt0a9KMK9mmybhv7N7Y13DcYuO2BCYv2m2mjPZmG0/A03OeqOiSj7SmtOWU0MaXkkXp+yd98KSZ4hF87cvCQmf2dImmXShA6q3v4ey963byLQatbrtNudMu1kObJfqn041hfKX6LB51ira3f2fzDL3JNilyHInVvnjiVm8ZnfXVbqfRQ7rrB5muXdo1cvpuvcb+Tki7r/LS03vS3QWd036BjDDru3gy2GdZx96DRXqfVcdMW2u0+q+0eADvaPXa2uwDsap/kzG1Psqt1EoA9rTXOarnVNnY319jbdBOW9zRPsLfh8n1KY6OAjw26w+3DA1fvI+tucvPh9R0cXd8BwKH+To703faxvluo9vFT2znWc9vHex3Weq5CJ3ptej13b/rdpqtXr4m6Lv/GKdHoumsX/2+6qtGIihN/9ml2N+9z+ze361ZW2pPVKkFaW5mE2//oXyaXIslkqQRhWR1cg06jtChURfNkv3Dn1ugOColCq2sjHV5RgSxyHdZXWpUKIowKSFFx8MtQ5Polr13yGrl8iolBFnlikF02MegYjVNiQAPrDBj0mvSBVmd92MGutnvDjjcWhcdPbQdgV+vksNM+q3V82Jnvba5FnbwTBV8EYEMIYEMMYnY31zi6vmMoMjHH+iucue1Jt91bGZZlrdceCldMv9uk0V5nAKjbYLDNRuvcFettJwKDjhOF+JrFwpAUgzQRgHQhyBKB5G8i6wHBpyoBmIblEgTcRataiReNOjrFopQRhXmxvrqttJUQE59XxmqI70WeMOSJQvPkoNCPPss6mJa4s/Tpd5tDS+GE1+Ee7zn12dU+WVoUIFsIDkfnAMPzk+xqnRwRhTRO9Nq02/1hHdJEYZTRaxoLQ7PnhDYWhaQYTCMCWceKiMO8WDpBgPmLQr+j2jvEqkShrOsDiotC2TyLMk9BHEdRYRifx6gotLrFrIQ066AosctoHL1ei3a7z4nIQtjR7pUWBRgVAigmBof6mxdS9UXh8VPb2dV2bqvjvQ47ItHyRQHc2h1xKxvQcBZRZCVsXIO8zjhyMyV+32VEIKutx7+pLHFIfsc8LIalFASYvyiUZRK30bQd4zQdV5WiUIe7bFIrYZJ4QpJx7qQirqOqyHMXFRGCLNZ67dKiAC6eAOWFIM4PNuIIsOGyitnZ7o6IQkyv16LVWXei0Gs6t1jCx5znRvKthQ0rLb6GLq+4z0kKQZE27qfJEwf/u2B24rC0ggDLJwqzYhohqJpJxKBq66AKEcgizWoo6zpKWgnNrk3lNioSVyhK7L/3OdZfGQaaD/V3DmMAhxMuoKqEIA4uw4YrK0m73R8RBSDHhQR51kIsDM2u0e+4+xT3Ob6XINnGs9pvso3EpIkDZFsPZcShTIwKllwQYOuLQlkroSoxqMI6qFsM8qyEOoUgSTJ4X7UouM5pfBkGU7iVYvzRRnkc6e9gT2tt2MHvbq5tEoWqhABGxSAOfvuMiIIXV4ANYUi6kcYFnTdItxZiB1Ur+pzX3rMsyzRxgGKupSxxKCsEw7JMdlogkE2jO5hoqG1ZyyDPZTRp4HlWJB9kkgHMTX7s5PH8PntI/JTsDzsFRoaeTkpylFDMLMTghLcdj5aqgjLXNjOPE6eGf5uOZbT1rN9NVswyO7ZRsJDJ7y+SSNJlku6VdJ+kN6Yc70j6aHT8dkkXeMeui/bfK+nF3v4PSfpbSXdNVnTHVrYOyjJv62AWQjDuR5aXPuuvKpL1SF6LeYzMKhNHSFoH/nwEYDjqx5+XMI7dzY3jvmj458V5+fn731m0zDAqarHQxcI3KcnRR5N2tEmyfqtZv70y81pc+snKlfvLl9QE3oN7OcjFwCskXZxIdjXwmJldBPxn3ItFiNJdiXuBxmXAe6P8wL1M/bLJij1fJv1hVz1yZ9EoWr9YBIoKQR2dd1b+0wrGPEdHNVI6q3hyVhpVWgdx57/Xm7Dm7/fTQnlR8GMZRdxZSWLLKDueUI5kcDnXXZRoR2ntpAox6HcmFwMoZiE8D7jPzB6I3nZ1E3BFIs0VwO9F2x8HXiRJ0f6bzKxrZg/i3sb0PAAz+wLuxeWBHGYZJK5btCaxBubNNFZGnpUwS7dREaqyDnxR8KlDFPKshDSGs7enGImVRty+q2q3aWKwvtKYepJjFkV+/efgXjgeczDal5rGzPrA47iXmRc5d6mYxuyvc7byIo0smoa6rYFZkOc6qv37KxaJmKQYpFkHaezO2F8lWaIAG26j2EpIo8qRWVnkWQlF2klZF1HZ0WoL78OQdK2kOyXdeepk/Q1rHPOcnXs6UEYE1D2V+zdLirgEYspaCZu+awI/duPU5o4hXgvIJ230zqQsousINtxGadZBci2oRSMtZpq1FtYkQ5eLCMI3gfO8z+dG+1LTSGoBZ+Jevl3k3LGY2fvN7Llm9txtK+lT3WfBtGIwr7WMyrAMMY4ynf2sxaJK19HIsZy2V5dFUISkdbCneWI4Oc3f76eF2bmOigSXY2FIE4F4qG/cuU7rlilrJUzS70wzh6VID/AV4OmSLpTUxgWJ9yXS7AN+Ptr+KeBWc69i2wdcGY1CuhB4OnDHxKWdE8EyKM40opI3b6DKjnwW1sU0rqO8WEIR0gLK8aqgPlUO14xZBFFIY5zbKKbIfI6ik8PqePgoG3cqQ26topjAa4FbcO9uvdnM9kt6q6SXRck+COyVdB/weuCN0bn7gZuBu4E/A15jZusAkm4E/hfwDEkHJV09SQVmMaV7kRejgq0TP1gkJhWLMm6vMlbCpu8p4DZKc4mkdYiTjDJKm3vgi8BWp87h7vMcrlyoJzGzTwOfTuz7dW/7JPDyjHPfDrw9Zf8rSpU0UCvL4C5aFNJEwTqzmxWdRrM33u+dNtwybRROml9+3KqjkB1M3tM8MVzTaG/KLOa6iWculyVt+O4kzDqOVQWhFyjIolsJgcXBtxIWcdXWNNJcLb6Lpuhw0yLM2m2UJnxVzUdII2klJoVhkUfQBUEIzJSt6t4q+jQ4TRxhVqQtaJeG37HH7qK9jcHwnQj+fsi2JJaJOu7JIj00bAlBWIQ3DeWxyC6ZRS5bnejJXubfRPl5ojDOSvBFYVb+4bRhp0VIcxdlrV1UhrJWQhWkjTQaN/T0dOT07AkmZJ5uo636ZL2oVCEUk7iOFnltrjR3UdZTf9VWQtWjjepgkZ70JyUIwpIThGK25AnFJK6jolZCkeGEVQVEJ8EXgeTntBFIs7AS0kZQpY20SrMU0oL0k85DWJY4QhCEkoTgcqAMRa0EXxR8K2Gczzpv6Om4he0mwV+2ehrqtBLmyXSvVp3vcifD753Lt1ZI8+RgZmZ2q2tT+Xyn9dXXYQ2cbvGDaeMEdbGoDxrJV1hWTdVWQtZyFuMWulsEFsXSP716gwmZVgjqYlEa0aJSVbC46HcNt2c0/jztDWp5AdG0mcow+qKZrNdULgJl3pmQF0eoc+jptMzrQW3pBaGuEUaxCFQtBKfDE/mkdazSr7poFkAZJmnTRRdj8/3naYvblcWfbHZ4MFru5OciZFkJeRRd9C7vhTmLuqjdrNgSvVOVorCo1kBMbBVUYR1sVXFaZjFIMs1iatOu9e+/xtJ/1aX/OsyixDOWJyVroloR8pbmyLtOaZbYIg91r3txu9OCuoUgfr3kMqx6WoR+Rwvp995KYpBkmh96GvPyq49bwqKs2JRZGtuvb95Io6rnIiSXNslbyHFebBlBWFTFrkMETkfroIhfft5ikPX9Wa6wrHaxyHMRfI56HXvSAijqLjo6RhyqGNU0bRzBdyFVLciLyHL1CkvEVrIGktRlGYx7aspbPG7eYpAkS8Cyhp5mWaf+sNNpljWeZwB1WncRjLqspmXaOILvxlsEK7mKV2fGbClBmMZKqOrGbmUhKMOg08i0QtZXWplWzqSiYNsXKxrol9Wvk19v//r47c9vx/6P3X9C9f3aixQIPTxobLIO0l6aE89F2N1cGzv0NA4s72qd3BRHGL7nuX1y6Dra2e4O/yZ9o9rpzJYbt7i+0pjY5O53tLDvTPZpnuxP7TaKyzqJ66jMNYrzT7s26yut1CfmuANNc7XEHW3qEtTb27OzFJ5MDHncXt0TbFU0umLQsZEJagMaWGfAoNekj/Or93ot2u3+cOjpjnaP470OO9vdYWD5zG1PDp/S97TWONTfyVmt40OXz97mGkfWV8e+E8E/dmR9dSgKh9d3ZL53+VB/51AUjvR3jIjCsf7KUBQeP7V9JJ5wrLeyacJavBT2oNccDr9tnNLw+sT/m72NGd/N3sYEwNhCa3U33Hqtrg3btt+Wk+1T3VNzXyK9CFtOEKZlGlEYdBq1i0Lc6KoQBZhOGMqQdW3iOkwiDFmiABO6kJKd/AToyZ4TppQOoG7roMj7gBuntEkUYsqIQkyaKIwjTTD8mcu+OBxd35G5kF5SHHxicfBHSfmxhD4wYHRORlVi0Dyx8RIlvw0m20TzxKmFDCxvCUGoOgi3qKKQ7DSrEgWoXxjyrkmeMExiLcAYYaig89+UX4aVEP/48+5XFWLgk7qSZxRLUHdDFBrt9eGcBF8UYtJEYVfr5DDo64tCTNbyFHmCkWU1+OLgWw0wKg6+1eCz1tusktYZMGCzKEBFYhC3se0rIw8K7ru3jYhClb/laShUAkmXAe8CmsANZvaOxPEO8PvAc4DDwE+b2YHo2HXA1cA68CtmdkuRPFMZzG4ExiKJwrg1cKpuSFULQ9nrkCUM07iRIEUY4s67SmGIRGGclRDT6A5qEd8066DoqzSBEffRars37Eh9UYhJE4XdXoc+KXGMwc8jy2oY51IaRx/nOvJx12k0lhiLwsRiACPtAkathbg9Jx8Y6mofeeT2JJKawHuAHwMOAl+RtM/M7vaSXQ08ZmYXSboS+C3gpyVdDFwJXAKcDXxW0vdE5+TlOXcWQRSKLKlbx9OFX/ZJGmZdVlLl8YWSwmAn0idFaTV9zR//h59lJbS6Rr8jmicHmwZG1GEd+IwbcZMlCrvaJ4drHPmiUJasuME4McmyGpIUEYfYdRRbCcnFADeusQADRuOTeWIQtxWtbi9sLcR5rq+0RkQhbiN1U6QXeR5wn5k9ACDpJuAKwO+8rwDeHG1/HHi3JEX7bzKzLvCgpPui/CiQ50KwCKJQBL8hVU1Rq6GquhYRwUmEYawbqYAwZInBCAkrIY1YFMY9Bba6zm3U7NrY8e+DnNFGRa2DtElqyTH8fkwhTRTiYLPPuJfpjJuDMAlJl5IvVFnrHiVFwREHmRNpO+AGZg5o0QBGB0WMxAy8trJJGEpYC/5vzxeFtIeHKijSe5wDPOR9Pgj8QFYaM+tLehzYG+2/LXHuOdF2Xp4LwbyXscgaiZNFnb7ILGGYlSWURlZ8ASYMPG9fyRQFrW7PFAU78eQmKyHpOkoLJMaiMM5KiGl205dRSDKNdZC2zMO44ZvJ1VCT6w9NssxFzKRvZks7L2vV1lZnfawo+J/HzQMZPoR4bWqwthFAb+xY3WgjU1gLcX/kt5V4rkoV8xHmH8XIQdK1wLXRx+6XP/GGu+ZZngo5Czg070JUyFaqz2R1OVp9QSoi3Jt5czTxf4NZ1Of8ogmLCMI3gfO8z+dG+9LSHJTUAs7EBZfHnZuXJwBm9n7g/QCS7jSz5xYo88KzleoCW6s+W6kusLXqs5XqAotXnyJOqK8AT5d0oaQ2Lki8L5FmH/Dz0fZPAbeamUX7r5TUkXQh8HTgjoJ5BgKBQGCG5FoIUUzgtcAtuCGiHzKz/ZLeCtxpZvuADwJ/EAWNj+A6eKJ0N+OCxX3gNWa2DpCWZ/XVCwQCgUBR5B7klwNJ10YupKVnK9UFtlZ9tlJdYGvVZyvVBRavPkslCIFAIBCojy212mkgEAgEJmdugiDpMkn3SrpP0htTjnckfTQ6frukC7xj10X775X04qJ51knV9ZF0nqTPSbpb0n5Jr1vWunjHmpK+JulT9ddi5HvraGu7JX1c0jck3SPpB5e4Lv8iamN3SbpR0syWbp20PpL2Rr+PJyS9O3HOcyT9ZXTOb0uayUsLqq6LpFVJfxq1sf2S8pf3mRYzm/kfLpB8P/A0oA38X+DiRJpfAt4XbV8JfDTavjhK3wEujPJpFslzyerzXcClUZqdwP+bRX3qqIt33uuBjwCfWua2Fh37PeCaaLsN7F7GuuAmij4IbI/S3QxctQT3ZgfwfODVwLsT59wB/H3c7LLPAJcvY12AVeBHvDb2F3XXZV4WwnA5DDPrAfHSFT5X4H504JbDeFGk9MPlMMzsQSBeDqNInnVReX3M7BEz+yqAmR0H7mFjlvdS1QVA0rnATwA3zKAOPpXXR9KZwD/Cja7DzHpmNotpabXcG9xow+1yc4hWgYdrrkfMxPUxszUz+yIwMq1c0ncBu8zsNnM96e8DP1lrLRyV18XMTpjZ56LtHvBV3Jyt2piXIKQth5Hs7EaWwwD85TDSzi2SZ13UUZ8hkWn5bOD2CsucRV11uR74NdwqAbOkjvpcCDwK/G7kArtBUrUL86RTeV3M7JvAO4G/AR4BHjezP6+l9JuZpj7j8jyYk2cd1FGXIZJ2Ay8F/sfUJR1DCCovOJLOAD4B/KqZHZt3eSZB0j8G/tbM/ve8y1IRLeBS4L+Y2bOBNWCmMauqkPQU3JPrhbgViXdI+pn5lirgE1luNwK/bdGCoHUxL0EosxxGfEHylsMokmdd1FEfJG3DicEfmdkf11LyzdRRl38IvEzSAZwp/UJJf1hH4VOooz4HgYNmFltsH8cJRN3UUZcfBR40s0fN7BTwx8A/qKX0m5mmPuPy9N0qs+oH6qhLzPuBvzKz6yso53jqDrZkBGBawAO4p5I4AHNJIs1rGA3A3BxtX8JocOwBXEAnN88lq49w/s/rl/3eJM59AbMNKtdSH1yA7xnR9puB/7iMdcGtMrwfFzsQzsf9y4t+b7zjV5EfVH7JEtflbbiHwsZM7sksviTjAr4EN3LmfuBN0b63Ai+LtleAj+GCX3cAT/POfVN03r14Ufe0PJe1PrhRBwZ8Hfg/0V/tDbuue+MdfwEzFIQa29qzgDuj+/NJ4ClLXJe3AN8A7gL+AOgsyb05gFsq5wmc1XZxtP+5UV3uB95NNAF32eqCszIMN6Ak7gOuqbMOYaZyIBAIBIAQVA4EAoFARBCEQCAQCABBEAKBQCAQEQQhEAgEAkAQhEAgEAhEBEEIBAKBABAEIRAIBAIRQRACgUAgAMD/B54DczKDxDI1AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "V_g = VectorFunctionSpace(mesh, \"Lagrange\", 1)\n",
    "w = TrialFunction(V_g)\n",
    "v = TestFunction(V_g)\n",
    "a = inner(w, v)*dx\n",
    "L = inner(grad(voltage), v)*dx\n",
    "grad_u = Function(V_g)\n",
    "solve(a == L, grad_u)\n",
    "plot(grad_u, title=\"grad(u)\")\n",
    "grad_u_x, grad_u_y = grad_u.split(deepcopy=True) # extract components\n",
    "plot(grad_u_y, title=\"y-component of grad(u)\")\n",
    "plot(grad_u_x, title=\"x-component of grad(u)\")\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The previous code can also be substituted by the method project, like:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "  Calling FFC just-in-time (JIT) compiler, this may take some time.\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "<matplotlib.quiver.Quiver at 0x7fc67efcbcc0>"
      ]
     },
     "execution_count": 29,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYQAAABJCAYAAAA9vFENAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJztnXd8HNW5979nd7V91XvfpmpbXbuqljEGQkIg3FwguTfvzQ2ENz33pgNJwCSkXW4SSAcDoYUSejHVNs2AG2Djgptsy92yZXVpV9LO+8cuRitmzWq0Wtsv8/18/LH2nJk5Z3ZnTnme33mOkCQJFRUVFRUVzamugIqKiorK6YHaIaioqKioAGqHoKKioqISQu0QVFRUVFQAtUNQUVFRUQmhdggqKioqKoDaIaioqKiohFA7BBUVFRUVQO0QVFRUVFRC6E51BaZDenq6VFxcPO3zJgISw6N++UwJDh7tIxAIX7GtT9BiNRuwmg2YjHqEgvoCDI/4P3Tt9xkYGqVvYDQsTaMRWEx6rJZg2Vqtsj57fHyC0dEx2TwpIHHoYB9TV6kbjAlYrQYsViNGY4KicgGGhn1EWgB/vH+YoZHw30KrEVhC37XVbECjUfZtj/nH8fnGZfPGxyc4crg/LE0ARlMCFqsRq9WI3qD8dRgc8kXMO3psEJ8v/LfQ6bRYLAasFgMWsx4hlN3z6Ng4Y+MTsnkjY2N09w+FpQkhsBgSsBkNWI0G9DqtonIhANJg5OyJAyBN+S2EHoQNNDYQFlD4Vg37xxgPBGTzBn0+eoZGwtI0QmA16LEZDdiMBnQahe/URICRCO8UksSB7v4Pv1N63aR2ZAbv1IifQISXqm9ghIHh8OdPqxFYTAZsZgPbt246KklSRjTliDMpdEV9fb20du3aaZ+3eechvvTTexWXm2gx4q0qprXGQWutE7NRH/W5X736H7y79YCicoWACncOLfVO5nvcFOWnRX3uy8s28/NrH1VULkBqmgVPsxtvs4sGr4uEhOgbjk9f+VeOHR/66ANl0Go1VJfn01LvYH6jm+yMxKjPffDulSz50zJF5QJk5ybjaSmhqa2E6np71B3T2NgEZ1/0v4rL1et11FUX0dzopK3JTUqyJepzr390GQ+8uUFx2c7MVDoqHJxV4aSqMCfqjkka70Q6ep7ichFWMLQhDAvAsBChsUV96pfveZRXd+xWXPSc3CwWlDg4u9xFaVZ61Oe9+c4uvvOLRxSXm5xoornGQUudg+ZaBwZ99AOQS35wB3sOHldU7up7vrtOkqT6aI5VTUZRkJeVRFFOKkW5qZgMynv56ZJkM1GUl0phbioZadG/MDNFoxXkFaRRUJBGQVE6Ol38HpP0ZAuFuSkU5qaSkmSOW7kJCVryC9MoKEojryBV8SxFCVkZNgrzUinMT8VmM8WtXJM+AXtmKvb0VPJTkxTPUhShLQCtA7T20GwhPiQZDdjTUrCnp5CTaI1buUJAQU4KhXnBZ1s/jQFWPPlYzBBG/WMc7O6XzRsbm+DrNzwYNuUy6HU0zimitdZBc7WDjBTlD87BI334/PJmjIeefovHnl8fluYoTKel3klLvZNyV7Zik9HQkI+j3QOyeX29w3z/2/cQmPjgt7dYDTR4nHhb3DR4nSQmKm+Y9h48zsSE/JT+z/e8wsp1nWFpFa7s4D3XOXAVZShumPr7huntkZ+Z7O7s5ufXPBSWlpxiwdPqxtNSQl2jA5M5+pnfZCRJYs/enoj5N9z4FNt2Hj7xWasRzK3Mp9njornRSUFeqqJyAbr7B+kfkTdXrdy2h18/9XJYWk6yjY5yBx3lDhqd+eh1ysxkkuSHia4ImQGk41+CwJFJiXowNIVmBAsQ2hxF5QIc6O1nZEzedPPw25u4/fV1YWmO9FQWlNhZUOqgOj8XncJ3anjUz5Gj8u/U0Iifr1/3IP6xD951szGBxqpiWuucNNXaSUlUPsDZd7g3omnwtsfe5IU3t4allRZn0lbjoLXGSYUjO+oZwhnlQ1CKUZ+APU/e3PLIi+sZGPaRkWI9YRKqqyzAqI/NTCAnM0k2fXjEz4o3tpGg01JTWUBLvYPmemfE46eLxWLAYjHI5i35y3ICExK5+Sk0tbjxtriZM68AnWJ7cjgFOSmy6Ye6+3nznd0YDToa5xXTUu+gqcZBWkpsRoiJSWYSI8wq7l4SbBgdriw8rW68rSWUVuTFZCYghKC4UP752rh5P9t2HsZmNeKpt9PU4MRT58BmM864XICMRCsZMiNdSZK4+sHnEALmFeSc6ATc2WkxmQkIoQedSzZPGnky2Blo0sHQgTCcBfpmhCY2M77cZHkzon98nKff3YpOo6G+KI8FJQ46ShwUpSXHpFyzUU9xBLPtnY+swj82Tk5GIi11TlrrHVSX56NPiE0Tm58lfw89fcO88tZODAlaGiqLaK1x0FJtJzNVmUXhY9EhnIy1W7r4xufbuWjBPKxm+QZ0NnjpjW14a+386yfrKHNmxa1cv3+czl1H+PYPz+ecc+bOyJE6XZ57eTMXnjWHSz5VR0Gu8lHxdOk5NsjgkI+rb/gX2uaXo42jCWzFS1v4t880cPFF9aSnx8/st2n/EdITLfzlnM/Q4i6Kqwns0PDrCNP3SbNcQoIuNgOcaFi2awc1Fdl8oaaW2qzcuJnAAoEAW9/dx/c/P5/zzp2HyaRspqmE51Zt4dyWMi47txZnXvT+kEh8LExGJ+OFN7fy4z89TYJOS215Pm01DlqqHeRmzO6D7B8b57Jv3sbhowMU5KTQXOegpd5JVVlezEbqkbjv3tdZcutLGI0J1NXbaWp24/E4SU2bXZvq8b5hLv2/tzAyOoazOIOWBifNDU7KXTmz3mD9/jdP8/Tjb2OzGWnwOvG0uGnwOLHNwDQWDbs6j/Dly5cAUFaeR1Ozi6ZmN3a7ctNYtHz5jkdYuWMP6VYz80sddJQ5aHIVYo7R7DcSB4dW8tqhbyHQkmGqJcfcRq6lHWtCwayWOxEIcM7Tt9LZ30OBNZmz8pwszHPjySxEr53dd+qJh9bwxxufQW/QUVNvx9NagrfFTXpm9KIIJQyO+Ljwh7fRNzSKIzeNtioHbVUO5jpz0IbUVEKIqE1GH4sOYWBolA3b5ZU+gYDE9bc8S/8U2aAzP43WGidtNQ4qnNknvtzpsnHrAfqHRmXzlq18j2df3hyWZrMY8FTbaal34q0uJlGhk/HYsUF2bD8kmzc87OeXNzzxITt/WXkuTU3BBsvhzFTcYK1dv4exMXl750NPr2P127vD0lKSzDTVO2hpcFJfVYxZ4QjrwL4e9nbJ2/KPHO7j5hufCUvTaAVz5xXibXXjbSkhv0DZrCUQkFizemfE/NuXvMyOHYfD0rKykvCGvuuq6kL001CcTGb74aMc6JW3a7+77xB/Xv5mWJpep8XrKKSjzE5HmYPsJGWzlvHAMN0jb0XIlVjb/XNGJ46GpdoSism1tJNjbifNOBeNUHbP648doGd0WDZv2f4d3Lv97bA0q05PW66DhXlOOnJdpBmVma56jw+xbYt8O+L3T/Crax/BP0X27CrNCT5fbSW4SpQPfNZt3cuIT95v8virG1nx1o6wtCSrkZa5dtqrnCxqLFU7hMls7jzEf177D8XlpthMNFc7aK2x0zTXPi098VeuUS471WoEc8vyaKlz0tboimibl+OlFZv52eLHFJULkJmZiLfZhbfJTW1t8bRkpxf9558Vy071CVpq5hYGncyNLjKnYWZ54N7XWfLn5YrKBSgoTMMT8qnMrSqcluz0vEW/VlyuyaSnvsGOt8lNU5OLpOToG6zFjy/jgdXKZadlORksKHOwoMxBZV5W1IOAAf9unt37L4rL1WuSyDY3k2tpJ9vcQoImej/SF5c/wMsHOz/6QBkEUJuex8J8N4vy3biSojezrHljB9f8t/J2JC3DRmOzG2+rm7pG57TMtf9yzR3sOaRMdrrudlV2GlPSki1kpljJTLFNSzs8U8xmAxmpNjLSrCTFUY4oBKSl20hPt5GRYYur7NRqMZCeaiU9zYrNGj+fjlarIS0jeL/pGba42tyTkkykpwfLVap0UoJBpyU70UpmopUMmzWuslOjNg2TLhOTNgOdiI2DPRqsCQayzDayTFZSDfGTNUPwncrItJGekUiC/gyWnQohzgNuArTAEkmSfjUl3wDcBdQBx4BLJUnaHcq7CrgcmAC+JUnSc6H024FPAUckSZoTTWUVy059Y3Qd6pXNG5+Y4Fu/fjhMdpqg01JfURDy2DvISVduB9x78DijEaZ6jz63nsdfCJedFuSm0FoXlJ3OLctTLJEbHBjl8OE+2bz+/hF+8L37wlZQG40JNDQ68Da58HhdpMxA+bOr62hE2emt977G62vDzStueybNDU5aGpyUOrMVN8bHjw/Rc1R+9WzXnqP8YspCPVuiicYmJ95mN/UeJ1aFyp9AQGJX55GI+f/zm6fZvu0D850QUFGZf8I8V1ScrrgxPtQ3QO+wvEnyjZ1d/M8zr4SlZdgsdJQGzUVeZyEmhf6EiYCPgbE9snkSE7x28L/CTEYCHRmmOnLNreRY2rEm5CsqF6Br4DiD4/KRBx7p3Mht760OSyu0JrMwz8XCfDcNGQWK/QlDQz4OHZAfpY8O+/nht+4JMxkZDDpqGhx4W0vwtLhJy1AuKth9sAd/BNnpHU+v5oU14bJTZ14a7VVO2qocVLnzYic7FUJogT8Bi4B9wBohxBOSJE02fl8OHJckySWEuAz4NXCpEKICuAyoBHKBF4UQJZIkTQB/B/5IsCOZVYyGBEqK5FduP7ZiAwPDPlISzbRU22mrcdA4p2haq5FPRiQzz8ion5dXbUOrEcwrz6e13klznYPCGKlvrDZjxAbutiUvEQhIZGYl0tTknrEdeyr2Qvlp+OHufla/swt9gpbaeUW0NDhpqnOQNY3VyCcjJcUSsSN74N7XASgoSqOpJejwq5iTHxPFkUYjcLrklWJbNu9n+7ZDmM166hve73CdJE9jNfLJyE6yyfoBJEli8ePBVdsVuZl0hMxC5TmZMZn9aDUGkg0lsnl7B59ndOIoek0SOeZWcixtZJu8JGhjo7IqtMm/U/6JCZ7peg+NEEGzUJ6bhfkuXImxkdpaLAac7mzZvAfuWonfN056hi3oUG51U11nxzCDcBWTKc6RbxeODwzz2oZOdFoN9WUFJ5zKuenKRDHRtACNwA5JkjoBhBD3AxcCkzuEC4HrQn8/BPxRBH+BC4H7JUnyAbuEEDtC13tDkqRXhBDFimodQ9Zs7OJbn2vnM2fNi1knEA3LVm6lucbBJRfU4S7OjFu5fv84O3ce4bvfO59F584hIUY66Wh47uXNXLCoiksvrCcvOzba8GjoOTbA4PEhfnLdRbQurECjUCCghBUvbORzl3q5+JLGWVdxTWbjvsOkmU3c+h8X0+wujKs5qG94OQ1pV5Jv+zy6GHUC0fDi/q00ZufyxZJGqtLy4lZuIBBgy4Yuvvmd8zjnwhoMcYxmsHT1Fs5pLuPfF9Ziz4k+tE0kPtJkJIT4LHCeJElXhD5/AfBIkvSNScdsDB2zL/R5J+Ah2Em8KUnSPaH024BnJEl6KPS5GHhqtk1GJ+PZlVu49q/PBFcnVxbSUu2gtWZmq5Ojwecf57KvLaG7ZxB7QRotDUEzUYU7R/Hq5Gi5956V3H7bK1gshhNmosZGB0mzHCqip3eIS/7vrfj845S7smlucNJc78BtV65oipbfXfsIzz2yjuRUCw1tpXg6yqhrdmGa5bUnndsP89XP/xWNVjCnqhBPWwme1hIKimeuGf8oLr/lIVbt2EtOso355Q46Khw0OPIxzPIgoG/4BXYf/RICA1ZjE4mms0k0nY1eN7uN9HhggguW/Yk9Qz3YrWl0ZJeyILuE6tQCdJrZtdk/ft+b/OXXSzGZ9dQ2ufC0l9LY5iY5dXbbkf7hUc5ffDsDIz5K8zKYP8fB/DkOKgqyTswEYyo7PdUdghDiSuBKgMLCwro9e+TtlidjcMTHls7DsnkTgQA/+dPTH5KdltmzaKt20FLjoKxYeYP13o5DDA7LhxZ4/uXNLF2xKSwtOdGEtzYoO22sLsaisMHq6Rli965u2bzhYR/XL34szM6v0Qgq5wTt2t4mF4WFyqfZ6zfvi7jM/v7H1rBqiuw0M812QnZaO6dA8Qjr8IHjHIggOz20v4ebrgtXXSUkaJnX6MDTHuwgsnKjV3FNJhCQWL92V8T8W296gZ3bwiXAeYWpQdNCeymVVcpXiXce6eFIn7zf5J09B/jj82+EpZn0CTS5C+mocNBeZifdpsx0NREYZtgfWXbadey/GJ8I96sYEypINC0k0bQIs74KIZQNfLb0HqR3bEQ2b9mBLfxj15qwtMQEI21ZbhZkl9Ca6SJRr0yg0d87zM6tB2Xz/L5xbvj+A2E+BCEEZXPz8bSX4JlfRrFLeTuyftcBRiOEwHngtfUsWx8uO01PNNNWEewczqpyxbRDaAKukyTp3NDnqwAkSfrlpGOeCx3zhhBCBxwCMoAfTT528nGhz8XEYYYwU9lpRoqFlmoHbTVO6qcZ1uKrVymXnep0GmoqC2iud9LaML2wFi+t2MLPrlcuO83NTaap2Y3X62LeNBusi770F8WyU6NBR/28otDswTmtsBb/vP0Vbvvdc4rKBSh2Z+HtKKOxvYyyeflRm5bGxib4VPPPFZdrtRmpb3LhaSuhodk1rcVy1z+yjAdnEO10bkE2HRUO5pc7KM2J3rk9OraTrQc7FJer02SQaDqLRNNCrMb5aKcR1uLK1+/htSM7PvpAGbRCUJdWREd2CQuySymyRm9mWbtyOz/++t2KygXIzEnGO7+UxvYSqhrs0zLXXvjzv7P7iDLZ6YY/fCemstM1gFsIYRdC6Ak6iZ+YcswTwH+E/v4ssFwK9jRPAJcJIQxCCDvgBlZzhmExGbCagnHrdbO84nEyCTptcG8Esx5jHO2SABaLMRgPyWqIq83doE8I7hNgNsRV4qvRCCw2IxarEYvNGF8JpikBszX4XcfKsR8NOq0Gq1F/4l8871mjsaDR2NBoEtGI+D3beo0Oq86ATWfErIufzxDAYjNgtgT3G5lts7BSopWdng/8nqDs9HZJkm4QQlwPrJUk6QkhhBG4G6gBeoDLJjmhrwG+BIwD/yVJ0jOh9PuADiAdOAxcK0nSbSerh9IZwvCon137j8nmjU8E+N5vHwszGWm1GqpL8oLRAmsdFGQpMyUA7Np7NOKmGk++sIEnX3w3LC0nM/FEtNOqSuXBsfr6hjl4QF5qOzAwytVXPRgmO9XrddTUFtHU5Mbb5CJjBhK5bZ2HGR+Xl53+/cE3eGNKtNPigjRa6oP+hMrSXMUvy7Hufo4ekpfa7t19lBuvDo92arYYqGtx4+koo6G1hCSFUttAQIq4ghXgD798ih1bw01GJeW5eNrceNpKcZVmK26M9/f00TMobz5Z07mX3y59LSwtxWKivcxOR4WDJnchVqMyk2QgMMLo2HuyeRISu49+eYrJSIPFUE+iaRGJprMx6JyK73nXwFEGxuWltk/vfZe7O1eFpWWbEunILqEjuxRPejEGrbIOaLB/hP175NuR0VE/P/lGuOw0IUHLvAZ70CTZXkpWrnIhxbb93fgimGHvWr6O59/eFpZWmJF8wp/QWFIY22inkiQtBZZOSfvppL9HgX+NcO4NwA0y6Z+LpuxYYDbqqXTKh9t98pWN9A/5SLQYaK6y01rjxDu3CJslNotl7AXyzsNR3xivrdmJEFBZkhvqBBzYC5Vr0ieTlGSO6CS+4/ZXCAQk0tKseL0uvM0uamuLZ7RL2mRKHPISzO5jA6zbsCfY4VbmhzoBJ3k5sVEcpWUkkhZBwvr4P4IhHHLyU/F0lNHYXsrc+uKYqKw0GkFZpbzDdNvmA+zYeiioSW90nHAqz0STPpm81CTyUj9sSpQkif8Jhb52ZaUxvyIY7XReofIwLJPRaEyYDTWyeb3DSxmfOIJG2LCZ5pNkWoTNuACdVvnAajJ2m/w75Q+M8901wU5/bnIuHdmldGSXUJakvMOdjDXRROlc+fUTD935Gn7fOEkpFhrbSvDML6XW68QcIeLwdCnJk5fN9w6NsHLzbrQaQY0jj/Y5DuZX2inOUiZf/9hHO31zw26++bl2Lo637PS1rXhq7Vz6qTrc9vjKTrdvO8R3vvcJzjln7rRCUsyUZ1Zs4vyFc7nswjrysmPTOERDz9EB+o4PcfX/Xkbr2ZVxNYEte+odLv33Ji76tyZS4xntdN9hksxGbrniMzS5i+JqDhodeY6C5B+QbP0immnshDZTXjr0HjVpefy28bPMTVa+8G26BAIBNr61h6/96HzOvaguZmsPouGp1Vs4t7aE/7OgFnt2HGSnpxOzITtd+tpmFv/tWczGBDxzi2mtDsYTn8lmFtHg841xydeXcOz4ECX2zBNmohJH1qyHTbjrrte48++vkphowuNx4m1yUV9vx2qd3RACR3sGueRrSxgfnwjNioJ7QNgLYrNw6GT8748f5oXH3yIt00ZjexnejjKqPc5Zf3l3bD7AN//lD+gStMxtsONZUI6no4zs/NkP//3Fv/yTtZ37KEhLoqMiKDutteeRMMt+MN/Is/QfvxyECb2hDb1hEXrjQrTa2Q3zPh6Y4OJXbmLfcA8uWxbtmWW0Z5ZRmZyHVqGiKVoevWslf/v1UqyJRupaSvB2lFHfVoItaXZDzvQPj/KJ625jyOdnTmE28+c4aJ/joCT3A0uDGu10CkMjfrZ1yYcWmJgIcNXNT4b5EISAOc4cWmuctNbYceYrN+Ns332E4WH5ZfZLV2zk6eUbw9LSUiwnQmHXzytS7Ezu7R2mq0ve3jk0NMq1P30kTHaq1WqoqioMRuD0usjNUz6C37TtIOMR7J13P7KKN98Kl2jmZCWd6BCrK/IVz1q6D/VyeL+83+RA1zF++9Pw/XANxgSqvc6gjXd+GWkKQxUHAgE2vxVZDv3XXzzFzik+hiJXFp4FZXgWlFM6r0Cx36TraC/d/fKy03W7DnDzsyvD0qxGPS2lxXSUO2grKybZoqzBkgIjjI1FUjdJDBz/GoFAuNRbl1CN3ng2euMidLpKxe/UjoHD9EeQnb54cCP37wmP8Jqit9CaUUJ7VhnedBcWnTIzzkDvMHt2yrcjvtExFn8z3Ieg0WqorC3CMz8oay6wR7XPvSybug7hG5OXnf7j5Xd44Z3tYWk5KTbaK4M+hNZKu9ohTGamstOc9MTgbmo1DmrLpufo/erVymWner2O+rmFNNcFR9KZ09hX+aWXZiY7LSxKCzqYvS4qK/Om1WBdeLly2anZpKexupiWeidNdXaSpzFTm6ns1F2Ri6ejDE97Ga6K6DdYGfOP8+mqnyguNzHFQkN7CZ4F5dS2uLFMY6Z2/cPKZacaIaguDu6mNr/CgSMzNep7Hh/fwfEj8xWVC6DR5IQ6h7PRG1oQIvqO6Rtr7uT17u0ffaAMCRotdan24Owhq5RcU/QDn7WvbuPHX7lTUbkAeUVpwedrfhmVtUXopjHwOZ1kpx97NBqBRiPQajRxtcUKgotbNBoNmnhugA5oNZrQPce3XI0IlqnViLjfs0arCX7X2jiXqxEflB3P50uARgTLjfezjdAAWgRagk96nIpFoBUatEKgiWO5EHy+hEag0Qri/GhHzcdihjA04md7l/yq3YlAgB/d9ESYyUgjBHPcOSc2qbbnRj9ymsr23UcYHolgMlr+YZNRRqqV5pDiqG5OoeJVu729w+zdK28yGhwc5ac/eThMdqrTBU1G78tOc2ag/Nm8/WDElcr3PLyKN6aYjPJzkk8ojqrKle8Y132oj8MRolEe6Orhtz95OCzNYEqg1us6oTpS6vQNBAJsfjvChvPAX2948kMmI3tpNo0dZXg6yimdxiK4qXQd7eXogPxsbN2u/dz0TLjJKNFkoLW0mI4KBy2lxSSZlfmNpMAI42PvyuchMXD8q1NMRiJkMlqEwbgIra5c8Tu18yQmoxdkTEZpeittmaW0Z5XhSXNiUrj+YKBvhD075CMeyJmMtDoNc2qL8XQEfVa5Rcqdvpu7DjMa0WT09odMRrmpiSHFkYOWiuLYyk7PdCwmPdWl8rLApa9tpn/Ih9mop2leMa01Dpqr7CTHaP+BSIHrfL4xfnrjkwCUOrM+cCrHKK5PcrKZ5Agbrdx156sEAhJJSSY8XhdNTS7q6uxYYiSRq3DLS3yP9gyybuNeNO9v/FMf9JUU5invcCeTkZ1ERrb8au5nHw4OJNKzkvB2lNLYXka1x4E+Bgv+NBoNc+qKZfN2bD7Azi0H0CVoqfI48HSU09hRRtYMfDSTKUxPpjBdvvO++dlghNei9OTgauQKBzXFuTFxKguNiQRDo2yeb+TZYGcgzOgN7cFOwLAQjVa5DX0yTpu8c3osMMG1G4Kdfoktm/asoFO5IikXTQycyrYkU8Tf+ZE7g9FOrYkmGtpK8HSUUd/qxhqjLVorCuXvuX94lDfe24MQMK8o54RT2ZWjTKjxsegQTsZr73TyrZDs1BTXaKfv0VhVxKUX1OOKs+z0vfcO8p3vxl92unT5Rj551hwu+3ScZafdA/QeG+Sa//0cLWfHN9rpskfXccnl7Vz0xVZS4ig73bj3EDaTnlu+fDHNJUVxKxdgeHQpBusPsVr/M66y01cOb2JeSh6/qrmEiqTZ3b95MoFAgA1rOvn6jy/g3IvrYjLIiJYnVm/m3NoSvrCgDrvCtQeT0V533XUzr1WcuOWWW6678sorY3rN4VE/v7l7OQ8vX8/2vd34xyfITJ39ndEKclO5+W/LeODxtaxat4ue3iGsZgMpyeZZteVqtRr27uvhr39bwdKl6+nqOkYgIJGRYZv1ziE3K4lf/OlZHlr6Nm9v7KJ/cJQkm2nWd4MzWQysW7mdO29+kRefeIeDXcfQaDWkZyfNegiBpFQLN37vfp68eyXvvdPFyNAoyek2LAo344mWjEQL961azy0rVvPshm0c7B3AkKAjK9E6636KCWmC7t4f0Dd4B/6xd5GkUXTabDTTiFekhEJLOrd0PsX9XS/z+tEtHPP1Y9YaSNPbZvWdEkKlAqegAAASDElEQVRwrHuAv/76aZb+czV7dhxhYiJAelbirIciKcxIZvH9L3Lvy2+zalsXvYMjJJoNJFtMJ+558eLFB6+77rpborqXj4MPYXjUT+fJQlf8fkroCo2gqiSPtmoHrdUOiiJsThENu7qOMhLBh/DUC+/y1PPhCpGsjESaGhy0NLionluguGPq6xvhQAR7+sDAKNf8+KEw2WlCgpaamqLgLl5eF5kKJZgAW08WuuKfr3/Ih1CYm3rCfDRnBrvEHTvST3eE0BX7dh/lxqvCQ1eYzHrqWtw0zi+jcX4pyalKQ1cE2LZhX8T8P/70EXZu3h+W5qzIw3NWOZ6FFbgq8xTPWvb19NEzKL/h/OrOffzu2fDQFUlmI+0h2WlLSTE2paErpBH8/i2Rcjl87AomwnwIGgz6OizGRZhN55CgK1HcSO8ZOsJQhNAVLx5+m4f2ht9zuiGRprRyWjIqqEtxzyh0xb7dR2XzRkf8/PTrd4X5EHS6YOiKxvmleOfPbO3J1v3dEWWnd69Yx/NTfAgF6UnMrwyaj7ylRarsdDKbOw/xxcXKZacFWcm01Thpq3ZQ5c6dltPzaz+4l43vKZOdmowJ1FcX09zgpKneQeo0Yu289NIWrv/Z44rKBXA6M090DqWlOdNaLHfhFcplpzarEW9NMPy3p6Z4WiFE/nnbK9z2W2WyUyEEZVUFIVlgKcXu6DecH/OP8+mKqxSVC5CamRh0MJ9VQXWLG6MpetPl9Y8u44FVymSnOo2GOnseHeXBsBaFadELCfxjO9h3uE1RuQA6bSFm0yLMxkWYDE0E42ZGx/ffWcKqY1s/+kAZDJoE6lJdNKdX0JxeTroh+gjCa1/bxo+/qlx2WuTKxDM/KDstm+bakwtvmIHs9GZVdhpTfGPj+PzjjPrHmAjErwMdHw/g843h84/hjzA6mC18vvHQvzHiOWgYH5/A5x/H5x+LuLhtNpAkCf/oGL4RP/7R+N6z3zeGb3QMn2+M8bH43fOEFAg+26F/8bxnSRpFkkZC/8fv2R6XJvBNjOGbGMMfiPM7NRL6nUfHCATkZ9Cnmo/FDGFwxMeWXfJysUBA4po/P/WhDXIq7FknZgXuwgzF09utOw4xOBRhg5yXNvPMsnDZaUqSmaZ6B82NTuqrizFPY7Q4mZ6eIXbvibBBzpCf6xY/GiY71WgEc+fm0+R10dTkpqBA+fR2/ZZ9EU1G/3h89Yc2yMlKt52QndZUKjeTHT5wnIMRNsg5uK+Hm66dskGOXke1x3FiVpChUGobCATYsKozYv6tv3iSzimy0wJnJo0LyvGcVUFFbRFahVLbXd09HI64Qc5B/vDC62FpZn0CLSVFwZXKpXbSrMrs+oHAML6TbJBzpOebU0xGoE+Yg9l4DhbTIvQJ8xRvkLN9YH9E2ekr3e/y6L7we05KMONNK6M5vYKGtBKsOuUb5HRG2CDH5xvjhu/cF75SWSMoqyrEO78MT0cZhQ7l7ciG3QcZ8ctHTX7wtfW8OGWDnIwkS3ClcqWDjrlOVXY6GavJQENFoWzes69voX/Ih1Gvo3FOUXCXtCoH6THaBL3UJb8pt88/zg2/CwaQdRSl09LoornBSXnJ9MwzkUhNtZAawSZ+9z0rCQQkrFYDjY1OmrxOGhud2GLk6Kwqlw8s1tM7xDub9iEElLtyTkhtnUWxifCalZsScdez5T9ZD0BKmpWG+aV4O8qobXJhNM9cWabRaKhucsnmdW45QOeWA2i0mmAso7Mq8JxVQW6MttC0Z6Riz5DvvG9ZEdx6JDc58YRZqMGRh14XiwivZkzGVtm8oZHnmQgcRmDAaGwN+g2Mi9DpcmdcLoDbJi8hHw9McON7QdlpkTmT5owKWtIrqEwqikkso8RkM9Uep2zeY/e+jt83jtlioLbZjbejjIY25SHVpzKvWF7K3T88yqptewGoKMg84Tcoz1cmX/9YdAgn49V1O/n2JW1ctLAKSzxlp69uob66iEs+XYcrQrjo2cDvH2fzpv189zvncc45c2IS/jlannphA5/smMOlF9aTlx2bkNfR0HN0kGM9Q1z9+8/RujC+stPn73+Df728nYuumE9qhNDcs8G7+w9hMiVwy+UX0+wqjOsq5F0DL2GxXE9+4mUkaOMnO33n+Dpa0lP4RM4lOG32uJUbCAR4e8Umvva9czn3Mi8GQ/zakSdWbWZRtZsvdNTiiFe0UyHEecBNBDfIWSJJ0q+m5BuAu4A64BhwqSRJu0N5VwGXAxPAtyRJei6aa8oxG9FOn1i2gV/97QWSbSaaax001znwzCtSvJdxtIyM+vnXr9xK/+DorCzSOhm33/Yy997zOpmZiXi9wWinNbXFsy6RO9Ldz+evuBWAmnmFNHucNDU6yZ7G1qBK+eXix1j+wibyC1LxNLtpanFTOS9f8aroaNm6rpNvn/0LjBYDtQsq8J5XRcOiuaTM8j1LksTn/v4Ab+87SFlWBgvcDs4qcTAnN2vWZafv9q7kwb03YtJacVtrKU2sx22rxaSNzWg5EhPSGPd2foHB8SNkGSsotjZRbG0mVV886+/UA795nCU/uofUnBQ859fivaCO2rPnYZzlduT44AifWHwbY+MT1LnyTswQCiYtVoxptFMhhBbYBiwC9hHcUvNzkiRtnnTM14B5kiR9RQhxGfAZSZIuFUJUAPcBjUAu8CJQEjrtpNeUQ2mHMDI6RtdBedvyxESA/77hEQaGPpCx6bTBvYxb6xy01DnJncHL23Wgh9EIO6Y99tx6nnghXCGSn5NMc52TloaZhXEYGBjhUAQJ5kD/KD/64QNhslOjMYHauuJQ6AonqalWReUC7NzdHXbtydx+96u8sSbc3u4ozqC50Umzx0WZO1vx2oDjPYMcOypvT+/ac4xfLg73IVitRhq8TrzNLhq8zmntZTyZQCBA58bIstPfffPv7NzwQWgLIQSldXY8587Dc14V9sp8xQ3Wwb4Bjg/L29Pf2N3Fb158NSwt3WJmvtvOWW4HzY4izNPYH3wyYwE/3T75e5aQuHf3DQyMf6CK0aClyFJOqa2essQG0gzKzUd9/v34A/JS2239z7P+eLi82JaQTbGliSJrE3mmKrQaZSP4of5hDu6U90WODI7yo3N/hn/Su643JlCzcC6eT9bh/VQdGfnKR/A7Dx3DH0FkcceLa3huyo5pjqzUE6Er6lz5Me0QmoDrJEk6N/T5KgBJkn456ZjnQse8IYTQAYeADOBHk499/7jQaSe9phyKZac7DnHF1fdO+7z3sRek0VrnpLXOQYU7Z1o7Tn31KuXRTq1mA56aYprrnXhr7dNawPXSiplFOy0ty6GpyYW3yYXLFb0EE+DiL/yJYz3KZKfJSWa8DQ5aGl3U1xRhnsYI64F732DJX5YrKlejFcyZW4C3xY23xU1BYfQv75h/nAuyvqKoXIDM/FQaz63Ce14V81pL0U9jj4Zrly7j/nXKZKd6rRZPcQELShwscNvJTYrepNXt28/N276hqFyAdEMepbZ6Sm31FFrK0YroBz5P7fshXUPKtmZPECYKLPUUW5sptHgw66JfMb/m2be5+vxfKCoXwFVjx/upYOfgrnNMy3R5OslO84C9kz7vC6XJHiMFNWR9QNpJzo3mmgAIIa4UQqwVQqzt7pZXzcw2ff0j9PYP09M3zFgcZYEjvjGO9w3T2z8cUak0W/T1Dp/4F2m0PxsMj/jpC91zpKCAs0FgQqKvd5je48P0HR8OU2DNNoN9w/QdG6D36ACjw/H7nf0TE/QMD3N8aJje4dG4yk6HxwcYGu9jaKKf8UD8fudxycfIRB8jE734AgNxKxegr7uf3iN99Hb3x1VePB2imSF8FjhPkqQrQp+/AHgkSfrGpGM2ho7ZF/q8E/AQnA28KUnSPaH024BnQqed9JpyKJadDvvYuE1+lD4RkLj+D8+EmYwA3EUZtNQ5aa13UObIVqz82bj1QMTGfNnK93hmxaawNJvVSFNtcGFW4zQXZk3m2LFBdkaa3g77+fnPHg9r9ISA8oq84GK0ZjfFxcqVP2+t3xOx43zo8XWsnrJSOTXFQlOjk5ZGJ3XVyvd2Prj/OPv2ypsGDx/q46YbnwlL02o1zKsuPDErULopUCAQ4K3lka2dd1z/MDvf3RuWll2Ujue8KjznVTG3uYQEhf6b7UeOcjDCBjkb9h/kD6+ER/406nQ02QtZUGKnw+0gy6bMNOgPjLJnSP6eJSQe2/enMJMRQKahgNLEBspsDeSb3WimMSuYzJHRrYxOyJtDdw2uZFPvE2Fpeo2FAksDxZZmiqyNGLXKTMC93X1sXycvL/aPjvGLz/8+zGQEUNrgxPuperwX1OGsUu7LeGvnPkb88usmHnr9XZZvCJedplhNtFXYmT/HwaLqkv8/TUZCiG4g8tZUJycdkF93fvpzptb9TK03qHU/Vah1jz1FkiRFFWo2mmHJGsAthLAD+4HLgM9POeYJ4D+AN4DPAsslSZKEEE8A/xBC/JagU9kNrCa4I8ZHXfNDRHtTcggh1kbbS55unKl1P1PrDWrdTxVq3U8tH9khSJI0LoT4BvAcQYno7ZIkbRJCXA+slSTpCeA24G4hxA6gh2ADT+i4B4HNwDjwdUmSJgDkrhn721NRUVFRiZYzKnTFTDiTe+8zte5nar1BrfupQq37qeXjFNwuqnjgpylnat3P1HqDWvdThVr3U8jHZoagoqKionJyPk4zBBUVFRWVk3BGdghCiPOEEFuFEDuEED+SyTcIIR4I5a8SQhRPyrsqlL5VCHFutNc8XesuhCgQQqwQQmwWQmwSQnz7TKn7pDytEOJtIcRTZ1LdhRDJQoiHhBDvCSG2hCTaZ0rd/zv0vGwUQtwnhIj5np5K6y2ESAs904NCiD9OOadOCPFu6JybhZidIEWxrrsQwiyEeDr0rGwSQnxk7LZTgiRJZ9Q/gqqknYAD0APrgYopx3wN+Gvo78uAB0J/V4SONwD20HW00VzzNK57DlAbOsZGMEbUGVH3Sed9B/gH8NSZ8syE8u4Ergj9rQeSz4S6E4wKsAswhY57EPjiaVRvC9AKfAX445RzVgNegtL1Z4BPnGbfuWzdATOwYNKz8ups1H2m/87EGUIjsEOSpE5JkvzA/cCFU465kODLCvAQsDA0krgQuF+SJJ8kSbuAHaHrRXPN07LukiQdlCTpLQBJkgaALUQIA3K61R1ACJEPfBJYMgt1nrW6CyGSgHaCkmskSfJLktR7JtQ9dJwOMIngQlIzoCzg1izUW5KkIUmSXgPCwgcIIXKAREmS3pSCLetdwEUxrves1F2SpGFJklaE/vYDbwHyG4ecQs7EDuGUxlaaIbNR9xOEpq01wKoY1vlD9YpUPsrq/nvgB8BsBkyajbrbgW7gjpC5a4kQYjbiO8e87pIk7QduBLqAg0CfJEnPn0b1Ptk1J4dYPR3f049ECJEMXAAsm3FNY8yZ2CGoyCCEsAIPA/8lSVL/qa5PNAghPgUckSRp3amuiwJ0QC3wF0mSaoAhQqFaTneEECkER7h2ghEELEKIfz+1tfp4EJqR3QfcLElS5H1XTxFnYoewHyiY9Dk/lCZ7TOgHSCK4cU+kc6O5ZiyYjbojhEgg2BncK0nSI7NQ77B6TS1f7pgo694CfFoIsZvgtPwsIcQ9Z0jd9wH7JEl6fzb2EMEOItbMRt3PBnZJktQtSdIY8AjQfBrV+2TXnGxmOR3f04/iFmC7JEm/j0E9Y8+pdmJM9x/BkVknwdHN+w6fyinHfJ1wh8+Dob8rCXeydRJ0IH3kNU/juguCttTfn2nf+5RzO5g9p/Ks1J2gY7A09Pd1wP+cCXUnGIl4E0HfgSBoC//m6VLvSflf5KOdyuefTt/5R9T95wQHbprZeM5jcu+nugIKf7DzCappdgLXhNKuBz4d+tsI/JOgE2014Jh07jWh87Yyycsvd80zoe4EFQ0SsAF4J/Qv5i/JbH3vk/I7mKUOYRafmWpgbei7fwxIOYPqvhh4D9gI3A0YTrN67yYYF22Q4GysIpReH6rzTuCPhBbXnu51JzjLkAiKPt5/T6+Yredd6T91pbKKioqKCnBm+hBUVFRUVGYBtUNQUVFRUQHUDkFFRUVFJYTaIaioqKioAGqHoKKioqISQu0QVFRUVFQAtUNQUVFRUQmhdggqKioqKgD8Pyg7nqFrE90HAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "grad_u = project(grad(voltage), VectorFunctionSpace(mesh, \"Lagrange\", 1))\n",
    "plot(grad_u,)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Heat Transfer Model\n",
    "The electrical heat generation is coupled to the stationary heat equation as the source term $Q$ on the right hand side, which is known as the electromagnetic heat source.\n",
    "\n",
    "$$\\nabla \\cdot(-k\\nabla T)=0$$\n",
    "\n",
    "Joule's first law states that the heat generated within an object is equivalent to the product of its conductivity $G$ and the square of the potential gradient:\n",
    "\n",
    "$$Q=G|\\nabla V|^2$$\n",
    "\n",
    "\n",
    "**Cambiar el valor de Q por uno numerico obtenido de las mediciones y también puedo tratar de calcular el valor de Q del Vector**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 45,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.quiver.Quiver at 0x7fc67e34d1d0>"
      ]
     },
     "execution_count": 45,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYQAAABJCAYAAAA9vFENAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJztnXmUG9Wd7z9Xu1pS7/uurVe3e29J3W27bfaEAMkkgeQN5/EGkskCk8kkk0AIwTBAZrINCWTgAYYQdsJqFrPZ7BiMjTE2Nl7a+97tdu+LpFa9PyScViMZdbUk24/6nNPntO6tqvsrqapu3d/ve39XSJKEgoKCgoKC6kQboKCgoKBwcqB0CAoKCgoKgNIhKCgoKCiEUDoEBQUFBQVA6RAUFBQUFEIoHYKCgoKCAqB0CAoKCgoKIZQOQUFBQUEBUDoEBQUFBYUQmhNtwEzIzs6WysvLZ7zfZEBidNwbuVKCA70DBALhM7Z1WjXmFD3mFD1Ggw4hw16A0THvZ479KUMj4wwMjYeVqVQCk1GH2RRsW62W12f7/ZOMj/si1kkBiYMHBpg+S11v0GI26zGZDRgMWlntAoyMThBtAvzRwVFGRsN/C7VaYAp91+YUPSqVvG/b5/UzMeGPWOf3T3L40GBYmQAMRi0mswGz2YBOL/92GB6ZiFrXe2SYiYnw30KjUWMy6TGb9JhSdAgh75zHfX58/smIdWM+Hz1DI2FlAoFJr8Vi0GM26NFp1LLahQBIw9GrJ/eDNO23EDoQFlBZQJhA5l016vXhDwQi1g1PTNA3MhZWphICs16HxaDHYtCjUcm8pyYD0e8pSWJ/z+Bn7ymt5ti9bJzNPTXmJRDlphoYGmNoNPz6U6sEJqMeS4qerZs/7pUkKSeWdsSplLqipaVFWr169Yz329h9kH/61QOy2001GXDXl9PZaKOzyU6KQRfzvt+/6kHWb94vq10hoMZZQEeLnQVuJ2XFWTHv+/ryjdxw7ZOy2gXIzDLhanfibnfQ6nag1cb+4DjvO7dx5OjI528YAbVaRUN1cfCcXU7yc1Jj3vfR+9/hrj8vl9UuQH5hOq4OJ57OChparDF3TD7fJKdf8HvZ7ep0GpobymhvszPP4yQj3RTzvtc/tZxH3vtIdtv23Ey6qm0sqrZTX1oQc8ck+bcj9Z4tu12EGfTzEPqFoD8NobLEvOt37n+SN7ftlN30nMI8FlbYOL3aQWVedsz7vbt2Bz+56QnZ7aanGvE02uhosdHeZEOvi/0F5Js/v4ddB47KanfVfT9ZI0lSSyzbKi6jGCjKS6OsIJOywkyMevm9/ExJsxgpK8qktCiTnKzYb5jZolILikqyKCnJoqQsG40meZdJdoaJ0qIMSgszyEhLSVq7Wq2a4tLg+RaVZMoepcghL8dCaVEmpcWZWCzGpLVr1Gmx5mRizc6kODNN9ihFFuoSUNtAbQ2NFpJDmkGPNSsDa3YGBanmpLUrBJQUZFBWlEFZYSa6GbxgJZMvxAhh3OvjQM9gxDqfb5If3vho2JBLr9PQNqeMziYb7Q02cjLkXzgHDg8w4Y3sxnjsuQ946sV1YWW20mw6Wux0tNqpduTLdhmNjEzQ2zMUsW6gf5R//9H9BCb//tubzHpaXXbcHU5a3XZSU+U/mPYcOMrkZOQh/f/c/wZvr+4OK6tx5AfPucWOoyxH9oNpcGCU/igjk53dPdzwy8fDytIzTLg6nLg6nDS32TCmxD7ym4okSeza0xe1/sbfPcuW7kPHPqtVgrraYtpdDtrb7JQUZcpqF6BncJjB8cjuqre37OK/nns9rKwg3UJXlY2uahtttmJ0GnluMknywuTuKJUBpKP/BIHDUwp1oPeERgQLEeoCWe0C7O8fZMwX2XXz+NqPufudNWFltuxMFlZYWVhpo6G4EI3Me2p0zMvhI5HvqZExL5df+yhe39/vdaNBi6u+nI4WO55G66xecPYe6o/qGlzy9Lu8/O7msLLK8lzmNdjobLRTY8uPeYRwSsUQ5GLQabEWRXa3PPHKOoZGJ8jJMB9zCTXXlmDQxWckUJCbFrF8dMzLq+9sQatR0zinJDiMbLFH3X6mmEx6TCZ9xLq7bltBYFKisDgDT4cTd4eTOXNL0Mj2J4dTUpARsfxgzwDvrt2BQa+hrb6cjmY7niYbWRnxeUNMTUshNcpNd99dbwBgc+Ti6qzA3VlBZXVhXEYCQgjKSyNfXxs27mNL9yEsZgOuFiueVjuuZhsWi2HW7QLkpJrJifCmK0kSv/jbiwgBc4sL6KoOdgLOvKy4jASE0IHGEbFOGnsm2BmoskHfhdAvAl07QhWfEV9hemQ3otfv57n1m9GoVLSUFbGwwkZXhY2yrPS4tJti1FEexW177xPv4fX5KchJDb7cNNtoqClGp43PI7Y4L/I59A2O8sYH3ei1alpry+hstNFRbyU3U55H4QvRIRyP1Zt2c/m353PBwrmYUyI/QBPBayu34G628o0vN1Nlz0tau16vn+07evjRz7/EmWfVoZuBH3O2vPj6Js5fVMc3v9JMSaH8t+KZ0ndkmOGRCX5xw9eYt6AadRJdYK++volvf62Nfzi/mezs5Ln9Pt5/mKxUE7ed8VU6HGVJdYEdGF2JMPyMLPM30WlijwHNluU7ttFYk8/FjU005RUmzQUWCATYvH4vP/3HLs4+sw6jUd5IUw4vvruJs9qruOisJuxFscdDovGFcBkdj5ff28zVtz2HVqOmqao4OMyqt1GYE5839Wh4fX4uumIJh3qHKCnIoL3ZRkeLnfqqori9qUfjwYdWctddr2MwaGluLqfd48DlspOZmVif6tGBUS785zsYG/dhL8+ho9VOe6udakdBwh9YN//meZ5buhaLxUCr246r3Umry4ZlFq6xWNixo4fLvrMEgOqqQjweBx6PA6tVvmssVi679wne3raLbHMKCyptLKy04bGXkhKn0W809o+8w2sHfoxATa6xkSJTJ0UpnVh0JQltdzIQ4Mxn72T7YB8l5nQWFdk5rdiJK7cUnTqx99TSJ1Zzy+9fQKfT0NRqxdXuwN3hJHsGogg5DI9NcP7PlzAwMo6tMIt59Tbm1duosxegDqmphBAxu4y+EB3C0Mg4H22LrPQJBCSuu+sFBqfJBu3FWcxrsNPZYKPWln/sy50pGzbvZ3BkPGLd8rc/4YXXN4aVWUx6XA1WOlrsuBvKSZUZZDxyZJit2w5FrBsb9XLTr5/5jJ+/qqoAj8dBu8eBzZYr+4G1et0ufL7I/s7HnlvDqrU7w8oy0lLwtNjoaLXTUl9Oisw3rP37jrJn95GIdYcPDvCn378QVqZSC+rmluIOuc2KS+SNWgIBiVXvb49af/fdr7Nt2+Gwsry8VNzu4HddX18qe6S25VAvBwYi+7XX7z3In199N6xMp1HjtpWysNJKV6WN/DR5oxZfYJTDY2uj1EqsOvxrxiZ7w0pTteXBzsHUSbahDpWQd87revfTNzEasW753m08sDXcLrNWx7wCG6cV2ekqcpBlkOe66j86wuZNByLW+Xx+fr34KbzT4oXOyvzQ9VWBoyJf9ovPms17GJuIHDd5+s0NvPrBtrCyNLOBjjor8+vtnNFWqXQIU9m4/SCXXP+g7HYzLEY66m101lvxzLXOSGn0vavly07VKkFdVREdzXbmtTmi+uYj8dprm7j+P56W1S5Abm4qHrcDt9tOU1P5jGSnF/yf/5EtO9Vp1TTWlQb9sG0OcmfgZnnkgZXcddsKWe0ClJRmHXuzq6svnZHs9Kyzfyu7XaNRR0tLOR53cPSQNoPg4+Kly3nkffmy0+r8HLqqbCyqslFbmBfzS8CgdxfP7v6m7HZ1qlQKU9opMnVSaGpHq4o9jnTJikd4fX/0Dvh4CKApp4jTipycUeLEkRa7m+X9d7v5xU8ektUuQFa25dj11dxqm9G8l3+4+h52HZQnO11ztyI7jStZ6SZyMszkZlrQxylIFAspKXpyMi3kZJlJS6IcUQjIzjKTnW0mJ8eSVNmp2aQnO9NMdpYZizl5MR21WkVWjoWc3FSycyxJ9bmnpRnJzraQnW1Jqv9Zr1GTl2YmL9VMjsWcVNmpUZNFiiYXoyYHtYhPgD0WzFo9eUYLeSlmMvXJkzUDZGWbyc6xkJ1jQas7hWWnQoizgT8CauAuSZL+c1q9Hvgr0AwcAS6UJGlnqO4q4FJgEvgXSZJeDJXfDZwLHJYkaU4sxsqWnU742H2oP2Kd3z/JFb97PEx2qtWoaakuYV6DjY56GwXZ8v2Aew4cZTzKUO/JF9fx9MvhstOSwgw6m4MSzLqqItkSueHhcQ4dGohYNzA4xs9+9kjYDGqDQUtrqxW324HbZSdjFsqfHbt7o8pO73zgLd6ZJjt1WnNpbw1KbSvt8ofVR4+O0Ncbefbs7l293LT4qbAyS6qRNo8dd7uTljYbZpnKn0BAYseOw1Hrf/Pb59m69e/uOyGgpqboWDyhvCxb9sP44MAQ/aORXZIrt+/mNy+8EVaWYzHRVRmUYLptpRhlxhMmAxMM+nZFrJOkAK8f+EmYy0iFhlxjE4WmDopN8zBri2S1C7B76CjD/siZB57YvoElm1aFlZWa0zmt2MFpRU5ac0tkxxNGRiY4uD/yc2RszMvPf/RAmMtIr9fQ2GLF0+mkzeMkO0e+qGDngT68UWSn9zy3ipffD5ed2ouymF9vZ169jXpnUfxkp0IINfBn4AxgL/C+EGKpJElTnd+XAkclSXIIIS4C/gu4UAhRA1wE1AKFwCtCiApJkiaBvwC3EuxIEopBr6WiNPLM7ade+4ih0QkyU1PoqLfS2WDDVVs2o9nIxyOam2ds3Mvr721BrRLMrS6ms8VOe7ON0jipb8yhdAyRWLLkdQIBKegW8szejz0da2nkYfihnkFWfbgDnVZN09wyOlrteJpt5MUp8JaRYYrakT3ywEoASsqyjklta2qL46I4UqkE9ihKsU2b9rN16yFSUnS0tFjxuIMB/PT0+Lyd5qdZIsYBJEli8TPBWds1hbksDAWVqwty4zL6Uav0ZOgrItbtGnqZscle9Ko0CkztFJvmkW90oVPHR7RQaol8T3knJ1m26xNUQtCUXcRpxU5OK3LgSIuP1NZk0mN3Rv6dH77/HbxeP9k5FlztTjydThqay9HHaSJreUHk58LRoVHe+mg7GrWKlqqSY0Hlwmx5ophYngBtwDZJkrYDCCEeBs4HpnYI5wOLQ/8/Btwqgr/A+cDDkiRNADuEENtCx1spSdIbQohyWVbHkVUbdvOjC+fz1YVz49YJxMLydzbjabJx4bnNOMtzk9au1+unu/swP/23sznjzDlok+gCe/H1jXzljHouPL+Fovz4aMNjoe/IEMP9I1yz+AI6F1WjkikQkMOrL2/g299o5WvfcJGZlbyZsRv2HiLLaOSui79Gu7M0qe6g/rFXcWddSknqt9Cqkye1fWX/ZtoKCrmkoo36LPkjkJkSCATYtG43V/zrWZx1XkPcOoFYWPbeJs7yVPG/Tm/CWhB7aptofK7LSAjxdeBsSZIuC32+GHBJknT5lG02hLbZG/rcDbgIdhLvSpJ0f6h8CbBMkqTHQp/LgWcT7TI6Hi+8vYnFty1Dr9PQWlsanNjROLvZybEw4fVz4RVL6OkbxlqSRUeznc4WGzWOAtmzk2Plgfvf5u6738Bk0tPaZsPtdtDWZptRMFMOff0jfON7dzLh9VMdmp3c3mLDaZWvaIqV/178FC8+uYb0TBOtnRW4FlTS3O7AmOC5J9u3HuT737odlVowp74U1/xKXJ0VlJTPXjP+eVz6fx/jvW17KEi3sKDaRleNjVZ7ccLjYP2jL7O991IEeiwGD2nG00kznoZOk9iHtD8wyXkr/syukT6s5iy68ivpyqugIbMEjSqxPvunH3qX2/7reYwpOpo8DlzzK2mb5yQ9wVLuwdFxvnzt3QyNTVBZnMOCOTbm19moKck7NhKMq+z0RHcIQojvAt8FKC0tbd61K7Lf8niMjE6wcXtkCWZACnDNrc99RnZaZc0LzlxutFFZLv+Btan7ICOjkVMLvPjGJp5/7eOwsnSLEXeTlc5mO231ZZhkPrD6+kbYubMnYt3oyATXX/9UmJ9fpRLU1hbj8ThwexyUlsofZq/buDfqNPuHnn6f96bJTnOzLLS3BOdhNM0pkf2GdWh/P/v3RJadHtx3lD9eF6660mrVzG21Bh/SC6rIK5Q3agkEAqxbvTNq/Z03v0T3loNhZUWlmbjmVeKeV0FtQ6nsuSfbD/dxeCBy3OTDnfu59aWVYWVGnRaPs5SuGhvzq61kW+TFiiYDo4x4P4hSK7HryI/xTYbHVYzaGtKMp5FmPJ0UXT1CyHvx2TRwgAHvWMS6Vw5s4qEd74eVpWoNzMtz0pVXQWeeg1StPIHGYP8o3Zsjy069E35u/PdH8E7JtiuEoKquGNf8ClwLqih3yH+OrNu+n3Ff5BQ4j76xjuXrwmWn2akpzKsNdg6L6h1x7RA8wGJJks4Kfb4KQJKkX0/Z5sXQNiuFEBrgIJADXDl126nbhT6Xk4QRwsbug/zTtfJlpzkZJjoagmktWmaY1uJ7v3xItuxUo1bRWFtCR7ONzhmmtXjttU38x/VPff6GUSgsTMfjceL2OJg7w7QW518qP9upQa+hZW5ZaPRgn1Fai7/d8yZLbn5JVrsA5c483AsqaZtfSVVdccyuJZ/Pz7meG2S3a7YYaPE4cM2voLXdOaPJctc/vpxH35UvO60ryaerxsaCGhuVBbEHt8d93Ww8sFB2uxpVDmnGRaQZT8NimI96Bmkt/nnl/bx1eNvnbxgBtRA0ZZXRlVfBwvxKysyxu1lWv72VX/7wPlntAuQWpIeurwrqW60zctde8B9/YechebLTdX/+t7jKTt8HnEIIqxBCRzBIvHTaNkuB/x36/+vACinY0ywFLhJC6IUQVsAJrOIUI8WoD+brN+rQJHjG41S0WjUmow5Tih5DEv2SACaT4Vg+pGT63PU6LaaUYLszSQ88W1QqgclsOPaXTJ+7waglJbQORTJTiWjUKswG3bG/ZJ6zWmVCrbKgVllQieRd2zqVBotGj0VrIEWTvJghgMmiJ8UU/J0T7RaWS6yy0y8BNxOUnd4tSdKNQojrgdWSJC0VQhiA+4BGoA+4aEoQ+mrgnwA/8K+SJC0LlT8EdAHZwCHgWkmSlhzPDrkjhNFxLzv2RXYl+CcD/PvvnwpzGanVKhoqi465jEryY58QNp0de44wFmVxnqXL1/PM8vVhZfFKjjUwMMqBKBK5oeFxfnHVo2GyU51OQ2NTWXBU4HaQMwuJ3Jbth/BHkZ3+5dGVvLMmfFJReUnWsXjCnIpC2TfLkZ4heg9Gltru2dnL76ZlO00x6WnucOCaX0Vrp5M0mVLbQCDAlo3RR4G3/PpZtm0OdxlV1BTi6qzANb8CR2Xs6xBMZ1/fAH3Dkd0n73fv4Q/PvxVWlmEyMr/KSleNDU9FKWaDPJdkIDDOmO+TaLVs7/3uNJeRCrO+JRRPOB29xi77nHcM9TLkjyy1fW7veu7f/l5YWb4hlQX5wRFBW3Y5erW8Dmh4cIx9uyI/R8bHvVxz+f1hLqMwl+T8StkuSYAt+3qYiDL7/74Va3jpgy1hZaU56SyoszF/jo22ytL4ZjuVJOl54PlpZb+a8v848I0o+94I3Bih/FuxtB0PUgw6au2R0+0++/oGBkcmSDXp8dRb6Wyy464rw2KKz2QZa0nkIen4hI+3VncjBNSGFsHpbLZjLYmPRC4tLSVqkPiee94gEJDIyjIH5x14HDQ1lc9qlbSpVNgiS/N6jgyx+qNdwQ63tvhYyut4KY6ycixkRenInn4omMKhoDgD14Iq2uZXUtdcFheVlUqlompOccS6LRv3sW3zwaAmvc0WCio7yYqT1LYoM42izM+6EiVJ4rfPBlNfO/KyWFATDCrPLZWfhmUqKpUBk74hYt3R0efxTR5GJSykGhcEOwHDQjRq+S9WU7FaIgfkvQE/P139GABz0guDQeX8CqpS8+NyT5lTjVTWRf6dH7v3LbwTftIyTLTNC4oWmtx2UqJkHJ4pFUWRZfP9w2O8vXEnapWgwVYU6gSslOfJk69/4bOdvrN+J1d8ez5fXZRk2enbm3E3lHPhl5txWpMrO9265SD/9pNzOPPMuhmlpJgty177mC+fXsdFX2mmaBajrpnS1zvEwNERfvHbC+k8vSapLrDlz67jwn9s54J/9JCZzGynew6RZjRwx2VfxVNRllR30NjYi5Sl/zsZ5ktQzWAltNny2sFPaMwq4g+tX2dOeuQHdyIIBAJs+GAXP7jyS5x1QTP6OL1YxcKzqzdxVnMFFy9qwpqXBNnpyUQiZKfPvbWR6+58gRSDFveccjobgvnEM1ITK8GcmPDxjcvv4sjRESqsuSE3kZ1KW17C0yb89a9vce9f3iQ11YjLZcftcdDSYo06kS1e9B4d5uuXL8Hvn2RORQEdzUHXmLU4PqOi4/H7ax7n5afXkpVroW1+Je4FVTS47Am/ebdt3McVX70FjVZNXZsN18IqXF3V5MtMpDcTLrntb6zevpeSrDS6QiOEJmsR2gTHwcbHltHfdylCGNHp56E3nIHecDpqdWLTvPsDk/zDmzezd7QPhyWP+blVzM+poia9GLVMRVOsPPnXt/m///U85lQDzR0VuLuqaJlXgSUtsSlnBkfHOef6JYxMeJlTms+CWhvza21UFP5dIKBkO53GyJiXLbsjpxaYDAS48pZnwmIIQsAce0Eo26kVe7H81AJbdx7+zKLyn7LstQ08u2JDWFlWhon2pqAEs3Vumexgcn//KLujZP4cGRnn2l89ESY7VatV1NeX4vY48LgdFBbJf4P/eOuBqLLT+55axcq1O8LKCnPTgkqqZjsNNcVoZUowew72c2hf5LjJ/j1H+MOvwteY1hu0NLjsuBYEfbxZufLcOIFAgI0fRJdD337jM3RPizGUOfNwLazGtbCayvoS2XGT3b399AxGlp2u2bGfP73wdliZ2aCjo7Kcrmob86rKSTfJe2BJgVF8vvXRaunv+z6BQLjUW6ttQG84Hb3hDDTaObLvqW1DhxjyRY6bvHxwA4/sCpfaZuhMdOZUMj+3Cle2A5NGnhtnqH+UXd2RnyMT4z6uuyI8hqBSq6htKgteX11VlFhjWuc+Ih/vPshEFNnpg298yMvrtoaVFWRYmF9rY0Gtjc4aq9IhTGXj9oNccp182WlBdiqdDTbmNdhoqppZoHc22U51Og0tc0qPrZWQO4N1lWcrOy0tyzoWYK6tLZrRA+u8f76dXpmy0xSjDld9OZ3NNjyNVtJnMFL7291vsuTmF2W1C+CsKcS1oArXgkoc1bEvsOLz+jmv7pey203NMNG6oBLXwmqaOp2YZjBSm43sVCUEDeXB1dQW1Niw5WbGfM5+3zZ6D8+X1S6ASl1wbOSg13cgROwd0xXv/4V3erd+/oYR0Ao1zVnW4Oght4oCY+wvPqvf3MIvv3evrHYBisqycHVV4VpQRW1TGZoZuGvPv+kv7DwsT3b60R/jKzv9wqMSArVKoFKpkuqLFYBQCdRJbhdArVKhUgXPO5mohEAV+q6TmXEUgsFhlVokNcYQbDd0zmqBKpnXlwCVUKFSqU7ANaZCoEKgJnilJwchBGqhQoUKkcR2IThiEKHfOcm3c8x8IUYII2Netu6JPGt3cjLAz29ZGuYyUglBnaOAzkYb8xrsWAtjf3Oaztadhxkdi+wyev7Vz7qMcjLNtIfSWDTPKZU9a7e/f5Q9UWbtDg+P86trHg+TnWo0QZfRp5PRCgrkK382bovuMvrrk591GRXnpx+LJzTMYsW4noMDHNof+S1q/54+/nDNE2FleoM2mGYgNBlNbtA36DKKsuA8cPuNSz/jMrJW5tP2qctobuyT4Kazu7ef3qHIo7E1O/bxx2XhLqNUo57OynK6amx0VJaTliIvbnR8lxH0931vmstIoNU2ojeeEXQZaapl31PdQ4cYnIHLKEtnpjM35DLKcmCUOf9gaGCMXVEWnYrkMlJrVMxpKsfVVYW7q4rCMvlB3417DkWdqfzg62s/4zIqzEhl/pygy6ijujy+stNTHZNRR0NF5Dwqz721kcGRCUwGHe66cuY12mifayU9TusPREtcNzHh45rfPwNApS2PzpAEsyJOeX3S01OiZtT8671vEghIpKUZcYUWZWlutmKKk0SuxhFZ4tt7dJg1G/YEO9yqwmCa72YbZbPocKeSk59GTn7k2dwvPLEGgOy8VNwLqmhbUElDmw1dHCb8qVQq5rSUR6zbtnEf3Rv3o9GqqXfZcS2qpq2rirxZxGimUpqdTml25M77Ty+8A0BZdvqx2ciN5YVxCSoLVQo6vSti3fjYMgKBQwiRgk4/f0pQWb4PfSp2S+TgtC8wybUfBeeaVFjymRdyC9WkFaGKQ1DZkmZkTnN5xLon7n0b74Qfc6qR1nkVuLqqaOl0Yo7TEq01JZHPeXB0nJWbdyEEzC0rOBZUdhTIE2p8ITqE4/H2uu38y0Xz+YeFczEmUXb6yluf4JpbxoXnteBIcrbTTzYfOCGy0+de3cC5i2r51rktFOUlMdtpzxD9vcNc/buL6Ei27PTJNXzzOwu44JJOMpIoO92w+yAWo447vvM12ivKktYuwPDYMgyWK0k1X4JKldg1hafyxuGPqc8o5D8bv0lNWmLXb55KIBDgo9U7+OEvv8JZX2uOy0tGrCxdtZGzGiq4eGEzVplzD6aiXrx48eytShJ33HHH4u9+97txPebouJff3reCx1esY+ueHrz+yeDKaAlOIVBamMktty/n0adW896aHRw9OoLZpCcjPSWhvly1WsWevX3cfserPL9sHbt3HyEgSeTkWBLeORTlpXHj/7zI35at5YOP9zA4PE6axZjw1eCMJj1r3t7KvX96hVeWfsiB3UdQqVVk56clPIVAWqaJ3/30YZ65720++XA3YyPjpGdbMMlcjCdWclJNPPTuOu54dRUvfLSFA/1D6DUa8tLMCY9T+KUA+45eSd/wvYx51yNJ42jUBahUif2dS1OyWbLjGR7d8yorezfS5x0kRWMgU2dJ6D0lhOAh/f/NAAASFElEQVRIzxC3/3YZzz/2Pru2HWZyMkB2XmrCU5GU5qRz3UOv8MBra3lvy276R8ZITdGTbjIeO+frrrvuwOLFi++I6Vy+CDGE0XEv24+TuuKnN09LXaES1FcUMa/BRmeDjbIoi1PEwo5dvVFTVzz70nqeeylcIZKXk4qn1UZ7m4OGuhLZHdPA4Bj790X2pw8Nj3P1NY+FyU61WjWNDWXBVbxcDnJlSjABNh8ndcXdj638TAyh9NNV4prt1FUWyl4l7sjhQXqipK7Yu7OX3131WFiZMUVHc4eTtpALKT1zFqkrPtobtf7WXz1B98Z9YWX2miJci6pxnVaDo7ZI9qhlb98AfcORF5xftX0v/70sPHVFWoqB+SHZaUdFORaj3NQVY4z7NkWrZW/vd/GHxRBUGHXNWIynYzaegV5TIfshvXvkEMNRUlesOLSWx/e+GVaWrUvFnV1De3YtTRnOWaWu2LuzN2Ld+JiPX11xX1gMQaNRM7el/Ni8l/xi+W7Czft6ospO71uxhpc+DI8hlGSnBVNhz7HhrixTZKdTma3stCQvnXmNduY12Kh3Fs4o6PnDnz3Ahk3yZKdGg5aWxnI8rXY8LTYyZ5Br57XXN3H9DU9//oZRsNtyj3UOlZUFM1L8zEZ2ajEZ8DSW09Fsx91QPqMUIn9b8gZL/iBPdiqEoKq+JCQLrKTcGfuC8z6vn/NqrpLVLkBmbiptXVW4FtXQ0OHEMIN1la9/cjmPyJSdalQqmq1FdIXWSijNit2NN+HbRvfBBbLaBdCqS7EYz8BsPB2T3k0wb2Zs/PzDO1nVFy2P0vHRq7Q0ZTjxZNfgya4hWx97BuHVb2/llz+Qv8BjmT332JyXqrkzm3ty/g2zkJ3eoshO48qEz8+E18+418dkIHkdqN8fYGLCh9frwxvl7SBRTHj9TEz4mfD6SOZLg98/GWzb68fnjzzKSASSJOEd9zEx5sU7ntxz9k74mBj3MTHhwx8lgVkimJQCwWvb72fC50/qOUvSOAFpDEkaJ7iibnLwS5NMBHxMBHx4A0m+p8Y//Z39BALJu7ZnwhdihDA8NsGmHVEWyAlIXP0/z35mgZwaa96xUYGzNEf28HbztoOMjERZIOfVjbywPFx2mpGegrvFRnubnZaGclJm8LY4lb6jx1kgZ9TL4uufDJOdqlSCujnFeNwOPB4nJcXy3WTrPtkb9WH+4NL3effDnWFledmWY4qjxhr5brJD+49yYHdfxLoDe/v447XhE/W0Og0NLtuxUUGOTKltIBDgo/e2R62/86Zn2D5tlFhizw3KThfVUNNUhlqm1HZHTx+Hoi2Qs+sAt7z0TlhZik5LR0UZXTU25lVayTLLS9ESCIwyFmWBHAmJ/Ud+NM1lBAbtHMzGM7AYz8CgrZO9QM62oX0M+iK7yd7sWc9T+6ZJbbUpuLKqac+upSWzArNG/gI526dlrf2UiQkfN/704fCZyipB1dwS3AuCCzCV2uQ/Rz7acYAxny9i3aNvruOVaQvk5KSZgjOV59joqrMrstOpmI16WmtKI9a98M4mBkcmMOg0tM0pY16DjY56G9np8vzJ06l05Ecsn/D6ufEPwQSytvJs2lsdtLfZqa6YmXsmGpkZpqgupvseeJtAQMJs1tPWasfjttPWascSp0BnfVXkxGJ9/SOs3bgXIYLS1E/TVdhL5acGmUpeYQZ5hZH9tCuuWQdARpaZ1gWVuLuqaPI4MKTMXlmmUqlo8Dgi1m3ftJ/tm/ajUquoa7XiWlSDa1ENhXFaQtOak4k1J3LnfceK4NIjhRmpQbdQtY1WWxE6TTwyvKZgMnRGrBsaewl/4BACPSZDR7ATMJyOVlM463YBHJbIEvLJwCR/2ByME5Wl5OHJrqE9u4aatPK45DJKTU+hwWWLWPfUAyvxTvhJMelpanfgXlBJa2eF7JTq05lrjSzlHhwd570tewCoKckNxg1qbVSXyJOvfyE6hOPx5gfd/Oib87jgtHpMSZSdrnhzE82NZXzzvGYcUdJFJwKv18/Gj/fxkx+fzZlnzIlL+udYefaV9Zw7v5aLzmuJW8rrWOjrHeZI3wi/+OO36FyUXNnpSw+v5BuXzuOCy7rIjFPK61hYv/cgRqOWOy77Gu2O0qTOQt4+9Bom83WUpV6IVp08qe3ao2vwZGVyTsGF2CyRH9yJIBAIsPbVDfzgJ2dy1oUu9PrkPUeWrtrIGY1OLu5qwpafpGynQoizgT8SXCDnLkmS/nNavR74K9AMHAEulCRpZ6juKuBSYBL4F0mSXozlmJFIRLbTpcs/4j/veJl0i5H2JhvtTTZcc+WvZRwrY+NevvH9OxkcHqeusoiO0JrCpUXxmaR1PO5e8joPPPAOubmpuN123G4HjU3lCZfIHe4Z5NuX3QlA49xS2l12PG128mewNKhcbrr+KVa88jHFJZm4PU7cHU7m1BXLnhUdK5vXbOdHi27AYNLTtLAW9zn1tJ45l4wEn7MkSXzrL4+wdu8BqvJyWOi0sajCxpzCvITLTjcMvMVje36LQW3GaW6iwtKKw9KEUZ3YBecnJR/3d1/MkL+HfGM15eZ2rGYPmbryhN9Tj/zmae668n4yCzJwfakJ91eaaTp9LoYEP0eODo9xzvVL8PknaXYUHZuYVjJlsmJcs50KIdTAFuAMYC/BJTW/JUnSxinb/ACYK0nS94QQFwFflSTpQiFEDfAQ0AYUAq8AFaHdjnvMSMjtEMbGfew+ENm3PDkZ4Mc3PcHQyN9lbJ+uZdzZZKOj2U7hLG7e3fv7GB+P7Pt76qV1LH05XCFSnJ9Oe2jWcn21/DQOQ0NjHIwiwRwaGufKnz8SJjs1GLQ0NZeHEtrZycyUf/N27+wJO/ZU7r7vTVa+H+5vt5Xn0N5mp93loMqZL3tuwNG+YXp7I/vT9+w+wk3Tkv2ZzQZa3cH0321uOxaZ8yECgQDb1++JWv/fl99D90d/T20hhKCy2Yrr7Hpc5zRgrS2W/cA6MDDE0dHIaRxW7tzNb16ZJsE0pbDAaWWR00a7rYyUGawPPhVfwEvvRGSprUSAB3fdwLD/76oYFWpKTTVUWFqptLSSpZfvPur37sMXiBxD+GTgZdYdDZcXW7T5lJvcWC0eioz1qFXy3uBHBkc50B05Fjk2PM6VZ/0H3in3us6gpfG0OlxfbsZ9bjM5xfLf4LsPHsEbJR3MPcvf58W14Sum2fIyj8UQmu3Fce0QPMBiSZLOCn2+CkCSpF9P2ebF0DYrhRAa4CCQA1w5ddtPtwvtdtxjRkK27HTbQS67+oEZ7/cp1uIsOpvtdDbbqHEWzGjFqe//Qn62U3OKHldjOe3NdtxN1hlN4JptttPKygI87Q7cbgcOR+wSTICvXfxnjvTJk52mp6XgbrXR0eagpbGMlBm8YT3y4EruvH2FrHZVasGcuhI87U7c7U5KSmO/eX1eP1/JkT9hMrcki7az5uI+p4G5nVXoZrBGw7XPL+fhNfJkpzq1Gld5CQsrbCx0WilMi92l1Tuxl1u3/lBWuwBZuiIqU1upsLRSklKNWsT+4rN0z5XsHpG3NLtWZaQkpQWr2UOZ2UWKJva5Ae+/sJZffOkmWe0COBqtuM8Ndg7OZtuMXJcnU7bTImDq68/eUFnEbSRJ8gMDQNZx9o3lmAAIIb4rhFgthFjd0xNZNZNoBobG6B8cpW9wFF8SZYFjEz6ODozSPzjKcBSlUqIYGBhloD/4F+1tPxGMjnkZCJ1ztKSAiSAwKTHQP0r/0VH6+0fDFFiJZrh/hIHeIfp7hxgfTd7v7J2cpG90lKMjo/SPjidVdjo2OcSIf4AR/wD+QPJ+Z39ggvHJfsYm+5mYHEpauwADPYP0Hx6gv2cwqfLimRDLCOHrwNmSJF0W+nwx4JIk6fIp22wIbbM39LkbcBEcDbwrSdL9ofIlwLLQbsc9ZiRky05HJ9iwJfJb+mRA4vpbl4W5jACcZTl0hEYFVbZ82cqfDVv2R32YL3/rE5a99nFYmcVswNNkpaPFTtsMJ2ZN5ciRYbqjDW9Hvdxww9NhDz0hoLq6CE97UHZaXi5f+fPBul1RO87Hnl7Dqg/CZypnZpjwtNnpaLPT3CB/bef9+4+yd09k1+DhgwPc/PtlYWVqtYq5DaXHRgVyFwUKBAJ8sOLjqPX3LH6c7vXh2VDzy7JxndOA6+x66joq0cqM32w93MuBKAvkfLTvALe88W5YmUGjwWMtZWGFlS6njTyLPNegNzDOrpFoHl6Jp/fdEuYyAsjRl1JpaaUitZViYwWqGYwKpnJobDPjk5HdoTuG32FD/9KwMp3KRKmpFavZQ6m5DaNangu4v2eArWsiy4u94z5u+vbNYS4jgMpWO+5zW3B/pRl7vfxYxgfdexnzRp438djK9az4KFx2mmE2Mq/GyoJaG2c0VPz/6TISQvQA0ZemOj7ZQOR55yc/p6rtp6rdoNh+olBsjz9lkiTFlGo2lteS9wGnEMIK7AMuAr49bZulwP8GVgJfB1ZIkiQJIZYCDwoh/kAwqOwEVhFcEePzjvkZYj2pSAghVsfaS55snKq2n6p2g2L7iUKx/cTyuR2CJEl+IcTlwIsEJaJ3S5L0sRDiemC1JElLgSXAfUKIbUAfwQc8oe0eBTYCfuCHUmieeqRjxv/0FBQUFBRi5ZRKXTEbTuXe+1S1/VS1GxTbTxSK7SeWL1Jyu5jygZ+knKq2n6p2g2L7iUKx/QTyhRkhKCgoKCgcny/SCEFBQUFB4Tickh2CEOJsIcRmIcQ2IcSVEer1QohHQvXvCSHKp9RdFSrfLIQ4K9Zjnqy2CyFKhBCvCiE2CiE+FkL86FSxfUqdWgixVgjx7KlkuxAiXQjxmBDiEyHEppBE+1Sx/ceh62WDEOIhIUTc1/SUa7cQIit0TQ8LIW6dtk+zEGJ9aJ8/CZGYJEXxtl0IkSKEeC50rXwshPjc3G0nBEmSTqk/gqqkbsAG6IB1QM20bX4A3B76/yLgkdD/NaHt9YA1dBx1LMc8iW0vAJpC21gI5og6JWyfst+/AQ8Cz54q10yo7l7gstD/OiD9VLCdYFaAHYAxtN2jwCUnkd0moBP4HnDrtH1WAW6C0vVlwDkn2Xce0XYgBVg45Vp5MxG2z/bvVBwhtAHbJEnaLkmSF3gYOH/aNucTvFkBHgNOC71JnA88LEnShCRJO4BtoePFcsyT0nZJkg5IkvQBgCRJQ8AmoqQBOdlsBxBCFANfBu5KgM0Js10IkQbMJyi5RpIkryRJ/aeC7aHtNIBRBCeSpgDyEm4lwG5JkkYkSXoLCEsfIIQoAFIlSXpXCj5Z/wpcEGe7E2K7JEmjkiS9GvrfC3wARF445ARyKnYIJzS30ixJhO3HCA1bG4H34mjzZ+yK1j7ybL8Z+BmQyIRJibDdCvQA94TcXXcJIeKzGkqCbZckaR/wO2A3cAAYkCTppZPI7uMdc2qK1ZPxPv1chBDpwFeA5bO2NM6cih2CQgSEEGbgceBfJUkaPNH2xIIQ4lzgsCRJa060LTLQAE3AbZIkNQIjhFK1nOwIITIIvuFaCWYQMAkh/vHEWvXFIDQiewj4kyRJ0dddPUGcih3CPqBkyufiUFnEbUI/QBrBhXui7RvLMeNBImxHCKEl2Bk8IEnSEwmwO8yu6e1H2iZG2zuA84QQOwkOyxcJIe4/RWzfC+yVJOnT0dhjBDuIeJMI208HdkiS1CNJkg94Amg/iew+3jGnullOxvv087gD2CpJ0s1xsDP+nOggxkz/CL6ZbSf4dvNpwKd22jY/JDzg82jo/1rCg2zbCQaQPveYJ7HtgqAv9eZT7Xuftm8XiQsqJ8R2goHBytD/i4Hfngq2E8xE/DHB2IEg6Au/4mSxe0r9JXx+UPlLJ9N3/jm230DwxU2ViOs8Lud+og2Q+YN9iaCaphu4OlR2PXBe6H8D8DeCQbRVgG3KvleH9tvMlCh/pGOeCrYTVDRIwEfAh6G/uN8kifrep9R3kaAOIYHXTAOwOvTdPwVknEK2Xwd8AmwA7gP0J5ndOwnmRRsmOBqrCZW3hGzuBm4lNLn2ZLed4ChDIij6+PQ+vSxR17vcP2WmsoKCgoICcGrGEBQUFBQUEoDSISgoKCgoAEqHoKCgoKAQQukQFBQUFBQApUNQUFBQUAihdAgKCgoKCoDSISgoKCgohFA6BAUFBQUFAP4fQI+dhweh+SsAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "import numpy as np\n",
    "Q = 1/rho*grad_u\n",
    "plot(Q,)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 47,
   "metadata": {},
   "outputs": [
    {
     "ename": "NameError",
     "evalue": "name 'assamble' is not defined",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mNameError\u001b[0m                                 Traceback (most recent call last)",
      "\u001b[0;32m<ipython-input-47-080bf2ec9c68>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m\u001b[0m\n\u001b[1;32m      1\u001b[0m \u001b[0menergy\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;34m(\u001b[0m\u001b[0mQ\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m----> 2\u001b[0;31m \u001b[0mE\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0massamble\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0menergy\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m",
      "\u001b[0;31mNameError\u001b[0m: name 'assamble' is not defined"
     ]
    }
   ],
   "source": [
    "energy = (Q)\n",
    "E = assamble(energy)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}