{
 "cells": [
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Copper wire heating - 'Hello World!' example\n",
    "***"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<img src=\"fig/copper_wire.png\" alt=\"Hello copper wire!\" width=\"600\"/>"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Physical formulation\n",
    "\n",
    "Let us consider a thermally insulated cooper wire $\\Omega$ of length $L$ and radius $R$ with a straight centerline whose direction coincides with the $+x$-axis.\n",
    "The wire is subject to a voltage drop $U$ connecting a heat pump control unit, which is located in a house cellar with the temperature $T_\\mathrm{c}$, with an external heat exchanger located outside at the ambient air temperature $T_\\mathrm{h}$.\n",
    "Our goal is to compute the equilibrium temperature distribution of the wire for the given parameters."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 167,
   "metadata": {
    "hide_input": true,
    "nbsphinx": "hidden"
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The tikzmagic extension is already loaded. To reload it, use:\n",
      "  %reload_ext tikzmagic\n"
     ]
    }
   ],
   "source": [
    "%load_ext tikzmagic"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 168,
   "metadata": {
    "alignment": "center"
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAlgAAAC9EAYAAAA/4N+HAAAJJmlDQ1BpY2MAAEiJlZVnUJNZF8fv8zzphUASQodQQ5EqJYCUEFoo0quoQOidUEVsiLgCK4qINEWQRQEXXJUia0UUC4uCAhZ0gywCyrpxFVFBWXDfGZ33HT+8/5l7z2/+c+bec8/5cAEgiINlwct7YlK6wNvJjhkYFMwE3yiMn5bC8fR0A9/VuxEArcR7ut/P+a4IEZFp/OW4uLxy+SmCdACg7GXWzEpPWeGjy0wPj//CZ1dYsFzgMt9Y4eh/eexLzr8s+pLj681dfhUKABwp+hsO/4b/c++KVDiC9NioyGymT3JUelaYIJKZttIJHpfL9BQkR8UmRH5T8P+V/B2lR2anr0RucsomQWx0TDrzfw41MjA0BF9n8cbrS48hRv9/z2dFX73kegDYcwAg+7564ZUAdO4CQPrRV09tua+UfAA67vAzBJn/eqiVDQ0IgALoQAYoAlWgCXSBETADlsAWOAAX4AF8QRDYAPggBiQCAcgCuWAHKABFYB84CKpALWgATaAVnAad4Dy4Aq6D2+AuGAaPgRBMgpdABN6BBQiCsBAZokEykBKkDulARhAbsoYcIDfIGwqCQqFoKAnKgHKhnVARVApVQXVQE/QLdA66At2EBqGH0Dg0A/0NfYQRmATTYQVYA9aH2TAHdoV94fVwNJwK58D58F64Aq6HT8Id8BX4NjwMC+GX8BwCECLCQJQRXYSNcBEPJBiJQgTIVqQQKUfqkVakG+lD7iFCZBb5gMKgaCgmShdliXJG+aH4qFTUVlQxqgp1AtWB6kXdQ42jRKjPaDJaHq2DtkDz0IHoaHQWugBdjm5Et6OvoYfRk+h3GAyGgWFhzDDOmCBMHGYzphhzGNOGuYwZxExg5rBYrAxWB2uF9cCGYdOxBdhK7EnsJewQdhL7HkfEKeGMcI64YFwSLg9XjmvGXcQN4aZwC3hxvDreAu+Bj8BvwpfgG/Dd+Dv4SfwCQYLAIlgRfAlxhB2ECkIr4RphjPCGSCSqEM2JXsRY4nZiBfEU8QZxnPiBRCVpk7ikEFIGaS/pOOky6SHpDZlM1iDbkoPJ6eS95CbyVfJT8nsxmpieGE8sQmybWLVYh9iQ2CsKnqJO4VA2UHIo5ZQzlDuUWXG8uIY4VzxMfKt4tfg58VHxOQmahKGEh0SiRLFEs8RNiWkqlqpBdaBGUPOpx6hXqRM0hKZK49L4tJ20Bto12iQdQ2fRefQ4ehH9Z/oAXSRJlTSW9JfMlqyWvCApZCAMDQaPkcAoYZxmjDA+SilIcaQipfZItUoNSc1Ly0nbSkdKF0q3SQ9Lf5RhyjjIxMvsl+mUeSKLktWW9ZLNkj0ie012Vo4uZynHlyuUOy33SB6W15b3lt8sf0y+X35OQVHBSSFFoVLhqsKsIkPRVjFOsUzxouKMEk3JWilWqUzpktILpiSTw0xgVjB7mSJleWVn5QzlOuUB5QUVloqfSp5Km8oTVYIqWzVKtUy1R1WkpqTmrpar1qL2SB2vzlaPUT+k3qc+r8HSCNDYrdGpMc2SZvFYOawW1pgmWdNGM1WzXvO+FkaLrRWvdVjrrjasbaIdo12tfUcH1jHVidU5rDO4Cr3KfFXSqvpVo7okXY5upm6L7rgeQ89NL0+vU++Vvpp+sP5+/T79zwYmBgkGDQaPDamGLoZ5ht2GfxtpG/GNqo3uryavdly9bXXX6tfGOsaRxkeMH5jQTNxNdpv0mHwyNTMVmLaazpipmYWa1ZiNsulsT3Yx+4Y52tzOfJv5efMPFqYW6RanLf6y1LWMt2y2nF7DWhO5pmHNhJWKVZhVnZXQmmkdan3UWmijbBNmU2/zzFbVNsK20XaKo8WJ45zkvLIzsBPYtdvNcy24W7iX7RF7J/tC+wEHqoOfQ5XDU0cVx2jHFkeRk4nTZqfLzmhnV+f9zqM8BR6f18QTuZi5bHHpdSW5+rhWuT5z03YTuHW7w+4u7gfcx9aqr01a2+kBPHgeBzyeeLI8Uz1/9cJ4eXpVez33NvTO9e7zofls9Gn2eedr51vi+9hP0y/Dr8ef4h/i3+Q/H2AfUBogDNQP3BJ4O0g2KDaoKxgb7B/cGDy3zmHdwXWTISYhBSEj61nrs9ff3CC7IWHDhY2UjWEbz4SiQwNCm0MXwzzC6sPmwnnhNeEiPpd/iP8ywjaiLGIm0iqyNHIqyiqqNGo62ir6QPRMjE1MecxsLDe2KvZ1nHNcbdx8vEf88filhICEtkRcYmjiuSRqUnxSb7JicnbyYIpOSkGKMNUi9WCqSOAqaEyD0tandaXTlz/F/gzNjF0Z45nWmdWZ77P8s85kS2QnZfdv0t60Z9NUjmPOT5tRm/mbe3KVc3fkjm/hbKnbCm0N39qzTXVb/rbJ7U7bT+wg7Ijf8VueQV5p3tudATu78xXyt+dP7HLa1VIgViAoGN1tubv2B9QPsT8M7Fm9p3LP58KIwltFBkXlRYvF/OJbPxr+WPHj0t6ovQMlpiVH9mH2Je0b2W+z/0SpRGlO6cQB9wMdZcyywrK3BzcevFluXF57iHAo45Cwwq2iq1Ktcl/lYlVM1XC1XXVbjXzNnpr5wxGHh47YHmmtVagtqv14NPbogzqnuo56jfryY5hjmceeN/g39P3E/qmpUbaxqPHT8aTjwhPeJ3qbzJqamuWbS1rgloyWmZMhJ+/+bP9zV6tua10bo63oFDiVcerFL6G/jJx2Pd1zhn2m9az62Zp2WnthB9SxqUPUGdMp7ArqGjzncq6n27K7/Ve9X4+fVz5ffUHyQslFwsX8i0uXci7NXU65PHsl+spEz8aex1cDr97v9eoduOZ67cZ1x+tX+zh9l25Y3Th/0+LmuVvsW523TW939Jv0t/9m8lv7gOlAxx2zO113ze92D64ZvDhkM3Tlnv296/d5928Prx0eHPEbeTAaMip8EPFg+mHCw9ePMh8tPN4+hh4rfCL+pPyp/NP637V+bxOaCi+M24/3P/N59niCP/Hyj7Q/Fifzn5Ofl08pTTVNG02fn3Gcufti3YvJlykvF2YL/pT4s+aV5quzf9n+1S8KFE2+Frxe+rv4jcyb42+N3/bMec49fZf4bmG+8L3M+xMf2B/6PgZ8nFrIWsQuVnzS+tT92fXz2FLi0tI/QiyQvpNzTVQAAAAgY0hSTQAAeiYAAICEAAD6AAAAgOgAAHUwAADqYAAAOpgAABdwnLpRPAAAAAZiS0dE////////CVj33AAAAAlwSFlzAAABLAAAASwAc4jpUgAAAAd0SU1FB+cHDw4ENRWXYfEAAHQGSURBVHja7Z11nJTV98c/M1ssbNEt3SlIK1IKSIctLd0gII1IKiiICAYpIaCIKGAAIqF0SEpJd+cuW78/zmeQ2d/sAn5Z2JXP+5/Pa3eeO8/M8zxz7z3nnnsOIIQQQgghhBBCCCGEEEIIIYQQQgghhBBCCCGEEEIIIYQQQgghhBBCCCGEEEIIIYQQQgghhBBCCCGEEEIIIYQQQgghhBBCCCGEEEIIIYQQQgghhBCJEeci01TTTf1OxnJgGpMkQaZeHXXthBBCCCGEEEIIIUQ8kuFJ0/camW58zvSLgqaZ07sfH/ys6Zuhpilm6hoKIYQQQgghhBCPH166BCL+SfKO6aAw02ejTc8OMS3QwjSv6//XTGtfNj3X2HTTPL7h77qmQgghhBBCCCHE44O3LoGIf/JeMT12wLTOCtOrX5sGvWhaYo1pm7dN96YwXTjfNHqUrqUQQgghhBBCCCGEiAcCXzUNaOn5dWcx04ovmbZjhFXKQF07IYQQQgghhBBCCPEI8WllWv+iaZd8pin26toIIYQQQgghhBBCPJY4GQmV5itTn/cezedIksu0UQXTlv6mgeV1j4QQQgghhBBCCBETpy7B40TwEdOX0piGpH245082ybQhtwaG1TSdwb+vrnQ/3qczNb/unRBCCCGEEEII8fiiKoSPFU9cMO3HLXorkpqeWR2/5w0uYfr6bp7vhun3GUxv/ey5XTUmdQ/YZnpys+6hEEIIIYQQQgjx+KEqhI8VTx02LRxlWoaOpO3D4+d8fi+YdqaDKtUQ09/Xm3pNNo14ztTJHFgluaWwVjnT98J074QQQgghhBBCiMcXObAeC7xKmj5/0zQgk+lzjUynVzINXf6Azpfd9OXGpmlym67heTpvNE3OCMCthXjcddNSX5uOTWn693bdQyGEEEIIIYQQQoj/NLnfMd3d2zQ62vTAVdOCex7s+fKWMe3L86T+IMbrJ0wX/2J6vrvptkOmLUub+r6veyeEEEIIIYQQQgjxWNCGSdAjW5i6HFjRpUy7+z3Y8wVfNk1TOe7j0pwzfYZbBDMxybszk+6ZEEIIIYQQQgghxGNBQBXTJRtMbzuuYujq30xT9tI1E0IIIYQQQgghhBAPkepMmn6eSdtjc2DdYI6qRmN1zYQQQgghhBBCCJHQUBL3/ySBjKhqssg0hcM0ur+pg8nRo6+Z+v9l+kaI6WJf0/O3dC2FEEIIIYQQQgjxqPHSJfgvkiLUNOQN0zVZTbeyCmHG8qbTlpr+UNt0Dx1b+1qZXpmtaymEEEIIIYQQQgghHiJZK5iu3m6au4quiRBCCCGEEEIIIRI6Tl2CxwrmunI05N8RuiRCCCGEEEIIIYRI6MiBJYQQQgghhBBCCCESNHJgCSGEEEIIIYQQQogEjRxYQgghhBBCCCGEECJBIweWEEIIIYQQQgghhEjQyIElhBBCCCGEEEIIIRI0cmAJIYQQQgghhBBCiASNHFhCCCGEEEIIIYQQIkEjB5YQQgghhBBCCCGESNDIgSWEEEIIIYQQQgghEjRyYAkhhBBCCCGEEEKIBI0cWEIIIYQQQgghhBAiQSMHlhBCCCGEEEIIIYRI0MiBJYQQQgghhBBCCCESNHJgCSGEEEIIIYQQQogEjRxYQgghhBBCCCGEECJBIweWEEIIIYQQQgghhEjQyIElhBBCCCGEEEIIIRI0cmAJIYQQQgghhBBCiASNHFhCCCGEEEIIIYQQIkEjB5YQQgghhBBCCCGESNDIgSWEEEIIIYQQQgghEjRyYAkhhBBCCCGEEEKIBI0cWEIIIYQQQgghhBAiQSMHlhBCCCGEEEIIIYRI0MiBJYQQQgghhBBCCCESNHJgCSGEEEIIIYQQQogEjRxYQgghhBBCCCGEECJBIweWEEIIIYQQQgghhEjQyIElhBBCCCGEEEIIIRI0cmAJIYQQQgghhBBCiASNHFhCCCGEEEIIIYQQIkEjB5YQQgghhBBCCCGESNDIgSWEEEIIIYQQQgghEjRyYAkhhBBCCCGEEEKIBI0cWEIIIYQQQgghhBAiQSMHlhBCCCGEEEIIIYRI0MiBJYQQQgghhBBCCCESNHJgCSGEEEIIIYQQQogEjRxYQgghhBBCCCGEECJBIweWEEIIIYQQQgghhEjQyIElhBBCCCGEEEIIIRI0cmAJIYQQQgghhBBCiASNHFhCCCGEEEIIIYQQIkEjB5YQQgghhBBCCCGESNDIgSWEEEIIIYQQQgghEjRyYAkhhBBCCCGEEEKIBI0cWEIIIYQQQgghhBAiQSMHlhBCCCGEEEIIIYRI0MiBJYQQD5ee7H638u8quiRCCCGEEEIIETfeugRCCPEwSPa8aeVtphkvmR7NYfpjdtOIzzy39x1rWqK86aHNpsdb6NoKIYQQQggh/usoAksIIeKVJP6mDWeZOiaZnj9rOqSXaYX5cb9PqRdNp/9sWr21rq0QQgghhBDicUEOLCGEiFdKrTcNvGL6XVPTk1lMM3c2zfVa3O/z1EzTtHRcHX9a11YIIYQQQgjxuCAHlhBCxAs+/UxLTjFdxwgqLDUpc4Z/1zM50NPz+/iz/VP5TI9dMN0nB5YQQgghhBDisUE5sIQQIn5grqrfkpjumGEaMt205m7T3cxltX6X57dJzq2FT35ruu+g6dH5usRCCCGEEEKIxwU5sIQQIl4IZ9L29TH+X/ZL04JM5v7hN6aXnvH8PjlHmmYsZTqbVQvDdukaCyGEEEIIIR4XtIVQCCEeKuULmHoxqfvKl+M+vhS3CnoxomvjZ7qGQgghhBBCiMcNObCEEI8rXtQd1Czxe7qAoqaFGWn1d5TpwemxfLzCpsVeMD09xPSvDLp1QgghhBBCiMcNbSEUQjwmOJhrKtVJ04JzTKPeNl3Vl3+3ip/zJxttmj7U9FwL0yu+no9/orZpYTqstrUxPb2CB6TWPRVCCCGEEEI8LigCSwjxX4FJ0f2GmaZtZlq5oum73Io3/7RpBzqQDjMXVXw5rlzcqGV6mpFXXm3ZDd9yP86rnGmNRqY5tpquS2d6XY4rIYQQQgghxGOHIrCEEIkM1w66zHQApaYjPh+r/RViv1a0tGmm86YhU00P5DLtscD00NGH87mv3jCd1pTnf8e0TgrTHWdM8/Y0fXGZ6U1uIdxWX/deCCGEEEII8bgiB5YQIpFRYKrpmEDTrNdNk2blAT1iNGBk1UXmunqHEVnrnn00n3/Bq6ZHL5sWHWiah9/nMCO1rmU3PbnQdGtVvsESPQNCCCGEEEKIxw1tIRRCJDLWMFLplyumfve4pW4mc2DNrWMa/c3D+by+jPCqx6TxL3ML444Jpp9v5Od7zpR+LeQfbLqEWx3PrdG9F0IIIYQQQgjxGJCVkShr8pjmrqBrIhIvT0WanmSuqOhoz7plsmmuRo/mc5b+wvQ4P8/Weabpjrkf51PZdPQU0/XMjZVni+61EEIIIYQQ4nFHEVhCiMTSXTF3VYkSpv0ZwXTuSeoG9+OvLzIdywinfV8+ms+d7CnTY9lMRzGJ11k6lH34ud5Ixu+Xit/vD9O/ntS9F0IIIYQQQjzuKAeWECKBE1DJtNlV01afm26mo+fV6qa18pkOqWD6C1+fX/PRfv5N3DI47TXTZIwcq9/NNHt60wzcKtgto+nGY7r3QgghhBBCCCEeQ7SFUCQmcq8wnfyM6U5W8evwnWlgJvfjUzPiajIjnqr+kLC+T7K5piX3m77wkmnhYaZJ+/PAQrr3QgghhBBCCCEeY+TAEgkZn3WmtVkl8I/MpkveNi3N7ObOvXG/T4qRfL9auqZCCCGEEEII8d9AWwiFEI+YdExu3p1b56qx+t7C903H3TQ9FXxv73ehl66pEEIIIYQQQvy3kANLCPGQcYw2LXvQtL+/acjHpn3Hmy5h1cDwfrpmQgghhBBCCPF4IweWEOIhEcBIqiZlTFszufq6I6YdqpnuL88GclwJIYQQQgghhAAgB5YQIt7Jn920zwbTgnNMJ7QxnZHN9OpmNtisayaEEEIIIYQQ4k7kwBJCPGD88prWOWD6FqsGnspv2srLdH1pNmihayaEEEIIIYQQIi7kwBJCPCAydDftzH6lNh1Vc4aYflzA9JyrwVpdMyGEEEIIIYQQ94IcWEKIf4nXq6Zl3zYdwCqBSemo6pHZ9KdvTcPz6ZoJIYQQQgghhPg3yIElhLhPQr4wbbnXtMnXpqs2mb7/hOnB3Wwgx5UQQgghhBBCiP8Jpy6BEOLeyD/cdPww09dbmY52mHZrbXrbcSWEEEIIIYQQQjwQFIElhIgF/+dN65w0fSul6dFipm/WMd30s2l0Rl0zIYQQQgghhBDxgRxYQogYZN5i2rOlaQVuFZyz2nTCVNPzAWwgx5UQQgghhBBCiHhFDiwhHnu86IAqv8q0f35TnzdMuzESa3lD08gAXTMhhBBCCCGEEA8TObCEeGxJUcS0ZTbTxtNNf51q+l410yPz2aCYrpkQQgghhBBCiEeBHFhCPD4UNCla3XTAb6ZPvGo68gXTeRNMb4bokgkhhBBCCCGESAjIgSXEfx7/7qYNgky7VTA94G/aJIfpzvRscE3XTAghhBBCCCGEEI+IrK1M1+QxzV1B1+Q/fb/nmI5nzqrd9U170XGVPLuukRBCCCGEEEKIxIBTl0CI/wpeP5hW9TL9cpdpnmdMOww3Hf2X6cWDumZCCCGEEEIIIRID2kIoRKIn5W7TNs+avlHL9OfvTN87a3o8j66VEEIIIYQQQojEiBxYQiQ++poUS2var5xpFkZUDVli+o2faWh9XTIhhBBCCCGEEIkZObCESDQkK2D64sem3Yqb/jnGtCkdWTv4/+jpumZCCCGEEEIIIYRIZCiJe+IkZxXTT9eZ7uT96/69adAVXSMhhBBCCCGEEP9lFIElRML7We4zrfqtaZ98pteXm7bbb7oqh2lUkK6ZEEIIIYQQQoj/tKWsSyBEQiH1FtM2jKh6Y7DpgnDTMeNMT0ayQX5dMyGEEEIIIYQQjwNOXQIhHhWOJqYlC5t+/oHpc3Rg9Qkx7VfT9GQrXTMhhHigJKWG6FIIIYR4BHBHBWbqUghxdxSBJR4Bvm+bls1uWqCZ6VU6cPD9f+SLOqiMmIoqaRpdzTTXStNGA02T1DWdtszUn1sJX/mJb3eZ7xeqZ0gIIe5twGG/W4v98Dv8P/vdVGv4+g7Ti1wwcBzSpRNCCPEA8Oc4k5JmwTOmzjJ8uYHp0oumf+uKCREHcmCJh4BfVdOSyU2b9DetVcH0egvTQ3Tc4OnE/oU5ULEqYPTvphkY8Zj9E9MrN/m9L/E6nOV1qmda2hVxlZn6lJ4lIYTwiGvB4Ab73QOmUZ1NA3aZ5l5rGtiU06DfTPf2MT0+mW/XWpdUCCHEvyCKYuNKkugvAcCRJ2QOAHi3yNsGAML7+tk4dZXDFurPMpUDS4i4kANLxANJGCH0FA2KpjQYSrc13UNHzI7GphtoMIx+kYZDksQ9YEW/Z+qbybT2VtP2hUyXcWQaw7//PMh2IXyfPdTn9SwJIYRHYjisQId/FCOoghnJ++xLpq+/ZXq6oemsrqbVGHk1iRGyXxQz9YnSJRZCCHEP4xBTf0S35DjUCgD896SuAADOftV/BoDoiq/nA4CIc/tHAgD819v4dfXZSXwfRmDhd11aIWJHDizxAPA1gyHjUzMAwJGheTcAiD5SrpvZC3+msuN6pDBdXdR00HbaH6dNzy77b1yPtOlN25Y2bWADFea/YjqGWwrPDovlDV7XMyWEEPdDUo4flf40bcotglVDTP24db0vcw+Oo8HwZBbTazQ8LlbStRRCCHH/BNwEgGQXaowGAO+LrcoCQNCZ3LsB4FTU5GwAED5wOu0EJ6uNF8tg6hioayiEEG5k5Qr1mjymuSv8u/fx5d7tctVNJ84CgFStT0YBQLq0v2y0/1fjCndgfgBwML7K52tvq5730QeH7T8DPknc19XZ1LQMHXKLfUz/6GJa5zVet0V6BoUQ4kGQbL3pC9yS/vUZ0+3MNfLpENOl35hGR5vuHmua+VfTlVxQaLtc11QIIcS9E8SF6fq5AMC5+4etAFBk47UeANB95KmWADAgQy1G/jqbA4DzDVf7nOlMNxY1rVxO11QIIdz4tw6sJLtNy71Kw6AKAKQ6ebwcALxMw2Ch40JPANgcWm4fAOyhW8o7l+mP503PZCgyEwAaJt1R0f7Tb33ivJ6BW0w7zDbdmdN00kRe3w565oQQ4kEQQAdVDUZKzc9teo797WfMdfUUTQMv5rBKR0fWxrSmLkfWGBbH2MIFhza7dY2FEELETnLm7G240PTbdwEg31cXOwHA4D42vmwffHwlANx01EwDAHurxHgXLuinOJenOAD4/vUnd2xUuqRrLIQQbtyrA8t3tGm5Rqafh5turwwAZWvu+xwAFv0c/QoAXF9zMczMgqbBNA/q3Wkm/KMFVgBA+HObNgFAt/muV/rmTyQXsIhJ/v2mXw7idSlh2n6VabIWetaEEOJ/IaCHac0ppvMmmB7mincoi1tEMxXC4EBT76Oe368Ni4hc62Ua2c99hGrZTNdcCCEEgFEmQSxC1XCz6XcsBrJ2LwC86r2mMgDsXhGRBwCi3j/+rY0qNRp4toNc6psKAH7cMfgTAMiV9CJzA1dcp0svhBBuxHRg5WJH7MWqd+WYI+TTTqa7mprOCgOAVqH9NgHA5fSHZ1sHfOUZ02YFTZ2FPHfUebeYbkgBAOErbxQBgNf67R9qR/S9mbCvmw8/ZwPmRlnHZL/f0RAqWdfU2U7PmBBC/BuSsv+szuIVc1abbg8y/YRbzkvOM+1dwDRsBjWfaX86tLzfcH9/r+GmHRhJfO1395GqpbZuCCHEY42rOm2D66YLF5uuf8J0wFQAmL6pQyMACD9/orCNIqcsSXt0zb6e7SCXeo807f8LAOz9NWwgAOTpe/mkHVFxvu6BEEK4kZVJA3+nAdCYWzEmMqfVDq40fM1cTs8fBIAd3Z79yDrc7WNNL081bUpHFmZ47qgLdDfdUNX0+kAAuNmvy0oA8Mk3hsnKB9RMmNcrAx1rI5jra89l00GHTNOm0DMlhBD/hoD2pjW4oDL/a9PLSU3PWBEQ1K1s6lXNvb1PUVOuu+Amt3SH8R8DWAXWd5p7O+cA03bs32/S4dWquu6JEEL898nO0SPFlylaAkCygS+xuMd3tEfWLTF9h0VBch0GgNDsz/c3e+bIANNTT5nWcou4cnFHxNV60wF7TUM/BIAdp3f/AABB6/9iCpJKS3V3hBDCJvpcSajZ3PQ4qy/9vc10LsuQ10hmmmScdbBP0mG1K8D0coRp0/XuHXNMLbDMdNMvpjfqmHaYCgBlj/vUsCM/6ENDY2LCuE7Od00rciD7hXvSV1wzrc09795F9EwJIcT9EMBk6tXZz85jUtudjKCazEjXg8HuI8qiOaaZFnh+Xy9uDezFceTmOdNwjjMDudXQt3mMcZFbv9dwy2HbubpHQgjxn6QrAKRdEbwJAIKPvTIZAJKe+3EPAPhu3xZqh71DeygvI4HR2EajyoGmxzKbnuxgWnNQ3PaQ71DTgcdNw4+Z/lIHAHq8V2mWHbme46GSuAshHlt8vgQAlCv9PgBg/oRjAADngZ72+oXBpo34d7JWAHDJ39XhPpnP1JKsR0df2W7avKepM8BzR533LdP10aY33zPtPMXUi1tBvNhhfzDO9FE7sIIYAdCZVah2c6vlp8ypkq2VnikhhLgf/JnL4wUbjzCPCwLnWY31S5YNL86t7D5Myl6BxT0OVXAfYZZ8ZJrhWCzjHnOI9GlrGnaV+pJpvxdMXWVFvGgwrGYyd1UhFEKI/xYBZQEgYPdLRwEgZNDSpgCQNvJQDQDwKjNkjR2Xl0VAvAoAwL6/XaNP5dGmR/xNT9OBVSsLR6ecnu0h71qm/Zkb+FZd0xXZTJ/YYUdm32C6kXZQ5ZK6Z0KIxwRfbvkrbY6pUhN6AEBKHG8JAME9Fne311t3MV3HFe+codaRJnvDtPZ10z9rmF7h/5sx5NXxrueOOt9k03UXTG8wiWGniaZeGwGgV2rX5/Vixz36F9MBkx7yBaPBUrC86XRGXP35rWnzjBz4KuvZEkKIeyHZV6bVuTVvHh1C27gl8FyA+8gxj8UwUv7l+f0qMULq76vu7RZx/HiG1QN95rm38zlr2vuQaSgjq0JZZKMfV7j9ucK+uqCpHFhCCJG4CeaW8drbAMDbyzZOFC12PQMAdNh1xAkAbVNW7W3HORvf2fofu+a5rKZHppme5sJ97Zc920EuzUuH1ZAnTcP+MF3OBf1sXgDw0g+uM+ZkEauNDUwrl9E9FEL8R/FjUsEyXFH+5CsASPnn8cMA0HBQ1AgAGOY9ZSMAtFwVzKTsWYpzwp4DAMq3L/aZdajjvEyjGTl16WPTZtPi7qgLMLR2I5MX3hhk2vGiqeOpO49f8Kbpsb25OwFArdQb+Ep/34dz3XxXmDbiCsemEabzmTyxZE89W0IIcS8EFDOtylwh85hLcTcjq8YzJ1VJ5l6sw616x593H0nmceKeeqjn81QYZnpgm3u76+zPu+wx9a7n3s6LVXTfZrXBm8dNb50yfYeRwBv4edpm1D0VQohEQWGTIFZJr/8lADjazm8PAHl7XC4NAANH2XixrfnBAwAQurZSSgA4XTQ2u6bKKtNjLFJ1qpJprXlx20NlZ5nu+8ndnloeZPrEmTuPD9hrGrgw9+sA4H1068f2n8reurVCiP8IPnRUle5g+klFAAjwtxxW9d+PWgsA378XXQ0AzqedXB8AIk8HtwKAcn6u98myFAACj61vAgCbKi/ZZh3qLYa6/rjc9MUSps4Onjvq21sFGXF1c6lpl+dNvTJ4bpd1FQDcurqiDwB0ez36or3St1T8Xr8MrOoxmjm+9jAXSx+G7oYs1DMmhBBx4U9D4YUkpl9dMj3X1b2nX0RHVZqznt+nJpOxn4hwbzeXWwRT/ea53dMnTGekNL3F5Lo36Hjq/JypM0ayd2++Xx+2D+PnCmO/H0XHVstRusdCCJGQCWC12nrsx799x3RdXQBoVGfLVgDYeTHqMgBEVjtc0kaZShYRFRrTLnFpJUZIHT5pejqvaS0uyCOV53alvjTd62d67pbpqIym2c95bucsCgDfftp9IgDkyHuekVgV9+oeCyESKT5XTMtwZXoC92j/WRsAUvde2hEA5m47PhQALq91efqnvGMakvnOjrLpVNO6L+SYAQBftf17OwCEJ3O1m0cHVuqmd7Zz8U+Hm4crDGu5knCTKxJdWph6v+S5o87S3nT5UAAI3xY2EwCafHuUVaf6rn2w18+LhlblVKY/FzVdOcH0eVcy9mF61oQQwhNJWfzj+SOmXzFSdjtzSo2/ZdqVkVbnWOwiuorpJOY4DIqlLHitH01PFXMfMeaw/065xHM7/7ymI5hsN5IRVNeYW6vj76aOye7tfJnzqi8dWmG+7udtqeS5QgiRoAjklvG63Dmx4EPT9Uzx8c4JAJh15e0tABA++3wy69WPFjet0s6zXeLSykNMD/c2PVPYtE7puNuVHmO653fTiyxy9dpmHnfdczvHCtNOPQBgb/HrbQEgb8kr/F4V5+ieCyESCX4HTEszmfgE5mLazhXjr1cDQOEkr+8EgPNvf8OqfpGnTadkMg329dxhBjxt+vEI0yh2oPNqmqbKfufxLv5pn3+qqSs5+41f2QFHskP+2vN5syc3/ZVVN8LeA4CbFwZuBIAkXuOcduSAwg/mOgaPNX17kOkuDnyf0EDKOUTPmhBCeCIZI6Cq/mE6j5FWF+qYfk/HT0muEPvUcm//CnNfnVvlPhJM5kJM8BTP5631mumJHu7t5jIpe6qNsYybw02H0XEWzq2LN6qbdmaydlfSdhferFrYmwso4dyy0cqpZ0AIIR4lwYyoqsf+fz53SmzgFrt3Ob7kzQIAt/LWqW+jximmPDk80rTyWM92iUur5Dc92sb0NJOr1z4etz1UlpFV+xgwcLGU6WsfxN3Ou7Jp5yOm13cBwK7FhzMAQPKT+7hwU2mangEhRALFh9XtSr1t+glzdOzkBP4bbvqr0QAA5gSmsyoYhyb5WMcX8Y3ppE2mQfs9d9BJ+fpHawEg+tfIRQDwU5vz7wBA//35cwNA3xSxdfB5WT1wPTv0mwyh7cy94c7YtgpeNl3RyPQWj+vbCgAG/52ESX1Hc+V+QKV/eSFXmBT5xPQrGlybGWnVOI2p72965oQQ4k6SjjGtxqq0swea7mDV1aMZ3Xv27Sz+UTA8ljfkFu1XmpueYw6qaOYY/IKRxUFbPDevsd70eIyk7XNouKT6ynM7XxYBGdbaNHyf6Y2cpp241dC5yb2d3wemf9AR1na0ngkhhHgolAaAVI0C3wCAFEVf2g4APlMWMsXHRvbjQ1eb5meVPserNjrUyGp6ijmnjjBpeqVNnu0Sl1b82/Qwt/adpj1T+2jc7UoxAGAv7aCLdJC9ziqF6Oe5nZP2U2c6uG4ONLWdLCPq1/KyI9czh2RlRQILcQ8oWdxDwY9Vjkqyg3qZ+vyvpqc44R/X33T5dACo+UkYAKBwgeGjAODkmqYFAOD6j4t+A4AxL7xXGQDGP5GUK8/msCre2a80AEzN8VYTAEg9pM0aADgQ+NcBAGgeGv0sAByfmPyCtUv7HAAMdcB0ZubvAOC1meMiAcAvtFhjADgzfHRXAHi9z9deALBzWZo27KJ3A0CPr9J2BIA2hz6cBwD+keXLAsCFWxMmAUDbJdMHAsCq9ametXYBzIESVcg0HasBOlhdCg5qGM+TjQYLc51U5t+dGXp7ibm0+rEM+xau3KesyfeZpWdRCPGY8TP7zUymhRix9Bqr/JXsa7qziGl/JkW/xfnBKEZWFXzddPq7ph244HIgs6mTDqgV50xdDqXBXJBpPojdOnMfDqFD7CbHwQ35eDxX3kdygeel9vwe/Lsfq0Zd3cX/s4rhZ1zASM+tgU0ZKTYsyNQrt+lsjjPe3ILoxQWlYI7Lqbi10eeCHh0hhHggzAAAR3Xv3AAQOL3UaAAInfTmpwAQ9nTa5QAQ+fePddndc+v53p0AcLxixBcAkGF5Lc7nPw82jVhm2oIOr+XFPZ++IvvzKczZmHSdaUtWw12Y2XO7kgtMp3KHTBouzLQPoVlB+wMxdng40vE45uwaznFwF8erZrPtW4Zxq7uT4ywm61ER4u44dAkeBs3owBrCUNiUNARuMLfVgQKmVxYDgNdCLy8ASNM1TS0AyN02z1EACJ0TdhoA/rq56wwAXNpxdSk7vq0A4JXReR4AUmwJeQoA8mXOWwwAbs0NjwaAvVmOdAWAC0XTMefUPjqKrpUBAN+xvqsBIEtYpvkAkHFb5kUAcCngUjgA7Azc2RQAwkdFcUse6gCAz6fexQEgU3T6PwAgy6dZ8wHAlVZXegPArvd3JQeA0IO3uOLtZDn1bNzTHskV9sMcEMCk9ZhtEsW/fXm9grnVJQsjraKfNt2fw/Qy6x06OYCgKjVCz6IQ4jFjvEkA+/sc3KIQTEfPZY5LeznBvuZqx0LfGf1Ns3Krty9zjxxnGfLDnU1vnWe/O5P9cjfTfDx/KtcEnyfYwwWGC+z3o+ioimYOq9zdaTAwcsuLDqu/apueYf8e9Tffl8nYk3CrYn5WLQzkVsebjATbzQWTG4wYK1yW78ecX8f5+Rwv6tERQoj/CVfErjmAageOA4Cse7IsBICItnNCAeBY1/HcSbHncwBImTpyIQCcO+N6m9p0MH3uKsbBSK0WjJz9ZRUAOGjV/rOVrzIXUqawersfcyu2ZuTxgrme27mqvE8dwHGI9kcHbomfmdXz1/WifdWe484I5trdQUdbk/4A8Gaf3QUAYNLwnBzPvuL5etGBtWy9Hh0hYkcRWPELtzDk5laMw9wa0ZNJbsPYQTlOAEAB3+QtAaBzup5VASDEJ8dTALDlzKzhADDi2wXVACA6k4PVCB0vm0YsAoAWqZ4+CABVhnUZAQB/N/xlMgD0DbGV6YiiaWmIdCpvOn0jALyYPXlBAKhTpcezAODs7pwDAAv/HLYYAFaccG3J83mL3+t5AHBOinoZAN7Z2WwSAORsWHYRAKzY+/FrAPBJ219bsp2rPPlkAPCd6FwAAN7PtNgDADcKXGBE1tc0JEAHlNfSO04HoC4HjB2sMvUxq1GdZe4SL1aRcrjex7Wi8aceRSHE40nkZ6ZNaSCcXW6alg6nolxQuc5cVxO5Jc/xnWkU+1c/bnXoyH62JA2BFIyIHcOFhU3sz72ZyzGCWwZzcOGiByfy+Tlhn8gFiKV0KHnz80YUNc3IcbRXDdpBrhV0JnMfyePOPGnqZJ2piFdMK9Jx15q5sPJ0Mh3C752SK/d/0EBa/As/x+d6doQQ4n9i/p3jhyP0rZMAkKbC110B4MOtnVsAwFffhfsBwEe1b7d7wqQmi3N8xoirsBamzTj/X7bqzpP944BKOt10gGsBhePVG5dMzXH1/9v5MuCgN8fLHFE0S7jVfGZ3z1/Ti5HB7V3jEiN6d0SaNp5quicJAEwudvvMC2J8gpl6ZIQQj5rUJsMrmM4KMfWvAQCfD3J1WwFpTT+dbxrZ2XQak5+HrHbvYN27u+jo/MyBtfk502unTF99CgA+7uZqkZVd5spZANDbu9p4O24lc2uFfmbafbCpd9+4z9vgqunZ86b7/zQtcMRTu/WlTa/mq94QAFqUPtHW/tO3uvuRyXnkW4wI+JNJd0dz4Hgiux4tIYS4H0ZzJaAfV6Sr0wCI4IIAd/ah0cq436cAt6xvzOI+ImxgZFe+J+JuP+RT93Y7uDWwWFQsDcqY1GGZ8dP8HtEHTWfQARYS4rl5UladXTbR/bzfc0vJXiaRbxOkZ0QIIeIDf0bcLuICRodWd74648eY9kVNVkc/aQv30cdoHz030v24mOq70nQoc11FOEy//9w0cK/ndl7NTXsnNQ3rZbqCdlGqWHJrOVKadqC9dCM17TEmdy/4jCd7KAM3sBecWiAUAAKr7GbxlEpf61kRQjxqXA4sTrBncwUhCatmBJ00nczkg1GTTafNNQ3x6Aj6p+N88k3TbayGcbWJacsy7FgHAUDq23F2WSYCQNG/tu4EgEP5jvI814eZdj1h6o24B4hXWK3jPD//Pg4UT78S9+etuRMAwmedmAwA3Wq5XulLQ+LJL0y/4V7zTYdMX+cWFX85roQQ4l8xmjkLB3KN2/990/FcKY6ig+cSV8ybMVLrnz2F7uTvYrqukntPv3mQaZH3PLfLxpXtX/3d2+3hQkWZHnF/j7rUU03c289eYJryfc/tnqPj7fAAzyNbSyXPFUKIeMH/iOlibs3u1OXOV0PPu3rjWjlMT9ERdHSn6fPvxW1fJFlgOvxZ04ispotZTCpjUs/2jM+Xpr3ZPqy16W90PGVf5Lmdkwv4HfrQjvrLdAPfr8CSuD9vyp0AsO7NOQUAIG+eK4zcqviDnhUhxKPG5cCigTDrFQCY9G26VdaBTZpnGkWH0zSWWQ3JFHfH9yQ76O2MvLrKKhwtXQ4kN4//DJoTy1aVPAYA+z45y/PcOGjadbGp46znjtqR1/TlA6YXFpru5cBSZpfndi6tQQfXqdoAEN7ryJsAUHXeKlaFWsWtJBtpMM1iksRCl/UICSHEg2B0BdOBc93/n6yu6TgaGJFcSLjMHCPNm5o6q3p+37xMDr9uj3vPv62ZaZFxnttl5Vb6FYHu7fYymXvppXF/nzr8XCffcm8/m0lyU27z3K5yL9NDMc4rB5YQQsQP/8+BVQMAgku7euGaOWkn/GF6xBVxFRy3fZGkjunwG6aRFU2X+JpmqO+5nQ/fvy8dZbdox/zGnSTZKntu50W7rSOrrd9MYrphs2m+CnF/3lR5TL+aDQD7FtorOV+5wB0mFZXyRAjxyKEDa1gjAEj+x/fDAeDyczNbsaNl5NV0Rl4lzxt3x1eIHv+tEabXGEHVuqWpc5DndrkYofVHbQCIGBvREgAi3+gx1f7vnTbu874cYnr+DdMDQ03LLY673Qvc+niCkWFHugDArJPtuWVkGXOOnD5k2pXVRYI76NERQogHyZhBpmPne349gFX8xrM/jmbS9UtMdtsshAdW99y+AKsYboxyHwk2M8luwSqe2+VgNcSVjdzb7eb4UOLnWL4QHWR1mCT4dEP39rO4QJJiiefmVViV9gyT2bcuqmdECCHig9sOLC5QtP/Beuvqb9FOoEPp2BbT547EbV/40o4ayh0h4VzI/5EOpkx7PLfzamD6dhvTsMamq5ua5izvuZ2jmWl7RlrdoF20ebtpwbp3Hu/ijogrpmqZTbst/F0A2DT5x7cBwHftn7wulVRsSgjxyKEDa8hKAKiX7Fo3ALjRMpIhptO3mSbvF3fHV2SU6dYUptdymbbi1j9nB88dbm6uXKzpCQBhWaOiAWDQhBM1ASBb0pxXAeCpHbENEC+FmZ7nHvQDPM8zveP+vDVmcEBaAwDhr59+HQAGFBzLcunLWGXqMnOxjGbklXd3PTJCCBEfjKEj6RKTsjeMpcpRILe6f8LIqegxphc5sW46lAem9Ny+YFLT9efcR4ZNzJlVaK7ndjm4xW9FB/d2u7qYljoa9/erRz3d3L39LCalTxFja4YXtziu/sm07XI9I0IIER+4HFg/1AGAJpELvrVe+tg6KlOSVH3Psz3iUr/dpkO5NTD8e9MlLsdVGc/tvGeZ9mJu4dArpiu/M81Z0XM7B3eetM9oev1l002M7CpUPW57KEWk6ez2plHjTMcvBoDnjj65047cQHuocik9K0KIR43LgXUaAF4tH7UJAC6/OX24dWAhr8Xd8RU9brqdKwxXw01bFeBx/p473LwlTNdOZYc7BAD+2P3RBgBwLvidK9O5MwNAytO327OjfmmN6YUxpvsGmJarE/fAUqOK6ck2AHD285tvAcBrTdeyOsc2pmkczesy8zfTfgf1qAghRHzyYXv3HvvsBtOGg3jAGPfjkzFSazz770hu6bvMXFpNmRvLWdDz+fKuNl2f1/28WxmJVaSk53ZZZ5iu6Obe7i/m6CrdKu7vWbev6akn3dvPogMsJZPCO1klUQ4sIYSIX1wOrEVXAGAc++Xo6MNMjv58h7jtiySdTIf7mUZ2NF0y2zTDOM/tfBjR1fum6S06slYy5Uq2rzy382Luqw4MMLjBLYYba5jmi4j786YcYfoVAw+i6ACbUMzUf4odmYMrLxs5vlUuqWdFCPGocW0hXAMA+ausTQcAX9RP1wUAOt2K2eEl+920TQvTbaNNr3GvdZtCps7knjvMXKwq+AdXBsLYcXabCgAZiuZMywk7q4HknG6vZ+aAMJjVO85wJeMAq32U+yPujrr6AtNjOwDg6EV7pcbFC8yZsoFbSxoWBYCsDf0/BgC/l8Yyme+AIXpUhBAiPhnF8egS+/8oFhc586xpg02e2wXkMv2EjqxojiPnT5mOGMXx5APP7QswWfumMu4jx2Zu4SsYS+6pHNxCuCqfe7tddIx1/Ns0eE6MhjxP3eOmp6u4t5/J9mnoeFu1wFQOLCGEiB/MgeX4eclHADD0wklW93t+FACcXB+bfVGaWwG/OGp6iw6gn74wzRTLVkMvbtV7m1v3wsaarmYgQM7snts5mLO3fVPTG9zJspkOq4IFY2lHe6q203QhzxtZ2nRibdOATQBQpIzruuR0mG70Mq38jJ4VIcSjxpXEvQInzlyx9n8OAMaMcnV8/h+ZfvCEacQg06vMedUmHTvkZJ47zpxlTdcwMiuMHWyP902d++zILAxVXdUJADpOqsQIryXLTV0rIgcY8lr++TvP4+KOHFesFnKsAwAcnUvHVeXrDIGdWsw0/7d3trvUvEIHAGh07DCrDfZro0dFCCHik1F1TcfONp3m596jn75u2qCz5/aBrPo34QzbvePefuUI0+yxbBEsxEisTcnd220sYFpgrOd2ORmJtbKce7tIOrgmhpoGnPXcvt5J0zPj3dvPyWi6i9+rzYt6RoQQIj5gBNaWxacBINulgXXZG1cAgPrbY9oXzzKJ+v7v3e2TH+m4ytzasz3k9ZJpT+bsDf3VdDV3vuSa6LkduPOkLYtsXU9luoU5hgs3jqUdtRE/14Vt7p/30zGmgRfuPD4Dl2Uy9cq7BwCSjtvZ1v5Taa2eFSHEo8blwOJK92zmFklCT31SeuY//M00qqbpAnagzzJU1SvIc4eZhx3177+wow407c5IKudmAMgw2PVxsowHgNwvbEkOAPtW7WX511Dm0nqPSeULfhN3R/0CVzJOvAoAh6/aK/WL/c0cJ+1YVj0oPwAsHO1qVyEVAISP/fs1AOj2kesd+2bRoyKEEPHJaC4UDGC/G8gFlS9XuPfwZ9uZNpju+X2CXjJ9kbml/oyRfH3VVtNcsZQDz88FjvW/urfb/L5p4Xc8t8s8yLQLt6Kcb+ve/lNGUAU5PLevxy3sp308j2yqQiiEEPHD7STuNUw7uvW3t273yhVeNz30rukZpkppx2qB6bt6tkt8Fpi+/YJp2A7TVSGm2T/z3M6rG9+fubhuMHJqEwMGCgzy3M7BXMCNTpheZk6uHawO/xIdZUFbPbcPehoAVpf/vBIA5KlwmeNbxR/1rAghHjUxHFizpgHAgEYp1lsHNprVByPzmy7gHuu0bg4kF3dsFeRKwR9FTEPp2e/Ban9exe48flxOmilLi/QGgC2rTv4MANGrwmface+c4wAwz3NH69JqXFk4XgIATu6OOAIAzav9XMuOeNpVLeoqAOxP5Wr3zBzTg0sAIPyjs/MAoE7+zTSo+r2sR0UIIeKT0dwqOHCm+/+TMxfhl5vde/wz3MpQ37W1sLHn9y3OrXw7w9zbr+b75azjuV3+Y6abSrq328zcVAUbxfJFmDukMZdmLnIciWZVQtuqAQQc99y8Hs97Otj9vHJgCSFE/HDbgcVI106vAUCr20Wknjll+jcX0s+yGFSDiXHbJd65TXuxOFYYAwLWMOVKjsue2zkZWdWOSeGv0w7awi2ABWPJMezIZ9rouuklBhzs5nlKFIj78wby9U+DAWBfn6gBAJDr54s37IiKv+lZEUI8alw5sEYCQEDLbz8EgIu/jude7giuFHzH5IPp5sbd8eWoZLq6Fzvq8+y4nzb1aei5XcYbpj+dBoDo7pGWlH3KkDH2f78tcZ+3GnNgHRsPADeyXeoDAKOKD+WWxPSv8/umAICbxW4PSKzOcYB/n/cHgOvzXykCAI70o1ktakA3PSpCCBGfjOGWwVGx5LpKQYfRzFTuI8BpJkOvXyHu9y951XRPTff2v7F6YbYNntsVsgECm1e5t9vALX/5Y4nQdVwwbeJteonJ36OCTCfw+wS8FqPhG/w+a0wv0rHXarueESGEiA9uO7Aqm7bPab19OaYsOUDH1XkWt2roKho1xrNd4pxq2uMp01BuEfz9VdPcA2Kxa5qbtvnB9NpfptsY4VUkIG576A1Whb9EB9leJpcvXfrO41z80y6AjriJTBUTfggAtjnXrAMA/zk7JtuRlU7oWRFCPGpcVQg3AkDtMpfGAsCNZpH1rANb0N80bb+4O77c1NX72FFXZcddztQ7veeONrNFWkX/1BUAwmtEnwaADyNOPQ0AxZw5DwLAc1zv9vrqtmHQAwDCrtbggHDCOurXTzQCgO1Pv1rXjvOZ4vnzPsvcWQcZSnvuOQ5InewIJ3OofMAy7QMm6lERQoj4ZAz73dPMOVV1iefjkt80nbHbvWc/U960/vdxn6ckI6F2ZnBv/1u4ac5XPLcr+JzphoXu7TayvHihu1QfbMwqiRe5hT2a1W0ndjANjJGby+s9U1UhFEKI+MXlwPrhRQB4qcGsqdbL7z9MO4EOpIZTPNszLvVhbqyebWkPcefKakZO5XzRczsnHVVts5hed5huTmtaaHws56Ud9QZ3vlxiqpfdDDgo9Ufc9lvgMdNPS5lGfWr65WUAePmJktvsyPWXTCuX0bMihHjUuBxYWwDgtTesS7vy9XdvWQeWdkrcHV8udny/s8xr6ArTHuxInVs9d7hPcAXilzymYV4AsPmNabsBwK/s2g/tyNxBAJAspevMgd8CwFeFO54DgPDMZ3i+YwdMq90llPeZD00PMuLrHLcyNtwPAO8WuG04cKVeDiwhhHg4fNjevcc+usj0+Yyej0/OyKgZe9zbnWF12XosDoL3PLd/qrDpzj/d269mhFTOdp7b5aejbVM293ab+5oWjCVHFugwa8KR5iKXZKIZcTaBda4C6Khy0lEnB5YQQsQvLgfWogMAMPaEK8n56Wu0E4bGbV94TzLtxQX7MOb6XcNcwTn2e27n5BbBdplNr88y3ZLUtGBaz+0cA03fOGt6iTtldrO4VolRcX/eAH6/iUtNI7nlcBaLZgXTDsrxhOlG2kGVS+pZEfFHDmalfm4r50Nc4ENLXRtxJ64thFMA4MnMm38EgIUlstQFgHe/i63jy84VhVVhpre4x/ttVt3wKea5Xfoepj89YxpBx9PgUAB4MjBvCk7Yh5vmcpU/nw0A3d+f3BQAzje8wZWCE9yDXv2vuDvqMty7vp8rIudnm744jcedB4AGQ12XxWuq6Whu/RgwTY+KEELEJ6OYJP3G2+49+GFuqavc23O7FKziN8vbvd2pgab1CsV93lLcMrLnpHv737hFMEtxz+0KcyV6y3n3dhuqm+aLpXqtM8C0KSO3LjEXViQNhvEdTQPp0FtFh5kcWEIIET+YA8vRdUkjABj540Vf2glccMdRz/aFk3bPW9xaGFrN9A8Wk8odcefxLv5p32qZ6bUlpttOmhZpGXe717iT5CIjrPZuNC39ftz2UEBF008WmkYuMv0qjWnKHwGg/O1avzmZs1EOLPEwqLzA9G8m+/mGlnkDplBIXcU18dO1erxxJXFnjo1ZFdmRv+S5w8zZynRVC3bUTA7YkxFNPunj7jjf5AqGa2VjKj3+3kwKmIVlzv/4xbTbfAB4Y+S6MgBw9JKr3THm0qq2O+4O3re+6RfT3M/b6oinz9eK8VafP5X0NQDItmJGHvtP/456VIQQIj55P9L0Ezp0vn7BvYee8xTHp7qe2yfnlrwvY1QPPM0JT735cZ+/JB1VOyu7t+/yetztCjLCa0NL93brGTlVIJacXuBKY6NvTC8yB9eVdKYVj5qufNNUDiwhhIgf/LcDgCP94p8AoGC9IQ3Ym/8GAK0+iGlfeA82fetj2kPcUbKKOadyTo3bPknHpOxbM5peXGxaYmPc7VLQYbaKAQE3b5lWqhp3u0BuXZxAuymKOb1mcudKivc9tQdyZjXd6GtaqZKeFRF/VGLgykV/9yf4OudZaxqa9mCV6mIbTYMTXG42b93MhwINB8dW/t0DAL7/GszJkYs5RaZyhbg4t04M+NJ0VKBp1MkY7zvbpAFdQyOSmR68ZfrFxwDQv3pEL56fOais+NILTYccBYARS/1/BwBn18jvAKDLiNVPAcDY3hlC7PiWnOgjHABSImI0AKwLKjUEAHJeadITANYVubAEADpFpW8IAOt3tGC1J2dKAAhKHjEWAFoULtIZAGqUqd8aAD7Gz9z60bIFP+d8PTJCCPFAhh+u7BYtanptr+neUxxX6NBqwNxUN9KYrizFaQInNFHvmu7khOZYUtNMa00ncEtiAc4rTnRl+y84fOTmeekwys/PN4DjYIpPTQ+NZDtWZYpgsvU1P5jma2pa4nnTGf043nFhJpQLI04WGQGTwx9nDq8CHEcnTud4yO9Z7pJp2O88fw09O0II8b8QzVyLfpZL92TmcQCw49tLdQHA4WjJhX1cBoBZb0RMAIDVa4sOAIAitVrNAICtGW6dAIC3Lv6eBQCWhZannVCJW9Wj9wNA+bnRCwHgx/cb9gYA73VFjgHAt7l3jgOAts1L5QGAi32L0pDHbgAo8nLkJwDwW6dakwEgsN/TbQBgUYsDvQGgbdkiRwHgRI4c6dluBABk/y2iBgCs7fpcAACk9nt1BgAsrXb8FgB0fGF7XQD4q2Q92jderEof8QkAFK+aPhwAjvySIR8AnI2qyeq42Vkt2NlBz5B4APPAjzn/GWLqG6NoTdIepmX5dxnOD8/0NN3Mqs8rWUV0Hedjm5gr9UrWh/2N5MB6yD05DYGSAFCzQXY6nCb7cELOKhQDOcH+gCvHUXk9v119TtAnsCO/RoOhOXOWrJ4IAN+M4eF5z9QAgHYD/9wBAD3m5PsYANKnPusNAF3X/HYCAD7qHcaJe0UaNuah9fvBJwMA9MlYeRoAZEub6iUAOLl8xw4AaLNpe2UA2OqTg5FeOZMCgNdBZ10A6F6lQnUAKFAg3RcA4DXU5cDLyu9fnsl/8YceFSGEeCDDDnMeZuRIEJ7T9DpDxJeUM63A15syWXptLlwsY/L1UNfKMHUFHV65OD6VpOPnHSZj31vNdD0dYM45pjcYAbWauT/KMkK4L6sfHuJEas3Xpg463KK4UPP9WxwvuWZflMd9OMx0bUHTgytpMHACtp0LRE46uvI1dr9OORmRZoYM4NQ4JIQQ/xs7aG4ygjYFcyfmYg7CiHEA4HXVazYAvHW1QgcAyDcuw0AAuJBhT1oA6Oi95VUAWB2avjvHs+Tsr3MCQJYLqVMAQO+RFWYCQNINPjcB4LdRK54DgNYDj04CgMuNS5XmeDIFANINCXkSAPqPrzQDAELm+GcDgE2ZVl8BgJaOQ2cB4EyOJ5n0BU0AIPn7yYIBYMDXlVsCQJqxtiV9T5N1jQCgZbL9vQHgcMn8rI7rYHGSyL4AUNaRoykAdDtUajYA9Gts1RbPTi1UwY5L9SrbFdcjJB7APJAR6elZts07edzHOzhPS0sHaik6mqM4b4qmv+EMA3B2PvRvJAfWw8HJjmsnAKx7L1cJ+3s8c1GV5Ir2gI9MRzPnR0ReAHA4+LjcDj2tT8/8p5yYX+eKeTOWFf9ttsd20U9zxbqceWAHnQoAgLCZrdMAwEcVf6Fn1XmWD+oyAGgT5fsTAAxp0ncHACT/KcU6AHAWXOYDANdWtB0IAFt9jjMk0bEDAKq84cwDADOf6HoAAFJ91LAFAES9sfE5APit5wl+3u2hpiO4kuJYqUdGCCEeBBF06LxPR88VJmF/lxN5HyZFbzqW/TAnzCF9TM9y4tOf48dN5phyhPA4LkCMZjmQVzjipMlq+jPP/zUn5E5GbgVwoWUYc129yWpMqf1MVzMZ+3RuEXQlXY/mymF+5rj6nIbBk/y8wdyK8XERji/12J4LQtlpiHxKB1o5bkX8crHpZwyd931Gz44QQvxPhjMjnZIyAnbmRY4LzwBA7WoTXwSALwp2mwAAKdO++DYAONf90RgAQuu0GgwAq48eYASvk/1yZG4AWJgy3QYAqPjGhF8BIKBSEluw2fFRWQDYtfndXgBwefbNprRrZgDAl3NSnACA2sk+OgEAQakDouz1SQ0B4LBPn8EAcOaNaz/ze6wBgI+cAVsAoFHT9/sAQPCWYIsUfmZOJgA4vavbjwBw2P8iAxGcHN8iygLApuuFKgNAoZyT6gDA4c7e0wHgZvVzdvxUrjdhOe0jr1f0DIkHMA/k76cK51VPneY8J8Zxlzif2shiOz8zsGZ5PtMDTAFx6RIbtNO1/W/CHFhD9wJAjoUrswHANf8NMwEgOjp0uGmveqZe2U3ZXf6j5UzrMzfWuS6mB781fXZyLO2oVSuYHr9ieuJ70xqp4m7nx6SDg7abhjMX11Iv08wVPbfzZtWN7qP5PaebrjkDAKem5B1Ag4oGx4DmelSEECI+Gcuk6zMyceLydIwD6OBqt8z0KlfKI+nQeZ9bvf27eX7/FDQsZpdyHxFOZzWt09Bzu0AunHzGiVAUk81f4MLO61zwccQyUSrIiK9N49zPu4EGS4EXPLfLye+7j+NRm6J6RoQQIj7wZ1XbReMBIOmBLmest+5WlnbCD6a2RTA6OudHcdsnaVjkaj6Pi/Q2/YD2iv+nntsl9zGdPtI0yrY2Rn9RxjTwmud2AcNMPynJdj1NvzptmsKtnYt/2hf4zXQjj7v6HAD8umiMBS78vomR0ZXW6VkR8UdlpnA4zxRFx1nEZwIXJqtynhXyvq7V4w0dWEPmAsDLfpEZAOBG0dDq1oH1/szUu0ncHXW9PqZnppgeesn02U/ibvccqwke/dn05LumNbLG3c5V5XAgB5RwlpldzuofWT733M5JB1l3lqcNfdN0XTLT3Kw26Gxg+gENjgET9agIIUR88iFzSUVzJa4/k7Z7f+9+nJNVaTrQ4LjOLYUR3Ho3ko4tv72ez5OSOU9mf+U+QpxkZHGddJ7bBXIF/HNGLEdPMb3Aif3r3KKIkZ7bF85ounWs+3k3cEt8vmLux3txq/tqRmApibsQQsQPVoUQ+KEvAPTMuIlVCK+PMl37K+2EM3HbJ2mWm85fYxpZxXQMk6z7H/PcLoRV2acdZ7u1pl+EmgYt8twu2XjT8UFs19t0zjemqUbG/XnzzTLdSLvvJpO7t78BAMnP5eHOlY2McFYVQhGflODOrRGMbC/JCEHfjLo24k5cEVgVAKDhT2HvAMC5En3oGPLNcGdH5+IOxxVXJs4MMj18ybSSb9wd5nNbTI8WNz3FdjXnxd3OJ6vpgJSmtz4wXXbZNOt0z+2cwaZd6TC72dZ0HSPK8qUGgNVBtw0HGgqjubd2wBQ9KkIIEZ982N69577JnFB9PzP17ud+vNcu044MHb/OJO3h3HI4gpFZSWJJRZCSDquv3nY/70lG3tZJ6rldELeYfDHMvd358aav/Rj39yxKB9e29O7t1zIHVx7mdnQwubscWEIIEb+4HFiLrgPAx8+7qpb/EWGad0/c9lDqZ0y/4Q6SKC6kj+XOkGQ9PdsnwYx4mvY12+01ndyer7/nuV2yLqYff2UayQX6uf1MUy2O+/Pm5QL+egYA3GT1xE6MHHNwXMzxp+lGLuTLgSXiE19uCfTNrmsh4oIOrOFJASDd5h/DAaDG34HDAKDE1ZgdXuY5ppPoADobafr3atOKtT13tC59nhFZx3aanmSEV82+cXe0z+QzXcKVjNCbpkt/N32ikOfzeTOyqxtDY0NrcEDKYZrbbUtkQ2Y8GdEnyc8AkP74ZDquBpTQoyKEEPHJ+4xA2r/Z9Dqr3d6qbtr/V05slri3c9Q3bc+qNNe4UneL1Wd+YfnlKus9nzcFHV1flXQfQU5zfKzTxnO7IFZj+pzVCKOZG+s8c4t8yZD3nDk9ty/M3Fybn3U/7/rOfH2e6aq0pnJgCSFE/GAOLGeZJRcAoNdrmxkplXskANzKFNO+cDLCqdl50w0BpuEDTMdwS1/SNp7tk+TpTKczBUvUJNpXC0yDvvbcLoB2zXhGgkUtMZ3Dv1N+47mdPx1w/bea/kX75zojsDpOMHXMAIDgy67rktMVgcUiJpWf1rMixMMjF3U1lRNb7HpA6tqqwOTliSb5vMuBVYEfn1sv/J8FgG9+d3V8GfeZLnnHNPKC6b78phU6eu4wXVqFHflRhtaepCOr5l1yXJXhSsC+1qYRlUyXrjDNMsxzOy+udHTjHvPQLqbruHc99w3P7RynAOD6+lbFAaDt9CtL7ZW+u/QTEkKI+GQ0HTrvcqW35xXTMBbPuMXqTv32m/rkdm/v4Ba/jqzSd45Jz6NfNz3MrXqV2no+f0r2819Ncx8ZTvJ9a7vmEZnc2wXVNP0szPTaWvf2K1aYZpvk+byFOV/Y2se93Wbm6jrG8tKtX9UzIoQQ8QEjsIov7gMA/vW7fXHnq0cP37YTuIDeijtUrjxheuOi6Zhapklj2YkSMsZ0GnNaRa4znUQHUrCX53bJGDjwMe2aSEZ4zWW7VFM9t/OjnTT8hOmtS6YXuVWxM7cMOkve2S4lS57Uz1HiLQBI8c5+Rg5Xevsh3RCXHT2B+tcDtttduoWaX78BER8P8P9ICpa9fo0roilcD65rRZYTVrDqUTTLioIrwfCiRlC5Uuukw8fJrQCHWLZxGqsPRY1PJNeZ38vJnCNRKwGgfpl0XGF2VT16vpDpKK4Ef8SkuSeuen7bKtyCN/lzU19epze5t/WHc57bleF9n8rytqloGDTjSvYiOiAv9IlhwPA8HZljZCjLq//JZLiu8ut7Y9ka0vw0APgeeT8YAJK0DGQZdOzRT1EIIeKVFhyOuCVwNFd6fVndb8Ax036MqIpicvSRc00jmZvqU05ElzFy9j2WV67B959ER1mLcqbLmXvrPNt1oAPMwRxUL3HCPvEtzg8qmH6fx/TKD6bdOK7MKcp591TTZxlhNYlbG5txHnKY4+qfHH+b8XtMS2H65Fr3yxM1W4+IEELEI5ssIOHmpojdAOBLa/Eka8Rm6tqC/ftoLqAcpgOkDfvvzYwAvnHL/Y2DWb5vzB+mjV5if/8bxw+OR5cj3dslvW46gvZTWxYjmZ/MtD1zP5772L2d7yGOmzzfW4wQ/uWSaQ8u6Owr7Gl8OXcuS1cA2Ob1XnEA2P1dmjoAcAHRPz+cG+HgglUdVgkunIMv0I677YCifQjae9HjY/EfcCHKye/taGp6i9WFJ9PePOXUj0AkIILZQbzOCXBnTijX0AHl8jnvZAfwFrPgt2WW+3asQteOjqxOdHB9zvLXp7OYbmGuilTPJ5IL44rAYoc4cygALOuX82W7Kov8TCMymg5fb5pk9J1XzcUdOa5WmR5lTqpT7UxrHYm7XTluJdw3xvRCI9NXq7sfF1O9mUurC0N5b75ouja3ad6knts5/jJ9M4nplQ0AEJ58+2QAKFNtCYesfrf0GxJCiPhkNKvxDdzp/n8fVgfsU8E0lONrWF7Tvn9zwl7H8/um59r5D93cR4DDzHpYeY3ndinZ788p7d7u5AnT2kPj/j5lG5vu+cK9/TIu3GSLJVl80cycT9R3b9ey3KO9P0W4pbEeDQefpHpmhRD/DVw5sBaz6mwnV2qSt2gn0K64+oLpNkZOPdk7brsmhDmqph40jZpqOoVVCYNne7ZPApqaftyf7Z41nccILfccV/9oEp5nOO2oiK6mi1klPsPKuD9vtvdNly8GgH2zoucBQLZlZ+ngqnjt4dwPR3nTKtVofzPQZArt1TDufLrC3JjvM3ChbUAMu92lXNAa1tJ0czPTm3R8VV2o34BIRHRd4f7TncsV2KT3uMfXyciehpygruaENF+HRHIB6MAalhQAMh1aWgsArpZZzr3X4dy7PbywqV8soa0urcLqgLerCv5tWqtw3O3KMAn7/uSmF8aavtou7nZeLA/blVsCXdUT1zEUNvdVz+2cdTkgsZzt1WWmf+YAgAvJi7Pqwfs0IAZ8qt+KEELEJ2OLm35BB49XDAePFyOd+rxsGspI3FtDTPs1pWPlE8/v73JkLfrTfUQ4zHaVCntul5Lt5sx2b3eKDrBaXElHXs/tS3NCvWe1e/tfGQGebYbndoVPmx7ixLr1jkd0Y6yMOsY5THe+ZpqlkJ5ZIcR/g9tJ3E8CgE+XTnQQtRhKO2E37YTGpkWaxW2fhDCp+1TmtoqknTN5uGlwOc/tkjU3/ZgpTyK58D9voGmqmZ7b+fG44YNMI7igv4RbCDP8eufxLv5pnz3QdAXtppsjAGDNd1PeBADnR5u5gFE5/aO9T8U3mJ5gzsujHBfzr76/98nNnUNLOX63WKvfgEhENBvq/hOeyxXFpAXv7338uLVtBLcaPN01kVwAOrCGLASAF/1v5QWAGw3CN9pVGcm93kl+jrujrsQqHYdbmp5m+dfa/jwun+d2paaZ/lXB9OJPpq9xYHC08tzOQQdXZx53g2VmN7ADz18ils/LahvNWV3kCnN77fzItBgNGWcj0w/GmQ6YqN+KEELEJ2O4xeEWk8V24Mqoo4z7cT78f39uNby1wtRyHQK9WR3Qq7/n82TiQtXiMe4jxEFuSa+wMpbhkhFQc0e4tzt+0LTmt3F/v7KM9Np3zL39Mub0eiJGknYv5hx51FUIszCVwFbOc8L5uV4bqWdWCPHfwOXA+qEpAHTqvrK89dKXmCNqJ3NbFZsetz0UXNN0CnNgRTFZ+9TXTUMyem6XlJFPY1nFMPJ50695vjSxtPN9xXQII7nCaVf9yJzDmXPF/XmzMJn78qVs/51pv7oAkOtUfgZobOT4/KirEOZmca1DzI11rKhp/un/7v1e43jcta9+AyIR0byq+095LldKk/7LlcU3mOS1xm+J5ALQgTU0MwC8tCj8JwC4MGrkl3ZV/OvdeXVc3BFxtcD0CENsT9NxVOcukVplOTD8RU//RR/T16yMa1Rs7bwamHaubHqDHfY6htjm2xXLeS+ZNqdj7grLxdpWwejo4l0AYG1x1zf0ZvpCObCEEOLh8GF79577WnLTDkze6tzufrwvHVv9uOUulCNHKLcS9mZuKZ+PPZ8vwyXTRTGSp//NiXGlbzy3S8UcjnOru7c7cdS01oC4v2c5tt9bxb39UkZ4ZeFE3Mky0o/agVWbOUdO7Xf/vJN9Tb0a6dkVQiRubkdgXQeAcWeiudD9Z1HaCQXjtodCmHJlChfEo+h4msrq5yEtPdsnyWiXjGOqk8g1pl8XMU39ged2SWhHDeXCf/gW0yUjTTMejvvzZg02XUb76BYjv/ozUsyLubdyzDLdSDvoUTuw8nAL4SHer2MMOMk/49+9X2FWO+7BFEFO5cBK1CTjjqkAzp8QSwogXy4oBjLC3ZE2kX3RB+3AqsAJ94ulEskFcG0hvAkAT2xcdhIA2owMSQ0A1SrE7PBcWpnlWI8yF9Ypau2Gno93aT5u8dsbZXoxxPS12nG38y5g2oVbBm9yYFk32DRPjbjbv8gQ2quuAemkaVG3LYqBzVyXxYtJgkfRcBjwpToFIYSIT95nkvRjP5qGnjW9wf63WRfP7by4gtqX41koI6FuMWKqLyOmfA55bs9aJLdzn7hGhEOMrHryC8/tUtGBM3epe7uT/Nx3i8gqU8l0zzX39r8k4fDMSLFVvC4P24Hlza2Dg2i4fEjH3rVFprtoOOQ8oGdXCJG4MQeWM+2SXQDQv/0e2hdFJwBA2L7Y7Ivgo6ZTuKMjqhr/Zi7g4D53Hu8iOtrsr+jo97hzJOpP0697mKZaG0s7ar9ZppE8z4+M+MpwKG57KBu3Rv76m+mt0nw/5gT2rgcA6bK4zpiTSc83ciHov+bASsncWT05v/BSZHGiIvlU+l1YPK8TFwD7FzFtyoAif24xLc7XOzxp+lagaVemnkibKZF88QftwArkXtpUWxPJBXAlca9gOpM5OfzLAUCWJ2N2fBW5MnA4jelp5riqvY3HVfLcYXqx7Os7LAPrWtkY1SfujtbBMq+ddpreZKTWRrbLXzLuDj4V93T/wL3o0Uy6WG+pp/N9xJ3dG55M+RsAlF62hB1jf0VgCSFEvDKKOSiGnzcdOc29h14+3zR9iOf2PnQ09Wf14FtMNh76nenbXFH2Gua5fabapr90dj/vKG5F9H0ilmGUDqx5MdodZxGYGndJ9l6OK7+H2C60jWljOpBWv2HadsLDvR85WbVxIiPgitekIdPTNIIGTeMmenaFEIkbRmC1XDwJAJIv7jnzzleL1olpXwTTATXpjGkkF/anMXVKSLe47Zsi3OJ3sKrpCUZw5Xgz7na5aW/teNr0IneeFF4cd7ss6U2XM6l7OO2ygZ+b+ly+83gv1sL1apJzPQA4Cm6iY6dy/kd7nx60A8vJBbBMrhyWl/RbSAwEcr7WkZHwtZmjMyUX+gowEn4lU06MYLVJ107R7N+bdqVr+hR3XLUISyQX4EE7sBIdLgcWk9bPZq4L/yEAsHy96+pUZtWNw3QEnWFS9jqhcXeYPtz73Yc5tMIYUruig2n2oZ7bOVn9o+N+0xtlTTdwpaDAq3ce7+L/O67mBplGDTIdP8U0IJaOPjgYAMK+n3QLALq8H8ky7X2bq7MQQoj4ZDQdN/25wpuHE9RDP7j31N/7mKZb5/l9fPn/AZzQhHFCcnOx6duMzPJ+x3P7V2mq3LhheouRWEPnmPrFkiw2dW/TeX3dP+9xOtJqPhP39x9zwr3decYE3+AWx1ZZH+79eHWyae8PONHnFs3R5d0/5yw6sHxe0zMshEic3K5CyOryHd0sjNDbG8aDubVvMnNdRX1mOo3J2kMuxG2fFG1vui2T6VVGULXMTPunnmf7JDcjtP5gCpWbjNzq9pyp9zjP7Z6g3bOU1QnDuRVy4A5Tv3yePq+Lca1f/BgAMnc6SUdRpeBHe58etANLJE5eP2T6ykZTx1j315NyvvQzI9iPciHymVXux02jg+swF07rVEskF0AOLJPbDiwmsU/C6hiVFpgeZe6p091N63T23FG61JfJ1/txZSKcOad+Zcee5bTndt7sSDt3ZAfNlYX1dETliY77vKkYAjuXKxSRYaYfs4NP9pzndsF9OSB9AQDhxyN+B4A2By6zCmHfl9RZCCFEfDKa5a0H1XT/f2U6hP4+7t5zf88VtAyxTDi8uSWwX0aaIL+ahtEB1acBHS+X3ds56MDqxCpH15lLIZIOmuF02CSJbWthWdN5X7l/3pNTTWsm89wu6A/TyTxv5G739i3LPZz7kITX872KpuVOu7/+PNfmL31mepD3Je94PcNCiMSJPw3ixZXZ/1cAgEZZbtsJH9NOOGV6O8cVF+iDv7yzt3bxj51RjFv8djA5+9Wepm/WiduuyVfKdB1z/t7kgn4XVk93bPDcLus50+WDTG/Rwdb/R9pbo+M+b82ZALDnrVMzASDfO1dpB1Wc82jv08NyYAVwvH8yt2n1ZabZPtFv5VGSvJtpd9rnaV6O5Tjer+2vmC5jpHvSGA6s7Hy9OCPx/BKLvS8HFifk/IHO9AOAQxlfYM6pQwxlPc2OsA6TBSLYc4fn/bJp3xTsMJmbamVr0+ylPbdz8jydmOTwJttvnGqaP1ncHW0qViOcM980cpTpeIbmJpvluV3QPNMv9ppGvA0AoX99+ToApJ0/hYbG0M9NA1mdKmjrXZTJ4WLLuXI3fEfxfZj7JWj7Xc53kR3uTRpgz97nCVlWPhkjFoKu3tv3DObQnjTbv/ueTl7PQBpKQSfu8bxcAfLt8e/O68fQ0+D093g/+bkC/ubnrnSfJ6QBnZQGaXC6ezwvI0CSMoky3rm/0zoYShvI6qhBR+9yPpYjDqKh79fyXz6/o/k+h+/xe17i5yxs6gpev2e4lcmfW8eCk93jefl7ScaiDY77zO3gYFLrAD6HQZH3eF7+3nxz/svrSwdD0NgY9y228/G+By7g9a13nydkZJQ/I2+Crt/j75S5nJKx/3RUv4sDi8nN37kYw7Dg96vPTe1nN8VwZHECkpnPbSAdV75ccfNmiPiAYqa36HgKY1Wlt7k1z2uhu+MriCHlgzihiWQkVDiTrQ9hRFIQt9oFcSLl5P1JzVxRX7d0/7yurYUvxKh65NvB1LWV4Zv07u06FeX1mO3+/AawKhNu3ufzy1yPgc9TmQzsSSY9ncR+LkdpnoefLysNmA2MLIvIY9rs4L2d17sCP//v9/h7cfUjIWzf5V/2+9wKGZz5PvuHurxer9zn9WXkYACvoytw4a7nZQ61JK7vme/+zuvNeWxQVupfdznfXirHRZ+e93lhaUgkYdGEoCv32D8wh1wyzm8cbe/z+nJ+lGyN+3h5r/1SEpdj9o37O68Xc6oEst8IOn2P95UGlc+P/+75TULHShB/Z0Hb7nI+3ocAzjsc9zvenOT1zXWf95XPnf/1fzkvpGPClZLlrvNCzrODaMd4L/p35w3huPwzf389/QFg1bspab988SkARLxj/fGZ9749BACTnsnY6M553tPBMe2MIjlN/2Tk1DWLsApvNRoACgX7u54fjksBWQDgQrl8jJxaWw4AwkKiTgDAkYKDnwSA2tOTrr7zPFGzXOfLUoGOK543nLmJBx4BgF8L+jXm9X3WfX6QrBUARL5Xg8nfT2QCgO21w/cCQPonjjMiuNLkx8OBlYXzjne4lfTUEtMuP8uJ9ChJSjs1Mxf6nL09H1eKZdlO0x4Y4rpv5x/WJ/XWzXoYRD0NACVyFvoaAFKX+OwMJ9RcSW7FAeK7qrEM6Ex23qMrDQXq2g2c2HIL3sEYye0dvL/t0piOYM6PXZwQNaVBsus6ADgc7gMDkIId/8dBpi+Gmn7a1LQXX7/u737eQDpoPhjDz1cUAE4/v7gFAFQq8XtlADg97E0mdyv/rmlZDuRocJcLSofBaE44frjP+1FnJ+1GV/UtOjywL5YBv6HpCRokbTkgXbrXDoH3711e2eKsBhZV+i4/T3rCV9CA73+f3zMbJ3QjLtDgYyho9LG7TCB/MR1LB80393nel9mxteSEO+L3u0yU6ag7QYOzFwfOo/c60Wa1scH8HZTKdG/ndc41XcOOdxAnHrfu8bypOcGbyAl3ClbZiPaLpQGvSzQnjNMYWTHtfp9fOmDbn+P7Zb7L9+RK12kaxgP4PfdcuEfDlAZJT0bOVO5Hh8Pdri8Nic2MyBlEh/7lVPd23rR87kfScMpKQzDqLucFHRNj6Ij4rsT9Xd8ajOxpTweWl+uBWBvL88fPd5kRM33osNt5j+fz4e+EKZ5QnY6cqDx3+Z2yH97BHFR96Ni8tCSW/oTPZyl+nyfZ//fg85ueET7XGQEVREO9Bu//r3xeznFc+YC5ruay/x7J4x10lPZj/zaQ7x/OHB9/MsdHN54ncCrfl5FhqehI6U6HTl2OB6dpyHSiI2gnk9+2czlGaRg1YEj7Z7xfrVawP6RD/w0aqBljTLTa8H7X5XXwYgd0hJFobfeYXrvH+5qC92ciHWap+DlT0xAM4oplTl63KF7nKP6+fekI9eJ5a7D6z7x2/ByxrFQX43gz0lUl6C7ly12O4usux2FH0z/us19qzPaNG91j/0uH5rYP+Zzwfl386h77X96nEZxg5/z13volBz/fN3SMjmeOtah7/J5F+fkG05ROynlXdCy53xyD6NB1VV+mI/6ne52mc77YwVWFszs/b5W79A90nO6uZdqb3/seu32k4nPxLh0z+ZhrLirjXc5Lh9V3fB7G0jCPuEcDuDD7hYF06ITQIRt9l/sazfsxiuP64vt8fptzp8RLHB+j+TuLbaLnMuyOc0t1G573SpV7O18yPv9D+DwV43WKqnyX52Gg6TIOMIPu83vm4ELDcI4vqWhHRB+P5fnlvPHGAtNhjLBYfZ/nbcDfWcmyABCSJVcvALjxZHPOE7JnB4DvSl6uAABtfsnyFACcXTWC7d4eDwDlUx2n/VKYVXKn0dDOznGk20YA+HTD5JwAsP3ySParxRsBQJFmfmcA4Gpkhp8AIKhqxo4AMKzq2Q8AYHjSKlkB4FZjx1wAcDQedMLsogxcSP6Cn6c8F4aH8vkc9isAVNxRgOPhEI6zSVsDQNMGwakAINzviR4AsL9mqqMA8FbG6GcB4OQRL+YUCvsJjwWH6RibSAd7HS5g+22Tz+BRcoN26o10cR+Xn+OdP8eZjVzoQEo5sP5TRJcHgFxpClnk2TOnh3BiTUNuQVHPDiSXg6tHFw7oHHh+p0HWjAbL3wtjTCA48WjHEWYkHQQ7aOA32X3nxOb/nzcV3388LaqGXHuYwCN70SF2LYbjKpie+w+yAMDFhs3PA8APQedfAIC39iQvBABnfmk03I77lntj13AFPYzVC3C3CCcOHAdcBuW393c/1nDgPu5amXSlj4xlRctBh08oLZdrXBlHkns7XxiTI0+axIkzIyCi77Ii6uT9P+dy6M2+v+95som7IePnWkO6y/s4OME+ePHfPe+/0tA6yOcj6i4rdg46WEI5gTrLlZnb5cvu9vuiI3QK798CztCj7nJ9HTSALnDiG+7yEA6+t/NeepMTQXb0vq7r+30s56NhHs2J1tFL/+76ruFzcWJj3IbT7fPS8Anjiu+x1/nClns7Xzgn6F9lNV3K5yo6/V3OyxX0y5ygXx94f9/zIifyY2hAJf3r3n43Dj53B/5lKP5aVtU7zZVQBw3s6Fh+Dw46UsLZXx5h/3qvHoCIuuwXeJ5V2e/xe7I63VVOfK93jPv4KPa7VRlRF0rD+GM+T9d78Tg6grLRYT+C1Wdy0HAAf29/0yELOozC2D8O4/2O4vH9+L5DOF59QsfMUPa/kYzsi+TzUokLJf3o+CnAyLLjT7F/cDms2P4Ms0m2Ye6uaD7XDelA+5y6lcVCctBQPEZD7COOBxu4Qh/e2t2xc5MO6Zt0SOEeqyhd5XlHsn3IPI7LPM/3NBx3juH5OD5HcH6Qi/3KRxx3y/F6ZOfv8c9YzruX/WBfV1XF3Hd5jpiLMoL3cS8NeAy4v9/Nz3QQ7ma/FrXp3vrfK3wurt1ndaJLdBh8xPEqmata5V0ieB38nZxkhFGUqwpy7Xs77346Yga7ihakivv36uD9jeLv7YDP/X3PSDqMFrDf/qP2PfYPnI9d4/zqqitC8917O+9lOnQm0GEW8Mo9npe/m1OcJ0WuvL/ve5Dz3GGckfrw9xY95S4N6Qg5sPDf9fuL6Xjf5ppPuhxZ4bF8T44HoVxYubH5/s4XyoXeSTzP143ucV7ISN2z7f/d9zxBg3MEF1z86MiNjsVx7KCDJpIOvb31/t15V3LBYHdSAGiz/6kjAPBcX69NABCde+Z1APgsy4TfAODsB9F0TIVRz24HgKH7i/J3P52/22xc6OhKD/TnSwCgTVnQQfJ5WgBo3XFXagDoln34XAB4oktAOABc2jdgKQCM9ln6MwDcyubN8eDcKQC4WuUJzgc+Z8Rwec7v3qVdNawuALwzOZz910E+94OGA8DoLBbZ/FS7t9IAwLDMPp0A4Ku3b6wAgL0319HhO4Xzlm2PmV0eyfEwkvcLJSASMN68TyVp157jvOyv4rE04LjnxXlMJOef97vj5aHzqLYQerGjTc5Ih7Sc8Kbhyrbf4Yd0ARhxMPQXAKj5zPnTAHDkmXpukT53VAX0Ng36wLQPc03dGmK6cqppdn/3di51MvdVR1YvvMEcVRvnmOaLJTeWV3PTjJtNv2IS+MgBpuPXmyY7fme7m0Vd7YOaA8DxsM9DAWBO3+jnAKB6VERzAEg2cTtXivouMM3JAcKxW52BEEI8DD6IWaicK9jT6VgJiSUCtQodu4dXu7dfMNU0Kw0wrxjtvWnQDaABeouOkps87i1GdHg9GcNAo+OsKx0hNxihdosO7ne5UhtCA84Rw4FbjOP/fjoQo/h5N3H86cYFkix0tOP7h3P989MQ/okT9o7cEtlkkGlTOgabcGTtQofZ33QYR9Jw7OyrZ1kIkbjwp0N70c8A8HERV7X06cwJnNzh2T7xZoqVYiw69Sf/f+1d09YZaD+V9Nw+J1OmrGF191Dm/O2+hHZTLDmuMtN+WsYcxeEpTQe9YupT7M7jO3G0+61lklcA4N1u1TYDQPd0x5MDQP6xoQMBwLl4BReMmnA8SJkkYd2nh53EPQ0dzpu4QNkrxlaT5IxATs2t135ciIVDv6n4wJHcNAkXEpy93F9Pz3nKZu7o+oWO/4BWnt8v/QrTGpwH+qZIJBfiYTuw0tIj3p2e+ZH01PfkyvBQrji8xxWp7K/E8wVwy4GVbOa3PgCQ95VkNQAge66YHWaJEaYbXjO99ZTpqnXsiGd47mhtr0p0dHsmDby+0nTzFtNCiKUdtQmTDx66ZRrRy3QiO+xAt5xc7Th9LlM+VSYAeKvd5F4AUPPbKAcAJI/ezfvQjxEquVwGTg51DkII8Sh4jxEDa7gF++QK95FgJleak8cSufY8J5DHGFkUzdwHR5ymbV1b4GNE7vgxwnEQI2XCudXrJidIPRjx4IyxF9ybW5m7c8HpJiOrwkI4geJWnRe5MNWeEUq/M5Lpd05wezDCMysjWJwFH9F8iBF9n3Ai15hbJZpxItiMEUnNGGnShCFW0953v0/f0ZAILKRnWgiROLAqhM7yS64AQPMp39OeSREKAFfHx7RL/FuYjmpjem6m6VVWV29TydQrvWe7Jieruq9ZbhrWw7QHq6T7eHlul4nFrX5hjuGIz00Hbzf1awcA1/7Z0FoaAEJq5Z8KABUbzT8OAFlHXvgbAIJqrusKAEnLtmJkdSqeyZFAI40emQOLuUDfYaRoNS6YdeDfg2l5fsjI0AIh+k3FB+WZG282d0TVPeX+egUuPF7j/PHj7+J+vzpckHw9sW0NbdbWvWuYwxwT/g84yVdyhnZOZq6HL5mUM+MI9+NKcoJ7hGW/mzwbzxfA5cBiyO8s157RlgDAeHdER5f6xXR3uOmVpKYL6JjKsdJzR+vF6hntz5pe32q64YZpge88twPP04RVDy+9Z3qc5WvHsKphYFMAOLPf9XVCaAC9mBcA/Pv+MhwAUlb5qxEAeK0d+A0AePfNvZUddB51BkIIkRAYzYifAdy6VodbM083dh8hZtIhlCKWrcZVORFZy5wHYYyIuspccm05DjtjJEnxZej4QG5ZCePWsZvM+dGTW419Kri386JjrDNz+p2s7f55XVuA19Kk6MHIrGzcKuX48dFed39ep4/p+Cvpd3/tSzPH3XmOp2eYNK9EJj3TQojEgT8dIouZ67BjuztfvZHP1asnvWn6Ie2X8EKmu1hlvVV+U0crz/ZN7t2maxi5dfOkaXfuaPGKpdr6E3Rw/ZLZ9NZV00GXTH0ZAOGk6ZbXImOfGVwGAPz+2PkLAAS/sGEGAASUanUKAEKGpvsJAJ5qnVjuU27azYeY0uAo7eX82eLnfC4H1gbOH5bx+XiWV8yP9nMwUwB8QwfbJFfS8Eb6bT0IHLwPnzI1SAR/kZ24Vd6HDsxBDETZzZQp45irGYzkB1PyFGGxnZ6H3NwhCRDukXdy5dObofsDqrt3ET8wZ0caJpN1cELq4BYBFPh3p3+d73+EWqGs5+Nc1Z0acM9xtrXxfF1cDix6jGdzxdrvvF2VEgxd/cvP9BLLtjZmxJMX/4+n3TtaL4a0dnzC9MYPphuXmuaLpYN2MAS3ySHTy/x7Jx1lxYsBQPsQ18cPYUfxku3Vn/f9YDMvtrYGAD/vQecAwGd+rtkAkHK8ow/wT80NIYQQCYXRXOEcFCNJcD2Ow6e2u48YMxm5lGJSjDdiBFUyOlDe54Q+gkn3r7EqYjtGHHnFyPnjTQfaQOb6CGOOwTC+Xy9GevkyUikDI3o7cUK7iZ9rCyOo32KkdVbXmninhHXdi3Ge8SEjv/3L3F/7IBZNWNjG/f50j9AzLYRIHNx2YLEKbacuAPB8NVev5h9q+uFA00hGSn3LBfd0I9n7TfFs3+Sh/sGtf6G1TN/Kburs7rldlvqmS7mDJdzXdMB8APj+V28WccpFB87gSwDgrP/ndgBIumVdKwAIWdy2MQAE/pVq3p3f+ulZCfaGJKX97aqiyQWuCowIPsOtfKcYyVyWC0zOolRXlfIJ/9vHcDmwNjIG70tGYieJkUPUi7nOxo4xXc6FN4eqFj5QXuPzsIwLga0YKdiROa9q0oHlipD7mpGFLZhDswmfi56cZ+YtnsC/cEpGFnVl0tAP+ZPd51o5dG094MR3NiOfhnCLQDP+kPzP3t95naz+MoEX6iBzLaXNkEAujMuBxQdgRn8A+Mvv6SvWQe7sYnp5oWkThsQ6/vbc0boip9qtNb0+1XTLbNNCVWNpR23sirjiCsduRloVp+MwgBP/FytwwnzJdL2tZD/7zssA4Nc8z9MA4JziVUw/diGESAyM5orqOwExXuCWtfqcsJ5xFQdgldkZzF2VvGsshglzU7mqEkbRoXSVuZvauNaeY2wR9OUEdSDH7VssznGdE6BJjABbTQfWer7v26zel5XJuJ0JfILUgZFmHUb9b+/Thbk/okNMf2Qy96Q99GwLIRI2MR1Yrly9fm1NR/1uGvGr6Xd5TdN9Fbddk+Mv0zX1TMPo+OrJiCvv7Z7bZfrU9GdGbN2ydiUG3wSA/O8WoMNmIMfHbSym9OswAHDUenMAAKS5kJb9siN/4rofDi4QNeQOpuFMarOcC1qXOD6fYnL/FYy0HsFx522OQ1nC/7fPcXsLISOV+/T0fJwXF8Y+YkDIci5kyYH1YPFmSoe8jNCvxirS+RkY5E2XM1jNOWs501r0wzzN40PmJZIvnITVH55i7M0zdISUY9LRot+YluCK69NckXyGSdYLuS7MfVYPcfBCj+UWg4NckUyXUPYWu5K4jwaAwn9ufAsArl3YH2Id5iUmSW98mB3rj547Wif3crfj8dfzmG5qb1qobywdPI9rXI3nY66rndcBYN17z3FiXY+Oxu8ZobaO1Y8GMdlsbnZcjjP6cQshRGJkLCOsFjJiOjCf5+MacGHoTHn3EWUGcyMkj6Vaa1IuXH3AHJSRXMG7ctm0DVfqHDEin3Nw4vMrk4C6qsL9yciungtMs7Pa16PKYXXXiR8NnXSscliHjr09XLEc94Fpen5Pr7vMU7wr8Hozkq3pCtNrRXl/OC7X58pnihdoKJ7Usy6ESFi4HFiL3gCANNl6TbfRZVQ/08jkpguYnD3tujtHHxf/2De5uLVw9TDT0J6mPVgEy/uCZ7soc1bTn6YDQORrUWkBYM1Tc4sBQNChocyZvIk5G1dwq2ML2j9pXVvrE3sxDRYvycMk3M9w4erpWrTXQ2i/05FXjgEoz3B8LnXFNDjZ//Yx/l8S9588HycHlvhP8jI9ycf5gL9wl9D6NOyAynGLnDO+yjnSgTVkOgC82si60OtbL5axDrTxJs8drEu9WH2wPbf63ehiuinQNP8Az+0crBbY+EPTy08AwPmfTswBgGabRtDROJ8T5A1M1jaAHvc87FAcqvIghBD/CT7s6z5STKxhGpjX8/H16Wg63cG93QxWpUmRJRZDhRPaUVyhi3jV9Crbd6lk2o4T6FXcsrCOuZ16cCtD1p2uAT5xXN98KU0/Z66OdYwY2xNquqsIrzuT6GfeEPf75eZK90e8/qvpcNzF+cp6Rrb9SsPqCya7f7mPnnUhRMLC5cD6YQUADOp64EsbVcJDTBfsMU1zOm67KPdPpr+fMg190/StfabOQZ7bPfEXAIQX/XkaAOxdbq/0WHs5FwAkSb+DW7xXMEKsBVPdpK6jexefpKlrummQaa9fPR/nxYWasdzKtpy5r+TAEomaEO6dncBy3IsY0VU4jakfHTMhnPA1oWe9xtB4/mCuCKyRAFA36GoKADhxsOlF61Adn3nuaB0tTduVNb3OKhhbmNuq0NOxdPAVASD6WqMlAHC5w4U/AeDrNXZEyfBDdNytpWEwmCvZuZkTzBmmZ0kIIf6LfJjbfcSIpGPkEzpeAtLFaMCk6w0YsXWGVfGiOU58yRDz5LFM8JPR4fJZsPt5r9GBs4GOnreZZD0bQ9Ed+xLn9fXlynEaOuRSc16SMoR/cwthaiZDtQiBON6PVXzTMFI9zWC+H7d0hnCek4YT+XTc+hFUQ8+6ECJhcXsLYVoAGHcymvbMd9wymD5ZLHYNNQdzBq/JahrGLYK9GLHl3fHO4y9x1Lk5MdMUANjY9OenAODdH+yIfEvDSgGAd6oVtAdbcvxLTQcJlM03XnEyR3Z+prDZ7ZqncF6QlDu5HJyHJGfKoelTOX9YZZqBC0GOK7qmIhGTko6Z9syK/wEjs95mCH9PJuGrzZwRvjfi+QO5cmBlBoCQagtbA0CJowGBAJDLK2YH7chn2pZ6jUndt9Q1LbztzuNdHOC3XDatXAMAmJ3kfG8AqPVntBMAkm/axa0MA7mim+dTdiBb9cwIIcTjwHvcGrH5E9Pzuzjy0IH09oi42zdkzsOztdxHri8LcYLJiN0Mz5u29Tf9Y4vpVkZk9RpkmoNFX5xddG+EEOK/jDmwnC8sCQeA9s8uZ7Gq9EsB4Eak+9F3OK6Yq3cVNaymaa8Opj7DAGDXVVdLRxsA8K2dpzwA1K245HcAePLHyAgA8Lm2mtVgWzDiKm1nNlyie/QwycUtg+8yUvlbVjeeucK0CR1Sab6jXc+UNrP4+nw+Mf259THdJl1T8V+AvvcArnRmYPK9QCaxhfMhfQ6XA4tJ62dxy6L/iwDw5gZXB+3Fjrgdq2HcqG66icnZC/QCgEiX22mV6+1D+L71MwFASL1l8wEg1Sf7p9r/B7CMeB5uzXCU0KMhhBCPI6MZmTOAydhbHXM3FXZXMC2xJ+73qc+cGKf2ubdfweSdKznObGT1wreZUysbF1CwRfdCCCEeJxiB5VxsW8LbdnKrbuuo5j6aREfn9jddM9E0dI1pD+5IcbI6oBeruubaBgCOwu9UBQDvOtu+BYCArRt6AUDA6dbnASDZ1XSFASBE8VWPGJc96sXIazCXtoMBF05XFUJGNnul4uv828Etnl5P8biNuqZCPDhcDiwmnZtNx5o/k7I7bMtfdNtfTK8PNt3KlYmC3wJAdDfX2wUuNm2wFwAcrReWBwDvHdu6AYBfW0sW7/NZ7ggASObtPAkAft11I4QQ4vFmdBvTAdxKmIELKms6upsOO2lKFIuxJo4xJmk54fz2Qox26017c6EoO1e2HTl17YUQ4nEmZhXCTl3ufPVWXddokuOa6Rof07AFpj2aAMDve5w37cgcjNgZwFQwWzsDgOP6qgUAkHxup/kAENwgzVEAKFkXgwHAUUh3Qggh7kZMB5YTAHzG+f9sHXLrgqbXWM51a1XTQjQQArhCUZ9JWhfQYFhnsVt/vLsdAJJE589q//fqoksuhBDi/zO6tunA0+7/z82cVmsWuUwI021c6azO2N9W3KSxeqzpWuZ0evsiDQo6xpxTdK2FEEL8w/9zYHUGgNEZXKNONqZIWTnNNPQgAETf6HUBAHK2zstxpz+rpm9h0a7fmESlFf9Oy0hgR21dcyGE+HfQgTXsEAB4j57jBQBnc3Q9ZB30NSZr3/gHAPz2Srn9dnxtlmf9jknu1vP/g5nMNR+TxDr66RILIYS4Ox9s53i0yPPrubLSMAh3mRSmZ1ktcD1zUrgirHKs4TikKkBCCCHi4LYDq75phwk2yuRg8vVVhQEgYkWEPwCs3PLFFwDgFz6IORs3ci/KSlbTbckquelu6NoKIcSDhQ6sIVEAUOXwsZoAcKP5jcMAcKnJgXAAeG3GQK5cf80qGJtY7WlwU9O86U2dDXRJhRBC3D9jZ5ru4NbAAh+YpuW404ZV7TasNN3Oqr69k5jm4BZEXNW1FEIIce+4HFg/fAIAzzg/3AoAUT+ZQ+pAJ1sw6ZnhTCcA8HlhWy87fiUX7luzKl3qVbqWQggRv7gcWFMBoO7n1kHP3h+VDABKpvuriL3+O1e6hzBpXV6WcfWaqEsohBDif+fD9qYxk7avZDWgDSz20Xe6aXbmXHQO0bUTQgjx7/H/23RRVwDotPz6zwAwPHO0LwAUHHvrZQDwnuwaj1q3Nk1X3NTxrq6hEMKFty5BvOKqdpgXAFYMsj+Wb8BCALhyKsk4+88xGgzRLOP6yi122Py/g/UKEaFLKoQQ4t6JouFQ5kX3/+dm8nY/OrbW+vPvdqaNGbHlnMMGLTQOCSGEuA9Wm/g0BgAszOMNALMz+/oBwOWTq38EgFudp7N4yA/M9Xvylxjvs0mXUgjhwqFLEK+w/l+n8aaVMpmeyWgaSQMB3OPtxZwiqEWN1CUUQgjx74lm9ducX5peYwTWKeYSiXJtCeRWda/3+HdlariuoRBCiH/BARNvC5iomfVpAAhMs7Q8AISlmrITAKL+PlkVACLm6YIJIe7O/wE1fprlk7kfeAAAACV0RVh0ZGF0ZTpjcmVhdGUAMjAyMy0wNy0xNVQxNDowNDo1MyswMjowMAV5oUQAAAAldEVYdGRhdGU6bW9kaWZ5ADIwMjMtMDctMTVUMTQ6MDQ6NTMrMDI6MDB0JBn4AAAAFHRFWHRwZGY6VmVyc2lvbgBQREYtMS41IAVcCzkAAAAASUVORK5CYII="
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "%%tikz -l patterns -s 600,320\n",
    "\\input{fig/formulation.tikz};"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## FEniCS Implementation\n",
    "\n",
    "Having clarified the basic formulation of the problem, we can start implementing it in FEniCS.\n",
    "First of all, we must import the ``fenics`` library into the ``python`` interface. Moreover, we import the ``pyplot`` library that is used for graph rendering.\n",
    "Finally, we will also need some features of the ``numpy`` package."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 169,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# Import relevant modules:\n",
    "import fenics as fe # 'import dolfin' does the same\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Problem parameters\n",
    "As a first step in the actual implementation, we need to specify electrical resistivity and heat conductivity of copper, our material of choice.\n",
    "Following table shows also some other materials for comparison.\n",
    "All the values are for the temperature around 20 °C and atmospheric pressure.\n",
    "\n",
    "<div align = \"center\">\n",
    "    \n",
    "| Material | Electrical resistivity $\\rho\\ [\\Omega \\cdot \\mathrm{m}]$ | Thermal conductivity $k\\ [\\mathrm{W} \\cdot \\mathrm{m}^{-1} \\cdot \\mathrm{K}^{-1}]$ |\n",
    "|--------|:-------------:|:-----:|\n",
    "|**Copper**|$1.68\\times10^{−8}$|   $384.1$ |\n",
    "| Silver |    $1.59\\times10^{−8}$   |  $429$   |\n",
    "| Stainless Steel  |    $69.0\\times10^{−8}$   |  $16$&ndash;$20$  |\n",
    "| Glass  |    $10^{11}$&ndash;$10^{15}$   |  $0.8$&ndash;$1.9$  |\n",
    "\n",
    "</div>\n",
    "\n",
    "It is recommended to gather all the problem parameters at the beginning of the solution script file, let us follow it by specifying the geometeric, the environment, and the material parameters."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 170,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# --------------------\n",
    "# Parameters\n",
    "# --------------------\n",
    "# Geometric parameters:\n",
    "L = 5.0  # Length of the wire [m]\n",
    "R = 0.005 # Radius of the wire [m]\n",
    "A = np.pi*R**2 # Cross-section area [m^2]\n",
    "\n",
    "# Environment parameters:\n",
    "T_c = 283 # Ambient temperature of the cellar [K]\n",
    "T_h = 298 # Ambient temperature of the outside air [K]\n",
    "U = 230 # Voltage drop [V]\n",
    "\n",
    "# Material parameters:\n",
    "rho = 1.68e-8 # Electrical resistivity of copper [Ω m]\n",
    "k = 384.1 # Thermal conductivity of copper [W m^-1 K^-1]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Geometry of the problem\n",
    "\n",
    "Next, we define the discretized geometry of our problem. There is always a possibility to use an external mesh defined in an *.xml* file. In our case, we use the predefined interval mesh provided directly by FEniCS environment.\n",
    "Initialization of the FEniCS object ``IntervalMesh(nx, x_1, x_2)`` creates an equidistantly discretized one-dimensional mesh with ``nx`` subintervals between the end points ``x_1`` and ``x_2``.\n",
    "\n",
    ">You can invoke a basic description of any function/object by ``help`` command. Try that for the object ``IntervalMesh``.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 171,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%capture\n",
    "# %%capture supresses the output of this notebook cell, which is quite lenghty for the help command\n",
    "help(fe.IntervalMesh)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 172,
   "metadata": {},
   "outputs": [],
   "source": [
    "# --------------------\n",
    "# Geometry\n",
    "# --------------------\n",
    "nx = 10 # Number of mesh subintervals\n",
    "mesh = fe.IntervalMesh(nx, 0.0, L)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can visually check the mesh using the FEniCS ``plot()`` command. The plot object is then displayed by adding the pyplot command ``show()``."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 173,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAWAAAAAmCAYAAAD6F0GJAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAAEuUlEQVR4nO3cT4jUdRjH8ffHbcTNkjmslO5a2yG8FGQsQhgdhNJKSjoV1CmoQ4ERGApdugVBdI7qEEUWaBEFmaAQRqW7/snUDAkr12KTWEwytPXpsFOo7cKOzvf3bPP7vGDYmR+zv8/zsMzDb7/fmVFEYGZm1ZuTXYCZWV15AJuZJfEANjNL4gFsZpbEA9jMLIkHsJlZkqvaeXJfX18MDg62FTD+xzl+OfUn5ybO0+iZw/UL5tG8utHWOS5HVm5mtnt2z92Ym5ndqdyRkZGTEbHw0uNtDeDBwUGGh4dn/PwP9o6yccsB+s5N/Hus0ejh+YduZe2y/nai25KVm5ntnqvLzcyuW25mdidzJf0w5fF2PogxNDQU7QzgFS9uZ3T8zH+Oz+2Zw7IbmjM+T7v2/jjO2YnzledmZrvn6nIzs+uWm5k9XW5/s5fPN6xs61ySRiJi6NLjRdeAT0wxfIEpm+qk6c5fOjcz2z1Xl5uZXbfczOzpzj/dXLscbS1BtGtxs3fKK+D+Zi/vPnlHsdzprrxL52Zmu+fqcjOz65abmT1d7uJmb8cyil4Br1+1lN5Gz0XHehs9rF+1tGRsWm5mtnuuLjczu265mdlV5Ba9Av5nofqlrUc4MX6Gxc1e1q9aWnzRPis3M9s9u+duzM3MriK3rU04Sb8CU+7mzUAfcPIyf/f/yj3XQ916rlu/cOU93zjV29DaGsBXQtLwVLuA3cw910Pdeq5bv1CuZ38SzswsiQewmVmSKgfwqxVmzRbuuR7q1nPd+oVCPVe2BmxmZhfzEoSZWZLiA1jSaklHJB2VtKF03mwg6Q1JY5K+ya6lCpKWSNoh6ZCkg5LWZddUmqR5knZJ2t/q+YXsmqoiqUfSXkkfZddSBUnHJB2QtE/SzL8MZybnLrkEIakH+A64GzgO7AYeiYhDxUJnAUl3AaeBNyPilux6SpO0CFgUEXskXQuMAGu7+e8sScD8iDgtqQHsBNZFxJfJpRUn6VlgCFgQEWuy6ylN0jFgKCI6/t7n0lfAy4GjEfF9RJwFNgEPFs5MFxGfAb9l11GViPg5Iva07v8OHAbKf0QqUUw63XrYaN26fkNF0gBwP/Badi3doPQA7gd+uuDxcbr8hVl3kgaBZcBXyaUU1/pXfB8wBmyLiK7vGXgFeA4o/zVos0cAn0oakfREJ0/sTTjrGEnXAJuBZyLiVHY9pUXERETcBgwAyyV19XKTpDXAWESMZNdSsTsj4nbgXuCp1hJjR5QewKPAkgseD7SOWZdprYNuBt6OiC3Z9VQpIsaBHcDq5FJKWwE80FoT3QSslPRWbknlRcRo6+cY8D6TS6sdUXoA7wZulnSTpLnAw8CHhTOtYq0NqdeBwxHxcnY9VZC0UFKzdb+XyY3mb1OLKiwiNkbEQEQMMvla3h4RjyaXVZSk+a2NZSTNB+4BOvbupqIDOCL+Ap4GtjK5MfNeRBwsmTkbSHoH+AJYKum4pMezaypsBfAYk1dE+1q3+7KLKmwRsEPS10xeaGyLiFq8LatmrgN2StoP7AI+johPOnVyfxLOzCyJN+HMzJJ4AJuZJfEANjNL4gFsZpbEA9jMLIkHsJlZEg9gM7MkHsBmZkn+Btfj/UVqu0p8AAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fe.plot(mesh)\n",
    "plt.show() # this forces the plot to show in a pop-up window when in command line/script environment"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Weak formulation\n",
    "\n",
    "As was said previously, any sensible nemerical simulation of a problem starts with its mathematical formulation, typically in the form of a PDE.\n",
    "This is sometimes called a classical formulation.\n",
    "Now more specifically, any finite element method computation starts with the so-called weak formulation of the problem, which stems from the PDE governing the problem, but should be understood rather as its generalization.\n",
    "One of the very nice features of FEniCS is that the implementation of the weak form closely imitates its (written) mathematical form.\n",
    "The weak formulation<a id='weak_form'></a> of our problem reads\n",
    "\n",
    "\\begin{align*}\n",
    "\t\\text{look for } T \\in V: a(T,\\vartheta) = \\mathscr{l}(\\vartheta), \\quad \\forall \\vartheta \\in V\\,,\n",
    "\\end{align*}\n",
    "\n",
    "where $T(x)$ represents our weak solution, that is, the temperature field, $\\vartheta(x)$ is the so-called test function, \n",
    "\n",
    "\\begin{align*}\n",
    "\ta(T,\\vartheta) = \\int_0^L T'(x) \\vartheta'(x)\\, \\text{d} x\n",
    "\\end{align*}\n",
    "\n",
    "is the bilinear form, that is a mapping $a: V \\times V \\rightarrow \\mathbb{R}$ linear in both of its arguments, and\n",
    "\n",
    "\\begin{align*}\n",
    "\t\\mathscr{l}(\\vartheta) = \\int_0^L  \\frac{A^2 U^2}{\\rho k} \\vartheta(x)\\, \\text{d} x\n",
    "\\end{align*}\n",
    "\n",
    "is the linear functional defined on $V$, which represents some linear space containing all the possible (weak) solutions to our problem[<sup>1</sup>](#fn1 \"Definition of Sobolev space\").\n",
    "Since the space $V$ might be infinite dimensional, a crucial step of any numerical method, which begins with such formulation, is to approximate $V$ using some smaller, i.e., finite-dimensional, space $V_h$.\n",
    "The hallmark of FEM is that it uses piece-wise polynomial functions for the definition of $V_h$.\n",
    "\n",
    "#### Function spaces definition\n",
    "Thus the first step in our FEniCS implementation of the finite element formulation is to define the proper function space approximation $V_h$ that contains the piece-wise polynomial solution and test functions. \n",
    "This is done by initializing the object ``FunctionSpace(mesh,type,deg)`` with three arguments.  \n",
    "The first argument is the discretized domain contained in the ``mesh`` variable.\n",
    "Second and third are the type ``type`` of the base functions and the order ``deg`` of these functions respectively.\n",
    "We will present various types of function elements throughout this course.\n",
    "You can check an encyclopaedia of finite elements e.g., [here](https://defelement.com/).\n",
    "In this simple setting, we define the conforming space of finite elements, i.e., $V \\subset V_h$, or more specifically the *Lagrange/Continous Galerkin* element (hence the acronym ``CG``) of order one, in other words, piecewise linear functions.\n",
    "\n",
    "<span id='fn1'><sup>1</sup>It is the result of functional analytic approach to the theory of PDEs that such object is the [Sobolev space](https://en.wikipedia.org/wiki/Sobolev_space#One-dimensional_case) $H^1_0 = \\{ v \\in H^1(0,L); v(0) = v(L) = 0\\}$.\n",
    "Further it holds that $H^1(0,L) \\subset AC[0,L]$ (which denotes the space of [absolutely continuous](https://en.wikipedia.org/wiki/Absolute_continuity#Definition) functions on $[0,L]$) meaning that it makes sense to use point-wise of the functions.</span>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 174,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# ---------\n",
    "# FE spaces\n",
    "# ---------\n",
    "V_h = fe.FunctionSpace(mesh, \"CG\", 1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Then, we need to initialize the solution and test functions (specified by the definition of $V_h$) as ``TrialFunction`` and ``TestFunction`` objects, respectively."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 175,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Definition of trial (=solution) and test functions:\n",
    "_T = fe.TrialFunction(V_h)\n",
    "Theta = fe.TestFunction(V_h)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Boundary conditions\n",
    "As this next step, we specify the geometric boundary of the domain through which we apply the Dirichlet boundary conditions.\n",
    "It is necessary to mark those entities of the mesh that constitute the boundaries $\\Gamma_{\\text{D}}$ and $\\Gamma_{\\text{N}}$.\n",
    "This can be done essentially in two ways."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 176,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# ------------------------\n",
    "# Boundary markings\n",
    "# ------------------------\n",
    "def left(x, on_boundary):\n",
    "    return fe.near(x[0], 0.0) and on_boundary\n",
    "# analogously, we could define the marking of the right boundary\n",
    "def right(x, on_boundary):\n",
    "    return fe.near(x[0], L) and on_boundary"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In the previous, ``on_boundary`` flag is ``True`` for all points that lie on the boundary of the mesh and ``near(x,x0)`` checks whether ``x`` is close to&nbsp;``x0`` (by default, within the tolerance given by the global constant ``DOLFIN_EPS``), which respectively gives the left and the right end of the beam.\n",
    "  \n",
    ">When you compare two real numbers in ``python``, always use certain **tolerance**. Value of the tolerance depends on the size of the numbers you compare. For unit sized numbers, ``DOLFIN_EPS = 3e-16`` is sufficient.  \n",
    "\n",
    "Try it for yourselves. Print the outcomes e.g., of ``0.1+0.1+0.1-0.3`` and ``100.1+100.1+100.1-300.3``:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 177,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "5.551115123125783e-17\n",
      "-5.684341886080802e-14\n"
     ]
    }
   ],
   "source": [
    "print(0.1+0.1+0.1-0.3)\n",
    "print(100.1+100.1+100.1-300.3)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "There is an equivalent way of marking the boundary. It employs the ``MeshFunction`` object that represents a numerically/Boolean valued function defined on the mesh. More precisely, the first argument in initialization can either be ``‘int’``, ``‘size_t’``, or ``‘double’``, which denotes the integers, unsigned integers and floats, respectively. The Boolean valued function is initialized with the ``‘bool’`` argument."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 178,
   "metadata": {},
   "outputs": [],
   "source": [
    "boundary = fe.MeshFunction(\"size_t\", mesh, mesh.topology().dim()-1,0)\n",
    "for v in fe.facets(mesh):\n",
    "    if fe.near(v.midpoint()[0], 0.0):\n",
    "        boundary[v] = 1 # left boundary\n",
    "    elif fe.near(v.midpoint()[0], L):\n",
    "        boundary[v] = 2 # right boundary"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The previous ``for`` loop iterates through facets of the mesh (vertices in our example) and assigns the value one to the left boundary and the value two to the right boundary of the bar.\n",
    "\n",
    ">note: it is preferable to set some default numerical value (e.g., zero) to all the other mesh facets, this can be done either by giving the last optional argument ``0`` in the initialization of the ``MeshFunction`` object or by calling ``boundary.set_all(0)`` **before** the ``for`` loop\n",
    "\n",
    "Though from the perspective of physics there are many boundary conditions (depending on what type of problem we solve), mathematically speaking all of the standard boundary conditions fall into one of the three following categories.\n",
    "\n",
    "##### 1. Dirichlet boundary conditions\n",
    "\n",
    "The first type of boundary conditions is the one that prescribes some value of the solution at the boundary.\n",
    "With heat transfer problems, this constitutes the case when the boundary is held at a constant temperature, e.g., by a large reservoir of water.\n",
    "\n",
    "> note: In the realm of variational, or weak formulations of PDEs, the Dirichlet boundary conditions are sometimes called **essential**, since they determine the definition of the solution function space.\n",
    "\n",
    "In FEniCS, Dirichlet boundary conditions are implemented through the object ``DirichletBC`` initialized with three arguments: the function space, the prescribed value at the boundary, and the corresponding subdomain specified by the previously defined ``left`` function (or equivalently by the ``boundary`` variable and its numerical value at the desired facets)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 179,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# --------------------\n",
    "# Boundary conditions\n",
    "# --------------------\n",
    "bc_left = fe.DirichletBC(V_h, T_c, left) # prescribe cellar temperature to the left end of the wire\n",
    "bc_right = fe.DirichletBC(V_h, T_h, right) # prescribe outdoor temperature to the right end of the wire\n",
    "\n",
    "# Group boundary conditions into an array:\n",
    "bcs = [bc_left, bc_right]\n",
    "\n",
    "# Equivalently, using boundary variable:\n",
    "#bc_left = fe.DirichletBC(V_h, T_c, boundary, 1)\n",
    "#bc_right = fe.DirichletBC(V_h, T_h, boundary, 2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "raw_mimetype": "text/restructuredtext"
   },
   "source": [
    "##### 2. Neumann boundary conditions\n",
    "The second type of boundary conditions specifies the value of the solution's normal derivative across the boundary.\n",
    "For our physical problem, this means that the heat flux over the boundary is given by prescribed rate $g$, i.e.,\n",
    "\n",
    "\\begin{align*}\n",
    "\t\\mathbf{q} \\cdot \\mathbf{n} = -k \\frac{\\partial T}{\\partial n} = g\\,,\n",
    "\\end{align*}\n",
    "\n",
    "where $\\mathbf{n}$ denotes the outer unit normal (in our case $-1$ for the boundary $x = 0$ and $+1$ for $x = L$), and $g>0$ some given heat inflow through boundary.\n",
    "\n",
    ">note: The category of Neumann boundary conditions is sometimes denoted as natural. The natural boundary conditions are reflected directly in the weak formulation of the problem.\n",
    " \n",
    "In the [weak formulation](#weak_form) above, the Neumann boundary condition would appear as a term from the integration by parts formula. To include this term in our implementation of the problem, we would first need to define the integrating measures that correspond to the boundary manifold of the mesh. To this end, we could initialize the ``Measure`` object with input arguments: string that specifies the integral type, ``\"dx\"`` for volume, ``\"ds\"`` for surface integral, and the respective geometry to which the measure corresponds, ``mesh`` for volume, ``boundary`` for surface integral."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 180,
   "metadata": {},
   "outputs": [],
   "source": [
    "dx = fe.Measure(\"dx\", mesh)\n",
    "ds = fe.Measure(\"ds\", subdomain_data = boundary)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### 3. Newton-Robin boundary conditions\n",
    "The last class of boundary conditions represents a linear combination of the two previous kinds.\n",
    "Therefore both the solution's value and its normal derivative are prescribed, mathematically\n",
    "\n",
    "\\begin{align*}\n",
    "\taT + b \\frac{\\partial T}{\\partial n} = g\\,,\n",
    "\\end{align*}\n",
    "\n",
    "where $a, b \\in \\mathbb{R}$ are arbitrary constants (weights).\n",
    "As the name might suggest, the convective heat transfer (Netwon's law of cooling) can be described by this type of boundary condition.\n",
    "Namely\n",
    "\n",
    "\\begin{align*}\n",
    "\t\\mathbf{q} \\cdot \\mathbf{n} = -k \\frac{\\partial T}{\\partial n} = r\\,(T_\\infty - T)\\,,\n",
    "\\end{align*}\n",
    "\n",
    "where $r$ denotes the coefficient of heat transfer and $T_\\infty$ is the ambient temperature."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Discretized weak formulation and problem form\n",
    "\n",
    "We are now ready to give the weak form of the problem. One of the very nice features of FEniCS is that the implementation closely imitates the [weak formulation](#weak_form) given above.\n",
    "Weak formulation<a id='weak_form'></a> is\n",
    "\n",
    "\\begin{align*}\n",
    "\tT - T_\\mathrm{b} \\in V_h: \\sum_{i=1}^n\\int_{x_{i-1}}^{x_i} T'(x) \\vartheta'(x)\\, \\text{d} x = \\sum_{i=1}^n\\int_{x_{i-1}}^{x_i} \\frac{A^2 U^2}{\\rho k} \\vartheta(x)\\, \\text{d} x, \\quad \\forall \\vartheta \\in V_h\\,,\n",
    "\\end{align*}\n",
    "\n",
    "where $V_h = \\{ v \\in C([0,L]); v(0) = v(L) = 0, v|_{[x_{i-1}, x_i]} \\in P^1\\}$, and $T_\\mathrm{b}$ is some particular solution that respects the boundary conditions, e.g., linear temperature distribution."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 181,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# --------------------\n",
    "# Weak form\n",
    "# --------------------\n",
    "#a = E*A*fe.inner(fe.grad(T),fe.grad(Theta))*dx\n",
    "# One can equivalently write the bilinear form in this one-dimensional problem as\n",
    "a = _T.dx(0)*Theta.dx(0)*fe.dx\n",
    "l = A**2*U**2/(rho*k)*Theta*dx"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We need to supply the subdomain measure ``ds`` with the value of the ``MeshFunction`` object that corresponds to the Neumann part of the boundary $\\Gamma_{\\text{N}}$, i.e., two.\n",
    "\n",
    ">*note*: it is possible to access individual elements of the gradient ``grad`` of a function ``u`` via ``u.dx(i)``, where ``i`` specifies the derivative with respect to the ``i``-th coordinate and the enumeration starts with **zero**, do not confuse this with the measure object ``dx`` that specifies the integration domain\n",
    "\n",
    ">``dot(u,v)`` vs ``inner(u,v)``: FEniCS makes distinction between the inner product and dot product of two elements ``u``, ``v``, ``dot(u,v)`` contracts (sums) over the last index of the first element and the first index of the second element, whereas ``inner(u,v)`` sums over all indices of the elements that must be of the same order\n",
    "\n",
    "* if ``u``, ``v`` are both vectors (rank 1 tensor), then ``dot(u,v)=inner(u,v)``\n",
    "* matrix-vector multiplication $\\sigma \\mathbf{n}$ is computed using ``dot(sigma,n)``\n",
    "\n",
    "Finally, we call the solver to compute the linear algebraic system generated by the finite element method. FEniCS is prompted to do that by calling the ``solve()`` function supplied with three parameters: linear variational equation, function that will store the resulting approximation and the problem specific Dirichlet boundary conditions.\n",
    "\n",
    "Let us first introduce the function ``u`` that will contain the finite element solution. In the code environment we define the object ``Function`` pertaining to the function space $V_h$. Then we can call ``solve(a==l, u, bc)``."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 182,
   "metadata": {},
   "outputs": [],
   "source": [
    "# --------------------\n",
    "# Solver\n",
    "# --------------------\n",
    "T = fe.Function(V_h)\n",
    "fe.solve(a == l, T, bcs)\n",
    "# Equivalent implementation:\n",
    "# F = a-l\n",
    "# problem = fe.LinearVariationalProblem(fe.lhs(F),fe.rhs(F),u,bc)\n",
    "# solver = fe.LinearVariationalSolver(problem)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Post processing\n",
    "\n",
    "\n",
    "There is more than one way of visualizing results from FEniCS computations, the readily available one is the tools offered by ``matplotlib`` library\n",
    "\n",
    "### Graphing with ``matplotlib``\n",
    "In the final stage of our simulation we visualize the computational results. This task is easily handled by the ``matplotlib`` library that we imported at the beginning and already used for mesh visualization."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 183,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYUAAAEGCAYAAACKB4k+AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAA2fklEQVR4nO3dd3hUZdrH8e+dDimUFEqAFHpIIEAoQggKioCIgoA0EdEFV7Gsri6W13VZdy1rryiKiCAqKsqKDWmhSUlASugkQEILhBrS87x/ZMgGBBJIZk6SuT/XNRczZ2bO+SUkueec55z7EWMMSimlFICL1QGUUkpVHloUlFJKFdOioJRSqpgWBaWUUsW0KCillCrmZnWA8ggICDChoaFWx1BKqSolISHhqDEm8GLPVemiEBoayrp166yOoZRSVYqI7L3Uc3r4SCmlVDEtCkoppYppUVBKKVWsSo8pXExeXh6pqalkZ2dbHUVdgpeXF40aNcLd3d3qKEqpC1S7opCamoqvry+hoaGIiNVx1AWMMRw7dozU1FTCwsKsjqOUukC1O3yUnZ2Nv7+/FoRKSkTw9/fXPTmlKim7FQURmSYiR0Rkc4llX4jIBtstRUQ22JaHikhWieemlHPb5Uyv7En/f5SqvOx5+Gg68DYw49wCY8zt5+6LyCvAyRKv322MibZjHqXszhjDziNniN+RTq0a7rRpWItmQT54uF3B56/lr0NwBwiL+9+y5HhIS4TYhys6slLnsVtRMMbEi0joxZ6Too+Kw4Be9tq+lVxdXYmKiip+/O2335KSksItt9xy3nH0l19+meuvvx4RYdSoUcycOROA/Px8GjRoQJcuXfj+++/LlcXHx4czZ85c9jWvv/4648ePp2bNmgD079+fzz77jNq1a5dr284iv6CQdXuPsyDpML9uPczeY2fPe97dVWge5Eubhn5ENPSjTcNatG7gi6/XJQbagzvAnLEwdHpRYUiO/99jpezMqoHmHsBhY8zOEsvCRGQ9cAp42hiz7GJvFJHxwHiAJk2a2D3o1ahRowYbNmw4b1lKSgo9evS46B95b29vNm/eTFZWFjVq1GDBggUEBwc7KG1RURg9enRxUfjhhx8ctu1K5Qo+oZ/JySd+Rzq/Jh1m0fYjnDibh4erC92a+TM+LpzrWgaRlVfAlgOnSDpwiqSDp1i8/QhzElKL1xHiX5OIBn5ENPCjTbAfEQ1qUc/PEwmLKyoAc8ZCzN2w7qP/FQil7MyqojACmF3i8UGgiTHmmIh0BL4VkTbGmFMXvtEY8wHwAUBMTEy1mTauf//+zJ8/nyFDhjB79mxGjBjBsmV/rItbtmzhrrvuIjc3l8LCQr7++muaN2/Oq6++yrRp0wC45557ePjhh89735IlS3j55ZeLi9LEiROJiYnh1KlTHDhwgOuuu46AgAAWL15c3D4kICDgoutNSUmhX79+xMbGsnLlSoKDg/nuu++oUaOGfb9J9lbKJ/TDp7JZkHSYBUmHWbX7GLkFhdSu6U6vlkHcEFGPHi0C8fE8/1eqaaAPA9s1BIoOLaWfzikqFAeLisWWAyf5cfOh4tf7e3sQ0dCPiIaBDG48jJbxL1HY4zFctCAoB3F4URARN2Aw0PHcMmNMDpBju58gIruBFkC5Ghv9479bSDrwh7pSLhEN/fj7zW0u+5qsrCyio6MBCAsLY+7cuQAsW7aseDnA119/TdOmTQEYPnw4kydPZsCAAWzcuJFx48ZdtChMmTKFhx56iFGjRpGbm0tBQQEJCQl8/PHHrF69GmMMXbp0oWfPnrRv377Ur+fBBx/k1VdfZfHixQQEBJz33KXWW6dOHXbu3Mns2bOZOnUqw4YN4+uvv2b06NGlbq9Su+ATuln3Eft6vcu8PQ1YMH85G1OLhsBC/Gsy5poQro+oR0xIHdxcyzZeICIE+XkR5OfFda2Cipefyclnq61IJB04xZaDJ0laMZ/xrjN5o2AQd8RP4dWNdTGhPYoPP7Ws50sND1c7fBOUs7NiT+F6YJsxpng/WkQCgQxjTIGIhAPNgT0WZKsQFzt8BFzy8BFA27ZtSUlJYfbs2fTv3/+S677mmmv417/+RWpqKoMHD6Z58+YsX76cQYMG4e3tDcDgwYNZtmxZmYrC5VxqvQMHDiQsLKy4wHXs2JGUlJRybauyyGsSy6GmI2gc/xLT3Ybyj68KgB20b1Kbx25sSZ+IejQL8qnQM6h8PN3oFFqXTqF1ixYkx2PmvMu+Xh/SSCL5fntvHt/1FH/5PYtZq1sB4CIQHuhTNE7RoKhQRDT0o663R4XlUs7JbkVBRGYD1wIBIpIK/N0Y8xEwnPMPHQHEAZNFJA8oBO41xmSUN0Npn+grm4EDB/LXv/6VJUuWcOzYsYu+ZuTIkXTp0oX58+fTv39/3n///TKt283NjcLCwuLH5b1OwNPTs/i+q6srWVlZ5VqflU5n57F0RzoLkg5zZtsiXjKf8I4ZzJjCXwjp2ZfI2AEE+Xo5LlBaIjJ0OiFhcYQAdLwDkkP4MC2R1IjrSDp4qnisYl3Kcb7bcKD4rQ1qedG9WQBP9GuFv4/nJTeh1KXY8+yjEZdYPvYiy74GvrZXlqpi3Lhx1K5dm6ioKJYsWXLR1+zZs4fw8HAefPBB9u3bx8aNG4mLi2Ps2LFMmjQJYwxz587l008/Pe99ISEhJCUlkZOTQ1ZWFgsXLiQ2NhYAX19fTp8+/YfDRz169Ch1vVXVgRNZ/Lq1aHzgtz3HyCsw9Km5gzdc32Br7Nvc1X0ANdNW0mvOWGgVBL4OPKZ/sdNOw+KQsDgaA43r1uTGNvWLnzqemctWW6HYfOAk8zYcYMn2I/x7UBR9SrxOqbKodm0uKrMLxxSefvpphgwZUvy4UaNGPPjgg5ddx5dffsmnn36Ku7s79evX58knn6Ru3bqMHTuWzp07A0UDwhceOmrcuDHDhg0jMjKSsLCw854fP348ffv2pWHDhixevLh4eYcOHS663qp4qMgYQ9LBU8UDxVtsY03hAd6M6x7G9RH16Lh/Jy6NZtLp3KDuuTGGtMRKfeZPHW8PujULoFuzoqJ+37WneeTLDYz/NIEhHRvxzM0R+F3q9FelLiDGVN0TeGJiYsyFk+xs3bqV1q1bW5RIlZUj/p8KCg0rdx/l16TD/Lr1CGknshCBjk3qcH1EPW6IqEfTQB+7ZrBKbn4hby3aybtLdlPfz4v/DGlbXDSUEpEEY0zMxZ7TPQVVLR3PzGXi7ERW7DqGl7sLPZoH8tD1zenVKogAJzjW7uHmwqN9WtKrVRCPfvk7Iz9czdhuofytbys9a0ldlhYFVe0kHTjF+E/XceR0Ds/dGsltHRo57R/C9k3qMP/BHrz40zamr0whfkc6r94eTXTj2lZHU5VUteuSqpzbf38/wOD3VpBfYJgz4RpGdw1x2oJwTg0PV54d2IZZ93QhO6+A295bySu/bCc3v7D0Nyuno0VBVQsFhYbnf9zKA7PXExVci/8+EEs7/TR8nu7NAvjpL3HcGh3MW4t2MejdFWw/dNrqWKqS0aKgqrwTZ3MZ+/Ea3l+6hzu6hjDrnq4E+lb/cYOr4eflzivD2vH+HR05dDKbm99azvtLd1NQWHVPOFEVS8cUVJW27dApxs9I4NDJbF4YHMXwzpWzSWJlc2Ob+nQMqcNTczfx/I/b+HXrYV4ZGk0T/5pWR1MW0z2FKubbb78lKSmp3OtZsmQJAwYMuOxrTpw4wbvvvlv8+MCBA+ddV2G1HzYdZPC7K8nOK+DzCV21IFyhAB9PpozuyKvD2rHt4Gn6vhHPZ6v3YYwp6hibHH/+G5Lji5aras25i0IV/MGvqKJQFhcWhYYNG/LVV185ZNuXU1Bo+M/P27hvViKt6vvy/QOxdGhSx+pYVZKIMLhDI37+SxwdmtThybmbuGv6WjJqRxY1Bjz3+3GuY2xwByvjKgdw7qJwrlVyBf/gz5w5k86dOxMdHc2ECRMoKChg7dq1tG3bluzsbDIzM2nTpg2bN2/mzJkz9O7dmw4dOhAVFcV3331XvJ4ZM2bQtm1b2rVrxx133MHKlSuZN28ejz32GNHR0ezevfu87c6ZM4fIyEjatWtHXFzRFbjZ2dncddddREVF0b59+/OuWD7n2Wef5eWXXy5+HBkZSUpKCpMmTWL37t1ER0fz2GOPkZKSQmRk5GXXO336dAYPHkzfvn1p3rw5jz/+eLm+lxc6mZXH3Z+s5Z3FuxnRuQmzx3clyM+BfYmqqYa1azBjXGcm39KG3/Yc47qvC1nR4ZWi34dF/zq/pbiq1px7TMEOk5ls3bqVL774ghUrVuDu7s59993HrFmzGDNmDAMHDuTpp58mKyuL0aNHExkZSX5+PnPnzsXPz4+jR4/StWtXBg4cSFJSEs899xwrV64kICCAjIwM6taty8CBAxkwYMBFD+NMnjyZn3/+meDgYE6cOAHAO++8g4iwadMmtm3bRp8+fdixY0eZvpYXXniBzZs3F3d8Ldne4nLr3bBhA+vXr8fT05OWLVvywAMP0Lhx46v+np6z8/Bpxn+aQOrxs/xrUCSjuoSUe53qf1xchDHXhBLbLIBH5/zOqF/zeLdhf/rHvwRxj2tBcBLOXRSg6Ac95m6ooB/8hQsXkpCQQKdOnYCiuRWCgop65z/zzDN06tQJLy8v3nzzTaCoJ8+TTz5JfHw8Li4upKWlcfjwYRYtWsTQoUOLm9TVrVu31G13796dsWPHMmzYMAYPHgwUtb9+4IEHAGjVqhUhISFlLgqXc7n19u7dm1q1agEQERHB3r17y10Uftp8iEe/3EANDzdm/6krMaGlfz/U1QkP9GHOhGv4/rsv6Pr7t3zoOoQxq6fiEdZDC4MT0KKQHF+0hxD3eNG/5fzBN8Zw55138vzzz//huWPHjnHmzBny8vLIzs7G29ubWbNmkZ6eTkJCAu7u7oSGhl51W+spU6awevVq5s+fT8eOHUlISCjT++zdVjs/P/+q11VYaHj91x28uWgX7RrX5v3RHalfSw8X2ZvbvuXcuuspkm+eypxl3ixMb8nUmaNxvf0TarS4zup4yo6ce0yh5HSLvZ7636GkCwefr0Dv3r356quvOHLkCAAZGRns3bsXgAkTJvDPf/6TUaNG8be//Q2AkydPEhQUhLu7O4sXLy5+ba9evZgzZ07xvAoZGUXTS5xrc30xu3fvpkuXLkyePJnAwED2799Pjx49mDVrFgA7duxg3759tGzZ8rz3hYaGkpiYCEBiYiLJycmlbqss6y2zSwz4Zy95hT/NWMebi3YxLKYRX4zvqgXBUdISYeh0wmL6Me+B7rTrMZDxWROZPucbVu+5+Fwfqnpw7qJg+8HnYq2Sr1JERATPPfccffr0oW3bttxwww0cPHiQGTNm4O7uzsiRI5k0aRJr165l0aJFjBo1inXr1hEVFcWMGTNo1apoZq02bdrw1FNP0bNnT9q1a8cjjzwCFE3b+Z///If27dv/YaD5scceIyoqisjISLp160a7du247777KCwsJCoqittvv53p06ef90ke4LbbbiMjI4M2bdrw9ttv06JFCwD8/f3p3r07kZGRPPbYY+e9pyzrLbOLDPgXfHEnT652Z+mOdCbf0oYXb2uLl7tzt6twqNiHi38vPN1cmdSvFY+Mv4fPPQYzfOpvPPd9Etl5BdZmVHahrbOVJf7w/3Rury3mbnJXT2VC9gNscm/LOyM70CXc37Kc6nyZOfk8/+NWZv62j+ZBPrw6LJqoRrWsjqWu0OVaZzv3noKqPMLiMB3HQfxLvJd5LccCuzBvYqwWhErG29ON526N4pNxnTmVncegd1fwxq87ySvQ5nrVhRYFVSmc3b6IMys+4I38QdztuYivbsynYe0aVsdSl9CzRSC/PNyTm9o24LVfd3DbeyvZdeSM1bFUBaiWRaEqHxJzBhf+/xxY/zO5n9/JvTkP4Nv373iP+hSPuePKNeCv7K9WTXfeGN6ed0Z2YH/GWW56cxmLtx2xOpYqJ7sVBRGZJiJHRGRziWXPikiaiGyw3fqXeO4JEdklIttF5Mar3a6XlxfHjh3TwlBJGWM4duwYXl5FZxEt2naYz7/7jsflEe4fdxfjYsOQ8J7lHvBXjnNT2wb8/Jc4mgX5MPGzRLYePGV1JFUOdhtoFpE44AwwwxgTaVv2LHDGGPPyBa+NAGYDnYGGwK9AC2PMZU9vuNhAc15eHqmpqeU+117Zj5eXF8HBwXywfC+vLNhBRAM/3r+jI43qaIfOquzQyWxueWc5riJ8O7E7Qb56+nBlZckczcaYeBEJLePLbwE+N8bkAMkisouiArHqSrfr7u5OWFjYlb5NOVBmTj4Pffk7P24+xK3RDXl+cFunnx2tOqhfy4uP7uzE0Cmr+NOMBL4Y31VPI66CrBhTmCgiG22Hl861tgwG9pd4Taptmapm9h7LZPC7K/l5yyGevqk1r90erQWhGokMrsXrw6PZmHqCR7/8nUKdvKfKcXRReA9oCkQDB4FXrnQFIjJeRNaJyLr09PQKjqfsaemOdG5+azmHT2czY1wX7ukRjohYHUtVsBvb1GdS31bM33SQ134tf58t5VgOLQrGmMPGmAJjTCEwlaJDRABpQMmOaY1syy62jg+MMTHGmJjAwED7BlYVwhjDe0t2c9fHa2hYuwb/nRhLbPMAq2MpOxofF87tMY15a9EuvklMtTqOugIOLQoi0qDEw0HAuTOT5gHDRcRTRMKA5sAaR2ZT9lFYaHj0y9958adt9I9qwDf3daNxXR1Qru5EhH/eGsk14f5M+noTa1MyrI6kysiep6TOpmiguKWIpIrI3cBLIrJJRDYC1wF/ATDGbAG+BJKAn4D7SzvzSFUNU+J38836NB7q3Zy3RrSnpoc25nUWHm4uvDe6A43q1GD8jHXsPZZpdSRVBtWu95GqPFbuPsroD1fTL6oBb49or+MHTir5aCa3vrOCAB8PvrmvO7VquFsdyelp7yPlcIdPZfPg7PWEBXjz4m1ttSA4sbAAb6aM7si+jLPcPytR+yRVcloUVIXLKyhk4meJZOYU8N7ojvh46iEjZ3dNU3/+NSiK5buO8vd5W7TjQCWmv62qwr300zbWphznjeHRtKjna3UcVUkMi2nMnvRMpizdTdNAH+6O1YtMKyMtCqpC/bT5IFOXJTPmmhBuidbrD9X5Hr+xJclHz/Dc/CRC/WvSu3U9qyOpC+jhI1Vhko9m8ticjbRrXJunbtKJjtQfubgIr90eTZuGfjwwez1JB7R5XmWjRUFViKzcAv48MwE3V+HdUR3wdNPWFerianq48dGdnfDzcueeT9Zy5JQ2r6xMtCiocjPG8NS3m9h++DSvD29PsE6Oo0pRz8+LD++M4fjZPP40Yx1ZuXpZUmWhRUGV2+w1+/kmMY0HezWnZwttPaLKJjK4Fm8Mj2Zj2kkenbNBm+dVEloUVLlsSj3Js/O20KN5AA/2bm51HFXF9GlTnyf6teKHTYd4dYE2z6sM9OwjddVOnM3lz7MSCPDx4I3h7XF10QvU1JX7U49w9qRn8vbiXYQFeHNbx0ZWR3JquqegrkphoeGRL3/n8Kls3hnVgbreHlZHUlXUueZ53Zr6M+mbjaxJ1uZ5VtKioK7Ke0t3s2jbEZ6+KYL2TeqU/galLsPd1YX3RnWkcZ2aTPhUm+dZSYuCumIrdh3llV+2c3O7hoy5JsTqOKqaqFXTnWljO2GAcdPXcjIrz+pITkmLgroih04WNboLD/ThhcFR2uhOVajQAG/etzXPu29WgjbPs4AWBVVm5xrdZeUVMGV0B7y10Z2ygy7h/vx7UBQrdh3jme+0eZ6jaVFQZfbCj9tYt/c4L97WlmZB2uhO2c/QmMb8+dqm+CW8ww/zvjz/yeR4WP66JbmcgRYFVSY/bDrIR8uTGdstlJvbNbQ6jnICj/VpiVdIDF0THyVhyXdFC5PjYc5YCO5gabbqTPf/Val2p5/h8a820r5JbZ7sr43ulGO4uAj3jr2L597O4dHF93PkVAJB22bC0OkQFmd1vGpL9xTUZZ3NzefPMxPwcHPhnZEd8HDTHxnlODU8XHnonrv51q0vQYlvkBl1pxYEO9PfcHVJxhiemruZnUfO8MbwaBpqoztlgaBjaxjjvpB3zW3kr/mQnB1LrI5UrWlRUJc0a/U+5q5P4+HeLejRXBvdKQvYxhDcbv+EFrc/z705D5D7+RgKdy+1Olm1ZbeiICLTROSIiGwusew/IrJNRDaKyFwRqW1bHioiWSKywXabYq9cqmw2pp5g8n+T6NkikAd6NbM6jnJWaYnFYwjXR9Sjd78hjM+eyLKlv1idrNqy557CdKDvBcsWAJHGmLbADuCJEs/tNsZE22732jGXKsXxzFz+PDORQF9PXr89GhdtdKesEvvweWMId8eGERrTjzt3dOOrhFTrclVjdisKxph4IOOCZb8YY/JtD38DtB1iJVNYaPjLlxs4crqo0V0dbXSnKhERYfItbejezJ8nvtnI6j3HrI5U7Vg5pjAO+LHE4zARWS8iS0Wkx6XeJCLjRWSdiKxLT0+3f8rqbPnrRcdsS/h27ue03DWNZwZEEN24tiWxlLocd1cX3h3ZkcZ1azJhZgIpR7V5XkWypCiIyFNAPjDLtugg0MQY0x54BPhMRPwu9l5jzAfGmBhjTExgoA5+lktwh6ILgWyFYdOy/9Jz42P4hndmdFdtdKcqr1o13fl4bCcEGPfJWk6e1eZ5FcXhRUFExgIDgFHG1tTEGJNjjDlmu58A7AZaODqb0wmLKxrEmzOW0z8+S6OF9/GCz98Yd8cYbXSnKr0Qf2+mjO7I/oyzPPzFeu2RVEEcWhREpC/wODDQGHO2xPJAEXG13Q8HmgN7HJnNaYXFUdBhHL6rX+NzcwP3jr2Lmh56obuqGrqE+zOpX2sWb0/n5y2HrI5TLdjzlNTZwCqgpYikisjdwNuAL7DgglNP44CNIrIB+Aq41xij0y85QnI82b9N5Y38QYzzXETTM4lWJ1Lqitx5TQit6vvyz++3kpVbYHWcKs+eZx+NMMY0MMa4G2MaGWM+MsY0M8Y0vvDUU2PM18aYNrZlHYwx/7VXLlVCcjw5s8dwz9n7OdHlMTxHzDhvjEGpqsDN1YV/DGxD2oks3lm8y+o4VZ5e0ezEju74jXuzJ5LTqDtP9Gv9vzGGNN1bUFVLl3B/bo1uyAfxe/RspHLSouCkMnPyGbGlC7+7t+OdUSUa3YXFFV0wpFQV80T/1ri7Cv/4r07MUx5aFJyQMYYn525iV/oZ3hzenga1tNGdqvrq+Xnx8PUtWLw9nYVbj1gdp8rSouCEZv62l+82HOCR61sQ2zzA6jhKVZix3UNpFuTDP77fQnaeDjpfDS0KTmbD/hNM/j6J61oGcv912uhOVS/uri5MHtiG/RlZTFm62+o4VZIWBSdyPDOX+2clEuTrxWva6E5VU92aBXBT2wa8t2Q3+zPOlv4GdR4tCk6isNDw8BcbSD+dw3ujO1C7pja6U9XX0ze1xkWEyd8nWR2lytGi4CRmrdnH0h3p/H1gBG0b1bY6jlJ21aBWDR7o3YwFSYdZvF0Hna+EFgUncCo7j9cW7KBreF1Gdm5idRylHOLu2DDCA7z5x7wt5OTroHNZaVFwAu8t2U1GZi5P3xShje6U0/B0c+XvA9uQcuwsHy5LtjpOlaFFoZpLPX6Wj5YnM7h9MJHBtayOo5RD9WwRyI1t6vHWop2knciyOk6VoEWhmnv55+0I8NcbW1odRSlL/N+ACIyBf83XQeey0KJQjf2+/wTfbjjAPT3CaFhbr1pWzqlRnZrcf10zfth0iOU7j1odp9LTolBNGWP41/ytBPh48Odr9SI15dzGx4UT4l+Tv8/bTG5+odVxKjUtCtXUL0mHWZOSwV9uaIGPp06ao5ybl7srf785gt3pmXy8QgedL0eLQjWUm1/ICz9uo3mQD7fHNLY6jlKVQq9W9ejdKog3F+7k0Mlsq+NUWloUqqFZq/eSfDSTJ/u3xs1V/4uVOueZmyPIKzT8+4etVkeptPQvRjVzMiuPNxbuJLZZANe2DLQ6jlKVSoi/N/fGhTPv9wOs2n3M6jiVkhaFauadxbs4mZXHk/1b64VqSl3En69tRnDtGjw7bwt5BTrofCEtCtXI/oyzTF+RwpAOjYho6Gd1HKUqpRoerjxzcwTbD59mxqq9VsepdC57WoqIdCjDOvKMMZsqKI8qhxd/2oari/BoH71QTanL6RNRj7gWgby+YAc3t2tAkK+X1ZEqjdLOVVwKrAUudxwiDAi92BMiMg0YABwxxkTaltUFvrC9JwUYZow5LkXHOt4A+gNngbHGGJ1BvowS9x3n+40HebB3c+rX0h9wpS5HRHj25ghufD2eF37cxqvDoq2OVGmUdvhorTGmlzHmukvdgD2Xef90oO8FyyYBC40xzYGFtscA/YDmttt44L0r/WKc1bkL1QJ9PZkQF251HKWqhPBAH/7UI5xvEtNYl5JhdZxK47JFwRjTq7QVXO41xph44MLv9i3AJ7b7nwC3llg+wxT5DagtIg1K276CHzcfImHvcR69oQXeeqGaUmU2sVczGtTy4v++20JBobE6TqVw2aIgIjeKyJCLLB8iIjdc5TbrGWMO2u4fAurZ7gcD+0u8LtW27MJtjxeRdSKyLj09/SojVB85+QW88OM2WtbzZaheqKbUFanp4cbTN0Ww9eApZq3WQWco/fDRMxSNK1xoCTC5vBs3xhjgisqzMeYDY0yMMSYmMFDPw/901V72ZZzlyZta46pzLit1xfpH1adbU39e/nk7x87kWB3HcqUVBU9jzB8+jhtjjgLeV7nNw+cOC9n+PTdXXhpQ8qNuI9sydQknzuby1qJdxLUIpGcLLZBKXQ0RYfItbTibW8BLP223Oo7lSisKfiLyh4PUIuIOXG0v5nnAnbb7dwLflVg+Rop0BU6WOMykLuLNhbs4nZ3HU/1bWx1FqSqtWZAv42LD+GLdftbvO251HEuVVhS+AaaKSPFegYj4AFNsz12WiMwGVgEtRSRVRO4GXgBuEJGdwPW2xwA/UHQm0y5gKnDfFX4tTiXlaCaf/pbCsJjGtKzva3Ucpaq8B3s3J8jXk2ecfNC5tFNVngaeA/aKyLlRmCbAR8D/lbZyY8yISzzV+yKvNcD9pa1TFXnxp224u7rwSJ8WVkdRqlrw8XTjqZta89DnG/hi7X5GdmlidSRLlHZKar4xZhJFx/rH2m5NjDGTjDF59o+nLmZdSgY/bj7EvT2b6pWYSlWgge0a0jmsLi/9vI3jmblWx7FEqb2PRKQFRYVgk+2ms19byBjDc/O3Us/Pk3t6hFkdR6lq5dyg8+nsfF7+xTkHnUu7TuEFiq4uvltE/uOYSOpy/rvxIBv2n+CvfVpS00MvVFOqorWq78eYa0L4bM0+NqWetDqOw5W2p9AbeIqi8YNSr25W9pWdV8CLP26jdQM/BndoZHUcpaqtv9zQAn9vD56Zt5lCJxt0Lq0o/JWiVhQzgCfsH0ddzicrU0g7kcXTeqGaUnbl5+XOpH6tWb/vBF8lplodx6FKG2heaowZboy53Rjzi6NCqT/KyMzl7cW76NUqiO7NAqyOo1S1N7h9MB1D6vDij9s4edZ5zqspbUyhfmkrKMtrVPm9uXAnZ3MLeKJfK6ujKOUUXFyEfwxsw/Gzubz26w6r4zhMaYePfijDOsryGlUOe9LPMPO3vQzv1Jjm9fRCNaUcJTK4FqO6hDBjVQpJB05ZHcchSisK7UTk1GVup/lfl1NlJ8//uA0vd1f+coNeqKaUoz3apwUPef3Al3NmUXSNrU1yPCx/3bJc9lLamIKrMcbvMjdfY8wf2lurivPbnmMsSDrMn69tSoCPp9VxlHI6tWt6ENX5Oh7I+BfLF9i6+yTHw5yxEFyWGYurFj3RvRIrLCyaUa1BLS/ujtUL1ZSyyrU33sYz247w6MqHyWELnhumw9DpEBZndbQKV9pAsxYNC837/QCb0k7y2I0t8XJ3tTqOUk7LxUUYOmQkn+b3xnPlKxBzd7UsCFD6mMIah6RQf5CdV8BLP20jKrgWt0brETqlrNYufyPjPBfxVsFg8td8WHQIqRoqrSjoFVIW+Wh5MgdOZvNk/9a46IVqSlnLNoZQOGQ6H7qN4N/ef8PMGVstC0Nph4cCReSRSz1pjHm1gvMo4OiZHN5bspvrW9fjmqb+VsdRSqUlwtDp+IXF8dcb9/J/3+Zx/fWv0C0tsdodRiptT8EV8AF8L3FTdvD6rzvIyivgif56oZpSlULsw8V//Ed2bkJksB+PrPEjs9NEa3PZQWl7CgeNMZMdkkQBsOvIaWav2c+oLk1oGuhjdRyl1AVcXYR/DIzktvdW8taiXUyqZl0GdEyhknn+h23UdHflod7NrY6ilLqEjiF1GNKxER8t30PK0Uyr41SosrTOVg6yctdRFm47wv29muGvF6opVak9fmNLBGHqsj1WR6lQpV3RnOGoIM6uoLBoRrXg2jUY2y3U6jhKqVIE+XlxW8dg5iSkkn46x+o4FabU6TiVY8xdn0bSwVM83lcvVFOqqrinRzh5BYXMWJVidZQK4/CiICItRWRDidspEXlYRJ4VkbQSy/s7OptVsnILePnn7bRrXJuB7RpaHUcpVUZNA324oXU9ZqzaS2ZOvtVxKoTDi4IxZrsxJtoYEw10BM4Cc21Pv3buOWNM9W3Jvfz18y56+XDZHsLOJPB2k6WI6Ni+UlXJhJ5NOZmVxxdr91sdpUJYffioN7DbGLPX4hyOFdyhqMNicjxHTmeTsPQ7PvB6m8aRsVYnU0pdoY4hdegUWoePlieTV1BodZxys7ooDAdml3g8UUQ2isg0EalzsTeIyHgRWSci69LT0x2TsqKFxRV1WJwzlk2fPs6r8jpnBn5Y7a6MVMpZTIhrStqJLH7YdNDqKOVmWVEQEQ9gIDDHtug9oCkQDRwEXrnY+4wxHxhjYowxMYGBgY6Iah9hcRxtPZreRz5hU8MhNIjuY3UipdRV6tUqiKaB3kxZuuf8iXiqICv3FPoBicaYwwDGmMPGmAJjTCEwFehsYTb7S47HY/10pnAbPU7Mq5aNtZRyFi4uwoS4pmw9eIplO49aHadcrCwKIyhx6EhEGpR4bhCw2eGJHCU5ntzPxzAheyKuvZ7GZdj04jEGpVTVdEv7hgT5evJBfNW+mM2SoiAi3sANwDclFr8kIptEZCNwHfAXK7I5QmFqIk+7/ZW02p0Y0y3kf2MMaYlWR1NKXSVPN1fGxYaxfNdRNqedtDrOVbOkKBhjMo0x/saYkyWW3WGMiTLGtDXGDDTGVP0Rm0uY4zWYL4+F8be+rfB0s12oFhZX1IlRKVVljezSBB9PN96vwnsLVp995HQyc/J55ZcddGhSm/5R9a2Oo5SqQH5e7ozs0oT5Gw+wP+Os1XGuihYFB5u6bA9HTufw1E0ReqGaUtXQXd1DcXURPqyijfK0KDhQZk4+05Ync2ObenQMuehlGEqpKq5BrRrcEh3MF+v2k5GZa3WcK6ZFwYG+SUzlVHY+4+OaWh1FKWVH4+PCyc4r5NNVVa9ZgxYFByksNHy8IoV2jWvToUltq+MopeyoRT1fercK4pNVKWTlFlgd54poUXCQpTvT2XM0k3HdQ3UsQSknMD4unIzMXL5KqFqN8rQoOMi05cnU8/OkX2SD0l+slKryOofVJbpxbaYuS6agsOq0vtCi4AA7Dp9m2c6jjLkmFA83/ZYr5QxEhHt7hrMv4yw/bT5kdZwy079QDvDxihQ83VwY0bmJ1VGUUg50Q0R9wgK8mbJ0d5VplKdFwc6OZ+Yyd30qg9oHU9fbw+o4SikHcnUR/tQjnE1pJ1m155jVccpEi4KdzV67j+y8Qu7qHmZ1FKWUBQZ3CCbAx4P3l1aNi9m0KNhRXkEhM1buJbZZAC3r+1odRyllAS93V8Z2C2XpjnS2HjxldZxSaVGwox83H+LQqWzGxYZaHUUpZaHRXUOo6eHK1CrQKE+Lgh19vCKZsABvrm0RZHUUpZSFatf0YHinJsz7/QBpJ7KsjnNZWhTsJHHfcdbvO8HYbqG4uOjFako5u7t7hGEoumapMtOiYCcfr0jB18uNIR0bWR1FKVUJBNeuwc1tGzB7zT5Ons2zOs4laVGwg4Mns/hh00GGd2qMt6eb1XGUUpXE+LimnM0tYObqytsoT4uCHXy6ai/GGMZcE2p1FKVUJRLR0I+4FoF8vCKF7LzK2ShPi0IFy8ot4LM1++gTUZ/GdWtaHUcpVcncGxfO0TM5zF2fZnWUi9KiUMHmrk/jxNk8xsXqxWpKqT+6pqk/UcG1mBq/p1I2ytOiUIGMMXy8IpnIYD86herMakqpPxIRxseFs+doJguSDlsd5w8sKwoikiIim0Rkg4issy2rKyILRGSn7d8q9Zd1+a6j7Dxyhru6hemcCUqpS+oXWZ/GdWvwfnzla5Rn9Z7CdcaYaGNMjO3xJGChMaY5sND2uMqYtjyZAB9PBrTTOROUUpfm5urCn3qEs37fCdbtPW51nPNYXRQudAvwie3+J8Ct1kW5MrvTz7B4ezp3dA3B083V6jhKqUpuaMfG1KnpzvtLd1sd5TxWFgUD/CIiCSIy3rasnjHmoO3+IaDehW8SkfEisk5E1qWnpzsqa6mmr0jBw9WFUV11zgSlVOlqeLhyZ7dQft16hJ2HT1sdp5iVRSHWGNMB6AfcLyJxJZ80RQfa/nCwzRjzgTEmxhgTExgY6KCol3fybB5fJaQyMLohAT6eVsdRSlURY64JxcvdhQ8qUaM8y4qCMSbN9u8RYC7QGTgsIg0AbP8esSrflfhi3T6y8gq4q3uo1VGUUlVIXW8PhsU05tsNaRw6mW11HMCioiAi3iLie+4+0AfYDMwD7rS97E7gOyvyXYn8gkI+WbmXruF1adOwltVxlFJVzD2x4RQUGj5eWTka5Vm1p1APWC4ivwNrgPnGmJ+AF4AbRGQncL3tcaX2S9Jh0k5kMU5nVlNKXYUm/jXpH9WAz37bx6ls6xvlWdKtzRizB2h3keXHgN6OT3T1Pl6RTOO6Nejd+g9j4kopVSYT4pry/caDzF69jwk9m1qapbKdklqlbEw9wdqU44ztFoarzpmglLpKUY1q0a2pP9NWJJOTb22jPC0K5fDxihR8PN0YFqNzJiilymdCz6YcPpXDdxsOWJpDi8JVOnIqm+83HmBoTCN8vdytjqOUquLimgfQqr4vU+P3UGhhozwtCldp5m97yS80jO0WanUUpVQ1ICLc27MpO4+cYfF2687G16JwFbLzCpi5eh+9W9UjxN/b6jhKqWriprYNCK5dg/eXWncxmxaFqzBvwwEyMnMZFxtqdRSlVDXi7urCuNgw1qRkkLjPmkZ5WhSukDGGaSuSaVXfl2vC/a2Oo5SqZoZ3akytGu58YNHeghaFK7RqzzG2HTrNuO46Z4JSquJ5e7pxR9cQfk46xJ70Mw7fvhaFKzRteQp1vT0YGN3Q6ihKqWrqzm6huLu6MHWZ41tfaFG4AilHM1m47TCjuzTBy13nTFBK2UegrydDOjbi68RUjpx2bKM8LQpXYPrKFNxchNFdQ6yOopSq5v7UI5y8gkI+WZni0O1qUSij09lFcyYMaNuQID8vq+Mopaq5sABvboyoz6er9pKZk++w7WpRKKMv16VyJidfu6EqpRxmQs9wTmXn8/na/Q7bphaFMigoNExfmUyn0DpENdI5E5RSjtG+SR06h9Xlo2V7yCsodMg2tSiUwa9bD7M/Q+dMUEo53r09wzlwsqjXmiNoUSiDj1ckE1y7BjdE6JwJSinHurZFEM2DfHh/6R6Kpq63Ly0Kpdhy4CS/7cngzm4huLnqt0sp5VguLsL4uHC2HTpN/M6j9t+e3bdQxX28IoWaHq7cHtPE6ihKKSd1S3Qw9fw8eX/pbrtvS4vCZaSfzmHehgMM6diIWjV1zgSllDU83Fy4OzaMlbuPsTH1hF23pUXhMj5bvY/cgkLu1DkTlFIWG9G5Cb6ebrwfb99GeVoULiEnv4BPf9vLdS0DaRroY3UcpZST8/VyZ2TXJvy46SD7jp2123YcXhREpLGILBaRJBHZIiIP2ZY/KyJpIrLBduvv6Gwlff/7QY6eyWFcrJ6GqpSqHCZ6zKebaxIfLi+xt5AcD8tfr7BtWLGnkA88aoyJALoC94tIhO2514wx0bbbDxZkA/43Z0LzIB9imwVYFUMppc7jG96ZKZ5vsXfdTxw7k1NUEOaMheAOFbYNtwpbUxkZYw4CB233T4vIViDY0TkuZ23KcbYcOMW/B0XpnAlKqcojLI6TAz7g1W/vYefnB/DP+BaGToewuArbhKVjCiISCrQHVtsWTRSRjSIyTUTqXOI940VknYisS09Pt0uuacuTqV3TnUHtK1WtUkopgtv3ZW3AILqmfgQxd1doQQALi4KI+ABfAw8bY04B7wFNgWiK9iReudj7jDEfGGNijDExgYGBFZ5rf8ZZfkk6xMjOTajhoXMmKKUqmeR4+mbNh7jHYd1HRYeQKpAlRUFE3CkqCLOMMd8AGGMOG2MKjDGFwFSgsxXZZqxKQUS44xqdM0EpVcmcG0MYOh16PVX075yxFVoYrDj7SICPgK3GmFdLLG9Q4mWDgM2OznYmp6hFbf+oBjSoVcPRm1dKqctLSzx/DCEsruhxWmKFbcLhA81Ad+AOYJOIbLAtexIYISLRgAFSgAmODvZ1Qiqns/MZ1z3U0ZtWSqnSxT78x2VhcRU6rmDF2UfLgYud0mPZKagAhYWG6StTaN+kNu2bXHSMWymlqj29otlmyY4jJB/N5C6dM0Ep5cS0KNhMW55CfT8v+kXWtzqKUkpZRosCsP3QaZbvOsqYbiG465wJSiknpn8BKZpZzcvdhRGddM4EpZRzc/qikJGZy9z1aQxq34g63h5Wx1FKKUs5fVGYvWYfOfmFehqqUkrh5EUhN7+QGatS6NE8gOb1fK2Oo5RSlnPqovDj5oMcPqVzJiil1DlOWxSMMUxbnkx4gDc9m1d8Yz2llKqKnK8oLH8dkuNJ3HeC31NPclf3UFz2LqvQmYuUUqqqcr6iENwB5oxl2YJv8PNyY4h/coXPXKSUUlWVFQ3xrBUWx9F+73PHV3cREzyUGt9+V+EzFymlVFXlfHsKwIl6XVleeyCxB6bZZeYipZSqqpyyKDTLXM8teT/ZbeYipZSqqpyvKDhg5iKllKqqnK8oOGDmIqWUqqqcb6DZATMXKaVUVeV8ewpKKaUuSYuCUkqpYloUlFJKFdOioJRSqpgWBaWUUsXEGGN1hqsmIunA3nKsIgA4WkFxqgJn+3pBv2ZnoV/zlQkxxly0PXSVLgrlJSLrjDExVudwFGf7ekG/ZmehX3PF0cNHSimlimlRUEopVczZi8IHVgdwMGf7ekG/ZmehX3MFceoxBaWUUudz9j0FpZRSJWhRUEopVcwpi4KI9BWR7SKyS0QmWZ3H3kRkmogcEZHNVmdxFBFpLCKLRSRJRLaIyENWZ7I3EfESkTUi8rvta/6H1ZkcQURcRWS9iHxvdRZHEZEUEdkkIhtEZF2FrtvZxhRExBXYAdwApAJrgRHGmCRLg9mRiMQBZ4AZxphIq/M4gog0ABoYYxJFxBdIAG6t5v/PAngbY86IiDuwHHjIGPObxdHsSkQeAWIAP2PMAKvzOIKIpAAxxpgKv2DPGfcUOgO7jDF7jDG5wOfALRZnsitjTDyQYXUORzLGHDTGJNrunwa2AsHWprIvU+SM7aG77VatP/WJSCPgJuBDq7NUF85YFIKB/SUep1LN/1g4OxEJBdoDqy2OYne2QykbgCPAAmNMdf+aXwceBwotzuFoBvhFRBJEZHxFrtgZi4JyIiLiA3wNPGyMOWV1HnszxhQYY6KBRkBnEam2hwtFZABwxBiTYHUWC8QaYzoA/YD7bYeIK4QzFoU0oHGJx41sy1Q1Yzuu/jUwyxjzjdV5HMkYcwJYDPS1OIo9dQcG2o6vfw70EpGZ1kZyDGNMmu3fI8Bcig6LVwhnLAprgeYiEiYiHsBwYJ7FmVQFsw26fgRsNca8anUeRxCRQBGpbbtfg6KTKbZZGsqOjDFPGGMaGWNCKfo9XmSMGW1xLLsTEW/byROIiDfQB6iwMwudrigYY/KBicDPFA0+fmmM2WJtKvsSkdnAKqCliKSKyN1WZ3KA7sAdFH163GC79bc6lJ01ABaLyEaKPvwsMMY4zWmaTqQesFxEfgfWAPONMT9V1Mqd7pRUpZRSl+Z0ewpKKaUuTYuCUkqpYloUlFJKFdOioJRSqpgWBaWUUsW0KCillCqmRUGpCiAioSKSZes7dCXvu93Wwl2vJ1CVghYFpSrOblvfoTIzxnwB3GOfOEpdOS0KSpVCRDqJyEbbJDbetglsLttozrbnsE1EpovIDhGZJSLXi8gKEdkpIhXWq0apiuRmdQClKjtjzFoRmQc8B9QAZhpjytJrphkwFBhHUduJkUAsMBB4ErjVLoGVKgctCkqVzWSK/rBnAw+W8T3JxphNACKyBVhojDEisgkItUtKpcpJDx8pVTb+gA/gC3iV8T05Je4XlnhciH4gU5WUFgWlyuZ94P+AWcCLFmdRym7004pSpRCRMUCeMeYzEXEFVopIL2PMIquzKVXRtHW2UhXANg/098aYK57+UkSuBf5qjBlQwbGUumJ6+EipilEA1Lqai9eAd4Hj9gil1JXSPQWllFLFdE9BKaVUMS0KSimlimlRUEopVUyLglJKqWL/DzoqDEX8ukZ4AAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# --------------\n",
    "# Postprocessing\n",
    "# --------------\n",
    "\n",
    "# Define exact solution:\n",
    "x_ex = np.linspace(0, L, nx)\n",
    "T_ex = [-(U*A)**2/(2*rho*k)*x_ex_i**2 + ((T_h-T_c)/L + (U*A)**2/(2*rho*k)*L)*x_ex_i + T_c for x_ex_i in x_ex]\n",
    "\n",
    "# Graph solutions:\n",
    "fe.plot(T-273.15)\n",
    "plt.plot(x_ex, T_ex-273.15*np.ones(len(T_ex)), \"x\")\n",
    "plt.xlabel(\"x [m]\")\n",
    "plt.ylabel(\"T [°C]\")\n",
    "plt.legend([\"FEM solution\",\"exact solution\"])\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Benchmarking results\n",
    "\n",
    "For many trivial problems, ours included, there is an exact solution. We can use this solution to benchmark our numerical results. The previous cell introduced the visual comparison of the explicit and discretized solutions.\n",
    "To make the benchmarking process quantitative, we use FEniCS' ``errornorm`` functionality.\n",
    "In order to benchmark the explicit solution, we first need to represent as an object that FEniCS can handle. Basically, there are two objects in the ``fenics`` library that are expedient to our goal.\n",
    "\n",
    "> **Expression** vs. **UFL**. FEniCS object ``Expression`` can be understood as an entity containing *precise* representation of some mathematical formula (e.g., the function $\\sin(x)$), whereas UFL obejct is essentialy representation of that said formula that belongs to the finite element space $V_h$.\n",
    "It is worth noting that anytime you want to work with expressions within the framework of FEM, you need to modify them so that they, again, belong to $V_h$.\n",
    "This can also be done in two ways.\n",
    "\n",
    "> ``project`` vs ``interpolate``. The result of ``interpolate`` is, simply put, a FE element function from V_h, that coincides with the original object in the nodes of the finite element (in our case of piece-wise linear functions, these are endpoints of every subinterval).\n",
    "``projection`` on the other hand does not have to coincide in the nodal points, but it is in some sense the best approximation to the original object that lies in $V_h$.\n",
    "More precisely it the result of the orthogonal projection of the object onto $V_h$, mathematically it is a function u_h \\in V_h that satisfies\n",
    "\\begin{align*}\n",
    "\\int_\\Omega (u_h - \\sin(x))v_h\\,\\mathrm{d}x = 0\\,, \\forall v_h \\in V_h\\,.\n",
    "\\end{align*}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 184,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "2.5805735030741737\n"
     ]
    }
   ],
   "source": [
    "# Define expression containg the explicit solution:\n",
    "T_sol = fe.Expression(\"-U*U*A*A/(2.0*rho*k)*x[0]*x[0] + ((T_h - T_c)/L + U*U*A*A/(2*rho*k)*L)*x[0] + T_c\",\n",
    "                      U = U, A = A, rho = rho, k = k, T_h = T_h, T_c = T_c, L = L,\n",
    "                     degree = 3)\n",
    "\n",
    "#T_sol_prj = fe.project(T_sol, V_h)\n",
    "\n",
    "# It is recommended to use a finer mesh for comparison:\n",
    "mesh_fine = fe.IntervalMesh(3*nx, 0.0, L)\n",
    "\n",
    "err = fe.errornorm(T_sol, T, mesh = mesh_fine)\n",
    "\n",
    "print(err)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Questions and excercises\n",
    "\n",
    "\n",
    "Q.1 Are all the properties and variables really constant over the cross-section? Is the one-dimensional simplication justified?\n",
    "<br/>Q.2 Consider the wire to be thermally insulated on both ends, i.e., $\\partial T/\\partial x(0) = \\partial T/\\partial x(L) = 0$, is such problem well defined? If not, how can we make it so?\n",
    "\n",
    "E.1 Study the role of finite element degree, possibly on some coarse mesh. What should the graph of the solution look like? Based on the values of the error is it more efficient to raise degree of the finite element, or the number of the mesh nodes?\n",
    "<br/>E.2 Assume that the control unit (and effectively the left end of the wire) is cooled by a fan blowing air at the temperature $T_0$, formulate the problem with this boundary condition.\n",
    "<br/>E.3 Try to the implement and test the formulation from Q.2.\n",
    "<br/>E.4 Study the difference between projection and interpolation on the step function $$f(x) = \\begin{cases} 0 & x <0 \\\\ 1/2 & x= 0 \\\\ 1 & x >0\\end{cases}$$\n"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Complete code"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Import relevant modules:\n",
    "import fenics as fe # 'import dolfin' does the same\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "\n",
    "# --------------------\n",
    "# Parameters\n",
    "# --------------------\n",
    "# Geometric parameters:\n",
    "L = 5.0  # Length of the wire [m]\n",
    "R = 0.005 # Radius of the wire [m]\n",
    "A = np.pi*R**2 # Cross-section area [m^2]\n",
    "\n",
    "# Environment parameters:\n",
    "T_c = 283 # Ambient temperature of the cellar [K]\n",
    "T_h = 298 # Ambient temperature of the outside air [K]\n",
    "U = 230 # Voltage drop [V]\n",
    "\n",
    "# Material parameters:\n",
    "rho = 1.68e-8 # Electrical resistivity of copper [Ω m]\n",
    "k = 384.1 # Thermal conductivity of copper [W m^-1 K^-1]\n",
    "\n",
    "# --------------------\n",
    "# Geometry\n",
    "# --------------------\n",
    "nx = 10 # Number of mesh subintervals\n",
    "mesh = fe.IntervalMesh(nx, 0.0, L)\n",
    "\n",
    "fe.plot(mesh)\n",
    "plt.show() # this forces the plot to show in a pop-up window when in command line/script environment\n",
    "\n",
    "# ---------\n",
    "# FE spaces\n",
    "# ---------\n",
    "V_h = fe.FunctionSpace(mesh, \"CG\", 1)\n",
    "\n",
    "# Definition of trial (=solution) and test functions:\n",
    "_T = fe.TrialFunction(V_h)\n",
    "Theta = fe.TestFunction(V_h)\n",
    "\n",
    "# ------------------------\n",
    "# Boundary markings\n",
    "# ------------------------\n",
    "def left(x, on_boundary):\n",
    "    return fe.near(x[0], 0.0) and on_boundary\n",
    "# analogously, we could define the marking of the right boundary\n",
    "def right(x, on_boundary):\n",
    "    return fe.near(x[0], L) and on_boundary\n",
    "\n",
    "boundary = fe.MeshFunction(\"size_t\", mesh, mesh.topology().dim()-1,0)\n",
    "for v in fe.facets(mesh):\n",
    "    if fe.near(v.midpoint()[0], 0.0):\n",
    "        boundary[v] = 1 # left boundary\n",
    "    elif fe.near(v.midpoint()[0], L):\n",
    "        boundary[v] = 2 # right boundary\n",
    "\n",
    "# --------------------\n",
    "# Boundary conditions\n",
    "# --------------------\n",
    "bc_left = fe.DirichletBC(V_h, T_c, left) # prescribe cellar temperature to the left end of the wire\n",
    "bc_right = fe.DirichletBC(V_h, T_h, right) # prescribe outdoor temperature to the right end of the wire\n",
    "\n",
    "# Group boundary conditions into an array:\n",
    "bcs = [bc_left, bc_right]\n",
    "\n",
    "# Define volume and surface measures:\n",
    "dx = fe.Measure(\"dx\", mesh)\n",
    "ds = fe.Measure(\"ds\", subdomain_data = boundary)\n",
    "\n",
    "# --------------------\n",
    "# Weak form\n",
    "# --------------------\n",
    "#a = E*A*fe.inner(fe.grad(T),fe.grad(Theta))*dx\n",
    "# One can equivalently write the bilinear form in this one-dimensional problem as\n",
    "a = _T.dx(0)*Theta.dx(0)*fe.dx\n",
    "l = A**2*U**2/(rho*k)*Theta*dx\n",
    "\n",
    "# --------------------\n",
    "# Solver\n",
    "# --------------------\n",
    "T = fe.Function(V_h)\n",
    "fe.solve(a == l, T, bcs)\n",
    "# Equivalent implementation:\n",
    "# F = a-l\n",
    "# problem = fe.LinearVariationalProblem(fe.lhs(F),fe.rhs(F),u,bc)\n",
    "# solver = fe.LinearVariationalSolver(problem)\n",
    "\n",
    "# --------------\n",
    "# Postprocessing\n",
    "# --------------\n",
    "\n",
    "# Define exact solution:\n",
    "x_ex = np.linspace(0, L, nx)\n",
    "T_ex = [-(U*A)**2/(2*rho*k)*x_ex_i**2 + ((T_h-T_c)/L + (U*A)**2/(2*rho*k)*L)*x_ex_i + T_c for x_ex_i in x_ex]\n",
    "\n",
    "# Graph solutions:\n",
    "fe.plot(T-273.15)\n",
    "plt.plot(x_ex, T_ex-273.15*np.ones(len(T_ex)), \"x\")\n",
    "plt.xlabel(\"x [m]\")\n",
    "plt.ylabel(\"T [°C]\")\n",
    "plt.legend([\"FEM solution\",\"exact solution\"])\n",
    "plt.show()\n",
    "\n",
    "# Define expression containg the explicit solution:\n",
    "T_sol = fe.Expression(\"-U*U*A*A/(2.0*rho*k)*x[0]*x[0] + ((T_h - T_c)/L + U*U*A*A/(2*rho*k)*L)*x[0] + T_c\",\n",
    "                      U = U, A = A, rho = rho, k = k, T_h = T_h, T_c = T_c, L = L,\n",
    "                     degree = 3)\n",
    "\n",
    "#T_sol_prj = fe.project(T_sol, V_h)\n",
    "\n",
    "# It is recommended to use a finer mesh for comparison:\n",
    "mesh_fine = fe.IntervalMesh(3*nx, 0.0, L)\n",
    "\n",
    "err = fe.errornorm(T_sol, T, mesh = mesh_fine)\n",
    "\n",
    "print(err)"
   ]
  }
 ],
 "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.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}