{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Cork stopper - 'Hello World!'\n",
    "***\n",
    "***"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "![caption](fig/cork_stopper.png \"Cork stopper\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## General formulation of the problem\n",
    "\n",
    "---\n",
    "Let us consider a beam $\\Omega$ with a rectangular cross-section that is constant along the centerline of the beam, which coincides with the $x$-axis.  \n",
    "The beam is loaded with the volumetric force density $\\mathbf{b}$ acting only in the direction of the $x$-axis. This results in extension or compression of the beam. Additionally, the right end, denoted $\\Gamma_{\\text{N}}$, is loaded with the constant surface force density $\\mathbf{g}$. The beam is fixed at the left end denoted by $\\Gamma_{\\text{D}}$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "hide_input": true,
    "nbsphinx": "hidden"
   },
   "outputs": [],
   "source": [
    "%load_ext tikzmagic"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "alignment": "center"
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAlgAAAC8EAYAAAD0vAwiAAAJJmlDQ1BpY2MAAEiJlZVnUJNZF8fv8zzphUASQodQQ5EqJYCUEFoo0quoQOidUEVsiLgCK4qINEWQRQEXXJUia0UUC4uCAhZ0gywCyrpxFVFBWXDfGZ33HT+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+QCGQsPKd+v0wMAAIAASURBVHja7J1lYBXH18afvTfuuAZ3d4oUKaW4FCtQJLgHLV4cSnEvTnAt7u7u7u6BkARC9N59P8yzvP+bJkiLpPT8vpzI3buzu7Mzc3QAQRAEQRAEQRAEQRAEQRAEQRAEQRAEQRAEQRAEQRAEQRAEQRAEQRAEQRAEQRAEQRAEQRAEQRAEQRAEQRAEQRAEQRAEQRAEQRAEQRAEQRAEQRAEQRAEQRAEQRAEQRAEQRAEQRCErxHXc0omSaukuW/MnzO7KenyA/9wUe6dIAiCIAiCIAiCIAiC8AnQ5ir5nbOSG/5U8vRzJTs/VdLZZHtc4XtK1urIP/wi91IQBEEQBEEQBEEQvj7McguEL09OOyWHRylpmadkMA1b3/or6bRXSa+XSn4/QclNY5V8UUvupSAIgiAIgiAIgiAIgvAJqM+UwapXlHQ15EMlMxRQsvUmJVfuVrJETbl3giAIgiAIgiAIgiAIwqekqRKel5U014v5Y/FPK9lqkpJlmiip7ZBbKAiCIAiCIAiCIAiCIHxBvHMp2Z2pg2XPKGnykXsjCIIgCIIgCIIgCILwnyRpNyXj7f6y7Uj/TMmeEUoWPs1/dJdnJAiCIAiCIAiCIAiC8J+mwhIla+b5zCdOoERW7kLYcYCSOUrG8vl8SrhalNTmy7MTBEEQBEEQBEEQhK8Xk9wCAcBmJQrfUbK8g5Lag89z+vw8f0VdyfUrlTy/J+bPx2uuZKPbSrq0kkcoCIIgCIIgCIIgCF8vYsASAHjkVrJUPyWLMrIpZblPe96sh5Xs/lLJqB/YLdcp6eSvpHmakvEOKNmKkVmOl5QMTSnPUBAEQRAEQRAEQRC+XsSAJQDI5kZ5TMnUrkrmnPdpzpe0mpJNryq5/4iS/lWUnFhGyfmPlRw9VMk/A5UsfJ2/31TSelWeoSAIgiAIgiAIgiAIwldNP0Y0WQsrqTOVb/wAJe2dP+75avH7azLiStvFf9RTosIgJY//qeSt2UrO81Ey60R5ZoIgCIIgCIIgCIIgCP8JkndQ8jBrShmGK0Nebq9kuksf6YQ9lEjA1EDH5LF8ro4S3juUzLKWn68iz0wQBEEQBEEQBEEQBOE/xU9NlAwroWR0A1ZEFyV9j8u9EgRBEARBEARBEARBED4jHruVXMsUvuiGq+jycC4lU/wu904QBEEQBEEQBEEQBEH4DFR6raT/EyUj7JWMmqSkpaaS4dyd8FU6JVt4yL0TBEEQBEEQBEEQBOFzo8kt+E+xRIn8F5XMXYx/T6ZEDRZRj39HyZkv2E0qKnn9sJJ7XsqtFARBEARBEARBEARBED4jGoumD2ZE1jwT/35O7o0gCIIgCIIgCIIgCF8ak9wCAUAldoe2/H0QpZPcGkEQBEEQBEEQBEEQvjRiwBIEQRAEQRAEQRAEQRDiNGLAEgRBEARBEARBEARBEOI0YsASBEEQBEEQBEEQBEEQ4jRiwBIEQRAEQRAEQRAEQRDiNGLAEgRBEARBEARBEARBEOI0YsASBEEQBEEQBEEQBEEQ4jRiwBIEQRAEQRAEQRAEQRDiNGLAEgRBEARBEARBEARBEOI0YsASBEEQBEEQBEEQBEEQ4jRiwBIEQRAEQRAEQRAEQRDiNGLAEgRBEARBEARBEARBEOI0YsASBEEQBEEQBEEQBEEQ4jRiwBIEQRAEQRAEQRAEQRDiNGLAEgRBEARBEARBEARBEOI0YsASBEEQBEEQBEEQBEEQ4jRiwBIEQRAEQRAEQRAEQRDiNGLAEgRBEARBEARBEARBEOI0YsASBEEQBEEQBEEQBEEQ4jRiwBIEQRAEQRAEQRAEQRDiNGLAEgRBEARBEARBEARBEOI0YsASBEEQBEEQBEEQBEEQ4jRiwBIEQRAEQRAEQRAEQRDiNGLAEgRBEARBEARBEARBEOI0YsASBEEQBEEQBEEQBEEQ4jRiwBIEQRAEQRAEQRAEQRDiNGLAEgRBEARBEARBEARBEOI0YsASBEEQBEEQBEEQBEEQ4jRiwBIEQRAEQRAEQRAEQRDiNGLAEgRBEARBEARBEARBEOI0YsASBEEQBEEQBEEQBEEQ4jRiwBIEQRAEQRAEQRAEQRDiNGLAEgRBEARBEARBEARBEOI0YsASBEEQBEEQBEEQBEEQ4jRiwBIEQRAEQRAEQRAEQRDiNGLAEgRBEARBEARBEARBEOI0YsASBEEQBEEQBEEQBEEQ4jRiwBIEQRAEQRAEQRAEQRDiNGLAEgRBEARBEARBEARBEOI0YsASBEEQBEEQBEEQBEEQ4jRiwBIEQRAEQRAEQRAEQRDiNGLAEgRBEARBEARBEARBEOI0YsASBEEQBEEQBEEQBEEQ4jRiwBIEQRAEQRAEQRAEQRDiNGLAEgRBEARBEARBEARBEOI0YsASBEEQBEEQBEEQBEEQ4jRiwBIEQRAEQRAEQRAEQRDiNGLAEgRBEARBEARBEARBEOI0YsASBEEQBEEQBEEQBEEQ4jRiwBIEQYgTuLxUMsUyJZ0rvueBrkpoqfl7f7mXgiAIgiAIgiB8bYgBSxAE4YuS9jslmzoo2ayxkiNpiEq59e3HJ0qo5MAVSn73Uu6pIAiCIAiCIAhfG2LAEgRB+CIkWKtkzQdK7ghVctV0JcuNUrL60rd/T5EMSnZepGSmtXJvBUEQBEEQBEH42hADliAIwheh+BUlXzBV8JKXkgndlUywQElz17d/T+58Skb+qeT5vXJvBUEQBEEQBEH42hADliAIwmfFnF3JRN2VPHTU9v/lnioZUVDJI9/G/D3uTDUsNFTJO1eVvHFF7rEgCIIgCIIgCF8bdnILBEEQPidWXcmVrFX1YoOSSS1Klv9FyWOXlTyTK+bviTdZyVzllNx3ScnHu+QeC4IgCIIgCILwtSEGLEEQhM+KflHJADfbvxdOq2Q6byVnr1Iy1D/m78nD1EOvI0oeSc/vzyr3WBD+tcsyGqJN4UpG5Pky7dCSK+lUWcnwECWtC+UZCYIgCILwpZAUQkEQhDhBweJUWFnM/aj57Z/Pt5o/bFTiZKTcQ0H4YLS40Qznk0rW3a5kzn1ftj3ma0pW6q1k2UAuG3NLlxEEQRAE4UshBixBEAQbXEoqmYgpfNq1T3s+t4dK5nit5D1GXN11jPnzrj8pmZ81sm6VVfL2EXl2ghAdryFKtmKKbr+5SvZvpORIX75PB79M+5wGKNniMMef0Uqezf5l71uUq5K78yv5fT8lyxWXPiUIgiAIwpdCDFiCIPxHcWympHdeJX9aruRg7gKYy15JrcwnHoaPK+nA4uvhVSkbxPz5DBmUzMeIq1Mcxx/LIxWEv+BQSsnME5Ssy9S4ATRktWuuZMZFX6Z9de4pmfOYkstOKxn5Xdy4f8+eK7kwlZL1arK9E6VvCYIgCILwuREDliAIXwtLldAKcHijQSoRd+krUV7JZkyJmcLIp600WE3OQoVtsZL7WGTdevfTNvvlASUPt1XSlb8n/I0faKdEMra/fRP+v6WSx1kU3hJfuoAgRMefkVW9+L4PZWRlSH2OF/uV1Id93nZl+UPJuiuVnMuapIFecfM+nk7G8WaWkq1o4HduLH1MEARBEITPhRRxFwThK6HgVCXr0MCThrvxZeAufimeKJlglJJGkWSDOVRoJ1FB+1zFk/XflZz5rZIeVKhbsGjylUdKOqZRMgmvL5i7Dx4vyi/aJ31AEP7yfnVX0njd/fvw/WaEo7MR6RT+edpjWqJkfaY23m+q5JH1/477uSYxx0sa9kueUHKzdDVBEARBEAThc6AxlWIoUwPmcWH9qWv/CMLHJCsjkU5y9z5dfz955KmSGb9QSozbFCU9qUA7sMZVhpxKZjmnZNJpSu70U3KvJ4/bIc9eEN6XH4KUDOB7H/pKyZ8CP8/508zkuHNTyare/677Z2LE6gga2Ocx9dJ+qfQtQRAEQRA++UpEboEgCF8Hl7ib2HBGCATVevvng92U/H2Fktc6fN72JmDkxSimNs2frWRihjJcp+HqMg1ZeTcomaEYP8/i7UF35dkLwkfnohLmF5RM2dX+YRHzkjQ865uUPNn6I7XXg3IZZaF3fN5IUW77YaexMrX6CCOvsrE2VgZv6TKCIAiCIHxqJIVQEISvZTjjLn7mNUqGDVLSM5bPL+auWptq/T1F7p+SoZSSP/ZX8sRRJSNZ1B1MZUyeSMn21ZXcOlnJJbw+ZJNnLwgfjOHAG6mEEw1VhRlR9J2XkvFYM8uJKcrPmLq7h6m9xxjJ+Tr5O87H97f4AiVvdVPyecA/uwxnpkbnrqBkqlZKWrmL6UGmHD+co2R6GsZz0WCv8/ovMRX5CncfxI9vP+/pdEq60KCVrQe/R3qWIAiCIAiffAEnCILwLyUla0f9zsir3sFKLvJXcv9R289fpoFrPHcXDE38Zdr97IKSaxIoOYWGKj0jFVIWoe+9TsmrVJT7HlHypRiuBOHvQ4OUmYafVoyI7LlTyXCOEzsqKnmMBqDSlZX8c5KSo/g9Kbq//XQJLEpmoYHoFiOmQqv9vea7M9X/59X8XqYk3mUqcrrdSo7m7oaV+fkfWevv5W0lv8/K8Sezkmk3vd/5H7Coe0A9JfOfVlKrKV1LEARBEARB+IRIDSzh34SJClMJpv5tZf/dx9osldrwcwmVLNWZih1T9HrVjiMXkkSJPIyYaM8aV78wtbFrHiWrNFPS5ao8e0H4p7ypgcXfjVp4mx8rmb7s24/PxvfyUgPb42fQ0Ow2IubjshZU8jzHq1YV/+Z8XUfJhg5K1l6rpLmo7eeyb1TyVm4l90zg308rmYLtvcTrvsXi8jlD368d9reUXE6D2ELuqmjfS/qYIAiCIAifCkkhFAThX4LnfCpujFRqyxose2YoOWIYFTGmyOCZEvsZGdGmrpLn4kqtFqb+nKbh6vJoJV14HSGM0Arn/42UQkEQ/jm6xfb3HYxcurHt7cddZA2sXTRcZeHf64xXcmUKJaMHMjlX57LrEIenLX+v3d4snp6G48UEGugtN2w/Z89xw5Pnu+Kj5KWGHGdY02s82xXUV8mrJ96vHVbep4D7SqYroKT5tJKR0sUEQRAEQfjoiAFLEIQ4Tk5GKv16XckMNOyMZjHzJaxJExJLEfYoprpsiOPXGdaVUh65IHx29Bbv+UGmDt5kzSvrGSU99ilZmjWnNqfl9zJSyS0DP7dQyUhK1P3AhjLlcH0NJYN+i/ljWWhIs+f/919he+Mpyc0XMTVaijXyvmc7GMFlNSuZhJtPmIOlLwmCIAiC8KmQGliCIMQxHKOUrMft5ufQYGVOqWST75Wc7atkyDS5Z4IgfF4imWqv17f9e4ofuLyKlqrs0Ifj21P+4fjfO+/ds0qeYu0uPIz5c/lo6Arx4udXfeQbYE/J1EgTU7jxvfQNQRAEQRA+FRKBJQhCHCFlGiV9OypZlQb2eUyBmUKDVWBuuVeCIHxhjO32Utj+2fIrf3Cw/XsQI0Gfs1aUqcPfPO+kt//beauS+ZjqeJ2pgvfyfeTrZy0xvaqShqFMbyldQxAEQRCET4VEYAmC8KWGn1NKlmFNlwXTlcyfVMnOzZUcuUxJMVwJghBXcGQEkma2/ftVphBaLtn+PZwpdjo/H2/Cp2lXZqYyZr6j5HHWpnoZy/jp3U7JDMk+7Dwadz10HKPkQ25GEfW79A1BEARBED6ZBim3QBCEz0u8QUp2ZdHj8Sz6dIa7X/kwAmsTFa7Il3LPBEH4tDi/p4HcgbWs8v3JZVR7Je+ydtQmS8zHBc5T8lV2JZN894ENvM3zMmWwPb8vub3tx3JvVzJxZo6rTGW0xI/2fXQclGcReO8P3NzCzHE8ZTYlHzP123JK+pIgCIIgCJ8KMWAJgvCp2alEXkZW/UEF7mfWZBlOha87I67uecgtEwTh06IVVVIfqGR1Flsvw938TD35QRZrN59UshoNWN9xV9SQx0qO56YSp47GfL7HTIV+OFTJdPyc2eH92htvnJLDWEvrd9bayr5ISScaqkrQUPWS42pgLLsK5i2kpHtlJY1IrffFI5WSKV2UPMddES2+0rcEQRAEQfhUSA0sQRA+EU7cJas2IwW6sjbL5URKNuNuVacaKWn1l3smCMKnRWutZCQN6xO4W+B5GmBa02BVaY+ST5hal4kRo0WLKHmD/5+yUcnlxi5/sRRVj+A4eKi8kt9zV9X4u5T0L/b2dkdy98DnHD9Psni6ToNYXUZU7fqJn6dhqgy3X72ZVckUbZXM+UDJ1Yx8fTnlw+5j9tO83t5Knu4ofUsQBEEQhE+NGLAEQfjIpOmmZFem1JR2VXLeRCWnUVEM0nnAZrlngiB8Hk7SkFR/tpKnGEEVwIj0bb8pmZ+RWWlZ0+ooI0gX5lTy0lwlH9XkF+96v/Nvn6VkAyvHy4xKvst+/6q/kr0TKFmkKY8vqeQVFlU/xBpbm9MrWZbHl+Curo/KKTn/iZIPD/+9+5iX4/sDjt8XJXVQEARBEARB+BxoXIgP5UJ8HrcH167JvRHejV1+KkqUuxiJsDNQyTKMxDL9KPdKEIT/Nm6MPF3CGlq9l//Nebu7kmYWU0fRWD7IFEO7yfxd/2ft9+yq5CpGfrXrIs9UEARBEARB+IyIAUv4OyRgRECP9UpeilByHLdr92aqC3rLvRIEQfhfKjESa8MjJVMP+ne0u3w1JTcygsvbXZ6lIAiCIAifCyniLgjCB5JvgZIzWUOmJmu3DApQsie3h79n5MQMk3smCILwv2x3U/Lq70rWMIqo346b7U3IlMq6rOE1l+2/J7vECoIgCILw2ZAaWIIgvAOXMkrWzqVkR9Y6ucRUGB/urnXxCA8YKPdMEAThbYSz2PpYeyW7c1fBb5lSuC+OtNP8VMmGHP+vpVFy5UZ5hoIgCIIgfG7EgCUIQiykZYRV13Alv2fx4BmXlJzJXbCCKsq9EgRB+DvcjVRyDmsE1uqr5NOWSl6Z/mXapbFd5Vjk3SWTklM7KxkpNQ0FQRAEQfjsiAFLEARiWkGF5YCSfbj9evASJdtWUHJvGiWj+so9EwRB+BicoKPgOVMK9R38xxcyYJlYnP1moJKHQpV8IYYrQRAEQRAE4Uvynyni7qmEOxUE09i42UwzFRe3C3wOnzjCKXFCJfulVvLiMiVHsshw0qkxH+fUQ0mXOJ5K4tpWScfDcfT9K8Tn/T2ff5842s4TvJ/P+fvCuNlOJ0YGur3iHxzjZjudW7F58+Poi8N2uY1R0i6O7jZnfsxx/SL7pVfcbKcDa1y5c5fWL2WYeme/zMB+eSNuj+suHH/szsTRftmIyw6+R6aacbSd59jOePxDRBztl9z90qVc3O6XzlwP2Q+Io/P4Oj7vFnz+HeJmO+2dOF4m4x/Kx812OnIcd94HQRA+G1LEXfgvEaREl2dKdv1FyZT1ObGfihvNTFRdyQGblGwcomT8h2znP3xvtfZKFnRVcgZTA6vQYDaAimBvnudx65i/x6g5PIkL4Ao0fLpMjVuPvRojyUasUbJoEy6Q4kgEqv1JJVtwF7Ke3J4+LRdupmpxRGFMx37BmmfNmOKUiIYiLU/caGdeFyUn0HBZgwqPeyV+oELcaOd3uZUcRYNbIRo2HB/EkeUB72PTUuyXo5XMoFHxiSOGwYR7lfxtm5Kt7imZhAoa5saNdmZke8ZXVbIO3xfPnPxAm7jRzhJMHRzho2SxF0o6LY5b43rdzUr23a5klqvsl3HEUBS/hpJD+F53nKFkCgvHyzRxo51JWcvyDwclf06hZDxGZGN23Ghn8aRK/r6Kv3M95Jg9bo3rVVliYdgfSmalI9AcR+6j834lf+Uu0l1Y2y4Vx3PtYNxoZ3K+NyO4CVBDvjcJ6GhFs7jRzjwcf37ne1SK63fnKsbAL6qXIAjCJ+E/E4FFFjPCyFpPyVOccAaMoyIZRYU925dpXzIuzA7/pmQYI532ctenDpk5L77mwug9DRwuTA1s0VvJc9z+fHZxLrR68bnnfb/vK8saWS/vKxlABXf5ZC6EWYzYm/cVm7/M/fRNpWQk+/fD8UpOq6NkOSo+8Qp+mfaxhjMm01CoJ1byLA1aQ2h4Lcp+6Xj5y7TTvaGSO7KyX9IQejCLkp1pKMzC8cPc+Mu0szSfoz/f82C2b21RLoRHccFuGIITfpl2NuACXOfC/D4N636MyKvESLyEPYyB+jPPC+eVHMNx0sLd5i5yXBw9UslvqYg7f6FImKSMUD3O+xjB8eYA35uuzZXMRUOR3aIv0868NAw9ouPgFRXbLXzvm9HAmoqefK3ul2ln7Z5Kht9R8omxLuC8UZW70CZmRDMafZl2DjjNcZ2G4Kt0wIzl/S1Bg5vLhC/TviSDldxPg0sU58dj7Ie/0CGQj+OQ/Rd6f9LQsPaU88srvjfb1irZkobX9Gyv6QuN6z9yk4NgjtdPkyi51KxkLRpak+79suvLTnzeuq7kdY6fYzm+f8f++CZC+DPjmpzLMUZSWxhBdIrzYj/qA/kMh0r3L9POtIxku7JHyVAaBHcy4rZdCSUzcBwwrf0y7SzFdeRzjpNPVyq5jOvon24pmdwYJ0eKzikIH2GFKrdA+H8D1hBOAN4ceBtzlyQ9jnkQHBvQ4GB4rvX3O85KQ8Y0LnirRzNQ6ZwIn3CCOVBSyY1ccBzkQulhNxpumvK4qFgWCjQ0uLOGFJze0b4cSib7TsnZVHDytbf9nIULjWtUzPaHKbmJisdJeoKe5Fcyxc9cMA+gIYQh4zNqKflnPioCVEhMjPgCFWlEuz7rLi7Edis5iwY0l6W2n4vgAuOSs5Lb6TnbxgXUeZ73KRfMkbEoQua7Ssan4mnmQho3334/LfQ0tuR1D0hPBdbL9nNBNGidYGTRTi5Ad3FBco0e1eeMnLAOivl8DvRcxtvK9yoR/xHLNvM6IwCN1KJBw5X0+SHaB2fyeVJR30/Dx1ZGbu1j+x8aBhu2V4/F0OUan/2Shh1wgY2QWNpZnZ9fzX5JBebbpNH6xRM+Fnp297ejgs7rPMH2PqQhNYSeX4RFO2E/JTxooHUxUtfmveN5c1wowZSd6RzP4kfzKEfw/b6YVsndHBe2cxfNM1yAPuX9i3gc8/lMPD6BkYphhCSefHs7I2k4rsuUiEnpYv7cS/ajU8X4fvN7D/J+XmaE5nOO05ZtsfRLfr8XxwGTYai9GEsD+T6a+P39+Z61DIr5409pwDpERWfLn+yXVDDv0RATREM5EsT8Pc4c/zyN1NRNb7+POiMqk9BgPp/jc65c0T7H8fSuv5J7mQq7lePDEUbqPeL48sowJDSNboHie8B5wbUU//4Og7KF58tTmO3k+ZI4R+sXxmPhcz/A8XwzHRqnGWHy5Bv2y1giY02MKI7HCEl7Ph8secf8Q4WxOt+7yTRY2/nafu7VI74nnB+2MMV0J+fJ6zR8PWNkrsUr5vPZM1UxXmm2mynUuPWOdlLh7s73p+usmMcRf743R2kA3sr27mUNr7uTlHzBSFw9T8znc+J5vOjogHHfrbG0j/NqEo77s3hdBZ9F65dsz30aLvezf22l4esAx/PH3AXypZHaVyLm83ryPXXhdei933Ef+VzSMPJvI1NGE6SN1i/5Pl33YDv53DYz8uQk2/mI/SK8Viz9kuOOF2tnOjDSC9vf8f7w+6pwfBjP87lVsf3c64lKnuV7tZnrwu1s7w2+YP4p+b3ZY1lvcB5JwHHLxHEDT94xrnNd2Y7rw4HnY/7cC37/Ud6vLVyX7J2j5G3e/xd8f62/xLIO5vvsxfWg1pH/iC2Clw4RF+4SPcuPBpjojmqWsLjHCPDDnBc3s/8d4vx+n47Ll0bEY3QDEudF93EcLwP493eUHtA5Dqfh+nBRoJLpVtl+LortvMpxeB8NXLvoWD3KfvuI1xl2Jxb9h+9TPL5nDlxvYf/b2xnFcaAE571ZfO5e0fSlMO7ifY6Gru00sG/l/b3C5/+M/SLywNeucU6YMGHChP93K3B+TM7nW4iO9GxcZ2p8buHUD07/ymmX43JurnuKcX1h4Ty7yNgtfbmvr6+vr69o+mLAEsSA9cUpQo/lUEMBeP1+x+lcyOZgaHzi940wo0JxkwrfQSrSG45y4ceFcvBR28MacWHemAOzluDtC1+dCwxHRljlpOLi8Z4esEhO3Jc4kV5hREIaXmcg79dAptAd4EKnHxcK3/eI0V7113Zyoo3PBWNOnscc//3aGU4P6gka5rZxAvqTC/nzjNTSGYngzRDy4VywJud91Z+/o508jzcNlRnYXi3n+7UzaAufNxWqLVxAL59CQ0y0lJqCXFgP5ILNyTCQzoqlffRwmji1ZuTCOsWM9+yXVNjvUFE6yNTDDVzYbaSi9qKn7WG1uNBrxf5iZuoIYolM0Wk4sWO7ctIAGe89a3VZaVC6xhpku6nIrGQ/2kXFKJJFou2ogHWgYaAaF3TWAu943qVs25WDCpl9y/drZwQNeqdo8NlBhWIpFYQL7D8WFtv2omF0CtuV7DTb8es77gfvfzIqUNmKf9j494rP9yDv62amKq2mInhrYjS7Cw12/Xg9Xs62z/Uv95EDgEYFOAPfa+/39cTzfb1PheoA+9kW9sc1HAcDntoeVo4pi70Y4Wd9hyFQ5zhs1CDJRcXKvcx7zgdUMK+MYDupSKxiqs8OGpTDjBRzRna0Ybtqc3zQ3xFZqtOQ504DUe7J0RSld43rGW0X7jto2FvF53mMCrZOQ2wCGiaGMeInMxV5a4t3tJPzfWKuB7LzuWnvWXsmhJFZh2gI38bx708aKm6kt/18Zr5nQ/keJaDBW39HZIrOfpSWDpO0nJ/xnrUin9Cwdojv7VYfjuu8zmc1bD9fmo6i3hzvzUbkZizt1Bk57cj+mIMGA88e79kvaWi+TgPHIYbo/knD6jYa0kI9bY/rTUN7ea4fooLfcR7WlHLmdRTg87Hv9n7tjKJB6Dzfl21UNNfyuRsOPyvfn/jsxwP4vuakoqm7vqOdNDAk4jyabRqfw4j3a+dLvqcn+D5s5/NdyvfwulGLjhG6yWnQnExHkKexTqz0jnGdhpJUvN6MH1iiwJ+lFw5yfNzC57eS9+dJPdvPf8txvQefu4uxucK3sdxHI8WW43tOXnf8Me/ZL2/TwMbr2sv3fA2fxxY+z9c0/JlpGG1LQ0JNOqas77iPOsc1Z84HeWi4cA59z/UG5/1zXNfs7MT3h/fvCB1IUXzfjVqUg+hgLsD1gp7jHeehQS8+1xs5uZ62W/B+7XzNefI459Vdx9lOrovPnfi6DVca1+flOM8MZP9KQEPxEhoI73Kc+4b9ojLns/ucPy/QMGitTr2L+tYUOurbBYsB6+tEdiEU/oU84EJyyShjhf+eExsV+Zb0yCd+lwGDkT9nuZDZy/Pto4J4kRE2r42FQzQD1nm2a6kREWDIWCLGrJzg4zHCKTm/zyO2iZ4L5Hv0FN3jxJyWBp6UVCCnMUJnAhcOj40Jlh68g0y9eGpMBFfecR/pEc3BhWDWeTb2lL9g4QR0g5FFh7mw2L9ayZM0RN3hQkg3IpBowAqmB3I9DW/ujJDBO6YkKxds319XMj0XCloskR3hXBBeMtrJBe5+GiyOMeIj4EjMxz/i55ZzQrY32pc7lufHBZAdF2Z1aIBK8a4FORWas1QY9tOgtJf37xwNuq9iiZi54s1+SQOFafw72snIAWcaThJT8YsXWwNpcHrA9+YUDXsH6HE8xAXuJUa4GKk1YGSKlQbeY9xk4TUNuXrudzxv/j8rDUQZqUjZx9Yvef23aKA5QsV9H38/yQdxmwt9i+GBp4EgjIrvOhrmPBhZEFskxxvFj5F8xZa9nwErgvfzCseZwzQY72QkzikuvJ/E8ryf0EC8ihEaTjS0IbYUNRqQNBoYatFA7/2u8ZLfe5EGsoN0hOyhPEMD9EtG9iBaraJbNOgt5vmN1MpYFR2OH568j9513mHAouH6MRe+Jzl+7qbCcILFjS/SABARbX2kU8E/yfctqkrM1/GXfsmFeZqySmYyDFixfZ5L/HtMxT2ywnYcOsj+eMOIHIq2iUYo+8OWu7btfef7w/G8IO9DNiqOWiyRmZGcr67x/wf5XHfzwk5ynHxcJ+bjA6igrqYh24UGfLxjswArx+dq7C9p873dgPWSBoCLVOyP0BG1m+PaWZ4vOJYaX3c47y/l/TMZ82d4LP2Sn/dk5FXS8e8wYDHV3p+OkeM0pO/j+36MitkFvlcRsTgUD/M+BDAi2nr2HfeR152I/TPn4HcYsNiv7tBge6I25x3OX4f5flzjdVuTReuXfM+20xB3YfB79ks+53z8vow/vX29EcV+fJ2G3cM0TO0xIp+OKfmQqWmIZrA0It1W8Xk4Wd4+P745L7+/HB0JGX94++df8XkaDrvjnBd20+B9hjWpAmMxIN/LpOQKzkMOq98xj/M6HHl9CTnfxup3NHZFZTtOGRHAXNceGsd2MuI2PFq/sfJ+nEjDfhv4fs9b5/uQiO9zOhqknGP7PPvRA/bPEyxlsIcBGof5/l2hA84SbUEdxfFuN8ef6zver50Wvl9ZOV9kPvx2tdpCQ/gtGmaOsl27+P4d43XfpuE61gv+SsjGSMNRXAdm53vZm+/r71z/W/n7Ui7kElJfqcLndp2GrT/4vFy5Lj5m1OJtCeGrRCKwBPz7IrD+KYu5oK5LD24EDVHXdtkusFbTwnWBCwt/TvgY82nbl4wLmpVU5L6hZ/11ESoGVHDWUbEJoaJdixO3PRcWYxlBsI4L4ojeH7edZTkxrKIhzJUK9HMu8I9zYlnOiWYPU7nuUdEPP/J5nncHKrhjqDCZaIB40Mp2Yfsnx8OTu/l/Ltwsaz5t++y4YBlHT287KlpRvE/XqfBuZQj+anroL3Kl9JTXo0/8tO00UqdW0VBahs/1NQ1n52m420DFdgMV+ys0SL36TDX1SrGI93K+Lwn5eyANrSdpWFzK938PDZe3qWCF//h52tmABrb5XJDpvK8PuMA+xPu5nAaoozQk3edxFqdPPC/QgD6KHv4uNNhYGAl0h+/9JhpKViayNWA9Yf/Qs37adialBraeC+H8fJ5hxvjNCI71HL830kN7hfOd4aj41OSl4XwDn2cyGm6DGPl4guPAMqbC7aej4wZTKsI+k8OxNvvhIsPAznY/YgrjAUYKrOZuvgc5H97nfBTZ9vO0cwD7WX8q8hZGfNymIrqd7/MqOjLOGOMlFUjrJ94FMgkjcVZwPizO+SWCBoBLHJc20uC6jhE3F+iQCqbihsBP2840HL+PGpFOXC8EcfOL07yvqzlfb2fq2A0qiKFun+d5V+d55zFyzCjN4M/3/yDXIyu5fjrAyDAjgiMyw+dpZyf2r7F0pFgY4XaXhq2dfI/+ZDtPMRL0KQ1q1m2ftn2uNDD/SY9EOTr4Irl+uMJ5ZxPbudbol3SkvDD0gU9cyyktDWKbaSjKxEi6VzSYnmJE3jquNzazf9zivPjK/Hmed0lGEq/i+xqP65xnXG8cZbuX01F7iHrVbY774dv/K5qmbQRWW67Lx/varoOr8D1dH8suuINoSP2Veup9Gmar05Fylu+Thf+3ukgE1teJRGAJ/0HCmUJwmoaBmfQkbqAh4A4nPj26B23MZ2ogFdWQ1UruYqjtNCqCR/yUrECPZAt6PE9wwTySC/qrLz5tM3UubB8btWUYmbOIBjXDI/zK9cs+bysXZDeoGPy5nooFDVXn6OmJqvFl2qcxtSeSKWEXGGk2m5Fuq6kQ3mbEiPVL7QbG1L/XeWj4oeIym57ELfRwPvH8wi84DS/PuLA0UikXUqE5SoPRXzzbOz5vMy1MLbpHA/ki3r/lfO7n+R6FP/1C/ZLjShT7nVG8fT4X6qupQFxlSoo1euRT1s/UUBoiQ/h+72Utu3k0aG3kwviRsc35pC/UL0+zX3bmuE7D4GK2/wANfi+iR+j++pnHS3qw71Hx2sCl/wKm6p1mvwg3IhIXfJnbGUUF8CLnlyU08K9i6s1FKsBWhy/0vFmDMrQTny+f6wK+7xv4//uUOiNZ8Lk3j2FEYhAjjnYxgmwBFcF9fH8Ch0c77jM/d50OqEfsd4v43iwwDMBUjUPP2N7/z94vabi6TQP6XCrkf/L5XmYt0MgnX6hfcv4O5TpoHzMTFhqpi3Tw3mEKnZ7rC7WTtUaDaDhbTweaH9/nXYwQDjDWuV9qnUnD32MadFfSwbjIiOSn4/i1kZK45aOdeOLEiRMnAtCgQbPPijzIgzyek7SSWkmtZNAqOMMZzpGT2r9q/6r9F9oyIHYc+X5oJaO9528izwzDk63hy0LHLuiwt6PB0I76R2SIGKz+G4gBS/gPMoWpN48DuVCnwqt7x432vaRnvj89tRep0ManR6I3FYniVDQmsN0LmCLz6sXnaedlphr+TEX1LBXI0DiWu2+Ehu+pzHbTsBYVEDfaZ+wCtLAU9Wve11uMDDEMcF+aCCoKvzEi8Brb/cwojuoZN9p5nZFXzfmcT1Znv2Qq1V+Kxn8hjlBxrMpIgvOMGImqFzfap9PQs4wROdOpsN6ix95aMo6Ml+yP/RlZdYGpuP5sJ9LEjXY+oCGyNVM3TzPi4vVaxClOMje4NhWC86xFF14ybrVzPefH5RzHjd3eLHrcaN8r1kwbymLEl1i03d/Y/GVr3GinkXLYguPkcdYsesXUbbSNG+08x9S6+rfZL+nICa8bt/rldjo6dzPl6RLfe0vyODKPc9wZzUjfazQAPDFqb3aPI/2ShqEuNKyep4Et0Mgcmh432nljnJJNaBA8w3WHUWLgMxCJSESmbYj6qI/641rrGfWMesYTJfEDfsAPCzNO2DRh04RNV6ugAzqgg3VM3DDwHPZT8iFThA31y9twJJ+3NVzZM8MjYz/b77nA77luzPfJRMf9byAphAL+eymE/xbMLAZbhYpENyNFjylmv7GI6VFGXlnzyz0TBEEQBEEQhK+XidMmTps4DUAoQhGaNYeeWc+sZ14zH6/wCq8y5oEFFlhu5METPMGTjf4qUms5DZanWQv45ZHPZdCyNUiZaaDsz9pgXfMoeZOO0NE0AN+hQfVbpjb7MvPjLlNbu7PW7/aGEnn130IisIQ4iBNz2FOzlkAEU4AimWKBUl/JhRoGZEZWGdu9ezDVqRFDzRvQ07idIfMTWfzwMSN1krOWhrFbIiKkDwmCILx13DV2r91u+3d7RnrYs4hxKFOC9IVy6wRBEISPgLFuX2b2NnubvbXI7598/+T7Jx6VivkU8ynmo2dNsCvBrgS7IkpZR1lHWUdpDf/yDc5whrP1JlZhFVal7o90SId0Dl1xD/dwDztghhnm9KeV4Yqb8qEhIy6PsFTKvkYTxkwYM2HMxuzadG26Nv1hMphggknz01voLfQWdg3V99w3dr2q9HEMRc5M+fRmiZRfjQjF6kp8w7/nY0S1HVPu+9KAtcVLyVsDpSv9NxEDlhAHsJ8PANrJvPEBQMvUpDoAWN1+4BAZzOLBb3aFivqXXzBrOejGbn+sieLKoplerDURnwN4MCPicrJ45WRGXpkYWg1je+GN0pcEQRBiXjFz3OVC2cpNGhzTKOnA1Bo3zi8W7n76nKnZ2C63UBAEQfgATJSMGLJyHtJWK+l0DamRGqkdmrrdd7vvdj/hec+fPH/y/Mlib71pvWm9GZkYTdEUTbWDMXyzCSZ9q35Jv6Rfsr+EZViGZQm7wg52MWr3OnToXmtUZFY5qM+VrQ0PeMCjy2j9D/0P/Y+bTM22zsQZnMGZtSxlMfQjb8KTiDUKvzV2t0yq5Gye35GOfY0ZJ1bWEgvnpiLWBhJx9d9GDFjCl+h2LAKbl0WUG/oDQIKaP84HgKRbUkYAwJWbd/oCQGS68RzoWCsc2r+t3xoTGIvHWxlJ5dVMyTqMLPuWHv71LLe4ljUzwlgDSzNq4xjbP2eTviQIgvC2cVfn7rJWFtf1ZE2nEqydVJmGKu+etocve6DkfDoSICnagiAIwtswdj9k8Xkr5yHNR8k0nEcq1gEAh4KlElimWqZmL+Ly6kT9E/UvrdfjJ52UdFKSXKPz57DPYZ8z+8WUUW5RblFu5il/OZMVVlj1elolrZJWKfksfbA+WB/c/k9cxEVcTFxRGbhAwxVeKEOV1QPf4Bt8E1YBR3AER+xmIB/yIZ+nPw7hEA4lWokwhCHMsac63oXF0lHh496mcO5S/pq72fbj5kXfctOYMEO/YTF/jfrSSxbDf1RUpSRe4aY8x4rz7zd53CUxcH3llgS5BcKnx/6Zknl7KNloDAC43S13AwB+6JnBCQCaTY8yA8DNk1fSAkCnVgEs5redRdfvZ/067kdm5nj3oiEqFQfsXjRwLSjMgd0r2oHToklBEAThrbiyuGtJpgLWYC5FRu5mFEwDl96IC2UWr77G4vrrGfGKp3IvBUEQhPeAuyWmYqRv7ZcJVydcnaBxfY+oJFFJPO3ibQ4sElik63cRG5EDOUpOcsFt3EYVF+vuyWcnn51yd4MyGN3d3SW2E0SrgVVVGbQaDlAyMSO0wugwfzIZ2ZEd2V+eUimFif9Qn3s2HLuxG7tXrVGfW9VI/b/bSPW7iX+H9ePenifMNFnF3eD73lKyYWy7Qa+M9jsNW2HzlDzG3T77M0JrlxNrbvUVQ9bXiRiwhE8xbnPAzVVOSR+mXpSOBwBOnW78AQDTguxuAEDVVrgPAG6uc14BQPbne0MAwHKpFS3umuu/+3440fBWOYWS3Xh/nnIXleZMITzJ7auNbY4FQRCED8NtnJIlyirZtI+S2egg2H1EydFMFbzH2hqTTytZnd+jcTc2uMo9FQRBEGJH4y6OqYYoWcdPyVqMxAod/F3J70qWXhyUKH3b9G3L/+Sdbvum7ZtKbPJyvFriaglTg/3dg5sFNwuonLG4vkRfAmiGw/5hrKdUkVVACEIQoi9CDuRAjpCGKIACKHDhOVZiJVY+W6Aiqjw24ziO47hdH7zGa7ye8VgZuBbnQUIkRMILdRGIQAQ6MYLJajjKjV3/pvzdO2NbvN0gb2bOy3QgPeLEe5G7nr7mrscw5HElzEzFSVFbyfQrlPyWm5AN5+6a1VnL8lFf6ZtfJ2LAEv42nsZmpzvNvQAgaHO+SgCAbxrtU//4ro6SV3cq2T0BABw6azkOAHnuzlcKRcULxdT/ezgAwMUfi+VWv5t+4Bn2/DvvUDKmBHbj9roVGSq7iBFY02jge2qEHH8rvUoQBOFDcDmrZOlMSjbm30s3UDLhcq6DKyvZi8Xagzgu47AS7bkrUjzOVyYHubeCIAhCDLCYuDfnk1qscVWH+k9YiJKzOY8svVR0QNEBpdK1zW+2M9u9flU64CeXn1x2pAPC9XA9omyiA/sb729cZKVTbvsm9k0crngXnfXNrG9mJbVbtnju4rmL50Yl+EsL5mM+5gN4iId4aDmDtmiLthcyYi/2Ym/2uYhABCK8WFNqy0+IQhSi5rnDHvawv5JRRVpFVVeGqzffuu3T3rYkT5T83dhdkPNvcxr6Vt9WUqfjCdzUC6xFqTFl0JGZKqVYxH0MHU25WRO4MEuwrJae+pUiBizhb2BvAoDwBCol0Fq/yQsAMKUqUREArLhKj3YfDjTb/QBA13MuU7/7cbc8UwklB88DgBTjXtDir+X6d94XUzwlSzAXu7eFA22Ykl0ZSryVnv2oUtKXBEEQPgT3TkoWZ40MHyoSWen53sd5ZyJTCbqXVLLuUSXv0KM9JB8VDS6YHx5Tci+363ZgagMSyD0XBEH4T3MjXd50edPlTR2c+Gjio4mP1jz8ZP2T9U/WNzgfVTiqcFThZFW8RnuN9hq95kGeTHky5dHXti3rXtY9v//tgRs7bOyQv3rhR+dCz4UGV2+QPU+/PP2uzwYcCzoWjJyipKZl7lIhbYW0FXMgQq+n19MfDmuxeN3idYuP1RpNvYEO7gOMPHoZqP+k/6T/VLC4isRqvBdzMRdzk+jq/1uyqwirhcnU75cYoRTl79vRt6NvR+w1LizmCKlPhQdTBJPTgBVFh/6VsUoGr2TK38nYvoHtpeFqT00l7xdQMikNhk4SOf2VIwYs4X/hwPdml78/AcC+koPyZGfPFQUApsnNxgMAfvw+AgAiD188DgB6wz4cKLcfBIAj2V8tAIBC54vQoz3/uZLx6Pluy5DQpasA4OmbCiP6MEo3/iGOe8LjsbaXD3fJaNaCihRDX4e3p+KURLqYIAjCe5FdCTfKEi2VbE7DVFmOu/f8lGxDT/h+OkqsjGx9wQX/b/Qs94qvpIkKwUAWkw1nJLBm7DY4QR6BIAjCfxGNm02lZO2lWtrLRS8XWVI0KFNrba21pTwTDXS54XLD6Zeonq45XXO65nzQ3qmpU1OnflmTm34y/YTymXvfPn77uKt3ijUpbqe4/dApeVLX4a7Dw2677sNN3NS6APogfRCG4v9TAXMiJ3LCQXuuPdeep/DGMzzDsxR+yhBVnpFIz2i4CeSmI9a5KpJqw3dKLqCD/DwDBPRZcasG1D3qiavpgOrIzVS631ayX11loLrJ+T6Ku9DrnPftKONTz2pCvStbGSW3VVRybxXpw183YsASAICeZ/2QkmYu5PM3AYDIDfVvAUCiKT83BIAc+80aAAR49kgMAGcGLTsAAP1nvxwFAL3pt3bIXLia+snvPAccKg5tXyu5dDHP7wQAUb+pX3JkizcXAG50d+oFAKE+6BLH7tdMtpQt7r1byZz3lRyzWcklHKhfV5UuJgiC8D64MJKqFA1PPnm4QOUuRXc5w5hLK5mOjo7iLDZ7gIYpnR7mKaf5eSokwwYo+ct8JS3MORwyng3YKs9AEAThPwUd7SmVAShhrTAAcAhVRcXzP8nbJKBBQIO0oxeOetzrca+QR/3GlS5ZumSm6y9rhQSGBL6OMGVIcSrFqcfXLeN/dfjVodmhb3LcrXC3QuK0M6dZbltumyJD49faVmvbHr+L/olLJC4R2LHgNL2SXklzArQJ2gQ9JUZHNYtqZq2Ky6Z9pn2mMverak20Jlo9+wt4iqd46kQ9I2yRkutYxGW+u5KXWXQ9KpwGq21x5cYaBjRlmArjvP37OCUftFOyxjWqV6l5PSzy/oK1snRuduVGvSoTI6xcOH9P5i7tMzjfP3wsXfrrRgxYAgB7Rlx50INd6kcA8AotFggAJe+mzgcAdZM+ywYAWYNbOQJAUvOfjgBwn+pE/ibG931DA44fPQaJqFC0PajkEmNXC6f/bYWuF8gKAAv/6NURADoOsp4DgFDorGWCFV/2PjnMVbIWDXyd01KhYpHDphxIjzMV8k0EmSAIghAjbqxtUZye2CbctTZ7UyUP3FCyBR0E5zgftLmi5CAuiPvVUtLkpeTwjkpG0jA1iZFcuo+Sv1Eh6MGaHBrnJ4fVShqGLWySZyQIgvBVQT0lFeebH9u5j3Uf6xHSyDlZxWQVs87KMf9B+Qfla5scsoXcCrlVodzsa5YTlhPZtnXpMq/WvFrh44KC52EenE6oTT50o9Q4CirHiT6SBibnSUq29Mp3JN+Ra3PTz7Qusi7SchacpvlqvrgVdD1gYsBExwahLZ4NfDbwepqkqVIEpAhI0NC+v+t21+1um59M0XPpufRcW5i5Mp+7sl9iraeo1v+WXfZsDVmBnTgvM5V/6Rglk5bic0mjpDsNdDodS5Gs6TWd+upTGrge8/ORq2TXwf8GYsD6T2L/p5K5uA2pz1Qla8QDgEL5k+QHgF5/mAMAoPQwfxMAePq3+UXpD396/u+3JSls/FSEuzvNu6xkAioIbekJX9Iv5vYUpMIwrw0AeP6SpRgARPQ9ozzwvXDky96vlCy+3pW1Vn5IxImEE8o0KkYvfpC+JQiC8DbcuCvQt3R5+HA+Ks6FbHKmAlykY2XcLC7YW9t+z2gWXddZs2oQa3z0a6ik6YGSv1FGcBfYCRzP9Q78Pz233WlAe8zdb+c/k2clCILwVXBVCW9vJWsz4qfmCQBImMBjjdtat7VdvNJUqTWx1sQkaxyyvS7xukS8dq9/cG3h2uJclpNVs2TIkmHtDvs1j6c/nu5ZWtvbe0XvFZEl9P6RXpFeQOEt1GNY8zYxU/3atu2ao2uO5XbLNE3XdH1c24rIgzw4sGu3dZh1mPbD01DLSstKx6oll7ukcknlcgiBkX0i+0QGr+5qHmgeaO4/qWBUrqhcUbkusoaudSwNNGP/rQ8imoGJqY7cXPCNPBf9sITShYX/xSS34L+APYvQ5mMNppEMlfUbCgDObhnLAcCPk64/B4CFk5Xhqnpn/0wA4Onftqj6/J/LY/7+wqn4faxileCFku3oEV/yMubj8rsoOZcKROJsAHDp0NKtAPCy10sf9Xft6ee9X3YZlCzDoux+VLgK0mPTiamWI4OUfDFd+pggCEJMuDBVvBxrG07lvDCKKYIvaMDqzwjXS3eVzEYH22Q6BjI2V9J8hgt5GqDG8u8DmIpgyahkH6Yi9GDRWHv+XR+l5CR6cPsYxW05z3gPZsP3yrMTBEH4N6LlUDJlWyV9ubnHQmZyVJkDAGlvri8FAMd/vl/aOsQ6pHUxVwfHmY4zTeesoRkzZcx0t4TLq6Rbk259emjUnhdXXlxxfbouyMXFxSUif8/tNY7UOLIvQ/3njpcdL0fm8KuKgziI5gm5q16bRbqu65q2yDVN3TR1H13xSGwqbypvrWl/B5GIhKdjxueNnzf2LFdx4lnHs46Jb8a7vCDJgiSrakc2H7xx8MYBaZb6d1jTYU2HNeepH1l3S2SRIPyPpi634GvERINRHhqqGnHIK8OB+8YZAMhsHfMzAOz+uVIwACTaWCQIAMy5nivDZlDb0wAAzxV+/OJ2tucpxNRAPxq4knDh3zZQySX0dGCH7XF5qbD48XMpaWpvkw0A+g57xfP8YpzX6/PcN48Itp+51Q0ZMbaVnpsJ3Ib9lpFaMlH6miAIwv/izJS/khzHm6RX8gemVFg5H7Th7korEysZxUjWC/S0zmXqQGl6sncyhWACDUujOS5bMik5lkXdTUzZGMjI3T6cR6z83O/2PF85JSex1oaJqQlDe/FCdlF+I89UEAQhTlMRG7ERG1MGlKpeqnqp6jX3aD21nlrPJutepXmV5lXBbDc9inoUdZ9z4lqBpgWaZhq1+srNBDcT1EuVc9/NPDfzRNT+oVR89/juN4ZE6O5t3du+vh9ZyzLEMsQU6HYQ7nDHIFY8CfkG1j3WPaZVBQ7kb52/9dXg8KaeOT1zhnR1DXud4XUGxwYbUxXqUKjD5QMBkZMyTso4IUODE9peba81om6ngDQBaTys3/x04sqJK5nsEqQ4fv74+SyNgxr6n/U/e69Nv8sR9SPqXyjf6jwAWHNqa+EHP2QAfOELsVsJgi1iwPo6HuMGAHAKzdMNAMKeNlAKROrSRo0merT7rQGAIeHHnQGgjzZARUA1qXYHAG4fDlkEAPPK7rgPAP3DV9oBwE/5MAgAlpwwzvcNDWNzRyiZiIpDWyomi6hYGCG7BgXG8TgWP08ZqGSHjADgmHrhOACIuFuFnhPtJx645xPdOHrgc1Gh6csUkwzcNXA4Fa7l9OCEBfM4MVwJgiAAANx8lCy2W0kfyuxM+btCw9ALRqymrcDPdVLy4EIlWeIKh5ga4UMHyG+VlCzI3QEHswahqbiSY+oqGcnUwFE0OFnpcR84m/NfHQ77Pyo5gvNXBHdzmmRlu5kibmosz1YQBCFOwnnDm5lENTdgMAZjcO0CLi9dXrq8dN9Z/mH5h+UDI8qY85vzY8Wp2iF/hvzp/Cqe9nTM0zGec0fvTPNNmm8e/WxxKpKySMrzj847O7RwaBE1SfvVOsE6QUuaeSVO4iQyAnCFK0IB5Ed+XAX0wfpgNLIb4VzXuW7EAvP8XClypbiZLMxeP6ofxcpyLeEOd82l2Dem+Kb41jH3V1/qf6l/6m+f7l9dfnX5Yutd6zzs+7BvwsNBJ6wLrQu1451yA0CSFus3A4D1l+Ys6QJn2+sUBMHG8iG34N+DPf3Fppb2fgCgd8oxGwA0h6apASCqVtlbAKCFXN8JADoGs0bTlpoAoOvIqn4fl5IKwj0AsG6ZqwNAybQhawHgXjt3ta34GP0+ACxZB9bKKpxfybmsTZKYhqx2TMFYdDzmlucvoOS8rkqm6KykLw1fc0sAgInqBtz1oeqHN0XQHT7unXRhbZNaVHx8ubvFHRYRbESD1vlCPCBEep8gCP9xuKB2Wa3ktzRQNbmoZJnuSu7i7n5tWBT9SBrOHz8r6XdWyQqMCJ7FBXpzpvbd66bkQUbq/rBfyW40KPUrouQg1hxBTSVGM+IripHBo41NQhg5O5CGsH7rlNRWKTm8k5IRjOy6f1NJR6bGI748ekEQhC+JVkpz0pyQP3lfh5QOKR0CauSI2B6x3RT68x49jZ4mrIO+C8dwDH2WLDha6mipoy1WnPmlxC8lfnjxNF2b39v8XnWOR86b+2/u9x43boh1knWStjxzmGm6abqeZoFD2gtpLzzq27GHb3Xf6isdUms4hVNahnVZUR7lcSQVtAXaAn0cEHk48rDdOT2ZuYe5h7XBvQtwgpPeNeyKnl5Pr+WMtwITMAFPgtrarbdbb7mxosNa77XeReNfXr6lwZYGhVy7n9XP6efwTaRygB/vEgIAmubH4KoERqmV9978SRVDf/NrIDVEFjW3Uk/Dc0k9FL5GxID1r8A0EwAsHXMHAIA1UcM6AOA5rc5tAAiPujwYAMKr9uXAt301AOh6IBfkHsoghbHcHc9nnJLzHgKA/96ugwDg7uiuVdTfPY3iuMolnrwAtzWfU1LJJKxl0o7fv3hHzO3Ok5vHHVAyJQ1EHR8qOZ//x0MACO+ifsnQ2qMuANwv4xQBAGG14fNx7qN3fSV7sgh9aSo0sy+xnV4c70dInxMEQQAAZxqESpTi/MEIplyMiEpBw5UnN7lwp6vlqoeSUYxoOkBHQTPWYvRj5NMPrNo6nTUGm3F+eJhcyXC2YyTnGzsarvpmU3IAl/CW3kqOMX7nkn3sBiXNrHE1gO3pw9TEKDp2Rg3gfMtdERElz14QBOGLwFq0ydMpWcMKf/g79GhUOJcll6Vky3zbPb09vb3LXH8QEj8k/v0pP/c7FHAoYCtOFAvsE9jnREHkKF25dOVU7d1nquN/U46Lq/UN/eN7tERLBHQKqh5SPWR/1yATQhGKjZE1tB+1H7FBvxPmEOZgfwq+j4Y+GpogELAGW4M1++CeaX5P8/sTx5DCWh4tjx6oF8VxHMfFeXvdh7kPC/1zTpnOLTq3aNvLoUbguMBxbj/MUo6eRqk4b3XmJld+YQDwHeOJdx7QWKMRj/7e7fLi/FjztJKnqbed/l4Zuiw/iiFL+JoQA1bcfCz0+GZnzY8G6wHAM3etOwBQ/HjqYgCQOfEsdwDY9rhXBQBoUezZUgBod9coPuv5PRf09FA3WabkXJqKOjUEgKTHgk6p3zUu6K2FASB02zc0ZM1hSl9SRm61oyKwsBMbfNi2/fno4fbj/oSpzyvpy/PMZTt0m6Os03N1AIB5t7q3AYAuP9mdAoAw6Hf+3n0091GylJeSvalY6UwpaUdFai+v23JF+p4gCP9tXFcqWYSGfR9akHLSsHOcC/CWTTnPMMJpCh0D5ZkqOJPjfSumZD9mbaq9T/i9TM2ewxpZ5Vk7a+YAJVsw4vcBk9cjVlMRGccZgzWz+jL3cDBTE7FUiXFVlYxkhPAozqtWRnINZApif0ZamejgsefuiDpTyD9ZCrsgCIKgYOpcSo7L1doBgGP/xp4AkLtQ7mJafi3/z0vsjmfMnjF7kvNBD4tdLXbV/2mKdnomPVO+fBOS/Lz7592ViywqeXXS1UnHHQ9PmFRlUhXv6q0Dreus67TVzRj5O4+RvZ1Y6zbQ9Gr+q/nOKYBLLpdcUq2xGxThH+Fvf95x4Ovkr5M7hgOJayauGXhRf5b8VPJTz5f4DzdHmCMshzZl0B/oD7Sl83tXK1et3IGL535J3SN1jyW7coRiHMYBc5Kq709PPakTHTazmZGiDwGAQnTv/9QwlRMAdDnkEgoAIdd1zptvSSGMQtT/u1mcKsAOdrBrzV15k3CeXM7NslY8V4asM9SnXueTLif8qy0lcgviAg6BSmZnylr9swDgmrj6PQAodjJDWgBoHKovBIAi36uQU925+wgAGGUKSg4A7WG3AwDSzvaYCwDflRu1GgDMDXyeKPPMonAAOHelmz0AFFrxkoYdRz82pBYAlNyXIyEAaAvm1QCAyFlJ7QAgdHjHQQDgmWZxAy70GYllHQUAD87kfQgACbfOHQYA2tXUgQAQlqFLLgDw6DFvhvq8mQYkfR4A3LuUMzMApGjkdwcAPPrmzgoAkZFn6YG35+fNr6lodIh2AyNtFY5ENJg1oazPVMd9NAgOp8J0jwqPiQqZvVGkPlz6pCAIXynG+Mbd/Vw3KlmIjo/G82hQuq5kfG6yMZSOjKH8XDjnD2sLJVvREDSDkbxVjVQ91kJsy+Mfc3fXfWuUbPGLkrNoAKtgGMBOKtma4/sDRnRZef6RVAjsGynZkwv9ITR4mcorOZHfF8nrGsvrtud1/cr57FdGmtFBj6WHbOc53JGuIwiC8I84n3ZJ2iVpl9ilTHAywckEJ5PUvNH7Ru8bvSuHhZ0KOxWWvV41/ZH+yC19sl56Fj1L52+SBGIiJtZr7NDWcYfjjmQBKzbnGJ5j+KLdYZmt5a3l9fAGDvpwfbg+uegLbZm2TG/yzfdJKyStENAleHiNIjWK7O3t9eD+ofuHEtfZG5A6W+psT44N37Zkz5I9rXcGX4h6EPXgliccBu8dvHdG7XjTryy/stz7cJtjCWsnrB1ULEnfHEVyFHnofucP1yKuRcIeLW9kGmYapqdacCRiQ8QG+zzn2vZx7OPYrIo1i28n304dpublZiTzGImcmg70TtzcatbgmG7Gb3rOqjDD/DjZzFzacG34qOzJ57rkdMl5s4B5wqtyr8pFLcbjiRMnTpw40b7PXw5+hEd4ZJ2MaqiGavbDcRqncdp8CxGIQIQ3HUtdOH82aqjkXka2reWu8vuOcd57xW8Noz5FB5VRK1gQ4hbSMb8o9jSgdOVzaM/aH4k2AIB3DUtaACg531l5Iha+XAwAG2avrgAAjw+/5MLa9IeSlpEAkDdtxkIAUOTm968BIGxWaCAAbNq+ugAAPFoSxIFMo0FHV8VyowplA4BymfKXAoBMB7XeAHCn65XLALBhws65AGBJAO7mBJWC1xBjAaD6yTLPAMC7eKYZAHAvz/UVALC+9vYmABCVUGcKCQwPdwIAqPK89HEASNs8awEAuBg//CYAbJ/14qj63DZut/6qJdvN7dZ1RgRYmPJnYkRVHu4iVYipjbep+Oxn5NlrGtJM9/g9r6QvCoLwH+E5x1GmdKcdqWRhRqB63ebnmKIHFlO/T8PXRu7WGsGUdTNrClpZu7Agv7cQi+sau8geZETuKR+Ov4wMtjKVPD8NTd8YFiQu2I8wEuwYDUgmGsC0INt5tEpeJZMwxfG1Rcl1dGQ8Z+1FjbVG7Lh7bmVeT7I/bG/TSXrAjzIyQH8tXUcQBOEf0T9LgiwJsiRwf1j2admnP5gKHw0/HX46dEz6VQ5nHM5oJ4PyO3k4eWTNFXgs8G7g3Uy50hY2pzSn1KPCuqVvnb71gyMH/3A+4nwk/NdUO/WV+kqMyJANznBGOIBXeAVnQNuv7UcbACmRUg8BdLNu1naEDNM8NU/8dPsERmKkpu1arEOHVuX6iIi9EXvNt8u3eu7/3N8j+fev9O56dy3tuuEpKqao+Gzl4F2ap+ap3zlX7cmRJ0fidYkaNXjz4M0NhxuXk5Pz0nxG+KanY6TLNRquON9ZbVNOkI0OnXll0QVdnifKkzvHphybfBZrXWu8rvE6d7n1lRJXSlwpWev726xZrFmsWbTbf7mTVlhhRVukRmqkdpmIV3iFVxUtCEAAAhK5QYMWg5bPmsbh1ZR8QP3pUGqYYILpxAYVyVXloTKEHaYjqTdLyOCWpCIKcQGJwPqiOHEA+X4rACSa/jwAAGq+cjsEAE4nPXwA4MCrzbUB4FjLS1x4m7lg12jxD38JAFWnl0kJAOkm5roPAIcrb30JAIc3XGDKhJkKiUYPuaU0ABRfnScrABRYkH8jAIRYApMCwJJ9q7ICgP/CYFVDZKE5ARfyHgDgNsAxPwDUv1FPKTZX9K0AsNzTrwUAPJ4W+AAA0MZMS76KkHK+bT8GAH4OrFcCAMwt7BICwJ+H/MYDwIPH8ZgrXqg522nsMsXaJvoLJd1OK1mKEVbJWYvlIHPm/2CKioWKk5mRV+AuiBYqOJgmfVEQhP8GOiOuMrH2VB4WRV/LTTie0WOrcZONfBw/SzLCtSUNXcd4/AFGTpm52cVBjqsnKCvTIFWUqQxJuZnHOo7bJkYeH2Fk1UmO75XMShZmSmNuGsjW04N8n8XjoxgptYgRYCnp+KhKh8VPLAK/lanuFxgZZsf5ZAmX4sk531RbxPvEFJPImTJPCIIgfAysqUOzhWYLLZWltp5KTxUZEBWSfkH6BUUzb01ybdm1Zcl+LZMv5M+QP7HGq2W6NunaPJy8e0iCJAmSBHs8nmbZZNlkKmc3VL+h30BObxc4wQlngTdmmniIh1eAXk4vh9EAkiCJ9jsAD3ggwqGjHqaHYaYpkbmuua7liXuZ+6/vv05UvM6GsKthVx16FvXYG39v/FwFlj4KbhXcyuVQ+xVuid0Shzo8CzjV8FTDTJEYBQD6/xuuuO3UPOo3GegI78x5btb3MRuunPgN48vSkFUc/dDvWf3BEXZp7NI83tKkoCXAEpDrgbULPOCBxZE0HJn+Ov+YYIJJH4e7uIu7low0WMWjPBfLA+inhD0DGcyskezWAZ7whKfTXKRHeqS37MUxHMMxhMRiCBOEL4p0yS+KUez2zw4AMDoiazwA6GKfdLz6u18KJbvkUfJlYwDQ+NT0N8NikVkcSKmAeHJg/J6pdmcfxHxcQVri57oBwPOOWTsDgF/W5YcAoFubuqcBoHlRa1sAmHHAOM6OxQE7cOE/nArBFnrca1DBiRpge71mtrMNa06NpKFqz24lK2ZR7axIhaknU0/q8XP3qFClLa1kd/69GBUXPw7EszoqGXRV+pggCEJM1OQ424uOkAZMHbi8zfZziZiKN5UGoRqcX+5zJmnG3fy21oz5PJXoOFjOcrVBrZRszsivDatjPq4kU/tmcdxPT4/1Bu5u2JIOmYfxbI+z48J8ZDMlO1EB2EqDW4vvlLx73vY4Rxq+ptFw9Yi7+PZiLUjclz4jCILwMRhiAoDkZ+JZAOBBLkfWSGz6mOv5XRy/myWrnqx6QKvgZr8W/7X44oXm0ZaVlpXWfaP99Bn6DEzseFarqlXFAUD/Tv8OIyxlUAzFEHSzFHIiJ3Ls06BDx9m1G137uvYNz3kme7PNzTZ3y58gpGjjoo0v/D7vds5vcn5zM3+q7/Ym3Jswd9my7U2pTKms7Q8nu1D4QuG0k6LrTXmYETKX+k861vTtzHlmZq2Yr9f8k5Kt6NAZzdIl+3wAoMLLpgsBYLPHTqbCt6bjZGesNbBsdyFMyl2B1/P8+cvwH9SXAozzcf7bEKrkYWaqXA2GPexhr6VWtbVmXoYOHfq905wHuQmWRGAJcQOT3II48BCSmTwAwGNwSqY0+DFVoysNMC8b/+/n/38gLcwBaC6LEyZgsVpfpvidzRzzcflv8ThGKiWbDgDDXU+nBYBf2ocwN1ovAgAzar8ZgJnC5xuo5G8sRnjyCNtLz0N0w5WJRRPbMJJqBAfY8zR0dR4NAFVvWznQalVtB14Hph6Wq8N2J1MyExWitmmUHMcUEjFcCYIgvB0rF9wZuftfU47TjkNsP+fP/7dJreQqGnZS0iUyq5OSZSfEfJ593F1pLWteJeU8N52Gr4qLYz5uzwK2iwvoW3ToVOLxM2m4SlHI9rio+UouZArgNUZ0/cDdB2dx3ZMqu+1xEXR83OG8qhm1r+ylrwiCIPxzHGa90UgyAECdMw0c1e+N/aj/MHKo8zdKBjfrnLVz1j+vAY/WPlobP4d20jLdMt20x/7HqJpRNc29AP91/us8nz2LCBsQNsCh56/DTKlNqbGzQpRplWmVtqFVXt8CvgU6+Kw5pAxXXowMnvvgdfHXxR3LZTuY2pLa8uTykb2NAhoFbNl7xaF2rtq5dlePrjflYgbKXEYEp69IvYcp9LOWxaLh0fDT6qaSI0speZF6VsdOALDZ4zUjtvDt37uzGh0wYCbKk91KTmamT20HJRtQT5zB1MBzwUqGA5GIRCRmKsOVVkN6qxCXkRTCOIB2xtQOAE62+pOGqG7K8t48WC24Z0Y/ohANRnO9lEzESKS2HIiW0sONA7bH5eNCfy49ASkGA4BlY8cXADA2e8ayAKBbvPn9UKmGnTQWUW/PiLGhNFyd5q5RTeghv1YtWkPZrtb0dP/OiK9L9Jw3YejtpQAAWHeXR63WlQe9r4synZ3vyhSXEvREbKIHYxxTUR7QMIak0psEQRDeC6OmFPcx6swUusjRSg6iISuckVlP6RBpw0hfjRFX1Tlyz6SBqzk3w9g2mQoIDWDtWWvKxAV8bS6wp9PT3pLF4DfOsG3mXkZUNaPDYg4NXxWoWkxjMffmVAgesx3HOe805Xnnsl3fn+ECnrs1NeX3P6Dhy8RdEv9/fydBEAThn+NAc1XkY9M+AHAL82AE8FzWkOrCCKygtcYxQxsObVg/GRAREBFgfzVvi8IlCpe4dLXCUesw6zDtQtBR/xH+Izwfdi6/IP+C/EOWLOw3+Pzg84236YNWT1s9rfivxrfkYUr6PM5XKbMmKZqk6IvJqwq6D3YfHDrTslbrqnW1PtbbbLTbaFf6jHFcdhp6/KjnZKCe1oWRyTONGonO0S61nhItOJ+OoOHrKg1VPkkAwO/cJepxCQ39xfz37mw4UwPnMmbs2G0lT/G+ht+NFjn1pvavbSSXIPw7EANWHMCS2OICAH9g810A+EMLSv2//9ffGLCK0EDkx4V5IioIbfcpuXhOzGcoyOLmfobhirtTqBRAu0rzqagMMTzhVwBglY8dF/jt+P2/MQLrFBUGH+7qdHWi7fnM9HC3SqPkSBqcjFpcPjzPxQCll/A6Z6ndBJvuzFETAJb3zf4TALyqa2Vo7ECefw3bGyGecUEQhL8Hi62/5q59jpwXujNSysTaUQOZUheWVckndGC02s6Rm6l5P9IgNpvF3Zv9ruRWLvyf0QPclp5rOClRmwaymYzoasEI4g3RIml30SHThAanWdzMoxJrisyg4akVv/chv3c/FQMfdx43SskfmFIxizHGLfl9VkaMSYkFQRCEj6rvMH/CfqKqvXQk3Va1223TjoxQCooH2JY8CcoWlM11YeG8aI3WWOd3Y2+9vfVyDY1XyNrL2kvb3r4kXuAFHBeE5ffL7zejMdRmHW/yVvLTIOTHkiqpGKHUoXol10quh1+GRujd9e5o8ZOf1zmvcyEB1qH7Cu4rmOuHnMz4mMddazNwt77OdMTPoB4W3XBl4v9bUj8bRYPXFW5K0thLyfM/A0CTXB/rzr6gY2cSM2z0JpLqJ3zNiAErbsBd8rRgADDxqVj+NP79DQfgudwNKiENN+3SKrm4ecxfm9/K45jjnIKhrL69ASB9xnlMVdSYkqGPAwCH1OanAFB5lu909fehHJBPUVFoTEXkWmvb85lZU6UNB/YRrGFygYpMY6b4XboHAMuYIDJihttoANjUrnVyAKhW59cMALB64BNOIB153QdK8EQ+0mUEQRD+CSbu4nc9jZJHxinZnL93X8oPcjOMgUzZC2NK3lN6klszgkrjOF2dnt1Z3CykaXIlt9GR8YxF5NsyRRzcBak2I6pmsNZVMxbJ3RQt1W8X550mTEFngDAq01EznTUkWzGF8AGL0u/jJig+PI8fI7bKUc5creR9GtweL7e9fkEQBOGfEGq4uyMtawFge/yTEwFAc1eGK52j9P/U6jUiaP0wFVNRJUk5K6zQqnRgaZEFYTGfKR8dL3NZMzEVXeUdi80YPmP4qHvzjoYeCD3gmLP2j+aF5oVWB0vaIe5D3Bseza0c61kncN7ISEd+F5rephslTqIFDJhoOGrFyORRjFi+7KVkQ7WpFS66/u9Rum4qDAAnplW/CQCVEiWpAQBPaupb3nUvbQ1U+lzpXcJ/agUrtyBOEQIALYw4qKoF6fGewyEqERfS7Zkyt2R+zF+Th55lP4bgpmSqXkcakuY9A4CW/x83tRcAzHWjtgNAt1Hf3QAAO9PgJerfZ+kxb8Ltya+1iXZCplq0ZIrG74wgu0zDV5NVABB5WhmuOjNjvH77dNMAoNvy0XMBIFGegfcBQE+quQFAaJ0nLM5+p550DUEQhI8KF9Y655Wp+ZXszUil8FNKduurZF8mGThEM+g8pSGqDSOg1qzmvMPNPWawlkaZaPPVMxq+2nHBv4Ie6mR0gEznAr788Zibv4cKQnN6nu8y8qoSI3WnMmUj6Qbb4w5SgWjGVMWb9JCX5W6K9fh5fZh0EUEQhI8IZwdsUJkeelc9BwDY/aUASL7X1GPoSE9KR70vN9NYWD7mE+Riyt+cTUqmbqJkZzpc5k58cv7J+fgVAftV9qss/Sw/P379+HX8sLST70fdj0r0Yhpr+WZkiZRfuLvuzEc8wY/RTsj5pzkdNSN5/DWm8jXpraSt4QoAHTWNfQEgld2oXwDANYcbHf26l3QWQYgdMWDFCfRTAJBms8cdAJhaqyoHTj9uE56UKQ3tKBcx5UGfbvs9+emBnk/DUSoWrfXlc55LBUF3B4DqDKR9oGUMBoBNDzOaAKBnx6RUAM6UUrIRF/qXM9qeT2OqSRsqPCNpcLtExaVRTgA4Xf38RABwyGtfDQAmLKh0FwD6Jl5UGgAqr2qpUkJaXOoMAOsP/zYMAMIGhhnFdA9IHxEEQfioMJIWTIGIosFoMg0/PTjuR9Dz3J0RtkMYyZst2gL7MWtUtWLNwjVc6Kemq2T2FCXLutse559Bybb0YK/opmRKlQqCGTRwNaTHO+FS2+N30UDW1HB40BNdmSmR01ncNula2+P2cbfFJowMu8nakU6MTMZz6SKCIAifBCNJ0B8AIq+5cn6oNoD6CjNNUqRRsgMjdefTkKTXtv26PExln99LybQ0JHXkcbNZC8q6N6JWRC276oB/U/+mnr3ckz/u8rhL/Nm5ylh7WntqWbxZI6sb561pnMesaW3Pl5rzYpdFSo5mpNY1BhA0oknu7I5o102HfBNGeI3rDAD6nKD2APA6MIR6D6ZJFxEEIY7izqLoW4cAwMhvbn0DALpu2avkcz8l699SksP2X2T+DUpezK5k0AIlfdzeflzuC0qeKQkAeg71nxe9LjgCQETSTJEAoEdFP85cT8m23yv52kXJ40WVzLodAE4mMa4zyU0AsM848A8A6LDsqTMAvLijvlHXTxXgccyBr0TP9wF+g/d16SuCIAgfk5qMyD3GhXfmutH0CzosOuRT8tVvtjPBLaZwlMoc8/cnpqKx+rLtcfeouHw/MebjEjKSdznnA72R7fErGWmVJH7Mx5emoermFNvj1tFDntwr5uO+ZarJbRbDHcYkd6SVviIIgvAxGZITABLnneisRulJPytp6AUB55Rs2PXtekzePUqeT6tk8Colmy56+3E5nxcvWLzguZynPVomb5l83eWwHIlvJL4RYGnbm5/LEfNx2YoreXSUbXtPVFBSFX3/63GmEko2WabkywTUv4oAwPru33OX3SsVlPyuzud6EkYRdyWdtim5sIKSvx3i/9NKsXchLiERWHEALdSUAQC8uqU5pP6yiSGnzZhMuHhXzEfmoYXej5b8lCmV7MQBcG7PmI/LwYW6n0rNyJ12MACsqXJhKwD82O7ESABweHxdeTBavWkpi420ZIrFSO5WeJm7HzZeCACadmkAAOR7UuAsAGijZz0FgMaZ+qwEgMH9E70GAK9U55ia6LMXACpXuZSO3ZKRXVKjTRAE4dNgpWsiA1P7fmbKhh33X9KvKfkHDUatOR7vKatkmgFKzqKh69uutt//lB7qNqyxuI7Hp2Rtxpn3uFCvYnvcMxa/7cAaVAPoyLjHFJIf6aGewnYkzmd7/C7WKGlGT/cdOkAqswbJVKbGJ71me9xBtnMZ5zetvvQRQRCEj4fDmz32rMcAoG7q+pls9YpNjORtznlhQa+YvynXMuox/D01UwS78Axz0sR8XDbWcPQr/6rAqwJOh3IHpY6fOv6Twwc92yxrs2zt8aWTezfs3XDhCpy3PS4zd8mdQ0dPNkb+TmUNY5+BSl6INh/BiMhiZsq43JyXLlFvmgcAlUecZKqjkdkC2aRKEN6CGLDiAOaEptkAcO72+m84oClL/MDVWwBA0/SmtkfkYeTWXO6Okaajkh05gM6lR1zva3tcdqZGGKG5GXMBgGVr9zsAUNtr7SEA2J0SNHzpaheqmajOCWa1kiOZ0nGVIbyNdqh23uIA346GrhlWAGjuUzotAIw6Y94KAJ4Xz7EdPmzHGRcA2LD+TUNpkEOE9A5BEIRPAms/2dHQ04O1C3uxtohdoJJRdHgs+EXJBkw53EyZjgvz2fx/iRy2p3mkiuKiFf++jgav1IwAm8V9mL57YXvcYxZTHziS88/PSj7kcTVoiJpCRSfhctvjd7HofHPuvnuP31OF8+c07mabmCkjFh4Xzl0WNQ/pIoIgCB8Ph6kcXU+YHwOA5/n4p9VfVnITkIaD+ftt6j8Jbb8hFyNp59IQlN5Hyc6MJJ7N3QL1orbHZWWNxrmHlcxUINGZRGeCmm3Y7/HC48XrMQEvE1dKXCmwviXbvo37NubaZhyXkRODXwKen5uK9GXJlA55lDxXKOYrbsyU9vGczx7weJ9mADDv0GnOhyZjv3kxXAnCeyAGrDhAVIglAgAmdFtbSQ3Yz5cCgDaAw/Cb3Tjy0YO8gPt0pGWtEV8W0Z1DA5j1WbQBn4rBAh6fuYySXX4FgMI/T28JAJE9dCPSKwUAzCpgdlS/tt6j5Bhuk36Z52ug2qmFsTbKpPwAYBrXsioA9NkQZQaA8ZUdwwDA894ZGtQa0aN/0qYmia7DBQDqxs8yAgA8fnLdzH/lll4iCILwUeGCOpQzjYmxtr/SgNOHqXsOobaH3afBqhlTHDYx5TsDI6X8WLy2ZAvb4x4xsqsFz7eWBrQ0NCTN8VKyzLGYm7v5O573MdtBw1tNft9URnQliub42M6ivk2oIN1m5FlVPyVnVFMyeS3ORJKyLgiC8AmwjlHSfqV5CQCcyLeL+kNb7h77vIzSK6LrP3mYkr6ADvAMJZXsyAinmazqax1se8YcjPhdwHkraykluz6s4lHF42Dnufv1vfpe5EA597LuZV9fsR7c93zf85zTMrMG4nzOF7m5yUmPq0qO53wWtc/2fBr3j2/MeWwii8nffUq9id9zfB0ANC4qfUIQ/g5iwIobGBb3HwFA+5UD96Y3Azc94HNZ5DZVYyU7sYj7nB4xf21Ofu9cFoXPSIWhGw1TM/YBwIkdug/PuBYAzJNNTwCgUaeWLOY+6rSSl/MBwJ3hzXOoCSYNJ5CFxQDA1DXJDgCYX+xQNwAYvCNrHQBwPniRCkFjKgino7VTK6Fkw9wA8JND33wAYFfYTM+EvlC6iCAIwked/rkwv85U9BlUDCK/VbIPU9l70lBlP9L2+IfzlGxBz/MWRjql5bwwmw6Wb1PbHvekgZKt0iu5ngYz1mDHLLaj9JyY2715nJLNRyt5n5+ryaLvf/CLEqe3PW4HDV7NWNT9DmtwVaWBayqvOyl3PbQelT4iCILw8XjNfJIIB8spANhc92g+ANASBXz3v5/7f8NVrtvUY+i4SEtDVyc13+iz6BB/E0JLshfkcXTAZ6E+1P20ruu6463pazAMwzAba82dzJ2suSy5R7wY8aJehqw0qM1tpmRu6i29aRibREe/niWaHsNNshozUnjCXiXvcn5qRD3t5JFo1xkBAIcOlGsJAImsiZhCr8v8IwhvQWoMxS0iAaCTsXl3VC4Ws/VjqKpRI6sTUwT9dvKTVW2/JhsNXHO4+0YmGsC60YA1Iw0ATEuu00OgqV2YmkVdA4COXqU0ALCrX4394+phAFjboLsbAFRbWJupID+yiO+K6wBwuH1oQwAomKa7St0ocKe1+n8TpqSciYz5sutT8Rg/FAC0pU7jASB0fjg95tpg6RqCIAgfFdZMtDIVfQJrEt5iRNIQLrh7X6B+wAX57zQMRdEx8oC1qFpwd8AZK5Qsx5oes1gTpClrhuxncfbHbkq2YurhNDoqKlMBmMVU9mYs1r4rwLb5W1j0tiV9+rM4L9Y0dqniPNiGBq1nVDx2GjVW6GGfQwNaFf4eQgVmfAfpIoIgCB+RX978tA0ArLDeAgDtLqbafjAHaxD6eSqZ/jgu4iJ8umSbUHVC1YlDd80x1zPXsw5InnqEywiXuvEe3r7T606vJCOyMBJrDmsmZuW81Z36zrTXAysNrDS9HpC4T+I+LwpZ9jyJeBIRr0cql1sLby1MemBaL9RGbSAzUwb7cDOPSYyVsk6N+cIa8vvHnVbyIVPafZgif8oc83H1vAEgfb/xUwHAfatHFAD4Qw951620LaruzHk6I2s7XuUu9GHFfH19fX19pesJAGBPA2k29pMHXEc9++XfdiUSgRU3YPFcFVw7Jio3FYL5xZVMx1S6TlQAZiXm5wfZfk1O1gKZzxTCLCwC39VPyakscmg5BgAB9EfUNNvvBoDlk5MWBYBfK6UIBAB93rVAAOhydlowAFRb2IkTQokmAGCe0iMLAAR3uLcDAAqm6cEIsXtUiBqyKPuJ6IarUkr8zKKNk7jL04PLALBVH/ErAITODDXKx++WLiIIgvAp0GYpGcH5YTTnm1+5ANZHKdmPikBPGqgcfrb9nnuMwGrG+Wsza1RlZAq4H33vJQrbHvfQqHHVTcl1TPlIy9SMOYwU+25JzO3f1EXJ5sY8wpT5Wqw9MoXbnCdysT1uOxWUpnQU3eXnXbfwA4+kbwiCIHxS/dMFAFrFM/6cm5Gy86k/ZGTqX+cLKVOmTOlvN+MYlmIpBlW/aF1oXWjat1Jv/rL5yw2NOjdIuj3p9oDJ84pqVs2qj8q2Qx3XpSHnATr8LW20p9pTvSnw4sSLE+5PXLc+8n/kn2BPnpnWDNYMWqtMdPz3Ym3E8ZyvoqIbrjhfNOK8OYGGrvt08DTkdRyPzXDFeXHSCQDQfSOCASCswGs6kvCBjnsv6nuj6BiaxcixSnOUoSt+lOwiKABOdOz1YgT7MmZA/UyDcTKjhum9uH4lEoEVJ7AmBoCqLXO3Vb/7ckGdlrWmOrNGyCyjY2WLNg+0UbIjd83IR0WjExWCN+Nusf89qtcdk1qoL2yiapYMaHgWAB6NetQRAH4ed6QHAGzM3ZIGsZNUQNqlBYCoNpah6vd9/H8wI6ka/KHk2SoxX28yhvz2p4c+NC+PGwYAU0qlTKN+L21sI1tE+oggCMJHpZMSOmsj4hTnIy7cR7KWlYWOksEsnt7PS0mN4/vv3DUp4pWSD/j/ptzsY9ZKJStw16iejAg+ys+FBSr5iCkizbkQn3FFyapBSs5m5FYT1tbaNcP2cjbS8NSUEVkzWRS3Ng1RO/txPoy2O+82zmtNuJCbS8cR1sr8IwiC8Ckw5h19PQD84ZaD4/I8RoSk5XzjuzD+tPjTXk6c06j2zdo3Jy8DtJ+1nxOmsMuij9RH6nULX/Lo7tH99alCWRunaJxiSwFt/pMuT7rE7322VO48ufPcGHGhutNEp4nhMz3SpJ+WftrDHEGpf1jww4KRLfWBtb1re+85UdDVLaFbwrDOTpVcF7sujjjT49dQhMLRMo6OGty1bbfGlPOGLAY/kUXl7zAiuSFrcp1pEvN1J2EKfT/qPRHlAODMveaBAPDw1FjWCjZ98647aEphSmFKAeAWbuGWNau1sLWwtXD8pjiO4zhe9iHMMMNchSVjjoxVcuULZcjaxvt8v7WWVEuqJQ0bo5XQSmglcF0vp5fTy+mVkA7pkA7Vqi+rvqz6MuhrXNe4rnGVnvsvf+/YP13pGCxNA21R1rg+xfXVWq7f1vL3W+OUfG28F42/9JWIASsucN80AwAqHG3FYrlhFZTsylC/OcbuTMejDaQc4FozF7sePconqEisoiEM6aKdkdumN1O7Fv48ogAAnG/p1AcAGlW05gGAU+1rcHfDbfze1XMAoPKZ9CUBwFKhrQ4A+pS0awBg773dWQGgYsuEqQEg/EppKkjW1gDwXQtzCwDYPL15UQAwjc64CQD2d9y/FACqfBOPIbC5+YJ4sZZKURronlDx0RJJpxEEQfg7WE8omYOGoyQ04BRnTcUUXOia6Wg4Q7mNEVBVfZTsTc+yE2tobWYNEPNFnicZ56HtSpaggawMI7cGMvJrPc9rpsfaSkPVMkYaZ6djJj0dODNYq6o/58W7nB/sWD0yNA/bQ0OWUUu+Bx1Cr3m+WzSM2Q3nwo4GscM03KWlwa204Yh5In1HEAThnxA12eVPlz+d46VxyXAtw7Vi0xKV0KZp0+40LOcd0iSkiXlRrt88U3mmujxxq31ha2Hr6ilhubQALUAPratH1Imog6FatZD4IfEdp+Z57gIXuOUEtL5aX/2Gdill6ZSl/V8D3hm8Mzxdn/2W9Yn1iWncqoDIqMgoO9fT4y7nv5w/yb6dlyrbV7bf08gui392/+zxfmx4OEG/BP1Ojn6ZwFzNXK3thMAxALBAK8F5wczYsMgJAHBha35/AMhWtv9WALi3G90BYPy5+VYAGO3rcZvzBlOyrL0AoERf0zkA2H6gWX0AMM3OMh8ADrU6OAZN0bRmoHOWhO0Stks4w7VVgWQFkhWo9F3ZTMczHc90PGF5S01LTUtNc3j0O3lxzcU1F9foXi5WF6uLNd6NFI9SPErxKP5Mc3VzdXN1AKdwCqfceRwDxVCiAXTo0O+nQWmURul9U4P3Be8L3nfUdMX3iu8V3wDXTHUz1c1UN5n9y6CXQS+DTFm9q3tX965eOj0iEIEImQf/pYYrllhwZaRVwrG2/3fcr6RhNy3A2tVtmWp4mBHu+/xoZ2Bpo3MMvAnmrtRoLgas/xDm2qZNAGDtHMkQ1K603M9mEXVrcDTDFRfardkhR3Jb8EtcyPsw9PTuiGin4u4YzWkIG8OUv2sbAGCc36omAHAqsjL/Hsligx7HASBZrx5LAaCNS3ZvADBv8kgCAOs3XigDAI26RIxVhquuLKao9wCAxPM8JgGAb+Vs1QHAHp41AGDrhIu/AkD9kNAtAPDyRDdGiCWl4Sw5Q4dbsZxjKK9PayO9RhAE4W8tZH7kOMtUB2MzjgmssXiZCxJ/biRuGkDJSKxTTMnIxlSFHpxPGtPBcMGBigpTRMx0TFxhRFU2eqq7einZcJyS5xkxHFmd7eRC6I6RSk7HTHru5jSdNULusJj8DRrYTHTQOLEG4zVeR0amNk7txeNYDPg6q4OYaQDLyFRDExd0XkY54XDpO4IgCP8Eq9mtsltl97ZZpuXalWtX8WPpm3m+8Hxx6lTYXc8Hng/uBV6r7TDWYezzQ0ma6Jn1zAkm9sqG+IiPMP08HuOx9p3mFD4tfJrnyPRHXa66XLVkBBCOcDwE9LX6Wi0C0C26RZtjvoOiKIp7CaB5aV76N2UQUjekrrt/mQPJBiQbEOwJez1Kj3qZO3Sq4zPHZxcyBS4MORFyYlVYLUaWRE5hgxMCQNtzWfsCQJqAlOcA4IXpeR8AaHHoWmIA2OJbko6O0lTsrf4AkKiO+y8A0OlZtsEAYL/T6y4AbN9zKSUA1GsU4onKqPyyQWuvdC/Tvcw6Ip5j3gV5F+QdWX+EQyKHRA6JqqRSkVTa2b/M5MP0Yfow7Db3N/c397frj/u4j/tJPeEIRzjGdv8dFsAd7nBPBz2LnkXPkna7na+dr51vuXmJXRK7JHa5HRpVLKpYVLHUk+52vNvxbsc0gfCHP/y9XGQe/Fev+7i+sWNJiPRp+I87sZiHmGzKUqmoxUj24mWU3EuD1gpGrm/qZbNMEwPWf6Rb1cQQAFg+YuI+AGhbe+YdAAiOQjAAeLx5SmZ2kNY04IzkLhkXmGLYaCsNWVtsz2BiKkQzWlzHUmG5xtojjTMAwKzIiyxmO4OWV90fALpmSXQIALpMGbcZAJIN9ToGANb+U9MAwOKtI3oBwIsrEQd4XDsAaBsU/xUA9Fo4eiYApDgVXwMAvfXMDQCwwn5YVwB43j+cHnvLVABwMVeKAoDw3O1NAGA52YOhjPeZyqHVlF4jCILwd7AwAqoyi5SPY8qc5xnbZcEwOkBu0rFiZmq4zpz0JpeU/JURXclmK7mKhqeROakI0HFipGI0YE2qATRoJaUhaR0XVMOoCETQAKbTEBWvE7+XC6lK/HvC00oOYcTU9q6c9zbyuriJyO+8zmo0TCWi4e03/r6VK7WerP1gl4TtYVF7PJC+IwiC8E+IahlgCbAENf01ydq+a/uGpG12MZFbIrf8v/bq9W3ub3PPnb4lXtTyqOUuHU3LtTvaHf2ofsryveV7u3AtfPeq3au6tW0xokDeAnmDO2fdWjZJ2SSnfAH9R/1H1ASgQ4cGoAIq4MjL8ViHdSh+cdL+Bvsb5OwQ3/VUu1PtMmbJePqxz2Of+EdONq0+qPqg7W22zcuaMmvKxPGyH86dM3fOXJ26Rx3AAZw+/uIMAOxe2GghAOQ+810+AHD8+epiALh2sV0gAGzJepmOFDNTIK1VAaCVnZcjAPT5dXR9AEh5KIEJAPSOc0oAwCqXIaMA4Nnw0PnYjM1omTDB2eZnm591WVs2+FzwueBzg78rvq74uuLrDqwJbRTaKLSRnftfTIHuVneru540aYGkBZIWSFTN84HnA88Hf7TEIzzCo2zzo32amT2PBsEZznA+XFWfrc/WZ2/zfj7n+Zznc85aLg2/NPzS8ICGpRKUSlAqwW/OOdPmTJsz7cPni4csHrJ4yLBWqIRKqCTz4L/U0uCnpDvXN1NYg61M9M/5KBnAiPmz/PxmGr4OMOPrMu0KL7iushqlFzZ+rivS5KF+SdwZObXytpJLmNowqyQAWGvzIS3T+P82RsQVI7Su0OPdkJ7lC6OjPV562pvxPGO5UL9Fg1FDDrxnrsTcvvjc/nUSI7p+uq/kH2xPj/VKhrS0Pc4rPxUj5no3ZC2SWdyNoytTVl6WjXY/pgLAb3cn+wDAqEX5UgNAQK/ylWjAmi19RhAE4WNQgwvakVyQujNlIhH/vp3zTVPmHtybYnu8HQ1Vv+SwNUjpnA8GMEVvJOcZC8dxs0pdR2fOL4NZI8vM+WIII56GMZUvKtoKKwVrI87geSvQs3iD5/Ohh3h/EtvjknIem7ZKyao0vN1mhFgzWu5K0CHkyPP26isGLEEQhH+OA7UX/cVvgQCQ5lqSYgBw9Xwz6jO6Z7RD1C7paM7I3VEeZQeVHXQii2v6qvmq5jswwa67flO/qV15OgjN0RwbtlY0LzIvsjqsujp0yNAhDfJmuP6o5aOW8Rv1u6H31ntr6W9T0a6fYvKjyY/GB2esYy1jLWN6UG+B7qw72xdtPsz3hO+J1i8qMEJ5ClPWHzN1rlEGJY+9jPkKPRmBPIa1iH2Yw+7H6+hMx0lwWgDIwPjghJ2VA+hYyL5CAGC93bK9+s/O3bHdy0n5JuWblA9QBqukea1zrHOsc9ZPxWVcxuV86WCCCaYLvdWnVzEzZg1TyC4ygCG0urZYW6wtRhsMwRAMcVqqZ9Az6BlmloMnPOF5r39qPbWeWu+14m79u/Xv1seDDps6bOqwSfryv9TukEnJxVWVrMR+cIvrsw20I6z0V/IkI+CDWYpI3xBXrkQisOIG3ObbFB8AHNg9tIqorH5qlV3JEfRkX/1NycbcdvxCzpi/tgkHLpYCwW3WMmlcV8kzlQFAoxlTNxIlwF3LMZmKyE88z9RCSvbkdubRDVeeVATGM2e2IYvEzWaqSTfuTvgXwxWLIo6tBwA5czVsBQB65vPcfh180SAGLEEQhI+CtpR2GRp0ZtKx0IeOiu9pkJpFg1ZTRm7d53wQRcfEaNZINNGA1I+GsP6BPJERkcV5y8JI4vGsLaLxvIM5D/aJ4IzEGlq/8XujqNg8YKRVC9bQmsFaIxW4EPPjbohNGem1lwv2xzR0tTIcNsM4vXCenJVHyRvcFfioseuht/QVQRCEf07EHOMniz0AXMPjrEoPUYar/9dDDFowdX1MKdRCLey9WSldRLqIh9N33NOGaEMwxSsdfsSP+qDZWwtULFDx6opjTkV2F9k9aXMdzhN98qE3egN3GEnceOnPi35etH3S2RqWFZYVJmvmB1oHrYNeQh84RB+i103ws44TOAEMnqg+/5SZLT6nlVSGq7/qTR50iIylgcuH84cf59cuzIwJTvu/V3e9rNtFALj1w4A+AFD6B+/6AHCngrUoP7I7tntp9bH6WH2M3/SEuIqruHrmrjJcLaOhYjnn67t06Fhm+/r6+vr6vvmaLW92J7yAC7iA+EpqrEmJfndwB3eAtaiACqgAYBM2QQxY/1L0o7QH0KA6nAbVedyM5zp3gY7MGu3AS3HtSsSAFTcwcooLAEB4RdAA1JqpEqO5oL7EDteYuwZEN1yZkirZJLWS4xjyd4Mph41Z/PxMpZgH4AQcMCdRkahDj/UfLC7fg+d9Fa2IXzwaqMbTg97AiwoBI666cKh7Gd/2OE8qJKM50DbuBQDBLlcHAUBQnVdqQjiqnZMuIgiC8FFhsQJzfxp6GPp9lfPGXC5oytIzPpuGHmOXwLtc4EY8V/J31s6ycFfaAVzID2REromqy0gaqCINTzX3NbLSwzeY88yvnO9MnB+Gc14LZ82GB0ypaE5P+Qw6ZirSMDeXM1wTGth2s52PGVnWkrUWwXmxKj2OafjnI3QcIVS6iiAIwsdE43iORwDg/Qu4eZOJuw82YwTSWNZcvKppjbRGekTjaxX2VNhzdNgd30enHp1KYIlI1OZwm8OdMoWpSJE9jVmFeiL1ntv83sY836kSRboX6X7xEmD/vf33UfGiCl0efXm0d/b8u1+6v3R37lh+JtqiLXDvmPq8DzNNjk6IWW/yYmDAWO5K2JiGLz8e34kO+OB2ttfvQUPWqO8BwL1dixYAYO5hCeQH/D7sfgbQkdOdgQ4BLBGgv6bB6pb0OQEIZYT5EOr3zxjJHnVD7o3wARgphFsaAUD5HKNPqoFx0jolgxcqeTKDktlL/e/A+f9Se6lks+1Kvkyo5NmUSuZ2juW4h0rW/F7JTZeUjCqv5OSuSro+j/l4r2VKzuXxlqFKznBR0iNxzMd51ODn7vG4aUrO3AgAOV409FCf3McJwPuu9BVBEISPSU0asPYzsjZdNAdDEdamulrIdgTfz9S+Joxksk9ve5xRtL0XHR1h9OiFsyZHLx5nd932OBMX2N24wAqlIS2MqRjzGKn1jSnm60lBx8umX23be5MKxLeWmI9LyqSWtQ1sjxvKSC6klb4iCILwMRmakPI7NepmukK9Y6CSQQuUPJ1JyZzesegxK5X0ecrjwpU876RkXi3m41C6dbrW6db69KzQrEazGhueWF1dXrq8DL1+KZf6f8E0sZzPTskqKlNEX1eYelNZJWddV9Izf8zHuxdTclpeHjcJAO5d2LgUADzSPWJETOk6n+tJGBFYSjptU3JhBSV/O8T/p30TqSUIcQCT3II4wApTCwCo+qoNi9K2Y+rgTXq0fXyUvLA75i/w2afk2KFK3mHEkg8VjzOxeJB/nqzkDG5XXp7F2WYwZLdHPyVDEtge50kFZixTARsyBdDP2F3qspLBT22Pcx+g5CgqNk0ZebaQqSSdWwDA+XgveB9M1XlghHQSQRCEj4kepmQ8Rt567LT9/yFGKDVlEfQbrAVVjJ7kSfRot+fugRpT/SyMXBrNSKchTOXTaLjqR0NVFxq+zNwvyUpD0Xh6sAdwnrEyUqshUwbn0KFRONr88oB/b8mqIpv597SMJJvNVMbim22Pe8zIsLbcHGQv22uSeUcQBOFjksv4Qfk1CqXMxEjb2bU4DlPvuMMUb5/SSp67F/MXNmKk1nh+7gFrFfpwd7RTeszH1Zt/bNOxTZkdu5d3veR6Kezpox4pJ6ec/KxTa/c6Y+qM2T3i2O2Yj6vJ+WMWI68qM/JqPvWWLimUDDphe5wbaz+OoKOlOfWhJX8CwB9lO+0DgOCbLzlPafbSWQQhdj5SCmE8Fv0q111JVy4ErcspWTxcY00jE6vbv6ktYXhGjaLyzw3LDGUXLqC5G9CeeF/TQzDXNamqVytAhaI7Q1aX0zN8u2+0Q7iwbsIaJOOZ2nGbioYPdyc8GX03gDxKNOAuhJNoS39QSsnGXMjvYa74q2gKjRc9AuM4kTTidrNzqKh0oaIQHK1WlQeLwI1msd5mzMVewM91nAYAJbK/5HPXqMDogYZmIq+qIAjCx0SnISoVHRlDuStNcyoUj5IruZ8OjYpcqPtw3u72gsexRhYYUTyRqX0RHL9/z8vzjVPy14pKDmQtRW2akmM47kdy4T+O25gc4N/7cd4vy5TCec6cBzm/HOT33GM7mzOVcAYjqSpwvvHj+qMpPd17ufvuA867ezhPOuzgdSWQviIIgvDPcaa7JPSsVhoASm+uwWpMoayx0/Oxksu46/itadG+gpHBjXjchIwc9zkP+LDo7/GtMbegHr9/0ll9tD5a83B57r3We+3TE4fm977c+/KCNaf1V91fdXdpDizDMpR6c1wtzndTOV++4K69rYspuZM1f4Oi7ZLuToPcSBrSWnLeWczSKh2cAWDYQ427xIPzHRpIb4mJPEwxzWcYCJmxY2Xghk4Do4n2BI16taFfv8EI4GGEt0692qg9tp2pmA/6yD2Pm3ykCCwHvpDefEEzM/Knm7H7AjtENxazy8yicmk7U3JgSMvtrzOzFsZPHNAmM8LoF3ZQx+CvSo04hjsAsKbpFC68Rw8AgKhkt8dFe1z0KDSlHM+aH7dpEGrIGljRDVfaPI6HTMGYxN0GHrJ2VSNOEOv4/+Dohiv6TMbRUNWQxdtnc+LoTM950HPb4zzoOR/N62rKAX8Ba3z50oD64jkAnOB05DzLXBAATHPMzM3FVXlVBUEQPskMxIVzRRZvn8HUvhR5bT93lbVF+rO4+jB6oO240Bu2XcmOLJpuxxockYzIGk4D0iDOA9oR/s6iot0YyWtfQ8lwOsD2c73QjBFfO+gIycTIYr8eShb5w7a9Dzg/tqDhawtlel6XH9cf3zJlUWfNSY0edFi+7HPxYE2TH/Lw9/jSVwVB+Ddj/o4ysUmNw2UtjMjtQwfKiD1K3poVTY+hntiY88dE6gX3vKiXpFHy+KRox3EeqEc9dQqLuz+unuVVlld3E0+IdP7J+aeIbSGHXkW+inQ2aUXaXWx3seObYvO1qBdNYwDGC6Y+NqZetZI1IgOP2p7XnfrVSOpBLQcoueiQku1HAcDI8Becd0zlpXe8D16cx1OzVEBJ6pETmEk0g3aCyrzvaTZEszMYkgEYhejo6sqMp1nVlSxRS+513OYjGbCecAAZxQGiNy2aa+rafo7vLfrTwNWPC9L+9LT2pwe0d0El6zKyaDg9rs492YEvfE0PwZrXehIAtgRcOgYAmmY9AQCmIm8+wkisRtzGfCwjr+6xiK4PB+4znWI+Q/3btgP+E3rAfegxP5E85uM8+f1jWIS3IYvyzqWC0ZUDRfCiaAO3sWviaiWNFJRF3CWqIz0ogYP+96iQVi6bAaBnwsoVAcCrjarxuAAAOt1JREFUlzu3qdXFAy4IgvBRMYroPjdzYc7xthK3DZ/GoubJMtseF8nahMM4X//OdYQ95+fB3HSkPReEJiomFjpCRhq1T+jI0Jgy348pgZ0ZAWZ2sz3vPRrMmvN7djFCKiMdKLP590LjbI97EKhkCxratjLCOC3XKbM4vxZnaoeFqRv6F3aUfcNdgCby/hTOLH1WEIR/M5EMazBNN00CgHUnlrD0ySTOG3o0vRF0ZDcsruR46i0PqYf4cJw8VTXmM/5EB8VkGiT8mYHi41tsa7Gt57Mc6Wodah2qpTRtNTmbnPU+OKxpmqbrNX5Un/uDkTlBHH+NTUEOpor5fG6MUB7OlPkWjPhdQr3Wl5tXvWgOAGM5y2YJVbUjtd4m7taO59JbYmIPUzH706E2iI6m29z8xcL7NpOOuH7NotkZDMnAmA7sfz8ZEdmM7PbuIfc6bvOxamDR86nTQh21VsnHpWw/ZslGyZB9namHVnY0K0P8LFzYvmBH/YO/X6PhKt7Yr+w5cCA0xQcAe/qF4Wj8uzE9EePpwb7HF6wxI6VOTov5a+sz9W7SbiWfMPSyMQ1RR/2A/99V4//x/F7JsYyI8+EAMZcTQGcW4w2KFqLrTs/7KCosLZnLvZj9wZepHC+q2R7nylRG1by8xZupYoh19TKGBiKvqiAIwsdEY2j8PUbk9qCj4VYvJSsxVXwaUyaS/mZ7fARdLEPoQf+N32fPXf2G0rHRYbWhslCF4QLyd0ZyDeZMZ6JndSBrKXaprqT5ou15bzPyt2mkkjup8GRhpLcf55lvonlQ73E+aU5D1hZGdmWkA2c2F76F2T493Zd5LmYazsqxxmQm7hb0A+d/8zTpu4Ig/BsJZ6XEyOpRtwDgYvL7pZUeYikd8xENGbM1gQagh5ynGjMV8Pj8mI+rS8fEFGasPDNS4LvU96/vvyP1seGWLpYupsemHqZOpk76c2vNQX8O+rNRtZrD1eemMqLnJY/34Xi8v3bMepMbS92MoAbHklhYytTG9geUDLCpbfUwpcsWAFg/pNczAPAumpL6kLWs9JaY0A19MESJQP4evITrGubyWMrZHvfGzmBIpnxG0vF1gYbSiazh7Ep7hTmx3PO4yScu4q7tj/Y7DVz4wFD456yxcYERO+6/f6XPYx0ARNI8Z1rsQ4PQRG4bfpeKREPej+OR0e4vB/b6tOlPYcjbI4bMNqTB6/Blm+HgTYlDT/5/HA2MPvS0+3Fb9E40RAXlsj2vOxWM0aeVbMHQzIWsidXeyBlPHW3AZ6rG71Q8WpcFgJBatyIBIDjsFWcA7ZC8qoIgCB8VLgDNNPRsY0pEky1K3hqhZBUafmYwtSN5tOKy4XTADKJjYhiX9vYzlBzOhXhH1ryyY2pcBGsq/ubH429yvGcK/GB6zrvRcGbf0va8tzkPNqFDZyfPk5Wp9PM4bxRJaXvcvdVKNtul5BZeV0bOk+WZmqi7fpnHkpClE4pFS4UpzVSTBOuk6wqC8O/mTar2RABI+iYBT6NjoSEdHpMaKnmfmScNOX8cjR6hRD2lLh0yf9AB8ZTzWSP+fnBHkctFLl8YCzjucNwR+TqqxK1St0ol7ZN3VPCq4FWum8Zw/A9mSnwj1rja2zVmvcmNkcsjOV+0ZqTWUtayasf5McDBtrmujBj7LR8AeNbxTQAAdvPt6ADCaOkj7wU3kUGfaHaNQn/v68624bqGEeZ2M+QW20D7ghPXXWant3/cgQZZu9ofuyH/kl0ILSyK559HSY8iX1d/0OcBgGmLXkMNjA1pYBrPhfc9hkA25vWfiqUmVE0u8CfzBXxMC7RhCDtWJebjPOn5HsNc4UbskEbtss4+HNDrReuYRqgsc4ibsZ2LGJHXsbqSgd/ZHmfHlIjB9Jy0YQjnstsAMNZl8BgAsEyNMnYdaSdjhiAIwieBKYQmRjrvYcpDMyoEd5jSUJlFZ6dxfkgWLcIpkpG9Q08r+TtT1e0G8O/8vVUZ2+OiOA+MoKIwNCfbQ1fOwNtKdurE5vrZHn+Xvzfj/Ld7t5IZWdNrDhWSXJdtj3vAyOYWdBBtixbphTRf5nEUZKTARRclrzAlPzNTLws5SJcVBOErYS8APH5jbmhAA/1ERtQ8pOLrwxTCE9Nj/pofGbk7hQ4Jf6aYGUXdD782Pnmzxs0ayc3AVY+rHt6/pAgJrBVYy21rOiu2Y7ve/ZmxiyHnlf3ZYz6fK+eZ36nPtGRR8aWNlOzA73mRJdp0y9I4Azm/tSsBABE/7CkPAM8XPDciidJI1/gSPGGAhpkpoHbX5Z4AQDK+dz/RENyUhuWOrG2ek/3VvFjJkozAb0K7QGvaIUpxt2jTjn/aItO/6wbuZo7ruSRfV8fQpwJAp821aJEfzwHxIYvN+dCTfPIIEGPKH4vfNlH70sKLn/+Fhr+jW2M+zoOegtFneTxzu+czkqsLt1EPqhNzu/MxAuwndtiLHLC70EIbEEvqXy5GctUbpeQ1Kkid7gPA0Q3+NFyZ6GmBbGcuCILwceF8ASNVrqDtv3exZmEzepLvsWhqZaa4TWWkbZIfbI+L4MJkKIvbjqCBzJmRWz70aCZvbHtcJOcNo6bWUDo2zJyXBtJH34nzjqmN7fG3OU81Y63I3YxkykxDXD2jJtcp2+PusR2tmHpwyCiq+5lrkNjzfhdjja6VnB/304HjRsPeD4xktuskXVgQhH+p3kMHg9o1TtfrsabVBI7fT+gg9wkwJTYl1jsfPzxhwoQJEydi3uJGixsNcXNa6p3QO6F/xfjUc5owgjgeI6K6D1Dy8Jzo+s9Y/7H+taq5Tbxw8cLFNGVK3dV76721LKYB8dvFb/fKvdepYimKpbjgvv/bmNvtSkfM7xyfW3OzkeWMQG5Px8jzWFIbczJgoH57JW/eB4BVRTpNAIAg+0CmwGtp33UHJ06cOHHiRGDiwIkDJw50vD3h9oTbE24nCJuQfELyCclxd8KECRMmTJCe9mGE0IC4lHpxeNf/9v1IOZX9lXaGZzRYreHfQ2iQGsV1VBsjEp7rp0O0M3xDx9sovt9pAv9py/5lBqzH9Pw+PfM1dQ9tnjkEALItLEeL7316whtwF6djof87AP9/6Go8RkRN5IBYgbUyZnD/jN0jbaaL/08V5EJ/bH8lmzJVYx5ra3RkbnDg7ugDv6IwQ2LnMOXDvIIDOlM/nh6I+UpZCx5zT1OhYS2S31er8xjFFbWXMogKgiB8UgWCBpHU9Hjni8XTvIMe5CacN+7w96p0tEznrn5JV9seF8aIpyk/K3mRnvECnC9mslh7ig22x0VwoTOMtbUG0YBmpoIyhLsOdaECES2jETepUEyhTz+cjpgu9MD35IzmEM0Rdpse/h2s4ahd/LzPw5uRVwnSKHmEKYSbmToTNU7JUjTApZgvfVgQhH8TpnHGT9YuAFA/QyUqwJOYAvaItasa0tF9JHjq5amXx9wFOnTo0KFHcYwLWhe0zvWXmn7lSpYreSxg86Xk+ZLne96gYrB9W/u2Fs2PetGOKTHrPx4sIz/m28fxH8ePf7rmb8leJHsR0OjkyILDCg67bDrUMUn6JOkDZkdvuZsRKcz5rC31lGWcgNqy3c+rx6w35aEjfi5rLrrREDD8KAB0+vUB50HN+sE31R72sPe+ihu4gRt/zMRkTMbkbotQH/VRP+3RSfkn5Z+U37Rr4rWJ1yZekz74jnXRZiWv0vAS9R+9D45Moa3D2m6H+X7uYI3sB7QfHWOEVU6WhqrKSMRVjIh8zVrYRRiYEsVa59aD/7SFdl/3A/Bg8XF3DiVR3H5bK8XL9+EHjeKB9NTqtDSGcCEZnPl/B9yPjXmEyR8AnhW6ygfqo16gP08qg1TN6ANwPBqsxjOV8GeG6s1kecSu9FC/zB3tfjCndzRTPXzY0eYy9aITa4wFJYl54C/I0Ni526iw8HPtuR36gmcxX2FeKhBzOVGl4sTSidc7pxYAuJTFFgB4vU03tjV35xdIyoQgCMLHXagxcsqLEU7jOS8E0/O4NZqBZwdD6w2Hx2wPLliM3fp2K9GaHrhHnCcecEHemql8c2gYqsBNQGbSMdWCpQHuM4I4kga231ibUWNtkD40pA2hx14zHDL0+EUxAmwl57P0DHUfwBTI/rxO8HtHsJaDhbVTMJey9ud9HvnZnhvcDck/PReILFp/gQa7HK+ULMD7fke6siAI/wqcqXW9nmAuDQBpr+ZmKvrFQCUbeyl5/IVxzKCUg1I2Pgq0fN3ydde88S82y98s/8bcdb7NmSpnqlvLCm5MlT9V/ic/Rix2uuZ0LaIFIt1yueUKvV1ol/my+bL1zhG3Dj079PRtacoXMShikN2VUTRQNCngtNZpbfgv11Mmupzo8ovmdzpVc6rmdM8SnOFoyNGQrG/Ua1ejliPH49Z0xCyjftmem54ETI1Zb8rNovDzeOVpWfurc24A0LTZLB6ewIj4jf/BN9UOdrCzT4woRCEq71LcwR3cqdUGIQhBiE9za0drR2vHDZ2xHuuxfsVxFZl1iYaIEO6KqMtE8vc0eAZ+JOKmAiY65qJo8HnFfv0619u/Jx7XPc6MFLdynWSmAy+YkX0vj3ye68rA81hpXzhaOObPuXM3Z09+fj8Db56whFEQI+tbcv31gnanO0btsr9d6+0rN2DlY6jnj0wJKMsaEsZufHtZxDaSu0aY+nEgoYHIK5OSp2hZX8gFpP/Dj9nKqE2WUQAwIf1UXwDoqZ1cDAAl+XR2vynV7kmL/1gOrT/T8DaHA2Q3plS8rBStgxlF1mnQamIYnLjLQueXMRmu/uc+GsXcacBKRo97e0ZQLTQi4qIZsHKxXX6MyEpNQ2AnHu9nNuYnAIhkBZR0R9xHA8DDLE7TACDsqD13o7Rj6oSxy8Rf3yS+cMZulowww4v3exJGUUk7o0ZLHmMuiuUAw/VPhSOKkW+6+cN6gInP02zU/Lr/rjUAz8NixVFGO6q85wk5UNrRtaAx8g9B7ziOiqiFmzFYy3zgdfL9M9OgCZd3HGAUUebno1hLTd/0nidkBIa5O89PBRpP33EcFWUr+7fF6N+337Mf0UBrZ+wWZhSfDovtAEoapKPo4dO/fc/r5H0009BrGvqO8xnQcGHs6hLFhQzav+d5WWzV7jQvg+Mlgt9xHMcZC5+jdeEH9iPeHzP7P5hqAMs75jkuACLnRXu/3wXHTTMdAKbN73mdbJeVhhmL94ddpzaMzTfGn3cUy3yzby1D4KOq8/l2e/v7ZeHnkzMFYgY9wy0YkbTduH5G3u5m8fM2fI6zOP5VLcXzcRxuxW2qn7Af7uM824zjjh/7e3n+f5ofz8v2PuEmIMYutCPpynHi33vSQDbQ2HWIqZDjWLPCQofTaM6TZhqC+jNVvx8XXhamRo7jOGoy3iNep72RCs+FmM75KypbtHHqXbDf2hkGOc7T9nx/8/O+7TQWwrzfj9n+Q1SAcrM/fc/ntpKeUn1WLP2Iz9WOHlFwV6xYU/ON+YSpjFGMuNPnved1clt7Y3MA4zre+b5wtyed41cUa3Bi9nue9wmvk5F/2tX3PK/xHnB+s7b9wPHoPK93Q7T7p7/9fDoX+FGMsMPQ9zwhI0jsqHBqfC/w6h3H8TlY+d5Ypn/geMT3xI5FscH1IMJjOcAYr6jgR3H80y98oILI+cVkPM977ziAhnUrU7EsVNRx9D1PyHHL3sd2/YHI2BpIST3DwnWG9QOTt0x0KJgb2I4377pOneeNosEfP77nCanP2PE8WraY1/GxrTes3OTD8vw93zOYHEwOLv0BPZmeLLK4/TZTTlPOZw/vNXSs5FgJKVqmC2kd0rpwwxN2K0NWhvSYZ3eg4IOCDy4m089UHV91fN2b7icfD3k8JF6dkdXvjb83PnGSyu55wvOEX98BxFsVb9WrPA719En6JPzo0yHqYNRBc45qbtZm1mZaof1Vavet3Xf3IcvF/S32t8i5r+pZ/yD/IK+0S8uVblW61cmc+yymu6a7lsHfXz1f43yN+3nst045MeVEgxuuynBQ8TcGPrShIr+iE5dJvA/PO8Z8tTn4PvttBwBLo4xeABAa1a00AIw6OJffax4Af/jD376YllBLqCXUejndc7rndM/c67X3a+/X3rg0YeWElRNW2vX7i+54O+p21G3rKLtLdpfsLtlNQC7kQi5jXaPlx3M8x/Psp6BBg5atJ8www9yYm2odTw0TTDCt26/11fpqfU/O16fr0/XpVwbgMR7jsUZHlNYh2vpCsMGN97sGA0JKc71ZmI6+RTRg9eW6JyoWu0tJ6rk+dOyloj5/hOPtPP79kOfnua7nnH9XUh+IjMXuYdjlIjhuHBth+/8w6lM7Yhl//j5fuQFrNx/8K+au1eEuSEtp0BhoWJwNRZwdzpUGhTz8+0CGwOXlAqkzPbQBjT9KM71RDQAe4tUNm9Ybyxp4sdjfeHaohsyNnk1DVBd6mIOj5Ux7UIEezQ7fjArGfE6wHanABuYFYkpRLEgDkB898imogLRncfn5xj6xc23Pm9fYzpwRV2l5/zryOcypCgANE6pItwXPzYUA4Jt0BaoDgP2hbuUA4H7bTFQcptOTHsaFuNYilgUWFZsdTIUZy+KN1vc0YKWn4jCSC1GHi+9YMHFhHc5t3Cdw4bv7AztAXRr4GvFFt4x4x0KSHpu7XDj1q/1+dhmDZEylGXROyZSckK1J3r5A05nSs4iGpEUfeJ3f0SDTkQs1U9V3XCcV8AAugIfRgHnxPQ1YjjQEdu7LCcIwmHzzjmUWDSobvJScmsxYLrzfefPzQfzC/u7eg/cvNkWXCmIkU2h/T2djZ3y3YkFDfZPxnEiZ2qW/w9BuYn+9Ts/iAHoQ37f0TyKOC78y1D8DJ1j9Hbu/mPhcJnPmW/+BBqyKVCRa8Hg7pjJjaiwHcBx4RUNtZ3pQ39cPYc8Uu15MTfiGKWr6OwxgRjHLg2zfoA98X/LQgDTYiDR6l2Gc5wnhvDacCurJ2J4D57NbxvhKT1hhGnb8qIBdouPhNQ3BGscBne9vOHcJjKLhqxoNAUaNk9bjaF/g/LqH646mVPRncUFdkeP6ZipUD+jxi+L/0VQJZ85rgZwvvXi+wYlsDdAT+N5H0uA0goYeKxXE/mk4fvK+Vue47835K5T9JjvnL40GzIc0hA2mIv6+jmtXKiC92U/zsZivHR0PWfjef8Pn4ct5HzQcaXzfIunAKMlxIz0V9euxGLCy0cD4KwcUdyNFJRaPKqhYhbB9g2lAOPe+7wsNMx04D3/Pfqznf8e4wI56mgaHoVQIXr2nASsVDQ6/8j1NQUeLXvEd4ycNMbN5/IoPfE/r8b2rx3WVZhi4Y0lB1biOesx+15njftB7GrDcaIDsx0jCHNQM9H7vGI/4vmzk+mHCBxqwcnOc7kXDupsRoZg8luukITiAnx/I9/BasQ87by++78W53rR0eMfzZL85SgPxMB4fkff9zpeUhrexPbmONopw14rlAPYfCzczmsn3de0HGrDKc53Zwcj4+OUd10nH8g0a9gZHfdg6MD7H5f4cxzIttDVMxXpe6kVrOL5Pp4NPz/auM1ojrBH+4wv1SNIpSadkraqurfa02lMHzXFVxsoZK4+oOGTgCM8RnjVPmqZgF3ZdS6e/unH3xl3nPHrgj4d+PLS/h3vVwpsKb7r0Xe6ViXsm7vnCx9Jaj9KjtJqmYCRHcvSBBw7hEDRMhDe84R8Pehm9jKlKFRRsULDBFS8gS9EsRe9VAsL3hu/1X5z9TMjUkKkPzmQ7cOPKjSsYnXjIrC2ztqzxXHQNAKrl8ywNACmRl/1ghbGrIN/fZzdj1pty01A75w4AhPjn7ggAfWo8NQPAxCUVXwOA9U45zptOVcyFzIXMhTzzlSpWqlipYikjyrYv275s+1Euk25NujXp1oueWiGtkFZI+8v63FzYXNhcGEH6d/p3+ncurXEUR3E0RXplmPrLcxsOHTr0JFAGrUqAFVZYKxzWF+oL9YWPhykD1429eIRHeJSxrPrcVMOh2VmMVTERxHnUyFg11vvZ6Kiqx3XqsjRvX4+t5rzhRn06Ldfn47geeen0ea/LcJzFOl9zvZif66LnHJ/PP/xcLbT7b3SwTNTI7KkAH/8mFoWHC14GZuEAazoN4sS3mAaiSpyYP1oJCkOx3Q4AZVpBeVIae3GAG0/PUQNaRmZxwdKVC7XgaAsXDxrWxnACbDqA7eVCqSMNLoH3YjFccYKaR89eMtbm6MAF4fwCMV9GXnqY5zHCLTUNWJ0yA0Ctu7MTAMCfqe34Ymaix6hRFwA41az6IgCwbnkQBgCRezrzBXpsKKg0uGFrLBMrn9d9Guqsxu4RKd/vMfgb28VzgW8yImhi8WxqXGhb6Am79jct48e5sHxFw6l16zsWELTUv2TEy8sPLDYcTAPHYr7/rkbkWL7YTkhJBeqqsRvlgQ877xUudKb9+vbn+Oa0tNCHsb88bvFh54uiZ3gbDZgX+P7r7zovPdV3jIWp4TrI+n7nvUtP1WwueBy4QNRjUxgY8Wnlfb51+sOuU2c/PMQi3E84Xul+77hOGhCCaJAI+UAP2ysa0Few/Z6r3u/+mmgYPD/y770v52l4mc4IOW01/1EjluukAh1BzTToAxULC3dV3cQF+qmV73mdnOAftfp713mf79lUGu7h+47nyeuMpGf43rvmd1qggzm+ded5unETj5qMQErE+/wHDSebaOgxM1IqihpTMhrsfqUiX53jtcYFTWtG5D6mQrCLhsgmNFDNZrHynIwMSkCDwiAaUu5xHjPxeRu7Hn5Hxc3YpXAoDZUmPqdxXFhF8j2ZSodPXs4rtWn4ycB5fQXfi52cV0N5fo2GzxCeLyDgw55nBOeXdfyeQ2x/BRaJvcGF7yoqDJoRUUEDgIWGnMSBNIAx8qIw/x7bZkmP+F7MocLr0Ovt/dcwiEdyHHzIyFnsfc/3hfPRTjoqri16z/eFEdD+7HfhHzi/vOBxi9kfXWu/53jP/nepyt97T49x/g9ysTU4Ipb+odE1+Zo1S0PLfdj5wrluXc1xZd+Z97y/fH/ujPub4xEVsdk0QDkUfkc/osIbRgfJ02x/77yb+b6enfye6yOOD09oIIzy+MD1ER3ffjyfg/aO+2tEKHL9d3HI37vOC1yXTcn5nv2W820QI65ePv6w84VwXF/hpqRXBZ73HQZGjfP2LY7f+gdusnQzZ0ibkDYvW59Ke3bn2Z1HrW5nE4QnCHf7fZm7HqQHRR5DzxfeL7wTnvcqetR61FqqYQ8t/Fz4uQRIsTNFoRSFng+evzbHohyLbvta2ujt9Hbo4uMBF7ggHNCua9f1hgBWY7WWEXjt/trdsavlhP0a+zVRWU+Ojn82/tngH7cMvt3qdqunOy/Xt+tt19uhaP4U9u3t28d7XnqnvkRfMs0r9T0AKJ08Ex34y2hYbmtEXMVmuOJ8OI/9Ne0qAHj0uGcNABif9LSav5em4XOuQYN+oT5WF6uL1eVu6aiTUSejTg675zHCY4THiDutLHkteS15cU3boe3QdsRQs9oe9rDX02IbtmFbikR6G72N3sZ7AZ7iKZ46+yoDFKAMVwDiIR7iWVfjJV7iZXhHlXKIrniBF3iRJBxJkRRJtYLKsLXCSxm0lq4RI9X74MF1sBv13iGG/k6HXA2uj04aNZ4Txvw9ienAOsRx6+X9uHm9Cbi+zMN1xVmjyPupWMYNI2OCdgn9HxfH1z7tBXbkgmccIy+mc0Jpz8iayGSf9vwmLmiNFILK9Dx9x9zUe0ve73tSUIHZQsPXcSrSPsf/WfvcuXBZuQUAzBOXVQaAqPYr2AHGUBFoxAWsHy1mnfkiBKeK9n00HI2iR6w5PYALuYDoSAXuxeZo94mhf3k5YM+lIc+bCx5fdrx5Rq50tFpguZizPZ8vZDp65jueBYD+F/0eA8Cg7Fm4cG5ABagKXdePOHHO4wS8gff3xbW3K6aCIAjC36MmDUtd6OFrQMXkFcfl6RzvqzMC6Rbnu2ZUmHYzUkCPFhv+PVO4ZnM+8KYhy/AwtqYD5Um0lOVSnF/nMKIvDSOp1tKj3YrHPeb3g54/B6Zw9KULtBcNGZFUKAbR8H6hie28nYmf30/D3BzOr6f5uYhMn/b+O/C6xnOefVadBhE6tMycz41Qe3s6aJpxPv+B6yc/3o+W/J7IpdK3BUGI2wyh48REg3nvvcog5E7HzQhmChiR1vPpUOx0emLLiS0nhLS8hNu4rSUd8aueTk+HGSgYcDDgoMckHLvgc8En9Z948CTbk2zxKy28kLdm3prXzV0HFGhUoNHlPc+OP0/5PKXXK0ukQ0KHhFGja1y5MPrC6DTh0/vMHjR7UPlsXpbQuqF1HZqvpCOnHQ3w/nui6U00NGTnvDSPAQLplSO6XtdUAJAix9pOAPCwb/UziVYmWploZcOgiCMRRyKOuNwJ2hy0OWjzIjPO4AzOLON8dM9wtH7zrjs4cc7EORPnAAhFKEKzDNUz6Bn0DGtn4yIu4mJGejQiZiA1UiP1QwuqoAqqvJ6MOZiDOW5d1XGWacpgtX6V+n2+DzzhCc8zPyoDV0RrX19fX1/fuNqPEjDSaQ0dmwXoIKuaR8mtn8kAZxgwK9MgPJkGrEmcjwvREFWZDrWr0RyMntTnu9EBN5sOgltnv9CNpWHKxIhzKw10Rmr8dyWUXEGD7Ug6zn+LJcAlD9eVVjpgzv72Txv4lUdguQSyQ/OGXaJl9OndD/wiegI1evaScmBxpoc9tOQ/a6e1DQC0zVmFkRodOIBn5QA5lwN8l2qxGK6Y2jGanp/maWi44uc7hsZsuHJl6PEQpsjU5bbjzowE8mUo9jx2VH1gNMOVH9vH1Id01QEgaoAKiY+/8PAZAHh5fTCL7FamgeopJ6YRNKBuouLwzMihl41fBUEQPi1GajAj6bTdXKhzIdZ6GT/H8bk6FypLGem36bSSXRlp84wRA9sZkdCc3zeLEU3VqRDojBhqTUfSUxpgdjPSt1l9LuBYO64qIwSz8vPj6DD5g5FKEUydHUrDjjND8LswwmUAI4mvM/VoHw1fI3idp3gfIphiiN8+z+1Pywjf+KwxdoXzc2qmLoIGQqOmjVHb8SznyaKMHCtO1SI1IyGuS88WBCGOo7Hmn55eGa5+pH4wmPpPZgYYLKC+0JnjedBs0yXTJd3P2jNkTcgax+L3Qnd57PLI8yxhr5NLTy7N5OiMgEMBh9wDFh2P7BLZxe64b/he7EUuBBzCH/jjzemdXFpVG1tt7IHvK01OkCFBhuCXCZbb97HvE3Vi+dNQhMIB7RgI4B+ttqYzI2n7MaW9YXUaIAoBwPGzY4IAoFgDuy4AEFF5BSOtzCFFBxcdXHSw/ZqCOQrmKJjD+Yq7m7ubu9uTiu6Z3DO5ZwptHW9hvIXxFuo/Bz8Ofhz8GGhcpXGVxm+LC41ABCIAhCEMYdoFFXEVPl9FXF1splIJb95UEVdJ42EJlmCJ46+quPv6g+r/c39T8nRruMIVrhFbfNv6tvVtiy3SRz+EjFz3PGHR/0BmgMxhLcmyzIyqTwPPIK6jrDQUJaYjUWPmhj9LteAzG7CMGoc1GIBUkpG+o9jeW8wIK8SIelfaQS7GEnnlSDtBGUYw7kj10Vr6aW+EHm0XPJ2RQ+9bDPmfkp4Wz7T0NG9iJFH4BxaRdGMEWWIOrKdpEAr9KBFk2iyzKwDkWlxFWfhL+jMUz7DA/sIIpODutkcaufmjuDBvzhdgIVP5OvA+B0YLaXajBfh3WoZbcUF/m4aonvRAzzUMStEizfJQUZnXAQCi/DLWAoAt3f4oBgCVn6elYaodDWPPqfiMoIdhE2sTBGSSQU8QBOFLoHMhnooGrBw0qBgZ20/oeGjFBVg4UxTLckndiCl6zky9a0uD1zMet7WUkk3pMJrJlJwfmSqocV5ozW2Yn9Ags5OGpCY870imyuVkDbuRPN5Ew9dCHl+UDpYMNKg9oCFrIyPL5vJ8xxi6HrXry97/b6nYHGJxqXE/R/vAiWi/MzUtAQ2A3zJyugCPL31ZDFiCIMRlTCzEYr2lapnVyVqBpQua01BlXq2kH/WTbkGOnR07R3YJmh0+Nnys/RjgZrmb5ZLl2uW7bOeynaW+r3Tj3m/3fkusef+MEziB9ouoZ7anAv5irG0LXBnQ8Fv8O7/f+T1pj6a1s1XMVvHOt2fqlHha4unZ7b/M1iZpk/QV/pZ17de1L3LJOM6FEcGDVQRKtY6tAODB/vBNAPDHyyu9AWBonhJp1DTjQn3zT6aGL0xRZmyZsWXGNjykPdQeag9/S6wn05Ppyaatf3Xq1alXp042eHXo1aFXh5Zzk46t/dVugbdp2As1DBpvivur7wH0vHpePW+Yl1ZPq6fVO/i9PlAfqA/M2RYWWGBJdVClCO5zUBFV80fDGc5wPjFR/R5VOG5HWL0TQ781AlOon+ssjfI3dnX8IMw06KRnjaoD423/f4CZVtu42Ud9RmgtZkbUFfaTjNT3H3P99CrFl7mdiRhp2J16flquo+YaJV1YCiUnSzI8ZaqjJVotRDPXL2X59xe0M5xPE9cMWEZxOe4iY6bHM4mxEs7DhS5ryDhy4fxmUw+jSOLlj/sg8jjYnu8IF9TG9tnvSz4uzJ2N2hVGrSDrR+n+PdXuKoHpbjOyqwEX4samKZHRIp/cmSM/gh9oxppCC5mS4UtFILButIGbL8RvTF1oTY/vcnqiO3PXgydGce0FtsfnKgQAlpl+1wDg2rWclwFg4s377QBg6vdl1PcsfUrD1ViGwq7jwvxFtugDsCAIgvBFDFgssh2fkUvjmeIdyNpPezlfPmXR80ZM1SvLCKnpXNDUZuSPzgVcO85Hz+gx30aPW3POS7OYklGdxSZ1pq63Yaj6E6aO72ZR4VKMQOrFEPUeu5UcToNYbaYqJuDve/l7HSpGJ/6vvSuPr+nc2s8+J8nJhBBJDFFDUiTVhIqhpi9+jYuaIlQVudfcyEBvW1+VutXWGLS+21a1ESIUEaHSpBRFCCGCXBqUIKJBmkkSmXPOvn+8z/ZVKvfjfm6T8j7/rF/OsLPPe/bZ71rPWutZfH+VNi1qbN2uuy0dPnee76a2j/b+PDrC8XSAe6y6n8D6mhnR0rfkNS4hIVGfYMX6kpIu+lYA8OwPXow/rlDiZRL3jeMcOlG1XiOueP/suvK9le+N6TuZlSL9mHDfwoA6hBUuBTWmMVqzonURn58x1hBsCK6ckZnmtN1pe0HitbFut9xuHTxeYHP2w7MfuvwZ+Bbf4kVYBgOA6dIHh8ymm03PnTNrveNwx+EHXUz+5qvMV61emG+HJCRl7HQstLWytXI/eV7nfsD9wHM/7nip34R+E/qNu3rQMMgwyDCo0b7yjPKM8ozmVWqhWqgWAshFLnKtmsEMZjDrzbipJ+Pln9l6n8j9LIEFBYe5b+d0MjUxNTE16TUDN3ETNyeOxDIsw7IOx1CKUpQmkDjYpBeVWaechWZWVQgrrHr8QS8hbTotO4UasVOpQYP743LbzbW8L5W29PGcjjUTS03pB2VoNXNMLJVTEieS/3cw/ZBxjIffZyLNha2yF2/X7fIW0Y84Tv8kM0PY1vy7F4fGrKG29GV+7hHsFFP5OZxZ8KMngRhFfqL6Qj0jsBx4wxlDralWdMjGaCJkrPjpT62JFSw9y+QXf56OVjzFvap++n9+LIpn9iTjl00mMP0RF86Sju4wthhk8Ya6v9njvFyMVqYTAPDZns+9AWC2cnQQALzMQsT4qfdu3Fy3UJb0TecNagsJqVkMCArCavzA6OguYW9woJ2wMW2EDeYFmJv7a+KqkjyrLrrTaAA4W7p+MQDs8vB8FgCifi63B4D0UxdTAMB0ayPLc+NZEZan3SB85bYtISEhUa/Afc3EzHIb5mDXUWtkMvfnw/xba9WLp32dlcFhdMzGsMVNpT8QSHHpfFZi72Ol1DRm9MLZmjGSCRqV+2oARa5z6BeYmME+w0TXTRJeVnQMr1Az6r9Zca1VSFdoJe3P1K9ld2FCzZLrfm31Ix6AodxB+lP5rKDrvY/+Fwmtn+QVLiEhUa+gNUiZH9N1B4C82Cxt2BQTIkcGP/idNowvlvH1AaxU2sa4LIRHzm9VgzJja/oi7ichbEWPyejn1c/r7PxDx/Xeem9jTP/bpdWl1VafK31XeK/wHhNjcRoAzn3wgQIAiVPe2NZQaajEL9DvGHdx3MXnV5Qet91ju6f/38zLLU5bnDZMKlVtFVvFdljH/1IylUwlY84P6rvqu+q7OufynPKc8hz1bRShCEXNL1BcveZQA97P9bxxM/9/z/px37jF+Lqsgsc5KiqpjrZDFapQNbWheDyVnS4VJ//gFVY10Jui+n3pd3hySJOWByplhdM79G+6sHOpiOv6Hadsng9+POfTjoUnZSSebtUyhe8EicckEm6v0r/axUIbJ2qDx16p2/UtYWX9e/TfnqdtTD5nF/2MDLZEptC/6krN1BasdD/Djq5zrKAvW/m4z/QxEVhlrAhKY4b1Km8kB0lcFPOLMzCDamAvpHKEX/g3dFTVx3M+dlxobcrQeTKBP697tOP4kDDqyR/EMk2r4tDj/BLUTLU1AFxvfIfi6YgCgHhtGCIaUCR3BSujprGlbzMDiRB+3oIaRJEtM+JL2IsbRC2RaH5fgSTo8uyB/52qAeitAKBvilc3AGi7eM18ADg1zX0hAFywOHAeAExr1jEQ2cuWjRytd2E5JCQkJCTqMbQWvAxWxFZzn+zCFowIapBMZkBxaPf974+jozONKZavmCB5VdNaZCYyiBVVeayw+p6OzxQmvNayAsyPiRQTM3tb2UIyliX1Hejgx2lTgFnZm8KWkCrHerbAidyHuW7OFKd/h5lIZ0oZtGAgdpcOdlVRLaEf18lATc4yBma3SAS60YMfydetYUBYwfUs06aLtZLXvoSERF3g3l1IVW0BYOuP268CwJqkw+J+NVd7+l4cw4A+lImLAHbyRDFQDvqMxFWNihsbtrxr2oghJD520gZ0aFPcpjjbpf9IY39jf8VV2bdh/IbxW/u0nAsA392cFgsAd6fPaAEAqS4XdpYnlCeo5zZP0EfoI9yVoxMcRzmOctpY+hWykY1DeM641LjUGKJMVEeoI9QRMMASlrDEKu6vFLGeygqcoFG/DQahQq3IEgRUwWjxYCX3BztOr7NgZfI+xnURrMg9vUK8r7Ldk0VY1cRtTt9MZSdVKivG17KCvJT+QAMSfjoSKUZWOhUUP97z6cgKwutM4BnvPPh1hdQCXU8/KewGr0tqVV/gdZq9r36s8x3+7o7Y1njipRq/aCYOE1kRh6Df6wyVJ/MC78ILdQ8JtRW8oJfb/h9vpEPf8xthPyJTe4hTAFcys1ye9XjO894UQpY8bvsYAHqODnseAJKibenoriCxNY1E0RZeMCFk8AtqiMjb2Am7lOO+AxcIG8MMRqARAK5Py80AgDZr9dHicVf2jPt3BQArh1fWAoDR/3YbAKj8OIzEZDwd6IKhcjuWkJCQ+CNiFImSWSRa5pIAmUMHaggraa9xmu1EEiaHpzz4eEOpafglM54t+LootvIFM3GVW8PB89OmD3I/tGUL4i3+3zimcjbQ4TvJYSrVwfV7fa0ZgMwkMTWAjndrEoQGBjRahnI59/uDibX4C67CzqDj7s1Ks47U+iihv1JJ/+QcNUH+kXn/91J6TV77EhISdYuFvK/rmNCYux34NXFlQ0IglPFGACtctnF/CmJFR759jfsuCZ/F3M9CWBmyM5D3z26xFbEVQE75Oc9znovd/IsL4griCv4+t9HXNl/bbOlZMrNBpwad7OY7F2b3zu6d4bP240KHQoeSYWGHfI74HAGuj9x/af+lR/mkQstK+2sO47YlTGAUfypa+9LfRHM0R/NfliATmch0IkFlxX3jB+6bG1lJlsJETvW9Hp0nm7iqb1A4RXMORdgTSCQeC/nX72vIxNsmVrD/idd5MFsO1y6Va/tweEKmEOrZO21ga11/MvSWrBBKo3q+gb3EOk7Vs6ajqE3v6UONqxdJVG1m5ngrj/O4iKta0QUAkqKtmWkOJSOrZbijmDmYSQa+wKXGjZu9p4s6CxvI1sEdPH+NyMr1B4A2azsyABjH62AERQJz5wNAWc5i0YLxcbydeDyP40jVYvnTkZCQkPhDgw6Yno59Flv9p7MlL5xE0SA63OEkkCZReyRx4P2Hi2MpfwD3qTCtIouOuIkVXrOpTeLBRMxYvr6A5xPN/TycFb1nmDgpv/nHWt5yJpAiuA6bSVhVcx11nC6k+46fX3tj0wcfr4Sv28jKcq0CwUgCrJLfowXXz4wVbuVv0trJS15CQqKeEABab9wtAHgxGYxXbJgIWco4LYAdJtHU4AkhUV+TuLJixdVH3F+CWfG7ky2DAe0AQFGUfcMNww3A60P6uPVx6+MWcMyjqUdTj56Ndzp3d+7eosmJT/y2+m31W2Y/xfqa9TXrW+6ZZi+YvWCW0sXN+KLxRWNq4c/D44bHIe6OM6IQhShgpmGmYabhYT93CRMRx2KEvTxKiK23Kkc2spHtTC2vRMZjEYxTz1A0uqKzJKrqEpasTO/BzrLerPw6RE1qhc+rtWhtFjF+30Di9nnu5+m35No+4h3kyfgYz7OEdAw/TzttLOVfhT1HZrSEBI9ChtSCYn6l5Pyv8hZ0hM9fZ6m9qfo/c95aBVZMGgC0ctv9BQBknnelSG4AA4goluwFUyskf4C4EfNnci9j8Q4zFUtY0RX9JQBkuYd0BQDnC3YMACawlHEYK7fySfhtZOuH1iP8ixxnJCEhIfFEYhRbB//K/fHPdMCuUnOqJTUjw1i5PJitaFcWCDuZGhuHuz/4+EP5+FckVpoxE36R+1QV/88Jqi2u5/58iomjyjfldyQhISHxJGKRRtRfFnHMwgjx53JWWs0gkRPNhEoQW8TzXB8c/7zNViyGL4k7rgHAeyM+uAwAi/K9mHjxZ0GD1bMDvAZ4DfDK3TNw+cDlA5dX2+qG6obqhq74wryXeS/zXhtfU+ep89R5bT2RilSkFhuhgw66MySWYl4RNoGE1GW2qpV2ePDntaVW0GgSb/4k5Oz5uU/oYIIJpg254v+ksMCg0ls7giSu6pQuIfHkz1bWl1gY0oya05dJnEZSezN5wb8+XiNO+ZtKHmAr4++sgXKtH/IbeTI+ho4/dDNqZZi0MZq8EPRsUVDoqEMTm6V4qpGPG7UpiHt/n/O+R2B1AYDQna4/AMBs31a8JUeTcAri+RU0efCN+wVNE2QuAFSj1RgAWL1yWiMAmPV2RwYevuxpvUMtrQj+YOI4bTBvjvxJSEhISDwNGEUx8HcG0DGj/anGFJwWbMkPZ+vZIE59urZM2Ils6ThMrUYbalv2ISEVyta1lhw68g01HNdRizGZ+3Z1Z/mdSEhISDwNWNgDAJpWN2wEADl6bdhG8Glht3GYRxD3pby8B8c/nVsAgHot8iKsYZ39zXNmnU91PrXm7Kbhz3z+zOfH/+ay/47XHa+yS10crLtbd7d+5spK12DXYNfg73e1d2jv0N7BLl85oBxQDnRppKar6Wr6gWCEIxzhk6vFdED7WggFlQRUFvevFBJjB5sLe5St3c24D/4lQlg3vu4oW8EjKTae3E3bByVRVS/RkjwDKwW1QhiTHXkISvEYlfsfrxVOtKxEh6bhWSCX+uGgyCWoSwgCSynYOR4Awjxf8geAKde3sFUghAFB3scPfr/HBQBQ/7JxHABcdvc4DQCfmjLeAYDVTnd1AGCakkOx902s4Ipl4JGridqPkN+FhISExNOEkdwH1nH67P9QbHwhW9irT97/+lbcJ8JYMj+QY7gvc3z158xQ9mBlVmdO/z1MEfjIbcKeovZJhdx3JCQkJJ4imLOvo8r1wzwACNkV6A0Afx/ecIF4ZicD+SB2iuT2/PX7SycIazmh0x0AuNV148sAsPti56NoiZYRM6qe8XbydrIMyJo3KHZQ7IA7OUaDp8HTsKpsj7mjuaO5Y8V4NERDNFTdTeNN403jW3RGO7RDO8e2Yopf5m5cx3Vc7xQqKqEsw37zIVSoUAFBcIEEl6mZaAUs4fAsI1v/8khI/ECN4UhqdKWwA6YiRhJWEhKPDjO5BHUPvYsuCgDO7Ymn6OBMoQmyJO8YACicynGXtztzX48iALiQtsEfAGJjBXG1xaksGwB+GnRNAQDTS5t4441rIuwvtMijlQGEhISExNMJZrrNKJY+55CwlWz1C6XmYTXFyG/w77dYidWWGh7tqYk4m9qRe9i6P437WQqHqVT4yCWXkJCQeHph2M9dJVlfDAB2r9oLLcDh0Y3FM0ERwub6AoA10yPJkcLGxHl2AwCLhPU9AeD7fh7HASC59bkPkIWssrSolJthN8Pau+yapN+q36o/n9vRUGQoMhSZphh9jD5GHyVQ7aZ2U7sZ3xAte74nUYACFAzbgQZogAZh3mKa32oOzao5tLVqGJqgCZpkrhKi61WLUIhCFDb8DjdwAzeqrQTxtZcFAhGUpklhh09FLgmrGHk1SEj8+5AVWHWKe1MIqbG17Q1hw4bd/zrdZgAYHOXVEQBaWK4tBoBkQ4ceAJBmmTQXAEze63qJ1+9mxjvHRq6xhISEhMRvMfKEsF9Sq6oJRXWNnJr3kU7Y1WwJ7E4iS2uF6MTUShI1KCPZup/MlvTKLXKNJSQkJCQ0WFPRt7r70iUAMMDQIxsA4t4fbQkA6t68pQCgGwiKtzudBQAlccQlXWNd4xZ7Z37ilemV2b7KNct+mv00M4vDTW/Puz3v9nvv/5w5JHPIzdnHk+3b27dvnKbaXdp/aX/6vt+ew/1TAV/RiYqq0a2EXXVWEFCb48TzTTk17rwf2qIt2mb/CRWoQIVTFvKRj3ybQJSjHOWJnPYauVnYFA4lqQyXFVYSEo8fsgKrfkDT7OIUDQQJ48ox5BOmA0BCv9c2AIDxbzmfAkBFWAADh1hqkuQxwECpXFIJCQkJidqhI0GVzmm2qWx5mEhNyfnUSBzE19mR2EoqEnYGWwST2eJR+Q+5phISEhIStaGUs2hRYkwFgD1vnigGAGVTnpBOSXdkC/twav+ONwGAOqTRCZ2vzveXIdvzBx8cfHCqrduzDrMdZju8PswX4QjH7Pe9TemmdFNoTDdBQH23E37wg9/tePjABz7GtbUQSZNFxRVchTWyUvg4Rdor2RLoVCkqrFyiRathor8gvDYWQQ899MkzxN+Vp7UDS+JKQuI/B1mBVae4J+I+jYHAl7yzczz2K8wR3KGo/CaO2fx2gbC5x4RVs+RaSkhISEg8PEa9LezsF4SNfEPYYCZQnDgtKZ77zZd8HXkrlO+VayghISEh8VBoiqZoBAgCaOFXWIZleNfptrJKWaV8feSkGqgGqoGTnsFd3MXdRhHwhS98Y+MVG8VG0W/+zK+HX4+Radcve4d6h/a/u6xYrVAr1IpZV1CGMpQBMMAAQ0W0aA28ehM90RM9d4YiC1nISjBX3BV3xf3ku+pidbG6uCAAHvCAx2h2qkwicbb1dWFfdRC2dRTjM9HqiHUcrnWKUw3Lx0iiSkLi94cksOoUDfyE3UHiqsd8Yc0ZWFyyFjaJ0yxuUoRdYUmrMpEHYgYdRrmmEhISEhK1Q10qrPtQYQdR00rHabYVC+mgcx86Te3EspV8nTal8LhcSwkJCQmJh4cpyvmi80Xniz5D+67uu7rv6s5vGcIN4YZw8wzlE+UT5ZOqCw5DHYY6DP1xbsuqllUtqy730W/Xb9dvN6qmWaZZpln655CDHOR0+wLbsR3b3ZbeE1V/MGJERVbRYDjAAQ4/dhEtgKlZsIc97KuOIx3pSB+7SFRg5Z8RbzvKDqVI7nfJ3C8rNkjCSkKi7iEJrDqFDcVzVzCAaHxD2CKK3ho5rUJ5i4GDNmazXK6dhISEhMS/AWpeNWNroO05Ya8lC1s9loEGWyf0lNFFEm2ZXEIJCQkJiUeHmuAywmWEy4hnbV7zf83/NX/r/c0tmls0t8hbqKxR1ihr1H6qqqqqqnQ1pZnSTGn6aDVJTVKT4AUb2MBGsRUEVgd/7MZu7G495OH/NVSoFcXwhCc8ExxxBVdw5dLLKEEJSuwoqh45VdiEBGEr07W3S+JKQqL+4J+C684vxeL5VgAAACV0RVh0ZGF0ZTpjcmVhdGUAMjAyMC0wMi0yNVQxMToxNTo0MSswMTowMCvyY44AAAAldEVYdGRhdGU6bW9kaWZ5ADIwMjAtMDItMjVUMTE6MTU6NDErMDE6MDBar9syAAAAFHRFWHRwZGY6VmVyc2lvbgBQREYtMS41IAVcCzkAAAAASUVORK5CYII=\n"
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "%%tikz -l patterns -s 600,320\n",
    "\\input{fig/formulation.tikz};"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Governing equations\n",
    "The problem for displacement $\\mathbf{u}:\\Omega \\to \\mathbb{R}^3$ of isotropic solid at small strains reads<a id='governing_eqns'></a>\n",
    "\\begin{align*}\n",
    "\t-\\text{div} \\boldsymbol{\\sigma} &= \\mathbf{b} && \\text{in } \\Omega, \\\\\n",
    "\t\\boldsymbol{\\sigma} \\cdot \\mathbf{n} &= \\mathbf{g} && \\text{on } \\Gamma_{\\text{N}}, \\\\\n",
    "\t\\mathbf{u} &= \\bar{\\mathbf{u}} && \\text{on } \\Gamma_{\\text{D}},\n",
    "\\end{align*}\n",
    "where\n",
    "\\begin{align*}\n",
    "\t\\boldsymbol{\\sigma} = \\lambda (\\text{tr} \\varepsilon(\\mathbf{u})) \\mathbf{I} + 2\\mu \\varepsilon(\\mathbf{u})\n",
    "\\end{align*}\n",
    "is the Cauchy stress tensor (with $\\varepsilon(\\mathbf{u})$ being the symmetric gradient of $\\mathbf{u}$, $\\mathbb{I}$ the identity tensor, and $\\lambda$ with $\\mu$ the Lamé parameters),\n",
    "$\\mathbf{b}$ is the volumetric force density,\n",
    "$\\mathbf{g}$ is the normal traction on the Neumann part of the boundary of $\\Omega$,\n",
    "and $\\mathbf{u}_D$ is the prescribed displacement at the Dirichlet part of $\\partial \\Omega$.\n",
    "When reducing the problem for an axial dilatation and a radial compression we suppose\n",
    "\\begin{align*}\n",
    "\t\\mathbf{u}\n",
    "\t=\n",
    "\t\\begin{pmatrix}\n",
    "\t\tu_x (x) \\\\\n",
    "\t\tu_y(y,z) \\\\\n",
    "\t\tu_z(y,z)\n",
    "\t\\end{pmatrix},\n",
    "\t\\quad\n",
    "\t\\mathbf{b}\n",
    "\t=\n",
    "\t\\begin{pmatrix}\n",
    "\t\tb_x(x) \\\\\n",
    "\t\t0 \\\\\n",
    "\t\t0\n",
    "\t\\end{pmatrix},\n",
    "    \\quad\n",
    "\t\\mathbf{g}\n",
    "\t=\n",
    "\t\\begin{pmatrix}\n",
    "\t\tg_x \\\\\n",
    "\t\t0 \\\\\n",
    "\t\t0\n",
    "\t\\end{pmatrix},\n",
    "\\end{align*}\n",
    "and that $\\mathbf{g}$ is zero on the walls and points in the axial direction at both ends.\n",
    "The system then partially separates; the equation for the component $u_x$ becomes independent.\n",
    "Its solutions than enters the equations for $u_y$ and $u_z$ via the traction on the boundary\n",
    "(here the Poisson's ratio comes into play); however, we are interested in finding solely $u_x$.\n",
    "When we integrate its equation over the cross-section with the area $A$, we obtain\n",
    "\\begin{align*}\n",
    "\t-A (\\lambda + 2\\mu)(u_x)'' = A b_x.\n",
    "\\end{align*}\n",
    "Since we want to describe the material properties rather with Young modulus $E$ and Poisson's ratio $\\nu$,\n",
    "we substitute\n",
    "\\begin{align*}\n",
    "\t\\lambda + 2\\mu = \\frac{E(1-\\nu)}{(1+\\nu)(1-2\\nu)}.\n",
    "\\end{align*}\n",
    "When supposing $\\nu = 0$,\n",
    "we end up with the reduced problem for $u:=u_x$\n",
    "\\begin{align*}\n",
    "\t-A E u'' &= b \\text{ in } (0,L), \\\\\n",
    "\tu(0) &= 0, \\\\\n",
    "\tAEu'(L) &= g,\n",
    "\\end{align*}\n",
    "where $b(x) = Ab_x(x)$ and $g = Ag_x$.\n",
    "Weak formulation<a id='weak_form'></a> is\n",
    "\\begin{align*}\n",
    "\tu \\in V: \\int_0^L A E u'(x) \\delta u'(x)\\, \\text{d} x = \\int_0^L b(x)\\delta u(x)\\, \\text{d} x + g \\delta u(L), \\quad \\forall \\delta u \\in V\\,,\n",
    "\\end{align*}\n",
    "where <a id='function_space'></a>$V = \\{ v \\in H^1(0,L); v(0) = 0\\}$[<sup>1</sup>](#fn1 \"Definition of Sobolev space\").\n",
    "Note that $H^1(0,L) \\subset AC[0,L]$[<sup>2</sup>](#fn2 \"Definition of Absolute continuity\"), thus $\\varphi$ is defined pointwise. Sometimes it is convenient to rewrite the previous variational formulation in the form\n",
    "\\begin{align*}\n",
    "\tu \\in V: a(u,\\delta u) = L(\\delta u), \\quad \\forall \\delta u \\in V\\,,\n",
    "\\end{align*}\n",
    "where $a(u,\\delta u) = \\int_0^L A E u'(x) \\delta u'(x)\\, \\text{d} x$ denotes a bilinear form given by the left-hand side of the formulation and $L(\\delta u) = \\int_0^L b(x)\\delta u(x)\\, \\text{d} x + g \\delta u(L)$ is a linear functional that contains the external loading of the problem.  \n",
    "\n",
    "<span id='fn1'><sup>1</sup>See the [definition](https://en.wikipedia.org/wiki/Sobolev_space#One-dimensional_case) of the Sobolev space $H^1$.</span>  \n",
    "<span id=\"fn2\"><sup>2</sup>$AC[0,L]$ denotes the absolutely continuous functions on $[0,L]$, see the [definition](https://en.wikipedia.org/wiki/Absolute_continuity#Definition).</span>  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---\n",
    "#### Geometric derivation*\n",
    "\n",
    "Displacement of the bar cross-section located at $x$ is given by $\\mathbf{u}(x)$. If we cut an element $[x,x+\\Delta x]$ of the bar then the left cross section is displaced by $u(x)=u_x(x)$ in the $x$ direction, and the right one is displaced by $u(x+\\Delta x)$, therefore the strain, or the relative deformation, of the element is \n",
    "\\begin{equation}\n",
    "    \\frac{\\Delta u}{\\Delta x}=\\frac{u(x+\\Delta x)-u(x)}{\\Delta x}\\,.   \n",
    "\\end{equation}\n",
    "Thus, for an infinitesimal element, the deformation at $x$ is given by\n",
    "\\begin{equation}\n",
    "    \\varepsilon(x)=\\lim_{\\Delta x \\rightarrow 0} \\frac{\\Delta u}{\\Delta x}=\\lim_{\\Delta x \\rightarrow 0}\\frac{u(x+\\Delta x)-u(x)}{\\Delta x} = u'(x)\\,.   \n",
    "\\end{equation}\n",
    "Using 1D stress-strain constitutive relation, i.e. Hooke's law, we get\n",
    "\\begin{equation}\n",
    "    \\sigma(x) = E(x) \\varepsilon(x) =E(x) u'(x)\\,,\n",
    "\\end{equation}\n",
    "where $\\sigma$ now represents only the $xx$ compononet of the stress tensor $\\boldsymbol{\\sigma}$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZAAAACaCAQAAAAtWtaVAAAAAmJLR0QA/4ePzL8AAAAJcEhZcwAAASwAAAEsAHOI6VIAAAAHdElNRQfkAhkLDypGpoK5AAAeK0lEQVR42u3dd3zf1X3v8edPe9uyZEkeeOKBbbDBbLDB7BhIAjQkhJGktI8mKUnaJk16097b5NJBRpvdJrcJvSHMhFECxoyYUTOSQIAABmyMGTa2PGRbkrWln/qHv/7y/claFlbA1nn5D/m3zvd8zznv8/l8zvqmHKiMdZ20MnnqNdhuo7dssEG9Jo3aBQL7gZwDNucp23zDm8qNU6lGlXkWKdeuRZN6G222wTaNdmnUFao6MLIEQlqznXZ6Lb6XXPkqValUrcY0YxVr165FnU02ect2DXbapTtUfeBgF0iqx+tOnVrstDZ6nS1HkSpjVKox3nRVcqR1arXJZm/ZbLudtmuVDk0hcLAJZCC6dGmzI36dJdsoFcpVGm+cWcbIkkKDWpu9Zavt6mzXIR1sTOBgF8jeLlnaNtsS7+QpN0q5KhNVm6NEnmxss1WtWlvV2aJRh84gmCCQkUe7zTYnXhcpU6rcOOONM1eRQvnabbPNZptttcUWLdp1BMEEgYw8mjWrTcQ4xUqUGGOiauPNVqpIoUY71Nmq1habbNeqTXsQTBDISKPbLrvAb6Ogv0ihEtXGqTbOoUYpk6Negx222GSzjRq0aAnzMEEgI48ujRptsS56nStfgTLjVakx1nSjlevQrNHOaJxsi0ZNdoV5mCCQkUeHDrtsiwSTkiNPvgrVxqpWbYYKJdq0abLdRrXRPExDmIcJAhmJDlmHDk22eyUSTLZchapUGKvaOLNUydGlXatatTbaZId6O7QEwQSBjDzB7J643B6H/FlylKlUbqxxxpmjQgpddsUO2Q7b7dAeJi6DQEaeYLp02Wpr/E5KrjFGK1dlghpHKJEjC3W22GiF1aHYgkBGsmTa1SaGlSlUHs/DnKMoCCQIJJCkRUv8/7TiUCDvPlmhCN6jpEIRBIEEAkEggUAQSCAQBBIIBIEEAoEgkEAgCCQQCAIJBIJAAoEDQiBZymX38VlKubwh5KDvFBk9hBR7UqDsXS/nMqV9fFK6n5aU5Bn1rt/lKAVD+NVoue9yvvOU7y+BnOCUPhdBdJvjQvn9/LrImXs11mOc1k8uZvnQoAu91AST9vp2kQtM3m9FOdr5Q5Lbnzmhj0+qXaJyP1Tx+WYNU/OZ2W8NJZnsI312BJDrNFU93pvr3GESSIFz1Qy69I7aHwI5zllW6uzz89+q8P5+fn+Wb6jOeGeh86zU0ecvnlLswkGtS8p3oR+4wbSMd3N8XMqq/Vbo5/uJM/f5VzXmRvsOsxzmoz7lgrjq1qp1xZB63qTt/ogKTw9LM8v1Bd8zdVDffV6HS/rxB47xbYdmvDPdZZ7QPCw5P8GPXTKob+6y0kW93+O+CKTCx9yrrp9vdLjNKQ7vs6/8nKmKMvrjT3jAln5S7HK74x05iNy1u8nz8nuktthsd+3D9qOcAUrgDKX+eJ9tyDG2eA0pS11gq9dd5BYnRZ/er8z73lFTONqJ7uin4+qd7EHV/nyHm+HiQaXY7U4zndynhf+swzIsTIE/8VS8x3//5rzA+xS53CGDSvE1v/GJ3tz5fRHIuZo9OcB3aj3lo702s5Sl2uVkmNhzpD0xQIrbPOHSQUQi3dJmeDnewQfFPuJejYO+wypXKezn8yUedqtTnbGP1bnYSl2Y5Rg/8oDlvijf1yMr0u4eF74DNyvXJVYmNmINlsssGITzscQ3rXW56YPsi3/lo31EVafLk85wek4wzgNDsJdXDsKdPNZ6/+5wFw4y1YdM6E3agxdIiaVWDKKXesQRZvTy/gLFbpeTcLEKnWdFP+7VHlaabfagmvdsv82wFkeq9Ot9KPoSx/RjQ8aa627XanflPtmQKcZGHcuJLo1sxUb3O8ox0Teek2XxkAUy2wyPDOF38/eKB/bmKJ3udKOZg7Qh/Ea1hb28P800N2UIJNsH/EbDEARylDEDfKfQye73Uxt8bJA2pNFTPri3e5jVRxYKlSpVqkxZ5B/PVu6FjH6xxgwlGG16osFsUuf4vVIsdob7vaHb2Pi9Gao8l5GXajOUYpTpiTGZLWpjd6Q/ZhjlBeNNjQcKllibOJuXAlNNki3beJN6tUr9HZ5wijW2esKyXmxIqUONk5JnsvE9ivlEr0T7BhuVmxZFVPVyVETfaPWcUwc55lcc10xpJOZFNmWcD5lvsilyZakxpZ/oZuCDIvKc4kGdbrHW5T2iuyyVZhiFEtONiePEndb2IvYc53jcGu0JUY4z21M93KIpJsmRbbzJ/Qz3DJzzYzVYY7WberEhWaqidlZmutHx+0+b3SNC7sPnLnOhc5TL165dtrv8UNrhmhPVkOcMo+WZ4V5TlVjoH6LHELR6xUI/7eH3L7LVGmO1JrQ/T7tNGSMclbId5h5TFTnWP0fng7R72VFyBrReC6UdocJ4Jb5nlwILLEsU5ninaXaYzTYqd4Qd/mUQ9itpP36EVtc615V+lej7Zjleo6M9KUu+xR51faLSj/Pz6P/LnOk13UiZojl+cAO/d4byDDH3Pkp0qWOVyNEirdO3PSTbUV5J3EelM3WYos3Lqs2Q7R8T+xT3jaO0ewGvuMHfu9g1ie5xkUN0OtIyk+VY7DueiRrvC85V2OOaC2V70nQNKqSiOjlUUaIEqHG6FrNts94Y8+zyjSEewFfoJHfo1u0Gl/iY22xItLMlxsoy1zJTFTreNdHW5rXyzLJxIAtS7h992Eo3ShnnFj/2sDRmWp/I7ily/NzPHeWvrZDr9ETPsMn4Hv1WjfnuRaPWhEZnektb/OpkxX7uFnN82SO6nZEwxpuM6zc22H3jx2jwqNvc6Gzno1yF9Qkb9j6/c7tH/I15ljvKCfs0vHiKNZEdeNwySxI2pMap7ne7tb4m268tzvDtZynw++j/zZ6ORDDNme5K9J/bFA0YhczxA7Pc4W5TNLvWtV5EifHeSHRcS612q/t82hJ3OcyiAUuu73HBU6yIjry72VpXJGzIQuPc6hZVrvaUBmeaFH+2WUUPF7TUIvfp0mKX6ti+zlAXnVy5ezh+qWfd7iF/bYHljnTCkOfAjtXgZfCSmxzhogx7XuoXbjHT31mp0xnGR5+02Ln3SFbOXq8/b7JPeR0bfdtTXoqFU5vw1Sf4L52KlHvYW+61KuEs1auQlxi6SznbM2rRoEml/EgW5bbH/XuRSe7Wochoz3jDCq/5XSLF0fIHCLYrzfFTT6NZt5PdpECO+oQPvM1LKFRgpR3+ybYojwVydCOtSI6iyPaldCTku8d+6MWGzLHKJpTo9IQ6fxVVzh5hPZ3IxZ6md5WNvqopfqdJfr8zCFT7J6tcbRdGOcEyrVFahQnLc4g2zyJfqcfU+07GAyDIURCXea5C+XKQku7lGSlHaY8HyDNtSK4ZHtEix2hvelmjP/FYwpsv6SHK02ywBk12KldiJxhjZ+LEyal2WIUChR61wzV2JMpn9zmW3ZH1ylWgQHaU856nihU62W3Re2k3uMTHYxtSaLLl2hUa7UWvediV0cGydGndOyrL2atQPugvvB55k/mqI4GkjE24Q+3uVo/Jyj2O1Rnnb3QrzfDC5zjbU6ahTJExCqOGNzbR83W4x3ZMUuUxrI0fg7MnxYGOl5ipNAqFRxtvNQqNSlT6m9FdLbDBWt0eikP7z6nWjW5lFrpGhxRSXvaD2FF4237stiH3uMDp7ojcoybkmO85m7W7N2Mg+3Df22tU61KHuCpyIPfcYeEAcyF/pManox53u0qlkUDylCfucrN7pXGE7V4g0Wx31+LFTouaTrfjzXBO5EO0+KEXe4h4seWJBnyzy1zhFq8h7X7bUW2yH0lbn7DUpJVk9P0TXWCNTyJfubSySCAVshNNe703o/rZaI1uD/dop1c4URrdUk5U5eIo542+32Og+Fj1ifb4kpt80UW+AzotV4eJavwHXvVqog6y9u6keja7M9THhn+ctLfiHzdkzDnsfsrGYdIZgfvbzb07oeezXa9WCrkuVh43hIbEtzqiYcqZsj3fS4oDP59jgR1Rzz3dGKvQrjkxxVgfNc6FVmU8I6TeLQojgUw02k8ju5LSELuUSfux24b8xLmutEID0cxQmblu28tnPlxzj8N7Us4x3+e9IUdW4vtd/Z7fW2CpX8f+8WRbYqvUaVfiLnfFLtC6RIf2dhN4zKtxSRb7rUejfrgrrune7MduG3K9r/iwa9AV1dZUZVHk0bO20onu4FzLrI2a81nmx4PAjYkBG1H7yrLQS70MWXd5yAtRzrOMca9nZEV2PnPmq8hJbk+0lrQbfDS2IXva2Qy5GQNEe8qns3+B5DrMq5FhS1nihUSx1Wf4lSndUhZZpxbFDrEmLpQ8WxNB40nqol6NHGucGI9O7J0ii7zprcjleiVuMnm2Jdyd3png9ajRnKjFI2jTnjEin9Jtgjm+rxuH6LQJbYmCqrPN73px5RZbnXF+FY+724WRDUmh22xjo65lhq1RH8liT/bI+Ynmuto2zFXgN3HJN/brQo41wX/F05XHuDuWVrsWJT3ustJR7tKOGnlRz7ybNxJWu9bqhBvb0wlc5N4ekr0pYUN219ZJ6qxDnunWxnWeZ0ciRD9Srjvj3K5yQjxVXK9YdqJJpnQbb57/J42J0omAudu62E5k2eLFHuNfb3NMhv2AF93oiy703Yx2tsEGFJpiTXSfWbozRgN7CdI7bImPvZzvaD9MRBLr1cQ91Vxfc5hyx3tJM+aZnegzxqmLG8UEx1gef5a23ei4OtcnArZZvuYIZU60WiMOMy9RPTV2RA5F39Rr1o1qS/3EqijiGRv3Yx/xN3IdocoLyHHqAD5/Mro50TqHmJT4V+kBXKlUgc+4EidIW4tSS+Jup9r0Hk7OfO+zQoFJJmZMS43RGYuq91H6HXFTOl+bW+JPmu2IhzNSzvMVJWaa7HmkLI6Hkvd1/Gq0XRl3PEmHB8z2YUzyj06Ub5HX1GGahYn6r9IUW7JRTrc8YSe3K44HVteriF2xLBf7sjzz1HgB2U4d0pq3Qmd6xYSMfE+wUp2Pm4gZrrFAiZOsUY/ZDk907BU9x7D2drHudJWZXneoP/UjKxKfPOdsFZFzcprLLTdVrWxFDnG4OxMm9VCrtKPCkT4p5eZItaVmmKDAqXZZL+05F6iMNLvYFVaothUlahzlrkTPMtOLAwrkQbPVaPWXnvYdXaj3SjzjWuiPpd1sjnXKFTnZroT32T9LfcwFexnjXHmWONOTPmWF6aptUq7BWV6OXbhjbc7wzmf6nsN8RBay7LA8Mdb1Vr9Lbna6w3wV2pxskb9LWIVWL5ohSxq5LjPFzxxpnTIFjpU1pHVoeT7hQh/usQauW7GUK9xsgT/xeym7dCpT5uR4tGv3vaxVjzJH+JhD3R3lrtB00+RY5HVv6vCSXBMit7jQJ+S40TzrjFboJM0ZMdpgOcEnXNrDHe+WZZRKF/tXJ/u4R1TYplupKkcnJgLGy4+HpPoUyKNyna1Rtv/0TIaJfV6XaVHVL1fkUOv8uQtcocH9Ca+x3EQ3ROPcJ3lVlgXW68REp1ltjSrH2KTdS1rNiATygHKTve4qF7hcg18lXJpRpvjlgEXzpBucJ8uz7onjpUd8SIFWNPmOORb7pcecpNom9/Xq86d69ahv6nWyqltKylb/bIzjfN9KZ9pubbx4JsdJHs3wasd61rPRNVJej0stZUG/i0DhJz7oEm3aXR0NN7y90uBzUeDb7keOdbKVnnCGS21zTz8zCal+Ip7X+rjntC6jPO5ryjX5K+e7VIOVifmMAnPcFo2oLVanwULrtKDS6ep8U7aT3abDmzY4PBJIs+863GLLPGaRsWr7qJ+Bcp7rtl4/TUtpkGOFfzHRGz7rgy7TaEXCZhzuzYzurJ/iKe91zOgrvpDoY0qjsLvnav6z/GeGT9wfX/a3iVsrk4WcvVI81XW97HWocl2PlaEppT2uPN6tiTmJ4mjwsUhZH0U82b8NYXdGSqm86G/y11Ndu/fMbK9M8Yu9lud8OlHab1+nt5yPdmNiHr4o8vELjRpgIdFXLTFUcpTJRrbRPWa8F7hlkIvM+WPfTrS0getnjzv29XiRzr6TbGd5GWOL33JFbxfbm9aEv5vkFvPidS3tGnWjw86M2egCZ7ojMf3TPz83I5586tAgjc4eKeY52y/3mknofQFCY48rb3Sf8+N7bIpCx+aM8bMktf5lQFeu9+u2R3+bMsLxV/cO+nrlHE9mDGv3fZ3ecr7THZbGld0cxY0t6gdYxfzjd7BAvlODLnTZmTEIkWWpB3oMaPTN3UrNT8wFDVQ/e6zBDzJmmvaNZDtLWtcjlCRc3n4F0hcvedAH+lntD2fb6b5Bp7jWfd4/wHz26VotG3Jx3KpoH3qbNq/ut4el5Tvefw/qm0c41M/e0SN07tE8hMWO6wfV7ewbu+epB8sWtzh3CDb7jX1YoT04Sl3o5t5WRO/bhqlf6HJGP+ZvrtmuHXBANsntmp3dT4qzzfeTIa8lYodrLdmPOwoHzyw58RKT/qhytht6mbPYF5r82NGDWvE8vEx2smv3SXYrvOID7/ozBnKd78V46vgdKu3IPlfIpMw2bp9TLHFkP+s2Z5nQZ8P66SD3KEzpEav8YZg0SMs1sdfNAXzK5/fpehPeAwKZZsoQLO2Rg45Zh4sS8/vdKj4MZJnnwgHcsXfGPMv62R+/v8k2y/EZM79DocBSMwZltYtc7ZsZuy/3v4A/NOCuisCwUOg413jdvcMokMV+rc3Lg97G807lcZVX1Vnh6HfYU/3SGtc4ZoAVtmP8mx12+o8hTvMNhuPVWeFPTRnWbuyAZ3/2wFlGmep4ZzpBlZTHfXpYHkbZrdA/Ox085/2JxRPDxRF+GcUxd/i/OodYat2KfMNi3bZ4wq884bU+Rpo+4Yfy0Omrbh+WBtzlSD+Rr9MaKzzo9zYP09EJBzj7MzzKc5pLnBr3elP85TA93rgg3tRZbcwfQCAT4pH9E3x+nw9HSJb2NKRU+6DFHnaTu3sdVp4SxXk5Lo53IO7vTqZKDnLMMdOZlrmul8V7gf0qkFZ3esQcpzvPXPle9oV30Jj6D+y+6TKwZq9VqMPBy1Y7At1u8g/vwCoWu9ZELVZZ7kGr+phv4nfqjcIu33THMFmQ450uW51HLbPSW4OeuwrsB7dtnMvd7/ZhHL6b4ad2WD6over7447e534v+P4gDwDoOwa52XIXqx4gUC/0GWut8/lhfJDncV71b07a5ysUDnmHYiCD0RYM68m/E9xk7h/wfor3wxGZueYmDgnon7/w5WG9n0pHDel+lrrOkpHzBN7h6+N3enbQscvb4+Dtkakvkq9RJ3KUxELr0iiNPEWytWtDaVzN3XbpQJYynXE6b+/Ta478/VLZ0XKDXCWxh99hl267t3zuvnJWtGpnz5Wbol+0ReFsct/cLu3IUiodp1OUcD6bI9dwfZSvnGgtG3RGy3byFWmK0imSI9so2dIadSFbqY5oKUtxYtR+912llEhF6QyObRnbxval47vcuR51m0e9OUwu9IgQyOA51d9HjTzbr3wpGse5yOc8j5n+VYVuZFnts+pwqq+6NWrE1zg+Wh7S4n95HDX+3XO+qhN/6rJ4J/O3XI8cXzfep2zEMb4WNbUsv/YlTbjYn/mMp1HlWmOjM0jWuyr6xdf93HeR58vOjq7c5avuxSjft9UXteMSfx7vbv9Z9Iu/tdCfewWzfcvo6K5e8meacYa/9RX3o8q3nITzpLR5SgNKHW2b56V0m2dSvHH2JW9IybFQylPD3mC7zMcY77fUG35tuYf23kMRBLK/C701amjZ8TLFzvgQgW5tWqOm1BY1iy6t8WBre/TrlLbEL5Lp7BFIZ/Rpe0Y6e7Zxtidysyed1ujzVI8r77F2LdE301H+u7XH6bx9B6k4Nx0ZedxzV+2JtLsSd5DSKqVNq1bkatOmVSqRrz3ppORE3+0Y5rpKJ+6l0Y5ejno46ALqd5/cOOxL6YhPGsnVHLkWRXEu05p0I1eB0f7B1dYqjs4k2e1udEYOSle0fqsgdoNSWqNVYsWyNEkjJ+EGdUanY+QpiNMpTly5Of5Fe+SqFSZ8+BYdSCmWjt2p/ITjuOcXOf3eVWGUTpYCn1ToO7J0R3nNUqwzuqtCefE9t0aSLCbaUzm8fNQNaj3oVk/YcrDL471hQTp66fda4/mBrl5WbnboUBg1sKZeerldvaTzNk0JUex98OXbNiDdy5WTv2jZaxFld+LKbb0s2mxJWM3GfsohrVmLVMbQazI3Lb0s32z6A9VWnW/72SB2eAaBBIbZtqfek/l62IPD7sgFgQQOWNpG1u2GZxQGAkEggUAQSCAQBBIIhCA9cCBQbr5K6z1jrkl2+s3BOPQbLEhgaIzxAR3W+pwvmabUd/t80HUQSGDEkXKW1R7zgnyL3a9dbcbjU4OLFRjRFKjzrN3Pt1qu0Z3uH4ZTtoJAAgcorR7UhSkq/VbvC2uCixUYsXRH644P1+ZlpOS+R5fGBIEE3gWyTVIk5RSv2YQxzjw4jw8KAgkMhdP80hLjLPCGNinHajs4dxeGGCQwtCC9Q6lTXW++E1QqTjxCKQgkMOK532aVVnvdEaba4Pl3cMB4EEjgoKMtfrr4M70+5zbEIIHASCAIJBAIAgkEgkACgSCQQCAIJBAIAgkEgkACgSCQQCAIJBAIAgkEAkEggUAQSCAQBBIIBIEEAkEggUAQSCAQBBIIBIEEAkEggUAQSOA9SXcogvcC4VST9w7ZSlUYpcIE453i/lAkQSAj2XZnKVRpjAo1JqhSI1+WLp022+IGK0MhBYGMJPuQLU+FKmNUG6fGWGXSOrXbbrPnLVen3g47o4OhA0EgB7UgcuUpU2OsKtVqVKiQ1qFJo1ovWm+zBg3qtYfiCgI5+F2mfPmKVBunSrUqY4yRq8UuDbZ42VveUq9J48F6TGcQSOBtUgoUKjLGBNXGqTRGqVJNdtmuzss222CbFs1apEOBBYEc7BQoUWSU8carUWW0EsU61auzzUu22KRWs1atIYoIAjn4yVWmVJkqh6g2Xpli+bLtsF2tF2yx1UYN2rUfnM/GCASBZEYRo4xWptIE49WoUCBXtl22qfWsLbbZYpt2HcFCBIEc/DFEliIVyo0xzgTVauTJ0a3DFps9p1adOls1SwdBBIEc/PZh9zxEpTGqjVejWgm6dKqz2YseUGenOvW6wkKPwMEukCy5co1SrVKVGjUqVUjr1KpBrdVW2KJevR06QkMIHHwC6e4hiDx5ilWrVqVatQoVsrVr0miztR62Ub1dGsI8RGAkCCTLKDXGmmCscSpVGqVYsyY71XnVozbYrtmuMA8RGHqoeqBS6cdyZWnXaLs6tbZ4y2bNWrSEoDow0gWSbyHe0qjtYH1GdyAQCAQCgUAgEAgEAoG3CWvVAoE+OcTVpoRiCAR656+0+1IohpFMOBerbyabodblpoaiGLlkhyLokz/ylI0+ZLPHQmEEAplM9b+VmGGNVcGGBAsS6MlFVnvRdmNcpLaHDRntSIdLazDPfAXqwv6RIJCRZj9O9gtt2GSpeZbbGX821vs1yvFZVQpM9Tcety0UWQjS31sUOFHpsKV+mqc0gDWuN9uHEl3K+/zeIx43xSL3KNQ5LBuuZjosNM/A0KlynUOHKe0p/o9R8asZXvFCHIeUO08uJlrlkxitZljWRH/aF0IlBwvy3uR0T6mPX73ieofFNqTBfTowS57fYafaEIEEgYwkppjQIyi/0TpXRDakK3KoFqizLpRhEMjIY4nfJuzHbhtyXRyH5KiUI88iq+zAOCcdwNvOAgMQFuP1ZLLT3ecDGe+lNWlymV94zfn+0mc0m+M/paUcZ2twsYJARg4fcJELe3k/z2Eu8K8mKlWt2p3GmeVQeZ4MhRYEMnL4bxf3ahG6Zdkg5WdeVeoJt1nkaK97MDzdIwhkJPGsZ/v9fKd7pHTjvuhvIATpgR7WJPNvIAgkEAgCOXBJmWD2AT3cmmOektAgg0D2P9mm+ay7XHJAuzy5vuH/W6o8NMoQpO+/3Fdb4BznOFSOh4Zx8eLwRzWFKpzjLE+6y8NeCU0zCOSdNqkiV1roOEXRO6cnFhgeiDUxGaVOs0St++QMMJoWCAIZIOposcwDljjd4QrwjH8/gGuiwEJV2OoJD3nMaaFxBoG8M9I2WOsBlU5xiTPV+t0BfDeFGqxzq//yjFYcE1vGQBDIEG3Ibldrq1vd4xTdB/TEXZfves7q8CSTIJDhoNnyA/wO2t0SmuN7jzBRGAgEgQQCQSCBQBBIIBAEEggEgQQCQSCBQBBIIBAEEgiMQA7cmfQ2z2g6iGvmdfmheb77/A+jlumzuUuY8AAAACV0RVh0ZGF0ZTpjcmVhdGUAMjAyMC0wMi0yNVQxMToxNTo0MiswMTowMBoaeRMAAAAldEVYdGRhdGU6bW9kaWZ5ADIwMjAtMDItMjVUMTE6MTU6NDIrMDE6MDBrR8GvAAAAFHRFWHRwZGY6VmVyc2lvbgBQREYtMS41IAVcCzkAAAAASUVORK5CYII=\n"
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "%%tikz -l arrows -s 400,200\n",
    "\\input{fig/equilibrium.tikz};"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "From static equilibrium of an infinitesimal element, see the previous figure, we get the following force balance\n",
    "\\begin{equation}\n",
    "    -\\sigma(x)A(x) + f(x+\\mathrm{d}x/2) + \\sigma(x+\\mathrm{d}x)A(x+\\mathrm{d}x)=0\\,,\n",
    "\\end{equation}\n",
    "which can be recast into\n",
    "\\begin{equation}\n",
    "    \\frac{\\mathrm{d}}{\\mathrm{d}x}\\left(\\sigma(x)A(x)\\right)=-f(x)\\,.\n",
    "\\end{equation}\n",
    "Using the stress relation we get\n",
    "\\begin{equation}\n",
    "    \\frac{\\mathrm{d}}{\\mathrm{d}x}\\left(E(x)A(x)u'(x)\\right)=-f(x)\\,.\n",
    "\\end{equation}\n",
    "If we suppose that $E(x)=\\mathrm{const.}$ and $A(x)=\\mathrm{const.}$, we get\n",
    "\\begin{equation}\n",
    "    -EAu''(x)=f(x)\\,,\n",
    "\\end{equation}\n",
    "which is the same relation we got in the previous derivation."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Implementation\n",
    "\n",
    "---\n",
    "First of all, we must import ``fenics`` library to ``python`` interface. Moreover, we import ``pyplot`` library that is used for graph rendering. And we will also need some features of the ``numpy`` package."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "import fenics as fe\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Problem parameters\n",
    "Let us set the material parameters of the problem. Following is the table that contains Young modulus and Poisson's ratio for some common materials.\n",
    "\n",
    "| Material | Young modulus $E\\ [\\mathrm{GPa}]$ | Poisson's ratio $\\nu\\ [-]$ |\n",
    "| -------- |:-------------:| -----:|\n",
    "| Rubber   | $0.01 - 0.1$  | $0.4999$ |\n",
    "| **Cork**     | $0.02 - 0.03$    |   $0$ |\n",
    "| Concrete | $30$      |    $0.1 - 0.2$ |\n",
    "| Steel    | $200$      |    $0.27 - 0.30$ |\n",
    "\n",
    "We decided to model a cork stopper, since cork material parameters justify the assumption $\\nu=0$ made in the governing equations. Thus, we choose appropriate dimensions of the cork stopper (radius 1&thinsp;cm, length 5&thinsp;cm) and some suitable magnitude of loading $b$ and $g$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# --------------------\n",
    "# Parameters\n",
    "# --------------------\n",
    "E = 0.025e9  # Young's modulus [Pa]\n",
    "A = np.pi*1e-4  # Cross-section area of bar [m^-2]\n",
    "L = 0.05  # Length of bar [m]\n",
    "n = 20  # Number of elements\n",
    "b = 0.03  # Load intensity [Pa*m^-3]\n",
    "g = 0.0005 # External force [Pa*m^-2]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Geometry of the problem\n",
    "\n",
    "Next, we define the geometry of the problem. There is always a possibility to use an external mesh defined in *.xml* file. In our case we use the predefined mesh provided directly by FEniCS. Initialization of the FEniCS object ``IntervalMesh(n, x_1, x_2)`` creates an equidistantly discretized one-dimensional mesh with ``n`` elements 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": 6,
   "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": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "# --------------------\n",
    "# Geometry\n",
    "# --------------------\n",
    "mesh = fe.IntervalMesh(n, 0.0, L)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can visually check the geometry of the mesh using the FEniCS ``plot()`` command. The plot object is then displayed by adding the pyplot command ``show()``."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAIQAAAD8CAYAAAC2NQwLAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4xLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvAOZPmwAABdhJREFUeJzt2L9P3HUcx/HXu4AtcZBEcChVr11qig4kRGc1ijpYYh2c/AN0baNNHWwXowwuDsbNrfVnVyZnzRE0jYkYimkEljamJjVYsX078G251xXKF77fg7N9PpJLj899vse78Mzd54jMFHDLnt0eAN2FIGAIAoYgYAgChiBgCAKGIGAIAqZ3OxcNDg5mo9GoeRR00vT09JXMHNps37aCaDQaajab27kUuyQiLpXZx1sGDEHAEAQMQcAQBAxBwBAEDEHAEAQMQcAQBAxBwBAEDEHAEAQMQcAQBAxBwBAEDEHAEAQMQcAQBAxBwBAEDEHAEAQMQcAQBAxBwBAEDEHAEAQMQcAQBAxBwBAEDEHAEAQMQcAQBAxBwBAEDEHAEAQMQcAQBAxBwBAEDEHAEAQMQcAQBAxBwBAEDEHAEAQMQcAQBAxBwBAEDEHAEAQMQcAQBAxBwBAEDEHAEAQMQcAQBAxBwBAEDEHAEAQMQcAQBAxBwBAEDEHAEAQMQcAQBAxBwBAEDEHAEAQMQcAQBAxBwBAEDEHAEAQMQcAQBAxBwBAEDEHAEAQMQcAQBAxBwBAEDEHA9FZ9gvMzi5qcmtXS1WXtH+jXifHDknTfrD37xJC+++Xyrs+x0drE6PCWfp+RmVu6QJLGxsay2Wzq/MyiTn5zQcsrN24/1rcnpJBWbuR9sdaum2br7+vRB689pYnRYUXEdGaObTh4odIrxOTUrMUgSSs37/xh3ctr3TLHemvLKzc0OTW7pVeJSmeIpavLVS7HDtjq76hSEPsH+qtcjh2w1d9RpSBOjB9Wf1+PrfXtCfX1xH2z1q6bZuvv67l92Cyr0hni1ntTt5yod2ONTxla+5SB/4+ynzL4wxQMQcAQBAxBwBAEDEHAEAQMQcAQBAxBwBAEDEHAEAQMQcAQBAxBwBAEDEHAEAQMQcAQBAxBwBAEDEHAEAQMQcAQBAxBwBAEDEHAEAQMQcAQBAxBwBAEDEHAEAQMQcAQBAxBwBAEDEHAEAQMQcAQBAxBwBAEDEHAEAQMQcAQBAxBwBAEDEHAEAQMQcAQBAxBwBAEDEHAEAQMQcAQBAxBwBAEDEHAEAQMQcAQBAxBwBAEDEHAEAQMQcAQBAxBwBAEDEHAEAQMQcAQBAxBwBAEDEHAEAQMQcAQBAxBwBAEDEHAEAQMQcAQBAxBwBAEDEHAEAQMQcAQBAxBwBAEDEHAEAQMQcBEZm79oojLki61LQ9KulLHUB3SzfPtxGyPZ+bQZpu2FcS6TxTRzMyxWp6sA7p5vm6ajbcMGIKAqTOIz2p8rk7o5vm6ZrbazhC4N/CWAVMqiIh4KSJmI2IuIt5d5/G9EXGuePz7iGi0PHayWJ+NiPH6Rq82W0S8EBHTEXGh+Pe5umerMl/L449FxLWION6J+e6QmXe9SeqRdFHSIUkPSPpJ0pG2PW9J+rS4/4akc8X9I8X+vZIOFs/Ts9n3LHurONuopP3F/SclLdY1Vx3ztTz+taQvJR2ve771bmVeIZ6WNJeZ85n5j6Szko627Tkq6fPi/leSno+IKNbPZub1zPxN0lzxfHXZ9myZOZOZS8X6z5L2RcTeGmerNJ8kRcSEpPlivh1RJohhSb+3fL1QrK27JzP/lfSnpIdLXltFldlaHZM0k5nXa5yt0nwR8aCkdySdrnmmu+otsSfWWWv/aLLRnjLXVlFlttUHI0YkfSjpxRrnKvW9N9lzWtLHmXmteMHYEWWCWJD0aMvXByQtbbBnISJ6JT0k6Y+S11ZRZTZFxAFJ30p6MzMv1jhXHfM9I+n1iPhI0oCkmxHxd2Z+0oE515Q4GPVq9X3soNYORiNte96WH4y+KO6PyA+V86r3UFlltoFi/7FOHdCqzNe2533t0KGy7H/sFUm/avXEfKpYOyPp1eL+Pq2ehOck/SDpUMu1p4rrZiW93IEf+rZmk/SepL8k/dhye6Rb5tutIPhLJQx/qYQhCBiCgCEIGIKAIQgYgoAhCJj/AM6TxpxOi8rkAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f96a3664080>"
      ]
     },
     "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": [
    "### Function spaces definition\n",
    "Now, we define the function space that contains the solution and test functions by initializing the object ``FunctionSpace(mesh,type,deg)`` with three arguments.  \n",
    "The first argument is the discretized domain contained in the ``mesh`` variable. 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 later in the course. In agreement with the [definition](#function_space) of the function spaces given in the first section, we define the conforming space of finite elements, that is the *Lagrange/Continous Galerkin* element (hence the acronym ``CG``) of order one, in other word, piecewise linear functions.  \n",
    "  \n",
    "Then we initialize the test and trial functions specified by the function space as ``TrialFunction`` and ``TestFunction`` objects, respectively."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# --------------------\n",
    "# Function spaces\n",
    "# --------------------\n",
    "V = fe.FunctionSpace(mesh, \"CG\", 1)\n",
    "u_tr = fe.TrialFunction(V)\n",
    "u_test = fe.TestFunction(V)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Boundary conditions\n",
    "As a next step, we specify the geometric boundary of the domain through which we apply the Dirichlet boundary condition. It is necessary to mark those entities of the mesh that constitute the boundaries $\\Gamma_{\\text{D}}$ and $\\Gamma_{\\text{N}}$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "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": 11,
   "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": 12,
   "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",
    "#### Dirichlet boundary conditions\n",
    "Dirichlet boundary conditions are sometimes called essential, since they determine the definition of the function spaces. 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": 13,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# --------------------\n",
    "# Boundary conditions\n",
    "# --------------------\n",
    "bc = fe.DirichletBC(V, 0.0, left)\n",
    "# Equivalently, using boundary variable:\n",
    "#bc = fe.DirichletBC(V, 0.0, boundary, 1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "raw_mimetype": "text/restructuredtext"
   },
   "source": [
    "#### Neumann boundary conditions\n",
    " In the [formulation](#governing_eqns) of the problem, we also prescribed a normal force $\\mathbf{g}$ at the right boundary of the domain $\\Gamma_{\\text{N}}$. Such boundary interactions falls into the category of Neumann boundary condition, which 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 appears as a term from the integration by parts formula. To include this term in our implementation of the problem, we first need to define the integrating measures that correspond to the manifolds representing the geometry of the mesh. To this end, we initialize the ``Measure`` object with input arguments that specify the integral type and the respective geometry to which the measure corresponds."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "dx = fe.Measure(\"dx\", mesh)\n",
    "ds = fe.Measure(\"ds\", subdomain_data = boundary)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Weak formulation\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."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# --------------------\n",
    "# Weak form\n",
    "# --------------------\n",
    "#a = E*A*fe.inner(fe.grad(u_tr),fe.grad(u_test))*dx\n",
    "# One can equivalently write the bilinear form in this one-dimensional problem as\n",
    "a = E*A*u_tr.dx(0)*u_test.dx(0)*fe.dx\n",
    "l = b*u_test*dx + g*u_test*ds(2)"
   ]
  },
  {
   "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$. Then we can call ``solve(a==l, u, bc)``."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "# --------------------\n",
    "# Solver\n",
    "# --------------------\n",
    "u = fe.Function(V)\n",
    "fe.solve(a == l, u, bc)\n",
    "# Equivalent implementation:\n",
    "# F = a-l\n",
    "# problem = fe.LinearVariationalProblem(fe.lhs(F),fe.rhs(F),u,bc)\n",
    "# solver = fe.LinearVariationalSolver(problem)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Post processing\n",
    "\n",
    "---\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.\n",
    "\n",
    "For many trivial problems, ours included, there is an exact solution. We can use this solution to benchmark our numerical results. The following cell introduces the exact solution and plots it next to the approximative result given by the finite element method."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXwAAAERCAYAAAB4jRxOAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4xLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvAOZPmwAAIABJREFUeJzt3XlcVXX+x/HXV0BR3BWXRATLUlFEJHfJsjTTrMxyL7OkZmqyaVqc6ldTMzUzNZXVaA7uZpqpOTVlpaaomOKW+66ggKG4AcrO/fz+uESaqCj3cLfP8/HgEffec873cyTffP2ec75fIyIopZTyfJWcXYBSSqmKoYGvlFJeQgNfKaW8hAa+Ukp5CQ18pZTyEhr4SinlJVwu8I0x04wxx40xOxx0vH8aY3YUfw12xDGVUsoduVzgAzOAOx1xIGNMPyASiAA6Ac8bY2o64thKKeVuXC7wRWQVcOr894wx1xtjvjPGbDLGrDbGtCzj4VoDK0WkUETOAVtx0C8TpZRyNy4X+JcQC/xBRDoAzwETy7jfVqCvMaaaMaY+cCvQ1KIalVLKpfk6u4ArMcZUB7oC840xv7xdpfizgcAbpeyWKiJ9RGSJMeZm4EcgHVgLFFpftVJKuR7jinPpGGNCgK9FpE3xmPteEWnsgOPOAWaLyOLyHksppdyNyw/piEgmkGiMeQDA2LUry77GGB9jTL3i78OBcGCJZcUqpZQLc7kevjFmLtATqA8cA14DlgMfA40BP+AzESltKOe3x/IHNhe/zASeEJEtFpStlFIuz+UCXymllDVcfkhHKaWUY7jUXTr169eXkJAQZ5ehlFJuY9OmTSdEJLAs27pU4IeEhLBx40Znl6GUUm7DGHO4rNvqkI5SSnkJDXyllPISGvhKKeUlXGoMvzQFBQWkpKSQm5vr7FLUJfj7+xMUFISfn5+zS1FKXYalgW+M+SPwGCDAduAREbmq5E5JSaFGjRqEhIRw3lw6ykWICCdPniQlJYXQ0FBnl6OUugzLhnSMMU2Ap4EoEWkD+ABDrvY4ubm51KtXT8PeRRljqFevnv4LTKmrFT8eEldd+F7iKvv7FrF6DN8XqGqM8QWqAUev5SAa9q5Nfz5KXYMmkTB/1K+hn7jK/rpJpGVNWjakIyKpxph/AUeAHGCJiFw0cZkxJgaIAQgODraqHKWUci2h0WT0n4z/3IfY3ngQUelfwAMzIDTasiatHNKpA9wDhALXAQHGmBG/3U5EYkUkSkSiAgPL9LBYhfPx8SEiIqLkKykpibi4OGrVqnXB+8uWLQPsPd6RI0eW7F9YWEhgYCD9+/cvdy3Vq1e/4jbjx48nOzu75PVdd93FmTNnyt22Uqp8CopsrE88xTvf7+Huj+KJ+CSXSdk9iTo8GVuH0ZaGPVh70fZ2IFFE0gGMMV9gX8hktoVtWqJq1aps2XLhJJtJSUn06NGDr7/++qLtAwIC2LFjBzk5OVStWpWlS5fSpEmTiiqX8ePHM2LECKpVqwbA4sU6/b9SzpJ8KptV+9NZtS+dHw+cJCuvEJ9Khsjg2rwblcmA/Sux3fw8lTZNg+bR7tnDxz6U07l4eUED9AJ2W9ieS+nbty/ffPMNAHPnzmXo0KGlbrdz5046duxIREQE4eHh7N+/H4D33nuPNm3a0KZNG8aPv/giTlxc3AX/YnjqqaeYMWMGH374IUePHuXWW2/l1ltvBexTVpw4ceKSx01KSqJVq1aMGTOGsLAwevfuTU5OjuP+MJTyIjn5RazYe5zX/7eT296No8fbK3h50Q52pGbSv11jJo2IZPP/3cH8PoUMPPgyvoNnUqnXK/bhnPPH9C1g5Rh+gjFmAfb56AuBn7CvTXvNXv/fTnYdzXREeSVaX1eT1+4Ou+w2OTk5REREABAaGsqiRYsAWL16dcn7AAsXLuT6668HYMiQIbzxxhv079+fbdu2MXr0aFavXn3RsSdNmsTYsWMZPnw4+fn5FBUVsWnTJqZPn05CQgIiQqdOnbjlllto3779Fc/n6aef5r333mPFihXUr1//gs8uddw6deqwf/9+5s6dy+TJk3nwwQdZuHAhI0ZcNAKnlCrFgeNnidt7nJX70klIPEV+oY0qvpXo1Lwewzs145Yb63N9YPULb3BI3XzhmH1otP116mbLevmW3ocvIq9hX8DErZU2pANcckgHIDw8nKSkJObOnctdd911yWN36dKFN998k5SUFAYOHEiLFi2Ij4/nvvvuIyAgAICBAweyevXqMgX+5VzquAMGDCA0NLTkl1eHDh1ISkoqV1tKeTIRYXtqBt/tSOP7nWkcTD8HQIsG1RnZuRm33BhIx9C6+Pv5XPog3Z+5+L1Qa4d0XP5J2/NdqSfuagYMGMBzzz1HXFwcJ0+eLHWbYcOG0alTJ7755hv69OnDlClTKMuiNL6+vthstpLXZbkP/nLHrVKlSsn3Pj4+OqSj1G8UFtnYkHSa73emsWRnGkczcvGpZOjcvC4Pdw2hV6uGNKld1dllXpZbBb67GT16NLVq1aJt27bExcWVus2hQ4do3rw5Tz/9NIcOHWLbtm1ER0czatQoxo0bh4iwaNEiPvnkkwv2a9asGbt27SIvL4/c3Fx++OEHunfvDkCNGjXIysq6aEinLMdVSv0qt6CINQdO8P3ONJbtPs6pc/lU8a1EjxaBPNv7Jnq1bECdgMrOLrPMNPDL4bdj+K+88gqDBg0qeR0UFMTYsWMve4x58+Yxe/Zs/Pz8aNSoEa+++ip169Zl1KhRdOzYEYDHHnvsouGcpk2b8uCDDxIeHk6LFi0u+DwmJoa+ffvSuHFjVqxYUfJ+ZGRkqcfV4RvlVeLH2x9uOn/oJHGVfey8+zOczStkxZ7jfLczjbg9xzmXX0SNKr7c1qoBfcIaccuNgQRUcc/odKk1baOiouS3C6Ds3r2bVq1aOakiVVb6c1Ju45cnWn+5YJq4Ctvno1jZ7h0+SQsmfv8J8ots1K9emTtaN6RPWCO6XF+PKr6XGY93ImPMJhGJKsu27vlrSimlrlXx3TAyfxSHmg2mwd5P+X3B06yO86dJ7SxGdmlGn7BGdGhWB59KnjVtiAa+Uspr2GzChqRT/HdLHZpl9+SJ3ROZbAZx/c19+VP7JrQLquXRc0Np4CulPN7etCz+uyWVr7YcJfVMDj0r7+Zlv2UktXySxxI/w7R9GJq6112A10IDXynlkX7OyOGrLUf575aj7P45E59Khh4t6vN2hzN03TwR8+BsqodGQ+KdF47pezANfKWUx8jMLeC77Wks+imVdYknEYGIprV5fUAY/cIbU796FYj/ER6cUaFPuLoKDXyllFsTERIST/HZ+iMs3pFGfqGN0PoBjO3VgnsjmhBSP+DCHZzwhKur0MB3sri4OCpXrkzXrl3LdZykpCT69+/Pjh07LrvdW2+9xUsvvVTyumvXrvz444/lalspZzh5No+Fm1P4bEMyh9LPUcPflyE3N+X+yCDCPfzi67XyrMC/wgMVriguLo7q1auXO/DL6reBr2Gv3InNJqw9dJI564+wZGcaBUVCVLM6PPnADdzVtjFVK7vmvfKuwuolDiuWRUuGzZ49u2QK48cff5yioiIOHz5MixYtOHHiBDabjR49erBkiX1Br3vvvZcOHToQFhZGbOyvE4R+9913REZG0q5dO3r16kVSUhKTJk3i/fffJyIi4qLZNFeuXFmyuEr79u3JyspCRHj++edp06YNbdu2Zd68eRfVO2PGDJ566qmS1/379ycuLo5x48aVzPw5fPhw4NcFVS513Li4OHr27MmgQYNo2bIlw4cPL9NcP0o50vGsXCasOEDPf8UxfEoCaw6cYGTnEJb+MZoFv+vK/R2CNOzLQkRc5qtDhw7yW7t27brovcs6tFLkn6EiP/zN/t9DK69u/1La79+/v+Tn54uIyO9+9zuZOXOmiIhMnjxZ7r//fnn77bclJiamZJ+TJ0+KiEh2draEhYXJiRMn5Pjx4xIUFCSHDh26YJvXXntN3nnnnVLb7t+/v8THx4uISFZWlhQUFMiCBQvk9ttvl8LCQklLS5OmTZvK0aNHJTExUcLCwkREZPr06fLkk0+WHKdfv36yYsUKEREJCAi4oI1fXl/quCtWrJCaNWtKcnKyFBUVSefOnWX16tWl/jkp5UiFRTZZvueYxMzaINf/+Rtp9uLXMvg/P8p/f0qRnPxCZ5fnMoCNUsaM9awhHbAP50Q9CqvehugXyn0h5ocffmDTpk3cfPPNgH1u/AYNGgD2uWjmz5/PpEmTLpg++cMPPyyZMz85OZn9+/eTnp5OdHQ0oaGhANStW/eKbXfr1o1nn32W4cOHM3DgQIKCgoiPj2fo0KH4+PjQsGFDbrnlFjZs2EB4eHi5zvNSx61ZsyYdO3YkKCgIoGSJx18malPK0X7OyOHzDSl8vjGZ1DM51A2ozKPdQxl8c1OaB155iU91aZ4X+ImrYONUe9hvnAqhPcoV+iLCww8/zN///veLPsvOziYlJQWAs2fPUqNGDeLi4li2bBlr166lWrVq9OzZk9zcXETkqi8ijRs3jn79+rF48WI6d+7MsmXLXGLq5MLCwiseT6mrISKsTzzFtDWJLN11DJtA9xvq89JdrbijdUMq+3rW6LOzWLmI+U3GmC3nfWUaY6y9cnr+pEi3veyQJcN69erFggULOH78OACnTp3i8OHDALz44osMHz6cN954gzFjxgCQkZFBnTp1qFatGnv27GHdunWAfaGTlStXkpiYWHIc+HUq49IcPHiQtm3b8uKLLxIVFcWePXuIjo5m3rx5FBUVkZ6ezqpVq0pmv/xFSEgIW7ZswWazkZyczPr160s+8/Pzo6Cg4KK2ynJcpcotfvwFfx9zC4qI+24h099+hsGx60hIPEVM9PWsfL4nsx/rRL/wxhr2DmTlEod7gQgAY4wPkAossqo9wJIlw1q3bs3f/vY3evfujc1mw8/PjwkTJpCUlMSGDRtYs2YNPj4+LFy4kOnTpzNs2DAmTZpEeHg4N910E507dwYgMDCQ2NhYBg4ciM1mo0GDBixdupS7776bQYMG8eWXX/LRRx/Ro0ePkrbHjx/PihUr8PHxoXXr1vTt25fKlSuzdu1a2rVrhzGGt99+m0aNGl0wxXG3bt0IDQ2lbdu2tGnThsjIXy9ax8TEEB4eTmRkJJ9++mnJ+/fdd1+px92zZ881/bkpVariGytO94tl+tFg9q79hreK3uWb6i/y94FtuTeiiV58tVCFTI9sjOkNvCYi3S63nU6P7L7056TKYntKBquXLGTI4Vf5pOh2RldeTtKtE2jTvb/eN3+NXHF65CHA3NI+MMbEADEAwcHBFVSOUqqiFBbZWLLrGNPXJLIh6TQBla8j7LpBjE2bDt1foG2Pu51dotewPPCNMZWBAcCfS/tcRGKBWLD38K2uRylVMTKyC/hswxFmrT1M6pkcmtatyiv9WjG0QRIBX37lsBsrVNlVRA+/L7BZRI5d6wGu5Q4XVXEqYlhQuY9D6WeZtiaRhZtSySkoonPzurx2d2t6tWqIz+HVMP+xX6+1hfbwmpkqXUFFBP5QLjGcUxb+/v6cPHmSevXqaei7IBHh5MmT+Pv7O7sU5WQ7UjOYGHeAb3ek4edTiXsjrmNU11BaX1fz140suLFClZ2lF22NMdWAZKC5iGRcafvSLtoWFBSQkpJSpnvJlXP4+/sTFBSEn5+fs0tRFUyKZ6qcGHeQVfvSqeHvy0NdmvFIt1D7VMTKci5z0VZEsoF65TmGn59fydOpSinXICIs33OcCSsOsPnIGepXr8wLd97EiM7NqOmvv/hdlec9aauUskxhkY1vtv/Mx3EH2ZOWRZPaVfnrPWE8ENUUfz+9f97VaeArpa4ot6CIhZtT+M/KQxw5lc0NDarz3oPtuLvddfj56JOw7kIDXyl1SWfzCpmTcJgpqxM5npVHu6BavNyvA3e0akilSnoThbvRwFdKXeT0uXym/5jEzB+TyMgpoOv19Xh/cARdr9e75dyZBr5SqkRGdgGxqw8yfU0S2flF9G7dkN/fegMRTWs7uzTlABr4SimycguYFp/ElPhDZOUW0i+8MWN7teDGhjWcXZpyIA18pbxYdn4hM388zH9WHeRMdgG9Wzfkj3fcSKvGNa+8s3I7GvhKeYP48fapiYufZs0tKGL5tws5sGUl72XfRc+bAnn2jhsJD9KhG0+mga+UNyieh75g4HQ+OxHC2mVf8NeCd9kc+AoLH+5Ch2ZXXnJTuT8NfKW8QEFwd+Lbvk3E7BGcKuzFP/x+ILn3JF7p1t/ZpakKpIGvlAcrsglfbU3lg2X7STpZlX/WuYuxOZ8hPZ4nTMPe62jgK+WBRITF29N4f9k+Dhw/S6vGNVl4ZyGR65dC9AuYjVOLpyfWGSq9iQa+Uh5mfeIp3ly8m63JZ2jRoDoTh0dyZ7V9VFr4R52H3stp4CvlIQ6mn+Uf3+5h6a5jNKrpzzuDwhkYGYRPJQPx83QeeqWBr5S7O3E2jw+W7WfO+iNU9fPh+T43MbpbKFUrnzd7ZfdnLt5Rh3S8jga+Um4qJ7+IqfGHmLTyEDkFRQzrGMzY21vowiPqkiwNfGNMbWAK0AYQYLSIrLWyTaU8XZFN+GJzCu8u2UdaZi53tG7IuL4tuT6wurNLUy7O6h7+B8B3IjLIGFMZqGZxe0p5tNX703lr8R52/5xJu6BafDAkgk7Ny7WonPIilgW+MaYmEA2MAhCRfCDfqvaU8mS7f87k79/uYdW+dILqVOXDoe3p37axzkmvroqVPfzmQDow3RjTDtgEjBWRc+dvZIyJAWIAgoODLSxHKfeTlpHLu0v2smBzCjWq+PLyXa14qGszqvjqcoLq6hkRsebAxkQB64BuIpJgjPkAyBSR/7vUPlFRUbJx40ZL6lHKneQWFBG76hAT4w5gs8FDXZrx1G03ULtaZWeXplyMMWaTiESVZVsre/gpQIqIJBS/XgCMs7A9pdyeiLB01zH++s0ukk/lcFfbRoy7sxXB9fTylyo/ywJfRNKMMcnGmJtEZC/QC9hlVXtKubuD6Wd5/X+7WLUvnRYNqvPpY53odkN9Z5elPIjVd+n8Afi0+A6dQ8AjFrenlNs5m1fIR8v3My0+EX9fH/6vf2se6tIMP59Kzi5NeRhLA19EtgBlGltSytuICF9tPcpbi3dzLDOPQR2CePHOlgTW0AenlDX0SVulnGDX0Uz+8tVO1iedIjyoFh+P6EBkcB1nl6U8nAa+Ulb6zdKCZ7LzWfTFXI7tWcuBKvfzj4FteTCqqd5PryqEBr5SVipeWrDo/ul8fjKU5d8u4B+29/iq5ZusGNiTWtX8nF2h8iIa+EpZKTSafbd8RMPZIzhe0It3Ky/n9N1TeCTqTmdXpryQ3gaglEXOZOfz4oJt9F4EC0xvxvouokb3GJpp2Csn0R6+Ug72y903f/16F6ezC3gz4jTDDq+Am3VpQeVcGvhKOdCRk9m8/N/trN5/goimtVl4ZxbNlr8KD87QpQWV02ngK+UABUU2psYnMn7ZPnwrVeL1AWGM6NwMnx8/0KUFlcvQwFeqnLYkn2Hcwm3sScuiT1hD/jIgjMa1qto/1KUFlQvRwFfqGmXlFvDukn3MXJtEwxr+/GdkB/qENXJ2WUpdkga+Utfg+51pvPblTo5l5fJwlxD+1PtGavjrPfXKtWngK3UVfs7I4bUvd7Jk1zFaNqrBxyMiaa9TIig3oYGvVBkU2YTZ6w7zzvd7KbTZGNe3JY92D9UZLZVb0cBX6gp2/5zJn7/YzpbkM/RoUZ83722rC5Iot6SBr9QlFBTZ+DjuIB/+sJ9aVf0YPziCeyKuwxid6Ey5Jw18pUqxJy2T5+ZvZUdqJgPaXcfrA8KoE6DrySr3ZmngG2OSgCygCCgs60K7SjlLQZGNSXEH+XC5vVc/aUQH7myjt1oqz1ARPfxbReREBbSjVLnsTcviuflb2Z6aQf/wxrxxTxvqaq9eeRAd0lFer7DIxn9WHeKDZfup4e/Lx8Mj6du2sbPLUsrhrA58AZYYYwT4j4jE/nYDY0wMEAMQHBxscTlKXWjfMXuvfltKBv3CG/PGgDDqVdc1ZZVnsjrwu4nIUWNMA2CpMWaPiKw6f4PiXwKxAFFRUWJxPUoB9l597OpDjF+6n+r+vkwYFkm/cO3VK89maeCLyNHi/x43xiwCOgKrLr+XUtbafyyL5xZsY2vyGe5q24g37mlDfe3VKy9gWeAbYwKASiKSVfx9b+ANq9pT6koKi2xMXp3I+8v2EVDZh38Pa0//8OucXZZSFcbKHn5DYFHxQyq+wBwR+c7C9pT6Vfx4+wLixdMQHzh+lhmfziLgxDZuazmGv97bhsAa2qtX3sWywBeRQ0A7q46v1GU1iYT5o5BB05mV1owfvl3AeJ/x7L71I7r2itSnZZVX0tsylWcKjeZ0v1h8Zo/kTP5tTKi8nML7Z9At7HZnV6aU0+hUf8ojLdt1jF5fCJ8U9mKs7yKqd4+hjoa98nIa+MqjZOcX8tKi7Tw2ayO9q+7j8WpxEP0CZuM0SNQbxJR30yEd5TG2p2Qwdt5PJJ44x5sRpxl25F3M4JnFa8j2gPmjLlxQXCkvoz185faKbMLEuAPcN3EN2XlFfPpoJ4YHncCcH+6h0fawT93szFKVcirt4Su3lnomh2fnbSEh8RT92jbmzfvaULtaZbjhmYs3Do3W3r3yahr4ym19tfUoLy/ajs0m/OuBdtwf2URvt1TqMi4b+MaYyDIco0BEtjuoHqWuKDO3gNe+3Mmin1KJDK7N+4MjaFYvwNllKeXyrtTDXwlsAC7XbQoFQhxVkFKXsyHpFM98toW0zFyeub0FT916A766kLhSZXKlwN8gIrddbgNjzHIH1qNUqQqLbHzww34mrDhAUJ1qfP54Fzo0q+PsspRyK5cN/CuFfVm3Uao8jp7J4em5P7Hx8GkGdQjiLwPCqF5FLz8pdbXK/LfGGBOOfeimZB8R+cKCmpQqsXzPMZ79fCsFhTY+GBLBPRFNnF2SUm6rTIFvjJkGhAM7AVvx2wJo4CtLFBTZ+Nf3e/nPqkO0alyTCcPa0zywurPLUsqtlbWH31lEWltaiVLFUs/k8Ic5m9l85AwjOgfzSr/W+Pv5OLsspdxeWQN/rTGmtYjssrQa5fWW7jrGc/O3UmQTXaBEKQcra+DPxB76aUAe9ts0RUTCr7SjMcYH2Aikikj/a65UebT8Qhv//G4PU+MTadOkJv8eGklIfb23XilHKmvgTwNGAtv5dQy/rMYCu4GaV7mf8hLJp7J5au5PbE0+w8NdmvFSv1ZU8dUhHKUcrayBf0REvrragxtjgoB+wJvAs1e7v/J83+1I4/kFWwH4eHgkfds2dnJFSnmusgb+HmPMHOB/2Id0gDLdljkeeAGocW3lKU+VV1jE3xfvYcaPSYQH1eLfQyMJrlfN2WUp5dHKGvhVsQd97/Peu+xtmcaY/sBxEdlkjOl5me1igBiA4ODgMpaj3NmRk9k8OWcz21MzeKRbCOP6ttQhHKUqgBERaw5szN+xj/sXAv7Yx/C/EJERl9onKipKNm7caEk9yjUs3v4zLy7YhjHwzgPt6BPWyNklKeXWjDGbRCSqLNtedtap4t73lRordRsR+bOIBIlICDAEWH65sFeeraDIxl++2snvP93M9Q2q883TPTTslapgVxrSGWeMOXGZzw32u3BiHVeS8jTHM3N5cs5mNiSd5pFuIfy5bysq++oMl0pVtLJMj3z3FbZZeqVGRCQOiCtbScqtxY+HJpElK0ttTDrF1E9m0qlgHyOGvKJz4SjlRFeaLfORiipEeYgmkTB/FDJoOjN/DmbZ4gX82+9DMu+ZTLCGvVJOpXPMKscKjSbv3mkUfDqSjLzb+LjKcsyDMwluqbNoK+VsGvjKoY6czObxb33om3cbY30XId2ex2jYK+US9MqZcpgVe4/T/6PVXHd6PU8ExEH0C5hN0yBxlbNLU0pR9vnwXy3tfRF5w7HlKHdkswkfLT/A+B/28WC9JN4q/Dc+g2fZL9yG9oD5o+CBGSUXcpVSzlHWIZ1z533vD/THPiGa8nIZOQU8O28LP+w5zsD2Tfhbo334BM/8NdxDo+1hn7pZA18pJytT4IvIu+e/Nsb8C7jqydSUZ9mTlsnjn2wi9XQOrw8I46EuzTAm4uINQ6M17JVyAdd60bYa0NyRhSj38uWWVMYt3E4Nf18+i+lMVEhdZ5eklLqCso7hb8c+WRqADxAI6Pi9FyoosvHW4t1MX5PEzSF1mDA8kgY1/J1dllKqDMrawz9/papC4JiIFFpQj3Jhp87l87vZm0hIPMUj3UJ46a5W+PnojV5KuYuyjuEftroQ5dr2pmXx2KwNHMvM4/3B7bivfZCzS1JKXSV98Epd0dJdx3jms58IqOLL5493IaJpbWeXpJS6Bhr46pJEhIlxB/nXkr20bVKL2JFRNKql4/VKuSsNfFWq3IIiXly4jS+3HOXudtfxzqBw/P10VSql3JkGvrrIscxcYmZtZGtKBs/3uYnf97weY4yzy1JKlZMGvrrA1uQzxHyykazcQmJHdqC3rkqllMewLPCNMf7AKqBKcTsLROQ1q9pT5fflllSeX7CNBjWq8MXvu9KyUU1nl6SUciAre/h5wG0ictYY4wfEG2O+FZF1FraproHNJvxryV4mxh2kY2hdPh4eSb3qVZxdllLKwSwLfBER4GzxS7/iL7n0HsoZzuYV8sxnW1i2+xhDOwbz+oAwXW9WKQ9l6Ri+McYH2ATcAEwQkYRStokBYgCCg4OtLEf9RvKpbB6buZED6WfPm/xML84q5aks7cqJSJGIRABBQEdjTJtStokVkSgRiQoMDLSyHHWedYdOMuDf8aRl5jLzkY483DVEw14pD1ch/3YXkTNAHHBnRbSnLm9OwhFGTEmgbkBl/vtkN7q3qO/skpRSFcCywDfGBBpjahd/XxW4HdhjVXvqymw24e+Ld/PSou10u6E+i57sRmj9AGeXpZSqIFaO4TcGZhaP41cCPheRry1sT11GbkERf5y3hW93pPFQl2a82r81vjrTpVJexcq7dLYB7a06viq7E2fzGDN8NIEcAAAOmElEQVRrI1uSz/BKv1Y82j1Ux+uV8kL6pK2HO3D8LI/MWE96Vh4fD+/AnW30yVmlvJX+m95TxI+HxFUXvLXrx6/538QXyMkv4rOYLhr2Snk5DXxP0SQS5o8qCf3VSxbS6PsnSKnWikW/76Zz2CuldEjHY4RGwwMzkPmjSKh3L62PfM6E+q/w6qOPUquan7OrU0q5AO3he5D8pt1ZWq0fnZOnsinwPl78XYyGvVKqhAa+h8jILuCfH8fSIX0RCU0f447sb6icHO/sspRSLkQD3wMkn8rmjY8+5skTf2NXtw/o9Oi7mAdmXDCmr5RSGvhu7qcjp7lv4hqCcvaS1mcSPXrfb/+geEyf1M1OrU8p5Tr0oq0b+27Hz4z9bAsNalbh7ph/ckOD6hduEBpt/1JKKTTw3ZKIMDU+kTcX76ZdUG2mPBxFfV2wRCl1BRr4bsZmE/72zW6mrUmkb5tGvD84An8/H2eXpZRyAxr4biS/0MZz87fy1dajjOoawqv9W1Opks6Jo5QqGw18N3E2r5AnPtlE/IETvHDnTfzulut1AjSl1FXRwHcD6Vl5jJ6xgV0/Z/LOoHAeiGrq7JKUUm5IA9/FHT55joemredYZi6TH+rAbS0bOrskpZSb0sB3YTtSMxg1fT2FNmHOmM5EBtdxdklKKTdm5RKHTY0xK4wxu40xO40xY61qyxPF7z/B4P+spYqvDwue6KJhr5QqNyt7+IXAn0RkszGmBrDJGLNURHZZ2KZH+GrrUf70+Raa16/OzNEdaVTL39klKaU8gGU9fBH5WUQ2F3+fBewGmljVnqeYFp/I03N/on3TOnz+RBcNe6WUw1TIGL4xJgT7+rYJpXwWA8QABAcHV0Q5LklEePv7vXwcd5A+YQ35YEh7faBKKeVQlk+eZoypDiwEnhGRzN9+LiKxIhIlIlGBgYFWl+OSCopsPDd/Gx/HHWRYp2AmDu+gYa+UcjhLe/jGGD/sYf+piHxhZVvuKju/kCc/3cyKvek8c3sLxvZqoQ9UKaUsYVngG3tqTQV2i8h7VrXjzk6dy2f0jA1sSznDm/e1YXinZs4uSSnlwazs4XcDRgLbjTFbit97SUQWW9im20g9k8PIqQmknM7h4xEd6BPWyNklKaU8nGWBLyLxgI5NlOJQ+llGTEkgK6+QT0Z3pFPzes4uSSnlBfRJ2wq262gmD01LQATmjulMmya1nF2SUspLaOBXoE2HT/PI9PUEVPHlk0c7XbxClVJKWUgDv4KsOXCCMbM20qBGFWY/1omgOtWcXZJSysto4FeAJTvTeGrOT4TWD+CTxzrSoIY+PauUqnga+Bb770+p/Gn+Vto0qcXMR26mdrXKzi5JKeWlNPAt9Mm6w7z65Q46h9Zj8sNRVK+if9xKKefRBLLIxLgDvP3dXnq1bMCE4ZE6VYJSyuksn0vHK8SPh8RVQPEkaN/tYfWSLxgftJJJI3VeHKWUa9AeviM0iYT5o7DdP52/7KjH/oTFTK46gap9Z+Hjo79TlVKuQQPfEUKjKbp/OjlzRlIv7zbGVVtB1WGzMM1vcXZlSilVQgPfAfIKi/jDmgDC8m5jrO8ipOvzGvZKKZej4w3llJ1fyKMzNpK1ZzmPV42D6BcwG6eVjOkrpZSr0B5+OWTkFDB6xgaqJMczo/pEqgydBaHRENoD5o+CB2bYXyullAvQwL9Gp8/lM3JaAnvTsvhfZB5VOsz6NdxDo+1hn7pZA18p5TI08K/BybN5jJi6noPpZ4kdGUXLlnddvFFotIa9UsqlaOBfpRNn8xg+OYGkk+eY8lAU0Td65zq8Sin3Y9lFW2PMNGPMcWPMDqvaqGjHM3MZEruOI6eymT7qZg17pZRbsfIunRnAnRYev0KlZdjD/uiZHGY8cjNdb6jv7JKUUuqqWBb4IrIKOGXV8SvS0TM5DI5dy/GsPGbpkoRKKTfl9PvwjTExxpiNxpiN6enpzi7nIsmnshkcu5ZTZ/OZ9WhHokLqOrskpZS6Jk4PfBGJFZEoEYkKDHStMfEjJ7MZEruOjOwCPh3TicjgOs4uSSmlrpnepXMJiSfOMTR2HbmFRczRxcaVUh5AA78UB46fZdjkdRTahDmPdab1dTWdXZJSSpWblbdlzgXWAjcZY1KMMY9a1ZYj7TuWxZDYddhEmDtGw14p5Tks6+GLyFCrjm2V3T9nMmJKApUqGeaO6cwNDWo4uySllHIYp1+0dRU7j2YwbPI6/HwqMS9Gw14p5Xl0DB/YnpLBiKkJBFT2YW5MZ5rVC3B2SUop5XBe38PfknyGYVPWUcPfl3mPd9GwV0p5LK/u4W9NPsPIKQnUCajM3JjONKld1dklKaWUZby2h78jNYORUxOoVc1Pw14p5RW8MvB3Hc1kxNQEavj7MXeMhr1Syjt4XeDvTctixNQE/H19mDOmE03rVnN2SUopVSG8KvAPHM9i+JR1+FYyejeOUsrreE3gH0o/y9DJCYBhzpjOhNbXsFdKeRevCPzDJ88xbHICNpswd0wnbmhQ3dklKaVUhfP4wE8+lc3Q2HXkFRbx6ZhOtGioT9AqpbyTR9+Hn3omh6GT13Euv4g5YzrRspFOhKaU8l4e28P/OSOHobHryMgpYPajnQi7TuezV0p5N48M/GOZuQybnMCpc/nMGt2RtkEa9kop5XGBn56Vx7DJ6ziemcvM0TfTXpclVEopwMPG8E+etYf90TO5zBzdkQ7NdMFxpZT6haU9fGPMncaYvcaYA8aYcQ5vIH48JK4C4PS5fIZPSeC60xv49uZNdAzVsFdKqfNZucShDzAB6Au0BoYaY1o7tJEmkTB/FGd3L2fE1AQCT65nSrUJhLTt4dBmlFLKE1g5pNMROCAihwCMMZ8B9wC7HNZCaDTZ90yh6LOH6Vt4OzHVVuA3ZBaERjusCaWU8hRWDuk0AZLPe51S/N4FjDExxpiNxpiN6enpV92I3w09WVf3Xp7y+YLKncZo2Cul1CVYGfimlPfkojdEYkUkSkSiAgMDr7oRvyPx9Mn5BqJfgI1TS8b0lVJKXcjKwE8Bmp73Ogg46tAWElfB/FHwwAy47WX7f+eP0tBXSqlSWBn4G4AWxphQY0xlYAjwlUNbSN1sD/lfhnFCo+2vUzc7tBmllPIEll20FZFCY8xTwPeADzBNRHY6tJHuz1z8Xmi0juMrpVQpLH3wSkQWA4utbEMppVTZeNzUCkoppUqnga+UUl5CA18ppbyEBr5SSnkJI3LRs1BOY4xJBw5f4+71gRMOLMcd6Dl7Pm87X9BzvlrNRKRMT626VOCXhzFmo4hEObuOiqTn7Pm87XxBz9lKOqSjlFJeQgNfKaW8hCcFfqyzC3ACPWfP523nC3rOlvGYMXyllFKX50k9fKWUUpehga+UUl7C5QP/SguhG2OqGGPmFX+eYIwJOe+zPxe/v9cY06ci6y6Paz1nY0w9Y8wKY8xZY8y/K7ru8ijHOd9hjNlkjNle/N/bKrr2a1WOc+5ojNlS/LXVGHNfRdd+rcrz97n48+Di/7+fq6iay6scP+cQY0zOeT/rSeUuRkRc9gv7tMoHgeZAZWAr0Po32/wemFT8/RBgXvH3rYu3rwKEFh/Hx9nnZPE5BwDdgSeAfzv7XCronNsD1xV/3wZIdfb5VMA5VwN8i79vDBz/5bUrf5XnnM/7fCEwH3jO2edTAT/nEGCHI+tx9R5+yULoIpIP/LIQ+vnuAWYWf78A6GWMMcXvfyYieSKSCBwoPp6ru+ZzFpFzIhIP5FZcuQ5RnnP+SUR+WUltJ+BvjKlSIVWXT3nOOVtECovf96eUpUNdVHn+PmOMuRc4hP3n7C7Kdc6O5uqBX5aF0Eu2Kf5LkAHUK+O+rqg85+yuHHXO9wM/iUieRXU6UrnO2RjTyRizE9gOPHHeLwBXds3nbIwJAF4EXq+AOh2pvP9vhxpjfjLGrDTG9ChvMZYugOIAZVkI/VLblGkRdRdUnnN2V+U+Z2NMGPBPoLcD67JSuc5ZRBKAMGNMK2CmMeZbEXH1f9mV55xfB94XkbMWdX6tUp5z/hkIFpGTxpgOwH+NMWEiknmtxbh6D78sC6GXbGOM8QVqAafKuK8rKs85u6tynbMxJghYBDwkIgctr9YxHPJzFpHdwDns1y9cXXnOuRPwtjEmCXgGeKl4CVVXd83nXDwcfRJARDZhvxZwY3mKcfXAL8tC6F8BDxd/PwhYLvYrHl8BQ4qvgIcCLYD1FVR3eZTnnN3VNZ+zMaY28A3wZxFZU2EVl195zjm0OBgwxjQDbgKSKqbscrnmcxaRHiISIiIhwHjgLRFxhzvRyvNzDjTG+AAYY5pjz7BD5arG2Vexy3CV+y5gH/bfbi8Xv/cGMKD4e3/sV+0PYA/05uft+3LxfnuBvs4+lwo65yTsPaKz2HsOrSu6/oo8Z+AV7D3cLed9NXD2+Vh8ziOxX7jcAmwG7nX2uVh9zr85xl9wk7t0yvlzvr/457y1+Od8d3lr0akVlFLKS7j6kI5SSikH0cBXSikvoYGvlFJeQgNfKaW8hAa+Ukp5CQ18pZTyEhr4SpXi/Klpr3K/wcXT3H5tVW1KXSsNfKUu7aCIRFzNDiIyD3jMonqUKhcNfOV1jDE3G2O2GWP8jTEBxpidxpjLzkVT3OPfY4yZYozZYYz51BhzuzFmjTFmvzHGHabeVl7O1WfLVMrhRGSDMeYr4G9AVWC2iOwow643AA8AMdjnSBmGfcGZAcBLwL3WVKyUY2jgK2/1BvbQzgWeLuM+iSKyHaB4LvofRESMMduxr06klEvTIR3lreoC1YEa2CevKovzF1axnffahnaelBvQwFfeKhb4P+BT7AunKOXxtFeivI4x5iGgUETmFM83/qMx5jYRWe7s2pSykk6PrFQpjDEhwNcictUrSRljemKfr72/g8tSqlx0SEep0hUBta7lwStgInDakqqUKgft4SullJfQHr5SSnkJDXyllPISGvhKKeUlNPCVUspL/D98h5gHcbYIlgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f96a35c5390>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# --------------------\n",
    "# Exact solution\n",
    "# --------------------\n",
    "x_ex = np.linspace(0, L, 10)\n",
    "u_ex = [-0.5*b*x_ex_i**2/E/A + (g+b*L)*x_ex_i/E/A for x_ex_i in x_ex]\n",
    "\n",
    "# --------------------\n",
    "# Post-process\n",
    "# --------------------\n",
    "fe.plot(u)\n",
    "plt.plot(x_ex, u_ex, \"x\")\n",
    "plt.xlabel(\"x [m]\")\n",
    "plt.ylabel(\"u [m]\")\n",
    "plt.legend([\"FEM solution\",\"exact solution\"])\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Complete code\n",
    "\n",
    "---"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import fenics as fe\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "\n",
    "# --------------------\n",
    "# Parameters\n",
    "# --------------------\n",
    "E = 0.025e9  # Young's modulus\n",
    "A = np.pi*1e-4  # Cross-section area of bar\n",
    "L = 0.05  # Length of bar\n",
    "n = 20  # Number of elements\n",
    "b = 0.03  # Load intensity\n",
    "g = 0.0005 # External force\n",
    "\n",
    "# --------------------\n",
    "# Geometry\n",
    "# --------------------\n",
    "mesh = fe.IntervalMesh(n, 0.0, L)\n",
    "\n",
    "fe.plot(mesh)\n",
    "plt.show()\n",
    "\n",
    "# --------------------\n",
    "# Function spaces\n",
    "# --------------------\n",
    "V = fe.FunctionSpace(mesh, \"CG\", 1)\n",
    "u_tr = fe.TrialFunction(V)\n",
    "u_test = fe.TestFunction(V)\n",
    "\n",
    "# --------------------\n",
    "# Boundary marking\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",
    "dx = fe.Measure(\"dx\", mesh)\n",
    "ds = fe.Measure(\"ds\", subdomain_data = boundary)\n",
    "        \n",
    "# --------------------\n",
    "# Boundary conditions\n",
    "# --------------------\n",
    "bc = fe.DirichletBC(V, 0.0, left)\n",
    "\n",
    "# --------------------\n",
    "# Weak form\n",
    "# --------------------\n",
    "a = E*A*fe.inner(fe.grad(u_tr),fe.grad(u_test))*dx\n",
    "l = b*u_test*dx + g*u_test*ds(2)\n",
    "\n",
    "# --------------------\n",
    "# Solver\n",
    "# --------------------\n",
    "u = fe.Function(V)\n",
    "fe.solve(a == l, u, bc)\n",
    "\n",
    "# --------------------\n",
    "# Exact solution\n",
    "# --------------------\n",
    "x_ex = np.linspace(0, L, 10)\n",
    "u_ex = [-0.5*b*x_ex_i**2/E/A + (g+b*L)*x_ex_i/E/A for x_ex_i in x_ex]\n",
    "\n",
    "# --------------------\n",
    "# Post-process\n",
    "# --------------------\n",
    "fe.plot(u)\n",
    "plt.plot(x_ex, u_ex, \"x\")\n",
    "plt.xlabel(\"x [m]\")\n",
    "plt.ylabel(\"u [m]\")\n",
    "plt.legend([\"FEM solution\",\"exact solution\"])\n",
    "plt.show()"
   ]
  }
 ],
 "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
}