tasks/task13.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Modelling discontinuous trapped concentration profiles\n",
    "\n",
    "It is possible to set the finite element type of trapped concentration in FESTIM.\n",
    "But why even bother?\n",
    "\n",
    "It basically helps have a better approximation of concentration fields in presence of discontinuities.\n",
    "To illustrate this, we will present two equivalent examples.\n",
    "\n",
    "A multi-material case with a trap defined in only one subdomain.\n",
    "And a single-material with a trap defined only in a subregion."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Multi-material\n",
    "\n",
    "Let's define a case with 2 materials on $\\Omega = [0, 1]$ with the interface at $x=0.5$.\n",
    "\n",
    "We will set the trap only in the first subdomain $x \\in [0, 0.5]$.\n",
    "\n",
    "Let's start with the default argument in ``F.Settings``: using CG (Continuous Galerkin) finite elements for the trapped concentration."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Defining initial values\n",
      "Defining variational problem\n",
      "Defining source terms\n",
      "Defining boundary conditions\n",
      "Solving steady state problem...\n",
      "Solved problem in 0.00 s\n"
     ]
    }
   ],
   "source": [
    "import festim as F\n",
    "import numpy as np\n",
    "\n",
    "my_model = F.Simulation()\n",
    "\n",
    "my_model.mesh = F.MeshFromVertices(np.linspace(0, 1, 100))\n",
    "my_model.materials = [\n",
    "    F.Material(id=1, D_0=1, E_D=0, borders=[0, 0.5]),\n",
    "    F.Material(id=2, D_0=1, E_D=0, borders=[0.5, 1]),\n",
    "]\n",
    "\n",
    "my_model.traps = F.Trap(\n",
    "    k_0=1,\n",
    "    E_k=0,\n",
    "    p_0=1,\n",
    "    E_p=0,\n",
    "    density=1,\n",
    "    materials=my_model.materials[0],\n",
    ")\n",
    "\n",
    "my_model.T = 300\n",
    "\n",
    "my_model.boundary_conditions = [F.DirichletBC(surfaces=[1, 2], value=1, field=\"solute\")]\n",
    "\n",
    "my_model.settings = F.Settings(\n",
    "    absolute_tolerance=1e-10,\n",
    "    relative_tolerance=1e-10,\n",
    "    transient=False,\n",
    ")\n",
    "\n",
    "my_model.exports = [F.TXTExport(filename=\"task13/trapped_concentration_multi.txt\", field=1)]\n",
    "\n",
    "my_model.initialise()\n",
    "my_model.run()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It seems to have solved without any issue. Let's have a look at the produced concentration profile for the trapped particles."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/home/remidm/miniconda3/envs/festim-workshop/lib/python3.11/site-packages/IPython/core/pylabtools.py:77: DeprecationWarning: backend2gui is deprecated since IPython 8.24, backends are managed in matplotlib and can be externally registered.\n",
      "  warnings.warn(\n",
      "/home/remidm/miniconda3/envs/festim-workshop/lib/python3.11/site-packages/IPython/core/pylabtools.py:77: DeprecationWarning: backend2gui is deprecated since IPython 8.24, backends are managed in matplotlib and can be externally registered.\n",
      "  warnings.warn(\n",
      "/home/remidm/miniconda3/envs/festim-workshop/lib/python3.11/site-packages/IPython/core/pylabtools.py:77: DeprecationWarning: backend2gui is deprecated since IPython 8.24, backends are managed in matplotlib and can be externally registered.\n",
      "  warnings.warn(\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAkMAAAGwCAYAAACq12GxAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAAA76UlEQVR4nO3de3hU1b3/8c8kM5MLkiAGIpcIAZWLgGCoEPhRpGpEaanaKpYeUAGP1HIQI0U4qFyOFG0toFWgeJAeLcV4wdZajhKpIHdLSNQWDlpEEyABApKLQC6T/fsjzOQKZse5MFnv1/PMU2bPnuSbXcx8+K6113JYlmUJAADAUBGhLgAAACCUCEMAAMBohCEAAGA0whAAADAaYQgAABiNMAQAAIxGGAIAAEZzhrqAC11VVZUOHz6s1q1by+FwhLocAADQBJZlqaSkRB07dlRExPl7P4Shb3D48GElJSWFugwAANAMeXl56ty583nPIQx9g9atW0uqvphxcXEhrgYAADRFcXGxkpKSfJ/j50MY+gbeobG4uDjCEAAAYaYpU1yYQA0AAIxGGAIAAEYjDAEAAKMRhgAAgNEIQwAAwGiEIQAAYDTCEAAAMBphCAAAGI0wBAAAjEYYAgAARiMMAQAAoxGGAACA0QhDAIxzutwjy7JCXQaACwRhCIBR/nW0RH3mvqtH//SPUJcC4AJBGAJglEWZn8pTZWn1ztxQlwLgAkEYAmCU/KIzoS4BwAWGMATAKPknCUMA6iIMATBKQTFhCEBdhCEAxiivrAp1CQAuQIQhAMbIPXGqzvMKD+EIAGEIgEE+P1Za5/npCk+IKgFwISEMATDGgcKv6zw/XU4YAkAYAmCQz48RhgA0RBgCYIzPC+sOk50iDAEQYQiAQRp0hpgzBECEIQCGKDpVoeNfl0uSOrWJkSSdIQwBEGEIgCH2nx0iS4yLUkLrKEkMkwGo5gx1AQAQDN4hsm4JF/mOMUwGQKIzBMAQB852hpLbtVKMO1KSdLq8MpQlAbhAEIYAGKGmM9RKMS5vGKIzBIAwBMAQ3jDUvd1Fvs7QKYbJAIgwBMAAnipLB46f7Qy1q+kMnaEzBECEIQAGOHzytMorq+SOjFDni2MV650zRGcIgAhDAAyw/+wGrV0uiVVkhEPRZztD3FoPQArDMLR06VIlJycrOjpaKSkp2rx58znP3bhxoxwOR4PH//3f/wWxYgDBVHSqQpZl1Tnm3aC1W7tWklRzNxmdIQAKszCUkZGhadOmafbs2crOztawYcN08803Kzc397zv27dvn/Lz832PK664IkgVAwgGy7K0bX+hxq3cqavnr9e8v+yp87p38nTy2TWGfMNkdIYAKMwWXVy0aJEmTpyoSZMmSZKWLFmid999V8uWLdPChQvP+b727durTZs2QaqyabK+PKHt+4/bfp/D4QhANfW/h83zZb+m5vwYdt/SvO9h701B+L/j7PexWVezvofN85v1Pey964vjXyvhoiiVVXh0prJKFZ4qtW8draS2Mep8caw6Xxyj7NyTev79fykn76Tvff+z/Qv9OKWz+nSKl1SzQau3M+QdJqMzBEAKozBUXl6urKwszZw5s87xtLQ0bdu27bzvHTBggM6cOaPevXvr0Ucf1YgRI855bllZmcrKynzPi4uLv13h57DzwAk9vf7TgHxtwERuZ4TGDExSftEZvbf3iJ746x6tuW+wHA5Hrdvqq8MQnSEAtYVNGCosLJTH41FiYmKd44mJiSooKGj0PR06dNCKFSuUkpKisrIyvfzyy7r++uu1ceNGffe73230PQsXLtS8efP8Xn99vTrE6SfXJtl6T71pEH4/X5Is2XtT875HM95j92dvznex/T2a8S2accHsviMY/58E4+f4v/xiXdUxXtGuCEW7IhXtipTDIR0pOqODX53Wwa9O60jJGbVyOzUutYsmDE1Wu9ZROnTytDZ/dkw7Pj+h9XuOaNgVCcovOiOpZiuOGDpDAGoJmzDkVb/NblnWOVvvPXr0UI8ePXzPU1NTlZeXp6effvqcYWjWrFlKT0/3PS8uLlZSkr3Q0hQjerTXiB7t/f51AZOUVXoU4XDIFVkz/bFTmxjdN6ybnnv/X1q4bq/aj+kvSbo41qWLW7kl1ZpATWcIgMJoAnVCQoIiIyMbdIGOHj3aoFt0PoMHD9Znn312ztejoqIUFxdX5wHgwhTljKwThLx+dl13tWsdpS+On9Lcs5Opu7Wr2aA1hlvrAdQSNmHI7XYrJSVFmZmZdY5nZmZqyJAhTf462dnZ6tChg7/LA3ABaRXl1C/SqrvCH52dWJ2c0Mr3urczdIZhMgAKs2Gy9PR0jRs3TgMHDlRqaqpWrFih3NxcTZ48WVL1ENehQ4f00ksvSaq+26xr16666qqrVF5erj/84Q9644039MYbb4TyxwAQBD9K6azfb/tCe/Krb4Lw3kkmSbGu6l99zBkCIIVZGBozZoyOHz+u+fPnKz8/X3369NG6devUpUsXSVJ+fn6dNYfKy8s1ffp0HTp0SDExMbrqqqv017/+VbfcckuofgQAQRIZ4dCj3++lsS/slCQlX1IThqLd1U3x0xWe8847BGAGh9Wc20IMUlxcrPj4eBUVFTF/CAhD8/+yR9v2F+qNnw1Rq6jqf/+VllWqz5x3JUn/918jfesOAWg57Hx+h1VnCADsevwHvRsci3bWTJc8Ve4hDAGGC5sJ1ADgL87ICLkja4bKAJiNMATASKw1BMCLMATASL5VqAlDgPEIQwCM5NufjGEywHiEIQBGivatQl0Z4koAhBphCICRYlmFGsBZhCEARvJOoGZ/MgCEIQBG8g6TMWcIAGEIgJFiubUewFmEIQBG4tZ6AF6EIQBGiuHWegBnEYYAGCnGxQRqANUIQwCM5A1D3FoPgDAEwEgMkwHwIgwBMBLrDAHwIgwBMBIrUAPwIgwBMBITqAF4EYYAGCnG7ZTEOkMACEMADMXdZAC8CEMAjMQwGQAvwhAAI3FrPQAvwhAAI8WwUSuAswhDAIwUe3aYrNxTpUpPVYirARBKhCEARvJ2hiSGygDTEYYAGCnKGSGHo/rPhCHAbIQhAEZyOBw1t9eXM0wGmIwwBMBYvtvrKypDXAmAUCIMATAWd5QBkAhDAAzm7QwRhgCzEYYAGCuWhRcBiDAEwGDRLsIQAMIQAIN55wyxPxlgNsIQAGN5h8nYuR4wG2EIgLGi2bkegAhDAAwWy631AEQYAmAw3wrUDJMBRiMMATBWDMNkAEQYAmCwGLdTErfWA6YjDAEwVoyr+lcgc4YAsxGGABgrls4QABGGABgs2rfoIrvWAyYLuzC0dOlSJScnKzo6WikpKdq8eXOT3rd161Y5nU71798/sAUCCBuxvu04qkJcCYBQCqswlJGRoWnTpmn27NnKzs7WsGHDdPPNNys3N/e87ysqKtL48eN1/fXXB6lSAOHAux3HGeYMAUYLqzC0aNEiTZw4UZMmTVKvXr20ZMkSJSUladmyZed93/3336+xY8cqNTU1SJUCCAe+FagrGCYDTBY2Yai8vFxZWVlKS0urczwtLU3btm075/tWrVql/fv3a86cOU36PmVlZSouLq7zANAy1axAzTAZYLKwCUOFhYXyeDxKTEysczwxMVEFBQWNvuezzz7TzJkztXr1ajmdziZ9n4ULFyo+Pt73SEpK+ta1A7gweRddPM0EasBoYROGvBwOR53nlmU1OCZJHo9HY8eO1bx583TllVc2+evPmjVLRUVFvkdeXt63rhnAhcnXGarwyLKsEFcDIFSa1i65ACQkJCgyMrJBF+jo0aMNukWSVFJSol27dik7O1tTpkyRJFVVVcmyLDmdTq1fv17f+973GrwvKipKUVFRgfkhAFxQvLfWV1lSuadKUc7IEFcEIBTCpjPkdruVkpKizMzMOsczMzM1ZMiQBufHxcXpk08+UU5Oju8xefJk9ejRQzk5ORo0aFCwSgdwgfIOk0msQg2YLGw6Q5KUnp6ucePGaeDAgUpNTdWKFSuUm5uryZMnS6oe4jp06JBeeuklRUREqE+fPnXe3759e0VHRzc4DsBMrsgIuSIdqvBYOl3hUZtQFwQgJMIqDI0ZM0bHjx/X/PnzlZ+frz59+mjdunXq0qWLJCk/P/8b1xwCgNqiXZGq8FSycz1gMIfFrMHzKi4uVnx8vIqKihQXFxfqcgD42aBfvqcjxWV6+z/+n/p0ig91OQD8xM7nd9jMGQKAQPDOGzrDZq2AsQhDAIwWc3bneobJAHMRhgAYLcZV/WvwNJ0hwFiEIQBGi/FtyUEYAkxFGAJgtBhX9TAZnSHAXIQhAEbzdoaYMwSYizAEwGix3E0GGI8wBMBozBkCQBgCYDSGyQA0azuOTz/9VBs3btTRo0dVVVVV57XHH3/cL4UBQDB4F11kAjVgLtth6IUXXtDPfvYzJSQk6NJLL5XD4fC95nA4CEMAwoovDJVXhrgSAKFiOww98cQTWrBggR555JFA1AMAQeWbM0RnCDCW7TlDX331le64445A1AIAQVczTFb1DWcCaKlsh6E77rhD69evD0QtABB0sW6GyQDT2R4mu/zyy/XYY49px44d6tu3r1wuV53Xp06d6rfiACDQohkmA4xnOwytWLFCF110kTZt2qRNmzbVec3hcBCGAIQV7zAZt9YD5rIdhg4cOBCIOgAgJLzDZGcIQ4CxvtWii5ZlybIsf9UCAEHHOkMAmhWGXnrpJfXt21cxMTGKiYlRv3799PLLL/u7NgAIOFagBmB7mGzRokV67LHHNGXKFA0dOlSWZWnr1q2aPHmyCgsL9dBDDwWiTgAICG9nqKyySlVVliIiHN/wDgAtje0w9Nvf/lbLli3T+PHjfcd++MMf6qqrrtLcuXMJQwDCSqy75tfg6QqPWkU1a5ciAGHM9jBZfn6+hgwZ0uD4kCFDlJ+f75eiACBYopw1vwaZNwSYyXYYuvzyy/Xqq682OJ6RkaErrrjCL0UBQLBERDgU7ar+VXiaeUOAkWz3g+fNm6cxY8bogw8+0NChQ+VwOLRlyxZt2LCh0ZAEABe6WLdTZyrK6QwBhrLdGfrRj36knTt3KiEhQX/605+0du1aJSQk6MMPP9Rtt90WiBoBIKBqdq4nDAEmatZMwZSUFP3hD3/wdy0AEBLcXg+YrUlhqLi4WHFxcb4/n4/3PAAIF97O0BmGyQAjNSkMXXzxxcrPz1f79u3Vpk0bORwN1+GwLEsOh0MeD79MAIQX9icDzNakMPS3v/1Nbdu2lSS9//77AS0IAIIt6uzdZOX8Yw4wUpPC0PDhw31/Tk5OVlJSUoPukGVZysvL8291ABAEzrOrTld42GsRMJHtu8mSk5N17NixBsdPnDih5ORkvxQFAMHkjKz+VVhJGAKMZDsMeecG1VdaWqro6Gi/FAUAweSK9HaGqkJcCYBQaPKt9enp6ZIkh8Ohxx57TLGxsb7XPB6Pdu7cqf79+/u9QAAINNfZzhBhCDBTk8NQdna2pOrO0CeffCK32+17ze126+qrr9b06dP9XyEABJgz4uwwWRXDZICJmhyGvHeR3XvvvXrmmWdYTwhAi+EdJqukMwQYyfYK1KtWrQpEHQAQMs5I7iYDTNas7Tj+/ve/67XXXlNubq7Ky8vrvLZ27Vq/FAYAweIdJmPOEGAm23eTvfLKKxo6dKj27NmjN998UxUVFdqzZ4/+9re/KT4+PhA1AkBAuZ3MGQJMZjsM/fKXv9TixYv19ttvy+1265lnntHevXt155136rLLLgtEjQAQUDWLLtIZAkxkOwzt379fo0aNkiRFRUXp66+/lsPh0EMPPaQVK1b4vUAACDQWXQTMZjsMtW3bViUlJZKkTp066R//+Ick6eTJkzp16pR/qwOAIHCd7QxVVtEZAkxkewL1sGHDlJmZqb59++rOO+/Ugw8+qL/97W/KzMzU9ddfH4gaASCgXGfnDJVX0hkCTGS7M/Tcc8/prrvukiTNmjVL06dP15EjR3T77bdr5cqVfi+wvqVLlyo5OVnR0dFKSUnR5s2bz3nuli1bNHToUF1yySWKiYlRz549tXjx4oDXCCC8OOkMAUaz1RmqrKzUX/7yF910002SpIiICM2YMUMzZswISHH1ZWRkaNq0aVq6dKmGDh2q3/3ud7r55pu1Z8+eRidvt2rVSlOmTFG/fv3UqlUrbdmyRffff79atWqlf//3fw9KzQAufC7mDAFGc1iWZeu//tjYWO3du1ddunQJVE3nNGjQIF1zzTVatmyZ71ivXr106623auHChU36GrfffrtatWqll19+uUnnFxcXKz4+XkVFRay6DbRQq3d+qdlv/kNpvRO1YvzAUJcDwA/sfH7bHiYbNGiQb5+yYCovL1dWVpbS0tLqHE9LS9O2bdua9DWys7O1bds2DR8+/JznlJWVqbi4uM4DQMvmYtFFwGi2J1A/8MADevjhh3Xw4EGlpKSoVatWdV7v16+f34qrrbCwUB6PR4mJiXWOJyYmqqCg4Lzv7dy5s44dO6bKykrNnTtXkyZNOue5Cxcu1Lx58/xSM4Dw4HJ65wwxTAaYyHYYGjNmjCRp6tSpvmMOh0OWZcnhcMjj8fivukY4HI46z73f93w2b96s0tJS7dixQzNnztTll1+un/zkJ42eO2vWLKWnp/ueFxcXKykp6dsXDuCCxXYcgNlsh6EDBw4Eoo5vlJCQoMjIyAZdoKNHjzboFtWXnJwsSerbt6+OHDmiuXPnnjMMRUVFKSoqyj9FAwgLNbvW0xkCTGR7ztCXX36pTp06qUuXLnUenTp10pdffhmIGiVJbrdbKSkpyszMrHM8MzNTQ4YMafLXsSxLZWVl/i4PQBjzdYYYJgOMZLszNGLECOXn56t9+/Z1jhcVFWnEiBEBHSZLT0/XuHHjNHDgQKWmpmrFihXKzc3V5MmTJVUPcR06dEgvvfSSJOn555/XZZddpp49e0qqXnfo6aef1n/8x38ErEYA4cd5tjNUUckwGWAi22HoXHN0jh8/3mAytb+NGTNGx48f1/z585Wfn68+ffpo3bp1vtv88/PzlZub6zu/qqpKs2bN0oEDB+R0OtW9e3c9+eSTuv/++wNaJ4Dw4vauM8Sii4CRmrzO0O233y5J+vOf/6yRI0fWmVfj8Xj08ccfq0ePHnrnnXcCU2mIsM4Q0PJ9eOCE7vzddnVLaKW/Tb8u1OUA8AM7n99N7gzFx8dLqu4MtW7dWjExMb7X3G63Bg8erPvuu6+ZJQNA6PiGyegMAUZqchhatWqVJKlr166aPn16wIfEACBYvIsucjcZYCbbc4bmzJkTiDoAIGR8nSHWGQKMZPvW+iNHjmjcuHHq2LGjnE6nIiMj6zwAINx4N2qtoDMEGMl2Z+iee+5Rbm6uHnvsMXXo0OEbV38GgAtdzaKLdIYAE9kOQ1u2bNHmzZvVv3//AJQDAMHnjGTRRcBktofJkpKS1MS78QEgLLgi6AwBJrMdhpYsWaKZM2fqiy++CEA5ABB83s5QlSV56A4BxmnWrvWnTp1S9+7dFRsbK5fLVef1EydO+K04AAgG75whqfqOssgIbgYBTGI7DC1ZsiQAZQBA6HjvJpOkSjpDgHFsh6G77747EHUAQMg4I2o6Q8wbAsxje86QJO3fv1+PPvqofvKTn+jo0aOSpHfeeUf//Oc//VocAARDZETtYTI6Q4BpbIehTZs2qW/fvtq5c6fWrl2r0tJSSdLHH3/M6tQAwpLD4fDNG2IVasA8tsPQzJkz9cQTTygzM1Nut9t3fMSIEdq+fbtfiwOAYPHOG2J/MsA8tsPQJ598ottuu63B8Xbt2un48eN+KQoAgs07b4id6wHz2A5Dbdq0UX5+foPj2dnZ6tSpk1+KAoBgozMEmMt2GBo7dqweeeQRFRQUyOFwqKqqSlu3btX06dM1fvz4QNQIAAHHzvWAuWyHoQULFuiyyy5Tp06dVFpaqt69e+u73/2uhgwZokcffTQQNQJAwDkjvDvXE4YA09heZ8jlcmn16tX6r//6L+3evVtVVVUaMGCArrjiikDUBwBB4XaeHSZj0UXAOLbDkFe3bt3UrVs3f9YCACHjm0BNZwgwju1hsh//+Md68sknGxz/9a9/rTvuuMMvRQFAsDmZQA0Yq1mLLo4aNarB8ZEjR+qDDz7wS1EAEGzeRRcrubUeMI7tMFRaWlpnsUUvl8ul4uJivxQFAMHmHSYrr6QzBJjGdhjq06ePMjIyGhx/5ZVX1Lt3b78UBQDB5ltniM4QYBzbE6gfe+wx/ehHP9L+/fv1ve99T5K0YcMGrVmzRq+99prfCwSAYGDRRcBctsPQ6NGj9ac//Um//OUv9frrrysmJkb9+vXTe++9p+HDhweiRgAIOBZdBMzVrFvrR40a1egkagAIV95FF1lnCDBPs9cZKi8v19GjR1VVb3z9sssu+9ZFAUCwuegMAcayHYY+++wzTZgwQdu2batz3LIsORwOeTwevxUHAMHinTNUwZwhwDi2w9A999wjp9Opt99+Wx06dJDD4QhEXQAQVN45Q5V0hgDj2A5DOTk5ysrKUs+ePQNRDwCEhIs5Q4CxbK8z1Lt3bxUWFgaiFgAIGW9nqLySzhBgGtth6KmnntKMGTO0ceNGHT9+XMXFxXUeABCOWHQRMJftYbIbbrhBknT99dfXOc4EagDhzLc3GROoAePYDkPvv/9+IOoAgJBycjcZYCzbYYhVpgG0RK4Idq0HTNWsRRdPnjyplStXau/evXI4HOrdu7cmTJig+Ph4f9cHAEFR0xkiDAGmsT2BeteuXerevbsWL16sEydOqLCwUIsWLVL37t21e/fuQNQIAAFXszcZw2SAaWx3hh566CGNHj1aL7zwgpzO6rdXVlZq0qRJmjZtmj744AO/FwkAgeb27VpPZwgwje0wtGvXrjpBSJKcTqdmzJihgQMH+rU4AAgW59k5QxUsuggYx/YwWVxcnHJzcxscz8vLU+vWrf1SFAAEm5POEGAs22FozJgxmjhxojIyMpSXl6eDBw/qlVde0aRJk/STn/wkEDUCQMC5mDMEGMt2GHr66ad1++23a/z48eratau6dOmie+65Rz/+8Y/11FNPBaLGOpYuXark5GRFR0crJSVFmzdvPue5a9eu1Y033qh27dopLi5OqampevfddwNeI4Dw44zgbjLAVLbDkNvt1jPPPKOvvvpKOTk5ys7O1okTJ7R48WJFRUUFokafjIwMTZs2TbNnz1Z2draGDRumm2++udFhO0n64IMPdOONN2rdunXKysrSiBEj9IMf/EDZ2dkBrRNA+HE5vcNkdIYA0zgsy7L1X35RUZE8Ho/atm1b5/iJEyfkdDoVFxfn1wJrGzRokK655hotW7bMd6xXr1669dZbtXDhwiZ9jauuukpjxozR448/3qTzi4uLFR8fr6KiooD+bABC638/ydfPVu/Wd7perNcmDwl1OQC+JTuf37Y7Q3fddZdeeeWVBsdfffVV3XXXXXa/XJOVl5crKytLaWlpdY6npaVp27ZtTfoaVVVVKikpaRDkaisrK2PzWcBAbMcBmMt2GNq5c6dGjBjR4Ph1112nnTt3+qWoxhQWFsrj8SgxMbHO8cTERBUUFDTpa/zmN7/R119/rTvvvPOc5yxcuFDx8fG+R1JS0reqG0B4qFl0kTlDgGlsh6GysjJVVlY2OF5RUaHTp0/7pajzcTgcdZ5bltXgWGPWrFmjuXPnKiMjQ+3btz/nebNmzVJRUZHvkZeX961rBnDhc0UwZwgwle0w9J3vfEcrVqxocHz58uVKSUnxS1GNSUhIUGRkZIMu0NGjRxt0i+rLyMjQxIkT9eqrr+qGG24477lRUVGKi4ur8wDQ8vlurWejVsA4tlegXrBggW644QZ99NFHuv766yVJGzZs0N///netX7/e7wV6ud1upaSkKDMzU7fddpvveGZmpn74wx+e831r1qzRhAkTtGbNGo0aNSpg9QEIbzWLLtIZAkxjuzM0dOhQbd++XUlJSXr11Vf1l7/8RZdffrk+/vhjDRs2LBA1+qSnp+u///u/9eKLL2rv3r166KGHlJubq8mTJ0uqHuIaP3687/w1a9Zo/Pjx+s1vfqPBgweroKBABQUFKioqCmidAMKPtzPECtSAeWx3hiSpf//+Wr16tb9r+UZjxozR8ePHNX/+fOXn56tPnz5at26dunTpIknKz8+vs+bQ7373O1VWVurnP/+5fv7zn/uO33333fr9738f7PIBXMC8iy6W0xkCjGN7nSHTsM4QYIbPjpToxsUfqE2sSzmPp33zGwBc0AK6zhAAtEQu5gwBxiIMAYBYZwgwGWEIAFSrM1RFZwgwDWEIACQ5I6o7Q54qS1UEIsAoTbqb7Pbbb2/yF1y7dm2ziwGAUPGuMyRVL7wYFREZwmoABFOTOkO19+qKi4vThg0btGvXLt/rWVlZ2rBhg+Lj4wNWKAAEkrtWGGISNWCWJnWGVq1a5fvzI488ojvvvFPLly9XZGT1v5w8Ho8eeOABbj0HELa8E6glwhBgGttzhl588UVNnz7dF4QkKTIyUunp6XrxxRf9WhwABIt3zpDE/mSAaWyHocrKSu3du7fB8b1796qKXyAAwpTD4fAFIm6vB8xiezuOe++9VxMmTNC//vUvDR48WJK0Y8cOPfnkk7r33nv9XiAABIsz0qHKKothMsAwtsPQ008/rUsvvVSLFy9Wfn6+JKlDhw6aMWOGHn74Yb8XCADB4oqM0JmKKjpDgGFsh6GIiAjNmDFDM2bMUHFxsSQxcRpAi8DCi4CZmrXoYmVlpd577z2tWbNGDkf1GPvhw4dVWlrq1+IAIJi8c4bKK+kMASax3Rn68ssvNXLkSOXm5qqsrEw33nijWrdurV/96lc6c+aMli9fHog6ASDg6AwBZrLdGXrwwQc1cOBAffXVV4qJifEdv+2227Rhwwa/FgcAweRda6iSOUOAUWx3hrZs2aKtW7fK7XbXOd6lSxcdOnTIb4UBQLB5O0MV3E0GGMV2Z6iqqkoej6fB8YMHD6p169Z+KQoAQsE7Z6iSNdMAo9gOQzfeeKOWLFnie+5wOFRaWqo5c+bolltu8WdtABBUNZ0hwhBgEtvDZIsXL9aIESPUu3dvnTlzRmPHjtVnn32mhIQErVmzJhA1AkBQeOcMMUwGmMV2GOrYsaNycnK0Zs0a7d69W1VVVZo4caJ++tOf1plQDQDhxhVx9m4ywhBgFNthSJJiYmI0YcIETZgwwd/1AEDIuJzMGQJM1KwwtG/fPv32t7/V3r175XA41LNnT02ZMkU9e/b0d30AEDTOCO4mA0xkewL166+/rj59+igrK0tXX321+vXrp927d6tv37567bXXAlEjAASFK5Jd6wET2e4MzZgxQ7NmzdL8+fPrHJ8zZ44eeeQR3XHHHX4rDgCCyembM0QYAkxiuzNUUFCg8ePHNzj+b//2byooKPBLUQAQCtxNBpjJdhi67rrrtHnz5gbHt2zZomHDhvmlKAAIBbdvbzI6Q4BJbA+TjR49Wo888oiysrI0ePBgSdKOHTv02muvad68eXrrrbfqnAsA4YLOEGAmh2VZtv6rj4hoWjPJ4XA0um1HuCkuLlZ8fLyKiooUFxcX6nIABNB/vvmJ/rgzV9NuuELTbrgy1OUA+BbsfH7b7gxV0T4G0EK5vHuT0RkCjGJ7zhAAtFRO795k/KMPMEqzwtCGDRv0/e9/X927d9fll1+u73//+3rvvff8XRsABJV3o1Y6Q4BZbIeh5557TiNHjlTr1q314IMPaurUqYqLi9Mtt9yi5557LhA1AkBQeBddZJ0hwCy25wwtXLhQixcv1pQpU3zHpk6dqqFDh2rBggV1jgNAOPEuulhOZwgwiu3OUHFxsUaOHNngeFpamoqLi/1SFACEgpPOEGAk22Fo9OjRevPNNxsc//Of/6wf/OAHfikKAELBN0xWRWcIMIntYbJevXppwYIF2rhxo1JTUyVVL7q4detWPfzww3r22Wd9506dOtV/lQJAgHknULNRK2AW22Fo5cqVuvjii7Vnzx7t2bPHd7xNmzZauXKl77nD4SAMAQgrTu4mA4xkOwwdOHAgEHUAQMh5F12kMwSYhUUXAeCsmkUX6QwBJrHdGZKkgwcP6q233lJubq7Ky8vrvLZo0SK/FAYAwcY6Q4CZbIehDRs2aPTo0UpOTta+ffvUp08fffHFF7IsS9dcc00gagSAoGAFasBMtofJZs2apYcfflj/+Mc/FB0drTfeeEN5eXkaPny47rjjjkDUCABB4fTOGWJvMsAotsPQ3r17dffdd0uSnE6nTp8+rYsuukjz58/XU0895fcC61u6dKmSk5MVHR2tlJQUbd68+Zzn5ufna+zYserRo4ciIiI0bdq0gNcHIHxxaz1gJtthqFWrViorK5MkdezYUfv37/e9VlhY6L/KGpGRkaFp06Zp9uzZys7O1rBhw3TzzTcrNze30fPLysrUrl07zZ49W1dffXVAawMQ/mpWoGaYDDCJ7TA0ePBgbd26VZI0atQoPfzww1qwYIEmTJigwYMH+73A2hYtWqSJEydq0qRJ6tWrl5YsWaKkpCQtW7as0fO7du2qZ555RuPHj1d8fHxAawMQ/rx7k9EZAsxiewL1okWLVFpaKkmaO3euSktLlZGRocsvv1yLFy/2e4Fe5eXlysrK0syZM+scT0tL07Zt2/z2fcrKynydL0nstwYYxO1kOw7ARLbCkMfjUV5envr16ydJio2N1dKlSwNSWH2FhYXyeDxKTEysczwxMVEFBQV++z4LFy7UvHnz/Pb1AIQPb2eIYTLALLaGySIjI3XTTTfp5MmTASrnmzkcjjrPLctqcOzbmDVrloqKinyPvLw8v31tABc275yhcobJAKPYHibr27evPv/8cyUnJweinnNKSEhQZGRkgy7Q0aNHG3SLvo2oqChFRUX57esBCB816wwRhgCT2J5AvWDBAk2fPl1vv/228vPzVVxcXOcRKG63WykpKcrMzKxzPDMzU0OGDAnY9wVgDu86QwyTAWax3RkaOXKkJGn06NF1hqe8w1Uej8d/1dWTnp6ucePGaeDAgUpNTdWKFSuUm5uryZMnS6oe4jp06JBeeukl33tycnIkSaWlpTp27JhycnLkdrvVu3fvgNUJIDz51hli0UXAKLbD0Pvvvx+IOppkzJgxOn78uObPn6/8/Hz16dNH69atU5cuXSRVL7JYf82hAQMG+P6clZWlP/7xj+rSpYu++OKLYJYOIAzULLpIZwgwicOyLFv/1efm5iopKanRicx5eXm67LLL/FpgqBUXFys+Pl5FRUWKi4sLdTkAAqiwtEwDn3hPknRg4S1+vTkDQHDZ+fy2PWcoOTlZx44da3D8xIkTQZ9UDQD+5Iqo+ZVIdwgwh+0wdK5b2UtLSxUdHe2XogAgFLy31ktSJfOGAGM0ec5Qenq6pOp1fh577DHFxsb6XvN4PNq5c6f69+/v9wIBIFi8c4YkOkOASZochrKzsyVVd4Y++eQTud1u32tut1tXX321pk+f7v8KASBIXLU6Q+xPBpijyWHIexfZvffeq2eeeYbJxABaHIfDocgIhzxVFmsNAQaxfWv9qlWrAlEHAFwQnGfDEJ0hwBy2J1ADQEvm25KDnesBYxCGAKAW77wh9icDzEEYAoBanGc7Q+xcD5iDMAQAtbjYrBUwDmEIAGpx+uYM0RkCTEEYAoBavKtQs+giYA7CEADU4vZ2hghDgDEIQwBQS01niGEywBSEIQCoxXl253rCEGAOwhAA1OJbZ4hFFwFjEIYAoBY6Q4B5CEMAUIvLyQRqwDSEIQCoxbvoIp0hwByEIQCoxXc3GXOGAGMQhgCgFt8K1HSGAGMQhgCgFvYmA8xDGAKAWlxnO0MV7E0GGIMwBAC1eIfJKirpDAGmIAwBQC01iy7SGQJMQRgCgFpqFl2kMwSYgjAEALX4OkPcTQYYgzAEALV4J1CzNxlgDsIQANTiXXSxnM4QYAzCEADU4mLRRcA4hCEAqMXJoouAcQhDAFBLzaKLDcNQWaUn2OUACALCEADUcq67yf5n2xfq9dg72vTpsVCUBSCACEMAUItvBep6YWjb/kJVWdKrf88LRVkAAogwBAC1eOcM1V908atTFZKkzZ8dY3I10MIQhgCglpp1huoGnq++LpckFZ+pVHbeyWCXBSCACEMAUItvAvU5OkOStGkf84aAloQwBAC1eBddrD1nyLIsnTxV7nu+8dOjQa8LQOAQhgCglpq7yWo6QyVllXW25/jHoWIdKykLem0AAoMwBAC11OxaX9MZ8s4XinFF6qqOcZKkD7jFHmgxCEMAUIt3mKx2J8g7X6htK7eGX9lOklhvCGhBCEMAUIu7kb3Jvjo7X6hNrEvX9WgvSfrgs2PysLM90CKEXRhaunSpkpOTFR0drZSUFG3evPm852/atEkpKSmKjo5Wt27dtHz58iBVCiAcORu5m8w7THZxrFsDLmuj1lFOnTxVoY8PngxFiQD8LKzCUEZGhqZNm6bZs2crOztbw4YN080336zc3NxGzz9w4IBuueUWDRs2TNnZ2frP//xPTZ06VW+88UaQKwcQLmoWXazdGaoeJru4lVuuyAj9vysSJEkb691if6bCo/X/LND+Y6WyLLpGQLhwhroAOxYtWqSJEydq0qRJkqQlS5bo3Xff1bJly7Rw4cIG5y9fvlyXXXaZlixZIknq1auXdu3apaefflo/+tGPglk6gDBRs+hiY50hlyRp+JXt9L//KNCmT4/poRuvlCSVnKnQxP/ZpQ8PnJAktW8dpdTul2hwt0vU9ZJWKj5ToaLTFSo6VaGSMxWqsqQIhySHQxEOKeLs/zocDjnOPncE8ec2jYOLe0Hpn3Sxrk1uG7LvHzZhqLy8XFlZWZo5c2ad42lpadq2bVuj79m+fbvS0tLqHLvpppu0cuVKVVRUyOVyNXhPWVmZyspqbpktLi72Q/UAwkVj6wx55wxdHOuWJA3vUT2J+qODJ3Xi63I5JN296kN9fLBI0a4IVVnS0ZIy/TnnsP6cczi4PwAQhh64rjthqCkKCwvl8XiUmJhY53hiYqIKCgoafU9BQUGj51dWVqqwsFAdOnRo8J6FCxdq3rx5/iscQFipmUBd+26yup2hDvEx6pHYWvuOlGjt7oN6dVeePj1Sqrat3HppwrW6vP1F2p37lXZ8fkI79h9XYWmZ4mNdio9xqU2MS3ExLkU4HLIsS1WWZOns/1qWqqqqn7P9WeBYYgjzQnNVx/iQfv+wCUNejnq9TcuyGhz7pvMbO+41a9Yspaen+54XFxcrKSmpueUCCDONdoa+rpkz5HVdj3bad6RET/x1ryQpMS5KqycN0uXtW0uShnRP0JDuCdKNwaocQHOFzQTqhIQERUZGNugCHT16tEH3x+vSSy9t9Hyn06lLLrmk0fdERUUpLi6uzgOAObyLLlZWWb5/PNUfJpPkW29IkpLaxuj1yUN8QQhAeAmbMOR2u5WSkqLMzMw6xzMzMzVkyJBG35Oamtrg/PXr12vgwIGNzhcCAO92HFLNJOrGwtDArm3Vt1O8+naK12v3D1FS29jgFgrAb8JqmCw9PV3jxo3TwIEDlZqaqhUrVig3N1eTJ0+WVD3EdejQIb300kuSpMmTJ+u5555Tenq67rvvPm3fvl0rV67UmjVrQvljALiAedcZkqrnDTkjrFq31tf8I8rtjNBbU4ZKOvewO4DwEFZhaMyYMTp+/Ljmz5+v/Px89enTR+vWrVOXLl0kSfn5+XXWHEpOTta6dev00EMP6fnnn1fHjh317LPPcls9gHOq3RmqqKpSVbml8srq+UO1O0MSIQhoKcIqDEnSAw88oAceeKDR137/+983ODZ8+HDt3r07wFUBaClcETWdoYrKKhVXeCRVd4Ji3ZGhKgtAAIXNnCEACIaIiOrFD6XqOUO+O8liXXSCgBaKMAQA9dTsT1bV6ORpAC0LYQgA6nGdbQ1VeizCEGAAwhAA1ONyetcaqqrZl6wVy3EALRVhCADq8S68WF5Z67Z6OkNAi0UYAoB6vLfXV1YxZwgwAWEIAOqp2Z+s9oKLhCGgpSIMAUA93rWGKj215gzFMmcIaKkIQwBQjyuyZrNWhsmAlo8wBAD1eIfJyj1VOskwGdDiEYYAoB7voouVHksnGCYDWjzCEADU41108euySp0+uzcZnSGg5SIMAUA93mGyYyVl1c8jHGodFXb7WgNoIsIQANTjnUB9tOSMJKkNm7QCLRphCADqqQlD1Z0h7iQDWjbCEADU44yoO0xGGAJaNsIQANTToDPEJq1Ai0YYAoB66k+gpjMEtGyEIQCox9sZKjrNgouACQhDAFCPd9d6LxZcBFo2whAA1OOMqPursQ3DZECLRhgCgHqc9TpDbQlDQItGGAKAerxzhry4mwxo2QhDAFBPwzlDdIaAlowwBAD11J8zRBgCWjbCEADUU7sz5HBIcTEMkwEtGWEIAOpx1poz1CbGpcgINmkFWjLCEADU46wVfhgiA1o+whAA1ON21vxqZPVpoOUjDAFAPbUnULP6NNDyEYYAoJ7aiy4yTAa0fIQhAKin9t1kDJMBLR9hCADqqT1M1oZhMqDFIwwBQD21t+NgXzKg5SMMAUA9tYfJ2LEeaPkIQwBQT+1FF9syZwho8QhDAFCPq86ii8wZAlo6whAA1ONwMEwGmIQwBAD1nCqv9P2Zu8mAlo8wBAD1nKmo8v259p1lAFomZ6gLAIALzfW92qtvp3h9p2vbUJcCIAgIQwBQT7QrUn/5j/8X6jIABEnY9H+/+uorjRs3TvHx8YqPj9e4ceN08uTJ875n7dq1uummm5SQkCCHw6GcnJyg1AoAAMJH2IShsWPHKicnR++8847eeecd5eTkaNy4ced9z9dff62hQ4fqySefDFKVAAAg3ITFMNnevXv1zjvvaMeOHRo0aJAk6YUXXlBqaqr27dunHj16NPo+b1j64osvglUqAAAIM2HRGdq+fbvi4+N9QUiSBg8erPj4eG3bts2v36usrEzFxcV1HgAAoOUKizBUUFCg9u3bNzjevn17FRQU+PV7LVy40DcvKT4+XklJSX79+gAA4MIS0jA0d+5cORyO8z527dolqe6KsF6WZTV6/NuYNWuWioqKfI+8vDy/fn0AAHBhCemcoSlTpuiuu+467zldu3bVxx9/rCNHjjR47dixY0pMTPRrTVFRUYqKivLr1wQAABeukIahhIQEJSQkfON5qampKioq0ocffqhrr71WkrRz504VFRVpyJAhgS4TAAC0YGExZ6hXr14aOXKk7rvvPu3YsUM7duzQfffdp+9///t17iTr2bOn3nzzTd/zEydOKCcnR3v27JEk7du3Tzk5OX6fZwQAAMJXWIQhSVq9erX69u2rtLQ0paWlqV+/fnr55ZfrnLNv3z4VFRX5nr/11lsaMGCARo0aJUm66667NGDAAC1fvjyotQMAgAuXw7IsK9RFXMiKi4sVHx+voqIixcXFhbocAADQBHY+v8OmMwQAABAIhCEAAGA0whAAADBaWOxNFkreKVVsywEAQPjwfm43ZWo0YegblJSUSBLbcgAAEIZKSkoUHx9/3nO4m+wbVFVV6fDhw2rdurXft/4oLi5WUlKS8vLyuFMtgLjOwcF1Dg6uc3BwnYMjkNfZsiyVlJSoY8eOiog4/6wgOkPfICIiQp07dw7o94iLi+M/tiDgOgcH1zk4uM7BwXUOjkBd52/qCHkxgRoAABiNMAQAAIxGGAqhqKgozZkzR1FRUaEupUXjOgcH1zk4uM7BwXUOjgvlOjOBGgAAGI3OEAAAMBphCAAAGI0wBAAAjEYYAgAARiMMBdDSpUuVnJys6OhopaSkaPPmzec9f9OmTUpJSVF0dLS6deum5cuXB6nS8GfnWq9du1Y33nij2rVrp7i4OKWmpurdd98NYrXhy+7faa+tW7fK6XSqf//+gS2whbB7ncvKyjR79mx16dJFUVFR6t69u1588cUgVRu+7F7n1atX6+qrr1ZsbKw6dOige++9V8ePHw9SteHpgw8+0A9+8AN17NhRDodDf/rTn77xPSH5LLQQEK+88orlcrmsF154wdqzZ4/14IMPWq1atbK+/PLLRs///PPPrdjYWOvBBx+09uzZY73wwguWy+WyXn/99SBXHn7sXusHH3zQeuqpp6wPP/zQ+vTTT61Zs2ZZLpfL2r17d5ArDy92r7PXyZMnrW7dullpaWnW1VdfHZxiw1hzrvPo0aOtQYMGWZmZmdaBAwesnTt3Wlu3bg1i1eHH7nXevHmzFRERYT3zzDPW559/bm3evNm66qqrrFtvvTXIlYeXdevWWbNnz7beeOMNS5L15ptvnvf8UH0WEoYC5Nprr7UmT55c51jPnj2tmTNnNnr+jBkzrJ49e9Y5dv/991uDBw8OWI0thd1r3ZjevXtb8+bN83dpLUpzr/OYMWOsRx991JozZw5hqAnsXuf//d//teLj463jx48Ho7wWw+51/vWvf21169atzrFnn33W6ty5c8BqbGmaEoZC9VnIMFkAlJeXKysrS2lpaXWOp6Wladu2bY2+Z/v27Q3Ov+mmm7Rr1y5VVFQErNZw15xrXV9VVZVKSkrUtm3bQJTYIjT3Oq9atUr79+/XnDlzAl1ii9Cc6/zWW29p4MCB+tWvfqVOnTrpyiuv1PTp03X69OlglByWmnOdhwwZooMHD2rdunWyLEtHjhzR66+/rlGjRgWjZGOE6rOQjVoDoLCwUB6PR4mJiXWOJyYmqqCgoNH3FBQUNHp+ZWWlCgsL1aFDh4DVG86ac63r+81vfqOvv/5ad955ZyBKbBGac50/++wzzZw5U5s3b5bTya+apmjOdf7888+1ZcsWRUdH680331RhYaEeeOABnThxgnlD59Cc6zxkyBCtXr1aY8aM0ZkzZ1RZWanRo0frt7/9bTBKNkaoPgvpDAWQw+Go89yyrAbHvun8xo6jIbvX2mvNmjWaO3euMjIy1L59+0CV12I09Tp7PB6NHTtW8+bN05VXXhms8loMO3+fq6qq5HA4tHr1al177bW65ZZbtGjRIv3+97+nO/QN7FznPXv2aOrUqXr88ceVlZWld955RwcOHNDkyZODUapRQvFZyD/XAiAhIUGRkZEN/oVx9OjRBonX69JLL230fKfTqUsuuSRgtYa75lxrr4yMDE2cOFGvvfaabrjhhkCWGfbsXueSkhLt2rVL2dnZmjJliqTqD23LsuR0OrV+/Xp973vfC0rt4aQ5f587dOigTp06KT4+3nesV69esixLBw8e1BVXXBHQmsNRc67zwoULNXToUP3iF7+QJPXr10+tWrXSsGHD9MQTT9C995NQfRbSGQoAt9utlJQUZWZm1jmemZmpIUOGNPqe1NTUBuevX79eAwcOlMvlClit4a4511qq7gjdc889+uMf/8iYfxPYvc5xcXH65JNPlJOT43tMnjxZPXr0UE5OjgYNGhSs0sNKc/4+Dx06VIcPH1Zpaanv2KeffqqIiAh17tw5oPWGq+Zc51OnTikiou5HZmRkpKSazgW+vZB9FgZ0erbBvLdtrly50tqzZ481bdo0q1WrVtYXX3xhWZZlzZw50xo3bpzvfO/thA899JC1Z88ea+XKldxa30R2r/Uf//hHy+l0Ws8//7yVn5/ve5w8eTJUP0JYsHud6+Nusqaxe51LSkqszp07Wz/+8Y+tf/7zn9amTZusK664wpo0aVKofoSwYPc6r1q1ynI6ndbSpUut/fv3W1u2bLEGDhxoXXvttaH6EcJCSUmJlZ2dbWVnZ1uSrEWLFlnZ2dm+JQwulM9CwlAAPf/881aXLl0st9ttXXPNNdamTZt8r919993W8OHD65y/ceNGa8CAAZbb7ba6du1qLVu2LMgVhy8713r48OGWpAaPu+++O/iFhxm7f6drIww1nd3rvHfvXuuGG26wYmJirM6dO1vp6enWqVOnglx1+LF7nZ999lmrd+/eVkxMjNWhQwfrpz/9qXXw4MEgVx1e3n///fP+vr1QPgsdlkV/DwAAmIs5QwAAwGiEIQAAYDTCEAAAMBphCAAAGI0wBAAAjEYYAgAARiMMAQAAoxGGAACA0QhDAADAaIQhAABgNMIQAAAwGmEIgFGOHTumSy+9VL/85S99x3bu3Cm3263169eHsDIAocJGrQCMs27dOt16663atm2bevbsqQEDBmjUqFFasmRJqEsDEAKEIQBG+vnPf6733ntP3/nOd/TRRx/p73//u6Kjo0NdFoAQIAwBMNLp06fVp08f5eXladeuXerXr1+oSwIQIswZAmCkzz//XIcPH1ZVVZW+/PLLUJcDIIToDAEwTnl5ua699lr1799fPXv21KJFi/TJJ58oMTEx1KUBCAHCEADj/OIXv9Drr7+ujz76SBdddJFGjBih1q1b6+233w51aQBCgGEyAEbZuHGjlixZopdffllxcXGKiIjQyy+/rC1btmjZsmWhLg9ACNAZAgAARqMzBAAAjEYYAgAARiMMAQAAoxGGAACA0QhDAADAaIQhAABgNMIQAAAwGmEIAAAYjTAEAACMRhgCAABGIwwBAACj/X+0LuBTT8cHQAAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "data = np.genfromtxt(\"task13/trapped_concentration_multi.txt\", delimiter=\",\", skip_header=1)\n",
    "\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "sorted_data = data[data[:, 0].argsort()]\n",
    "plt.plot(sorted_data[:, 0], sorted_data[:, 1])\n",
    "\n",
    "plt.xlabel(\"x\")\n",
    "plt.ylabel(\"trapped concentration\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "These oscillations around the interface discontinuity overshoots and undershoots are due to the [Gibbs phenomenon](https://en.wikipedia.org/wiki/Gibbs_phenomenon).\n",
    "\n",
    "This phenomenon occurs in finite element when trying to approximate a discontinuous function with continuous elements (see [ten Eikelder et al](https://doi.org/10.1142/S0218202524500040) for more info).\n",
    "\n",
    "The best way to get rid of this behaviour is to use discontinuous elements.\n",
    "Here, we'll change the ``traps_element_type`` in ``my_model.settings`` to ``\"DG\"`` (Discontinuous Galerkin). "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Defining initial values\n",
      "Defining variational problem\n",
      "Defining source terms\n",
      "Defining boundary conditions\n",
      "Solving steady state problem...\n",
      "Solved problem in 0.00 s\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjcAAAGwCAYAAABVdURTAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAAA3bklEQVR4nO3de3xU9Z3/8fdcSAJIYjUaQQIGUAG5mlRI+FGKlyhS8VZNtQteYFdElkukCIsWYcVoazF4CYqL+MAixgta67KViAIBlEoMSgurLKIJEMSAJsRLYibn9weZgcxMxjlhkmHOeT0fj3nUOXNm8sl50OSdz/leHIZhGAIAALAIZ7QLAAAAiCTCDQAAsBTCDQAAsBTCDQAAsBTCDQAAsBTCDQAAsBTCDQAAsBR3tAtoaw0NDdq/f786deokh8MR7XIAAEAYDMPQkSNH1KVLFzmdoXsztgs3+/fvV2pqarTLAAAALVBeXq6uXbuGPMd24aZTp06Sjl6cxMTEKFcDAADCUV1drdTUVN/v8VBsF268t6ISExMJNwAAxJhwhpQwoBgAAFgK4QYAAFgK4QYAAFgK4QYAAFgK4QYAAFgK4QYAAFgK4QYAAFgK4QYAAFgK4QYAAFgK4QYAAFhK1MNNQUGB0tLSlJCQoPT0dBUXFzd77rp16+RwOAIe//u//9uGFQMAgJNZVMNNYWGhpk2bpjlz5qi0tFTDhw/XqFGjVFZWFvJ9n3zyiSoqKnyPc889t40qBgAAJzuHYRhGtL74kCFDdOGFF2rx4sW+Y3369NE111yjvLy8gPPXrVunkSNH6uuvv9app57aoq9ZXV2tpKQkVVVVRXTjTE+DoYqq7yP2eQDMST4lXgntXL7n39d51D7OFeIdAGKJmd/fUdsVvK6uTiUlJZo1a1aT49nZ2dq8eXPI9w4ePFg//PCD+vbtq3vvvVcjR45s9tza2lrV1tb6nldXV59Y4c049G2t/t/D77bKZwP4acmnxGnd70bqlHi3/rm/Stc+uVn/9osemnH5+dEuDUAbi1q4qayslMfjUUpKSpPjKSkpOnDgQND3dO7cWUuWLFF6erpqa2v1/PPP65JLLtG6dev0i1/8Iuh78vLyNG/evIjXH0y8O+pDmABbqq1vUGVNnfZ+/Z16n5WoTw4cUZ2nQdvKv4l2aQCiIGrhxsvhcDR5bhhGwDGv888/X+eff+yvsMzMTJWXl+uRRx5pNtzMnj1bubm5vufV1dVKTU2NQOVNndkpQZ88MCrinwvgp1204G0dPFIrT0PTu+z+zwHYQ9RaDcnJyXK5XAFdmoMHDwZ0c0IZOnSodu3a1ezr8fHxSkxMbPIAYC0u59E/iAg3AKQohpu4uDilp6erqKioyfGioiJlZWWF/TmlpaXq3LlzpMsDEEOaDTfRmy8BIIqielsqNzdXY8eOVUZGhjIzM7VkyRKVlZVp4sSJko7eUtq3b5+WL18uScrPz9c555yjCy64QHV1dfrzn/+sV199Va+++mo0vw0AUdZcuKmncwPYUlTDTU5Ojg4dOqT58+eroqJC/fr10+rVq9W9e3dJUkVFRZM1b+rq6jRjxgzt27dP7du31wUXXKD//u//1pVXXhmtbwHASaC5cNNAuAFsKeoDiidNmqRJkyYFfe25555r8nzmzJmaOXNmG1QFIJa4HHRuABzD3GUAMc/XuTHo3AAg3ACwAG+48e/U1Dc0RKMcAFFGuAEQ87zhxr9TQ+MGsCfCDYCY1/xsKTo3gB0RbgDEvOYGFJNtAHsi3ACIec0NKGaFYsCeCDcAYh6L+AE4HuEGQMxrdhE/tl8AbIlwAyDmNTsV3MOgG8COCDcAYp6bqeAAjkO4ARDznA4W8QNwDOEGQMxzuxo7NwHbL0SjGgDRRrgBEPN8nRsPnRsAhBsAFuAbc+PfuTEkgxlTgO0QbgDEPGczs6UkFvID7IhwAyDmuZtZ50YKXLUYgPURbgDEvOYW8WvuGABrI9wAiHnOZjbObO4YAGsj3ACIec0NKJYIN4AdEW4AxDwGFAM4HuEGQMxrbvsFiQHFgB0RbgDEPDo3AI5HuAEQ80JOBSfcALZDuAEQ81zMlgJwHMINgJjnch79URZsfA3hBrAfwg2AmOdq/Enm8RBuABBuAFhAyM4Ns6UA2yHcAIh5vs5NkC5NfZBuDgBrI9wAiHm+zk2QcBNs1WIA1ka4ARDzXEcnSwXv3DDmBrAdwg2AmOdyhejcEG4A2yHcAIh53nVugnVp6NwA9kO4ARDzQu0KTucGsB/CDYCY5wyx/QKdG8B+CDcAYl6oqeCscwPYD+EGQMwLNRU82KrFAKyNcAMg5oXcOJPODWA7hBsAMc/lHXPDgGIAItwAsABvuGEqOACJcAPAAnxTwdl+AYAINwAswBmqc8OAYsB2CDcAYl6ozg0DigH7IdwAiHlO3/YLDQGvBZtBBcDaCDcAYp7b5d1+IfA1wg1gP4QbADGPzg2A4xFuAMS8Y2NuAl8j3AD2Q7gBEPOOrXND5wYA4QaABfhWKA7WuWG2FGA7hBsAMc8bboIt2EfnBrAfwg2AmOe7LRWkdUO4AeyHcAMg5nl3BQ+WY9hbCrAfwg2AmBdqQDG7ggP2Q7gBEPNcoaaCM6AYsB3CDYCYx1RwAMcj3ACIecdmS0mGX6eGcAPYD+EGQMzzDiiWAsMM4QawH8INgJjnch0XbujcALYX9XBTUFCgtLQ0JSQkKD09XcXFxWG9b9OmTXK73Ro0aFDrFgjgpBeqc8NUcMB+ohpuCgsLNW3aNM2ZM0elpaUaPny4Ro0apbKyspDvq6qq0rhx43TJJZe0UaUATmbeMTdSYLhhKjhgP1ENNwsXLtT48eM1YcIE9enTR/n5+UpNTdXixYtDvu+OO+7QzTffrMzMzJ/8GrW1taqurm7yAGAtocINnRvAfqIWburq6lRSUqLs7Owmx7Ozs7V58+Zm37ds2TLt3r1bc+fODevr5OXlKSkpyfdITU09oboBnHxC3ZYKtt8UAGuLWriprKyUx+NRSkpKk+MpKSk6cOBA0Pfs2rVLs2bN0ooVK+R2u8P6OrNnz1ZVVZXvUV5efsK1Azi5OJ0OefMNnRsA4SWEVuQ47i8u6egaFf7HJMnj8ejmm2/WvHnzdN5554X9+fHx8YqPjz/hOgGc3NxOh370GAGzpRhzA9hP1MJNcnKyXC5XQJfm4MGDAd0cSTpy5Ii2bt2q0tJSTZ48WZLU0NAgwzDkdru1Zs0aXXzxxW1SO4CTj9PhkGQE6dwE2ZMBgKVF7bZUXFyc0tPTVVRU1OR4UVGRsrKyAs5PTEzU9u3btW3bNt9j4sSJOv/887Vt2zYNGTKkrUoHcBJyNw4qDlzELxrVAIimqN6Wys3N1dixY5WRkaHMzEwtWbJEZWVlmjhxoqSj42X27dun5cuXy+l0ql+/fk3ef+aZZyohISHgOAD7cTYbbkg3gN1ENdzk5OTo0KFDmj9/vioqKtSvXz+tXr1a3bt3lyRVVFT85Jo3ACCF6Nww5AawHYfhv8ucxVVXVyspKUlVVVVKTEyMdjkAIiTjgSJV1tTpb9OGa8f+auW+9JEkaViv07ViwtAoVwfgRJn5/R317RcAIBKcjbMs6z3sLQXYHeEGgCV4b0v5L9rHkBvAfgg3ACzBO6DYf9E+poID9kO4AWAJvs4NA4oB2yPcALCE5jo3TAUH7IdwA8ASmu3ckG0A2yHcALAE32wpOjeA7RFuAFiC29W4iJ/BVHDA7gg3ACzB1di58bDODWB7hBsAluByNtO5sdci7ABEuAFgEa7mBhQzFxywHcINAEtwNTcVnM4NYDuEGwCW4Gpm+wXG3AD2Q7gBYAku59EfZ2ycCYBwA8ASGmeCB9yG8r9NBcD6CDcALMHbufHv1PgPMAZgfYQbAJbgavxp5h9uGFAM2A/hBoAl+Na5Cdh+gXAD2I27JW/69NNPtW7dOh08eFANfvu2/P73v49IYQBgRnO3pQg3gP2YDjfPPPOM7rzzTiUnJ+uss86So3HJc0lyOByEGwBR4RtQ7D/mxpAMw2jyswqAtZkONw888IAWLFige+65pzXqAYAW8XVugoyx8TQYvo01AVif6TE3X3/9tW644YbWqAUAWqy5AcUS08EBuzEdbm644QatWbOmNWoBgBZrbsyNFLhqMQBrM31bqlevXrrvvvv0/vvvq3///mrXrl2T16dMmRKx4gAgXN7OTbAuDZ0bwF5Mh5slS5bolFNO0fr167V+/fomrzkcDsINgKhwN3Zugi3ax0J+gL2YDjd79uxpjToA4IQ4G2dDBRtQTOcGsJcTWsTPMAwZ3MsGcBLwzoYKOuaGcAPYSovCzfLly9W/f3+1b99e7du314ABA/T8889HujYACJuvc8OYG8D2TN+WWrhwoe677z5NnjxZw4YNk2EY2rRpkyZOnKjKykpNnz69NeoEgJDczWy/0NwxANZlOtw8/vjjWrx4scaNG+c7dvXVV+uCCy7Q/fffT7gBEBVOwg2ARqZvS1VUVCgrKyvgeFZWlioqKiJSFACY5e3cBLsFxc7ggL2YDje9evXSSy+9FHC8sLBQ5557bkSKAgCzvLuCBxs8TOcGsBfTt6XmzZunnJwcbdiwQcOGDZPD4dDGjRu1du3aoKEHANqCK1TnhnAD2Irpzs3111+vLVu2KDk5Wa+//rpWrVql5ORk/f3vf9e1117bGjUCwE9yNc6WCrbVAuEGsBfTnRtJSk9P15///OdI1wIALeakcwOgUVjhprq6WomJib7/DsV7HgC0JXeoMTcMKAZsJaxw87Of/UwVFRU688wzdeqpp8rR2P49nmEYcjgc8ng8ES8SAH7Ksc5NQ8BrdG4Aewkr3Lzzzjs67bTTJEnvvvtuqxYEAC1xbBG/wNcIN4C9hBVuRowY4fvvtLQ0paamBnRvDMNQeXl5ZKsDgDC5fNsv0LkB7M70bKm0tDR99dVXAccPHz6stLS0iBQFAGZ5p4J7guQYwg1gL6bDjXdsjb+amholJCREpCgAMItF/AB4hT0VPDc3V5LkcDh03333qUOHDr7XPB6PtmzZokGDBkW8QAAIh4sBxQAahR1uSktLJR3t3Gzfvl1xcXG+1+Li4jRw4EDNmDEj8hUCQBiOdW4CXwu29g0A6wo73HhnSd12221atGgR69kAOKmE6twEW7UYgHWZXqF42bJlrVEHAJwQ32ypIDmGzg1gLy3afuGDDz7Qyy+/rLKyMtXV1TV5bdWqVREpDADMcLmanwoebJAxAOsyPVvqxRdf1LBhw7Rjxw699tpr+vHHH7Vjxw698847SkpKao0aAeAnHVvnJvA1OjeAvZgONw8++KAeffRRvfnmm4qLi9OiRYu0c+dO3XjjjerWrVtr1AgAP+nYCsV0bgC7Mx1udu/erdGjR0uS4uPj9e2338rhcGj69OlasmRJxAsEgHA4feEmMMjQuQHsxXS4Oe2003TkyBFJ0tlnn61//OMfkqRvvvlG3333XWSrA4AwuUOEG3YFB+zF9IDi4cOHq6ioSP3799eNN96oqVOn6p133lFRUZEuueSS1qgRAH6Sr3MTJMhwWwqwF9Ph5oknntAPP/wgSZo9e7batWunjRs36rrrrtN9990X8QIBIBy+AcVB5oJzWwqwF1Phpr6+Xn/96191+eWXS5KcTqdmzpypmTNntkpxABAuF50bAI1Mjblxu9268847VVtb21r1AECLuBhQDKCR6QHFQ4YM8e0zBQAni1ADitl+AbAX02NuJk2apLvvvlt79+5Venq6Onbs2OT1AQMGRKw4AAhXyKngwfZkAGBZpjs3OTk52rNnj6ZMmaJhw4Zp0KBBGjx4sO9/zSooKFBaWpoSEhKUnp6u4uLiZs/duHGjhg0bptNPP13t27dX79699eijj5r+mgCsh6ngALxMd2727NkTsS9eWFioadOmqaCgQMOGDdPTTz+tUaNGaceOHUFXO+7YsaMmT56sAQMGqGPHjtq4caPuuOMOdezYUf/2b/8WsboAxB6no/kBxcFWLQZgXabDzRdffKGsrCy53U3fWl9fr82bN6t79+5hf9bChQs1fvx4TZgwQZKUn5+vt956S4sXL1ZeXl7A+YMHD27SHTrnnHO0atUqFRcXNxtuamtrmwyArq6uDrs+ALHD7QrRuSHbALZi+rbUyJEjdfjw4YDjVVVVGjlyZNifU1dXp5KSEmVnZzc5np2drc2bN4f1GaWlpdq8ebNGjBjR7Dl5eXlKSkryPVJTU8OuEUDsOLZxJp0bwO5MhxvDMORo/CFyvEOHDgUMLg6lsrJSHo9HKSkpTY6npKTowIEDId/btWtXxcfHKyMjQ3fddZev8xPM7NmzVVVV5XuUl5eHXSOA2OGdCt5gSP53pujcAPYS9m2p6667TpLkcDh06623Kj4+3veax+PRxx9/rKysLNMF+Ael5sLT8YqLi1VTU6P3339fs2bNUq9evXTTTTcFPTc+Pr5JrQCsyRtupMBxN3RuAHsJO9wkJSVJOho+OnXqpPbt2/tei4uL09ChQ/Wv//qvYX/h5ORkuVyugC7NwYMHA7o5/tLS0iRJ/fv315dffqn777+/2XADwB6ahBu/W1PMlgLsJexws2zZMklHB/HOmDHD1C2oYOLi4pSenq6ioiJde+21vuNFRUW6+uqrw/4cwzBYMRlA6HDDCsWArZieLTV37tyIffHc3FyNHTtWGRkZyszM1JIlS1RWVqaJEydKOjpeZt++fVq+fLkk6cknn1S3bt3Uu3dvSUfXvXnkkUf07//+7xGrCUBsItwA8DIdbr788kvNmDFDa9eu1cGDB2X439v2eML+rJycHB06dEjz589XRUWF+vXrp9WrV/umk1dUVKisrMx3fkNDg2bPnq09e/bI7XarZ8+eeuihh3THHXeY/TYAWIzLESrctHU1AKLJYfink58watQolZWVafLkyercuXPA4F8zt5Siobq6WklJSaqqqlJiYmK0ywEQIYZhKG32aknSnCv7aMHqnb7XrhnURfm/Mb+COoCTh5nf36Y7Nxs3blRxcbEGDRrU0voAIOIcDoecjqNTwQNmS3FXCrAV0+vcpKamBtyKAoCTgauZ/aWYCg7Yi+lwk5+fr1mzZunzzz9vhXIAoOWaDzf8QQbYienbUjk5Ofruu+/Us2dPdejQQe3atWvyerCtGQCgLTS3BQPhBrAX0+EmPz+/FcoAgBNH5waA1IJwc8stt7RGHQBwwnzhxm9cYD3hBrAV02NuJGn37t269957ddNNN+ngwYOSpL/97W/65z//GdHiAMAMl/PojzT/Tk0DkyAAWzEdbtavX6/+/ftry5YtWrVqlWpqaiRJH3/8cURXLwYAs1yNP9H8w009c8EBWzEdbmbNmqUHHnhARUVFiouL8x0fOXKk3nvvvYgWBwBmuOncAFALws327dubbHTpdcYZZ+jQoUMRKQoAWsLZXOeGMTeArZgON6eeeqoqKioCjpeWlurss8+OSFEA0BLezo1/mGkg3AC2Yjrc3Hzzzbrnnnt04MABORwONTQ0aNOmTZoxY4bGjRvXGjUCQFi8G4P7hxk6N4C9mA43CxYsULdu3XT22WerpqZGffv21S9+8QtlZWXp3nvvbY0aASAs/p0b776+rHMD2IvpdW7atWunFStW6D//8z/14YcfqqGhQYMHD9a5557bGvUBQNicja0b7wBit9OhHz0G4QawGdPhxqtHjx7q0aNHJGsBgBPibgw33s6N0+GQZAQs6gfA2kzflvr1r3+thx56KOD4H//4R91www0RKQoAWsLXuWk41rk5/jkAe2jRIn6jR48OOH7FFVdow4YNESkKAFrC1TjGpr6hQdKxsMOAYsBeTIebmpqaJov3ebVr107V1dURKQoAWuLYIn7e53RuADsyHW769eunwsLCgOMvvvii+vbtG5GiAKAlvIv4eQcUu+jcALZkekDxfffdp+uvv167d+/WxRdfLElau3atVq5cqZdffjniBQJAuPyngrv8Zk8BsAfT4WbMmDF6/fXX9eCDD+qVV15R+/btNWDAAL399tsaMWJEa9QIAGHxH1DsctC5AeyoRVPBR48eHXRQMQBE07Gp4E0HFLPODWAvLV7npq6uTgcPHlRD4w8Rr27dup1wUQDQEk6Ht3Nz9LmbcAPYkulws2vXLt1+++3avHlzk+OGYcjhcMjj8USsOAAwg84NAKkF4ebWW2+V2+3Wm2++qc6dO8vh3bwFAKLMO4DY05hl6NwA9mQ63Gzbtk0lJSXq3bt3a9QDAC3mCzfezo3DG3YIN4CdmF7npm/fvqqsrGyNWgDghBwLN0efuxuXLDYMFvID7MR0uHn44Yc1c+ZMrVu3TocOHVJ1dXWTBwBEi3/nxnXcbXO6N4B9mL4tdemll0qSLrnkkibHGVAMINq8Ycbjt4if91g7V1TKAtDGTIebd999tzXqAIAT5nKFDjcA7MF0uGEVYgAnK5ffAOIm4YbbUoBttGgRv2+++UZLly7Vzp075XA41LdvX91+++1KSkqKdH0AEDbfmBtPkHDjIdwAdmF6QPHWrVvVs2dPPfroozp8+LAqKyu1cOFC9ezZUx9++GFr1AgAYTm2zs3RIONkQDFgS6Y7N9OnT9eYMWP0zDPPyO0++vb6+npNmDBB06ZN04YNGyJeJACEw38quPeYp8FgzA1gI6bDzdatW5sEG0lyu92aOXOmMjIyIlocAJjhPxXce4xwA9iL6dtSiYmJKisrCzheXl6uTp06RaQoAGgJ/6ngzR0DYG2mw01OTo7Gjx+vwsJClZeXa+/evXrxxRc1YcIE3XTTTa1RIwCExRVkL6lgxwBYm+nbUo888ogcDofGjRun+vp6SVK7du1055136qGHHop4gQAQLv8BxccfqyfcALZhOtzExcVp0aJFysvL0+7du2UYhnr16qUOHTq0Rn0AELZQnZsGZksBtmE63FRVVcnj8ei0005T//79fccPHz4st9utxMTEiBYIAOEKFW7qWecGsA3TY25+85vf6MUXXww4/tJLL+k3v/lNRIoCgJZw+7o0x455BxTTuQHsw3S42bJli0aOHBlw/Je//KW2bNkSkaIAoCWOX7TPizE3gP2YDje1tbW+gcTH+/HHH/X9999HpCgAaAm3q/lww2wpwD5Mh5uf//znWrJkScDxp556Sunp6REpCgBaIljnxk24AWzH9IDiBQsW6NJLL9VHH32kSy65RJK0du1affDBB1qzZk3ECwSAcLmdgeHGSbgBbMd052bYsGF67733lJqaqpdeekl//etf1atXL3388ccaPnx4a9QIAGFxBgk3dG4A+zHduZGkQYMGacWKFZGuBQBOSNDOjSNwYT8A1ma6cwMAJytXsM5N4yDjBjo3gG0QbgBYRrBw4+3cMBUcsA/CDQDLcDFbCoAINwAsJNiAYmZLAfZDuAFgGcEGFLsYUAzYTlizpa677rqwP3DVqlUtLgYATkTQqeAub+emoa3LARAlYXVukpKSfI/ExEStXbtWW7du9b1eUlKitWvXKikpqdUKBYCfEnIqONkGsI2wOjfLli3z/fc999yjG2+8UU899ZRcLpckyePxaNKkSUpMTGydKgEgDKEHFJNuALswPebm2Wef1YwZM3zBRpJcLpdyc3P17LPPRrQ4ADAj6FRwJ50bwG5Mh5v6+nrt3Lkz4PjOnTvV0IK/jAoKCpSWlqaEhASlp6eruLi42XNXrVqlyy67TGeccYYSExOVmZmpt956y/TXBGBNQRfxo3MD2I7p7Rduu+023X777fq///s/DR06VJL0/vvv66GHHtJtt91m6rMKCws1bdo0FRQUaNiwYXr66ac1atQo7dixQ926dQs4f8OGDbrsssv04IMP6tRTT9WyZct01VVXacuWLRo8eLDZbwWAxYTu3DBbCrALh2GYmx/Z0NCgRx55RIsWLVJFRYUkqXPnzpo6daruvvvuJrerfsqQIUN04YUXavHixb5jffr00TXXXKO8vLywPuOCCy5QTk6Ofv/734d1fnV1tZKSklRVVcUYIcBiPt77jcY8scn3fPi5yTqtY5z+sm2/7h3dRxOG94hidQBOhJnf36Y7N06nUzNnztTMmTNVXV0tSS0KCXV1dSopKdGsWbOaHM/OztbmzZvD+oyGhgYdOXJEp512WrPn1NbWqra21vfcWzMA6wnWufEOMm5gnRvANlq0iF99fb3efvttrVy5Uo7GHxz79+9XTU1N2J9RWVkpj8ejlJSUJsdTUlJ04MCBsD7jT3/6k7799lvdeOONzZ6Tl5fXZCp7ampq2DUCiC1Bw42TvaUAuzEdbr744gv1799fV199te666y599dVXkqQ//OEPmjFjhukCHH5TNw3DCDgWzMqVK3X//fersLBQZ555ZrPnzZ49W1VVVb5HeXm56RoBxIagKxQ72RUcsBvT4Wbq1KnKyMjQ119/rfbt2/uOX3vttVq7dm3Yn5OcnCyXyxXQpTl48GBAN8dfYWGhxo8fr5deekmXXnppyHPj4+OVmJjY5AHAmpxB/jByMRUcsB3T4Wbjxo269957FRcX1+R49+7dtW/fvrA/Jy4uTunp6SoqKmpyvKioSFlZWc2+b+XKlbr11lv1wgsvaPTo0eaKB2BpbmfgjzQXU8EB2zE9oLihoUEejyfg+N69e9WpUydTn5Wbm6uxY8cqIyNDmZmZWrJkicrKyjRx4kRJR28p7du3T8uXL5d0NNiMGzdOixYt0tChQ31dn/bt27P1AwAFyTbHwg0DigHbMN25ueyyy5Sfn+977nA4VFNTo7lz5+rKK6809Vk5OTnKz8/X/PnzNWjQIG3YsEGrV69W9+7dJUkVFRUqKyvznf/000+rvr5ed911lzp37ux7TJ061ey3AcCCgnZuHAwoBuzG9Do3+/fv18iRI+VyubRr1y5lZGRo165dSk5O1oYNG0IO7j0ZsM4NYF0Hj/ygixYcG/s3/Nxk9e2cqKc3fKZ/HZ6mOaP7RrE6ACeiVde56dKli7Zt26aVK1fqww8/VENDg8aPH6/f/va3TQYYA0BbC7ZxJlPBAfsxHW6ko2Ncbr/9dt1+++2RrgcAWizUgGKmggP20aJw88knn+jxxx/Xzp075XA41Lt3b02ePFm9e/eOdH0AELZQA4rp3AD2YXpA8SuvvKJ+/fqppKREAwcO1IABA/Thhx+qf//+evnll1ujRgAIS6gBxWy/ANiH6c7NzJkzNXv2bM2fP7/J8blz5+qee+7RDTfcELHiAMCMoJ0bV2PnxkO4AezCdOfmwIEDGjduXMDxf/mXfwl7TygAaA2hOjescwPYh+lw88tf/lLFxcUBxzdu3Kjhw4dHpCgAaIkgW0sdt0Ix4QawC9O3pcaMGaN77rlHJSUlGjp0qCTp/fff18svv6x58+bpjTfeaHIuALQVh8Mhl9PRJMgQbgD7MR1uJk2aJEkqKChQQUFB0Nekoz9kgm3TAACtyeVwyKNjQcZNuAFsp0V7SwHAycrldEjH/V3lJNwAtmN6zA0AnMxcfgNv6NwA9tOicLN27Vr96le/Us+ePdWrVy/96le/0ttvvx3p2gDANP9w42S2FGA7psPNE088oSuuuEKdOnXS1KlTNWXKFCUmJurKK6/UE0880Ro1AkDYAjo3Ljo3gN2YHnOTl5enRx99VJMnT/YdmzJlioYNG6YFCxY0OQ4Aba3Zzg3hBrAN052b6upqXXHFFQHHs7OzVV1dHZGiAKCl/HcGZyo4YD+mw82YMWP02muvBRz/y1/+oquuuioiRQFASzGgGIDp21J9+vTRggULtG7dOmVmZko6uojfpk2bdPfdd+uxxx7znTtlypTIVQoAYWBAMQDT4Wbp0qX62c9+ph07dmjHjh2+46eeeqqWLl3qe+5wOAg3ANocA4oBmA43e/bsaY06ACAiGFAMgEX8AFiK/4Bi707hhBvAPkx3biRp7969euONN1RWVqa6uromry1cuDAihQFASwR0bhr/hCPcAPZhOtysXbtWY8aMUVpamj755BP169dPn3/+uQzD0IUXXtgaNQJA2AJnS9G5AezG9G2p2bNn6+6779Y//vEPJSQk6NVXX1V5eblGjBihG264oTVqBICw+Ycbl7dzw2wpwDZMh5udO3fqlltukSS53W59//33OuWUUzR//nw9/PDDES8QAMwIDDdHf8zVewg3gF2YDjcdO3ZUbW2tJKlLly7avXu377XKysrIVQYALRAQbhoHGDfQuQFsw/SYm6FDh2rTpk3q27evRo8erbvvvlvbt2/XqlWrNHTo0NaoEQDC1tz2C/WMuQFsw3S4WbhwoWpqaiRJ999/v2pqalRYWKhevXrp0UcfjXiBAGCGd9E+L2+4aSDcALZhKtx4PB6Vl5drwIABkqQOHTqooKCgVQoDgJZwNrdxJrelANswNebG5XLp8ssv1zfffNNK5QDAiXEHDChuDDcMKAZsw/SA4v79++uzzz5rjVoA4IQ5mxlQTOcGsA/T4WbBggWaMWOG3nzzTVVUVKi6urrJAwCiKaBz42JAMWA3pgcUX3HFFZKkMWPGyHHcvW3DMORwOOTxeCJXHQCY1FznhgHFgH2YDjfvvvtua9QBABHR3JgbOjeAfZgON2lpaUpNTW3StZGOdm7Ky8sjVhgAtERz69xIR7s3/p0dANZjesxNWlqavvrqq4Djhw8fVlpaWkSKAoCWCrgtddxzujeAPZgON96xNf5qamqUkJAQkaIAoKWauy0lsQUDYBdh35bKzc2VJDkcDt13333q0KGD7zWPx6MtW7Zo0KBBES8QAMzw79y46dwAthN2uCktLZV0tHOzfft2xcXF+V6Li4vTwIEDNWPGjMhXCAAm+Hdujl+x2EO4AWwh7HDjnSV12223adGiRUpMTGy1ogCgpfy3Xzg+7BBuAHswPVtq2bJlrVEHAEREQOeGcAPYjukBxQBwMnMFmertDTyEG8AeCDcALCVYuHGyMzhgK4QbAJYSsnPDzuCALRBuAFhKsHDDzuCAvRBuAFiK//YL0nG3pRhzA9gC4QaApbhcDCgG7I5wA8BS6NwAINwAsBSmggMg3ACwlKBTwRlQDNgK4QaApfivUCxJbpe3c9PQ1uUAiALCDQBL8d8VXDpuKjjZBrAFwg0ASwk2oNh7q6qezg1gC4QbAJYSdBG/xmNkG8AeCDcALCVUuKFzA9gD4QaApYTs3DBbCrAFwg0ASwnZuWHjTMAWoh5uCgoKlJaWpoSEBKWnp6u4uLjZcysqKnTzzTfr/PPPl9Pp1LRp09quUAAxIdhUcO8gYzo3gD1ENdwUFhZq2rRpmjNnjkpLSzV8+HCNGjVKZWVlQc+vra3VGWecoTlz5mjgwIFtXC2AWOAMOVuKcAPYQVTDzcKFCzV+/HhNmDBBffr0UX5+vlJTU7V48eKg559zzjlatGiRxo0bp6SkpDauFkAscAfZONPF9guArUQt3NTV1amkpETZ2dlNjmdnZ2vz5s0R+zq1tbWqrq5u8gBgXaE6N4QbwB6iFm4qKyvl8XiUkpLS5HhKSooOHDgQsa+Tl5enpKQk3yM1NTVinw3g5ON2Bv5YI9wA9hL1AcUOv7+yDMMIOHYiZs+eraqqKt+jvLw8Yp8N4OQTJNswoBiwGXe0vnBycrJcLldAl+bgwYMB3ZwTER8fr/j4+Ih9HoCTW6jODQOKAXuIWucmLi5O6enpKioqanK8qKhIWVlZUaoKQKxzBevc+LZfINwAdhC1zo0k5ebmauzYscrIyFBmZqaWLFmisrIyTZw4UdLRW0r79u3T8uXLfe/Ztm2bJKmmpkZfffWVtm3bpri4OPXt2zca3wKAk4yLzg1ge1ENNzk5OTp06JDmz5+viooK9evXT6tXr1b37t0lHV20z3/Nm8GDB/v+u6SkRC+88IK6d++uzz//vC1LB3CSCrUrOAOKAXuIariRpEmTJmnSpElBX3vuuecCjhkMCAQQQqjtFwg3gD1EfbYUAERS0HDT2M3x8McRYAuEGwCWEmxAsXfVYg8bZwK2QLgBYCnBBhQ76dwAtkK4AWApwQYUuxlzA9gK4QaApbiCbJzpJNwAtkK4AWApdG4AEG4AWEqw2VJ0bgB7IdwAsJRQU8FZoRiwB8INAEsJFm68t6XYFRywB8INAEvhthQAwg0AS3GH6NwQbgB7INwAsBRnkNlSdG4AeyHcALAUOjcACDcALMUZbMwN2y8AtkK4AWA5/t0b73OmggP2QLgBYDn+3RvvDKoGwg1gC4QbAJbj37nx7hRO5wawB8INAMvx31/K1fiTjs4NYA+EGwCWE3hbis4NYCeEGwCWE3hb6uj/sv0CYA+EGwCW02znxkO4AeyAcAPAcgI6N6xzA9gK4QaA5fhvweC9LcUKxYA9EG4AWI7bFfy2FOEGsAfCDQDLaW4qOOEGsAfCDQDLcTUzoJhwA9gD4QaA5QSEm8ZODlPBAXsg3ACwnMDODRtnAnZCuAFgOc2FG7ZfAOyBcAPAcujcAPZGuAFgOYGzpRoX8SPcALZAuAFgOf6dGzfhBrAVwg0Ay/EPN062XwBshXADwHICOjcuOjeAnRBuAFhOs50bwg1gC4QbAJbjP6CYMTeAvRBuAFhOc1PBCTeAPRBuAFhOwG0pwg1gK4QbAJbT7FRwZksBtkC4AWA5oQYUGwQcwPIINwAsp7nOjSRxZwqwPsINAMvxny3lPC7cMO4GsD7CDQDL8S7a53tOuAFshXADwHKczWycKTGoGLADwg0Ay3E3s86NJHk8hBvA6gg3ACzH6R9uHHRuADsh3ACwHP/OjdPpkDff1Dc0RKEiAG2JcAPAcvw7N9Kx7g3ZBrA+wg0Ay/Hv3EjHxt3QuQGsj3ADwHL817mRjoUbsg1gfYQbAJYT9LYUnRvANgg3ACwn1G2pBmZLAZZHuAFgOaEGFNezQjFgeYQbAJYTqnPD9guA9RFuAFiO//YLEuEGsBPCDQDLoXMD2FvUw01BQYHS0tKUkJCg9PR0FRcXhzx//fr1Sk9PV0JCgnr06KGnnnqqjSoFECtchBvA1qIabgoLCzVt2jTNmTNHpaWlGj58uEaNGqWysrKg5+/Zs0dXXnmlhg8frtLSUv3Hf/yHpkyZoldffbWNKwdwMnM5A3+0EW4A+3BH84svXLhQ48eP14QJEyRJ+fn5euutt7R48WLl5eUFnP/UU0+pW7duys/PlyT16dNHW7du1SOPPKLrr7++LUsHcBJzBfmzzTtb6kD1D9r79XdtXBFgLy6nQ52T2kft60ct3NTV1amkpESzZs1qcjw7O1ubN28O+p733ntP2dnZTY5dfvnlWrp0qX788Ue1a9cu4D21tbWqra31Pa+uro5A9QBOZqE6N1Nf3NbG1QD2c2aneP19zqVR+/pRCzeVlZXyeDxKSUlpcjwlJUUHDhwI+p4DBw4EPb++vl6VlZXq3LlzwHvy8vI0b968yBUO4KR30TmnKS25o67od5bv2FUDu+iLQ//HIn5AG4hvF90hvVG9LSVJDr8pm4ZhBBz7qfODHfeaPXu2cnNzfc+rq6uVmpra0nIBxIBup3fQuzN+2eTYXSN76a6RvaJTEIA2FbVwk5ycLJfLFdClOXjwYEB3xuuss84Ker7b7dbpp58e9D3x8fGKj4+PTNEAAOCkF7W+UVxcnNLT01VUVNTkeFFRkbKysoK+JzMzM+D8NWvWKCMjI+h4GwAAYD9RvSmWm5ur//qv/9Kzzz6rnTt3avr06SorK9PEiRMlHb2lNG7cON/5EydO1BdffKHc3Fzt3LlTzz77rJYuXaoZM2ZE61sAAAAnmaiOucnJydGhQ4c0f/58VVRUqF+/flq9erW6d+8uSaqoqGiy5k1aWppWr16t6dOn68knn1SXLl302GOPMQ0cAAD4OAzDXlMHqqurlZSUpKqqKiUmJka7HAAAEAYzv7+jvv0CAABAJBFuAACApRBuAACApRBuAACApRBuAACApRBuAACApRBuAACApRBuAACApRBuAACApUR1+4Vo8C7IXF1dHeVKAABAuLy/t8PZWMF24ebIkSOSpNTU1ChXAgAAzDpy5IiSkpJCnmO7vaUaGhq0f/9+derUSQ6HI6KfXV1drdTUVJWXl7NvVSviOrcNrnPb4Dq3Ha5122it62wYho4cOaIuXbrI6Qw9qsZ2nRun06muXbu26tdITEzk/zhtgOvcNrjObYPr3Ha41m2jNa7zT3VsvBhQDAAALIVwAwAALIVwE0Hx8fGaO3eu4uPjo12KpXGd2wbXuW1wndsO17ptnAzX2XYDigEAgLXRuQEAAJZCuAEAAJZCuAEAAJZCuAEAAJZCuDGpoKBAaWlpSkhIUHp6uoqLi0Oev379eqWnpyshIUE9evTQU0891UaVxjYz13nVqlW67LLLdMYZZygxMVGZmZl666232rDa2GX237PXpk2b5Ha7NWjQoNYt0CLMXufa2lrNmTNH3bt3V3x8vHr27Klnn322jaqNXWav84oVKzRw4EB16NBBnTt31m233aZDhw61UbWxacOGDbrqqqvUpUsXORwOvf766z/5nqj8HjQQthdffNFo166d8cwzzxg7duwwpk6danTs2NH44osvgp7/2WefGR06dDCmTp1q7Nixw3jmmWeMdu3aGa+88kobVx5bzF7nqVOnGg8//LDx97//3fj000+N2bNnG+3atTM+/PDDNq48tpi9zl7ffPON0aNHDyM7O9sYOHBg2xQbw1pynceMGWMMGTLEKCoqMvbs2WNs2bLF2LRpUxtWHXvMXufi4mLD6XQaixYtMj777DOjuLjYuOCCC4xrrrmmjSuPLatXrzbmzJljvPrqq4Yk47XXXgt5frR+DxJuTLjooouMiRMnNjnWu3dvY9asWUHPnzlzptG7d+8mx+644w5j6NChrVajFZi9zsH07dvXmDdvXqRLs5SWXuecnBzj3nvvNebOnUu4CYPZ6/w///M/RlJSknHo0KG2KM8yzF7nP/7xj0aPHj2aHHvssceMrl27tlqNVhNOuInW70FuS4Wprq5OJSUlys7ObnI8OztbmzdvDvqe9957L+D8yy+/XFu3btWPP/7YarXGspZcZ38NDQ06cuSITjvttNYo0RJaep2XLVum3bt3a+7cua1doiW05Dq/8cYbysjI0B/+8AedffbZOu+88zRjxgx9//33bVFyTGrJdc7KytLevXu1evVqGYahL7/8Uq+88opGjx7dFiXbRrR+D9pu48yWqqyslMfjUUpKSpPjKSkpOnDgQND3HDhwIOj59fX1qqysVOfOnVut3ljVkuvs709/+pO+/fZb3Xjjja1RoiW05Drv2rVLs2bNUnFxsdxufnSEoyXX+bPPPtPGjRuVkJCg1157TZWVlZo0aZIOHz7MuJtmtOQ6Z2VlacWKFcrJydEPP/yg+vp6jRkzRo8//nhblGwb0fo9SOfGJIfD0eS5YRgBx37q/GDH0ZTZ6+y1cuVK3X///SosLNSZZ57ZWuVZRrjX2ePx6Oabb9a8efN03nnntVV5lmHm33NDQ4McDodWrFihiy66SFdeeaUWLlyo5557ju7NTzBznXfs2KEpU6bo97//vUpKSvS3v/1Ne/bs0cSJE9uiVFuJxu9B/vwKU3JyslwuV8BfAQcPHgxIpV5nnXVW0PPdbrdOP/30Vqs1lrXkOnsVFhZq/Pjxevnll3XppZe2Zpkxz+x1PnLkiLZu3arS0lJNnjxZ0tFfwoZhyO12a82aNbr44ovbpPZY0pJ/z507d9bZZ5+tpKQk37E+ffrIMAzt3btX5557bqvWHItacp3z8vI0bNgw/e53v5MkDRgwQB07dtTw4cP1wAMP0FmPkGj9HqRzE6a4uDilp6erqKioyfGioiJlZWUFfU9mZmbA+WvWrFFGRobatWvXarXGspZcZ+lox+bWW2/VCy+8wD3zMJi9zomJidq+fbu2bdvme0ycOFHnn3++tm3bpiFDhrRV6TGlJf+ehw0bpv3796umpsZ37NNPP5XT6VTXrl1btd5Y1ZLr/N1338npbPor0OVySTrWWcCJi9rvwVYdrmwx3qmGS5cuNXbs2GFMmzbN6Nixo/H5558bhmEYs2bNMsaOHes73zsFbvr06caOHTuMpUuXMhU8DGav8wsvvGC43W7jySefNCoqKnyPb775JlrfQkwwe539MVsqPGav85EjR4yuXbsav/71r41//vOfxvr1641zzz3XmDBhQrS+hZhg9jovW7bMcLvdRkFBgbF7925j48aNRkZGhnHRRRdF61uICUeOHDFKS0uN0tJSQ5KxcOFCo7S01Dfl/mT5PUi4MenJJ580unfvbsTFxRkXXnihsX79et9rt9xyizFixIgm569bt84YPHiwERcXZ5xzzjnG4sWL27ji2GTmOo8YMcKQFPC45ZZb2r7wGGP23/PxCDfhM3udd+7caVx66aVG+/btja5duxq5ubnGd99918ZVxx6z1/mxxx4z+vbta7Rv397o3Lmz8dvf/tbYu3dvG1cdW959992QP29Plt+DDsOg/wYAAKyDMTcAAMBSCDcAAMBSCDcAAMBSCDcAAMBSCDcAAMBSCDcAAMBSCDcAAMBSCDcAAMBSCDcAAMBSCDcAAMBSCDcAAMBSCDcAYt5XX32ls846Sw8++KDv2JYtWxQXF6c1a9ZEsTIA0cDGmQAsYfXq1brmmmu0efNm9e7dW4MHD9bo0aOVn58f7dIAtDHCDQDLuOuuu/T222/r5z//uT766CN98MEHSkhIiHZZANoY4QaAZXz//ffq16+fysvLtXXrVg0YMCDaJQGIAsbcALCMzz77TPv371dDQ4O++OKLaJcDIEro3ACwhLq6Ol100UUaNGiQevfurYULF2r79u1KSUmJdmkA2hjhBoAl/O53v9Mrr7yijz76SKeccopGjhypTp066c0334x2aQDaGLelAMS8devWKT8/X88//7wSExPldDr1/PPPa+PGjVq8eHG0ywPQxujcAAAAS6FzAwAALIVwAwAALIVwAwAALIVwAwAALIVwAwAALIVwAwAALIVwAwAALIVwAwAALIVwAwAALIVwAwAALIVwAwAALOX/A9ZFrSQNrb8NAAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "my_model.settings.traps_element_type = \"DG\"\n",
    "my_model.initialise()\n",
    "my_model.run()\n",
    "\n",
    "\n",
    "data = np.genfromtxt(\"task13/trapped_concentration_multi.txt\", delimiter=\",\", skip_header=1)\n",
    "\n",
    "sorted_data = data[data[:, 0].argsort()]\n",
    "plt.plot(sorted_data[:, 0], sorted_data[:, 1])\n",
    "\n",
    "plt.xlabel(\"x\")\n",
    "plt.ylabel(\"trapped concentration\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we can see the oscillations have disappeared! 🎉\n",
    "\n",
    "**Note:** The **new** weird behaviour at the interface is due to the fact that ``TXTExport`` projects the solution onto a new DG functionspace and this projection isn't perfect. This will not show up in Paraview with ``XDMFExport``"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Traps non homogeneously distributed\n",
    "\n",
    "This phenomenon occurs for all discontinuities, even in single material cases.\n",
    "To illustrate this, let's consider a case with one material on $\\Omega = [0, 1]$ and a trap with a density of 1 for $x \\in [0, 5]$ and 0 otherwise.\n",
    "\n",
    "We do this using ``sympy.Piecewise``."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Defining initial values\n",
      "Defining variational problem\n",
      "Defining source terms\n",
      "Defining boundary conditions\n",
      "Solving steady state problem...\n",
      "Solved problem in 0.00 s\n"
     ]
    }
   ],
   "source": [
    "import festim as F\n",
    "import numpy as np\n",
    "import sympy as sp\n",
    "\n",
    "my_model = F.Simulation()\n",
    "\n",
    "vertices = np.concatenate([\n",
    "    np.linspace(0, 0.5, 50),\n",
    "    np.linspace(0.5, 1, 50),\n",
    "])\n",
    "my_model.mesh = F.MeshFromVertices(vertices)\n",
    "my_model.materials = F.Material(id=1, D_0=1, E_D=0)\n",
    "\n",
    "my_model.traps = F.Trap(\n",
    "    k_0=1,\n",
    "    E_k=0,\n",
    "    p_0=1,\n",
    "    E_p=0,\n",
    "    density=sp.Piecewise((1, F.x < 0.5), (0, True)),\n",
    "    materials=my_model.materials[0],\n",
    ")\n",
    "\n",
    "my_model.T = 300\n",
    "\n",
    "my_model.boundary_conditions = [F.DirichletBC(surfaces=[1, 2], value=1, field=\"solute\")]\n",
    "\n",
    "my_model.settings = F.Settings(\n",
    "    absolute_tolerance=1e-10,\n",
    "    relative_tolerance=1e-10,\n",
    "    transient=False,\n",
    ")\n",
    "\n",
    "my_model.exports = [F.TXTExport(filename=\"task13/trapped_concentration.txt\", field=1)]\n",
    "\n",
    "my_model.initialise()\n",
    "my_model.run()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Again we can plot the trapped concentration profile and observe similar oscilliations."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjcAAAGwCAYAAABVdURTAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAAA7eElEQVR4nO3de3iT9f3/8VeaNG1BWgWkcqilgArISVsRyphDsYpszB2EjQ08wO8rMqZQGcIXFeGr1jlXwM2iOGGXJ+xEdM4v30lFgQKKoxbFwdQBWg5FLIyGcmhpev/+aJMmbYHcMW2aO8/HdeW67Cd3knfvMfLic7QZhmEIAADAImLCXQAAAEAoEW4AAIClEG4AAIClEG4AAIClEG4AAIClEG4AAIClEG4AAIClOMJdQEurqanRgQMH1K5dO9lstnCXAwAAAmAYho4dO6YuXbooJubsfTNRF24OHDiglJSUcJcBAACCsHfvXnXr1u2s10RduGnXrp2k2puTmJgY5moAAEAgXC6XUlJSvN/jZxN14cYzFJWYmEi4AQAgwgQypYQJxQAAwFIINwAAwFIINwAAwFIINwAAwFIINwAAwFIINwAAwFIINwAAwFIINwAAwFIINwAAwFIINwAAwFIINwAAwFIINwAAwFIINwAi1qnTbhmGEe4yALQyYQ83eXl5SktLU3x8vNLT01VYWHjW6ysrKzV37lylpqYqLi5OPXv21LJly1qoWgCtxSHXKWU8/I5m5G8LdykAWhlHOD88Pz9f06dPV15enoYNG6ZnnnlGo0aN0o4dO3TxxRc3+ZqxY8fq66+/1nPPPadevXrp0KFDqq6ubuHKAYTbjlKXKiqr9fG+8nCXAqCVCWu4yc3N1aRJkzR58mRJ0qJFi/T2229ryZIlysnJaXT93//+d61fv167d+9W+/btJUndu3dvyZIBtBLlJ09Lkk67a8JcCYDWJmzDUlVVVSoqKlJWVpZfe1ZWljZv3tzka958801lZGTo8ccfV9euXXXppZdq5syZOnny5Bk/p7KyUi6Xy+8BIPK56sKNu4Y5NwD8ha3npqysTG63W8nJyX7tycnJOnjwYJOv2b17tzZu3Kj4+Hi9/vrrKisr09SpU3XkyJEzzrvJycnR/PnzQ14/gPBynaodjj7tJtwA8Bf2CcU2m83vZ8MwGrV51NTUyGaz6aWXXtLgwYN10003KTc3V3/+85/P2HszZ84clZeXex979+4N+e8AoOWVe3tuGJYC4C9sPTcdO3aU3W5v1Etz6NChRr05Hp07d1bXrl2VlJTkbevTp48Mw9C+fft0ySWXNHpNXFyc4uLiQls8gLDzDEtV03MDoIGw9dw4nU6lp6eroKDAr72goECZmZlNvmbYsGE6cOCAKioqvG2ff/65YmJi1K1bt2atF0Dr4um5qWbODYAGwjoslZ2drT/96U9atmyZdu7cqRkzZqikpERTpkyRVDukNHHiRO/148ePV4cOHXT77bdrx44d2rBhg37zm9/ojjvuUEJCQrh+DQBhUB9uGJYC4C+sS8HHjRunw4cPa8GCBSotLVW/fv20evVqpaamSpJKS0tVUlLivf68885TQUGBfv3rXysjI0MdOnTQ2LFj9fDDD4frVwAQJq5T9NwAaJrNiLK9y10ul5KSklReXq7ExMRwlwMgSMMff1d7j9QuJNj96E2KiWl6IQIAazDz/R321VIAEAzXyfqdyU8zNAXAB+EGQMSpqTG8w1ISG/kB8Ee4ARBxKqqq5TugzkZ+AHwRbgBEnPITp/1+pucGgC/CDYCI41kG7lHN4ZkAfBBuAEQc3/k2EsvBAfgj3ACIOK5GPTeEGwD1CDcAIo7vMnCJXYoB+CPcAIg4DefcMKEYgC/CDYCI0zDcsBQcgC/CDYCI03BCMT03AHwRbgBEnEY9N8y5AeCDcAMg4jRcLUXPDQBfhBsAEafxnBt6bgDUI9wAiDiuU/5Lwem5AeCLcAMg4nh6bhwxNkls4gfAH+EGQMTxhJsL2jolcfwCAH+EGwAR5dRpt6qqa+fYdKgLN25WSwHwQbgBEFE8K6VibFJiQqwkNvED4I9wAyCieDbwS0yIldNe+1cYE4oB+CLcAIgonvk2ifGxstdNKGYpOABfhBsAEcVzInhSQqxi7bXhhp4bAL4INwAiiqfnJinBp+eGcAPAhyPcBQCAGd5hqQSHbLa6nhuGpQD4oOcGQERx+fTceDfxo+cGgA/CDYCI4juh2BFT+1cY4QaAL8INgIjiuxTc03PDhGIAvgg3ACJK/ZybWDnsLAUH0BjhBkBEKW9izg09NwB8EW4ARBTffW7sdXNuOH4BgC/CDYCIUj+h2OGziR/DUgDqEW4ARBTPhGK/TfzouQHgg3ADIGK4awwdO1U7LFU7oZiDMwE0RrgBEDEq6oKN1HATP4alANQj3ACIGJ75Nm2cdsXaY7xLwasZlgLgg3ADIGL47k4siaXgAJpEuAEQMXwnE0vyHr/AqeAAfBFuAEQM3xPBJXmHpVgKDsAX4QZAxPA9EVwSS8EBNIlwAyBi+J4rJUmxMSwFB9AY4QZAq3Ts1GkZhn9oaTihuL7nhmEpAPUINwBanec27tHA+WuUt26XX3ujCcV2VksBaCzs4SYvL09paWmKj49Xenq6CgsLz3jtunXrZLPZGj3+9a9/tWDFAJrT68X79D9v7VCNIf3t4wN+z5WfrN+dWKpfLcU+NwB8hTXc5Ofna/r06Zo7d66Ki4s1fPhwjRo1SiUlJWd93WeffabS0lLv45JLLmmhigE0p/Wff6PfvPqJ9+fPvj6moyeqvD83nFDs3cSP1VIAfDjC+eG5ubmaNGmSJk+eLElatGiR3n77bS1ZskQ5OTlnfF2nTp10/vnnt1CVgXHXGCotPxnQtTabLbDrAvzsAN9OtgDfMfD3C1Arr6+1/+8RyGWt+Z4Yqj024cjxqtrHiSpVVdfoyovPV48Lz/Ne9/Heo7rrxSJV1xgaM7CLPt1frt1lx7X1y/9oZN9kSfVzbur3uWFYCkBjYQs3VVVVKioq0uzZs/3as7KytHnz5rO+9oorrtCpU6fUt29f3X///RoxYsQZr62srFRlZaX3Z5fL9e0KP4PDxyv1nd++1yzvDVjVJZ3O07V9OumKlPM19/VPdaLKre/06qgnbhmoB//6qXaXHdeHXx7xhhvPnJvEeM8+N3Wb+DEsBcBH2MJNWVmZ3G63kpOT/dqTk5N18ODBJl/TuXNnLV26VOnp6aqsrNQLL7yg6667TuvWrdN3v/vdJl+Tk5Oj+fPnh7z+psQ5zj3KF/BfwQFeaAR4oRHw+wV4XYBvGPj7BXghIpI9xqYL2sSqfVunjhyvUoLTrgNHT+mLQxX64lCF97rLuyRqyS+vlNMRo8Fp7fXKP/bqwz1HvM97h6Xa0HMD4MzCOiwlNe7+NgzjjF3il112mS677DLvz0OHDtXevXv1xBNPnDHczJkzR9nZ2d6fXS6XUlJSQlC5v07t4vXZw6NC/r4ITsDhyyKhL5CQG+oAaaa2eIddMTH+/78uP3Fa67/4Ru/u/FrvffaNOrWL059vH6x2dcu8B6e1lyR9ur9cxyur1cZpP/NScObcAPARtnDTsWNH2e32Rr00hw4datSbczZDhgzRiy++eMbn4+LiFBcXF3SdiEwBzxkJeOJQwJ8c6je0rKQ2sRozsIvGDOwiqfE/bLpd0EZdz0/Q/qMnVVxyVOmpF3iHnzxzbmJZCg6gCWFbLeV0OpWenq6CggK/9oKCAmVmZgb8PsXFxercuXOoywPQwpoKpJ7emw/3HPb22thjbGrjtNf9N0vBATQW1mGp7OxsTZgwQRkZGRo6dKiWLl2qkpISTZkyRVLtkNL+/fv1/PPPS6pdTdW9e3ddfvnlqqqq0osvvqjXXntNr732Wjh/DQDNZHBae71evF9b9hzR9+t6eJISYr1ByDPnhqXgAHyFNdyMGzdOhw8f1oIFC1RaWqp+/fpp9erVSk1NlSSVlpb67XlTVVWlmTNnav/+/UpISNDll1+u//3f/9VNN90Url8BQDPy9NwU7z2qb47Vrnr0DElJPvvc0HMDwEfYJxRPnTpVU6dObfK5P//5z34/z5o1S7NmzWqBqgC0Bj06tlXH85wqq6hS4RdlkuqXgUs+OxQz5waAj7AfvwAAZ2Kz2by9NwU7ahcfJPr23LAUHEATCDcAWrWruteGm13fHJfUINzYORUcQGOEGwCtmqfnxsNvzk3dsBQ9NwB8EW4AtGq9L0pUO595Nr7hxu5dLWUEvIEiAOsj3ABo1ewxNu/QlFS/O7FUv4mfRO8NgHqEGwCtnu/QVFM9NxIrpgDUI9wAaPV8w01iQv0QVay9/q8wwg0AD8INgFavX5ckJcTWHrlwxp4bVkwBqEO4AdDqOR0xmpiZqks6nacB3c73tjsYlgLQhLDvUAwAgZgzqo/mjOrj12az2WSPscldYzChGIAXPTcAIpqn94aN/AB4EG4ARDSOYADQEOEGQERz1K2YOs3J4ADqEG4ARDR6bgA0RLgBENHszLkB0ADhBkBE82zkR88NAA/CDYCIVn94Jj03AGoRbgBENEfd4ZnVTCgGUIdwAyCiMaEYQEOEGwARzRFTtxSccAOgDuEGQETzDEu5mXMDoA7hBkBEqz9+gZ4bALUINwAimmdYijk3ADwINwAiGpv4AWiIcAMgotXPuaHnBkAtwg2AiOaZc8M+NwA8CDcAIprnVPBqem4A1CHcAIho9Zv4MecGQC3CDYCI5um5YSk4AA/CDYCIxvELABpyBPOizz//XOvWrdOhQ4dU06Ar+MEHHwxJYQAQCO8mfgxLAahjOtw8++yzuuuuu9SxY0dddNFFstls3udsNhvhBkCL8i4FZ1gKQB3T4ebhhx/WI488ovvuu6856gEAU7yb+DEsBaCO6Tk3//nPf3TLLbc0Ry0AYFr98QsMSwGoZTrc3HLLLVqzZk1z1AIAprGJH4CGTA9L9erVSw888IA++OAD9e/fX7GxsX7P33333SErDgDOhU38ADRkOtwsXbpU5513ntavX6/169f7PWez2Qg3AFpUfc8Nw1IAapkON3v27GmOOgAgKJ7VUvTcAPD4Vpv4GYYhw+AvFADhwyZ+ABoKKtw8//zz6t+/vxISEpSQkKABAwbohRdeCHVtAHBOHL8AoCHTw1K5ubl64IEHNG3aNA0bNkyGYWjTpk2aMmWKysrKNGPGjOaoEwCaxMGZABoyHW7+8Ic/aMmSJZo4caK37Yc//KEuv/xyPfTQQ4QbAC3KwSZ+ABowPSxVWlqqzMzMRu2ZmZkqLS01XUBeXp7S0tIUHx+v9PR0FRYWBvS6TZs2yeFwaNCgQaY/E4B12OuGpTh+AYCH6XDTq1cv/eUvf2nUnp+fr0suucTUe+Xn52v69OmaO3euiouLNXz4cI0aNUolJSVnfV15ebkmTpyo6667ztTnAbAe71JwhqUA1DE9LDV//nyNGzdOGzZs0LBhw2Sz2bRx40atXbu2ydBzNrm5uZo0aZImT54sSVq0aJHefvttLVmyRDk5OWd83Z133qnx48fLbrfrjTfeMPsrALCQ+nBDzw2AWqZ7bn7yk59oy5Yt6tixo9544w2tWrVKHTt21Icffqgf/ehHAb9PVVWVioqKlJWV5deelZWlzZs3n/F1y5cv165duzRv3ryAPqeyslIul8vvAcA6vPvcMCwFoI7pnhtJSk9P14svvvitPrisrExut1vJycl+7cnJyTp48GCTr/niiy80e/ZsFRYWyuEIrPScnBzNnz//W9UKoPXyHJzJsBQAj4B6bnx7Oxr2gnzbXhGbzeb3s2EYjdokye12a/z48Zo/f74uvfTSgN9/zpw5Ki8v9z727t1rukYArReb+AFoKKDujwsuuEClpaXq1KmTzj///CbDhyeUuN3ugD64Y8eOstvtjXppDh061Kg3R5KOHTumrVu3qri4WNOmTZMk1dTUyDAMORwOrVmzRtdee22j18XFxSkuLi6gmgBEHjbxA9BQQOHm3XffVfv27SVJ7733Xkg+2Ol0Kj09XQUFBX5zdQoKCvTDH/6w0fWJiYnavn27X1teXp7effddrVy5UmlpaSGpC0BkoecGQEMBhZtrrrnG+99paWlKSUlpcjjJ7JBPdna2JkyYoIyMDA0dOlRLly5VSUmJpkyZIql2SGn//v16/vnnFRMTo379+vm9vlOnToqPj2/UDiB6eCYUn+ZUcAB1TE8oTktL8w5R+Tpy5IjS0tICHpaSpHHjxunw4cNasGCBSktL1a9fP61evVqpqamSajcMPNeeNwCim52eGwAN2AyTx3rHxMTo66+/1oUXXujX/tVXX6lv3746fvx4SAsMNZfLpaSkJJWXlysxMTHc5QD4lj7cc0Rjn3lfaR3b6r2Z3wt3OQCaiZnv74B7brKzsyXVrm564IEH1KZNG+9zbrdbW7Zs4SgEAC3Ou88NS8EB1Ak43BQXF0uqnVuzfft2OZ1O73NOp1MDBw7UzJkzQ18hAJyFd4diVksBqBNwuPGskrr99tu1ePFihnQAtAr1m/gRbgDUMj2hePny5c1RBwAExTMsxYRiAB5BHb/wj3/8Q6+++qpKSkpUVVXl99yqVatCUhgABMIzLMVScAAepg/OfOWVVzRs2DDt2LFDr7/+uk6fPq0dO3bo3XffVVJSUnPUCABn5BmWoucGgIfpcPPoo49q4cKFeuutt+R0OrV48WLt3LlTY8eO1cUXX9wcNQLAGXEqOICGTIebXbt2afTo0ZJqz206fvy4bDabZsyYoaVLl4a8QAA4G+9qKZaCA6hjOty0b99ex44dkyR17dpVn376qSTp6NGjOnHiRGirA4Bz8OxQXGNINQxNAVAQE4qHDx+ugoIC9e/fX2PHjtU999yjd999VwUFBbruuuuao0YAOCPPqeBS7XJwZ4ztLFcDiAamw80f//hHnTp1SlLtwZaxsbHauHGjfvzjH+uBBx4IeYEAcDYOnzBTXVMjp/kOaQAWYyrcVFdX629/+5tuuOEGSbXnTM2aNUuzZs1qluIA4Fw8E4olNvIDUMvUP3EcDofuuusuVVZWNlc9AGCKZym4JLlZMQVAQUwovvrqq73nTAFAuNljbLLVdd6cZsUUAAUx52bq1Km69957tW/fPqWnp6tt27Z+zw8YMCBkxQFAIBwxNp12G2zkB0BSEOFm3LhxkqS7777b22az2WQYhmw2m9xud+iqA4AAOGJidNrtZiM/AJKCCDd79uxpjjoAIGj1G/kRbgAEEW6++uorZWZmyuHwf2l1dbU2b96s1NTUkBUHAIGoP4KBOTcAgphQPGLECB05cqRRe3l5uUaMGBGSogDADHvdiil6bgBIQYQbz9yahg4fPtxocjEAtATvsBRzbgDIxLDUj3/8Y0m1k4dvu+02xcXFeZ9zu9365JNPlJmZGfoKAeAcvMNSLAUHIBPhJikpSVJtz027du2UkJDgfc7pdGrIkCH6f//v/4W+QgA4B0/PDUvBAUgmws3y5cslSd27d9fMmTMZggLQangOzzzNsBQABbFaat68ec1RBwAEjZ4bAL5MTyj++uuvNWHCBHXp0kUOh0N2u93vAQAtzTPnhuMXAEhB9NzcdtttKikp0QMPPKDOnTs3uXIKAFqSZyk4B2cCkIIINxs3blRhYaEGDRrUDOUAgHmxMayWAlDP9LBUSkqKDIN/HQFoPewcvwDAh+lws2jRIs2ePVtffvllM5QDAObVH79AuAEQ5KngJ06cUM+ePdWmTRvFxsb6Pd/U0QwA0JwcHL8AwIfpcLNo0aJmKAMAgle/FJw5NwCCCDe33nprc9QBAEHzLgVnWAqAgphzI0m7du3S/fffr5///Oc6dOiQJOnvf/+7/vnPf4a0OAAIhGdYik38AEhBhJv169erf//+2rJli1atWqWKigpJ0ieffMLuxQDCor7nhmEpAEGEm9mzZ+vhhx9WQUGBnE6nt33EiBF6//33Q1ocAATCzvELAHyYDjfbt2/Xj370o0btF154oQ4fPhySogDAjFhWSwHwYTrcnH/++SotLW3UXlxcrK5du4akKAAww84+NwB8mA4348eP13333aeDBw/KZrOppqZGmzZt0syZMzVx4sTmqBEAzsrB8QsAfJgON4888oguvvhide3aVRUVFerbt6+++93vKjMzU/fff39z1AgAZ8UmfgB8md7nJjY2Vi+99JL+53/+Rx999JFqamp0xRVX6JJLLmmO+gDgnDyrpZhQDEAKItx49OjRQz169AhlLQAQFM+wFEvBAUhBDEv99Kc/1WOPPdao/Xe/+51uueWWkBQFAGY4WAoOwEdQm/iNHj26UfuNN96oDRs2mC4gLy9PaWlpio+PV3p6ugoLC8947caNGzVs2DB16NBBCQkJ6t27txYuXGj6MwFYi8Ne+1cZxy8AkIIYlqqoqPDbvM8jNjZWLpfL1Hvl5+dr+vTpysvL07Bhw/TMM89o1KhR2rFjhy6++OJG17dt21bTpk3TgAED1LZtW23cuFF33nmn2rZtq//6r/8y+6sAsAg7B2cC8GG656Zfv37Kz89v1P7KK6+ob9++pt4rNzdXkyZN0uTJk9WnTx8tWrRIKSkpWrJkSZPXX3HFFfr5z3+uyy+/XN27d9cvf/lL3XDDDWft7QFgfbHscwPAh+memwceeEA/+clPtGvXLl177bWSpLVr12rFihV69dVXA36fqqoqFRUVafbs2X7tWVlZ2rx5c0DvUVxcrM2bN+vhhx8+4zWVlZWqrKz0/my2dwlA62dnKTgAH6bDzZgxY/TGG2/o0Ucf1cqVK5WQkKABAwbonXfe0TXXXBPw+5SVlcntdis5OdmvPTk5WQcPHjzra7t166ZvvvlG1dXVeuihhzR58uQzXpuTk6P58+cHXBeAyOPtuWFYCoCCXAo+evToJicVB8Nms/n9bBhGo7aGCgsLVVFRoQ8++ECzZ89Wr1699POf/7zJa+fMmaPs7Gzvzy6XSykpKd++cACthmfODcNSAKRvsc9NVVWVDh06pJoG/1JqaiJwUzp27Ci73d6ol+bQoUONenMaSktLkyT1799fX3/9tR566KEzhpu4uDjFxcUFVBOAyMRScAC+TE8o/uKLLzR8+HAlJCQoNTVVaWlpSktLU/fu3b2hIxBOp1Pp6ekqKCjway8oKFBmZmbA72MYht+cGgDRx3P8wmnCDQAF0XNz2223yeFw6K233lLnzp3POYR0NtnZ2ZowYYIyMjI0dOhQLV26VCUlJZoyZYqk2iGl/fv36/nnn5ckPfXUU7r44ovVu3dvSbX73jzxxBP69a9/HXQNACJf/fELzLkBEES42bZtm4qKirwB49sYN26cDh8+rAULFqi0tFT9+vXT6tWrlZqaKkkqLS1VSUmJ9/qamhrNmTNHe/bskcPhUM+ePfXYY4/pzjvv/Na1AIhc3p4b5twAkGQzDMPU3wZXXXWVFi5cqO985zvNVVOzcrlcSkpKUnl5uRITE8NdDoAQ+PunBzXlxSKlp16g1+4KfFgbQOQw8/1tes7Nb3/7W82aNUvr1q3T4cOH5XK5/B4A0NLqN/FjWApAEMNSI0eOlCRdd911fu2eJdxutzs0lQFAgLxLwZlQDEBBhJv33nuvOeoAgKDF1h2cyT43AKQgwo2ZXYgBoCXU99wwLAUgyE38jh49queee047d+6UzWZT3759dccddygpKSnU9QHAObGJHwBfpicUb926VT179tTChQt15MgRlZWVKTc3Vz179tRHH33UHDUCwFk57CwFB1DPdM/NjBkzNGbMGD377LNyOGpfXl1drcmTJ2v69OnasGFDyIsEgLOh5waAL9PhZuvWrX7BRpIcDodmzZqljIyMkBYHAIFwcCo4AB+mh6USExP9dg322Lt3r9q1axeSogDADAdLwQH4MB1uxo0bp0mTJik/P1979+7Vvn379Morr2jy5MlnPJkbAJqT5/gFloIDkIIYlnriiSdks9k0ceJEVVdXS5JiY2N111136bHHHgt5gQBwLiwFB+DLdLhxOp1avHixcnJytGvXLhmGoV69eqlNmzbNUR8AnBOb+AHwZTrclJeXy+12q3379urfv7+3/ciRI3I4HBxGCaDF+R6/4DkKBkD0Mj3n5mc/+5leeeWVRu1/+ctf9LOf/SwkRQGAGZ4JxRLLwQEEEW62bNmiESNGNGr/3ve+py1btoSkKAAww7MUXGLFFIAgwk1lZaV3IrGv06dP6+TJkyEpCgDM8KyWkui5ARBEuLnqqqu0dOnSRu1PP/200tPTQ1IUAJjh13PDpGIg6pmeUPzII49o5MiR+vjjj3XddddJktauXat//OMfWrNmTcgLBIBz8Z1zw3JwAKZ7boYNG6b3339fKSkp+stf/qK//e1v6tWrlz755BMNHz68OWoEgLOy2Wx+K6YARDfTPTeSNGjQIL300kuhrgUAgmaPscldYxBuAJjvuQGA1ijW03PjZlgKiHaEGwCWwLAUAA/CDQBL4AgGAB6EGwCWwOGZADwINwAswbMcnE38AAS0WurHP/5xwG+4atWqoIsBgGA56oalTjMsBUS9gHpukpKSvI/ExEStXbtWW7du9T5fVFSktWvXKikpqdkKBYCzoecGgEdAPTfLly/3/vd9992nsWPH6umnn5bdbpckud1uTZ06VYmJic1TJQCcg+cIBpaCAzA952bZsmWaOXOmN9hIkt1uV3Z2tpYtWxbS4gAgUPa6wzNZCg7AdLiprq7Wzp07G7Xv3LlTNaxSABAmsXZWSwGoZfr4hdtvv1133HGH/v3vf2vIkCGSpA8++ECPPfaYbr/99pAXCACB8C4FZ0IxEPVMh5snnnhCF110kRYuXKjS0lJJUufOnTVr1izde++9IS8QAAIRy7AUgDqmw01MTIxmzZqlWbNmyeVySRITiQGEHccvAPAIahO/6upqvfPOO1qxYoVsttq/UA4cOKCKioqQFgcAgfKslnIz5waIeqZ7br766ivdeOONKikpUWVlpa6//nq1a9dOjz/+uE6dOqWnn366OeoEgLPy7HPDJn4ATPfc3HPPPcrIyNB//vMfJSQkeNt/9KMfae3atSEtDgAC5VkKziZ+AEz33GzcuFGbNm2S0+n0a09NTdX+/ftDVhgAmBHLJn4A6pjuuampqZHb7W7Uvm/fPrVr1y4kRQGAWUwoBuBhOtxcf/31WrRokfdnm82miooKzZs3TzfddFMoawOAgMXWHZzJPjcATA9LLVy4UCNGjFDfvn116tQpjR8/Xl988YU6duyoFStWNEeNAHBO9NwA8DAdbrp06aJt27ZpxYoV+uijj1RTU6NJkybpF7/4hd8EYwBoScy5AeBhOtxIUkJCgu644w7dcccdoa4HAIJCzw0Aj6A28fvss880bdo0XXfddRo5cqSmTZumf/3rX0EVkJeXp7S0NMXHxys9PV2FhYVnvHbVqlW6/vrrdeGFFyoxMVFDhw7V22+/HdTnArAWB0vBAdQxHW5Wrlypfv36qaioSAMHDtSAAQP00UcfqX///nr11VdNvVd+fr6mT5+uuXPnqri4WMOHD9eoUaNUUlLS5PUbNmzQ9ddfr9WrV6uoqEgjRozQD37wAxUXF5v9NQBYjHcTP3YoBqKezTAMU//M6dGjh375y19qwYIFfu3z5s3TCy+8oN27dwf8XldffbWuvPJKLVmyxNvWp08f3XzzzcrJyQnoPS6//HKNGzdODz74YEDXu1wuJSUlqby8nDOxAAvJ+b+demb9bk3+Tpru/37fcJcDIMTMfH+b7rk5ePCgJk6c2Kj9l7/8pQ4ePBjw+1RVVamoqEhZWVl+7VlZWdq8eXNA71FTU6Njx46pffv2Z7ymsrJSLpfL7wHAejgVHICH6XDzve99r8l5MRs3btTw4cMDfp+ysjK53W4lJyf7tScnJwcckn7/+9/r+PHjGjt27BmvycnJUVJSkveRkpIScI0AIkf9hGKGpYBoZ3q11JgxY3TfffepqKhIQ4YMkSR98MEHevXVVzV//ny9+eabfteei+dUcQ/DMBq1NWXFihV66KGH9Ne//lWdOnU643Vz5sxRdna292eXy0XAASyofik4PTdAtDMdbqZOnSqpdpVTXl5ek89JtaGlqWMaPDp27Ci73d6ol+bQoUONenMays/P16RJk/Tqq69q5MiRZ702Li5OcXFxZ70GQOSzMywFoE5QZ0sF8jhbsJEkp9Op9PR0FRQU+LUXFBQoMzPzjK9bsWKFbrvtNr388ssaPXq02fIBWBSb+AHwCGoTv1DJzs7WhAkTlJGRoaFDh2rp0qUqKSnRlClTJNUOKe3fv1/PP/+8pNpgM3HiRC1evFhDhgzx9vokJCQoKSkpbL8HgPBjEz8AHkFt4rd27Vp9//vfV8+ePdWrVy99//vf1zvvvGP6fcaNG6dFixZpwYIFGjRokDZs2KDVq1crNTVVklRaWuq3580zzzyj6upq/epXv1Lnzp29j3vuuSeYXwOAhTjsbOIHoJbpfW7++Mc/asaMGfrpT3+qoUOHSqqdULxy5Url5uZq2rRpzVJoqLDPDWBNKz4s0ZxV2zWyT7L+dGtGuMsBEGJmvr9ND0vl5ORo4cKFfiHm7rvv1rBhw/TII4+0+nADwJo8w1JuloIDUc/0sJTL5dKNN97YqD0rK4sN8gCEjXdCMcNSQNQzHW7GjBmj119/vVH7X//6V/3gBz8ISVEAYJZ3KTj73ABRz/SwVJ8+ffTII49o3bp1fnNuNm3apHvvvVdPPvmk99q77747dJUCwFnEskMxgDqmJxSnpaUF9sY2m6lDNFsKE4oBa1rzz4P6rxeKdMXF5+v1qcPCXQ6AEGvWCcV79uwJujAAaC6xdoalANQKap8bAGht2MQPgEdQOxTv27dPb775pkpKSlRVVeX3XG5ubkgKAwAzHHaWggOoZTrcrF27VmPGjFFaWpo+++wz9evXT19++aUMw9CVV17ZHDUCwDk5WC0FoI7pYak5c+bo3nvv1aeffqr4+Hi99tpr2rt3r6655hrdcsstzVEjAJwTw1IAPEyHm507d+rWW2+VJDkcDp08eVLnnXeeFixYoN/+9rchLxAAAsGp4AA8TIebtm3bqrKyUpLUpUsX7dq1y/tcWVlZ6CoDABPouQHgYXrOzZAhQ7Rp0yb17dtXo0eP1r333qvt27dr1apVGjJkSHPUCADn5F0KTrgBop7pcJObm6uKigpJ0kMPPaSKigrl5+erV69eWrhwYcgLBIBAeHtuGJYCop6pcON2u7V3714NGDBAktSmTRvl5eU1S2EAYEZsDD03AGqZmnNjt9t1ww036OjRo81UDgAEx86p4ADqmJ5Q3L9//1Z5ZhSA6OY5ONNNuAGinulw88gjj2jmzJl66623VFpaKpfL5fcAgHCw+4Qbk+cBA7AY0xOKb7zxRknSmDFjZLPZvO2GYchms8ntdoeuOgAIkGeHYql2aMqz7w2A6GM63Lz33nvNUQcAfCsOnzBT7TYUaw9jMQDCynS4SUtLU0pKil+vjVTbc7N3796QFQYAZniGpSSpuqZGEukGiFam59ykpaXpm2++adR+5MgRpaWlhaQoADDLs4mfxOGZQLQzHW48c2saqqioUHx8fEiKAgCzfDpuWA4ORLmAh6Wys7MlSTabTQ888IDatGnjfc7tdmvLli0aNGhQyAsEgEDYbDbF2m067TbqhqUARKuAw01xcbGk2p6b7du3y+l0ep9zOp0aOHCgZs6cGfoKASBA9pi6cMOwFBDVAg43nlVSt99+uxYvXqzExMRmKwoAghEbE6NTqmEjPyDKmV4ttXz58uaoAwC+tfojGBiWAqKZ6QnFANBaOTg8E4AINwAsxFG3ZIo5N0B0I9wAsAzPRn703ADRjXADwDI850lVu5lzA0Qzwg0Ay6DnBoBEuAFgIZ4jGJhzA0Q3wg0Ay6jvuWFYCohmhBsAluGo67lhEz8guhFuAFiGZyn4aYalgKhGuAFgGZ5wQ88NEN0INwAsw8HxCwBEuAFgIfYYVksBINwAsJBYVksBEOEGgIWwiR8AiXADwELYxA+A1ArCTV5entLS0hQfH6/09HQVFhae8drS0lKNHz9el112mWJiYjR9+vSWKxRAq2f3LgVnWAqIZmENN/n5+Zo+fbrmzp2r4uJiDR8+XKNGjVJJSUmT11dWVurCCy/U3LlzNXDgwBauFkBr18ZplySdrHKHuRIA4RTWcJObm6tJkyZp8uTJ6tOnjxYtWqSUlBQtWbKkyeu7d++uxYsXa+LEiUpKSmrhagG0dkkJsZKk8pOnw1wJgHAKW7ipqqpSUVGRsrKy/NqzsrK0efPmkH1OZWWlXC6X3wOANSXWhRvXKcINEM3CFm7KysrkdruVnJzs156cnKyDBw+G7HNycnKUlJTkfaSkpITsvQG0Lon03ABQK5hQbLPZ/H42DKNR27cxZ84clZeXex979+4N2XsDaF0S4x2SCDdAtHOE64M7duwou93eqJfm0KFDjXpzvo24uDjFxcWF7P0AtF6eOTeuk9VhrgRAOIWt58bpdCo9PV0FBQV+7QUFBcrMzAxTVQAiGROKAUhh7LmRpOzsbE2YMEEZGRkaOnSoli5dqpKSEk2ZMkVS7ZDS/v379fzzz3tfs23bNklSRUWFvvnmG23btk1Op1N9+/YNx68AoBVhQjEAKczhZty4cTp8+LAWLFig0tJS9evXT6tXr1Zqaqqk2k37Gu55c8UVV3j/u6ioSC+//LJSU1P15ZdftmTpAFohT8/NsVPVctcY3k39AEQXm2EYUbVPucvlUlJSksrLy5WYmBjucgCEUFV1jS69//8kSR8/mKWkNrFhrghAqJj5/g77aikACBWnI0YJsbW7FDPvBohehBsAlsKkYgCEGwCWkphQO5WQScVA9CLcALAUem4AEG4AWEpivGcjP8INEK0INwAshZ4bAIQbAJbC4ZkACDcALIVdigEQbgBYSv2wFIdnAtGKcAPAUhLj65aCMywFRC3CDQBLYUIxAMINAEthzg0Awg0AS/H03DAsBUQvwg0AS/EdljIMI8zVAAgHwg0AS/EMS512Gzp1uibM1QAIB8INAEtp67TLHmOTxKRiIFoRbgBYis1mq18OzqRiICoRbgBYDsvBgehGuAFgOd5wc4JwA0Qjwg0Ay2GvGyC6EW4AWA4ngwPRjXADwHIS4z0b+XF4JhCNCDcALIcJxUB0I9wAsJzEhNql4IQbIDoRbgBYThITioGoRrgBYDkMSwHRjXADwHLqJxQTboBoRLgBYDneYSnCDRCVCDcALKd+Ez+WggPRiHADwHI8PTcVldWqdteEuRoALY1wA8ByPKeCS/TeANGIcAPAchz2GLV12iUx7waIRoQbAJbEcnAgehFuAFgSJ4MD0YtwA8CSOBkciF6EGwCW5NnIj3ADRB/CDQBLqt/Ij9VSQLQh3ACwJCYUA9GLcAPAkhITave6YUIxEH0INwAsiZ4bIHoRbgBY0tlOBqc3B7A2wg0ASzrTyeDrPjukAQ+t0b1/+Vg1NUY4SgPQzMIebvLy8pSWlqb4+Hilp6ersLDwrNevX79e6enpio+PV48ePfT000+3UKUAIklSm6aHpVZvL5UkvfbRPj2yeqcMg4ADWE1Yw01+fr6mT5+uuXPnqri4WMOHD9eoUaNUUlLS5PV79uzRTTfdpOHDh6u4uFj//d//rbvvvluvvfZaC1cOoLXzDks1ODjzwz1HvP/93MY9Wrphd4vWBaD52Yww/rPl6quv1pVXXqklS5Z42/r06aObb75ZOTk5ja6/77779Oabb2rnzp3etilTpujjjz/W+++/H9BnulwuJSUlqby8XImJid/+lwDQKh0sP6UhOWtlj7Hp34+Mks1m09euU7r60bWKsUnTRvTSk+/+W5KUO3agfnxlN+9ryyoqtenfZXLXGLqgrVMd2jp1QRunEhNiVVNj6HRNjU67DVW7a9TUyJatiXpstobXNL6q4TVApLLH2NTl/ISQvqeZ729HSD/ZhKqqKhUVFWn27Nl+7VlZWdq8eXOTr3n//feVlZXl13bDDTfoueee0+nTpxUbG9voNZWVlaqsrPT+7HK5QlA9gNbOsxTcXWPoRJVbbeMc3l6bvl0SlZ11mU5UufWnjXs0a+UnOu2u0TfHKrX2X4e0be9RMVoFBK9Tuzh9OHdk2D4/bOGmrKxMbrdbycnJfu3Jyck6ePBgk685ePBgk9dXV1errKxMnTt3bvSanJwczZ8/P3SFA4gICbF2xdptOu02VH7ytF+4Gdy9gyTpv2/qo7KKSr2x7YDue2273+v7dk5Uh/OcOlxRpf+cqNLh41Wqqq6RVNvDEmuPUWyMTTExDbpbmghFDZua6jBvKksRsBCp4mPtYf38sIUbD1uDfljDMBq1nev6pto95syZo+zsbO/PLpdLKSkpwZYLIELYbDYlJcSqrKJK5SdPq8v5CfXhJu0CSVJMjE2P/3Sgjle5tfGLMg3r1UHX9k7WiN4XqnOSf5e6YRiqctfIERMje8NAA6BVCVu46dixo+x2e6NemkOHDjXqnfG46KKLmrze4XCoQ4cOTb4mLi5OcXFxoSkaQERJjK8NN66Tp/Wf41X67OtjkqSrurf3XuN0xOjZiRnnfC+bzaY4R3j/NQogMGFbLeV0OpWenq6CggK/9oKCAmVmZjb5mqFDhza6fs2aNcrIyGhyvg2A6Jbos0vx1q/+I0nq1ek8dTiPf/AAVhbWpeDZ2dn605/+pGXLlmnnzp2aMWOGSkpKNGXKFEm1Q0oTJ070Xj9lyhR99dVXys7O1s6dO7Vs2TI999xzmjlzZrh+BQCtmCfcuE5V68M9hyVJg9Pan+0lACwgrHNuxo0bp8OHD2vBggUqLS1Vv379tHr1aqWmpkqSSktL/fa8SUtL0+rVqzVjxgw99dRT6tKli5588kn95Cc/CdevAKAV8z1fyjPf5mrCDWB5Yd3nJhzY5waIHnNf366XtpRo8nfStHzzl3LXGNo8+9qQ778BoPmZ+f4O+/ELANBcPD037312SO4aQyntEwg2QBQg3ACwLE+42fXNcUn1+9sAsDbCDQDL8kwo9mC+DRAdCDcALCupQbhhpRQQHQg3ACzLczK4JF3YLk6pHdqEsRoALYVwA8CyfHtuBqe1P+vRLgCsg3ADwLJ8ww3zbYDoQbgBYFmJCfX7lDLfBogeYT8VHACaS1JCrEZcdqEMSZd2ahfucgC0EMINAMuy2WxafvvgcJcBoIUxLAUAACyFcAMAACyFcAMAACyFcAMAACyFcAMAACyFcAMAACyFcAMAACyFcAMAACyFcAMAACyFcAMAACyFcAMAACyFcAMAACyFcAMAACyFcAMAACzFEe4CWpphGJIkl8sV5koAAECgPN/bnu/xs4m6cHPs2DFJUkpKSpgrAQAAZh07dkxJSUlnvcZmBBKBLKSmpkYHDhxQu3btZLPZQvreLpdLKSkp2rt3rxITE0P63qjHfW4Z3OeWwX1uOdzrltFc99kwDB07dkxdunRRTMzZZ9VEXc9NTEyMunXr1qyfkZiYyP9xWgD3uWVwn1sG97nlcK9bRnPc53P12HgwoRgAAFgK4QYAAFgK4SaE4uLiNG/ePMXFxYW7FEvjPrcM7nPL4D63HO51y2gN9znqJhQDAABro+cGAABYCuEGAABYCuEGAABYCuEGAABYCuHGpLy8PKWlpSk+Pl7p6ekqLCw86/Xr169Xenq64uPj1aNHDz399NMtVGlkM3OfV61apeuvv14XXnihEhMTNXToUL399tstWG3kMvvn2WPTpk1yOBwaNGhQ8xZoEWbvc2VlpebOnavU1FTFxcWpZ8+eWrZsWQtVG7nM3ueXXnpJAwcOVJs2bdS5c2fdfvvtOnz4cAtVG5k2bNigH/zgB+rSpYtsNpveeOONc74mLN+DBgL2yiuvGLGxscazzz5r7Nixw7jnnnuMtm3bGl999VWT1+/evdto06aNcc899xg7duwwnn32WSM2NtZYuXJlC1ceWcze53vuucf47W9/a3z44YfG559/bsyZM8eIjY01PvrooxauPLKYvc8eR48eNXr06GFkZWUZAwcObJliI1gw93nMmDHG1VdfbRQUFBh79uwxtmzZYmzatKkFq448Zu9zYWGhERMTYyxevNjYvXu3UVhYaFx++eXGzTff3MKVR5bVq1cbc+fONV577TVDkvH666+f9fpwfQ8SbkwYPHiwMWXKFL+23r17G7Nnz27y+lmzZhm9e/f2a7vzzjuNIUOGNFuNVmD2Pjelb9++xvz580NdmqUEe5/HjRtn3H///ca8efMINwEwe5//7//+z0hKSjIOHz7cEuVZhtn7/Lvf/c7o0aOHX9uTTz5pdOvWrdlqtJpAwk24vgcZlgpQVVWVioqKlJWV5deelZWlzZs3N/ma999/v9H1N9xwg7Zu3arTp083W62RLJj73FBNTY2OHTum9u3bN0eJlhDsfV6+fLl27dqlefPmNXeJlhDMfX7zzTeVkZGhxx9/XF27dtWll16qmTNn6uTJky1RckQK5j5nZmZq3759Wr16tQzD0Ndff62VK1dq9OjRLVFy1AjX92DUHZwZrLKyMrndbiUnJ/u1Jycn6+DBg02+5uDBg01eX11drbKyMnXu3LnZ6o1Uwdznhn7/+9/r+PHjGjt2bHOUaAnB3OcvvvhCs2fPVmFhoRwO/uoIRDD3effu3dq4caPi4+P1+uuvq6ysTFOnTtWRI0eYd3MGwdznzMxMvfTSSxo3bpxOnTql6upqjRkzRn/4wx9aouSoEa7vQXpuTLLZbH4/G4bRqO1c1zfVDn9m77PHihUr9NBDDyk/P1+dOnVqrvIsI9D77Ha7NX78eM2fP1+XXnppS5VnGWb+PNfU1Mhms+mll17S4MGDddNNNyk3N1d//vOf6b05BzP3eceOHbr77rv14IMPqqioSH//+9+1Z88eTZkypSVKjSrh+B7kn18B6tixo+x2e6N/BRw6dKhRKvW46KKLmrze4XCoQ4cOzVZrJAvmPnvk5+dr0qRJevXVVzVy5MjmLDPimb3Px44d09atW1VcXKxp06ZJqv0SNgxDDodDa9as0bXXXtsitUeSYP48d+7cWV27dlVSUpK3rU+fPjIMQ/v27dMll1zSrDVHomDuc05OjoYNG6bf/OY3kqQBAwaobdu2Gj58uB5++GF61kMkXN+D9NwEyOl0Kj09XQUFBX7tBQUFyszMbPI1Q4cObXT9mjVrlJGRodjY2GarNZIFc5+l2h6b2267TS+//DJj5gEwe58TExO1fft2bdu2zfuYMmWKLrvsMm3btk1XX311S5UeUYL58zxs2DAdOHBAFRUV3rbPP/9cMTEx6tatW7PWG6mCuc8nTpxQTIz/V6DdbpdU37OAby9s34PNOl3ZYjxLDZ977jljx44dxvTp0422bdsaX375pWEYhjF79mxjwoQJ3us9S+BmzJhh7Nixw3juuedYCh4As/f55ZdfNhwOh/HUU08ZpaWl3sfRo0fD9StEBLP3uSFWSwXG7H0+duyY0a1bN+OnP/2p8c9//tNYv369cckllxiTJ08O168QEcze5+XLlxsOh8PIy8szdu3aZWzcuNHIyMgwBg8eHK5fISIcO3bMKC4uNoqLiw1JRm5urlFcXOxdct9avgcJNyY99dRTRmpqquF0Oo0rr7zSWL9+vfe5W2+91bjmmmv8rl+3bp1xxRVXGE6n0+jevbuxZMmSFq44Mpm5z9dcc40hqdHj1ltvbfnCI4zZP8++CDeBM3ufd+7caYwcOdJISEgwunXrZmRnZxsnTpxo4aojj9n7/OSTTxp9+/Y1EhISjM6dOxu/+MUvjH379rVw1ZHlvffeO+vft63le9BmGPS/AQAA62DODQAAsBTCDQAAsBTCDQAAsBTCDQAAsBTCDQAAsBTCDQAAsBTCDQAAsBTCDQAAsBTCDQAAsBTCDQAAsBTCDQAAsBTCDYCI98033+iiiy7So48+6m3bsmWLnE6n1qxZE8bKAIQDB2cCsITVq1fr5ptv1ubNm9W7d29dccUVGj16tBYtWhTu0gC0MMINAMv41a9+pXfeeUdXXXWVPv74Y/3jH/9QfHx8uMsC0MIINwAs4+TJk+rXr5/27t2rrVu3asCAAeEuCUAYMOcGgGXs3r1bBw4cUE1Njb766qtwlwMgTOi5AWAJVVVVGjx4sAYNGqTevXsrNzdX27dvV3JycrhLA9DCCDcALOE3v/mNVq5cqY8//ljnnXeeRowYoXbt2umtt94Kd2kAWhjDUgAi3rp167Ro0SK98MILSkxMVExMjF544QVt3LhRS5YsCXd5AFoYPTcAAMBS6LkBAACWQrgBAACWQrgBAACWQrgBAACWQrgBAACWQrgBAACWQrgBAACWQrgBAACWQrgBAACWQrgBAACWQrgBAACW8v8B5kPeeV2PG8EAAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "data = np.genfromtxt(\"task13/trapped_concentration.txt\", delimiter=\",\", skip_header=1)\n",
    "\n",
    "sorted_data = data[data[:, 0].argsort()]\n",
    "plt.plot(sorted_data[:, 0], sorted_data[:, 1])\n",
    "\n",
    "plt.xlabel(\"x\")\n",
    "plt.ylabel(\"trapped concentration\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "By changing the trap element type to DG, we can greatly mitigate the Gibbs phenomenon. As before, the DG projection in ``TXTExport`` isn't perfect and one should use ``XDMFExport`` for a cleaner plot."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Defining initial values\n",
      "Defining variational problem\n",
      "Defining source terms\n",
      "Defining boundary conditions\n",
      "Solving steady state problem...\n",
      "Solved problem in 0.00 s\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjcAAAG2CAYAAACDLKdOAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAAA2OUlEQVR4nO3de3RU9b3+8WdymQSQxGI0cgkxggoYbiZVEn7UohhFWmptJZUWvMBZTZUiRERy0CocNdraEKwGiwU8WMQoorWetBJRIIBSidHSQtUCmoATY0LNJAoTMtm/PzAjQy7MDpPszOb9WmuWzDd7z3xmL8w8fG/bYRiGIQAAAJsIs7oAAACAYCLcAAAAWyHcAAAAWyHcAAAAWyHcAAAAWyHcAAAAWyHcAAAAWyHcAAAAWyHcAAAAWyHcAAAAW7E83BQUFCgpKUnR0dFKSUlRSUlJm8fefPPNcjgcLR4XX3xxF1YMAAC6M0vDTWFhoebMmaOFCxeqrKxM48aN08SJE1VeXt7q8UuXLpXL5fI9Kioq1KdPH91www1dXDkAAOiuHFbeOPOyyy7TJZdcomXLlvnahg4dquuuu065ubknPf/ll1/W9ddfr/379ysxMbEzSwUAACEiwqo3bmhoUGlpqRYsWODXnpGRoe3btwf0GitWrNCECRPaDTYej0cej8f3vKmpSYcOHdJZZ50lh8PRseIBAECXMgxDdXV16tevn8LC2h94sizcVFdXy+v1Kj4+3q89Pj5elZWVJz3f5XLpL3/5i5599tl2j8vNzdWiRYtOqVYAANA9VFRUaMCAAe0eY1m4aXZi74lhGAH1qDz99NM688wzdd1117V7XE5OjrKzs33Pa2trNXDgQFVUVCgmJqZDNQMAgK7ldruVkJCg3r17n/RYy8JNXFycwsPDW/TSVFVVtejNOZFhGFq5cqWmTZsmp9PZ7rFRUVGKiopq0R4TE0O4AQAgxATSAWLZaimn06mUlBQVFxf7tRcXFys9Pb3dczdv3qx///vfmjFjRmeWCAAAQpClw1LZ2dmaNm2aUlNTlZaWpuXLl6u8vFxZWVmSjg0pHTx4UKtXr/Y7b8WKFbrsssuUnJxsRdkAAKAbszTcZGZmqqamRosXL5bL5VJycrKKiop8q59cLleLPW9qa2v14osvaunSpVaUDAAAujlL97mxgtvtVmxsrGpra5lzAwBAiDDz/W357RcAAACCiXADAABshXADAABshXADAABshXADAABshXADAABshXADAABshXADAABshXADIGQdOerVabYPKYAAEG4AhKSquiNKfeB1ZT//vtWlAOhmCDcAQtK/XHWq9zTq/QNfWF0KgG6GcAMgJNV7Gq0uAUA3RbgBEJLqjxBuALSOcAMgJNXRcwOgDYQbACGJnhsAbSHcAAhJ9Z6jVpcAoJsi3AAISUwoBtAWwg2AkORmWApAGwg3AEISc24AtIVwAyAkMSwFoC2EGwAhiZ4bAG0h3AAISfTcAGgL4QZASKo7wlJwAK0j3AAIOYZh0HMDoE2EGwAh5/BRr5oMq6sA0F0RbgCEHCYTA2gP4QZAyOGmmQDaQ7gBEHLouQHQHsINgJDDZGIA7SHcAAg5dfTcAGgH4QZAyKHnBkB7CDcAQk49G/gBaAfhBkDIoecGQHsINwBCDkvBAbSHcAMg5LAUHEB7CDcAQg7DUgDaQ7gBEHLouQHQHsINgJDDnBsA7SHcAAg59NwAaI/l4aagoEBJSUmKjo5WSkqKSkpK2j3e4/Fo4cKFSkxMVFRUlAYNGqSVK1d2UbUAugPm3ABoT4SVb15YWKg5c+aooKBAY8eO1e9//3tNnDhRu3fv1sCBA1s9Z8qUKfrss8+0YsUKDR48WFVVVWps5BcdcDoh3ABoj6XhJi8vTzNmzNDMmTMlSfn5+Xrttde0bNky5ebmtjj+r3/9qzZv3qx9+/apT58+kqTzzjuvK0sG0A0wLAWgPZYNSzU0NKi0tFQZGRl+7RkZGdq+fXur57zyyitKTU3Vr3/9a/Xv318XXnih5s2bp8OHD3dFyQC6AU+jVw3eJqvLANCNWdZzU11dLa/Xq/j4eL/2+Ph4VVZWtnrOvn37tHXrVkVHR+ull15SdXW1brvtNh06dKjNeTcej0cej8f33O12B+9DAOhy9NoAOBnLJxQ7HA6/54ZhtGhr1tTUJIfDoTVr1ujSSy/Vtddeq7y8PD399NNt9t7k5uYqNjbW90hISAj6ZwDQdZhvA+BkLAs3cXFxCg8Pb9FLU1VV1aI3p1nfvn3Vv39/xcbG+tqGDh0qwzB04MCBVs/JyclRbW2t71FRURG8DwGgy9XRcwPgJCwLN06nUykpKSouLvZrLy4uVnp6eqvnjB07Vp9++qnq6+t9bR9++KHCwsI0YMCAVs+JiopSTEyM3wNA6KLnBsDJWDoslZ2drT/84Q9auXKl9uzZo7lz56q8vFxZWVmSjvW6TJ8+3Xf81KlTddZZZ+mWW27R7t27tWXLFt1111269dZb1aNHD6s+BoAuxJwbACdj6VLwzMxM1dTUaPHixXK5XEpOTlZRUZESExMlSS6XS+Xl5b7jzzjjDBUXF+uXv/ylUlNTddZZZ2nKlCl64IEHrPoIALpYc8+NMyJMDY2smgLQksMwDMPqIrqS2+1WbGysamtrGaICQtAzb3+ie1/+h87q5VTNlw06/+xeeuPO71pdFoBOZub72/LVUgBgRvOw1BnRlnY8A+jGCDcAQkq956gk6Ywowg2A1hFuAIQUX88N4QZAGwg3AEJK3dcTinszLAWgDYQbACGFnhsAJ0O4ARBSmpeCM6EYQFsINwBCii/cREVaXAmA7opwAyCkNA9LMecGQFsINwBCSp2HOTcA2ke4ARBSmFAM4GQINwBCRqO3SYePeiUxoRhA2wg3AELGlx6v78/03ABoC+EGQMio+/rWC1ERYYoM59cXgNbx2wFAyKhnd2IAASDcAAgZTCYGEAjCDYCQUcfuxAACQLgBEDLouQEQCMINgJDBrRcABIJwAyBkcOsFAIEg3AAIGXVHji0FZ1gKQHsINwBCBhOKAQSCcAMgZDChGEAgCDcAQgab+AEIBOEGQMgg3AAIBOEGQMioO8JScAAnR7gBEDK+2eeGnhsAbSPcAAgZ7HMDIBCEGwAhg54bAIEg3AAICU1Nxjfhhp4bAO0g3AAICV82NPr+TM8NgPYQbgCEhOZem8hwh6Ii+NUFoG38hgAQEo7fndjhcFhcDYDujHADICRwXykAgSLcAAgJ9WzgByBAhBsAIcF36wUmEwM4CcINgJDg67lhWArASRBuAISEOjbwAxAgwg2AkEDPDYBAEW4AhIR6z1FJzLkBcHKEGwAhgftKAQiU5eGmoKBASUlJio6OVkpKikpKSto8dtOmTXI4HC0e//rXv7qwYgBWqGNYCkCALA03hYWFmjNnjhYuXKiysjKNGzdOEydOVHl5ebvnffDBB3K5XL7HBRdc0EUVA7AKPTcAAmXpb4m8vDzNmDFDM2fOlCTl5+frtdde07Jly5Sbm9vmeeecc47OPPPMLqoyMN4mQ67aw1aXAdjWB5V1kqTe9NwAOAnLfks0NDSotLRUCxYs8GvPyMjQ9u3b2z139OjROnLkiIYNG6Z77rlH48ePb/NYj8cjj8fje+52u0+t8DbUfOnR/3vkzU55bQDfYIdiACdjWbiprq6W1+tVfHy8X3t8fLwqKytbPadv375avny5UlJS5PF49Mwzz+jKK6/Upk2b9J3vfKfVc3Jzc7Vo0aKg198a7lQMdK6BfXpqZEKsJCn86//dvE2GhRUB6I4s79898e6+hmG0ecffiy66SBdddJHveVpamioqKvToo4+2GW5ycnKUnZ3te+52u5WQkBCEyv2d0ztaHzwwMeivC6B1zvBwSVJDY5PFlQDobizraoiLi1N4eHiLXpqqqqoWvTntGTNmjD766KM2fx4VFaWYmBi/B4DQ5/y6p5RwA+BEloUbp9OplJQUFRcX+7UXFxcrPT094NcpKytT3759g10egG6OcAOgLZYOS2VnZ2vatGlKTU1VWlqali9frvLycmVlZUk6NqR08OBBrV69WtKx1VTnnXeeLr74YjU0NOiPf/yjXnzxRb344otWfgwAFmgONx4v4QaAP0vDTWZmpmpqarR48WK5XC4lJyerqKhIiYmJkiSXy+W3501DQ4PmzZungwcPqkePHrr44ov1f//3f7r22mut+ggALOIM/6bnpr25egBOPw7DME6rpQZut1uxsbGqra1l/g0QwmoPH9XIRRskSR8+MNHXkwPAnsx8f/PbAEBIOn7rhQaGpgAcp0PDUh9++KE2bdqkqqoqNTX5/1L51a9+FZTCAKA9zcNS0teTiqMsLAZAt2I63Dz11FP6xS9+obi4OJ177rl+49wOh4NwA6BLhIU5FBHmUGOTwYopAH5Mh5sHHnhADz74oO6+++7OqAcAAuaMCFNjg5dwA8CP6Tk3//nPf3TDDTd0Ri0AYIpvrxuv1+JKAHQnpsPNDTfcoA0bNnRGLQBgSvO8Gw89NwCOY3pYavDgwbr33nv19ttva/jw4YqM9L9D7+zZs4NWHAC0h12KAbTGdLhZvny5zjjjDG3evFmbN2/2+5nD4SDcAOgyhBsArTEdbvbv398ZdQCAab5ditnnBsBxTmkTP8MwdJptcAygG4mi5wZAKzoUblavXq3hw4erR48e6tGjh0aMGKFnnnkm2LUBQLsYlgLQGtPDUnl5ebr33ns1a9YsjR07VoZhaNu2bcrKylJ1dbXmzp3bGXUCQAvfLAUn3AD4hulw87vf/U7Lli3T9OnTfW0/+MEPdPHFF+v+++8n3ADoMiwFB9Aa08NSLpdL6enpLdrT09PlcrmCUhQABIJhKQCtMR1uBg8erOeff75Fe2FhoS644IKgFAUAgXBGhEsi3ADwZ3pYatGiRcrMzNSWLVs0duxYORwObd26VRs3bmw19ABAZ2EpOIDWmO65+dGPfqQdO3YoLi5OL7/8stavX6+4uDj97W9/0w9/+MPOqBEAWsWwFIDWmO65kaSUlBT98Y9/DHYtAGAK+9wAaE1A4cbtdismJsb35/Y0HwcAnY2l4ABaE1C4+da3viWXy6VzzjlHZ555phwOR4tjDMOQw+GQ1+sNepEA0BrfnBt6bgAcJ6Bw88Ybb6hPnz6SpDfffLNTCwKAQDX33LDPDYDjBRRuLr/8ct+fk5KSlJCQ0KL3xjAMVVRUBLc6AGgHE4oBtMb0aqmkpCR9/vnnLdoPHTqkpKSkoBQFAIH4ZodihsMBfMN0uGmeW3Oi+vp6RUdHB6UoAAgEPTcAWhPwUvDs7GxJksPh0L333quePXv6fub1erVjxw6NGjUq6AUCQFtYLQWgNQGHm7KyMknHem527dolp9Pp+5nT6dTIkSM1b9684FcIAG1gnxsArQk43DSvkrrlllu0dOlS9rMBYDmWggNojekdiletWtUZdQCAaQxLAWhNh26/8M477+iFF15QeXm5Ghoa/H62fv36oBQGACfDhGIArTG9Wuq5557T2LFjtXv3br300ks6evSodu/erTfeeEOxsbGdUSMAtIphKQCtMR1uHnroIS1ZskSvvvqqnE6nli5dqj179mjKlCkaOHBgZ9QIAK1ih2IArTEdbvbu3atJkyZJkqKiovTll1/K4XBo7ty5Wr58edALBIC2MOcGQGtMh5s+ffqorq5OktS/f3/94x//kCR98cUX+uqrr4JbHQC0g6XgAFpjekLxuHHjVFxcrOHDh2vKlCm644479MYbb6i4uFhXXnllZ9QIAK1yhodLItwA8Gc63Dz++OM6cuSIJCknJ0eRkZHaunWrrr/+et17771BLxAA2sKwFIDWmAo3jY2N+vOf/6yrr75akhQWFqb58+dr/vz5nVIcALSnOdx4mwx5mwyFh7W87x2A04+pOTcRERH6xS9+IY/H01n1AEDAmsONxNAUgG+YnlB82WWX+e4zBQBWat7nRiLcAPiG6Tk3t912m+68804dOHBAKSkp6tWrl9/PR4wYEbTiAKA9keHfDEN5vF5JkdYVA6DbMN1zk5mZqf3792v27NkaO3asRo0apdGjR/v+a1ZBQYGSkpIUHR2tlJQUlZSUBHTetm3bFBERoVGjRpl+TwD24HA4uAUDgBZM99zs378/aG9eWFioOXPmqKCgQGPHjtXvf/97TZw4Ubt37253t+Pa2lpNnz5dV155pT777LOg1QMg9ESFh6mhsYlwA8DHdM/NJ598ov79+ysxMdHv0b9/f33yySemXisvL08zZszQzJkzNXToUOXn5yshIUHLli1r97yf//znmjp1qtLS0syWD8BmWA4O4ESmw8348eN16NChFu21tbUaP358wK/T0NCg0tJSZWRk+LVnZGRo+/btbZ63atUq7d27V/fdd19A7+PxeOR2u/0eAOyDYSkAJzIdbgzDkMPRci+JmpqaFpOL21NdXS2v16v4+Hi/9vj4eFVWVrZ6zkcffaQFCxZozZo1iogIbEQtNzdXsbGxvkdCQkLANQLo/gg3AE4U8Jyb66+/XtKxCXw333yzoqKifD/zer36+9//rvT0dNMFnBiU2gpPXq9XU6dO1aJFi3ThhRcG/Po5OTnKzs72PXe73QQcwEaal4MTbgA0CzjcxMbGSjoWPnr37q0ePXr4fuZ0OjVmzBj913/9V8BvHBcXp/Dw8Ba9NFVVVS16cySprq5OO3fuVFlZmWbNmiVJampqkmEYioiI0IYNG3TFFVe0OC8qKsoviAGwl+aeGw9zbgB8LeBws2rVKknSeeedp3nz5pkagmqN0+lUSkqKiouL9cMf/tDXXlxcrB/84Actjo+JidGuXbv82goKCvTGG29o3bp1SkpKOqV6AIQmhqUAnMj0UvBAJ/IGIjs7W9OmTVNqaqrS0tK0fPlylZeXKysrS9KxIaWDBw9q9erVCgsLU3Jyst/555xzjqKjo1u0Azh9MCwF4ESmw81nn32mefPmaePGjaqqqpJhGH4/93q9Ab9WZmamampqtHjxYrlcLiUnJ6uoqEiJiYmSJJfLpfLycrMlAjiN0HMD4EQO48R0chITJ05UeXm5Zs2apb59+7aY/NvakFJ34na7FRsbq9raWsXExFhdDoBTNPN/39Hre6qUe/1w3Xhp25t/AghtZr6/TffcbN26VSUlJdz2AEC3QM8NgBOZ3ucmISGhxVAUAFiFOTcATmQ63OTn52vBggX6+OOPO6EcADCH2y8AOJHpYanMzEx99dVXGjRokHr27KnIyEi/n7d2awYA6Cy+fW7ouQHwNdPhJj8/vxPKAICOcYaHS2JYCsA3TIebm266qTPqAIAOYUIxgBOZnnMjSXv37tU999yjG2+8UVVVVZKkv/71r/rnP/8Z1OIA4GS+mXMT+B5bAOzNdLjZvHmzhg8frh07dmj9+vWqr6+XJP39738P6u7FABCIKHpuAJzAdLhZsGCBHnjgARUXF8vpdPrax48fr7feeiuoxQHAybAUHMCJTIebXbt2+d3ostnZZ5+tmpqaoBQFAIFiKTiAE5kON2eeeaZcLleL9rKyMvXv3z8oRQFAoJhQDOBEpsPN1KlTdffdd6uyslIOh0NNTU3atm2b5s2bp+nTp3dGjQDQpuZhKfa5AdDMdLh58MEHNXDgQPXv31/19fUaNmyYvvOd7yg9PV333HNPZ9QIAG2i5wbAiUzvcxMZGak1a9bof/7nf/Tuu++qqalJo0eP1gUXXNAZ9QFAu5hzA+BEpsNNs/PPP1/nn39+MGsBANPouQFwItPDUj/+8Y/18MMPt2j/zW9+oxtuuCEoRQFAoKJYCg7gBB3axG/SpEkt2q+55hpt2bIlKEUBQKAYlgJwItPhpr6+3m/zvmaRkZFyu91BKQoAAsWwFIATmQ43ycnJKiwsbNH+3HPPadiwYUEpCgACRbgBcCLTE4rvvfde/ehHP9LevXt1xRVXSJI2btyotWvX6oUXXgh6gQDQHm6/AOBEpsPN5MmT9fLLL+uhhx7SunXr1KNHD40YMUKvv/66Lr/88s6oEQDa1Nxz42HODYCvdWgp+KRJk1qdVAwAXe34YSnDMORwOCyuCIDVOrzPTUNDg6qqqtTU5P+vpYEDB55yUQAQqKjwcN+fj3oNOSMIN8DpznS4+eijj3Trrbdq+/btfu3N/2Lyer1BKw4ATqa550Y6thz8+OcATk+mw83NN9+siIgIvfrqq+rbty9dwAAs5RduGpukKAuLAdAtmA437733nkpLSzVkyJDOqAcATAkPcyg8zCFvk8GKKQCSOrDPzbBhw1RdXd0ZtQBAh7AcHMDxTIebRx55RPPnz9emTZtUU1Mjt9vt9wCArvbNLRiY8wegA8NSEyZMkCRdeeWVfu1MKAZgFd9eN/TcAFAHws2bb77ZGXUAQIcxLAXgeKbDDbsQA+huori/FIDjdGgTvy+++EIrVqzQnj175HA4NGzYMN16662KjY0Ndn0AcFLfzLkh3ADowITinTt3atCgQVqyZIkOHTqk6upq5eXladCgQXr33Xc7o0YAaBd3BgdwPNM9N3PnztXkyZP11FNPKSLi2OmNjY2aOXOm5syZoy1btgS9SABoD3NuABzPdLjZuXOnX7CRpIiICM2fP1+pqalBLQ4AAsGwFIDjmR6WiomJUXl5eYv2iooK9e7dOyhFAYAZLAUHcDzT4SYzM1MzZsxQYWGhKioqdODAAT333HOaOXOmbrzxxs6oEQDaxbAUgOOZHpZ69NFH5XA4NH36dDU2NkqSIiMj9Ytf/EIPP/xw0AsEgJOJigyXRLgBcIzpcON0OrV06VLl5uZq7969MgxDgwcPVs+ePTujPgA4KV/PDXNuAKgD4aa2tlZer1d9+vTR8OHDfe2HDh1SRESEYmJiglogAJwMS8EBHM/0nJuf/OQneu6551q0P//88/rJT35iuoCCggIlJSUpOjpaKSkpKikpafPYrVu3auzYsTrrrLPUo0cPDRkyREuWLDH9ngDshR2KARzPdLjZsWOHxo8f36L9u9/9rnbs2GHqtQoLCzVnzhwtXLhQZWVlGjdunCZOnNjqaixJ6tWrl2bNmqUtW7Zoz549uueee3TPPfdo+fLlZj8GABthKTiA45kONx6PxzeR+HhHjx7V4cOHTb1WXl6eZsyYoZkzZ2ro0KHKz89XQkKCli1b1urxo0eP1o033qiLL75Y5513nn72s5/p6quvbre3B4D9sVoKwPFMh5tvf/vbrfaUPPnkk0pJSQn4dRoaGlRaWqqMjAy/9oyMDG3fvj2g1ygrK9P27dvbvZmnx+OR2+32ewCwF/a5AXA80xOKH3zwQU2YMEHvv/++rrzySknSxo0b9c4772jDhg0Bv051dbW8Xq/i4+P92uPj41VZWdnuuQMGDNDnn3+uxsZG3X///Zo5c2abx+bm5mrRokUB1wUg9DChGMDxTPfcjB07Vm+99ZYSEhL0/PPP689//rMGDx6sv//97xo3bpzpAhwOh99zwzBatJ2opKREO3fu1JNPPqn8/HytXbu2zWNzcnJUW1vre1RUVJiuEUD3xlJwAMcz3XMjSaNGjdKaNWtO6Y3j4uIUHh7eopemqqqqRW/OiZKSkiRJw4cP12effab777+/zd2Ro6KiFBUVdUq1Aujevum58VpcCYDuwHTPTbA4nU6lpKSouLjYr724uFjp6ekBv45hGPJ4PMEuD0AIYVgKwPE61HMTLNnZ2Zo2bZpSU1OVlpam5cuXq7y8XFlZWZKODSkdPHhQq1evliQ98cQTGjhwoIYMGSLp2L43jz76qH75y19a9hkAWC+KpeAAjmNpuMnMzFRNTY0WL14sl8ul5ORkFRUVKTExUZLkcrn89rxpampSTk6O9u/fr4iICA0aNEgPP/ywfv7zn1v1EQB0AywFB3A8h2EYhtVFdCW3263Y2FjV1tZyqwjAJjbu+Uwz/nenRg6I1Z9m/T+rywHQCcx8f1s25wYAgoV9bgAcL6Bhqeuvvz7gF1y/fn2HiwGAjmApOIDjBdRzExsb63vExMRo48aN2rlzp+/npaWl2rhxo2JjYzutUABoC6ulABwvoJ6bVatW+f589913a8qUKXryyScVHh4uSfJ6vbrtttuYwwLAEoQbAMczPedm5cqVmjdvni/YSFJ4eLiys7O1cuXKoBYHAIFgKTiA45kON42NjdqzZ0+L9j179qipiV8sALqe8+t/bHmO8jsIQAf2ubnlllt066236t///rfGjBkjSXr77bf18MMP65Zbbgl6gQBwMk56bgAcx3S4efTRR3XuuedqyZIlcrlckqS+fftq/vz5uvPOO4NeIACcTHO48TYZ8jYZCg9r/+a7AOzNdLgJCwvT/PnzNX/+fLndbkliIjEASzWHG+nYpOIezvB2jgZgdx3axK+xsVGvv/661q5dK4fj2L+QPv30U9XX1we1OAAIRPM+NxIrpgB0oOfmk08+0TXXXKPy8nJ5PB5dddVV6t27t37961/ryJEjevLJJzujTgBoU2T4N8NQHq9XUqR1xQCwnOmemzvuuEOpqan6z3/+ox49evjaf/jDH2rjxo1BLQ4AAuFwONjrBoCP6Z6brVu3atu2bXI6nX7tiYmJOnjwYNAKAwAzosLD1NDYRLgBYL7npqmpSV6vt0X7gQMH1Lt376AUBQBmsRwcQDPT4eaqq65Sfn6+77nD4VB9fb3uu+8+XXvttcGsDQACxrAUgGamh6WWLFmi8ePHa9iwYTpy5IimTp2qjz76SHFxcVq7dm1n1AgAJ0W4AdDMdLjp16+f3nvvPa1du1bvvvuumpqaNGPGDP30pz/1m2AMAF2peTk44QaA6XAjST169NCtt96qW2+9Ndj1AECHNPfceJhzA5z2OhRuPvjgA/3ud7/Tnj175HA4NGTIEM2aNUtDhgwJdn0AEBCGpQA0Mz2heN26dUpOTlZpaalGjhypESNG6N1339Xw4cP1wgsvdEaNAHBSDEsBaGa652b+/PnKycnR4sWL/drvu+8+3X333brhhhuCVhwABIqeGwDNTPfcVFZWavr06S3af/azn6mysjIoRQGAWVHscwPga6bDzXe/+12VlJS0aN+6davGjRsXlKIAwCx6bgA0Mz0sNXnyZN19990qLS3VmDFjJElvv/22XnjhBS1atEivvPKK37EA0BWYcwOgmcMwDMPMCWFhgXX2OByOVm/TYDW3263Y2FjV1tYqJibG6nIABMn8de/r+Z0HdNfVF+n28YOtLgdAkJn5/jbdc9PUxL+KAHQ/vn1u6LkBTnum59wAQHfkDA+XxLAUgA6Gm40bN+p73/ueBg0apMGDB+t73/ueXn/99WDXBgABY0IxgGamw83jjz+ua665Rr1799Ydd9yh2bNnKyYmRtdee60ef/zxzqgRAE7KF2664Vw/AF3L9Jyb3NxcLVmyRLNmzfK1zZ49W2PHjtWDDz7o1w4AXSWKnhsAXzPdc+N2u3XNNde0aM/IyJDb7Q5KUQBgFkvBATQzHW4mT56sl156qUX7n/70J33/+98PSlEAYJaTHYoBfM30sNTQoUP14IMPatOmTUpLS5N0bBO/bdu26c4779Rjjz3mO3b27NnBqxQA2sGEYgDNTIebFStW6Fvf+pZ2796t3bt3+9rPPPNMrVixwvfc4XAQbgB0meZhKfa5AWA63Ozfv78z6gCAU0LPDYBmbOIHwBaYcwOgmemeG0k6cOCAXnnlFZWXl6uhocHvZ3l5eUEpDADMoOcGQDPT4Wbjxo2aPHmykpKS9MEHHyg5OVkff/yxDMPQJZdc0hk1AsBJRbEUHMDXTA9L5eTk6M4779Q//vEPRUdH68UXX1RFRYUuv/xy3XDDDZ1RIwCcFMNSAJqZDjd79uzRTTfdJEmKiIjQ4cOHdcYZZ2jx4sV65JFHTBdQUFCgpKQkRUdHKyUlRSUlJW0eu379el111VU6++yzFRMTo7S0NL322mum3xOA/TAsBaCZ6XDTq1cveTweSVK/fv20d+9e38+qq6tNvVZhYaHmzJmjhQsXqqysTOPGjdPEiRNVXl7e6vFbtmzRVVddpaKiIpWWlmr8+PH6/ve/r7KyMrMfA4DNEG4ANDM952bMmDHatm2bhg0bpkmTJunOO+/Url27tH79eo0ZM8bUa+Xl5WnGjBmaOXOmJCk/P1+vvfaali1bptzc3BbH5+fn+z1/6KGH9Kc//Ul//vOfNXr0aLMfBYCNcPsFAM1Mh5u8vDzV19dLku6//37V19ersLBQgwcP1pIlSwJ+nYaGBpWWlmrBggV+7RkZGdq+fXtAr9HU1KS6ujr16dOnzWM8Ho+vp0kS978CbKq558bDnBvgtGcq3Hi9XlVUVGjEiBGSpJ49e6qgoKBDb1xdXS2v16v4+Hi/9vj4eFVWVgb0Gr/97W/15ZdfasqUKW0ek5ubq0WLFnWoRgCh4/hhKcMw5HA4LK4IgFVMzbkJDw/X1VdfrS+++CJoBZz4CyjQX0pr167V/fffr8LCQp1zzjltHpeTk6Pa2lrfo6Ki4pRrBtD9RIWH+/581GtYWAkAq5kelho+fLj27dunpKSkU3rjuLg4hYeHt+ilqaqqatGbc6LCwkLNmDFDL7zwgiZMmNDusVFRUYqKijqlWgF0f809N9Kx5eDHPwdwejH9f/+DDz6oefPm6dVXX5XL5ZLb7fZ7BMrpdColJUXFxcV+7cXFxUpPT2/zvLVr1+rmm2/Ws88+q0mTJpktH4BN+YUbJhUDpzXTPTfXXHONJGny5Ml+w0fNw0lerzfg18rOzta0adOUmpqqtLQ0LV++XOXl5crKypJ0bEjp4MGDWr16taRjwWb69OlaunSpxowZ4+v16dGjh2JjY81+FAA2Eh7mUHiYQ94mg3ADnOZMh5s333wzaG+emZmpmpoaLV68WC6XS8nJySoqKlJiYqIkyeVy+e158/vf/16NjY26/fbbdfvtt/vab7rpJj399NNBqwtAaHKGh+lwk5dwA5zmHIZhmJp5V15eroSEhFYnAldUVGjgwIFBLTDY3G63YmNjVVtbq5iYGKvLARBEIxdtUO3ho3o9+zsafE5vq8sBEERmvr9Nz7lJSkrS559/3qL90KFDpzzJGABOhW+vG3pugNOa6XDT1lLt+vp6RUdHB6UoAOgIdikGIJmYc5OdnS3p2L409957r3r27On7mdfr1Y4dOzRq1KigFwgAgYri/lIAZCLcNN+c0jAM7dq1S06n0/czp9OpkSNHat68ecGvEAAC5NulmFswAKe1gMNN8yqpW265RUuXLmUyLoBuhzuDA5A6sBR81apVnVEHAJwy5twAkDowoRgAuiuGpQBIhBsANsJScAAS4QaAjTAsBUAi3ACwESYUA5AINwBsJCLs2AajTebuKgPAZgg3AADAVgg3AADAVgg3AADAVgg3AADAVgg3AADAVgg3AADAVgg3AADAVgg3AADAVgg3AADAVgg3AADAVgg3AADAVgg3AADAVgg3AADAVgg3AADAVgg3AADAVgg3AADAVgg3AADAVgg3AADAVgg3AADAVgg3AADAVgg3AADAVgg3AADAVgg3AADAVgg3AADAVgg3AADAVgg3AADAVgg3AADAVgg3AADAViwPNwUFBUpKSlJ0dLRSUlJUUlLS5rEul0tTp07VRRddpLCwMM2ZM6frCgUAACHB0nBTWFioOXPmaOHChSorK9O4ceM0ceJElZeXt3q8x+PR2WefrYULF2rkyJFdXC0AAAgFloabvLw8zZgxQzNnztTQoUOVn5+vhIQELVu2rNXjzzvvPC1dulTTp09XbGxsF1cLAABCgWXhpqGhQaWlpcrIyPBrz8jI0Pbt2y2qCgAAhLoIq964urpaXq9X8fHxfu3x8fGqrKwM2vt4PB55PB7fc7fbHbTXBgAA3Y/lE4odDoffc8MwWrSditzcXMXGxvoeCQkJQXttAADQ/VgWbuLi4hQeHt6il6aqqqpFb86pyMnJUW1tre9RUVERtNcGAADdj2Xhxul0KiUlRcXFxX7txcXFSk9PD9r7REVFKSYmxu8BAADsy7I5N5KUnZ2tadOmKTU1VWlpaVq+fLnKy8uVlZUl6Vivy8GDB7V69WrfOe+9954kqb6+Xp9//rnee+89OZ1ODRs2zIqPAAAAuhlLw01mZqZqamq0ePFiuVwuJScnq6ioSImJiZKObdp34p43o0eP9v25tLRUzz77rBITE/Xxxx93ZekAAKCbsjTcSNJtt92m2267rdWfPf300y3aDMPo5IoAAEAos3y1FAAAQDARbgAAgK0QbgAAgK0QbgAAgK0QbgAAgK0QbgAAgK0QbgAAgK0QbgAAgK0QbgAAgK0QbgAAgK0QbgAAgK0QbgAAgK0QbgAAgK0QbgAAgK0QbgAAgK0QbgAAgK0QbgAAgK0QbgAAgK0QbgAAgK0QbgAAgK0QbgAAgK0QbgAAgK0QbgAAgK0QbgAAgK0QbgAAgK0QbgAAgK0QbgAAgK0QbgAAgK0QbgAAgK0QbgAAgK0QbgAAgK0QbgAAgK0QbgAAgK0QbgAAgK0QbgAAgK0QbgAAgK0QbgAAgK0QbgAAgK0QbgAAgK0QbgAAgK1YHm4KCgqUlJSk6OhopaSkqKSkpN3jN2/erJSUFEVHR+v888/Xk08+2UWVAgCAUGBpuCksLNScOXO0cOFClZWVady4cZo4caLKy8tbPX7//v269tprNW7cOJWVlem///u/NXv2bL344otdXDkAAOiuHIZhGFa9+WWXXaZLLrlEy5Yt87UNHTpU1113nXJzc1scf/fdd+uVV17Rnj17fG1ZWVl6//339dZbbwX0nm63W7GxsaqtrVVMTMypfwgA3cac58r08nufatb4wfrJpQlWlwOctsLDHOob2yOor2nm+zsiqO9sQkNDg0pLS7VgwQK/9oyMDG3fvr3Vc9566y1lZGT4tV199dVasWKFjh49qsjIyBbneDweeTwe33O32x2E6gF0Z4+/+W89/ua/rS4DOG2d0ztKf1s4wbL3tyzcVFdXy+v1Kj4+3q89Pj5elZWVrZ5TWVnZ6vGNjY2qrq5W3759W5yTm5urRYsWBa9wAN3WhGHxeuNfVfI0NlldCnBai4q0dkqvZeGmmcPh8HtuGEaLtpMd31p7s5ycHGVnZ/ueu91uJSTQXQ3Y0fdG9NP3RvSzugwAFrMs3MTFxSk8PLxFL01VVVWL3plm5557bqvHR0RE6Kyzzmr1nKioKEVFRQWnaAAA0O1Z1m/kdDqVkpKi4uJiv/bi4mKlp6e3ek5aWlqL4zds2KDU1NRW59sAAIDTj6WDYtnZ2frDH/6glStXas+ePZo7d67Ky8uVlZUl6diQ0vTp033HZ2Vl6ZNPPlF2drb27NmjlStXasWKFZo3b55VHwEAAHQzls65yczMVE1NjRYvXiyXy6Xk5GQVFRUpMTFRkuRyufz2vElKSlJRUZHmzp2rJ554Qv369dNjjz2mH/3oR1Z9BAAA0M1Yus+NFdjnBgCA0GPm+9vy2y8AAAAEE+EGAADYCuEGAADYCuEGAADYCuEGAADYCuEGAADYCuEGAADYCuEGAADYCuEGAADYiqW3X7BC84bMbrfb4koAAECgmr+3A7mxwmkXburq6iRJCQkJFlcCAADMqqurU2xsbLvHnHb3lmpqatKnn36q3r17y+FwBPW13W63EhISVFFRwX2rOhHXuWtwnbsG17nrcK27RmddZ8MwVFdXp379+iksrP1ZNaddz01YWJgGDBjQqe8RExPD/zhdgOvcNbjOXYPr3HW41l2jM67zyXpsmjGhGAAA2ArhBgAA2ArhJoiioqJ03333KSoqyupSbI3r3DW4zl2D69x1uNZdoztc59NuQjEAALA3em4AAICtEG4AAICtEG4AAICtEG4AAICtEG5MKigoUFJSkqKjo5WSkqKSkpJ2j9+8ebNSUlIUHR2t888/X08++WQXVRrazFzn9evX66qrrtLZZ5+tmJgYpaWl6bXXXuvCakOX2b/PzbZt26aIiAiNGjWqcwu0CbPX2ePxaOHChUpMTFRUVJQGDRqklStXdlG1ocvsdV6zZo1Gjhypnj17qm/fvrrllltUU1PTRdWGpi1btuj73/+++vXrJ4fDoZdffvmk51jyPWggYM8995wRGRlpPPXUU8bu3buNO+64w+jVq5fxySeftHr8vn37jJ49exp33HGHsXv3buOpp54yIiMjjXXr1nVx5aHF7HW+4447jEceecT429/+Znz44YdGTk6OERkZabz77rtdXHloMXudm33xxRfG+eefb2RkZBgjR47smmJDWEeu8+TJk43LLrvMKC4uNvbv32/s2LHD2LZtWxdWHXrMXueSkhIjLCzMWLp0qbFv3z6jpKTEuPjii43rrruuiysPLUVFRcbChQuNF1980ZBkvPTSS+0eb9X3IOHGhEsvvdTIysryaxsyZIixYMGCVo+fP3++MWTIEL+2n//858aYMWM6rUY7MHudWzNs2DBj0aJFwS7NVjp6nTMzM4177rnHuO+++wg3ATB7nf/yl78YsbGxRk1NTVeUZxtmr/NvfvMb4/zzz/dre+yxx4wBAwZ0Wo12E0i4sep7kGGpADU0NKi0tFQZGRl+7RkZGdq+fXur57z11lstjr/66qu1c+dOHT16tNNqDWUduc4nampqUl1dnfr06dMZJdpCR6/zqlWrtHfvXt13332dXaItdOQ6v/LKK0pNTdWvf/1r9e/fXxdeeKHmzZunw4cPd0XJIakj1zk9PV0HDhxQUVGRDMPQZ599pnXr1mnSpEldUfJpw6rvwdPuxpkdVV1dLa/Xq/j4eL/2+Ph4VVZWtnpOZWVlq8c3Njaqurpaffv27bR6Q1VHrvOJfvvb3+rLL7/UlClTOqNEW+jIdf7oo4+0YMEClZSUKCKCXx2B6Mh13rdvn7Zu3aro6Gi99NJLqq6u1m233aZDhw4x76YNHbnO6enpWrNmjTIzM3XkyBE1NjZq8uTJ+t3vftcVJZ82rPoepOfGJIfD4ffcMIwWbSc7vrV2+DN7nZutXbtW999/vwoLC3XOOed0Vnm2Eeh19nq9mjp1qhYtWqQLL7ywq8qzDTN/n5uamuRwOLRmzRpdeumluvbaa5WXl6enn36a3puTMHOdd+/erdmzZ+tXv/qVSktL9de//lX79+9XVlZWV5R6WrHie5B/fgUoLi5O4eHhLf4VUFVV1SKVNjv33HNbPT4iIkJnnXVWp9UayjpynZsVFhZqxowZeuGFFzRhwoTOLDPkmb3OdXV12rlzp8rKyjRr1ixJx76EDcNQRESENmzYoCuuuKJLag8lHfn73LdvX/Xv31+xsbG+tqFDh8owDB04cEAXXHBBp9YcijpynXNzczV27FjdddddkqQRI0aoV69eGjdunB544AF61oPEqu9Bem4C5HQ6lZKSouLiYr/24uJipaent3pOWlpai+M3bNig1NRURUZGdlqtoawj11k61mNz880369lnn2XMPABmr3NMTIx27dql9957z/fIysrSRRddpPfee0+XXXZZV5UeUjry93ns2LH69NNPVV9f72v78MMPFRYWpgEDBnRqvaGqI9f5q6++UliY/1dgeHi4pG96FnDqLPse7NTpyjbTvNRwxYoVxu7du405c+YYvXr1Mj7++GPDMAxjwYIFxrRp03zHNy+Bmzt3rrF7925jxYoVLAUPgNnr/OyzzxoRERHGE088YbhcLt/jiy++sOojhASz1/lErJYKjNnrXFdXZwwYMMD48Y9/bPzzn/80Nm/ebFxwwQXGzJkzrfoIIcHsdV61apURERFhFBQUGHv37jW2bt1qpKamGpdeeqlVHyEk1NXVGWVlZUZZWZkhycjLyzPKysp8S+67y/cg4cakJ554wkhMTDScTqdxySWXGJs3b/b97KabbjIuv/xyv+M3bdpkjB492nA6ncZ5551nLFu2rIsrDk1mrvPll19uSGrxuOmmm7q+8BBj9u/z8Qg3gTN7nffs2WNMmDDB6NGjhzFgwAAjOzvb+Oqrr7q46tBj9jo/9thjxrBhw4wePXoYffv2NX76058aBw4c6OKqQ8ubb77Z7u/b7vI96DAM+t8AAIB9MOcGAADYCuEGAADYCuEGAADYCuEGAADYCuEGAADYCuEGAADYCuEGAADYCuEGAADYCuEGAADYCuEGAADYCuEGQMj7/PPPde655+qhhx7yte3YsUNOp1MbNmywsDIAVuDeUgBsoaioSNddd522b9+uIUOGaPTo0Zo0aZLy8/OtLg1AFyPcALCN22+/Xa+//rq+/e1v6/3339c777yj6Ohoq8sC0MUINwBs4/Dhw0pOTlZFRYV27typESNGWF0SAAsw5waAbezbt0+ffvqpmpqa9Mknn1hdDgCL0HMDwBYaGhp06aWXatSoURoyZIjy8vK0a9cuxcfHW10agC5GuAFgC3fddZfWrVun999/X2eccYbGjx+v3r1769VXX7W6NABdjGEpACFv06ZNys/P1zPPPKOYmBiFhYXpmWee0datW7Vs2TKrywPQxei5AQAAtkLPDQAAsBXCDQAAsBXCDQAAsBXCDQAAsBXCDQAAsBXCDQAAsBXCDQAAsBXCDQAAsBXCDQAAsBXCDQAAsBXCDQAAsBXCDQAAsJX/D4PeY4WRfx9zAAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "my_model.settings.traps_element_type = \"DG\"\n",
    "my_model.initialise()\n",
    "my_model.run()\n",
    "\n",
    "data = np.genfromtxt(\"task13/trapped_concentration.txt\", delimiter=\",\", skip_header=1)\n",
    "\n",
    "sorted_data = data[data[:, 0].argsort()]\n",
    "\n",
    "plt.plot(sorted_data[:, 0], sorted_data[:, 1])\n",
    "\n",
    "plt.xlabel(\"x\")\n",
    "plt.ylabel(\"trapped concentration\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Conclusion\n",
    "\n",
    "Now you may wonder, why not use DG elements all the time then?\n",
    "\n",
    "Well the thing with DG is that it duplicates degrees of freedom on the edges. Here's an example on three triangular cells with CG elements vs DG elements (order 4).\n",
    "\n",
    "![Alt text](https://canada1.discourse-cdn.com/free1/uploads/festim/original/1X/2e63978aac2a2f81a3e73e9612b4ad7904b88aaa.png)\n",
    "\n",
    "This means that problems using DG elements have overall more degrees of freedom than when using CG.\n",
    "Therefore, since most problems are mono-material, we decided to set the defaults to CG."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "festim-workshop",
   "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.11.9"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}