commit 5aada84c7c014ae273c88bc93283467593b1e3f4 Author: Gerardo Marx Date: Wed Apr 21 12:51:10 2021 -0500 Heat Equation, Dirichlet boundary condi diff --git a/.gitignore b/.gitignore new file mode 100644 index 0000000..763513e --- /dev/null +++ b/.gitignore @@ -0,0 +1 @@ +.ipynb_checkpoints diff --git a/Poisson.ipynb b/Poisson.ipynb new file mode 100644 index 0000000..39d0971 --- /dev/null +++ b/Poisson.ipynb @@ -0,0 +1,219 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Introduction\n", + "The Poisson's equation is a second-order partial differential equation that stats the negative Laplacian $-\\Delta u$ of an unknown field $u=u(x)$ is equal to a given function $f=f(x)$ on a domain $\\Omega \\subset \\mathbb{R}^d$, most probably defined by a set of boundary conditions for the solution $u$ on the boundary $\\partial \\Omega$ of $\\Omega$:\n", + "\n", + "$$-\\Delta u =f \\quad \\text{in } \\Omega\\text{,}$$\n", + "$$u=u_0 \\quad \\text{on } \\Gamma_D \\subset \\partial\\Omega \\text{,}$$" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "here the Dirichlet's boundary condition $u=u_0$ signifies a prescribed values for the unknown $u$ on the boundary.\n", + "\n", + "The Poisson's equation is the simplest model for gravity, electromagnetism, heat transfer, among others.\n", + "\n", + "The specific case of $f=0$ and a negative $k$ value, leaves to the Fourier's Law." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Comparative analysis\n", + "Along this example, the fenics platfomr is used to compare results obtained by solving the heat equation (Laplace equation) in 2-D:\n", + "\n", + "$$\\frac{\\partial^2 T}{\\partial x^2}+ \\frac{\\partial^2 T}{\\partial y^2}=0$$\n", + "\n", + "the problem is defined by the next geometry considerations:\n", + "\n", + "![](physicalproblem.png)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "The resulting contour of temperature, solving using finite diferences, is shown next:\n", + "\n", + "![](resulteq.png)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Solving by Finite Element Method with Varational Problem formulation" + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "metadata": {}, + "outputs": [], + "source": [ + "#1 Loading functions and modules\n", + "from fenics import *\n", + "import matplotlib.pyplot as plt" + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "metadata": {}, + "outputs": [ + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAQgAAAD8CAYAAACLgjpEAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJztnXt0lHWa5z9PBRIqgEgSCNcQJRC5hYCAoCDYEVoRETx4CYjx7MxxZ7f7nJk9vWe3154dZ3p25kzvbk/39KxnbMemRbmJCAK2tiJeuNtAuDWXAGKAhJgYRAMBpCHP/lEXK2UKivf3i6Te/D7n5OStt9566ve+VfnkvVR9H1FVHA6HoyUCN3oADoej7eIE4XA4EuIE4XA4EuIE4XA4EuIE4XA4EuIE4XA4EuIE4XA4EuIE4XA4EuIE4XA4EtLhRg+gJXJycjQ/P/9GD8Ph8C07d+6sV9Ue11quTQoiPz+fHTt23OhhOBy+RUSOJ7OcO8RwOBwJcYJwOBwJcYJwOBwJcYJwOBwJcYJwOBwJuaYgRKS/iHwgIgdEZL+I/GV4fpaIrBORI+Hf3RM8viy8zBERKbO9Ag6Ho/VIZg/iMvAjVR0KjAd+ICJDgR8D61V1ELA+fLsZIpIFPAvcAYwDnk0kEofD0fa45ucgVLUGqAlPnxWRg0Bf4CFgSnixhcCHwH+Pe/j3gXWq+gWAiKwD7gOWmgx6z549vPHGGwwfPpysrCyTUgBs2LABgEmTJiEi1urdfffdxrVUlY0bN1qrd+zYMaqqqsjKymL48OHG9SLrescdd5CRkWGtno11tV2vtraWiooK0tPTGT9+vHG9yNhGjBhB9+7m/zcPHTrE+fPn+dGPfmRcK4JcTyaliOQDG4DhwAlVvTk8X4Azkdsxy/9XoJOq/q/w7f8JXFDV/9tC7aeBpwHy8vJuP3488ec4fvGLX9DQ0JD0uB2O9sSzzz57zWVEZKeqjrnWckl/klJEugCvA3+lqg2x/2lVVUXEKP1WVV8AXgAYM2bMVWs98MADLF0a2gm55ZZbKC0tpWPHjp6f+7333mPz5s0ATJgwgalTpxrtSfz617/ms88+Q0SYM2cOQ4cO9VwL4O/+7u8AyMzMpKysjJ49e3qudfLkSRYsWABA3759eeKJJ+jUqZPneps3b+a9994DYNSoUTz44ING227RokV88sknADz44IOMHj3acy34ZttlZGQwf/58+vbt67lWfX09zz33HAA9evSgrKyMzp07e65XXl7O2rVrARg6dCgPP/wwaWlpnuutXr2aY8eOeX58SyR1FUNEOhKSw2JVXRmeXSsivcP39wbqWnhoNdA/5na/8DwrjBs3jsrKSpYuXcqf/vQnz3VEhEAgwNixY9m6dSvr1q3DJO07PT2dPn360L9/f1asWMGBAwc814qMb9iwYQQCARYuXEhdXUub+voYO3YsNTU1LFq0iIsXLxqNDWD8+PHs2rWLtWvXGm27Dh06kJ2dTUFBAWvXrqW8vNxzLYBgMEhhYSGZmZm88sorVFebv/1Gjx7NmTNnWLhwIY2Njcb1xo8fz4EDB1i5ciVXrlwxrmeTZK5iCPAb4KCq/nPMXWuAyFWJMmB1Cw9/B5gmIt3DJyenhedZoaioiFmzZlmRBMD9999vVRJz5861JomsrCyeeuopa5IYPHgwjz76qBVJAEyZMoW7777bmiQee+wxa5K46aabKCsrsyaJ/Px85s2bZ00S48ePZ9q0aW1SEsnsQdwFzAe+JyK7wz/TgX8CporIEeDe8G1EZIyIvAgQPjn598D28M9PIycsbWFTEiJiVRIZGRlWJZGdnW1VEoWFhdYkISJtWhLdunVr05KYMGFCm5TENQWhqptUVVS1SFWLwz9vqeppVS1R1UGqem/kD19Vd6jqn8c8foGqFoR/ftsaK+Ek4R0nCe+0B0n45pOUThLecZLwjt8l4RtBgJOECU4S3vGzJHwlCHCSMMFJwjt+lYTvBAFOEiY4SXjHj5LwpSDAScIEJwnv+E0SvhUEOEmY4CThHT9JwteCACcJE5wkvOMXSfheEOAkYYKThHf8IIl2IQhwkjDBScI7qS6JdiMIcJIwwUnCO6ksiXYlCHCSMMFJwjupKol2JwhofUmY0NqS+Pzzz43qOUl4p7Ul0dTUZFSvJdqlIKC5JD7++GOjjRsviRMnThiNLV4SJn800FwSkYASE2IlEYnD80q8JA4fPmxVEhcuXDAaX6wk3njjDaNa0FwS77//vnG9WEns3bvXuiSuK3Luu2LMmDF6td6cP//5zzl37hwQSpQy4dNPP41O5+bmkpmZ6bmWqlJZWRm9nZ+fb5Su9PXXX3Pq1KnobdN1PX36dLOoPpvbLi0tjby8PM+14rfdgAEDCAS8//+6fPkyJ0+ejN42XdevvvqKL774JqnAtF5lZWUzEdp8Lf76r//6mslUyUbOpeQexOXLl6PT9fX1XLlyxfNP7Ju6traWixcveq7V1NRE//7fBGhVVlZy+fJlz/U6dOjQLMz0008/NVrXm29uFhlKTU2NtW135coVGhsbrW2748ePG207EWkWzWe67bp06dJs2504ccKoXuy6AjQ0NFh7LWxyzUxKEVkAzADqVHV4eN6rQGF4kZuBL1W1uIXHVgJngSvA5WSMlQyzZ8+OZlIGAgFmzZpllAq8fv16Nm3aBHyze2+ScblgwYLof69+/foZZ1xGchVFhHHjxjFs2DDPtWIzKQOBAA888IBRxuWWLVui5106duxIaWmpUcblsmXLqKioACAnJ8c44/If//Efo+eYRo4caZRxGZtJ2aFDB6ZNm2aUcRmbSRkIBHj00UeNMi5Xr17NJ598YpRrGU8yexAvEYqqj6Kqj0XCYwhlVa5s6YFh7gkva0UOscyYMYNLly6xcOFCzpw5Y1QrEAgwe/ZsKisrWbJkifGJy/z8fGtXN0SEO+64g379+vH666+zf/9+z7Ui3HfffVYzLh9++GFr8XW5ublMmjTJ2onLUaNGWTtxCVBSUmI143LWrFnWTlzaaNsQSzKJUhuAFmPiwnmVj2LY58IrvXr1Yv78+dYkUVRUZE0SYD/jct68edYk0ZYvgQLcc889ViXRnq5u2MT0HMQkoFZVjyS4X4F3RWRnuO+FdXr37t1mJdEan5NoL5IQESeJNoCpIEq5+t7DRFUdDdxPqGVfwvZGIvK0iOwQkR3Xe63eScI7ThLeaQ+S8CwIEekAPAy8mmgZVa0O/64DVhHqz5lo2RdUdYyqjunRo8d1j8dJwjtOEt7xuyRM9iDuBQ6palVLd4pIZxHpGpkm1BPjjwbPd02cJLzjJOEdP0simcY5S4GtQKGIVInIn4Xvepy4wwsR6SMib4Vv5gKbRGQP8Afgd6r6e3tDbxknCe84SXjHr5JI5ipGqar2VtWOqtpPVX8Tnv+Uqj4ft+wpVZ0enj6mqiPDP8NU9R9aZxW+jZOEd5wkvONHSaTkJymTwUnCO04S3vGbJHwrCHCSMMFJwjt+koSvBQFOEiY4SXjHL5LwvSDAScIEJwnv+EES7UIQ4CRhgpOEd1JdEu1GEOAkYYKThHdSWRLtShDgJGGCk4R3UlUS7U4Q0L4l4feMy9aUxK5duzzXgtSURLsUBDSXxP79+42z/GIlUVNTw/nz5z3XipeEKbGS2LZtm3G9WEns3r3buF6sJI4fP85XX33luVa8JC5evGj02sZKYsuWLZ7rRIiVhI1silhJHDlyxLokUjKTMpKwBBin51y50jwuvC3VU9Vmb+5AIGAUCNKW19V2vfa87Z555plrJqIlm0l5zci5tsiQIUM4ePAgAD179uTWW281qrd58+bo9G233fat7EaTerfffrtRfJ2qRv9zNTU1cdddd1kbW7du3RgyZIi1egMHDsTLN3ET1SsuLjaKr4ut19TUxJ133mkkiT179kTDkjMyMhg1apSVsQHk5eXRp08fK/VM3m/xpKQgiouLo4Koq6ujpKSEgQMHeq4nImzatIlgMEhVVRUlJSVGGZcnT57kxIkTiAiff/45paWlRi9a7K5tU1OTUcZlYWFhNJPyzJkz9OnTh6FDh3oeW2ZmJuvWrSM9PZ3q6mpKSkqMMi7r6+upqKggEAjw2Wef8cQTTxhJory8PBp9f+HCBaOMy+Li4mgm5fnz58nKyjLKuMzKymLt2rUEAgFOnTpFSUmJUcZlY2Mjx44d8/z4lkjpcxDz5s2jR48eLF26lE8++cSoViAQsHriMj8/32pznkmTJlk7cQnwyCOPWGvOA/DUU08hIrz88svGVzdyc3OtxteNGTPGWnMegJkzZ1rNuHzqqaesZlzaJKUFEQwGmT9/vjVJtMbVDVuSAPsZl63RwcuWJNpTB6+bbrrJ6tUNm6S0ICC0i9teJJEqvUCdJK4f25dAbZHyggAnCSeJ5HCSuH6SSZRaICJ1IvLHmHl/KyLVIrI7/DM9wWPvE5EKETkqIj+2OfB4nCScJJLBSeL68NQ4J8wvIs1zVPWt+DtFJA14jlCi9VCgVES8ny5PAicJJ4lkcJJIHqPGOddgHHA0HD13CVgGPOShznXhJOEkkQxOEslhcg7ihyKyN3wI0tKHBvoCJ2NuV4XntTpOEk4SyeAkcW28CuLfgIFAMVAD/Nx0ICaNc1rCScJJIhmcJK6OJ0Goaq2qXlHVJuDfabkhTjUQ2+O8X3heoppGjXNawknCSSIZnCQS40kQItI75uZsWm6Isx0YJCK3iEg6oT4aa7w8nwlOEk4SyeAk0TJeG+f8bxHZJyJ7gXuA/xJeNto4R1UvAz8E3gEOAstV1bxvvQecJJwkksFJ4tt4apyjqvNVdYSqFqnqTFWtCS8bbZwTvv2Wqg5W1YHfZeOclnCScJJIBieJ5vjik5TJ4iThJJEMThLf0K4EAU4SJjhJtD9JtDtBQOpJwuSNmWqSaGsZl7GSePPNN43G1tYzLlsiJQNjIpmFx48fN0oIKikpYfHixdTU1ABw6tQpo3GVlJRE30Rnz541qpeTk8O4ceP4+OOPAdi4cSO33Xab53ojR47k1KlTnDwZ+uzawYMHyczM9FxvypQpLFmyJPqmPHnyJMFg0HO9adOmsXLlymimosm269q1KxMnTmTDhg0AbN++neLiYs/1Bg8ezKlTpzh69CgA+/btIzs723O9SZMmcerUKaqqqgCorKw0SuKaOnUqy5cvN5ZrS6R8JqXD4WhOu8+kvOeee/jggw8A6NevH5MmTTKqt3Tp0uj02LFjKSgosFbvvvvuM4qvi683Z84co/i6vXv3Rvtj5OTkcO+99xrthcWObeTIkUbxdfH1TOPr4uvNnj3bKL6usrIymjLeuXNno/g6gOXLl0fDZm+77TbjjMvIurb7TMpevXoBoePXqqoqamtrjSQxceJENm3aRG5uLuXl5QwePNhIEnl5edTV1SEibN26lbKyMiNJiAjdunXjyy+/ZMeOHcydO9fzmyAYDLJ//36ysrKor6/n+PHjRhmXU6dOZd26deTm5rJ3714GDx5sJInCwkI++eQTOnXqFN12JpIIBoN07NiRc+fO8Yc//MEo4zIrK4utW7fSvXt3zpw5Q0VFhZEkpk+fztq1a+nRoweHDh1i0KBBRhmXxcXFLpMyllmzZjFixAjef/99Nm7caFQrEAjw5JNPkpOTw7Jly6LHm17p1auX1ROXI0aMsNacB+zG1wHRvhu2TlyWlZW1yQ5egNXmPIDVqxu2SWlBiIhVSWRmZlqVRFvu4AV2JRHbnMeGJHJyctqNJGxfArVJSgsCQv/5nSS80Zpt/vwuibbeC9QWKS8IcJJwkkgeJ4nrwxeCACcJJ4nkcZJIHt8IApwknCSSx0kiOXwlCHCScJJIHieJa+M7QYCThJNE8jhJXB1fCgKcJJwkksdJIjFeG+f8HxE5FE61XiUiNyd4bGU4eWq3iCT+ckUr4SThJJEsThIt47VxzjpguKoWAYeB/3GVx98Tbq5zzS+GtAZOEk4SyeIk8W08Nc5R1XfDmZMA2wglVrdZnCScJJLFSaI5Ns5B/Afg7QT3KfCuiOwUkaevVsR2X4x4nCScJJLFSeIbjAQhIj8BLgOLEywyUVVHE+rP+QMRuTtRrdboixGPk4STRLI4SYTwLAgReQqYAczTBGurqtXh33XAKlpusPOd4iThJJEsThLeG+fcB/w3YKaqnk+wTGcR6RqZBqbRcoOd75x4SWzatMmoXrwkTpw4YVQvXhKmXyeOlURsgIoX4iXx3nvvGdWLl0RFRYVRvXhJXLhwwahevCRMiJeE7YzL3bt3G9Vr8TmutUC4cc4UIEdEqoBnCV21yADWhcMytqnqX4hIH+DFcG+MXGBV+P4OwBJV/b2NQR8+fBiAd999lyFDhniu06tXL/bt2wdAU1MT27ZtMxrX4MGDqa2tBULpQzbq7dmzBwhlUppkSEJIPJH8zcWLF/P973/fc634AJwNGzbQpUsXz/UGDhwYzcusra013naFhYXs3LkTCGVSZmVlGdXLy8ujsrISoFl2pheCwSDp6elcunQJCL2P+/fvf41HJWbAgAHRPdeGhgbPdVrCZVI6HD6j3WdSzpkzhxUrVgCQn5/PnDlzSEtL81xvzZo1HDx4EIA777zTOOPyZz/7WbOxDhw40Eq9jIwM4/i6Tz/9lOXLlwOhPYrS0lKjDMP169cTkfmoUaOM4usA/uVf/iV6rP/ggw8aZ1xGtl1aWhpPPvmkUXxdXV0dv/3tb4FQ6pVJfB3Atm3b+OijjwAYOnQoM2bMMNp2y5Yt4/jx4y6TMrIBRo4cyZ49e3jzzTeNJJGdnY2IMHz4cLZs2UKnTp2MJJGXl8fXX3+NiLBq1SpKS0uNJCEiFBQUUFVVxauvvmokichhQFFREfv27YuOz+ubKjKOkSNHsmvXLjp16mQkiQEDBlBXV0eXLl1488036dSpk5EkgsEgvXv3pq6ujuXLlxtJInKIN3z4cA4cOMCKFSuMJHHTTTcBoSzJ3bt3k5GRYZRxGcnKtElKfxdj7Nix3H///Rw6dIgVK1ZEE4K9YDu+LhgMWm3O0xoZl7Y6eIHd+Lr09HTrGZc2O3gNHjzYanydzQ5etklpQQCMGzfOmiRa4xJoKnXw8vMl0Lbc5g/ariRSXhDgJGGCk4R32nIvUFv4QhDgJGGCk4R3/C4J3wgCnCRMcJLwjp8l4StBgJOECU4S3vGrJHwnCHCSMMFJwjt+lIQvBQFOEiY4SXjHb5LwrSDAScIEJwnv+EkSvhYEOEmY4CThHb9IwveCACcJE5wkvOMHSbQLQYCThAlOEt5JdUm0G0GAk4QJThLeSWVJJCWIBL0xskRknYgcCf9u8euFIlIWXuaIiJTZGrhXnCS84yThnVSVRLJ7EC/x7d4YPwbWq+ogYH34djNEJItQAtUdhPIon00kku8SJwnvtLYkTGhtSZimraeiJJISREu9MYCHgIXh6YXArBYe+n1gnap+oapnCDXciRfNDSFWElu3bqWpqclzrXhJmGZSxkvC9IWPlcSqVauMakFzSXz44YdGteIlYZpJGS8J00zKWElEgnZMiJWEqRDjJbFnzx7rkjAJjMlV1Zrw9GeEMijj6QucjLldFZ5nxPr16wF48cUXjTIp41m8eLFRGk/si1NZWWn8hurUqVN072bjxo3U19cb14v8wSxevNh428Wu7/PPP0+vXr2M6kWora013nbp6enR8W3fvp1z584Z1cvMzIzWWLlyZTSBzAa//OUvrb2Pz549i6oaJVPFYiVRSlVVRIzUFW6s8zSEEpmuRuzxYGVlpVFYao8ePaK7jkePHiUrK8sovi47O5vTp08DcPDgQUx7fKSlpUUlYVqvQ4fmL/eRI0eM4utit92ZM2dQVSPBxtazse2CwWBUiDbqxWJaLzs7u9khS1VVlVF8Xey2u3Llyrdea6+YVKkVkd6qWiMivYGWzuJUE0rEjtAP+LClYqr6AvAChEJrr/bEpaWl0fj2yO54165dr3sFIqxfvz4afd+jRw8eeeQRI0m8+OKLVFdXA6HkJtOMy0hIb1paGtOmTaOgoMBzrZMnT7JgwQIAOnfuTGlpqZEktmzZEt1V7t69u3HG5dKlS6Op5QUFBcYZl5FtF9kdN4mvq6+v57nnngNC77s5c+YYZVyWl5ezdu1aIBQ/Z5pxuXr1ao4ePWpNDmB2mXMNELkqUQasbmGZd4BpItI9fHJyWnieFaZNm8bZs2dZuHAhZ8+eNaoVCAS4//77qaio4LXXXjM6cZmWlsaAAQMoKiqycuJSRBgzZgw9evSw0pwHoKSkxNqJS4AHHnjA2onL3Nxca1c3gsEgI0eOpH///lZOXAJMnjzZWnMegOnTp1tLpgoE7H5yIdnLnEuBrUChiFSJyJ8B/wRMFZEjwL3h24jIGBF5EUBVvwD+Htge/vlpeJ4V8vLymDdvnjVJRE5c2pCEiPDQQw9Zk0RsxqUNSfTq1Ysnn3yyTV7dAPsZl3PnzrUmiciJS1uSsJ1xaZNkr2KUqmpvVe2oqv1U9TeqelpVS1R1kKreG/nDV9UdqvrnMY9doKoF4Z/f2l6BtiyJQCBgVRKxVzf8LonW+JxEW5aE7YxLW/jik5ROEt5xkvBOe5CELwQBThImOEl4x++S8I0gwEnCBCcJ7/hZEr4SBDhJmOAk4R2/SsJ3ggAnCROcJLzjR0n4UhDgJGGCk4R3/CYJ3woCnCRMcJLwjp8k4WtBgJOECU4S3vGLJHwvCHCSMMFJwjt+kES7EAQ4SZjgJOGdVJdEuxEEOEmY4CThnVSWRLsSBDhJmOAk4Z1UlUS7EwR8WxKmaUOpJAnTjMtUk4QJrS2JtpRxmYh2KQhoLondu3cbZVJCc0lUV1cbvdHjJWGaMxgriY8++sioFjSXxI4dO4zrxUri2LFjRjmS8ZK4cOGCtT2JDz74wHOdCLGS2LZtm3G9WEkcOnTIuiTkRrUVvxpjxozRq73xIilBgFHcHPCtvQeb9Tp16mSU7tPU1MT58+etje38+fPNRNiWt13Hjh3JyMjwXEtVaWxsjN7u3LmzUTLV119/3Uz6bXnbPfPMM9dM9RKRnao65lp17WVTfYfk5+dTWVkJhP479uvXz6heeXl5dLpPnz7GL1ak3sWLFykuLjZK+WlqamL37t1A6E0wevRoK2ODUPLVwIEDrdXr2bMnN998s5V6f/rTnxgyZIiRYFWVXbt2AdDY2Gh1212+fNkovi6+Xvfu3Y0zMyP1bEbOea4kIoXAqzGzbgX+RlV/GbPMFEJRdJ+GZ61U1Z96fc4IEyZMiAriyy+/ZMaMGfTv399zvczMTDZt2kTHjh05ffo0M2bMMMq4rK+vj0bfX7hwwTjjMiIIgG7dunH33Xd7rlVcXBzNpDx37hxDhgwxyrjMzs5m3bp1pKWlRbedScZlY2NjNPq+oaGBuXPnGmVcHjp0KHrIkpGRYZRxOWHChGgm5ddff82tt97KsGHDPI+tb9++0UzKyLYzybhsamri2LFj1hKtweAchKpWqGqxqhYDtwPngZaaLmyMLGdDDrE89thjdOnShUWLFnHy5MlrP+AqBAIBnnjiCRoaGqxc3cjPz7caXzdx4kRGjBjBBx98wIYNG4zGBvDwww9bzbi02ZwnNzeX2bNnU1lZyZIlS4zj68aMGWPt6gbAjBkz6NevH6+//jr79+83qgWhbWcz49Imtk5SlgCfqOpxS/WSomvXrjz11FPWJJGXl2dVErYzLiPNeWxIwnbGZWt08LIlCduXQNPT06PNeWxIwvYlUJvYEsTjwNIE900QkT0i8raIeN8fS0B7kkRsBy8bkrD9OYn2JInYDl5+loSxIEQkHZgJvNbC3eXAAFUdCfwr8MZV6jwtIjtEZMf1Xh92kvCOk4STxNWwsQdxP1CuqrXxd6hqg6qeC0+/BXQUkZyWiqjqC6o6RlXHeDmb6yThHScJJ4lE2BBEKQkOL0Skl4RPqYrIuPDznbbwnC3iJOEdJwkniZYwEoSIdAamAitj5v2FiPxF+OYc4I8isgf4FfC4tvIns5wkvOMk4SQRj5EgVLVRVbNV9auYec+r6vPh6f+nqsNUdaSqjlfVLaYDTgYnCe84SThJxOLb72I4SXjHScJJIoJvBQFOEiY4SThJgM8FAU4SJjhJOEn4XhDgJGGCk0T7lkS7EAQ4SZjgJNF+JdFuBAFOEiY4SbRPSbQrQYCThAntXRImpKok2p0goPUlYTvj0uS/V7wk2lrGZXuWRFvLuGyJlEyUisSwff755wSDQc91Zs6cycKFCzl+PPQt9S+++MJzrS5dujB9+nRWr14drWVSr6CggPHjx0dzCw8dOkRxcbHnepMnT+bMmTMcOXIEgOrqarKysjzXmzFjBi+//DLHjh0DQiE5nTp18lQrIyODGTNm8Npr33zfz2Tb9evXj8mTJ0fzN/fs2cP48eM91xs/fjxnzpzh8OHDAJw4cYK+fft6rjd9+nQWL14clevnn3/ueU9RRJg5cyZLlixpFk9oi5TPpHQ4HM1p95mUd911F5s3bwZC/y3GjLnmel6VN9745lvot99+u1F8XXy9e++91zjjMrbegw8+aBRft3PnzughVU5ODhMnTrQ2tqKiIm699VZr9SZPnmwUXxdf74EHHjCKrzt27Bh79+4FQiG4U6dOtTa2wsJChgwZYqWeyTrGk5KCyMvLY/PmzXTp0oWqqiqKiooYO3as53r19fVs2rSJrKws9u3bx8iRI40kUV5eTlVVFR06dGDXrl2UlZUZZVyuXr2ajIwMLl68yOHDh40yLrOysliwYAGZmZnU19fT0NDApEmTPI+tsbGRdevWkZWVxf79+ykqKjIKwj148CAVFRUEg0F2795NWVmZkSTeeecdmpqauHTpEgcOHKC0tNTzH1Dfvn3Zu3cvwWCQxsZGamtrjTIur1y5wtq1a+nWrRuHDx+mqKjIKAg30jbAJil9kvKRRx5h8ODBvPXWW2zfvt2oViAQsH7i0mYHr7Fjx1q7ugEwa9Ysa815gOiJy6VLlxqfuMzNzbV+4tJWcx7A6iVQgHnz5llrzmOblBZEWloajz76qDVJtMbVjbba5k9EWq2Dlw1JtMZHSJmGAAANl0lEQVTVjbYqifT0dKsdvGyS0oIAJ4m22ubP75Jo671AbZHyggAnCSeJ5HGSuD5shNZWisg+EdktIt+6NikhfiUiR0Vkr4iYtTdKgJOEk0SyOEkkj609iHvCjXFaut54PzAo/PM08G+WnvNbOEk4SSSLk0RyfBeHGA8BL2uIbcDNItK7tZ7MScJJIlmcJK6NDUEo8K6I7BSRp1u4vy8Q+1dVFZ7XDJO+GPE4SThJJIuTxNWxIYiJqjqa0KHED0TEU2dZ074Y8ThJOEkki5NEYowFoarV4d91hJr3jotbpBqI/Vhiv/C8VsdJwkkiWZwkWsa4L4aIdI1MA9OAP8YttgZ4Mnw1YzzwlarWmDzv9eAk4SSRLE4S38Z0DyIX2BRujPMH4Heq+vu45jlvAceAo8C/A//Z8DmvGycJJ4lkcZJojmnjnGPhpjgjww1y/iE8P7Z5jqrqD1R1oKqOUNXE3+NuRZwknCSSxUniG3zxScpkcZJwkkgWJ4kQ7UoQ4CThJJE8rS0JE74rSbQ7QcC3JXG19KpkaG1JmH6dOFYSK1asMKqVapK4cOGCUb14SZiQipJIycCYSIbk1q1bGTRokOc6hYWFHD58mIsXLwKh7EITRo0axfr164FQeIdpveHDh1NeXg7Axo0byc7O9lwrIyODvn37Ul0dusK8Zs0aSkpKPNfLz89n79690W23Y8cOOnfu7LleUVER7777LgC1tbXG266oqIiPP/4YgO3btxtlSALccsst0TCWlStX0tTU5LlWnz59SE9Pj267LVu20KdPH8/1hg0bxokTJ7h06RKXLl3yXKclXCalw+Ez2n0m5UMPPRRNj87Ly2PmzJkEAt6PlhYvXszp06cBmDBhglF8HcCvfvWr6PTs2bONMy4j9QKBgHF83ZEjR3j77bcB6Nmzp1F8HcCbb74Z/c86atQoo/g6aL7tpk+fTkFBgbV68+fPN4qvq6mpiSZvd+vWjblz5xrlP3700UfRPaUhQ4Zw7733eo6vA3jppZdoaGhwmZSZmZkADB48mMOHD7Nx40YjSQwZMoRNmzYxePBgtm7dSvfu3Y0kkZeXR2NjI6rK7373O5544gkjSYgIAwYMoLq6mjVr1hhJonfv0PfkItvuvffeM5LEwIEDOXbsGIWFhezatYvu3bsbSaKwsJCamhqCwSDvvPMO3bt3N5JEMBgkJyeH+vr66LbzKonICd5BgwZFRWsiiby8PPbs2UNhYSEHDx7k5ptvNsq4vPXWW10mZSx33303U6ZMYc+ePaxZs8bouDAQCFjNuLR94rJ///5WO3i1RsbliBEjrJy4DAaDPPnkk+Tk5Fjp4NWrVy+rVzdGjBhhrTkP2M+4tElKCwJC0ei2JNGhQ4c2LYlUafNn6+qGTUm05Q5e0HYlkfKCACcJE5wkvNOWGwbbwheCACcJE5wkvON3SfhGEOAkYYKThHf8LAlfCQKcJExwkvCOXyXhO0GAk4QJThLe8aMkfCkIcJIwwUnCO36ThG8FAU4SJjhJeMdPkvAsCBHpLyIfiMgBEdkvIn/ZwjJTROSrcFOd3SLyN2bDvX6cJLzjJOEdv0jCZA/iMvAjVR0KjCeUaN1S7/KN4aY6xar6U4Pn84yThHecJLzjB0l4FoSq1qhqeXj6LHCQFvpdtBWcJLzjJOGdVJeElXMQIpIPjAI+buHuCSKyR0TeFpFhV6lhrXFOIpwkvOMk4Z1UloSN5r1dgNeBv1LVhri7y4EBqjoS+FfgjUR1bDfOSYSThHecJLyTqpIw7YvRkZAcFqvqyvj7VbVBVc+Fp98COopIjslz2qA9S+LcuXNG9ZwkvNPakmgNPOdBSOhL678BDqrqPydYphdQq6oqIuMICem01+e0yeTJkwH48MMPAYyCOiKSeO2113jrrbeAUCSbVyKSeOmll1i0aJHxf4eIJBYtWsTy5cuNakFIEgBvv/02FRUVRrUikgB4//33AcjNzfVcLyKJl19+mWXLlhl/jT0iiVdeeYUlS5YY1YKQJABWrVpFZWWlUa2IJCAUvwhw0003GdWMxyQw5i5gPrBPRHaH5z0D5EGoNwYwB/hPInIZuAA8rhb2hbZs2QLAiy++yOjRo03LAaCqrF692iiZKhgMRqcrKytZu3at0Ziys7P54osvgFAmZWNjo1G9bt26UV9fD4RStGxtO4AlS5YYZWbGBtbU1tYab7usrCxqa2uBUCalqSi6d+8eDcBduXKl8R93hw4duHz5MgDPP/88Q4e2dAEwOWL/pBoaGlBVo394saR8JmVaWlqzP0wvxO52d+7c2WjjXr58ORpGCtClSxdrY2vr9YLBoFF8nao2k2BbXlfb9TIyMoyj4iL12n0mZWlpaTSCvGfPnsyfP99IEuvXr2fTpk0AFBQUGGdc/vrXv+azzz4DQqlXphmXESGmp6fz6KOPGsXXnTx5kgULFgChPRTTjMstW7ZEj38HDBjAnDlzjCSxePHi6LmDcePGGWdcRrZdWloas2bNYuDAgZ5r1dfX89xzzwGhvTGT+DqA8vLy6J5Snz59KC0tNZLE6tWrqaiosJpJmdIftZ4yZQp1dXW88sorxv0PAoGAtROX6enp9O/f39qJSxGhuLjY2olLgEmTJllrzgNQUlLCoUOHWLFihdHufFpaGj179rQaXzd8+HBrfTcA7rzzTmsnLgG+973vWWnOA1iVA6S4IAoKCnj88cetScLm1Y322sHLhiRExHrGpc3mPLYzLm128LJNSgsCnCRMaMuSaI1LoO2lzZ9NUl4Q4CRhgpOEd9qDJHwhCHCSMMFJwjt+l4RvBAFOEiY4SXjHz5LwlSDAScIEJwnv+FUSvhMEOEmY4CThHT9KwpeCACcJE5wkvOM3SfhWEOAkYYKThHf8JAlfCwKcJExwkvCOXyThe0GAk4QJThLe8YMk2oUgwEnCBCcJ76S6JNqNIMBJwgQnCe+ksiRMI+fuE5EKETkqIj9u4f4MEXk1fP/H4XDbG4qThHecJLyTqpIwaZyTBjwH3A8MBUpb6IvxZ8AZVS0AfgH8zOvz2aS1JWESwuMk4SSRLN+FJEz2IMYBR1X1mKpeApYBD8Ut8xCwMDy9AigRW1lYhsRKYufOnUb/+aG5JE6ePGlVEqapX7GSsBFuGiuJbdu2GdeLlcTRo0etSsJU/rGS+P3vf29UC5pLIhJSZEKsJPbv329dEp4j50RkDnCfqv55+PZ84A5V/WHMMn8ML1MVvv1JeJn6q9W+VuTcunXrormUphH5sT04unbtSqdOnazVy8nJMYqvu3LlSjSTEszX9ezZs83i8GxuO9v1srOzjVK9mpqaOH36m3xk07GdP3++WRxeW952P/nJT+jQ4ephcSkXOSciTwNPQ+g/1NXIysoCQsEdkWmv9OjRgwMHDgCh2C+TuLT4ej179jSqFakXSY42fRPFjs12vf79+xtF18XX69mzp3Hwao8ePTh06FB02pSjR49y6dIlK/Vi17Vbt25W65m+h2MxEUQ1EBuO2C88r6VlqkSkA9CNBLH3qvoC8AKE9iCu9sS33347t99+u8dhOxyOZDE5B7EdGCQit4hIOvA4sCZumTVAWXh6DvC+jdh7h8Px3eB5D0JVL4vID4F3gDRggaruF5GfAjtUdQ2hxjqviMhR4AtCEnE4HCmC0TmIcDu9t+Lm/U3M9EXgEZPncDgcN4529UlKh8NxfThBOByOhDhBOByOhDhBOByOhDhBOByOhLTJ7t4i8jlw/BqL5QBX/ch2iuCH9fDDOoA/1iPZdRigqtf8+GabFEQyiMiOZD5L3tbxw3r4YR3AH+thex3cIYbD4UiIE4TD4UhIKgvihRs9AEv4YT38sA7gj/Wwug4pew7C4XC0Pqm8B+FwOFqZlBTEtcJyUwURqRSRfSKyW0QSR2i1IURkgYjUhdPCIvOyRGSdiBwJ/+5+I8d4LRKsw9+KSHX4tdgtItNv5BiTQUT6i8gHInJARPaLyF+G51t7PVJOEEmG5aYS96hqcQpdXnsJuC9u3o+B9ao6CFgfvt2WeYlvrwPAL8KvRXH4m8ptncvAj1R1KDAe+EH4b8Ha65FygiC5sFxHK6GqGwhle8QSG068EJj1nQ7qOkmwDimHqtaoanl4+ixwEOiLxdcjFQXRF4jNbq8Kz0tFFHhXRHaGMzlTlVxVrQlPfwbk3sjBGPBDEdkbPgRp04dJ8YR7zowCPsbi65GKgvATE1V1NKHDpR+IyN03ekCmhCMFU/HS2L8BA4FioAb4+Y0dTvKISBfgdeCvVLUh9j7T1yMVBZFMWG5KoKrV4d91wCpCh0+pSK2I9AYI/667weO5blS1VlWvqGoT8O+kyGshIh0JyWGxqq4Mz7b2eqSiIJIJy23ziEhnEekamQamAX+8+qPaLLHhxGXA6hs4Fk9E/qDCzCYFXotwE6rfAAdV9Z9j7rL2eqTkB6XCl6B+yTdhuf9wg4d03YjIrYT2GiCUDbokFdZDRJYCUwh9a7AWeBZ4A1gO5BH6Fu6jqtpmTwImWIcphA4vFKgE/mPMcXybREQmAhuBfUCkNdwzhM5DWHk9UlIQDofjuyEVDzEcDsd3hBOEw+FIiBOEw+FIiBOEw+FIiBOEw+FIiBOEw+FIiBOEw+FIiBOEw+FIyP8HEMx3uv8tpFkAAAAASUVORK5CYII=\n", + "text/plain": [ + "
" + ] + }, + "metadata": { + "needs_background": "light" + }, + "output_type": "display_data" + } + ], + "source": [ + "#2 Create mesh and define function space\n", + "mesh = RectangleMesh(Point(0,0),Point(20,20),10, 10,'left')\n", + "V = FunctionSpace(mesh, 'Lagrange', 1) #Lagrange are triangular elements\n", + "plot(mesh)\n", + "plt.show()" + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "metadata": {}, + "outputs": [], + "source": [ + "#3 Defining boundary conditions (Dirichlet)\n", + "tol = 1E-14 # tolerance for coordinate comparisons\n", + "#at y=20\n", + "def Dirichlet_boundary1(x, on_boundary):\n", + " return on_boundary and abs(x[1] - 20) < tol\n", + "#at y=0\n", + "def Dirichlet_boundary0(x, on_boundary):\n", + " return on_boundary and abs(x[1] - 0) < tol\n", + "#at x=0\n", + "def Dirichlet_boundarx0(x, on_boundary):\n", + " return on_boundary and abs(x[0] - 0) < tol\n", + "#at x=20\n", + "def Dirichlet_boundarx1(x, on_boundary):\n", + " return on_boundary and abs(x[0] - 20) < tol\n", + "\n", + "bc0 = DirichletBC(V, Constant(0), Dirichlet_boundary0)\n", + "bc1 = DirichletBC(V, Constant(100), Dirichlet_boundary1) #100C\n", + "bc2 = DirichletBC(V, Constant(0), Dirichlet_boundarx0)\n", + "bc3 = DirichletBC(V, Constant(0), Dirichlet_boundarx1)\n", + "bcs = [bc0,bc1, bc2,bc3]" + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "metadata": {}, + "outputs": [ + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAQsAAAD8CAYAAABgtYFHAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJzsvWeMXFl75/c7N1XqwCabqZlzzuREcmbIya9mpJXWK0tYrCXbC8mAFrD9xd6FP1hYY421jbUBw8Ya8q68a0BaaaHVC+mdN5Cc4SSGGYYZhmEmh+QwNFOn6q5w4/GHe2/Vrap7K3Q3J5l/oNBVN5w+94TfeU64zxFSSp7qqZ7qqVpJ+a4j8FRP9VQ/DD2FxVM91VO1paeweKqneqq29BQWT/VUT9WWnsLiqZ7qqdrSU1g81VM9VVtqCQshxCIhxIdCiAtCiPNCiP8yOD5TCHFQCHE1+NuXcP/vBddcFUL83nQ/wFM91VN9OxKt1lkIIeYD86WUXwghuoFTwN8Bfh8YllL+cyHEPwb6pJT/bd29M4GTwE5ABvfukFKOTPuTPNVTPdUTVUvLQko5KKX8Ivg+DlwEFgC/Afzb4LJ/iw+Qer0JHJRSDgeAOAi8NR0Rf6qneqpvV1onFwshlgLbgM+BuVLKweDUfWBuzC0LgNuR33eCY3Fh/wHwBwC6ru/o7+/vJGpP9VRP1YEGBwcfSylnd3JP27AQQnQB/wH4r6SUeSFE5ZyUUgohprRuXEr5J8CfAAwMDMj/5Hf+IRePC6yuFE6XitWl4qbAyQrsHDg5cLIesstF0VyyOYu+bJH+TIGBzBgLjSGWGEP0KCUen9vAxUvrAdiw6hp7dn2BaBaZNvWnf/V3MC0DgDd2H2XF4jvTECr8yz//bQCE8Pidd/Yzo3t8ymHefTibv31/LwD9fSP85huH0FR3yuF+enIbX11ZBcCmNVfYveP0lMME+JO/+Lu4ngrAWy8dZtnCe9MSbpi2iuLyu+/8ip6uwpTD/GZwHj//8CUA5swa5jde/xBNmXrafvjZLi59vQyAhWsuoawbJO+muWPN4l6pl8elHCPFLMWCgeeoiAkVraigFUAvgFaUqCYYEy7ahIs2ZqLkJ1j3Wob/6X//Z7c6jU9bsBBC6Pig+DMp5V8Hhx8IIeZLKQeDcY2HMbfeBV6J/F4IfNTq/7mmR7pbYcVmjwtnHGxUbEVgaQLXFNg62ClwXIHrgpLxsNCQQkdTNXKaxoyUxrihYmiQ6fWHSJYvu875qyvpzhXYs+tLxBSJsWj+fR4OzSSbKfPB0Wfp7ipOCzAM3WJu/zCPhmfw8w/38Dvv/IoZPRNTCnNGyb9/1dJbXL25mPePPOcDQ5taoe7vGwVgzfIbnLu8mu5ckd07pw6MJQsGGR7twTAcDh55nt9684NpAYaquiyc94D7j/p578OX+N13pw6MGT0+zJcuucHNW8vYf2Qnb7x2AFX1phRud98jYBmLll3n9uW19KUF6tr7SNXBkR6mJyk5ULQEHgKhKmhCQQccCborUV2QtsQzPbwyKEWJOe5MKj4tYSF8E+JfAxellP9r5NTfAr8H/PPg79/E3L4f+B8jMyVvAP+k1f90yh7XP8yzYm8P67Z5fHXRxTMVXFMFJJ4l8EzwdIFnKXi2igOUDZ2CbTBmZxjRc/SqRQB06T/m7hcPk0qZfHZ6C6Z02bXzxJSA4eCiG2Xeeutn/PyXv8ZPD+zljdcOsGTxN5MPFADJjFkPeebZw7z3y3f48/fe5N2f/IyenslbGBPSAmD5qousXHKbX378Ij89sG9agAHw1ktHQTM5+sVWTOmwc8epKYVnSw89XebtN3/Be794h7/ev483X9/PooVThbGkb9YDduw8xnu/fIc/+9mbvPtrP6O7a/IwnvBsANavu0D/gtucPPoSv3j/DV7cNzVglANbfdWLRymgMXxmDRnPgBWlSYc5FbWzzuJF4B8A+4QQp4PPT/Ah8boQ4irwWvAbIcROIcS/ApBSDgP/A3Ai+PzT4FhLPb5qcv1wid5ZsHGDjWa7qKZvVilW8LFBsQSYPjBMS6NgGUzYBnk7zZibBaAk/a7CuJRsfv5TVqy5wJdntnP4xC5GXRjzJvexJbhAKmXxa2//nFkzhzjw/hvc+mZx+znQRP39Q7zz9nvYts7PfvEu+Xz3tIS7eNV5Xn7pY27cGeDfH3iFIcthzDMn9SlLv5Ual2Ve3vMxa1Zf4tSXOzl5ase0xDWdNnnnJ+8xo3eU/Qff5PadhdMS7uzZj3nn7fcomyn+5r13uZfvmnQ5KASVuuDByrUXWff8Me7dXsKHH7zFIzvLkJee1KcYNHJjZOh99gJdS+9SOrcM4/J3M57XzmzIYSmlkFJullJuDT6/kFIOSSlflVKuklK+FkJASnlSSvkPI/f/qZRyZfD5fzqJ3OOvbR8YMyWbV1sYlt8CqpZEtQjAIRCW/xiOpVK2I9aFk2PMzVLwUgCMeCmGZZplz51gxZoLXDi7nXOndjHVt/THPCjrFnve/DkzAmBcvLV40oVPAmYQpycFjLWrL/PySx9z585C9r//Jo6jTjlMIXhiwHjprffo7h3lVwff5NLthVNO2zEPjFmP2fvWe1hWikO/eJfCRNeU4zrkpVm45gpLnjvF4zuLOHnoVYbtLka9bMefstQBGHa6yMssxq5riMVDpC7Oof/r3JTj2qk6mg35LvT4axt0jRXPGWxdbnLyrgIouEatdeGaCh741oVuMKEH1oWWxfB8y2LMy6B5/iMve+4EZaly4ex2Smis2Nb5GIYlVRypMOSlmaWUMVIWr7z5cz7a/2sc/uANdr96gIFFk++SjAUWrD5ziFfeeo8Pf/UOf/Pzd9n39s/o6rBLMhGEVZTVcOevvMwzHhw//DI/P/gme17dj9phlyQ0lcc80INwt774MZaEU1/upCxh0/bOuyQO4Ebimkqb7A3S4NP332TPa/uZv2Dq40Ny5gTb3jjIF/vf4P1f/Do73tpPpsMxjLGgfOWlQT8w6mWZs/prAG59toOLH+5h7p4vER12SYpBIzfmZhGKJC+z5Lc9QHUM5lzrYcI2GZrZUZBT0g9iuffjWx5XT3nM6PHYsbCEbnkdWRehhp0uhoLPmMwy/9mzzF51nRtnN3P9y21TsjBC03FcV9j8+gd09Y3w6QdvcOnWio7NTxkzV9M3a4i9b72H4+gc+uW7TEyThbF89WWe2f0x9+8u5NMP3sSdJgvjmd0fs2zVJc6f3sm5L6bHwgiB0dM7yqcH3+LK7eUdpy1SUJJa9TfQ2z/E9jcPYJspTv3qTUoTk2u1C16KUc/v+g45XajLH9K/6yuK9+Zw55OdjFq+pdvuJ7QswAfGiJMj76S5uWGCx3NNlt9Ksfb+1NO1Xf0gYAHw+A5cPSvo6/bYvriMbsm2xy7Knp/oeTdTGccYcroYdrtY8twXFWCc/+IZRtz2zUQbFRel5hiAnrLY/sZBuvtGOPPhXh7djl1a0lTlSIEOP15fgW1vHPz/FTDq02DCEGx5432yvWOc/mAfQ3cHphzXUS+LnFli9eufYJkpTvzqLR7kZ7VdDiZkqhrfoDECv4LL5UPkdlzBHpzFyOHNHQHDDMptBRR2mjE7w4RjcG6lxZ1ZLhsHFTYNT32Auh1977shUT0aFHiGwpq1Lrvmlvh8NOPPhkRmRiTg6GrNzMgsT5AG8l4G4TYmbNfOK5SlzuC5dZQ8nb7NV9vqklhSw5W1vA2BgQ4rXjvM9fd3c+bDvWzZ+yGzF92dchr0zBr2zeYDr/P+L32zOdvdeiQ/74WDvHqlRY2qd+Ut1sujXDjyAofef5st+w6haq3N5kLYr/bSaF7jlNyKFz7HlCrnT++kJDVWbDvTMkwASyo4CW2ZkTbZ8eYBTu1/g9Mf7GPrq4eYtaCzadVKPkWU6x9hzeufcOngy1zYv4+B1z5Hy5VbhlV2/LDyboYshUqDBH5FZ9kYiryJ98VSzKNrKD1zB9TWZqzhGaSAUSdL3g1AYRsUgs9niwQ7HNg+5CFdh3Pqk63OPxjLItTDByoXr+nMzLk8M7NUtTBC6wIarIuy6yfiqJOtDHpGP3kvi779Bt0rv2H0/EpGzq7qqEsyFOnehK0KgJayWfHaYTJ9Y5z5cC83v1neVkslEZhSazgeqmfWMNvfOIjraJz61ZsUx6c+MAewYNU11r94lKF7A5w5tA/XmXrxEALWv3iUgZVX+frMVq5/uWXSYUXTomiorHr9U1K943z5wT5u3VnWZtpSY96HCvOuPMNm/r7jeJbOvfefxSk0grWZGkABfplbXKS85R7ag270zxeTL6fJ280/YbmtWBQhKCyDsq1j2hqfzta43i3YMeqwqWhNOm3b0Q/Ksgj14LGK1AXrl1g8pxQ5WsrimSLRuii7fgs5bqeRIrm17Nl+g5SnMXp+Jaank9l0s6mFYUsVD4UxN1tZ0wHUAAMV1rz2CZfff4lrH73AyleOMmPh5DqaNa1hX5lVr33K5YMvceJXb7HmjY9JdycPzBWkX+gLXgo9plUFmKEUWbDqGgAXjrzAmUP72rYwEuMZaOD501hS4+szWylLnQVbLzQNw0HFlUpsWKG0tMXa1z/m0sGXuXroRVbtO0LvwIO24liTR3Uqz3DofuUs+Y+2cOf95+jZewY1ZyZeb3vVtDUo1gAjOmY2tLBM1n1E31ez6Tk5wJ2tI8gmPNbDboidYcKtA4Wl4VgqnqPy6UwVbIcdJQtFFVwhlRzoFPSDhAXA/WENVxNsWmDygijyqZXFtQWKJfCIjF3oBmZgHuedNF6zVekayG13SQGli0sAWgIjVLSAADXwGFNTzN57ivuHdnH1oxeZt+cLsgseJYblASVPrynQs7TGrkZu1ihrXv+Eywdf4vKBl1sCo5XCiplbcY+l8iQ3j+7k5KHXWbX3CErCSH4paKXHvAyql9x3FgKWvnASgHtnNwC0BEYrDTldoMGcfScZ/OAZrh56kbkvf0F2/uOm95WkTj2C6vNPmzlBzytnyH+0hfyHW1oCA3xYlCLWRKi8XbVOBgccSu4YAxd7WXi6ryUwAB8SbgwobH9sSVoqh7MawjHZ5pioClx7AsD4wcICYHBMx9MEW+aW2a0U+cTM+taFIcFUKtaF6fqAKNgGLlVY9OqlSkb26OVKBs/YdhchNUoXl1CWOsqGu7HAcKWKlEpNQQsh0QAPo8i8fSe4f2gX9z/d3hIY9UpsCXudaQVGqNkrbwJw8+hOrn74YlNgtKvJAqOZFQCgpmzmv3qcwQ+e4cHH29sCBjTmUcO5Xg9lz2XcT9cw+uFW1JcvI7KNpr70qhUzCRJjdqby/c48B9OdYNmVLuZ/OYurG8djgdHt+jAo2AYFJwYUpn+TsBSELTimZBCU2ayaqGn4esJomQad6AcNC4B7EzqeDttmltnjFfnIyvrLwSPWheX4hbxo66R0iwm7jUTc/JiU1DEuDVD0dKx1j6if0Ux7KopUahIxycIYc7OgQvbl8zgfbeL+p9vpfvE8xkDcglaBJfWG7k2Syr0O8/adZPDQLi7u38v81z5H76pdEmy61UE4t0XlC6Usfcxs9xyPPt/EhUMvMfelLxqAEa4FGHa6UER7o/Ldz1yiLHXund1A0TOYuflawzWWp+JJpSUoQk1oOrlXzpH/aAv3P95B956vMOY1uk2RCExPjwVF3DExs4i65zLOJ2uxP1pHcfdNZLZ2IFdxM+SACTcNFBIhES13N+dILLfMmutplp3t5at1pQZg2MGLdAXLoOwmg0KxfItaOHDCTqMI2JAxUbolN8zpszB+8LAAuFv0TeFtM8u84hY5ZGYxDSrWheP6CVu0DRy7g4DXjzAbSF3x3+SNAwbUtiYAfVq1ZY+zMHpeOUf+o02MH9nQBBjx99eEFQFJamae+ftOMHhoF4PvPxsLjMmoe4U/g/Po8008+GR7LDA6lRAw+9lzAIx+5b+xGgeMODVLDyXlVLoO459uTARGu+FV8rUHlBdvkT2yhOzhpbHACBUHiigkCpHvV2cJTNdm802ddeeznFrl4EWAYQWWRdnxBzOjoAjXFYWgUGx/zZGwBKesNEoa1s2wUFy4k5+eav6jgAX4wPA02NFTZq9V5AMzgx1YF47rv4BWtnXsTgeM148ww1PJXZmN5Wrk14xUgKFJFRVB3k7To1en2FrCQwX2XINP1zB+ZAPK89dQ5o9VrpESTE+j1dKghoLe69D9yjnyH23mbjgw1+XHy3H9Qlz0UjjBfe1YLTD9wAjjbey8TkqqjH61ClPqZDdW35p2qA4ed6IQGCMfbmP8040oL1xFmZevucaUeiUNoqrPt6hGu6Dw7CD9n88n/ekyHj9/DzfjW1KGY5ADiq4BlBKtiSgoCoFrgyt9KrbrseO2yrYrCseWyQownMBQMy0Vx24NCsUC1fKPn55IocyQrJlloZYl9/JTd8rwo4EFwJ2ygdRgR7bMvkKJQ2SwDAXXFYDEslTEJJ54bE2BJUD3df/l2SgwQkVbFKAlPIThogb9Ye/YSqgDBrRvVUSl9U3Q88pZ8h9trg7MdcWvFeioIi4dIeddpnBiDXc/3kX37q8QqqxZ8CaUzgAiBOR2XQGgdH4pQA0w2lHsM2igvnQJ95O1eEdXQQwwQjUDRH2e2jMsHgfA6D82UAMMgIJjUAgs11aQACjbftrd6AfHETw7CM9eF3yySOApAseVgKzMerQHiuDdKVNw5kEa4ZZZucBGNVUGRxMftS39qGABcNsy8FTYlS6zr1TifSONRAE8n85Wc8LmDKuSubnI+MatNQXmeyp91/soezqPVo7T7akIqdS0JL26b/q3BQ8FeP423ccWJQIjSU1BEgOM6VB6uT8tWTixhvHDG+ne/dWUw2wXGJ1aGCLlJgKj7OlYCZCoz7eaONgZyMHEjhGWnJrJzGMLuLlzCMtRCV1OdQIJ8N9lArg2EzxX8PxDj9234KP5Cp7rAdK3KCylMj4BtAAFqMGbiBe/1lEWeSxb7qKUYfB4+2lYrx8dLAC+cQ2kA89Q5rWxMnezPo09W8VR/cTWDLeSUSnDqcnAJF1dWWappzLnRhe2p+BKBVfWwicKDmgDHobH+PO3yR5dgjy2koln7pAJBjjNmAId7dIkaczNQo+HsudKZSRf2ex7Nyx4RmV6r11FrZh6YGj9yXDrpILL7Xf8GajzS31rJWamqZUarAQVePE22cNLkEdXMfHcbTIx97UERJ3KvTa3dgyz+NRMFp/o585yP33Kjg7YHUEC/PeZAK70gnQVXhhyeOWex5ARlC1TQdidg0IzPRRTcuW8irLaY8l6BVFKAZN7CepHCQuAW9LAE/CcW2b+eNVUDOemJ+crCM4v8+3MgVt+IS7mnIbZlS69OjDSFjwEjD/zgLnH55I5vhDhJU+8NzOb60Ei+orVrs5nK5s8VXM1VNgl4yjyJvbJpdj3/dcex7wswp38OIYQoOy4iQfIi8G7NH3xy9ibpUGDDJfi7ltkDy8h89miStp2CohQYV5PZKG0Jc+aMz2svFB9RycOFK0gAdVyebkbsAUv5G0WBePT4dRoBRZtg8JDNT3Ussu14x5iq2Txjs5WpEb1vYRFqldj5b4eUIJKo6pIVYDqIdVgxFj1+3VSpTLlJJXgowKKP1AYrXa/cdfigRGGWaWrUKLfq3eoihH5Xk2q6CrQbEFjyVc9qE3663rNqtEZteeUKrZsXZIOCrNxfRZOyai7tnllTJr7UFIualBIvaOrkYuGmobTvqpWlfvhOkTPNHtwGumi9NnqhsPtFnc7Al075ZHO+79TV/vpK6Tqrq3mb/1krR2Z05wTOe54CpYOWpCFa8720Dvbxo3839rvAn/JHUhPEDZZ/ncX3Pgu8htjZcqBBSs8/4MAoYFQQBjBcRcULyjLrkTx8AHuSoQrQbr45tbk9L2EBUD/qjSlfGARCOGXSwEIDykECIkMjlV6AsF14XGAvFDokX4GdXkSYXp+nkTzJfG7iHytZrqCoJCS5Ez//NxHOsV0bUVWmqwUrfVtXNv9kYpEeH642XvdONn25nqVJsvY6+XdnonMTcN7BDkTJax0Y1lcU4cpTqsCkLVQij4o9Tu9eC1WTnoJSyDrV9NIJIJq2pqR6c9oFZJ1XcsoVjwZX6FDzXmsUTIkXqQg1YQnI18qx2VwXFZ+5hVBT1Dx57suBXw/m+ElQgSRViTCHwelUqwkCBmEKWXwAaRCedQmPaN1lztO31tYmOMOl38+RtnUIJ3CSxvItIqb1nBTwRhESsFJKbgpgZsC1xB4BrgG/t8gl1cpJrtcv8DlVYX3e9K4QiCNSMFOVb8rerXbohnV7ymjWrjSusWuKxpzR/1i9rDf5sYSK3YdRk5PrphdMefWH5jvP5+QjK0fojy3vSnOekUHVZXRNLmPlgMgMzbFF75B5jpZdBIv/eos0uf9XSC8HpPSs7fbeqOylTJHFqM98tt4e2Eea92jpl2HdjRmZ6ppq0jurx2j0F9N/2aL9QoJ5wqWwYxxwe6L/vmSAR+vhKLhp0Grbod/UXWBVSjFEmwyTbZZfvxGUflUZILuSHz3QzP9MlztfjiIsoNStqBsQtlEFkus/62OdgCoxqnVBUKIPxVCPBRCfBU59pcRf5w3hRCx7pyDc+eC6062GylrwkFRBevenUG6uxpFUXZQyw6qmdx61SRkXYN0PJVinuPyWr6MKmVN5oQZBrUZGc3gaMaXbR0J5DMet2a7LL2TYuHXGQpmY6EKXymO04RtNBRSV/MYHShidtv0nZqLe28GY3am8mlX0TcYC07gzWn1MDgq6cNLKYxNz9uqAOVN91Ef5sh8vijRnG43rnk7jYOC1WNSWDhO6vJs5Pn5bYcVTav6dPMUyciCIlbOYdHpmSgPcrF5ANV8i8u7gmXUDF4CXFrgoDvw0lVB1ooHhWdX3+nAVBJBEdVxLc0CHPbIEsokByejarULYZLasSz+DfB/AP9v5J/9x+F3IcS/AJrN9+2VUrZeqB+RZ0su/NUD1v9Hc1n3ZpaL+4uUy+Cl4yucP6jjr6VwjeSCelMamGnB7rI/S/J+bxo3yCRpeH7GBRaGZ6sVC8Ox1IqFYVpaxcJwpe/K79xSv4VeNegn5+WFBgh/GjaqaKGrtzbigHFrxzBLTs2svHA0MdunXxwwwoHTVrJ7ImsFPpvP4+cGcRNWI7ajLlcjDQwPlMiIx8w4249+bAlDOx9M3cIQktHN/vszPVf9NS7jq6srMjsBZ1Su7nFhyzhrTvew+lw3VzaNk5/p52ES1IEGONTrUVrh7grJnmsKe64IPlgqGYtUsThrApJBIQLoXlV0EPCMLPOiWuIYGWJN2Cesdhz2fgLErkcOtgn4beDfTXO8KD62ufBXD3wL460c6a5q4oTWhWJ6FdMrTqF1Ee3O35IGh9Np5rkur435FgZEMixC+2grkGRhgD9ldmKB4NZsl1WDGmvuqCCrrU9cIWtmbUjpvxeQR+fWjmHMbpuFp/voepS8zr9Za1ovu9cHhuIo9H82H7U4Pb3R4qJxRjc/JvU4w6yTcydlYTRIwOjmR4wOFOm52odxcXbHFlZoOYRAtjwVR5dc3pqnlHVZdbYb42GmqQURl4dlW6ds61huNf3u6xofLBWkXHj1hiRn++WrU1DUfLfhmjA46aYZUByeT02PhdGppurdZA/wQEp5NeG8BA4IIU4F2xN2pOJjmws/G0FRYN1enaxmI8rJraBq+n24sCuSpJbAgNhuiWOpFWiYlobnRSqDgBMLBF/PkjXACNUKGkngyKNzYfMEpS6Xhaf70AY78w8ZVqwJxwdNwfG9hz3OqtzcMfydAyMJcrb017GM2RnGnAz3NowxOlBk9tfdzL4W332KQqEeELH/G50vNpQoZj02XcjQN1J9casVIOLW5diBO8LhjODAApWUB2/ecciU4stVO6AIpZpw3TM4ZaYZ0ByezX37wJgqLH6X5lbFbinlduBt4I+EEC8lXSiE+AMhxEkhxEnTrQ7MlYZdHxgqrNtnkI7Ulah1Ec4vxymcnVQi43mTBQbUDVRJUbU0BHyxSFaAseIbo2H9S3vWRm1hcXXJpS15Sl0uK7/qRhvMdVQpklTucbi5wx/DmHlsgOJYV2LlTfqUIw5aws/gPIfBDWOkHmfoPTFAvpx8f9sS1ACj90rvpJ/f9pQKnB0dTm8sVoCRibHemgHCtDRsu3Y60rFUhtIKBxaoGC68fd8i58iOQRFaxGHDp1jwtWPwZSHNfMPhmRnfLjAmDQshhAb8FvCXSddIKe8Gfx8CPwWeaXLtn0gpd0opd6bU2lHv0rDLhV8VUBTJhhcFGdUf6ExS1LpQ6iyMJwEM8AuNaWk1wFj3QLDiG4OyFT9VlQQN8LshUaujHhgzHjeGORmAlHscvtkxhOIqLDk5E704dWe9AKMLSgxuGCM3ZLDo9EzafHs9URO2wYRjcGV1iUfzyiy4mWXBjeawiaZf+Ikb2xuTBkdXO0xkJLuu6PSP+RU2CRAQye86RcvFI1XnwDwDw5W8PWiRC1wltNv1iCosx6oFNyyD0/k081IOO/u/PWBMxbJ4DbgkpYzdvEEIkRNCdIff8bcunPTLBKVRLwAGbNgtSGdrEyi0LupnQOJUD4yjWhUYenB/O8CQngjX2FQvtTRMW6sBxoZBQdlKbp2aWRuVa2yDfGA2T+S8RGDUK4RG0fGvLbtaA1C+K2AkdRsmbAM3WE7fAD0BN9YWaoARB4VmA5Whomlu6/DZGpvxNOy6otM7FG9FJEECwLUbx7oepxQOzExheJK3h02iXgNqlm9HVmeG5VM0GXe+WfKBMTfjsn1+iSnuS96W2pk6/XfAMWCNEOKOEOI/D079DnVdECHEgBDiF8HPucBhIcQZ4Djwcynlr6YS2dKox4UPLRQFNu7yyCl2pSuSJNWqJroamU6NZsoNtQqMfaVSDTCSBj6jio5lhDJtjWNzFa72UQFG2AA0a7EKloEEbDc+axwNzmwsVoCRfpBpu3IkacI2eJxRuLQlj+IoLD4xCzufblqZw49Z8eakx56/M8flxtoCuSGDgS9mUShPvtsUquAYfLXCZnCOzYKbWZbeai+sChhctQFyGjGlAAAgAElEQVTMZVtnHJ1PVkrG0/DC14K5+eaACPM9hATEDGSaCo8NlQM9GVKe5K18ie6yrAFFqGgjFm30olZF9PftMZ2zwynm5Fy2Lio/cWC0Mxvyu1LK+VJKXUq5UEr5r4Pjvy+l/L/qrr0npfxJ8P1rKeWW4LNBSvnPpiPCpTHJhUMBMJ6RpDO1CRRaF2FXJE41GdECGBBvZXi2GizTraoBGkJwfL6IBQY07wtDrcURLdxRYGy8mGHWcMT92iRa2FDFbpdLW/OormDtlz2kSlMd0vL1eL7JjbUFekZ0Vp/r7rhLEvs8Ai6tKjM4x2bZ7VQDMOrTLslqq09/S4NPVkrGUoIXvlaYN1FbvsI8rm8cADw3eSBzSFPZHwDjjWKRnOe1BYrKsTpQqJYfrzvDBuceppjT7bJ5uflEgfGD2woAfGCcPyJRFNiyxSanuC2nUQFEBB7TAQzwl/I2tTTqgLH2joJpNha0ZtAIFS34Y54RC4zY+2yDUtgNcbSmQPk2gLHibC+lcjzYKmM0UsGTojnsAmDc7ndZdjvFgq8zLbtzoZoNWI57Gh8sEYyl4JVvJHNGRDIg7MDXRE0g8QOZI57GwWyWlJS8USiSC15DSAJFCNUkUITX3nukc/5eiv5e74kC4wcJC4BiHs4f9nxgbLXJ1L2bEbUulEhLpjwBYEDdyrxAlVbI1irA2PgYtjyUmKYaa9qWbR0kuJ7aEh5jnsGx1Q75jGTDxQy5B5NfDh2tqI/SKl9uLKG4Cmu+7MUdTyVW6tD1W9FpXvkLtsGtfsmllWX6RlU2XczU5EtHcY1aDLbBmWUOt/tdVt/TWH0nHppRODhu4zX1XQ1LExwY0BjTYd+gy0ChWr7CvK7Pbz+g5jMeQ6rKB3qOFJLXrAJdVjXcuK5HK4WzgPfvq1z6Rqe/12PjGvuJAOMHCwvwgfHVcYEQsG29RU44LadRYRqB4YmGqdGkQuTYGkdnaTXAQMqWg2bRrkrs3L4Gn6+1Gc9IdlzTyD1It2V+t9JEl8fpjUVUV7DtXJZ0eXpWDN6f53QMjJbPI6gBxvJvjJbpBsnjESHkTVVwYGEVGPPyjeNVANhB2tgRb9tNZjyGlAgwZIGs9FqCIsmqiOruY53L1zX6+zzWb3SmHRg/aFgAFCcE5wNgbF1vkUk1dkUqA5wRc68ZMBQ7HhhhptcXhujgZ6jY1kcIjs7SuNyj1ACjEkwTaISKqwT1wJgzUhuXcL8J8FebxvXn4yrhkwLG9ZkKZ5c59I2qrDuf87skdXFxPeF3Q1rArpIWjs7nCwU3Z0rW3xesH0yOq+uJpgOWUZmq4FdzU4xqglcf2AyUInSLyfd6SCTNeIy4KodEDgPJa16BbDCtFh3MDEGaBIqo7wrwXyC791Dj6mWVWbMka7dJxDTW8B88LMAHxtkzGooC21eZdOPUDHRGFV3ZqUTWYfhv7kXOxQBDlZK4RTMVxRQeqLM2hOCzOUoFGJsGBfWT/xKB44qO4DEudT5eQSIw2lEcQB4YGsfW2CiOYMvZHDJfez5c6pwEoCQo3Z7tcXaZw+y8YOcVnXbddza1GAScXCxrgBG1Hvy0bIRIHCSisDdVwf55BmN6AIy8jM1nYSesmUgYyByzVD52fGDsVQt02X4ixK0+bgcU/nGPwUGVa+cFM+fAmmfEtAHjRwELgGJB4expDUVItq8yyRqNpS9K7XpohOoYGFLUWhlQhUbSuEYEGJtHPDYNCpyYwlcJrqHAx8vW4OMVMJaGHdc0+h43N8PbVT4n+WytjebCc5cMsq33Cm5L7QKjVXciKtPWODpP4foMWH9fsOlhMoWaQaIhXEdl/8w0Y5rg1WGTAdNv9n0vVrV51w4ownI2LFU+KecwhOTldIEupxrfpFmjZqAA/92pB7cF1097zJwnWP2iPi3A+NHAAnxgnL5goAjJrqUluqXrux9LGuCcJmBATNckVBI0HK0GGNuGPBxTiR1xbwiyCTxsDT5dKRlLw3M3BKH/X9Pxr7NdtaF1bqcyTicwov/v6gyVU4sls/OC7ZcNLNM/Hs6GtAJEXFpIIfh8QHB9Bmx+RAMwPLcxnZsOWAb5ZyqC/bPS5FWF14ZNFkzENUjtg8K/Bka8ABhIXuoqkFG8unKaPE5Rr3Blsyg7PLgJX5+w6VugsnpvdsrA+FHBAqBQVPjiahoFHxhZPdnCgKkBA9k4wNkRNByNo306l7vUCjDCLon0lKZz+jVB11WYJGC0o2YgeWhofLJSorrw7EUDdUKvbODU7L5WQLo5C04tlswd9xdDNeuStGtlSSE40q9xrVuw+RFsGGxcog/tQSKUsBQsx18vMaYqvDpeZsByEE4wFuHEj0+0WmwFkC+rHJnIoQvJSzkfGFG10/2I08Prrg+MRTqr3+xFTGFx7o8OFgCFcgAMIXl2QZGuYHBBtas1ezLAuOVVgfFqoYQWVOxo4QgVQiOxi1K5MBj0jAIjRlFwtAOPCU/jgyUKoynBczeUjoDRTKNZH0Sa5zt5ybXYWqFd1QMj9M7XLhygMY2kEBydq3KtW7B12GPLUK3j5k4gEc1HUxE+MJQAGE51XXYza6KyetiKL38jrsrRkRyGItndVyCHV7nX/9u6+wFU3sxWyn7AD88X+fpoib4lKda8OxuhTi7Pvrdu9aaqQlnhxM0Mu5aWWDqjmmuqVXWQo1i++z3/uO+OLzwO/rkwo0IXfbc8AzR4wSkjXCiIRrPTM2rNjbCg1bjxCwtkyqsAA2DziF+gZYtXu+uBEXX/F8pSBYeWwr6bkg33/f9nOypmiwoedR8YpxAYe64JVjzuvOAlVfrL3WAPSJ67JxFAvo1Z31bgDIEBLluHA2e5cWkbN2AZZyHi57ON4GA2y+vFIlsD13fCEZXmtx1rIq6RGkXlyEiOF/sKLOlqdHvYKSgo+//84dlxsG2Wv9wT+0ztSEzWxdaT1MDAgPxP/37wCkpYGRv+hlf7X6pOe2vPpfTa5yvHFZS6QzLpXOR7JnKVA1hx26wnKWH+WwCZCE+KSmPcOlEYJcOtdVpVmqYmIlPHlCcRrq2AE6mzky2uQtalbciYuPBaOOWNKiVljbPfUrPCFPkukq4JlK7zMmZWZlpkzZ+KA18Zd042/DW6/Ez64z/+41NSyp2N/zlZ31/LQghGrhdBC6KoBiVG9bNGBn/DJ5DBeS80sZSglVdh0cxq6ZuwFQpB6WvYniOS6zXnouuwIteswq5E4Ruh4RJsQxAjqcSUiJhrhZSsDubysx5c7mrScsaFmaA1+WpNkcDd7uRrO9GqyL7DoymYmPz7YbHh6h7cygncdmDsNb9GAKsngrR14XK2rvgnzD6ImHArMxUCVnlVC8AFBoNCWbPgLAKqmuEIN/64cGFxbzXcfFnBCte5BF6/FTdwqRDu1+KE9waBVv6G5/0L5m6eXOZ/b2GhpRVGbpQYHQwOpFM1f710YNqn/crkpv1HiXr+BnBSChOmxbr5FlJCTvM4O5KuACPqs9OrK+hu5HflnFftkqQ9ySIcPCAjPT7Ss5VCXd8VSToGdd0TYFm5SGgQlVX4YoZWsy1BkqJeyet1pVfh3dt+YUm5cCujcjdXS8u4rkwrjRuS7Q8kroAeC44PCCaa+EGNqlkXImM6LCxKJJBx4dB8BU9pDDd23CFJpsKKiWKF0UVFcLrbH3CN63LEeayqHI9cft3Tecv2PbBnkNwWOo8iXa12xsf8c9XyoZpQcgRrZlk4HuRSHpdvpjBtpY0uSECNoAsS/pVF//34VPfkhiq/lwOc5pjNxGCZ1e/MZkbo1LnuwcM+WTM3e/X6/G4GRcDzs4vktHABTPygZzNF+6AjKHwm0syVLq/YxYoDnbiCFjcQCo0FVQIXshqXsxpbxly2P3ahjdWTTd9ZCPTpXJURQ7B30GVBoRZS9QOEzT71+mCJQPPgtZuSdEGZVBj1epwSHJ2jMlCU7Bt0KxvotPOcFdWteZHAuZzGlYzKtgmbraNO56CoPxYUwSMiQx6Fl2WRuc2cUcSoHhRRnbyZQVclO1abpGJm974tfS9hIT248Bd3KDywWP3ObPqWtrZtQ7qqdQNA0TdRxy2V4x0Ao1krENVNYSQCIwka9aofcZcCjvboPjAKDtvHbR8YCatE65VUoSwFDi5IBka7Ciu8F4wBPVJVDgxoqK7vd7LLbr+L1EzXepUKMPbe9RBtPHvSorhQEjjSa3A1pbGtZLGlWJu5nYA+OpBpIfhAZMlLhT2qD4y2Zt2agAIgX1I4eTODFgIjwUJ90vpewgLANT0u/vUDCg8sVr3e2xYwWkm15LQBoz5Tk4Dh3zdJK0OIRmCE4baoEFF5tooMul3SVbBUMS3AqNdwWnBwgYbmTQ0Ynq2CFEjpf7+S1TkyS2NB2WPfQxvViwm3g/QQrkDYKkdyqRpgdAL3+vUToVxL4SMnS54AGOrkt1mI7GxJcUzw5bUUmirZtsH6ToDxvYUFgGvJWmDMCzvy8V2RZn45Q6nm9AGj/t2QVsBotyAC4IpgKbGaDIxQHVQU8Ctg2dPYP8dgWBfsvecyMFZrjbRt5tepGTDiwm/3f17t1hqB0QkgwnUvNQdFDTA2m7UtQDvWRNxxCx8Y457CiykfGJ1aFXFvTo8XVR8YWgCM1A/Lu/cTVwUYjxwfGItaj8nWd0UUp9GbVhIwEsNM6pLI2kyOAmOvVQsM/972rYxQwlY5lk5xOdUEGKE6qECWKjgw12DYEOx7aLOw2DjA2apiS1dpuO6xqrF/ruED47ZDttjoIGgyutqtcaTX8IHxwG5I23olLoyLSLEVjulprmkaW60qMDqxJuKc6zqWwsflKjDmaJO3MKJd6cKI4MwFA02TbN5qf6vA+N7DAnxgXPrFqA+MV7JtAaMdxQFjMgOeULtCLwTGHJKB0TE0hOBYLgKMMaf1ooM2wNEOMCaj4ZTiA0P6rvDDNyonpchzXM0GwDA99o2YDWnbDiCgLq2F4Gi6Coz6MQyIh0T9m8p+uNXvFgqHx7OMuwrPd/nAmKxVAdXGb7ygcO6Mhq7D5q02Rvp74t07Ya/TPxZC3I3sd/qThHvfEkJcFkJcE0L846lEtAYYe7P0zfYLdVJXpN66qFeYIZ0Ao50BzyRg6FZjhibCQSa8cRgFRslm+5iDMGPeeo1T+GakLWohYvrvO3ynwKiLD6bir0uQIhZ09cDQgjRoKx2IXzeBEHymZLim6GxxTTY5gYWRYE3EDUTWT4CoFlhS4fBEFRizg9Wx9a4TOtXEuMLZABibnpHfCjDaSd1/A7wVc/x/k1JuDT6/qD8phFCB/xN/g6H1wO8KIdZPJbIVYDx2fWAs6Mwwqn/ZZqrAEF6MTwuqLU4UGC9LHxixZmyTaboGoNQBY1vJgmCT504qTL0sR+VAX5phXfGBMSbjK3G7Cq4fRmP/zDSaB2/ft+kqMPkwAwlL4ZpmcCSX8oExXm7aJYlNx5jzCMFnWroCjM1W/Guebb39GbUwpMKR4SwTjsJzM6rAiAsvalWEi66gttELy/HEuML54wItBMbktn5tWy1zqtlepy30DHAt8PJtAX8B/MYkwqmRa0kuHSz4wHhR7xgY9U59pwqM+nO1YTcCQ5XJwGg2blFT4BOAEWqy4LAUwYGZKR8YIyYLyzEWRn1FD960bAaAYV1h/6w0Wrh3htN5lyTuma6ldY7kUgzYbgMwWgEiek3tPxIcl2muo7NZmmyU1Zoc1+0AEHX5GVc+bKlweCQCjFRoAbd48BaayAvOn/CBsXG3eKLAmMqYxT8SQpwNuil9MecXALcjv+8Ex2KVtH1hnFybGmDMmlW7gq2+KyKc9ky0VoOerYARBw3VhNtWPDA6sTLqr1Fshc/1dCIwQlUqmd06XGgTGJNQx8AINnBqBb0aYIyV0cw20tBtYs3ZgBB8LmqBEVep4/K8WUNiS4VjD7NM2Aq7ZpcqwKjcG7Eqoo2aYjVaFVAt5xN5wfkjHpoRACMb+2hT1mRh8S+BFcBWYBD4F1ONSLPtC+NUAcaw9C2Mee3/ryTrAhqB0eN1VlmSrIzblsFxpxYYkDwNJ7zmraJ/keBzPc0VXW8KjBo5Sk0rHfexHfWJAeNAT8YHxpBJd4nEOLSj0DL4WjE4lk4z4LrsLZVQYtKgpZVRD+8AGF97PjA2KPXTqvFhNLtGtSS2Jzj2KABGf4n+TNiwxYOiXRVG4fwRiWbA+n3GEwHGpGAhpXwgpXSllB7wfxO/h+ldYFHk98Lg2LTJteHSxxaFYcmaZwR989pf/j1ZYDSss4ipm0lWxi2ZDIzEuftW5rQQfJZKVYCxfcJGaad1bSHbUTnYlWFE9YGxaMJrWqmF3RpCwlIY1vzduTTp787V5XZWMZLS45quNwCjra5IkwFM1RKccH1gbFR9YCTlbWM868KLdGltT3D8TpaCrbBjfhUYrRRnVUC1zBdG4eIhC80QPjBy0+NrJNSkYCGEmB/5+ZvE72F6AlglhFgmhDDwtzv828n8v2aKAwY0dkWg1pyD6QEGNJqe0evqr60AQ9YCA5KBUQ0vofBHgLHZsnz/Cm1WlmayFMGBHh8Ye8fLLLRiCnULvxtx6hQY7T7HNV3nM80Hxr5CvIVRCbMJoGu7HLXAWK/XWhiqRY3XeD++deHVO402fWB8frcKjFk5P22jZTLafU7yhFWvwojk4ocBMN7qmlZgtDN1GrfX6f8shDgnhDgL7AX+6+Dayl6nUkoH+EfAfuAi8O+llOenLeYRVYAx6nsznjlneqaRxi2VY4/aB0YzaER1SxocdycHjGqYorYSJQAjVDhmoTiiI4i0BYxJqB4YPeU6uHm+j4a2x3CCz3XV4DMtzfxgFW0cMJId4cYPOCqW4FQpzQ1bZ4NRBUZcfrcDilC2Jzh5I0PBVNi2uFwBRivFWRVQXUZQGJFc3F9ASwXA6Jqe5VTtzIY07HUqpfwHUspNUsrNUspfl1IOBtdW9joNfv9CSrlaSrliuvY6TZJrw4Wj0gfGVo+Zc2Ts8u/6dRfNrAuA4kQtMHJqtaTFTpsmQKPeyogC4xW3ceGWcKutXzsACQc9j4t0IjAS72vxcRxl0sBoFu6op/F+Jose7P/Z5bXfJWkGvHpgaFbygDI0g0Q0zwQnrSow6scwoPZdDmgOCv+3xHYFJ29VgTGzO1g/FDNVCq1BAUDZpDDkVoHx633TAowfxArOduU6UWDIioVRv9x7KsDYPaMOGAl1sZ2uSQiM2SIeGDX3tQsOITgu0lxV2gdGO3IchfczWUYUHxiLii6KJRBu6KIwHgitNKyqHMxm0WVrYHQS7g3P4HPhA+MlGW9hJEHC/19xRwVfTqS5aeqsz5isTSfPe7YDilC2K/jySopiWbB5hVkBxlRVGHK5+LORKjC6p7bk/kcFC6gCYyLvA2PWrPjCN53ASJw2baNrUg8Mo8XWi9AGOITgc60KjG0lC6XNqdNmsoTg/awPjFdKJRY409QlaQKMTgBRny5fiwAYNAKjGSSa56Xgi2IVGOt0M2blZvugAL/sOa7gi6vpCjD6ehuB0a5VEf1eeORw8b1RNEOw4e/NnRIwfnSwgAAYJwUTeVi3wWFW//S8gj1u1wGjyTqMqFp1TaLA2KMW0YLl3u0s2EkERwQYm1yLrY6/dkXYrbsezVQPjIWTBEZDl8TVeF/PoUt4o1Ck25Px/jHbff5AUWC87PpesuL2V+0M+ILTo2lulXTWdZmszUUWbk0CFKEcV3D6vEGxJNi01mbmjGmyMB45XPgPD9DSKhv+3lxSPZPbdOp76VZPNRQyWQORqkZPZCJETEUYF3HhJlPUfL9xUbBhp6Svz8+gHt2lrEX56OEZtb/d6G9Z61oPz3d/eG4kxc7+csVpb1fEygh9Knpx+eEkHLdhWFe56BlsCErmSsXiuqdDpLDVu/1rUKRQe0HSXVINeqTHgPQjNk+6jDc4H61TGxsInVTTvOoVmRf4eZzheVgwJQvGBU5oaXY7pUrXKW6dS/3YQDMpFgyjckkxWBek7TLF5kaQERXIxERbsYltThUbUOFawaBL9Zib8uM4W3UwIzVKsYBIfqsWYER/e5CKXu+BBldv6GxZbzGjx0+DbFb6lkQuiKrpQlf43QY9iLxlVeuGaUHgclKWyrgW3PhwmFVv9ycnVgt9b717/+Ef/uF3HY2neqofrX5U3r3NvM2tT0crv4URQbIRwbVefQQv8h29aoms2VoF4sOHCkOPBZ5W22RIvbZpcbXa317deU+DHf3VJvjSmMF4gmtvmZDKXsLx59KlyvfjThqnxX4Asg2rcom0WRi4fx5HcEZLt2Plt9QepxrXG4rOHWXqRUpx4UVZDfe0SFFo0WOufz8jSS9okbQ103gRt//NLJYkl5qKI1ncZTM77VsXE7bC1aFaE1Cp8xim1g24i7rzSrAUft36qlX14A6MPhZgR8bK7OieCZHvVjUxpGVFvvvH1/zmQPzDtND3EhbSg1SPTm6myjfHxoOjRUQ2fEvGrHr7xo58t/DSYUa5yMDjt34Rlq+TFIvQ3+/x6KHGw6G6vS5TtYXRqfvtpupgkoIFRYe5GQfLEyzJ2Rx7pJNX4na5SuhCBHnq1p3bIUt4CAwkSxSbT90sbjNgWNU4JSmPwkLpUETQjWSG43JWpBK9hsd2l2KUlR47XJMxobDYs7mpaNxVm9/cekYHlkmbeTjYwCpp8b7IUYjZrLOTF7EUC1y1hI0ghWSxanPEzOIhEgeiw/viFI5PlFyF2ekiZVfQpXt0GR5Xh6uZ0WycAhoH21XTz5N0GpYtdymOw+z58HgQxiJroJVy5L5yJFEjg5yhR28Ar+Q3bnO3FhKetLm+lwOcZt5m8OQIC56fyZKXZkxbuF+d1ZmYEKxb7zCnu7ap6HR2JCykeVtJXLhVE36TZcJxBfWWo3PcTDNHuLykFFHbsAPCqcBmFeiESHMVnY1YbJZm4rRqdOAw7lOvD/Qso0LhJafEQsvu6N4kjaBwSOTQkLwmC+RkdTu/Vs9ZeQ6rMe1v2jonrTRzVZcXU/G+RqL3xinOH8WZ4TR38hqrZlqsmhks3OoYFI2D8V+dEBQLsG67pHd2fHy+DX0vYYGEr/c/ZPDkCAM7e6cNGI4D585o0waMcOOY6CxJCIwk5yadAOMb1+C4mWa24gNDt2Rb7yX4cU2uTO0Co5nCih8uTHNswSFyjKKwR5YYkB1QoYlGhFoFhlegp81lz63e4bjpGBVgPN9VRKmDcbN7k/JW2HD2YSMwpirHhvMhMJ4T3xkwvp+wCFQPjKhJVT+fHCo651z/UplqOriuqAFGf9/UpqfCRVn1wKj3hxHVZIGxO+1bGO2+zASRVjh6vRDTAox62UJwSEwvMFQT8pbKR7YPjL1adcPgOLWbNqoFt4sGXxTTzNFqgTEZUFQlaoCxYna1bHZqVSiR1/gd2/dbUZr47oDxvYYF+MC4fzo/rRZGFBgbVts1wOjYuqBaIb9NYIRhdAIO8K2B8I3KL+w0193vHhjRbkXFIvKo2fJvFJWPnGRgdJIO0WX6t6xaYGgJ+aVaMjkvGx5RcP52irsjGivn2DXAmIrccZfzh2UVGPO+3er7vYcFwI1DIxVgLH4uN6WwwlVw0w0M/3gyMOIKWrMCLtzaBUFJwGgnrGQJTnpVYGx1TFRTxlfeDlQPjIWWHRtmp2HXA6Pb9jqyJJIGMaPAeHZGY5ekmTURF3+/fAi+ulcFxqqZbSxeaUOOTQUYa/bo3yowfhCwgAgwtuaqwJhEVySqZsCo11SBAcmFrlVhDwt6K2CEYXUGjiow1qsWmxSTuKWTSRU9nG5UrdrjniX42M4x5im8qJaY3+7cZgvlLZVPyjl0JC+nC+Ti3uaLxrsJJEIpFtyeMPgyn2aOUQuMpDxLAl1tuagCY/l8h2XzqxFp1QWpfwkyWo4rwMjLbxUYPxhYQAIwJqFoRiQBI8kreFvhmzFvq04RGP69cLdkcLLYHBjRMNuDRnvA6FQ2go/cHGNy6sCIPsuop/JxOYcukoHRDiTCcEN9U64Fhm7Fl4HOLC3BpRs694bUBmBMRd64xcUPrSowBp78KogfFCwgAMZXxe8EGPXWhYjbRi+iyQCjncp924oAw2g9rapY1UVFyfv1fn+AoViB53SZnCZxwAgB0S4k4sKNAmNXf6mxS9IEFEpMvvplRnDxllEBxvK6MYzWVkW8xetYVIHxavaJA+MHBwuAm4cnqsDYGkmghK5IdNVbkiYLDLzkLgk0B0ZiX7iNOloBhuYDw+hgH4poV6X2I/iinOa6/eSBkRyH9sOtAAPJK6kCWaU9azDJ+Q34+XI3r3NmJM3sdC0wkt9UjU+j2rJSBcayRS5LF07ByoqU7QowxrwnDowfJCygDhg7Wzv4jVM9wUNgjBdaj2E0hNUCGElewxPn7N3WlScKjOe7fGC006o2l+ALqwqMzdJEsWRHFToOAK4l+KQUGcOo90XXoUILYrys8ulEDk1IXupKBkY78Y7mxe1CLTCSFm4l5Xu8013B5Ssagw+rwGhlVbQjJ29ycX+hCoyFU99EPE4/WFhABBgbU5MGRr1cV/DV6VpgtLQuAjUDRrNtBlrN3TerpPXAUJEdmePxqgJjnWGxUW+0MML4VLo3bYDERvBJ2QfG86nOgZH0XGNuMjDahURcHkSBsWNeCUXUd0k6AUVopQouXa8CY8nSzuBQYzFH5JiyCoy3ep8IMCa7feH/IoS4FOwb8lMhROwCCCHEzcBX52khxMnpjHioBmAkdUUSFEdy1xWcuWhMGhhJhWgqwAgVB444YFTDnSw0WgNjMuoUGO2Crx4YXa7X5qBx82nRe8M65x6m6c/WAmNyoAjlA+P+oMKSpV7HwKhRpLxXgDHqPBFgTHb7woPARinlZuAK8E+a3L832OKwo9dhO9HNwxPcv2i2tDCSBopiTb+ibABGJ5osMDrZAzMKjWbAACoeqBWntiSN8U8AACAASURBVAI2r4TfHjDq4yRk5/9qvKRyZMQHxu6+Alk1eQyjVVpHxyfujNcCI3GWxI4PL25mTTUlVy6rNcBoKIdObZlrp/Fzxspc/NloFRhLpsfihjZgEbd9oZTyQOC9G+Az/D1BvlPd/LxcBcaW5EGedvcVARqAMSdb2woqHogmqx4nA4xq2J1D4+5Ec2AkxtNq9hGcmZgaMOLC9SzBkfEcedcHxjx9cmMY9ZbWmNMaGE0HOBPWT0SBsW1RuaFLkhi/WFCEx0QNMBat7HC6vhw/4uqYsgKMtb8xZ9qAMR1jFv8Z8MuEcxI4IIQ4JYT4g2aBdLJ9YZIqwFinVYDRDo0heWCpvktSDwxovoNUUrekHWB0ammAD4xT+c6B0Vw+MG6YPjA2KWYQNyou6poBJ0m2FBye8IHxbK59YLQaaI0DRjtp2WxaVDUlg481zt9LMaurFhhqwu7wzUERygfGgzuweCUVYNQ3am2V4/ouyc9GKQ3b0waMKcFCCPHfAQ7wZwmX7JZSbsffSf2PhBAvJYXV6faFSYoDRkO8O7AuVNNrHMPobby/1ZZzkwUGVKHRLjjulKvAeCE3fcD4sugDY23GYn16mrokbQKj06nVMUfl2MMsGpLdM1p0SVosO4/m3d1RvQYYhpXgkqAtUATHyy7XvhINwGiqBKsiKseUXPirBxVgzFg2tT0NJw0LIcTvA+8Af18m+OaTUt4N/j4Efkr8NofTrpufl7l/xakAo12vkM2mraLA2LTMSgRGKyujXu0CoxJGm+AIgdGv+8DQ29wcurnigDF1xQKjxaKsJEXTJ2+rHHuURROS5+cUG4DRDiTi8iwKjM3LTdS6AtYZKMJyJGqBsTbiwatN67hesljCKXsVYKz77QVTAsakVnAIId4C/hvgZSllMeGaHKBIKceD728A/3TSMe1QN7/wM2FgnYbnNma4KDsVT1qtpJoebkqpAGPLOotNyywUBcaLjSjSTK/B01Y1LD8uUc9bITCeWVDi+dlF9DYRHgWGazTG407ZHw3f0VPmuRl+Nil268rX3DGwDwyAtRmLmKRtS/VxcBEctXK80Ffg2VwJVcCI3V5CNANnCIznZxd5fk6xUrFbLdluNg2umR4PHqgojsG6xRa9XZEZuEmBIpTg+kkX4QoWrRV0sOdSU0sjBMa6vzubdb+9AEWbnI3QsrYE2xe+AvQLIe4A/z3+7EcKOCh8t2yfSSn/CyHEAPCvgl3J5gI/Dc5rwJ9LKX/VTqRSvTrzllVnY4VeG80af5yVJ6l7FF2llJdIKVGCErJgtcAsRguBg9Tq3eDZeLG11cUL/HI+HlbpDRzndGclC/rjTed6P5+N52t/35/QWNLrh7W822bC6TxTZX0zBwzbKrMM31x+rq/EmfzUrYwxVwXsSuVbnTYxpdJkOXmdMvGHB02dvmB7wD7dY1kmnmwiSqk2ZggHSzpLuvy0XdVjYZbj7c2KH86EBlhxPOiu/h4vKvTk/Fq9daXJFbWxStVvcgWgxIxziGD2Y2JMMgeBovhxnLdapcYBm+3Q8NB1ZcX3vdlVc2zo0gRdLz/17v1UT/VUHehH5d27+Njkwl/exbM8RLrWE63IxAyA1l1DqkreRZs05q70H/XhN5Jb52XN2JxMNzrZdVPxSRPdV2T3rqrpd/GWwaOx5N2e4roJtf+vev61Zb5DVU/C53czFCLmeL1z33bk6oJZhsuzM3xPY3lH4bORLNMxjLGx22Rxxm+xb5Z0Lk6kpsOHDj+ZM1H5fno4zYNSZztpJc3ChGnrevD53SxFRzTtcvhhNT8/K+2wYbWfBuMFwblLBp6Hvy9IfVhmvOkVtwZo+UbJrMX+c9+/6nDnK8ffD6RedV0QWWqcTZSRa3b80fLkh2mi7yUsnLJHZqbBqnfnceEv7vjuvqMqRT19h6rL0HLVA3hpzD/3+JbLnCUqtokPjPBOJaYiWy5uzJiGG7n20bBCV1Zi2bBmkYXlGDweS0hSp9FjeI2KVWDYLoyZKt0pj23zynx2N0sxBEbgWbCZJ+8GmWB6fth3ixoDGYdtvSU+G2nhNbwN5QPz93ZJY2nGxvQEFydSxO7a04EelFT+v/bOLEiO47zzv+xjbmAuYDAYYHCQAOYASIAUSJEiTZMSJVGySJmm6JVseeWwvZJ35Qg71g/r9YOtkMPh3Yddb9i7YYe89tq7lmRblkSRFiXxEGWKIkUSPEACMzhnBiCGmPvsu6sq96G6Zqq7q7rr6pnuYf8jJqaPqqzsrMxfZX755ZetEQ1FCm7qTJFWm5lJla+uhi3CzpVC1WAhFaa9UeWW3iSvjjeTVO3vSyStoZX4LaG0RioXIX16LkRPt8bgjVnePh9Fs9hwSbPZ8U1YMCC+pNG9L8zshErv4QhKGq6dsrBNpPLrvkwWQ0pLrpfI6mSy6Hsnqsq1IUpS5cJj19m+t5nhT+8lFHVQ8RxMJU28lmXqgsKew4L9R9fTtJtKtZodKTRWKSplZ0kMuZleXc2EeGWyGSHgjj0JWqLWlny3Uawm41HeWNBnSe5sT9CQ1VxPzVrp9Goz44koA60ZhtqcTauar1v4B5DVBD+daWE1G+LkjiQ9TfZl66YsVtJhXplsISwktx1I0hy1vi/l7lehMXNqJsy5yxE62zVuOpwhFHJWnuWm8sdeyTIzprD3WIS9J9w8JYJVVcICYG50dQ0Yg5/c4QwYZSTRZ0muj0nHwLBSITAKp1WtHLcMuQNGuCQw1s9xB43JhA6M7kaV23ckCZu8Ecs1YDtJ4PRqUxEwvKZnKCtLA8MLMAHiqyFOXWkmHLIGhltQGJqajejA6JQMH1PygOFmRWnedKmEsVcUHRgnmjYNGFULC8gHxtDDPXnAyIv0bSeb3saVV9KWwLCS0xvsxNPTUDl/DFh3SXYKDHDXcEoBwzZ9iwYvVPN3cGa+kSsxHRhHm4Nz3CoEhldIAIicsWY1FbYEhldQGJq9FuLC+bAlMIry4gIgY68ozFzMbBowqhoWsA6MbX2NRcDwo/G3ZBEw/AxHwB0wwFkvI5yWJFZLD0mKz3PWkLwAo7wEby3qwDi8PcNge4YggRFL54DR4m6lpq3bfQEw2mTpdMuBwqgX01Nhx8BwrFSasZ8kNw0YVQ8LcAEMB3YLQ6FUxhUwrGQVVi9oYBhKrIY4Nd7kGBjgzK5RC8BYCwKcErz8rg6MW3Y7A0apkAGGDGBEQpJbD6dpanC+1iPv+4JVp2ZgHB3OOgZGOY/NzQJGTcACdGBcfHIuDxh+hiKGrIBhJbvhiFXbqhQwYukwp8abCCG5o88ZMAwZe1sY+4aYVW3ACKfXY3AW5lXRnAHDCSTMSi4LXr/YRCRsDQyvAZwNYHR0w+CtMg8YrlZAFygPGCf9bY/hVDUDC4D5C4kiYPiRQfBCYLgZjoDzIUmpCucGGK9ONBMSOjC2SW87qhWG9J9ajHJ6OgeMrsoDw88+IqWA4RYSZvtRLBmyBEb5HkXp7+cmJJfOCEtgFKVl16uweOitAeNk64YAo6ZgARbAiBQAw8VQxCynwLBTOK2VnCWxi7hllhdg3HYgyTapumogdno3lgNGsw6MaMb7hkPrxwvOTjVydVkHxlBrJpC8FgKjN5r1DAmzCoHRQul6YLfuY+37XD2amRSOgeFGY88uMXMuuSHAqDlYQD4wBj/eXgyMQin2N9xM8uIhSfENLTc7EgQwnPSXCoHR0qC5fqpayQyMk7uL405CPjjWNhkqCRTBmVkdGIe6MhzpCsboKZNwaryZWCrEif4UO9ucAT5UZvXbGjBCkluOZmhqtJuydgYKQ3nAOK4Rctr6yjwAx360uiHAqElYgAkYvVFnwHAoMzAO3KixGcAAffFRuZ6GFTD060tf4HACDPcKBhiFv03RBKcmmllNlweG054bQGIJ3hxpIBzCEhhuQWFoZlJw+Q2Njh4YeL/IA4brpegmiOQB4852d+k4VM3CAkoAw+VQpPAmOQFGOQUBDFjvLttVdDtgrOfDGzSqCRjl4FcOGE78WgyF0trafYklQpbA8AoKQzNX4fIb0hIYfmQAo//OjooAo6ZhATowLj27UtkexkBx9FgnzlpBAcOQHTjKAQPWF0OFss57HZsFjHBa6rMhWum4EmZZAcMrJMwqBEZrqQCelK8Xhi3MNzBsHohjP1pl5kysIsCoeVgAzF9O+waGVRdwDRgHsQSGcBD5pdDwaQcMt1NzhQ3BCTCK8yZL/kFlgBFOw+i1Bt5ZiHCoK8NgezpXTv5sLoomePNCI7GkDozu7c5misqVvRkYN5/I0tRkE8W7zM53hUZzMzAGb5MIdwtrbXX5qfmKAGNLwAIsgKFYLeUtPSb0AgwcRuEuBwzwNpdv7m14AUbpPOuNd3o+wtuTjXQ3q9y2Sw+Fb7+LfGkArZ8nGLneyDsLEW7YmeVwjz+jp1EGiip445IOjJtvSJcEhhtIJxfhrdMRwuHSwLCT3ezazFXdjbu9N8TA3dF8YCjup8UN36NKAGPLwAKsehjBpFsOGE7Xj1QKGIYiaY3UiggUGIauL0d5e7KRrlaVW/fpka3DabkWBSqc8dIr8AcMu2GZE2A4LWdzzzAeC1kCI2wTtNdQuWn42XHVHhhWcmiTMwOj/2e6HZ1TSlsKFlAAjPtbi4FRpqDtLNIbCQyksyGOnVIrgtcv6J6elQaGf7kDRjmDr6FSwPAyfDRUCIyWMrvClwOFUd9cA8NCVh7NBjD23bPDNzAcwcJmC8MuIcTTQoiLuf+dNud+LnfMRSHE53zl1qHWgNETzgEjYKPnGjDy5QYYRkW0BQbr3WQvvY14SvcVCCG5bX/lgBEO5HFjBYx1OQVEoQqBsbPZ2b4k5WY71oARgmO3S5ps4ok6BYWhPGB8IEzIChgenA4vfmeS6dPLvoHh9Fb/LcVbGP4e8KyU8jDwbO59noQQXegBft+Pvg3AH9pBJWiZgbHvfQVh+Dz2LiAHjCuw5yB07yqOBxVOWWxDZyM7YBTGHgZ342tDa8AQkuP9+m8OK94anllmYBzqcekbYCsdGJOzOjB2blMR0p1vRKFCaQ0tITl9NkosITg2kC0JN7vehJVScxpnXxWEw3DkuHX0+JJ5s6ljs+MqYz9J0t4XYe9x5/vnlFsndelfptaA0XHQm+OWo1G9lPJ5IcSBgo8/iR71G+DvgB8B/6ngmI8CT0spFwCEEE+jQ+frpa7X1Bll4J6+os+F3Z0O2/fZhCkM3i2PtLF83TS+LFVzrELtrSUqMTDRuh0GTlhVsIxlpG0rGWH9VFWsuQH396k0NpbpLpfKo0lZVdAQ1dM6fmOGmUXnhjOtRBEJ0+U/cEOClVTpZ0/ISTs0pbm9VeOmgy6fpDYGZ3PZ7t+j0lxgoLRaQWyZvYIhjKJANBcX9ehJydyUzB1T4t5oGhC1/k7VgDDCVLjDH28lE5egquhB9a3OswaAVNcLXfjsYPsxAe6SUl7PvZ5CD/1fqD3AO6b313KfFSm3veHnAdrb2+keaCMxl09fYdc4SkxSJxZVWjp1mEQaQ7TtANVYSlyy9GTJxphYhZZcSPgdvfr74upWOg3zcTJ3WFYBY+eDnm6NmMW+JHnpQ9lwl4UBdHe0qyTS/odmsaSgrVlPvKVRItBQAxjtxJOC1ly6PZ0qsWS5H+gsXUUFY+cHo2xdmV0sgKIVcLdzhyRtuZOOOY0Sv0fqddlcb7ftjJBqVNEUm/M0DbtBgizIc2I2TctOb8vaA5kvkFJKIfxZu6SUXwG+AvpWANmEyvlvvUtyfh0YIauo3jkVB/DNqamR3qEGDrxf/z4d0zj3TBxNWf/eTlpTqVDaksEPhOjq0d8tzsHEeYFVRbAK/GsntTHEfXfq0Zk1DcavRphbdG7t0myCAm9r1rh9SE83owjeGmskmfZvcNjXk+XwXt0WkMwITl9qRJP+QXTiUIru7Tp5ZpfDjL0bpbBsvdhyzGU7cSnMwoLDjYxSStH1DW3rkNx8R879PAOjP7UHRkmX7oLhcd/Njey7Va/zqVWN808uYrXIuNQQRLOI9H30l73tY+6ntkwLIXYD5P7PWBwzCfSb3u/NfVZSmZjeko99tp/mbg+x7200/lLSZPQsf3w5X32pSuIr8G7OhmHph4G72IvhtIaiwPWZsKXRs2yey9g2Lk1GCQl9RWWzzQIpLzp3NUpnm8bxQ+lAZkmkhJV4iMnZCAd7FW7oy2KUrVejL+i9/MnpMPGYYPiYQleXA+c1h/fvylmNUASO3S1otNikyOs2hGMv6jaMgY+2B+O45fH2+IHF44Axu/E54DsWx/wA+IgQojNn2PxI7rOS0rKSM1/VRy9mYFhRsqxMtJ4bz3Lp+QJg+DB2AiAl46PCETDcQENRsJ0lcSK7BhVPiTWjZ5DAmFqIMHKlIVhgAOfeia4B48aeDKG0t/gdhgSgpiVvvxVxBAw39yy+DCM/kSWBYasS9XD2YoaxH63S3t8QHDA8yOnU6deBl4ABIcQ1IcSvA/8F+LAQ4iJwf+49QoiTQoj/DZAzbP4R8Gru78uGsbOcknMZS2D41fyEe2CUk0ipjoABziufUCUkJGfejHgGBlhDwzxLUu3AAMGFS2Emp8Mc2KtysF/B68I+80yHooiywHAMChPA7IDhtVdhaPZ8atOB4QgWUsrPSCl3SymjUsq9Usq/llLOSyk/JKU8LKW834CAlPKUlPI3TOf+jZTyUO7v/7jJnBtgOAqxl5MlMErIyY2uBDBAt+KfeTNCbNU7MABELpqTyEpCaY3kMrx5tqHqgGHATQ+rJ3OgE1wYi3gGht10aClgOIa6xXG+ehhWSusPs80GRtV7cBYCo7HFQ6XOFt/QImAoPocj6MCYeEOrCDDefmsdGD3bvMduNCueDK0D41CKVqH4cgQzVA4Y5msU/tnLPTCc+E1YAcMPKAwVAaO1hNHXRc92M4FR9bCAYmA0dQaz6MMtMJyq0sAYGlbo2aY4diAqpTVghODEcIbmpvU03Tfodc1cDzF6STd6nrghTSTrbXVtvpwBw41zFRQDo3NnMCHv1oARhuEPNpQGhgttFjBqAhaQD4yjj+6yBIaboYghN0MSN+POjQBG9w7NdcOwUilgWMkMDmPDHiugTM+FGb0UpWO7xs2DzkPhl5Y9MPyUhaIIzr4C8VUYvEWWBYbTGK3J6Qyjz2UIRSoHjMFP9mwIMGoGFmAChhC2wPCiPGDcHdkUYLjZ+awQGOCvoYB7YDjVRgDjxr4sYR+zJEb5q4rg7KuiLDCcgsKoK4klaQ8MH8b1NWDsb9oQYNQULCAHjL+/6g0Yafsb4wYYbuQUGOBuq0QrYIA/aNQSMMJpyeXRENffDbFvv8aBgypuZ0msIF0OGG5BYagkMHxo9nyKy0/Nbwgwag4WoANj5BvTlsDwMhQxlAeMexpsgeF2GmziDa10AB2TggAGeIdGNQPD+E3rv0tw8ULYEzBKlbMdMLyCwpBrYJSISm9IJpLMno1vCDBqEhYAyYWsLTBKqky3bw0YO0SgwCgbccssCSEHiyzKAQNMDSzjvNFXEzCKAVEol8CQEFIclG0BMLra/YHCUB4w7ovS2BZQD2MDgFGzsAAfwCgjp8BwK1fAYL2bXPIp6AAYZukBe7Wiv0JtNDAK8yM0CY4vWR4Ybj1oYR0YiWXJwO2CTqulkuZcZJyln1iSjP4gpgPjgTZrYDiwZRT2oisNjJqGBVgDo+xQxMGNcAIM4WGJpVtgGCpV2d0CwzJ9C4CkluCtNyOEhH9gmNOdmxScPxemY7vG8SMZomUC3TqTNTC8QGItxZSCFlM5+6IksYIjYDhSKk1iQWP0qXhpYHhQMTCCSRe2ACxgk3sYDuMgmOUVGGDf2wgCGFZKxPWoUCEBtwxnaA2peQ3f2N0rlCmGTamey8x0mPPnwrR3SI7dpAQ2S2IGxsH+9cVnrlMyla+aJThgmB5UfoBR6oFoBsbQo3sCA8aWgAVYAKOjTB/M4ZSVE2CEUpnK2jBsVAiNSgNDCDh+IktTczBOS0EDQy8PlbG3JFPvQP+NsO+wu7IVKcXSkFkKGKGMs3B9VnXOEhg+1yrBOjA6bmgJDBhbBhaQD4zhT7SXB4ZDzU9kufRSdnONniVk7m2814Bh3dMSXD4rXAHDDhJmWQHD7wIxsADGtmCa5fSpeS59dyowYGwpWEABMB7sKA0MFwSfv5isamAYCqcUiKuMvMqWBYYTw68bYDidEoViYHTsdtiEytS1PGA81FkSGG7cA2ZOrwQGjC0HCzADAx0YnQH1MN7R8oARrhQwAhhiqorQgbFCRYHRXClgFNwyZ4AolAUwTGXrpDdhJbmaYfSHaRJLkiN3R+noC2bIm1jQGH1iiVBUMPxQJ03bvNdbc+yXQmB4jXa/JWEBOjBGn1gCAb1HS6wRdjM+TKXzgNHeGybaZF3wbu0YZmBYRff2IlURjJwSxFZg/4HgomKZgbG7L7h0zcDo6JA0RPzNZujKB0YoBEKRniAB6w8CNQujP8qQWJLsPVrihrmsX4l5ZQ0YPUPFoSK9Oh2agbG939uaeSELo7lWgfr6+uSv/cpvlD7ICRyFINqUz8NsymflzgX5jTauZ0BqEqWUjcsFyKMN6wdnM2WCuzpUJALCVAzZUgxzk9eCANVZh3Y+oORoK2oKW6JpoAawIl8IiJjym824rPc2h0cayIvEnTXvyuajbdnWWydplirbFr238qUvfek1KeVJN3kK6BkWvKQmmRtdtf1eRMp30UTuEd170zpJ43MKqUWL2ucgvTXlth7oPaKnL0KCxUllPQiwhaTj9CW7b8gBqUFwfUy/8063FbBT3/7115oK81YRUx1KmrZQ2LN3Hb6rq4JkLhq5Fx8Us4z8hkIwO5OLgu9RRvh+AfRalG3Z88vsOWrUA9Ajxy+9q7nLsE365nobm86SWnDWS5Fl8tt3m7ete6oWFtGWMCtXE8yfi1l+XyrSt1mipZnUssqBu7ehqZKWrjATL6ySWrIo0BKRvq2Oa2gRdO0NIzVJU1uIc89nSgBDKRMtfF07+yES1Su1psKVs+vGOekiUrhZM5OSEx/Q04g2wvK8YGHGf68lExMcHJSoKrS2SMbPQCoh8DvCbWrR6NqpP0gbm2H0NYGmucuv1VCj5wCEjH1asnB1tDQwnAwl566oHPtwrk40C5amNJbGHMaLLTFMySQ19t3ehpLRaOmKMP7MHOnl8t2scrFqW3Z4C1Hp+Y4KIQaEEG+a/laEEL9TcMy9Qohl0zF/4CTt9HKWlWtJBh7uo3uwzWsWizT6xBIAww+VmSUpJ9MNji9qjmZJwJ0d493LORvGYcH+o+uNxDDKeR1zXzit2zAGTki6eoIbgp59VRAytvNrCSbd1SW4+JagvQuG3icdT6uWLB8JkxclU+OSvQOCfUP2AHJrqL74YobEsuTIXVE69gb3HB79F92GcfTRXTS2b97z3TMspJTnpZQnpJQngPcBCeDbFof+2DhOSvllR2lrMPqP10oCw0uk78SiUhoYLo1RhgpnSYKI6Yk0GT0LgGHICzQUhTWjZ5DAiK/CmVeCB8bsdeEYGE7LQ0oYO10aGF78J5QMjH5/lcSiypH7WgIDRmJu3ei5mcAIajbkQ8BlKeWVgNJDzciywHCiQutxclENDhja+ri8IsCgPDDAPTTMsyRBAiMRE5sCDK89LTtgeHa0ymRRMzD6VNwZMFzUtcS8wsg/T28qMIKCxaex37/0TiHEaSHE94QQR+0SEEJ8XghxSghxKq3qvYaggFGossBwI5N1ejOBAe6GKFsBGH6GY4YKgRGER6YrYDiUTCRJzGY3FRi+YSGEaAAeAr5h8fXrwH4p5XHgz4HH7NKRUn5FSnlSSnmyMbxuvCwFDE+bDuVUEhhuffMrPSTBOTAMOXJfrkVgvKbR3iUZvkUrctzyKjMw+m/y0AAt1oZUAhjApgIjiJ7Fx4DXpZTThV9IKVeklLHc6yeBqBBih9sL+O1hSBsHgGoChhNouAUGlIdGLQDD3GOauwYXX5Ns3wFDd4hAgBFKZZh4Oc30JYU9RyPugFGinpQEhov6VTiU3ixgBAGLz2AzBBFC9Iqcx4oQ4vbc9ea9XKQmhiSp9FolcAsMcNbL8AIMAJELaCuyal7jE7lAtdUAjLx8qRKktadlUMAohPT4KcUdMBw0+K3Uw/AFCyFEK/Bh4Fumz35TCPGbubefAs4IIU4DfwZ8WvpwGbUChp+hiCFbYHhdKlylwLCTEeRl9AWN+JJk4IRGd6e/KVpDhcBoDitFsPJyHT/AEIpqW86OgeGibhQBIygYbzAwfMFCShmXUnZLKZdNn/2llPIvc6//p5TyqJTyuJTyDinli34zvOE9jE0ABhKEal+hIXhggO5WPfKiJL4ER24TdO3WP7dr3IUN3O6Y5JzK2Rc0QkJy9G5BU2sg2XUNDMf2Ibc9DAfKA8ZH2+nYF8zevRsJjJpcSFYIjK7DDhfGlOmFVAQYBYvP3Mb0NLrKVhV9I4FhJ8MVWqRKuxgnVuDsTyShEBsODE/BiUoBw2N9UDMw+vgiiXmlssDoiJY/yYNqEhaQD4wjP7fDOTDKKFAbhiGfwDBkBY5qAIZTbTQwvEDCLEtg+IlilUrr9fZflpwBw2kELvKBceyz/RUBRs3CArwBw8kSX0tg+A11lko7DqDjRGZw1IFRAIzbJZGSS2udKw8YQz4SMtUfV8BwocRslrNffYdwNFQRYNQ0LGAdGKvvpivfwwggNmKQwDAUSmW48kr6PQ+MhUsZLv80y/aeEAP3RAPzwxh/Ic70+TR7jjfRf4vDxYZmWdSbcsDwGrciPp3m7NcqA4yahwXkCv7bM46B4fRGVGRIQmWAAXDllTRTF5T3DDDMvStjuDF3RQsOGKap8PGXUt6AUdIPI9gehjEzWClgbAlYAGhZd8BwqiJgNAUQiYUcMCqwkdHE68oaMA4MBhPFaTruHwAAF9lJREFUqpqAUcrgaygQYFg0ctfAcOSHUZkhSSWAsWVgAe6A4aabVzFgVGjnMwMYfYMRDgxqwax32ERgOAFEoTwDw9SbsJJjYLjyw8gHRruHfUms/I2CBsaWggXUgWHIDIx9JyJre1uIjFLUfXfaCCsJjJEfZggJOHaXoCWSBU3fwtAP6FwBowwkzCoLDA+2LTMwBh7cScfB4vibXhQkMKo2UpYfGcAYeriHIz+3g9V3/RsmYR0YQw92MPxQB2pWj2TlV/MTekM+dE9zyajhbjXxug60vsEIrZ2lbRhOGqXW1LAGjOEPCI7cJogtlj7HaWNPpmDkuQzD9zUw/MEGVFWiBjChMXdFA7LceEc0B4yCcvBotB5/SX+S7zneRFuP6cHhwwiuZiQj35hi+JEeBh7cSXzGWQGU82I2gHH0l/o59tl+T7voQZXCItwYoqXNndU51FhMzKsvLDH4cA/b9+qrWNt6G8jEClu3gnAYos/QlRdjHL6/fe19S6f/DlpyWeXdM2n23KznZfdAhNlx/ySaGVNp6QzRvkt/rHb2hcgkvLobr8/7X30DBn+2ge3deuPb3qqgWAXBtYl+bqcrb2Q5dGcDRuTglg7/htrEsuT6OZW+Ib269w2Gmb+Y61U2e7930+cztHSFad+tp9uxW5Bd9t6kZDIFhLny/CJDj+xiW5/eBlp7GtAU+3umpcvfTylh4tkZDj/ovTtYtdG9v/CFL2x2Nuqqa8tqS0X3Ti9nGXvaXQjqUIP1eGzgoZ6113Pn4sxfiBcf1ODeCj3wQMfa66svx0jGg5muHPjgupXv0vMJVOH/NnXvC7Njv967SK1qXDmt+N38DICBn1kvt9lxlYVrQeyInp/ulTezpFZ9ZNYUosBcthefT5R8YjtSJsvOgWa6Duq9gMSCwjuvxt1tA5CxHm6Y6+3026ssjRfb2DQXXp6Ghj61x/U5UKWwkKqksT1K+74Wxl0Awy7i9/hzCxy8r4vEXIbuIy3MX0ywcDFRcFQS0eLOqLQwnqZjXwNKWqP3WDMjTyyRSvkvUiUj0RRJtEnQc6SBc8/E16OGO41AXqBMQrJjf5h0XNK0LcS27hBXT/s30l55I8v+W6IkljS694dYuKayOOl/ynbxXZX2nhCqoofaH/1hhlTMZcO2sB9oqiSbkjQ0C3oONXD+2bh3u1Mu/WxSo+tgI+lVlZauCG07I7zzisUDyUKljOxXX1hk392dxGcz7BxqY+FSsggYXlZdL407y1uhqnI2JL2qMPnyAn23d3Lwwz3lT3CocrMkXrzmkosKI48vrW+V2KQE4+k5npsl6QkzeH/r+iyJYbX3eI2J17L6LMlQhH3HA4xA/a8Z4guSw3dF6dwTTLVKLEtGnssQCsHQBxtoanPYcytTPnNjGS6/kGR7b5iBD7UG54fx/CrTI0n23NpK/+0BuaUCI/88TWIuE+gsiRdVJSyQMPHMrGtglKNspRy3UktqPjCCcg2fsAHG2oXdQ0NimlYNEBhqFs5VABhJN8BwUR5zY1nvwLC5hpSVAYaa1hj55symA6M6YZGTF2CUUzlgePXJtwWGT2iUBQZ46m1sGWD46Gl5AoaD6zgFhpu6Vg3AqGpYgHtgOBnDbSgwYGOAsZYJ542npoERzQbSe3MFDBfX24o9jKqHBWyBHgb47mW4AsZaZspfs2aAkUqTnE4x8v0YoZBk6IFWmrYFU33LAsPjvSsFDK/1qxAYnYeCs42UUxBbAUwIId7ObU94yuJ7IYT4MyHEJSHEW0KIW71cZ0sAAzYeGMY1jem5TLao8lcdMFJpfQMnTSvKa3JJY+T7cUJhNgYYPnsvle5hDD7St2HACKpncV9ue0IrJ4+PAYdzf58H/sLrRZwCw8100nsGGJaZXB+yTLwYZ2o0rQPjmAjE3mILDPNQyeqvjDYMGEowywQKgeG1Tpm1BozZjQPGRgxDPgn8X6nrp0CHEMKzz+lm9DC8qiwwPDbGQIFh0sTLKR0YxxrZd7IpP592f4pS8jh1Nc25H6wSn1N1YPQE47S1IcB4oCOwss0Dxl0d5U9woOxSgjNfe2fDgBFECUvgKSHEa0KIz1t8vwd4x/T+Wu6zPFltX2injQaGnydBSWBAbQDDp9QsnHs6rgPj3hY69wWT2YoAI5VmbiTG5edW2N4XDR4Yb62y9/3tvoFh9J7VlLZhwAgCFndLKW9FH258UQhxj5dE7LYvtFM5YHjxbNtUYHiARh0YAQHDYgg0dzEdKDBkIolMJBl7ZiEwYBjaKGD4hoWUcjL3fwb4NnB7wSGTQL/p/d7cZ761pXoYUAeGR3kGRjlPz4CAUVhv/ALD6kG4EcDwvSOZEGKb8Rr4CHCm4LDHgX+bmxW5A1iWUl73c12zSgHD625lmw4MlyuB68BwCQxFde7p6RMYdvXFKzBK1elKA8Nvz2IX8EJue8JXgO9KKb9fsIXhk8AYcAn4K+A/+LxmkbaM0dMsVXU1Q1AHhgNgpNKeVtp6BUa5B0utDUn8bl84ltua8Hhum8I/zn1u3sJQSim/KKW8UUp5k5SyyBcjCG0kMPxOfTkGRt5J5cFRB4YFMBqVQKaB3QDDsE84kRtgOO0pVwoYNeHB6VRWwPC7cXJVAWPtZHu/hDowIDmVZOTxRd01/MEOmtqD2cLBCTC81AsnwHBbjysBjC0FC9jgIYnqz2fAFzDyEsqHx3sOGBbwTC6ojDyxSCgsNgwYMu2991ILQ5ItBwsoBobMuo8mVChbYGiy8kZP14mmmT8X49KzKyZgBBPFa9OB4cLjc+OAEUzZFgFD6On66R0HCYyqjJQVhCaemQVgz/u7AkuzMGq4CIm1CMwGMNxG24J1YAw/1MHwgx1EGoNh+PzlNLDCoQ9tZ/D+4AxdEy/rlbfvmLeoXVYygDH44VYO39tCKCSIzfoPJGQAY/jBToYeDK5hz13Uy/bG+7bTtqu97PFONfbMAgB73x9cmgYwjv1SP4OP9BGKeKtfVQmLpo4ou2/o9J1OZkVBSonIEXrvne2klvyHklu6klqLGN7a00DvLdvyvhc2sUDLaXEiw65hHTa7b24hvRqMa3RsRmFbr56nwftbmfjJavmToqWrRmo1P3TenuONZJMuphqy1vdh6UqKbT1tALTtjNB7UzDLsBcn0vQM6WntuaWVbMJf6D+ZyQIq8ZkMbb06NIcf2cX4cwt+s0pyIb8nvOeOLtSs/4CpC5ditO3e4fn8enTvuup6D2pLRfdOzKY5+w/XUNP+g7/2393Nnjv04cjsSIzx5xYJIrT1rb++h0iTPg6+/PS8ddRwQDS7ezre9ms7AdAUycjjiySX/PcwtvVGGfy4bjiLzyucf3IpkKdV/+2t9B7T7Tczo0muvBTznSbA+z63Y21DoMvPrbAwHswKUKNs1azGyONLpJadla1M2tul2vub1iJxx6bTnHtsFk3xX2/3393OruP6cOT6qUWu/OtcIBHZb/3Ng57Oq0pYKCmNpq4GBh7uY+Tr11Az/go+vax366ZPL7PreDvZeJYrP17xnc+Va2ladkTJxFRu+FAXSkqziBoOpPXPnNozlLRGbCZLS3eEIw+065XaJzCMbvfshRQ7DjVy6EPbOffksu9Q+MZQafZ8kp6hZpS05OpP/QNjeTJD0/YwSlpy8J5tKGmNxQn/W5RpimR5MkNbT5QjH21n9InywChnwM7Ec2UwEmPHUCuHP7aDc4/N+CpbLZkitajXl+nTy+w+2YmS1rj6oznPaRpKzHoDb1XOhihJlfPffpe23U0Mf2Yv4YZgsjnx7Iw+S3JbJ/t/ZnsgaaoZd3urOp05KdpbNYhZEmD+ckqfJemNMvjx9sAMfuMvxJg6k6DvRAv77mgLJE01Izn33SXicwqHP9xO54FgdhhPLCiOZknc3C+A2dE4l74/z/b+RgZ/vsdz2RbOflz+3jRTry/Rf1c3++71bnPwq6qEBcDC+VjwwDBHDQ8QGG5dw51WwsoBIx08MCRM1BAwyk2rep0OnwsIGPmZkVUBjKqFBVQIGFQHMMAZNGoKGNQ+MNz2JqzkBxilfCo2GxhVDQvY+sCA8tCoA2ODgPGJdpo6gjHjVaSHweYCo+phAe8NYEBpaNSBETwwjPJOXIsx8o0pQmHB8KO7Ng0YTj01NwsYNQELeO8Aw5DMKkXwqAPDHzDWyjJXtmYl5rKM/PP0pgHDrUv3ZgCjZmAB7z1gGDKgIRNJEpMxRv5pGmQdGE6AYS67ctpsYLjVRgOjpmAB711gmJVcyDLyjRwwHuygsdH/QjnYOsBwA4hCbQYw/CwU20hg1BwsYOOA4TcWBmwAMICjv9hLY2PWVyMxVMvA6NjtP84I1HsYdqpJWMDGAOPgR3pqChhNneuV2gwOtw2oWoFR+JtkIomylGD0G9eJz2Q48omddN4YzMKzjQIG2WDc2DcCGJ5bmBCiXwjxnBBiRAhxVgjx2xbH3CuEWM5tbfimEOIP/GU3X3VglAaGWdJY6p3OWDa6wr+5t5c2FBhFcFP17QudAE/NSEa/OV2TwBj6N3tqpofhp3UpwO9KKYeBO9D3DBm2OO7Hua0NT0gpv+zjepaqA8M5MNxq7u0lLj45pwPjgW2IbKoYLBl9qblMlgeQ8Tf+1CxTb67owDjZGMx2fjUEDC2ZYub1eS4+PkX7/paaAYbnliWlvC6lfD33ehUYxWKnsY1QHRiVA8b8hYQOjL5Ghh7uIRQNaC3JDxd1YJxsZ/89AYWRq3JgaMlUXh2aPbNSU8AIpFUJIQ4AtwAvW3x9pxDitBDie0KIo0Fcz0p1YNSBAdUHDAMQdvWmloDhu0UJIdqAbwK/I6UsXPf9OrBfSnkc+HPgsRLpON7r1E51YNSBAdUBjFKAKFStAMPvjmRRdFB8VUr5rcLvpZQrUspY7vWTQFQIYZlrt3ud2qkOjDowYPOA4QYSZtUCMPzMhgjgr4FRKeV/tzmmN3ccQojbc9eb93pNpyoCRkABcAuBITUVKTVf4NgoYDR3eYsLWqhCYIQb6sBo6tDLVmb8BecpAkY0mHobFDD8PHLuAn4FeFsI8Wbus98H9oG+KxnwKeDfCyEUIAl8Wm5Q0E8DGAMP97F9bzCVBIqjhsendVCYgRFqdtczsooaHoQMYAw/uovDHw/O0DV/IQHMcfjjO9YCFweh8R8uAtB3Ug8lF5vy74NgAGPokV0c+cTOtVB9fmUAY/iRHg5/LLiynT2jj+QPP9RLx4HgIrJf/p7+4Oi/q9tzGlUbsPfR+345kLS6jqzP5WuKxtKYRdg7txLQdXg93YWLMfvYiGHnT4dwQ4j2/vXGt3A5gLwCLd3RtadfkOl23bjeC8rEFGLT/sPeFaYbVF4jjaE8sHlOV80P8djS00BTx7qr+cKFYGKQmuttfDq9FhrSl0z1dksF7O24oYXEnP/KF5tK0darV5JQJETzjgbfMT0L0+063EZsysFQRJR/qmViCg1t+m3purGF2LT/p6tSEPS4fV9TUbh5L4pNp2nbpYfBb2iL0JyVgZRtfCZDa4/eALtubCE+kyGIh1o2oRJtCa+l66hsi66bD38lVVC2+1tILgZbb1t3NSLC+I6XChCbTtG2y1tvsGphoaQ1LnznOskAgLH7ZAc3fHQXoFeYIIIAAww+2kf3EX3PkJWrScafnnF8bqmhyp3/cT8AUpNMvrJiHQTYpVp7Grj5s7sBUNMaF5+cI7Xofw+V3e/bzoGf1fd4ycZVRr89gxZA1PChX+ih44A+fFy+muTK80u+04T1stVUybWfLrN4udghzK0NatueJm7+VT1dJa1y/lvvklr0D+O9d3Wx/149GnlmVWH0nyYDAcbRX97r6byqXBuSWVVAg2Of7ad5RzARkQAufXcq8Jie8emUp71VzfPvhZVTSWtMv70auNET4MqPF0EIjj66K7BZEoDLT80HOksipW6zCNroqSmSqdOrxKfzjZ7l/CGcaPyZGUJhwbHP9tPUGYxBGeDiE9d1o+cvBjRL4pE3VQkLTZG8/fdXAwfG/LnVwKdVJcFsxpxXWaVETWQZ+fpVVq4lg50lmc/NkgQMjLnz8ZqZVlXiGc5+7Sqx6ymOfGInHf3BlEFiNs2Zr74TODBmz6xw8YmpYIHhQVUJC4DUQrYiwKi4H0aAu7erGcnoP15bA0bn/kjJHolTrU2rBgyMavPDKCwrvbz0x6qa1hj5h2vErqcY+IW+PIOiHyVmM1sWGFULC6gDA/KBMfBwH92D+ZXaqkE4AclWAoaX3w91YLhVVcMC6sCA8sCwU17jSeuWfy2TWfssPrnKmb+/CsDRR3fR2KIF0nNxCwyrhi5VFanl5+fyd69z/dQifSfb2XfnNt82BqgDw42qHhZQBwZ4B0Y5Jef0Sg25su3OL9tSPRctmUJms5bHzZ5e4MJj19nW18jgJ3eAkvbcAzBr7AczXD+1yJ47uzjwwZ2BlEEdGM5UE7CAOjBg84DhVXOjq1x47Drb9zYz/Om9gQ1J6sDYHGDUDCygDgyoA8NQHRgbD4yaggXUgQF1YBiqA2NjgVFzsIA6MKAODEN1YGwcMGoSFlAHBtSBYagOjI0BRs3CAurAgDowDNWBUXlg1DQsoA4MqAPDUB0YlQVGzcMC6sCAOjAM1YFROWBsCVhAHRhQB4ahOjAqA4wtAwuoAwPqwDBUB0bwwNhSsIA6MKAODEN1YAQLDL9bATwghDgvhLgkhPg9i+8bhRD/mPv+5dxmRBVXITBadjYGkm7Fo4YbwAigrRQCY8fQNv+JUgyMloBgXAiMoGBc28AIpmyLgOExarifrQDCwP8CPgYMA5+x2Ov014FFKeUh4E+B/+r1em5lBkbvrcEETYF8YHQcaCVSgR5GpDEcODD2fsB7VOdCmYGx60RwZWsGRvu+FsIBhcI3AyMUCSZNMzCCLFszMHYdbw8sXTMwvEa79xzdWwhxJ/AlKeVHc+//M4CU8k9Mx/wgd8xLQogIMAXsLLcdQF9fn/yNX/13nvJVqMKnv5r2H3uzMF01q0EQyQrynqiVyGutpVsv28qku9HRvfcA75jeXwPeb3eMlFIRQiwD3cBcYWJCiM8Dn8+9Tf/Rn3z5jI+8Vat2YPHbt4C26u+CrfvbBtyeUDXRvaWUXwG+AiCEOOWWerWg+u+qPW3V3yaEOOX2HD8DuEmg3/R+b+4zy2Nyw5B2NmD7wrrqqit4+YHFq8BhIcRBIUQD8Gng8YJjHgc+l3v9KeCHG7V9YV111RWsPA9DcjaI3wJ+AISBv5FSnhVCfBk4JaV8HH3j5P8nhLgELKADxYm+4jVfVa7676o9bdXf5vp3VeVep3XVVVf1act5cNZVV12VUR0WddVVlyNVFSzKuY/XsoQQE0KIt4UQb3qZtqoWCSH+RggxI4Q4Y/qsSwjxtBDiYu5/52bm0YtsfteXhBCTuXv2phDi45uZRy8SQvQLIZ4TQowIIc4KIX4797nre1Y1sHDoPl7ruk9KeaLG5+3/Fnig4LPfA56VUh4Gns29rzX9LcW/C+BPc/fshJTyyQ3OUxBSgN+VUg4DdwBfzLUr1/esamAB3A5cklKOSSkzwD8An9zkPNVVICnl8+gzW2Z9Evi73Ou/A35+QzMVgGx+V81LSnldSvl67vUqMIruWe36nlUTLKzcx/dsUl4qIQk8JYR4LefavpW0S0p5Pfd6Cti1mZkJWL8lhHgrN0ypueGVWblV37cAL+PhnlUTLLa67pZS3oo+zPqiEOKezc5QJZRzutsq8/F/AdwInACuA/9tc7PjXUKINuCbwO9IKVfM3zm9Z9UECyfu4zUrKeVk7v8M8G30YddW0bQQYjdA7v/MJucnEEkpp6WUqpRSA/6KGr1nQogoOii+KqX8Vu5j1/esmmDhxH28JiWEaBVCbDNeAx8BttKqWrNb/+eA72xiXgKT0ZhyepgavGdCCIHuST0qpfzvpq9c37Oq8uDMTU39D9bdx/94k7MUiIQQN6D3JkB3sf9arf42IcTXgXvRl25PA38IPAb8E7APuAL8opSypoyFNr/rXvQhiAQmgC+Yxvk1ISHE3cCPgbdZjwzy++h2C1f3rKpgUVdddVWvqmkYUldddVWx6rCoq666HKkOi7rqqsuR6rCoq666HKkOi7rqqsuR6rCoq666HKkOi7rqqsuR/j/0pglfNVsx2AAAAABJRU5ErkJggg==\n", + "text/plain": [ + "
" + ] + }, + "metadata": { + "needs_background": "light" + }, + "output_type": "display_data" + } + ], + "source": [ + "#4 Defining variational problem and its solution\n", + "k =1\n", + "u = TrialFunction(V)\n", + "v = TestFunction(V)\n", + "f = Constant(0)\n", + "a = dot(k*grad(u), grad(v))*dx\n", + "L = f*v*dx\n", + "\n", + "# Compute solution\n", + "u = Function(V)\n", + "solve(a == L, u, bcs)\n", + "\n", + "# Plot solution and mesh\n", + "plot(u)\n", + "plot(mesh)\n", + "\n", + "# Save solution to file in VTK format\n", + "vtkfile = File('solution.pvd')\n", + "vtkfile << u" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Results after editing color-map on paraview\n", + "![](paraview-results.png)" + ] + }, + { + "cell_type": "code", + "execution_count": 1, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "\u001b[33mWARNING: The directory '/home/fenics/.cache/pip/http' or its parent directory is not owned by the current user and the cache has been disabled. Please check the permissions and owner of that directory. If executing pip with sudo, you may want sudo's -H flag.\u001b[0m\n", + "\u001b[33mWARNING: The directory '/home/fenics/.cache/pip' or its parent directory is not owned by the current user and caching wheels has been disabled. check the permissions and owner of that directory. If executing pip with sudo, you may want sudo's -H flag.\u001b[0m\n", + "Requirement already satisfied: pandoc in /usr/local/lib/python3.6/dist-packages (1.0.2)\n", + "Requirement already satisfied: ply in /usr/lib/python3/dist-packages (from pandoc) (3.11)\n", + "\u001b[33mWARNING: You are using pip version 19.1, however version 21.0.1 is available.\n", + "You should consider upgrading via the 'pip install --upgrade pip' command.\u001b[0m\n" + ] + } + ], + "source": [ + "!sudo pip install pandoc" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.6.7" + } + }, + "nbformat": 4, + "nbformat_minor": 2 +} diff --git a/Readme.md b/Readme.md new file mode 100644 index 0000000..2d49bf6 --- /dev/null +++ b/Readme.md @@ -0,0 +1,101 @@ + +# Introduction +The Poisson's equation is a second-order partial differential equation that stats the negative Laplacian $-\Delta u$ of an unknown field $u=u(x)$ is equal to a given function $f=f(x)$ on a domain $\Omega \subset \mathbb{R}^d$, most probably defined by a set of boundary conditions for the solution $u$ on the boundary $\partial \Omega$ of $\Omega$: + +$$-\Delta u =f \quad \text{in } \Omega\text{,}$$ +$$u=u_0 \quad \text{on } \Gamma_D \subset \partial\Omega \text{,}$$ + +here the Dirichlet's boundary condition $u=u_0$ signifies a prescribed values for the unknown $u$ on the boundary. + +The Poisson's equation is the simplest model for gravity, electromagnetism, heat transfer, among others. + +The specific case of $f=0$ and a negative $k$ value, leaves to the Fourier's Law. + +## Comparative analysis +Along this example, the fenics platfomr is used to compare results obtained by solving the heat equation (Laplace equation) in 2-D: + +$$\frac{\partial^2 T}{\partial x^2}+ \frac{\partial^2 T}{\partial y^2}=0$$ + +the problem is defined by the next geometry considerations: + +![](physicalproblem.png) + +The resulting contour of temperature, solving using finite diferences, is shown next: + +![](resulteq.png) + +# Solving by Finite Element Method with Varational Problem formulation + + +```python +#1 Loading functions and modules +from fenics import * +import matplotlib.pyplot as plt +``` + + +```python +#2 Create mesh and define function space +mesh = RectangleMesh(Point(0,0),Point(20,20),10, 10,'left') +V = FunctionSpace(mesh, 'Lagrange', 1) #Lagrange are triangular elements +plot(mesh) +plt.show() +``` + + +![png](output_6_0.png) + + + +```python +#3 Defining boundary conditions (Dirichlet) +tol = 1E-14 # tolerance for coordinate comparisons +#at y=20 +def Dirichlet_boundary1(x, on_boundary): + return on_boundary and abs(x[1] - 20) < tol +#at y=0 +def Dirichlet_boundary0(x, on_boundary): + return on_boundary and abs(x[1] - 0) < tol +#at x=0 +def Dirichlet_boundarx0(x, on_boundary): + return on_boundary and abs(x[0] - 0) < tol +#at x=20 +def Dirichlet_boundarx1(x, on_boundary): + return on_boundary and abs(x[0] - 20) < tol + +bc0 = DirichletBC(V, Constant(0), Dirichlet_boundary0) +bc1 = DirichletBC(V, Constant(100), Dirichlet_boundary1) #100C +bc2 = DirichletBC(V, Constant(0), Dirichlet_boundarx0) +bc3 = DirichletBC(V, Constant(0), Dirichlet_boundarx1) +bcs = [bc0,bc1, bc2,bc3] +``` + + +```python +#4 Defining variational problem and its solution +k =1 +u = TrialFunction(V) +v = TestFunction(V) +f = Constant(0) +a = dot(k*grad(u), grad(v))*dx +L = f*v*dx + +# Compute solution +u = Function(V) +solve(a == L, u, bcs) + +# Plot solution and mesh +plot(u) +plot(mesh) + +# Save solution to file in VTK format +vtkfile = File('solution.pvd') +vtkfile << u +``` + + +![png](output_8_0.png) + + +# Results after editing color-map on paraview +![](paraview-results.png) diff --git a/output_6_0.png b/output_6_0.png new file mode 100644 index 0000000..0672190 Binary files /dev/null and b/output_6_0.png differ diff --git a/output_8_0.png b/output_8_0.png new file mode 100644 index 0000000..0c85d5e Binary files /dev/null and b/output_8_0.png differ diff --git a/paraview-results.png b/paraview-results.png new file mode 100644 index 0000000..a661427 Binary files /dev/null and b/paraview-results.png differ diff --git a/physicalproblem.PNG b/physicalproblem.PNG new file mode 100644 index 0000000..5318246 Binary files /dev/null and b/physicalproblem.PNG differ diff --git a/resulteq.png b/resulteq.png new file mode 100644 index 0000000..8bcc2ed Binary files /dev/null and b/resulteq.png differ diff --git a/solution.pvd b/solution.pvd new file mode 100644 index 0000000..1ecfa5f --- /dev/null +++ b/solution.pvd @@ -0,0 +1,6 @@ + + + + + +