{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "d12b1a16",
   "metadata": {},
   "outputs": [],
   "source": [
    "method_type = \"HDG\" # available choices \"HDG\", \"DG\" and \"HDG_Newton\"\n",
    "\n",
    "order = 3 # from 0 to ....\n",
    "ictype = 1\n",
    "\n",
    "smax = 1\n",
    "nel = 100\n",
    "\n",
    "hmax = 1./nel\n",
    "dt = hmax \n",
    "\n",
    "\n",
    "from ngsolve import *\n",
    "from netgen.occ import *\n",
    "from netgen.geom2d import SplineGeometry\n",
    "from netgen.meshing import IdentificationType\n",
    "import numpy as np\n",
    "from ngsolve.solvers import Newton\n",
    "\n",
    "\n",
    "from ngsolve.meshes import Make1DMesh\n",
    "mesh = Make1DMesh(nel, periodic = True)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "de4e7132",
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "fesrho = L2(mesh, order = order, dgjumps = True)\n",
    "feshat = Periodic( FacetFESpace(mesh, order  = order))\n",
    "\n",
    "gfrho = GridFunction(fesrho)\n",
    "gfrho_old = GridFunction(fesrho)\n",
    "\n",
    "\n",
    "\n",
    "dS = dx(element_boundary = True)\n",
    "n = specialcf.normal(mesh.dim)\n",
    "h = specialcf.mesh_size\n",
    "\n",
    "\n",
    "myfes = fesrho * feshat\n",
    "(rho, rho_hat), (rho_test, rho_hat_test) = myfes.TnT()\n",
    "\n",
    "wind = CF(1.0)\n",
    "\n",
    "def evolveIE_HDG_Newton(gfrho_old_vec):\n",
    "\n",
    "    gfrhoold = GridFunction(fesrho)\n",
    "    gfrhoold.vec.data = gfrho_old_vec\n",
    "\n",
    "\n",
    "    bf = BilinearForm(myfes, condense = True) # create bilinear form\n",
    "\n",
    "\n",
    "    bf += 1/dt * rho*rho_test * dx - 1/dt* gfrhoold*rho_test*dx \n",
    "    bf += - wind * rho * grad(rho_test)*dx \n",
    "    bf += wind * rho_hat * n* rho_test * dS \n",
    "\n",
    "    bf += (rho_hat - rho)*rho_hat_test*dS \n",
    "\n",
    "    bf.Assemble() # assemble bilinear form\n",
    "    gfA = GridFunction(myfes)\n",
    "    Newton(bf, gfA,maxit=20, dampfactor=1, maxerr = 1e-12)\n",
    "    \n",
    "    return gfA.components[0].vec\n",
    "\n",
    "def evolveIE_HDG(gfrho_old_vec):\n",
    "\n",
    "    # standard implicit Euler HDG \n",
    "    gfrhoold = GridFunction(fesrho)\n",
    "    gfrhoold.vec.data = gfrho_old_vec\n",
    "\n",
    "    bf = BilinearForm(myfes) # create bilinear form\n",
    "\n",
    "\n",
    "    bf += 1/dt * rho*rho_test * dx \n",
    "    bf += - wind * rho * grad(rho_test)*dx \n",
    "    bf += wind* rho_hat *n* rho_test * dS \n",
    "\n",
    "    bf += (rho_hat - rho)*rho_hat_test*dS \n",
    "\n",
    "    bf.Assemble() # assemble bilinear form\n",
    "\n",
    "    rhs = LinearForm(myfes)\n",
    "    rhs += 1/dt * gfrhoold*rho_test*dx \n",
    "\n",
    "    rhs.Assemble() \n",
    "\n",
    "    gfA = GridFunction(myfes)\n",
    "    gfA.vec.data = bf.mat.Inverse(freedofs = myfes.FreeDofs())*rhs.vec\n",
    "\n",
    "    return gfA.components[0].vec\n",
    "\n",
    "def evolveIE_DG(gfrho_old_vec):\n",
    "\n",
    "    # standard implicit Euler DG \n",
    "    gfrhoold = GridFunction(fesrho)\n",
    "    gfrhoold.vec.data = gfrho_old_vec\n",
    "\n",
    "    (r, dr) = fesrho.TnT()\n",
    "    bf = BilinearForm(fesrho) # create bilinear form\n",
    "\n",
    "\n",
    "    bf += 1/dt * r*dr * dx \n",
    "    bf += - wind * r * grad(dr)*dx \n",
    "    bf += wind * 0.5*(r + r.Other()) * n* dr * dS \n",
    "\n",
    "\n",
    "    bf.Assemble() # assemble bilinear form\n",
    "\n",
    "    rhs = LinearForm(fesrho)\n",
    "    rhs += 1/dt * gfrhoold*dr*dx \n",
    "\n",
    "    rhs.Assemble() \n",
    "\n",
    "    gfA = GridFunction(fesrho)\n",
    "    gfA.vec.data = bf.mat.Inverse(freedofs = fesrho.FreeDofs())*rhs.vec\n",
    "\n",
    "    return gfA.vec\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "7be9f0dd",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjgAAAGdCAYAAAAfTAk2AAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjEsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvc2/+5QAAAAlwSFlzAAAPYQAAD2EBqD+naQAAVf1JREFUeJzt3QlYlFXbB/A/MOyyimyyiRuaC4pKuLSor2ulrWqLaS6V2ltZufSWli2WlV+vZWbupqVZaqa+pLlkKorivqGIyA6Csss6813nDENiqKgMzyz/33U98czwzHAzIXNzzrnvY6HRaDQgIiIiMiGWSgdAREREVNeY4BAREZHJYYJDREREJocJDhEREZkcJjhERERkcpjgEBERkclhgkNEREQmhwkOERERmRwVzJBarUZqaiqcnJxgYWGhdDhERERUC6I3cX5+Pnx9fWFpefMxGrNMcERy4+/vr3QYREREdAeSkpLg5+d302vMMsERIze6F8jZ2VnpcIiIiKgW8vLy5ACF7n38ZswywdFNS4nkhgkOERGRcanN8hIuMiYiIiKTwwSHiIiITA4THCIiIjI5THCIiIjI5DDBISIiIpPDBIeIiIhMDhMcIiIiMjlMcIiIiMjkMMEhIiIik6PXBGfXrl14+OGH5aZYouvg+vXrb/mYnTt3omPHjrC1tUWzZs2wdOnSf1wzd+5cBAUFwc7ODuHh4YiOjtbTd0BERETGSK8JTmFhIdq3by8Tktq4cOECBg4ciAcffBBHjhzBa6+9htGjR+P333+vumb16tWYOHEipk+fjkOHDsnn79u3LzIzM/X4nRAREZExsdCIvcfr4wtZWGDdunUYPHjwDa+ZPHkyNm3ahBMnTlTdN3ToUOTk5CAyMlLeFiM2nTt3xtdffy1vq9VqufHWK6+8gilTptR6sy4XFxfk5uZyLyoiIiIjcTvv3wa12WZUVBR69+5d7T4xOiNGcoTS0lLExMRg6tSpVZ+3tLSUjxGPvZGSkhJ5XPsC6cOx5Bz8GJ2EJh4OCGroiCYejgho6ABblZVevh4REZEhKK9QIzWnGPFZBUjIKkRCdhG6NHHHgLY+isVkUAlOeno6vLy8qt0nbouE5OrVq7hy5QoqKipqvObMmTM3fN6ZM2fi/fffh74dSRIJTmK1+ywtgCAPR9zXvBEeDPFEeBN32Fkz4SEiIuNVodbgaHIOdpzJxM7YSziTnoeyiuoTQldLK5jg6JsY8RHrdnREwiSmtepaez9XvNKzGS7I7LUQCVlFKCgpR/ylQnks3ZsAe2srdGvWEA+390X/Nj6wUbGQjYiIDJ9Go8GeuGz8HJOEP89ewpWismqfF+9nQQ3/nsEID3aHkgwqwfH29kZGRka1+8RtMc9mb28PKysredR0jXjsjYiKLHHoW3t/V3lc+8OQVVCKQ4lXsDM2EzvOXEJ6XjH+OJ0pj4+dT+O5ewMxrEsAGjbQf3xERES3S4zErDucgqV7L+BsRkHV/U52KtzXohF6tvSU01GNXe1hKaYtDIRBJTgRERHYvHlztfu2bt0q7xdsbGwQFhaGbdu2VS1WFouMxe0JEybA0IiF1Y2cbNH3Hm95iITndFo+Ik+my6msjLwSfL7lLOZsj8OjoY3x797N5Q8IERGR0vKLyzD/z3is2H8ROZWjNQ42VngizA8PtfNFxwBXqKwMdxZCrwlOQUEB4uLiqpWBi/Jvd3d3BAQEyKmjlJQULF++XH7+pZdektVRkyZNwgsvvIDt27fjp59+kpVVOmKq6fnnn0enTp3QpUsXfPnll7IcfeTIkTB0IuFp7essjwkPNsOm46lYsicBx5JzsfpgEn49moKX7m+KF+9rCnsbrtMhIiJl1tf8HJOEz36PlbMQgr+7PZ6PCMJTnf3hbGcNmHuZuGjaJ3raXE8kKKKB34gRI5CQkCCvu/Yxr7/+Ok6dOgU/Pz+8++678rpriSTos88+k4uSQ0NDMWfOHFk+XluGVCYuXv6DF6/gs8hYRCdclvf5uthhyoBWeLidj0yKiIiI6sOBhMt4/7eTOJGirTYO9nDEW31bos893rAygOmn23n/rrc+OIbEkBIcHfG/YdPxNMzcfAYpOVflfRHBDfHFU+3hy2krIiLSo8KScpnY/HQwuWp9zau9mmN4RJBBFcMwwTHCBEenuKwC3+2Kx7yd53G1rAIu9tb45LG26K9gqR0REZmuo0k5eHXVYdm7RkwaDO0cgDf6tICHARa/MMEx4gRHR5Saix84sT5HGNrZH9Mebg0HG4NaF05EREa81mb+rvOYveUsytUauTzi/4aEIjy4IQwVExwTSHCEsgo1/m/rWcz78zzE/yUxFzrv2TC09HZSOjQiIjJi2QUlmPDDYUTFZ8vbA9v54OPBbeHiYG0y79+GM7FG/2BtZYlJ/ULww+h74e1sh/isQjwxb69ssERERHQn4jLzMfibPTK5EWXfs55oh6+HdTD45OZ2McExAhFNG+J/r/aQ2zzkl5TjhaUHsGLfRaXDIiIiI7MnLguPfrMXSZevIrChAzZM6I6nOvmbZMUuExwj4eZog+9HheOxjo3lvOk760/gg42n5DkREdGt/BidiOcXRyO/uBydAt2wblw3NPNsAFPFBMeIiFK9L55sjzf7tJC3F+2+gJdWxKCkvELp0IiIyEBpNBp89vsZTF17XC4mHhTqixWjw+HuaANTxgTHyIhhxAk9m2POsA4y4dl6KgMvfh8jy8uJiIiuT24+3nwac3ecl7dFb5svh4TCztr0u+UzwTFSj7T3xdIRnWFnbSm3qh+z/CCTHCIiqpbczNh4Cgv+uiBvfzDoHrz+rxYmud6mJkxwjFjXZh5YMqIL7K2t8Ne5LIxedlDu+kpEROZNo9Hg/d9Oyf0OhY8ebYPnIoJgTpjgGDlRYbV0ZGdZ6rc7Lgujlh1AUWm50mEREZFC1GoNpv16Ekv3JsjOxKIb/jPhgTA3THBMgOg6ufyFLnC0scLe89l4ecUh2SSQiIjMz6eRZ/D9vosyufn08XYY2iUA5ogJjonoFOSO5aO6yDU5ohHg5F+OySFKIiIyH6K6dv6ueHn+6WPtZI8bc8UEx4SEBbpj7tMd5Zb2aw+l4LPfY5UOiYiI6slvR1NlfzRhUr+WeKqz+SY3AhMcE9OrlRdmPtpWnn+z8zyW7dUuMCMiItO193wW3vjpqDx/PiIQL9/fFOaOCY4JEln7xH9pmwG+99tJ/O94mtIhERGRnpxKzcOLy2NQWqHGgLbemPbwPWZTCn4zTHBM1Cs9m+Hp8AC5C/mrq4/gaFKO0iEREVEdy8wvlvsTin0KuzRxx+ynQuUyBWKCY7JE9v7BoDboFeKJ0nK17HZ8Kb9E6bCIiKiOiN/t41YcQnpeMZo2csSC5zqZRYfi2mKCY8JEFv/l0FD5gy/+AYxbGSP/QRARkfF7/7eTOHjxCpzsVFgwvBNcHKyVDsmgMMExcU521vIH38lWhQMJVzBj40mlQyIiorv0w/5ErNyfKHvdzBnaAcGNTHdX8DvFBMcMiB/8/w4Llf8QVuxLxI/RiUqHREREd+hgwmVM33BCnr/ZpyUeDPFUOiSDxATHTPQM8ZL/EIRpv55AzMXLSodERES3KT23GC/JbvUaWTE17gGWg98IExwzIv4hiH8Q4h/GhB8O40phqdIhERFRLVWoNfj3qsPIKihBiLcTPnuiPcvBb4IJjhkR/xDEP4hgD0ek5RZjErdzICIyGl9tP4foC5flvoPzng2Do61K6ZAMGhMcMyP+QcwZ1gE2VpbYeioDy6MuKh0SERHdwr74bMzZdk6ef/xYWzTxcFQ6JIPHBMcMtWnsgqkDQuT5R5tO42RqrtIhERHRDVwuLMVrq45ArQGeDPPDoNDGSodkFJjgmKkRXYPQu5WnbO39yg+HUVhSrnRIRER0HbGM4K01R2Uvs+BGjnh/0D1Kh2Q0mOCY+Xocb2c7xGcVYvoG9schIjI0S/YkYNuZTNioLPH1sI5wsOG6m9pigmPG3Bxt8N+hoRDblvwck4xNx7gpJxGRoYhNz8cn/zsjz98d2AqtfZ2VDsmoMMExc+HBDTH+wWby/J31x7lfFRGRASirUGPiT0fkMgKxp+Cz9wYqHZLRYYJDeKVnc7T2ccaVojK8ve44S8eJiBQ2d0ccTqbmwdXBGjMfb8t+N4aa4MydOxdBQUGws7NDeHg4oqOjb3jtAw88IP9HXn8MHDiw6poRI0b84/P9+vWrj2/FJIm53S+eag9rKwtZOr7ucIrSIRERma0TKbn4enucPJ8xqA08neyUDsko6T3BWb16NSZOnIjp06fj0KFDaN++Pfr27YvMzMwar1+7di3S0tKqjhMnTsDKygpPPvlktetEQnPtdT/++KO+vxWT1srHGa/1biHPxYLjtNyrSodERGR2Ssor5NRUuVq7FcPD7XyUDslo6T3BmT17NsaMGYORI0eidevW+Pbbb+Hg4IDFixfXeL27uzu8vb2rjq1bt8rrr09wbG1tq13n5uam72/F5L14XzDa+7siv7gck3/hVBURUX378o9zOJtRgIaONvhgUBtOTRlqglNaWoqYmBj07t377y9oaSlvR0VF1eo5Fi1ahKFDh8LRsXrXxp07d8LT0xMtW7bEyy+/jOzs7Bs+R0lJCfLy8qod9E8qK0t88WQ7OWW16+wlrDqQpHRIRERm43DiFcz/87w8/+jRtmjYwFbpkIyaXhOcrKwsVFRUwMvLq9r94nZ6evotHy/W6ogpqtGjR/9jemr58uXYtm0bPv30U/z555/o37+//Fo1mTlzJlxcXKoOf3//u/zOTFczTydM6qvddfzjTaeRkVesdEhERCavtFyNyb8ck92KB4f6ol8bb6VDMnoGXUUlRm/atm2LLl26VLtfjOg88sgj8nODBw/Gxo0bceDAATmqU5OpU6ciNze36khK4sjEzYzs1kQ7VVVSjvfYAJCISO++23VeTk25O9pg+sPsVmzwCY6Hh4dcIJyRkVHtfnFbrJu5mcLCQqxatQqjRo265dcJDg6WXysuTrvq/HpivY6zs3O1g27MytICMx9tKz/+70Q6/jhV/f8fERHVnQtZhZhTWTU17aHWsgkrGXiCY2Njg7CwMDmVpKNWq+XtiIiImz52zZo1cu3Ms88+e8uvk5ycLNfg+PhwtXldER0zR/doIs+n/XqCe1UREemBKOb4z7rjcoqqR3MPDAr1VTokk6H3KSpRIr5gwQIsW7YMp0+flguCxeiMqKoShg8fLqeQapqeEtNPDRs2rHZ/QUEB3nrrLezbtw8JCQkyWRo0aBCaNWsmy8+p7rzWqwX83e2RmluML7acVTocIiKTs/ZQCvaez4atyhIfDmbVVF3S+65dQ4YMwaVLlzBt2jS5sDg0NBSRkZFVC48TExNlZdW1YmNjsXv3bmzZsuUfzyemvI4dOyYTppycHPj6+qJPnz744IMP5FQU1R17Gyt8OLgtnl8cjaV7L2BwB1+083NVOiwiIpNwubAUH246Jc9FH7LAhtWrhenuWGjMsNmJKBMX1VRiwTHX49zaq6sO49cjqXI7hw0TuslyciIiujuioZ8YwQnxdsJvr3SHNX+31un7N19NuqV3H2oNF3trnErLw/Koi0qHQ0Rk9PbFZ8vkRsxIzXysLZMbPeArSrfk0cAWk/uFyPP/23qWO44TEd2F8gp1VQuOp7sEoEMAO/HrAxMcqpUhnf3RtrGL7I0zK/KM0uEQERmtFfsu4kx6vtwp/M0+2saqVPeY4FCtiJ447w/SNp9aE5MsW4oTEdHtySoowRdbtVWpIrlhzxv9YYJDtdYxwA1PhPlV7TiuFj3FiYio1j6LjJUbGt/j64xhXQKUDsekMcGh2yLW4jjZqnAsORc/HeSWF0REtXUkKQerK39vzhh0jxwZJ/1hgkO3pZGTLV7t3Vyez/o9FrlFZUqHRERk8MSI9/RfT8jzxzo0Rligu9IhmTwmOHTbnu8ahOaeDWSTqtlbY5UOh4jI4K2JScLR5Fw0sFVhSn9tVSrpFxMcum2iX8P7j2gXHK/Yn4i4zHylQyIiMlgFJeX47HftwuJXezWHp7Od0iGZBSY4dEe6NvPAv1p7oUKtwcebWTZORHQj3+48L6ungho6yBFwqh9McOiOTe0fApWlBbafycSeuCylwyEiMjipOVex4K94eT6lfyvYqPi2W1/4StMdC27UAM/eGyjPP9x0Wo7mEBHR3z7/PRYl5Wp0CXJH33u0m0xT/WCCQ3dFzCc726lwOi0Pv8QkKx0OEZHBOJacg7WHU+T5Ow+1goXYeIrqDRMcuiuiC+e/e2nLxj/fEovCknKlQyIiUpxGo5Ej28KjHRqjnZ+r0iGZHSY4dNeeiwhEgLsDMvNL8N0u7VwzEZE523IqA9EXLsNWZYm3+nK/KSUwwaG7ZquyqurrMH/XeaTnFisdEhGRYkrL1Zi5WTt6M6ZHMHxd7ZUOySwxwaE60b+NNzoFuqG4TI0v/9D2eyAiMkerDiQiIbsIHg1s8dIDTZUOx2wxwaE6IRbPTR3QSp6LPariMguUDomIqN6JdYhztp2T52JbG9G5mJTBBIfqTFigm2z+J6rFRWkkEZG5Wbz7ArIKShHY0AFDO/srHY5ZY4JDdUosphMb5EaeTJc75xIRmQuxP9/8ykKLN/q0lNvakHL46lOdauHlhMc6+snzT/93RpZKEhGZg292xMl9p+7xdcZDbX2UDsfsMcGhOvda7+awsbJEVHw2/jrHLRyIyPSl5FzF8qiL8nxSvxBYiqFsUhQTHKpzfm4OsjeO8GnkGai5hQMRmbgvt55FaYUa9wa7477mHkqHQ0xwSF/GP9hMVg+cTM3DxuNpSodDRKQ35zLy8csh7VY1k/uFcEsGA8EEh/TC3dEGY+8LludfbIlFWYVa6ZCIiPRCbFMjBqrFZpodAtyUDocqMcEhvRnVvQkaOtrgYnYR1lb+dUNEZEqOJ+fi95MZsnr0zT7cksGQMMEhvXG0VeHlyi6ec7bFyfblRESmZPZWbc+vQaGN0dzLSelw6BpMcEivnr03EJ5OtrLCYPXBJKXDISKqMzEXr2BH7CVYWVrg1V7NlQ6HrsMEh/TKztpKLjgW5m6PQ3FZhdIhERHVif/bqt137/GOjRHk4ah0OHQdJjikd0M6+8PHxQ7pecX4MTpR6XCIiO7a/vhs7I7LgsrSAq/05OiNIWKCQ/UyijOhZ+Uozo7zuFrKURwiMl6iQ/vsytGbpzr7w9/dQemQSKkEZ+7cuQgKCoKdnR3Cw8MRHR19w2uXLl0qewhce4jHXf/DNW3aNPj4+MDe3h69e/fGuXPa3VvJMD0Z5g8/N3tkFZTg+30JSodDRHTH9p7Pxv4Ll2XH9gmVU/BkhgnO6tWrMXHiREyfPh2HDh1C+/bt0bdvX2RmZt7wMc7OzkhLS6s6Ll7Utr/WmTVrFubMmYNvv/0W+/fvh6Ojo3zO4uJifX87dIdsVJb4d+UivG//jJf7tRARGRvxB7bo7SU8HR4AX1d7pUMipRKc2bNnY8yYMRg5ciRat24tkxIHBwcsXrz4ho8Rozbe3t5Vh5eXV7Ufri+//BLvvPMOBg0ahHbt2mH58uVITU3F+vXr9f3t0F14rENjBDV0kDvuLtvLURwiMj5/nr2EQ4k5sFVZYlxlGwwywwSntLQUMTExcgqp6gtaWsrbUVFRN3xcQUEBAgMD4e/vL5OYkydPVn3uwoULSE9Pr/acLi4ucurrRs9ZUlKCvLy8agfVP5WVJV7trR3FWfhXPAo5ikNERkT8gf3fbef+boHhXH35BJlRgpOVlYWKiopqIzCCuC2SlJq0bNlSju78+uuvWLFiBdRqNbp27YrkZG0nXN3jbuc5Z86cKZMg3SESJ1LGw+185SjOlaIyrNhXfeqRiMiQ7YnLxuHK0ZsXK7eiIcNlcFVUERERGD58OEJDQ3H//fdj7dq1aNSoEebPn3/Hzzl16lTk5uZWHUlJbDin5CiOri/Od7viWVFFREY0eqOtnBrWJYCjN+ae4Hh4eMDKygoZGRnV7he3xdqa2rC2tkaHDh0QFxcnb+sedzvPaWtrKxcuX3uQcgZ3aAx/d3tkF5Zi5X6O4hCR4dsXfxkHEq7IyqmX7ufaG5h7gmNjY4OwsDBs27at6j4x5SRui5Ga2hBTXMePH5cl4UKTJk1kInPtc4o1NaKaqrbPScqyFqM4D2hHcebvimd3YyIyeHMq196IxqXeLhy9MQZ6n6ISJeILFizAsmXLcPr0abz88ssoLCyUVVWCmI4SU0g6M2bMwJYtWxAfHy/Lyp999llZJj569OiqCqvXXnsNH374ITZs2CCTH/Ecvr6+GDx4sL6/Haojj3X0Q2NXe1zKL8EqdjcmIgMWfeEyouKzYW1lUbWBMBk+lb6/wJAhQ3Dp0iXZmE8sAhZrayIjI6sWCScmJsrKKp0rV67IsnJxrZubmxwB2rt3rywx15k0aZJMksaOHYucnBx0795dPuf1DQHJcIm+OOIXxTvrT2Den+cxtEuA7HhMRGRovtquHb15spM/+94YEQuNWDllZsSUlqimEguOuR5HOSXlFbh/1k65R9UHg9vguXsDlQ6JiOgfO4Y/Pm+v3HNqx5sPcFsGI3r/NrgqKjIftiqrquHeeTviUFquVjokIqIaR28e7+jH5MbIMMEhRYkFe42cbJGaW4z1R1KUDoeIqMqJlFzsjL0ESwtg3INce2NsmOCQosS6mzE9msjzeTvPo0JtdjOmRGSgvtmpbU/ySHtfBDZ0VDocuk1McEhxT4cHwsXeGheyCvG/E2lKh0NEhLjMAvzvhLY7/suVbS3IuDDBIcU1sFVhRNcgeT53x3nZMZSISEnf/il+FwH/au2Flt5OSodDd4AJDhmEkd2C4GBjhdNpeXLOm4hIKclXirD+sHZNIHcMN15McMgguDrYyN15ha93xHEUh4gUs2BXPMrVGnRr1hAdAtyUDofuEBMcMhijuzeR+7yIvhOicygRUX2T3dUPaDdk1m0pQ8aJCQ4ZDLE775Od/OT53J3nlQ6HiMzQ4j0XUFKuRqi/KyKaNlQ6HLoLTHDIoLx4X1NYWVpg19lLOJ6cq3Q4RGRGcq+W4fuoi/J8/IPN5N6HZLyY4JBBCWjoIHtOXNuDgoioPnwflYCCknK09HJCrxBPpcOhu8QEhwzOS/drqxYiT6Yj/lKB0uEQkRkoLqvAkj0J8vylB4JhKdoXk1FjgkMGR/Sc6BniKXtQLPjrgtLhEJEZ+DkmGdmFpWjsao+H2mlHkcm4McEhgx7F+eVQMjLzi5UOh4hMmNgiZsFf8fJ8dI8msLbiW6Mp4P9FMkidg9zQMcBV7jCuGzYmItIHsUXMxewiuDlYyw2AyTQwwSGDJKoXdKM4K/ZdRH5xmdIhEZEJEk1FxbYMwvAI0VFdpXRIVEeY4JDB6t3KC00bOSK/uBw/7E9UOhwiMkF74rJxIiUPdtaWeL5yTzwyDUxwyGCJKgbRF+fv5lsVSodERCZm/i7t6M2QTv5wd7RROhyqQ0xwyKAN6uALL2dbZOSV4NfDqUqHQ0Qm5ERKLv46lyWbi47uEax0OFTHmOCQQbNVWWFU9yby/Ntd56FWcxNOIqoburU3D7Xzgb+7g9LhUB1jgkMGb1iXADjZqRB/qRB/nM5QOhwiMgFJl4uw+XiaPNdNhZNpYYJDBs/JzhrPhAfK84Vs/EdEdWDR7gsQA8I9mnugta+z0uGQHjDBIaMwomsQrK0sEJ1wGYcTrygdDhEZsZyiUqw+kCTPx97HtTemigkOGQVvFzs80r6xPNd1HCUiuhMr9yfialkFQryd0L2Zh9LhkJ4wwSGjMeY+7WLjyBPpuJhdqHQ4RGSERLsJXXd0MXojmoqSaWKCQ0YjxNsZ97doJOfNF+/mWhwiun2i3URWQQm8ne3wcHtuqmnKmOCQUdHNl/90MBlXCkuVDoeIjIhoM/Fd5RT3C93Fuj6+BZoy/t8lo9K1aUO09nGW8+cr919UOhwiMiJ/nr2EuMwCNLBVYWiXAKXDIT1jgkNGRcyX60Zxlu69iOIybt9ARLXz3S7t6M2wLv5wtrNWOhzSMyY4ZHQGtvOBj4udnEf/9UiK0uEQkRE4npyLqPhsqCwtMLKbtmCBTBsTHDI6Yt78hcpfUAv+usDtG4jolnTtJcS2DL6u9kqHQ6aS4MydOxdBQUGws7NDeHg4oqOjb3jtggUL0KNHD7i5ucmjd+/e/7h+xIgRcqri2qNfv3718J2QoRjaxV/Oo4v59D/PXVI6HCIyYKk5V7GpclsGbqppPvSe4KxevRoTJ07E9OnTcejQIbRv3x59+/ZFZmZmjdfv3LkTw4YNw44dOxAVFQV/f3/06dMHKSnVpyJEQpOWllZ1/Pjjj/r+VsjAtm8Y0tlfni/i9g1EdBPL9iagQq3BvcHuaNPYRelwyFQSnNmzZ2PMmDEYOXIkWrdujW+//RYODg5YvHhxjdevXLkS48aNQ2hoKEJCQrBw4UKo1Wps27at2nW2trbw9vauOsRoD5nf9g2WFsDuuCycSc9TOhwiMkCFJeX4ITpRno/uztEbc6LXBKe0tBQxMTFymqnqC1payttidKY2ioqKUFZWBnd393+M9Hh6eqJly5Z4+eWXkZ2dfcPnKCkpQV5eXrWDjJ+/uwP6t/GR5xzFIaKarDmYhPzicjTxcETPEE+lwyFTSXCysrJQUVEBLy+vaveL2+np6bV6jsmTJ8PX17dakiSmp5YvXy5HdT799FP8+eef6N+/v/xaNZk5cyZcXFyqDjHtRaZhVA/tYuNfj6QiM79Y6XCIyICIaanFldsyvNC9CSzFkC+ZDYOuovrkk0+watUqrFu3Ti5Q1hk6dCgeeeQRtG3bFoMHD8bGjRtx4MABOapTk6lTpyI3N7fqSErS7iJLxq9jgBs6BriitEKNFVFs/EdEf9t6KgOJl4vg6mCNJzr6KR0OmVKC4+HhASsrK2RkZFS7X9wW62Zu5vPPP5cJzpYtW9CuXbubXhscHCy/VlxcXI2fF+t1nJ2dqx1kOnRVEd/vY+M/Ivrbot3a0vBnwwNhb2OldDhkSgmOjY0NwsLCqi0Q1i0YjoiIuOHjZs2ahQ8++ACRkZHo1KnTLb9OcnKyXIPj46Ndj0HmpU9rL/i52eNKURnWHmLjPyICjiTl4EDCFVhbWWB4RKDS4ZApTlGJEnHR22bZsmU4ffq0XBBcWFgoq6qE4cOHyykkHbGm5t1335VVVqJ3jlirI46CggL5efHxrbfewr59+5CQkCCTpUGDBqFZs2ay/JzMj8rKsqozqfiLjY3/iGjRbm3hwSPtG8PT+e8lDmQ+9J7gDBkyRE43TZs2TZZ+HzlyRI7M6BYeJyYmyj42OvPmzZPVV0888YQckdEd4jkEMeV17NgxuQanRYsWGDVqlBwl+uuvv+RUFJmnpzr5wclWhfOXCuWGekRk3o39Nlc29hvVndsymCsLjUZjdn/uijJxUU0lFhxzPY7p+HDjKSzcfQHdm3lgxehwpcMhIoXM3Hwa83fFo2vThvhhzL1Kh0MKvX8bdBUV0e14/prGf7Hp+UqHQ0QKNfb7sbKxH0dvzBsTHDKpxn9979FW5y3Zw8Z/ROZo7aFk5BWXI6ihAx5sycZ+5owJDpkU3V9saw+nILugROlwiKgeqa9p7CcKD9jYz7wxwSGTEhbohnZ+LigtV+OH/dphaiIyDzvPZuJCViGc7FR4IoyN/cwdExwyKRYWFnihsmR8+b6LKCln4z8icysNH9YlAI62KqXDIYUxwSGTM6CtD7ycbXEpvwSbjv3dgoCITNeZ9DzsicuWhQZs7EcCExwyOTYqSwyPCKr6i84MOyEQmZ3FlaM3/dv4wM/NQelwyAAwwSGT9HSXANiqLHEyNQ/RFy4rHQ4R6VFWQQnWH0mV5y901/5xQ8QEh0ySm6MNHqvcPVg3L09EpmnlvkRZWNDe3xUdA9yUDocMBBMcMlmjKv+S23o6A4nZRUqHQ0R6IAoJvt93UZ6/0C1IFhoQCUxwyGQ183RCj+YeEEtwlkVpe2MQkWkRhQRiikoUFogCAyIdJjhk0l6obPz304EkFJSUKx0OEdUhUUCwpLKxnygssLbiWxr9jT8NZNLub94IwR6OyC8pxy8xyUqHQ0R1KObiFRxPyZUFBaL3DdG1mOCQSROt2kd0067FWbo3QbZyJyLTsLhyz7nBoY3h7mijdDhkYJjgkMl7vKOfbN0uWriLVu5EZPxScq7i95MZ8nwkS8OpBkxwyOSJlu1DO/vLc918PREZt+VRCahQa9C1aUOEeDsrHQ4ZICY4ZBbEAkTRwv2vc1k4m5GvdDhEdBeKSsuxKjqpatdwopowwSGz4O/ugD6tveU5R3GIjNvaQynIvVqGwIYO6BniqXQ4ZKCY4JDZGFm52Hjd4WTkFJUqHQ4R3WFpuCgYEJ6PCIKVGJolqgETHDIbXZq4o7WPM4rL1PixcnibiIyLmGaOyyxAA1sVnuyk3Y6FqCZMcMhsiBbuulGc76MSUF6hVjokIrpNSypLw58IE9WR1kqHQwaMCQ6ZlYfb+6Khow1Sc4ux5ZS2xJSIjINo9bAj9hLEdlMjurI0nG6OCQ6ZFTtrKzwdHlDtL0EiMg7LKtfe9GzpiSAPR6XDIQPHBIfMzrP3BkJlaYEDCVdwIiVX6XCIqBbyi8uw5qB27ZyuOznRzTDBIbPj5WyHge20uw6zZJzIOKw5mIzC0go092yA7s08lA6HjAATHDJLuvn7346m4lJ+idLhENFNiI7Fy6ISqkZvRMEA0a0wwSGz1CHADaH+riitECXjiUqHQ0Q3sTM2Exezi+Bsp8KjHRorHQ4ZCSY4ZLaqSsb3XURpOUvGiQyVbip5WJcAONiolA6HjAQTHDJb/dv4wNPJVk5RbT6epnQ4RFQDsXfc7rgsuZfccxGBSodDRoQJDpktG5UlnrtX+wuTJeNEhkm3LYPYS87PzUHpcMiIMMEhszYsPEAmOkeTc3Eo8YrS4RDRNcSecWsPJVebUiYyqARn7ty5CAoKgp2dHcLDwxEdHX3T69esWYOQkBB5fdu2bbF58+Z/bLY2bdo0+Pj4wN7eHr1798a5c+f0/F2QKfJoYItH2vvK86UsGScyKKsPJMm948QecmIvOSKDSnBWr16NiRMnYvr06Th06BDat2+Pvn37IjMzs8br9+7di2HDhmHUqFE4fPgwBg8eLI8TJ05UXTNr1izMmTMH3377Lfbv3w9HR0f5nMXFxfr+dsiES8bFOpz0XP4MERkCsVfc8qiL8pyl4XQnLDRiOESPxIhN586d8fXXX8vbarUa/v7+eOWVVzBlypR/XD9kyBAUFhZi48aNVffde++9CA0NlQmNCNfX1xdvvPEG3nzzTfn53NxceHl5YenSpRg6dOgtY8rLy4OLi4t8nLOzc51+v2Scnvx2r+xs/ErPZnijT0ulwyEye5En0vDSikNwd7TB3ik95TYrRHm38f6t1xGc0tJSxMTEyCmkqi9oaSlvR0VF1fgYcf+11wtidEZ3/YULF5Cenl7tGvHNikTqRs9ZUlIiX5RrD6JrjezWRH78YX8iissqlA6HyOzpSsOf7hLA5IbuiF4TnKysLFRUVMjRlWuJ2yJJqYm4/2bX6z7eznPOnDlTJkG6Q4wgEV2rT2sv+LrYIbuwVHY3JiLlnEzNxf4Ll+WecWLvOKI7YRZVVFOnTpXDWbojKUm7YRuRjsrKEs9FBFWVpep55paIarFreP+2PvB2sVM6HDJSek1wPDw8YGVlhYyMjGr3i9ve3t41Pkbcf7PrdR9v5zltbW3lXN21B9H1hnb2h63KEidT8+R6HCKqf9kFJVh/JLVaAQCRwSU4NjY2CAsLw7Zt26ruE4uMxe2IiIgaHyPuv/Z6YevWrVXXN2nSRCYy114j1tSIaqobPSdRbbg52lTtc7N0Lxv/ESlh1YEkuXVKOz8XdAxwVTocMmJ6n6ISJeILFizAsmXLcPr0abz88suySmrkyJHy88OHD5dTSDqvvvoqIiMj8cUXX+DMmTN47733cPDgQUyYMEF+XpQKvvbaa/jwww+xYcMGHD9+XD6HqKwS5eREd0OUowq/n8xASs5VpcMhMitlFWp8X1kaLhr7sTSc7obedy0TZd+XLl2SjfnEImBR7i0SGN0i4cTERFlZpdO1a1f88MMPeOedd/D222+jefPmWL9+Pdq0aVN1zaRJk2SSNHbsWOTk5KB79+7yOUVjQKK7EeLtjIjghoiKz8byqARM7d9K6ZCIzEbkiXSk5xXLBpwD2vooHQ4ZOb33wTFE7INDN7PlZDrGfh8DF3tr7JvaC/Y2LFElqg+PfbMHhxJz8Gqv5nj9Xy2UDocMkMH0wSEyRr1aecHf3R65V8uw/kiK0uEQmYWjSTkyubG2ssAz9wYoHQ6ZACY4RNexsrTA85Ul42KXcTMc5CRSbNfwh9v5wtOJyw3o7jHBIarBk5384WBjhbMZBdh7PlvpcIhMWmZeMTYeqywN567hVEeY4BDVQKy/ebyjX7WW8USkHyv3J6KsQoOwQDe082NpONUNJjhEN6D7S3LbmQxczC5UOhwik1RSXoGV+/8uDSeqK0xwiG6gaaMGuL9FI4glOMv2an8BE1Hd2ng0DVkFpfB2tkPfe2ruRk90J5jgEN2E7i/KNQeTUFBSrnQ4RCZFLODXLS5+LiIQ1lZ8S6K6w58mopu4r3kjBHs4Ir+kHL/EJCsdDpFJibl4BcdTcuUecMO6sDSc6hYTHKKbsLS0qFqLI/7SVKtZMk5UV3QL+AeHNoa7o43S4ZCJYYJDdAuimsrJToULWYX48+wlpcMhMgmpOVcReTJdnrM0nPSBCQ7RLTjaqjCkk788X7yHu4wT1YXlURdRodbIvd9a+XDLHKp7THCIauH5rkGwtAD+OpeFcxn5SodDZNSullbgx+hEec7ScNIXJjhEteDv7oB/tfaS50sqqz6I6M6sPZws93oLcHeQe78R6QMTHKJaGtmtify49lAycopKlQ6HyHhLwysXF4uRUbH3G5E+MMEhqqXwJu5yrUBxmRqrDiQpHQ6RUdodl4VzmQVwtLHCk52026EQ6QMTHKJasrCwqFovsHxvAsor1EqHRGR0Fu++ULWhrbOdtdLhkAljgkN0Gx5p74uGjjZIzS3G7yczlA6HyKjEXyrAjthLsLDQTk8R6RMTHKLbYGdthWfCtR1Xl7BknOi2LKtcoN+zpSeaeDgqHQ6ZOCY4RLfp2XvFnjkWOHjxCo4l5ygdDpFREFVTayq3O9Et2CfSJyY4RLfJ09kOD7XzrbaegIhuTmxYW1RagRZeDdCtWUOlwyEzwASH6A7oFhtvOp6GjLxipcMhMmhiQb5u3ykxeiMW7BPpGxMcojvQzs8VnYPcUFahwfdRF5UOh8igbT2VgZScq3BzsMajHRorHQ6ZCSY4RHfohcp1BCv3X0RxWYXS4RAZLN0ebs+EB8qF+kT1gQkO0R0SWzc0drXHlaIyrD+conQ4RAZJLMQ/kHBFLsx/LiJQ6XDIjDDBIbpDKivLqrU44i9U0YKeiKrTLcQXC/O9nO2UDofMCBMcorvwVGd/2XL+bEaBbEFPRH8TC/A3HkurNqVLVF+Y4BDdBdFqXrScFxaxZJyomuVRCShXa9AlyB1t/VyUDofMDBMcors0omuQbD2/M/YS4jILlA6HyCBcLa3AD/sT5fkL3bktA9U/JjhEdynIwxG9Qrzk+dK9HMUhEtYfSZEL8P3c7PGv1t5Kh0NmiAkOUR0Y1V27vuDnmGTkFJUqHQ6RosSCe93iYjHCaWXJxn5U/5jgENWBe4Pd0crHGcVlavwQrR2WJzJXu85l4VxmARrYqjCks3aNGpFJJTiXL1/GM888A2dnZ7i6umLUqFEoKCi46fWvvPIKWrZsCXt7ewQEBODf//43cnNzq10n2nxff6xatUqf3wrRTYmfwdGVozhix+TScrXSIREpZuFf8fKjSG6c7KyVDofMlF4THJHcnDx5Elu3bsXGjRuxa9cujB079obXp6amyuPzzz/HiRMnsHTpUkRGRsrE6HpLlixBWlpa1TF48GB9fitEt/Rwe194OtkiI68Em46nKh0OkSJi0/Px17ksiFkpMT1FpBSVvp749OnTMjk5cOAAOnXqJO/76quvMGDAAJnA+Ppqd2O+Vps2bfDLL79U3W7atCk++ugjPPvssygvL4dK9Xe4YkTI25sL18hw2Kgs8XzXIHz2eywW/nUBg0Mbc1NBMjuLdmtHb/q38YG/u4PS4ZAZ09sITlRUlExCdMmN0Lt3b1haWmL//v21fh4xPSWmuK5NboTx48fDw8MDXbp0weLFi9lFlgzC010CYGdtiZOpedgXf1npcIjq1aX8Eqw/rB29HNWDjf3IREdw0tPT4enpWf2LqVRwd3eXn6uNrKwsfPDBB/+Y1poxYwZ69uwJBwcHbNmyBePGjZNre8R6nZqUlJTIQycvL++OvieiW3FztMETYX5YsS9R/iUb0bSh0iER1Zvv911EaYUaHQNc0THATelwyMzd9gjOlClTalzke+1x5syZuw5MJCEDBw5E69at8d5771X73Lvvvotu3bqhQ4cOmDx5MiZNmoTPPvvshs81c+ZMuLi4VB3+/lzVT/qja0n/x+lMxF9i4z8yD8VlFVix76I8H90jWOlwiG4/wXnjjTfk+pqbHcHBwXJ9TGZmZrXHinU0olLqVmtn8vPz0a9fPzg5OWHdunWwtr75Kvzw8HAkJydXG6W51tSpU+VUl+5ISkq63W+bqNaCGzVA71aeVZtwEpmDtYdScLmwVDb269Na2/iSyKimqBo1aiSPW4mIiEBOTg5iYmIQFhYm79u+fTvUarVMSG42ctO3b1/Y2tpiw4YNsLO79e6zR44cgZubm3xMTcT9N/ockT6M6h4sR3BE4783/tVSTl0RmSq1WlO1uHhktyZQWbHFGilPbz+FrVq1kqMwY8aMQXR0NPbs2YMJEyZg6NChVRVUKSkpCAkJkZ/XJTd9+vRBYWEhFi1aJG+L9TriqKiokNf89ttvWLhwoSwjj4uLw7x58/Dxxx/L/jlEhtT4r01jbeO/lfu1w/ZEpurPs5dw/lIhnGxVeKqTn9LhEEl6TbNXrlwpE5hevXrJ8vDu3bvju+++q/p8WVkZYmNjUVRUJG8fOnRIVlgdP34czZo1g4+PT9Whm1YS01Vz586VI0ShoaGYP38+Zs+ejenTp+vzWyG6g8Z/2nUIy6IuoqRcm6ATmaIFbOxHBshCY4b11WJkSCw21pWgE+lDWYUaPT7dgfS8Ysx6oh2e6sTF7WR6TqTk4qGvdsv9pnZNehCNXe2VDolMWN5tvH9zopRIT6ytLPFCd20n1wW74tmriUx69Oahdj5MbsigMMEh0qOhXQLkhoNi48GdsZeUDoeoTqXkXMXGY2nyfAxLw8nAMMEh0iNnO2sM66Kdmvpul/YvXSJTsWT3BVSoNejatCHaNHZROhyiapjgEOmZLJu1tEBUfDaOJ+cqHQ5Rnci9WoYfoxPl+dj7OHpDhocJDpGe+bray53Gr12vQGTsVkUnorC0Ai29nHB/i1v3RiOqb0xwiOrB6MqNBzcdT0PyFW1bBCJjVVquxpI9CVU/26ItApGhYYJDVA/u8XVB92Yecr3C4t3aNwYiY/Xb0VTZ/sDTyRaPhGpHJ4kMDRMconoypnKdwqoDicgtKlM6HKI7Itod6KZaR3QLgq3KSumQiGrEBIeontzX3AMh3k4oKq3Aymhu30DGuy3DmfR8ONhY4ZkugUqHQ3RDTHCI6olYp6CrNhHTVMVl3L6BjM/8P7WjN8O6BMDFgdsykOFigkNUj0Q1la+LHbIKSrD2UIrS4RDdlqNJObLdgWh7MKq7duE8kaFigkNUz9s3jKrs+PrdrvNy0TGRsfj2z/Pyo1hYLNofEBkyJjhE9WxoZ3+42FsjIbsIW06mKx0OUa3EXypAZOXP60v3N1U6HKJbYoJDVM8cbVV4PiKw6i9ibsJJxmDBXxcgflR7hXiihZeT0uEQ3RITHCIFDO8qymstcTQ5V65pIDJkmfnF+OVQsjx/kaM3ZCSY4BApwKOBLZ7q5F+tKoXIUImuxaJ7cccAV3QOclM6HKJaYYJDpJAxPYJhaaHtK3IqNU/pcIhqlF9chhX7LlatveG2DGQsmOAQKSSgoQMGttO2uZ+/S1udQmRoftifiPzicjRt5IjerbyUDoeo1pjgECnoxcrGfxuPpSExm5twkmERzSgX7b4gz1+8ryksxZAjkZFggkOkoDaNXXB/i0ayHw5HccjQiIXFmfkl8HGxw+AOjZUOh+i2MMEhUtj4B5vJj2sOJiMzr1jpcIik8gp1VWM/scWIjYpvF2Rc+BNLpLAuTdxlZUpphbpql2YipYlp06TLV+HuaIOhnQOUDofotjHBITIA4ypHcVbuT8SVwlKlwyEzp1Zr8M3OOHku9pyyt7FSOiSi28YEh8gAPNCiEVr7OKOotAJL9yYoHQ6ZuT9OZ+BsRgGcbFV49l5t120iY8MEh8gAiN4iurU4IsEpKClXOiQyU2LrkLk7tWtvnosIlPumERkjJjhEBqJfG28Eezgi92oZftivbaxGVN/2ns/G0aQcuZXIC92bKB0O0R1jgkNkIKwsLfDSA02rNjYUPUiI6tvcHdq1N8O6BMgtRYiMFRMcIgMyOLQxfF3scCm/BGsOJikdDpmZQ4lX5AiOytICYyqbUBIZKyY4RAZE9BrR7dY8b+d5ucEhUX2Zs+2c/PhYx8Zo7GqvdDhEd4UJDpGBGdLZH55OtkjNLZadZInqg1h3szP2kpwq1S14JzJmTHCIDIydtVXVKI5YD1FWwVEc0r+vtmtHbwaF+iKwoaPS4RAZdoJz+fJlPPPMM3B2doarqytGjRqFgoKCmz7mgQcekCWz1x4vvfRStWsSExMxcOBAODg4wNPTE2+99RbKy1lWS6bjabnA0wbJV65i3eEUpcMhE3ciJRd/nM6E2EuTozdkKvSa4Ijk5uTJk9i6dSs2btyIXbt2YezYsbd83JgxY5CWllZ1zJo1q+pzFRUVMrkpLS3F3r17sWzZMixduhTTpk3T57dCVK9E51ix/49uFEfsC0Sk79Gbh9v7ommjBkqHQ2TYCc7p06cRGRmJhQsXIjw8HN27d8dXX32FVatWITU19aaPFSMz3t7eVYcYAdLZsmULTp06hRUrViA0NBT9+/fHBx98gLlz58qkh8hUiA6yYh+gi9lF2HD05v9miO7U6bQ8/H4yAxYWwCs9OXpDpkNvCU5UVJSclurUqVPVfb1794alpSX2799/08euXLkSHh4eaNOmDaZOnYqioqJqz9u2bVt4eXlV3de3b1/k5eXJ0SIiU+Fgo8KYHtpRnK+3x6FCrVE6JDJB4mdLGNjWB808nZQOh6jOqKAn6enpcn1MtS+mUsHd3V1+7kaefvppBAYGwtfXF8eOHcPkyZMRGxuLtWvXVj3vtcmNoLt9o+ctKSmRh45IhoiMgWiVP3/XecRnFWLjsVQMCm2sdEhkQs5m5GPziTR5/krP5kqHQ6TsCM6UKVP+sQj4+uPMmTN3HJBYoyNGZMQojVjDs3z5cqxbtw7nz2v3RrkTM2fOhIuLS9Xh7+9/x89FVJ8a2KowurJd/lccxaE6Jn6mNBqgfxtvtPTm6A2ZeYLzxhtvyPU1NzuCg4Pl2pnMzMxqjxWVTqKySnyutsT6HSEuTjuMKh6bkZFR7Rrd7Rs9r5jmys3NrTqSktghlozH8K5BcLZTIS6zQI7iENWF2PT8qp+nCVx7QybotqeoGjVqJI9biYiIQE5ODmJiYhAWFibv2759O9RqdVXSUhtHjhyRH318fKqe96OPPpLJk24KTFRpiYXIrVu3rvE5bG1t5UFkjJztrGVF1edbzuK/f5yTayVUVmxhRXfnyz/OVo3e3OPronQ4RHVOb78lW7VqhX79+smS7+joaOzZswcTJkzA0KFD5foaISUlBSEhIfLzgpiGEhVRIilKSEjAhg0bMHz4cNx3331o166dvKZPnz4ykXnuuedw9OhR/P7773jnnXcwfvx4JjFkskZ0awI3B2u5Fmf9EY7i0N05mZqL/51Il5VTr/VuoXQ4RHqh1z8DRTWUSGB69eqFAQMGyFLx7777rurzZWVlcgGxrkrKxsYGf/zxh0xixOPEdNjjjz+O3377reoxVlZWsqeO+ChGc5599lmZBM2YMUOf3wqR4mtxdN2NxX5B7G5Md+P/tp6VHx9q58u1N2SyLDQaMUhpXkQVlVhsLNbjXNtjh8iQFZWW475ZO5BVUIpPHmuLoV0ClA6JjHTPqUFz98iuxVsn3s/GfmSy79+cyCcyor44Lz/QrKr6paS8QumQyAjNrhy9GdyhMZMbMmlMcIiMyDPhAfBytkVKzlX8dIDVgHR7Yi5exp9ntTuGv9qLfW/ItDHBITKyncYnVG6G+PWOOBSXcRSHau+LLdrRmyfD/LhjOJk8JjhERuapzv7wdbFDRl4JVuy7qHQ4ZCSizmdj7/lsWFtZsO8NmQUmOERGxlZlhX9XTi98s/M88ovLlA6JDJyoJfk0UtthfmjnAPi5OSgdEpHeMcEhMkJPhPkhuJEjLheWYsFfF5QOhwzcllMZOJKUA3trK7zSi6M3ZB6Y4BAZIdHJ+K0+LeX5wr/ikVXw92ayRNcqr1Djs99j5fmo7k3g6WSndEhE9YIJDpGR6tfGG+39XFBUWoGvt2v3aiO63trDKXIfM1cHa4y9P1jpcIjqDRMcIiNlYWGByf1C5PnK/ReRdFnbEZxIR1TZfVnZ92b8A83kvmZE5oIJDpER69rMAz2ae6CsQlPVwI1IR1TZpeYWw8fFDs9FBCodDlG9YoJDZOQm9dWO4qw/koLTaXlKh0MGIq+4TPZKEl7v3UL2UCIyJ0xwiIxcWz8XDGznA7GrnG4xKdGCXfHIKSpD00aOeKxjY6XDIap3THCITMCbfVrK9vvbz2TKhm5k3jLyirGwsn3AW31byqo7InPDn3oiE9DEwxFPV+4u/uGmU1CrNUqHRAr6/PdYXC2rQFigG/re4610OESKYIJDZCJe690cTrYqnEzNk6XBZJ5Opubi50PJ8vw/A1vJajsic8QEh8hENGxgi/GVewzJv+BLuRGnOW7J8NGm03I91sPtfdExwE3pkIgUwwSHyISM6BqExq72SM8rxoK/4pUOh+qZWIMlNtS0UVliUl9tp2sic8UEh8iEiFLgKf21ZePf/nkemXnFSodE9aSsQo2PN5+W5y90awJ/d26oSeaNCQ6RiXmonQ86BLjKLRy+2MLmf+bix+hEnL9UCHdHG4x7sKnS4RApjgkOkYkRi0rfGdhKnv8Uk8Tmf2bS1O/LP87J89d7N+eWDERMcIhMU1ige1Xzvxm/nZKLT8l0fbXtHC4XlqKZZwMMq2wXQGTumOAQmagp/UJgq7JEVHw2Nh9PVzoc0pO4zHws2ZMgz8XIHZv6EWnxXwKRiRKLTF+6X7sW46NNp1BUWq50SFTHxMjcextOoVytQe9WXnigpafSIREZDCY4RCbs5Qeaws/NXu4o/c2O80qHQ3Us8kQ6dsdlybLwaQ+1VjocIoPCBIfIxMvG3xmofeP7blc8ErIKlQ6J6oho5PjhJm1Z+Ev3BSOgIcvCia7FBIfIxPW9xws9mnugtEKNDzaeUjocqiPzdsYhJeeqbOz48gPaDtZE9DcmOERmUDY+/eF7oLK0wLYzmdh+JkPpkOguJWYX4dtd8VULi+1trJQOicjgMMEhMgOifPiF7k3k+fu/nUJxGfepMmYzNp5Cabka3Zt5oF8b7hZOVBMmOERm4pWezeDpZIuL2UX4ZicXHBur30+m44/TGXJE7r1HWnO3cKIbYIJDZCac7Kwx7eHWVes3RP8UMi75xWWY/utJeT72vmA083RSOiQig8UEh8iMDGzrg54hniir0ODttSegVrPDsTERe4uJneIDGzrg372aKx0OkfkmOJcvX8YzzzwDZ2dnuLq6YtSoUSgoKLjh9QkJCXK4taZjzZo1VdfV9PlVq1bp81shMgni38qMQffA3toK0QmX8dPBJKVDolo6nHgFy6K0HYs/GtxWtgAgIoUSHJHcnDx5Elu3bsXGjRuxa9cujB079obX+/v7Iy0trdrx/vvvo0GDBujfv3+1a5csWVLtusGDB+vzWyEyGX5uDnijTwt5/vHm08jML1Y6JLqFsgo1pq49LvcWe6xDY3Rv7qF0SEQGT6WvJz59+jQiIyNx4MABdOrUSd731VdfYcCAAfj888/h6+v7j8dYWVnB27t6RcC6devw1FNPySTnWmJE6Ppriah2RnQNwvojKTiRkocPNp7GV8M6KB0S3cTCvy7gTHo+3Bys8Z/KneKJSKERnKioKJmE6JIboXfv3rC0tMT+/ftr9RwxMTE4cuSInNq63vjx4+Hh4YEuXbpg8eLF3C2Z6DaIDRk/eawdLC2A346mYkdsptIh0U163vx321l5/vaAVmjYwFbpkIjMO8FJT0+Hp2f1jd9UKhXc3d3l52pj0aJFaNWqFbp27Vrt/hkzZuCnn36SU1+PP/44xo0bJ0eHbqSkpAR5eXnVDiJz16axC17opu2N85+1x5FXXKZ0SHQdsQh80i9HUVymRkRwQzwR5qd0SESmm+BMmTLlhguBdceZM2fuOrCrV6/ihx9+qHH05t1330W3bt3QoUMHTJ48GZMmTcJnn312w+eaOXMmXFxcqg6x1oeIgIl9WiDA3UFuxvkht3EwOMujErAv/rJcFD7zsbbseUOkzwTnjTfekOtrbnYEBwfL9TGZmdWHvcvLy2VlVW3Wzvz8888oKirC8OHDb3lteHg4kpOT5UhNTaZOnYrc3NyqIymJlSNEgoONCp8/2R7iffOng8ncxsGAXMgqxCeR2j8Wpw4IQZCHo9IhEZn2IuNGjRrJ41YiIiKQk5Mj19GEhYXJ+7Zv3w61Wi0TktpMTz3yyCO1+lpinY6bmxtsbWuemxb33+hzROauSxN3jOrWBAt3X8CUX45jy+tucHWwUToss1ah1uDNNdqpqa5NG+LZ8EClQyIyOnpbgyPWzvTr1w9jxoxBdHQ09uzZgwkTJmDo0KFVFVQpKSkICQmRn79WXFycLCkfPXr0P573t99+w8KFC3HixAl53bx58/Dxxx/jlVde0de3QmTy3uzbEk0bOSIzvwTTN2g75ZJyFu2OR8zFK2hgq8KsJ9rBUqwGJyLD6YOzcuVKmcD06tVLlod3794d3333XdXny8rKEBsbK6eiriWqovz8/NCnT59/PKe1tTXmzp0rR4hCQ0Mxf/58zJ49G9OnT9fnt0Jk0kTTODFVJd5Hfz2SisgTaUqHZLbOZeTj8y1nq3YKF32LiOj2WWjMsL5aVFGJxcZiPY7oskxEWrMiz8iNOBs62uD31++DB0uS672h3xPz9uJoci4eaNkIS0Z05sJiojt8/+ZeVERU5dXezRHi7YTswlK88dNR7lVVz2ZvPSuTG2c7lexTxOSG6M4xwSGiKrYqK/x3aAfYqizx59lLWLg7XumQzMaus5cwb+d5ef7J4+3g7WKndEhERo0JDhFV09LbCdMfvkeez4qMxZGkHKVDMnliP7CJPx2R58+EB2BAWx+lQyIyekxwiOgfhnXxx8C2PihXa/DKj4fY5ViPxDSgmA7MKihFSy8nvPtQa6VDIjIJTHCI6B/E2o+PH2sLPzd7JF2+irflTtZcj6MP83fF469zWbCztsTXT3eQFW1EdPeY4BBRjVzsrTFnWAeoLC2w8VgaVh9gB/C6dijxCj7fEivP33/kHjT3clI6JCKTwQSHiG6oY4CbbAIoTNtwEseSuR6nrlzKL8G4FYdk1+KH2/viqU7cI4+oLjHBIaKbGtsjGL1beaK0XI0Xv4+Rb8x0d8RrOW5lDNLzimUH6Y8fbcOScKI6xgSHiG5KbBMwe0goghs5Ii23GONXHpIN6ejOzdh4EgcSrsDJVoXvhneCk5210iERmRwmOER0S8521lgg3ohtVYhOuIwZv51SOiSj9WN0IlbsS5Q7uP93WCiaNmqgdEhEJokJDhHVingjFm/I4o35+30XsSo6UemQjE7MxcuY9usJef5mn5boGeKldEhEJosJDhHVmnhDfuNfLeT5tF/FNMtlpUMyGqk5V/HSCjG9p8GAtt4Y90BTpUMiMmlMcIjotox/sBn6t/FGaYUao5cdlLtf083lFJVi+OJouUBbNPP77In2XFRMpGdMcIjotog35tlPhaJjgCtyr5bh+cXRSM8tVjosg1VcViETwbjMAng722HJyM5wtFUpHRaRyWOCQ0S3zd7GCoue7ywrq1JzizFiSbRMdqg60ePm3z8exsGLV+Bkp8KyF7rA19Ve6bCIzAITHCK6I26ONlg2sgsaOdniTHo+xi4/KEcrSEtsbTF9wwlsOZUBGytLWYUmNjIlovrBBIeI7pi/uwOWjuyMBrYq7L9wGa+vPoJy9siR5myLqyoH/3JoKO4Nbqh0SERmhQkOEd2Ve3xd8N1zYbC2ssD/TqTj1VVHzL4R4H//OIf/++OsPJ/+UGsMaOujdEhEZocJDhHdta7NPDDvmTA5FbPpeBpe+eGw3I7AHKelZm89W5XcTO4XghHdmigdFpFZYoJDRHWid2svzH8uDDYqS0SeTMf4Hw6ZVZIjkpsvtpzFnG3n5O23B4TgZfa6IVIMExwiqjMPhnjKxbQiydl6KkNuKFlSXmEWyc2nkbH4ekecvP3OwFYYex+TGyIlMcEhojp1f4tGWPR8J9iqLPHH6UwMXxQtG92ZKjFKNeWX4/j2z/Py9nsPt8boHsFKh0Vk9pjgEFGd69G8EZaM+Lu66tFv9iIhqxCmJreoTPYAWn0wCZYWwIeD23DNDZGBYIJDRHpbePzLy13R2NUeF7IKMfibPYi+YDp7VyVmF+GxeXuw93w2HG2ssPD5Tnj23kClwyKiSkxwiEhvRGO7deO7or2fC3KKyvDswv1YeygZxu5gwmWZsJ2/VAgfFzuseakrdwYnMjBMcIhIrzyd7LBqbETVBp0TfzqKt9YcRWFJOYxx64Wvt5/DkO/24XJhKdo2dsGv47uhta+z0qER0XWY4BBRvexdNffpjvh3r+ays++amGQMnPMXjiblwFik5FzFsAX78PmWszLReaS9L1a/eC88ne2UDo2IamChEfWNZiYvLw8uLi7Izc2FszP/8iKqT/vis+WWDmm5xVBZWmBinxZ48b6msBKrdA3UpmNpmLr2GPKKy+V6mxmD2uCxjo3lzupEZJjv30xwmOAQKVJ99Pa647LrsdDOzwXTH26NsEB3GBJR+fXR5tOyp4/Q3t8Vc4aGIrCho9KhEZmlPCY4N8cEh0h54lePmKqa8dspFFSuxxkU6osp/UPg42KvaGz5xWWyad+S3Qly3ZAYXXr5/qZ4tXdzWFtxZp9IKUxwboEJDpHhuJRfgi+2xMpeMuK3kZ21pewCPKJrENwdbeo1luKyCvxyKBn/t/UcsgpK5H33tWiEdwe2QnMvp3qNhYju7v1bb3+KfPTRR+jatSscHBzg6upaq8eIXGvatGnw8fGBvb09evfujXPntPu66Fy+fBnPPPOM/MbE844aNQoFBQV6+i6ISN8aOdnik8fb4bcJ3dElyB3FZWq5n1PEzG2Y/PMxnEnP03sMablXMSvyjPya/1l3QiY3wR6OWDyiE5aN7MzkhsgI6W0EZ/r06TIBSU5OxqJFi5CTc+tqiU8//RQzZ87EsmXL0KRJE7z77rs4fvw4Tp06BTs7baVC//79kZaWhvnz56OsrAwjR45E586d8cMPP9Q6No7gEBkm8eto8/F0ue3B8ZTcqvsjghvi8TA/PNCyETwa2NbJ1xLTYnvisvDb0VT870S6rIwS/NzsMap7EzwTHij31CIiw2FQU1RLly7Fa6+9dssER4Th6+uLN954A2+++aa8T3wDXl5e8jmGDh2K06dPo3Xr1jhw4AA6deokr4mMjMSAAQNkIiUeXxtMcIgMm/h9EHPxCpbsSZA7k+uSD1G01M7PFQ+2bIQezT3QzNMJLvbWtXrOq6UViM8qQNT5bOyMvYT9F7JRVvH3r7/wJu4Y2a0J/tXay6AruojMWd5tvH+rYCAuXLiA9PR0OS2lI76J8PBwREVFyQRHfBSjQrrkRhDXW1paYv/+/Xj00UcVip6I6pIov+4U5C4P0X9mdXQitp3JxMnUPNk7Rxxf/qGdvm7oaIMgD0cENXSEs73qH2tqErKKkJBdKMvSrxfg7oCeIZ54spMf7vF1qbfvj4j0z2ASHJHcCGLE5lritu5z4qOnp2e1z6tUKri7u1ddU5OSkhJ5XJsBEpFxEHtZTezTUh4ZecXYGZuJ7WcycSgxRy5Qzi4slYcY8bkVMdrTprEzHmzpiQdDPOU6G/ayITJNt5XgTJkyRa6TuRkxjRQSEgJDItb1vP/++0qHQUR3ycvZDkM6B8hDt45G9KoRIzTi49WyimrXqywtEdjQQY7wNGnoCLd6rsoiIiNJcMT6mBEjRtz0muDg4DsKxNvbW37MyMiQVVQ64nZoaGjVNZmZmdUeV15eLiurdI+vydSpUzFx4sRqIzj+/v53FCcRGY4Gtiq0aewiDyKiO05wGjVqJA99EFVTIknZtm1bVUIjEhGxtubll1+WtyMiIuRi5ZiYGISFhcn7tm/fDrVaLdfq3Iitra08iIiIyDzorQYyMTERR44ckR8rKirkuTiu7VkjprLWrVsnz8U8uKi2+vDDD7FhwwZZHj58+HBZGTV48GB5TatWrdCvXz+MGTMG0dHR2LNnDyZMmCAXINe2goqIiIhMn94WGYuGfaKfjU6HDh3kxx07duCBBx6Q57GxsbLUS2fSpEkoLCzE2LFj5UhN9+7dZRm4rgeOsHLlSpnU9OrVS1ZPPf7445gzZ46+vg0iIiIyQtyqgX1wiIiIjIJBbNVAREREpBQmOERERGRymOAQERGRyWGCQ0RERCaHCQ4RERGZHCY4REREZHKY4BAREZHJYYJDREREJocJDhEREZkcvW3VYMh0zZtFR0QiIiIyDrr37dpswmCWCU5+fr786O/vr3QoREREdAfv42LLhpsxy72o1Go1UlNT4eTkJHcxr+vsUiROSUlJ3OdKj/g61w++zvWDr3P94Ots/K+1SFlEcuPr6ys33L4ZsxzBES+Kn5+fXr+G+B/Kf0D6x9e5fvB1rh98nesHX2fjfq1vNXKjw0XGREREZHKY4BAREZHJYYJTx2xtbTF9+nT5kfSHr3P94OtcP/g61w++zub1WpvlImMiIiIybRzBISIiIpPDBIeIiIhMDhMcIiIiMjlMcIiIiMjkMMG5A3PnzkVQUBDs7OwQHh6O6Ojom16/Zs0ahISEyOvbtm2LzZs311us5vI6L1iwAD169ICbm5s8evfufcv/L3RnP886q1atkp3ABw8erPcYzfF1zsnJwfjx4+Hj4yMrUVq0aMHfHXp4nb/88ku0bNkS9vb2svPu66+/juLi4nqL1xjt2rULDz/8sOwmLH4HrF+//paP2blzJzp27Ch/lps1a4alS5fqP1BRRUW1t2rVKo2NjY1m8eLFmpMnT2rGjBmjcXV11WRkZNR4/Z49ezRWVlaaWbNmaU6dOqV55513NNbW1prjx4/Xe+ym/Do//fTTmrlz52oOHz6sOX36tGbEiBEaFxcXTXJycr3Hbsqvs86FCxc0jRs31vTo0UMzaNCgeovXXF7nkpISTadOnTQDBgzQ7N69W77eO3fu1Bw5cqTeYzfl13nlypUaW1tb+VG8xr///rvGx8dH8/rrr9d77MZk8+bNmv/85z+atWvXiipszbp16256fXx8vMbBwUEzceJE+T741VdfyffFyMhIvcbJBOc2denSRTN+/Piq2xUVFRpfX1/NzJkza7z+qaee0gwcOLDafeHh4ZoXX3xR77Ga0+t8vfLyco2Tk5Nm2bJleozSPF9n8dp27dpVs3DhQs3zzz/PBEcPr/O8efM0wcHBmtLS0nqM0vxeZ3Ftz549q90n3oS7deum91hNBWqR4EyaNElzzz33VLtvyJAhmr59++o1Nk5R3YbS0lLExMTI6Y9r97USt6Oiomp8jLj/2uuFvn373vB6urPX+XpFRUUoKyuDu7u7HiM1z9d5xowZ8PT0xKhRo+opUvN7nTds2ICIiAg5ReXl5YU2bdrg448/RkVFRT1Gbvqvc9euXeVjdNNY8fHxchpwwIAB9Ra3OYhS6H3QLDfbvFNZWVnyF4z4hXMtcfvMmTM1PiY9Pb3G68X9VHev8/UmT54s54ev/0dFd/c67969G4sWLcKRI0fqKUrzfJ3FG+327dvxzDPPyDfcuLg4jBs3Tibtojss1c3r/PTTT8vHde/eXe5SXV5ejpdeeglvv/12PUVtHtJv8D4odhy/evWqXP+kDxzBIZPzySefyAWw69atkwsNqW7k5+fjueeekwu6PTw8lA7HpKnVajlK9t133yEsLAxDhgzBf/7zH3z77bdKh2ZSxMJXMTL2zTff4NChQ1i7di02bdqEDz74QOnQqA5wBOc2iF/qVlZWyMjIqHa/uO3t7V3jY8T9t3M93dnrrPP555/LBOePP/5Au3bt9Bypeb3O58+fR0JCgqyeuPaNWFCpVIiNjUXTpk3rIXLT/3kWlVPW1tbycTqtWrWSfwmLqRgbGxu9x20Or/O7774rk/bRo0fL26LKtbCwEGPHjpUJpZjiort3o/dBZ2dnvY3eCPy/dxvELxXx19S2bduq/YIXt8V8eU3E/ddeL2zduvWG19Odvc7CrFmz5F9ekZGR6NSpUz1Faz6vs2h1cPz4cTk9pTseeeQRPPjgg/JclNhS3fw8d+vWTU5L6RJI4ezZszLxYXJTd6+zWKt3fRKjSyq5TWPdUex9UK9LmE20DFGUFS5dulSWu40dO1aWIaanp8vPP/fcc5opU6ZUKxNXqVSazz//XJYvT58+nWXienidP/nkE1ke+vPPP2vS0tKqjvz8fAW/C9N7na/HKir9vM6JiYmyCnDChAma2NhYzcaNGzWenp6aDz/8UMHvwvReZ/H7WLzOP/74oyxl3rJli6Zp06ay+pVuTPxeFS05xCHSiNmzZ8vzixcvys+L11i81teXib/11lvyfVC09GCZuIESNfwBAQHyDVWUJe7bt6/qc/fff7/8pX+tn376SdOiRQt5vSiV27RpkwJRm/brHBgYKP+hXX+IX2BUtz/P12KCo7/Xee/evbKlhHjDFiXjH330kSzRp7p7ncvKyjTvvfeeTGrs7Ow0/v7+mnHjxmmuXLmiUPTGYceOHTX+vtW9tuKjeK2vf0xoaKj8/yJ+npcsWaL3OC3Ef/Q7RkRERERUv7gGh4iIiEwOExwiIiIyOUxwiIiIyOQwwSEiIiKTwwSHiIiITA4THCIiIjI5THCIiIjI5DDBISIiIpPDBIeIiIhMDhMcIiIiMjlMcIiIiMjkMMEhIiIimJr/BwFGm6pkLNVTAAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# set initial condition\n",
    "#\n",
    "if ictype == 1:\n",
    "    # smooth solution\n",
    "    rho0 = CF(cos(2*pi*x))\n",
    "\n",
    "    tend = 1\n",
    "else:\n",
    "    # nonsmooth solution\n",
    "\n",
    "    rhomin = 0.1\n",
    "    rhomax = 2\n",
    "    rho0 = IfPos((x - 0.75)*(x-0.25) , rhomin, rhomax) # 0.1 2\n",
    "\n",
    "    tend = 0.1\n",
    "\n",
    "\n",
    "\n",
    "gfrho_old.Set(rho0)\n",
    "\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "xvec = np.linspace(0,1, nel)\n",
    "aux = gfrho_old.vec.FV().NumPy()\n",
    "if order == 0:\n",
    "    plt.plot(xvec, aux)\n",
    "else:\n",
    "    plt.plot(xvec, aux[::order+1])\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "4916bcdf",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.01\n",
      "0.02\n",
      "0.03\n",
      "0.04\n",
      "0.05\n",
      "0.060000000000000005\n",
      "0.07\n",
      "0.08\n",
      "0.09\n",
      "0.09999999999999999\n",
      "0.10999999999999999\n",
      "0.11999999999999998\n",
      "0.12999999999999998\n",
      "0.13999999999999999\n",
      "0.15\n",
      "0.16\n",
      "0.17\n",
      "0.18000000000000002\n",
      "0.19000000000000003\n",
      "0.20000000000000004\n",
      "0.21000000000000005\n",
      "0.22000000000000006\n",
      "0.23000000000000007\n",
      "0.24000000000000007\n",
      "0.25000000000000006\n",
      "0.26000000000000006\n",
      "0.2700000000000001\n",
      "0.2800000000000001\n",
      "0.2900000000000001\n",
      "0.3000000000000001\n",
      "0.3100000000000001\n",
      "0.3200000000000001\n",
      "0.3300000000000001\n",
      "0.34000000000000014\n",
      "0.35000000000000014\n",
      "0.36000000000000015\n",
      "0.37000000000000016\n",
      "0.38000000000000017\n",
      "0.3900000000000002\n",
      "0.4000000000000002\n",
      "0.4100000000000002\n",
      "0.4200000000000002\n",
      "0.4300000000000002\n",
      "0.4400000000000002\n",
      "0.45000000000000023\n",
      "0.46000000000000024\n",
      "0.47000000000000025\n",
      "0.48000000000000026\n",
      "0.49000000000000027\n",
      "0.5000000000000002\n",
      "0.5100000000000002\n",
      "0.5200000000000002\n",
      "0.5300000000000002\n",
      "0.5400000000000003\n",
      "0.5500000000000003\n",
      "0.5600000000000003\n",
      "0.5700000000000003\n",
      "0.5800000000000003\n",
      "0.5900000000000003\n",
      "0.6000000000000003\n",
      "0.6100000000000003\n",
      "0.6200000000000003\n",
      "0.6300000000000003\n",
      "0.6400000000000003\n",
      "0.6500000000000004\n",
      "0.6600000000000004\n",
      "0.6700000000000004\n",
      "0.6800000000000004\n",
      "0.6900000000000004\n",
      "0.7000000000000004\n",
      "0.7100000000000004\n",
      "0.7200000000000004\n",
      "0.7300000000000004\n",
      "0.7400000000000004\n",
      "0.7500000000000004\n",
      "0.7600000000000005\n",
      "0.7700000000000005\n",
      "0.7800000000000005\n",
      "0.7900000000000005\n",
      "0.8000000000000005\n",
      "0.8100000000000005\n",
      "0.8200000000000005\n",
      "0.8300000000000005\n",
      "0.8400000000000005\n",
      "0.8500000000000005\n",
      "0.8600000000000005\n",
      "0.8700000000000006\n",
      "0.8800000000000006\n",
      "0.8900000000000006\n",
      "0.9000000000000006\n",
      "0.9100000000000006\n",
      "0.9200000000000006\n",
      "0.9300000000000006\n",
      "0.9400000000000006\n",
      "0.9500000000000006\n",
      "0.9600000000000006\n",
      "0.9700000000000006\n",
      "0.9800000000000006\n",
      "0.9900000000000007\n",
      "1.0000000000000007\n"
     ]
    }
   ],
   "source": [
    "# time step\n",
    "t = 0\n",
    "\n",
    "from ngsolve.webgui import Draw\n",
    "\n",
    "\n",
    "\n",
    "gfrho.vec.data = gfrho_old.vec\n",
    "\n",
    "while t< tend- 1e-8:\n",
    "    \n",
    "    if method_type == \"HDG_Newton\":\n",
    "        gfrho.vec.data = evolveIE_HDG_Newton(gfrho_old.vec)\n",
    "    elif method_type == \"DG\":\n",
    "        gfrho.vec.data = evolveIE_DG(gfrho_old.vec)\n",
    "    elif method_type == \"HDG\":\n",
    "        gfrho.vec.data = evolveIE_HDG(gfrho_old.vec)\n",
    "    else:\n",
    "        raise Exception(\"Method undefined\")\n",
    "\n",
    "    \n",
    "    gfrho_old.vec.data = gfrho.vec\n",
    "\n",
    "    t += dt\n",
    "\n",
    "    print(t)\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9d9d6df8",
   "metadata": {},
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "e95bb64f",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "err L2 = 0.6455159615766916\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x1ecea14e3c0>"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjgAAAGdCAYAAAAfTAk2AAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjEsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvc2/+5QAAAAlwSFlzAAAPYQAAD2EBqD+naQAAYJxJREFUeJzt3Qd4FNUaBuAvvQAJhBBCIHTpJZBAqIKAVBUUFZAuRZoNlaIiKiqKyEUUQeko0hSUovQOgUAg1FBCkZaEmt43e59zJpsCSUggm9md/d77zN3Z2dnNv4PJ/nvKf6z0er0eRERERBpirXYARERERIWNCQ4RERFpDhMcIiIi0hwmOERERKQ5THCIiIhIc5jgEBERkeYwwSEiIiLNYYJDREREmmMLC5SWloabN2+iRIkSsLKyUjscIiIiygdRmzgmJgZeXl6wts67jcYiExyR3Hh7e6sdBhERET2Ga9euoUKFCnmeY5EJjmi5MVwgFxcXtcMhIiKifIiOjpYNFIbP8bxYZIJj6JYSyQ0THCIiIvOSn+ElHGRMREREmsMEh4iIiDSHCQ4RERFpDhMcIiIi0hwmOERERKQ5THCIiIhIc5jgEBERkeYwwSEiIiLNYYJDREREmmPUBGfPnj14/vnn5aJYourgX3/99cjn7Nq1C40bN4aDgwOqV6+OxYsXP3TO7NmzUblyZTg6OsLf3x+BgYFGegdERERkjoya4MTFxaFhw4YyIcmPy5cvo1u3bnjmmWcQHByMd955B0OHDsXmzZszzlm5ciXGjh2LyZMn4+jRo/L1O3XqhFu3bhnxnRAREZE5sdKLtceL4gdZWWHt2rXo0aNHrueMHz8eGzduxKlTpzKO9e7dG5GRkdi0aZO8L1psmjRpgh9//FHeT0tLkwtvvfnmm5gwYUK+F+tydXVFVFQU16IiIiIyEwX5/DapxTYDAgLQoUOHbMdE64xoyRGSk5MRFBSEiRMnZjxubW0tnyOem5ukpCS5Zb1ARnHjKHB0CeBWDShdTbl1qwLYOhjn5xEREZkCXSoQdQ24exG4d1G5rdQCqJt7o4axmVSCEx4ejrJly2Y7Ju6LhCQhIQH379+HTqfL8ZyzZ8/m+rpTp07FZ599BqO7EQQEPTBmyMpaSXSqtweeehao1AqwczR+LERERMaSplO+1F/YDFzYCkScBtJSsp+TEs8Ex9hEi48Yt2MgEibRrVXoyjcGnv4gSwZ7CUiOAe5eULZDcwE7Z6BKG6D+y0DtFwBb+8KPg4iIqLCJES2XdgHBvwOh24CEe9kft3EA3Kqm92BUBSq3gppMKsHx9PREREREtmPivuhnc3Jygo2NjdxyOkc8NzdiRpbYjK68r7Jl/Y8h7jZwLRC4sEXJcmNuAuf/VbYS5YAmQwDfwUAxd+PHR0REVFDJ8cCJlcChn4HbIZnHHVyB6u2Apzop3VGu3mLcCEyFSSU4zZs3xz///JPt2NatW+Vxwd7eHr6+vti+fXvGYGUxyFjcHzNmDEyOlRVQ3AOo/ZyyiYQn4hQQsl7pyooJA3Z8Aez+FmjwKtBmPFDSCC1LREREBZUYDez/HjiyAEi4rxyzKwb4vAbUewmo0BSwMak0IhujRhYbG4vQ0NBs08DF9G83NzdUrFhRdh3duHEDS5culY+PGDFCzo4aN24cXn/9dezYsQOrVq2SM6sMRFfTwIED4efnh6ZNm2LmzJlyOvrgwYNh8kTC41lf2Vq/D5xeCxyaA9w8Bhz7FTj5B9DqHaDFW4C9s9rREhGRpY6vCV4GbP9c6YUQSlYC/N8AGvUDHF0BS58mLor2iZo2DxIJiijgN2jQIFy5ckWel/U57777Ls6cOYMKFSpg0qRJ8rysRBL07bffykHJPj4+mDVrlpw+nl8mNU1cXP6rB5X/kK4eUI65VACe/Qyo11NJioiIiIrCfwHApvFA2HHlfunqQPtPgFrPAdY2akdXoM/vIquDY0pMKsExEP8MokVn6yfKVDuhcmvgxbmAawW1oyMiIi1LilUSm2O/ZY6vaTseaDLMpCbDFOTz23RGA1k60VIj+jTHHAae+UiZbXVlLzCnJXDmb7WjIyIirboRBPzcOj25sQJ8BwFvBgHNR5tUclNQTHBMjZ0T0GYcMGIf4NUISIwEVg0A1r0JJMepHR0REWlprM3eGcCCjsC9S8rwiEEbgee/B4qXgbljgmOqRB2BIVuBVqJ+jxVwdCnw89NAxBm1IyMiInMXdwdY2h3Y/hmQlgrUfREYuQ+o3BJawQTHlNnYAR0mAwPXAyW8gLuhwMJOSoElIiKix3H7HDCvnTIMQkz77j4beHkR4FQKWsIExxxUaQ2M3K8s85AUDSx7FTi8QO2oiIjI3FzaBcx/Foj8DyhVBRi+S5n6rcEZu0xwzIWzG9B/LdCwD6DXARvHAps+VPpQiYiIHiVoMfBbTyApCvBuBgzdDpSpAa1igmNOxGj2HnOAdh8r9w/OBlb2B1IzV0onIiJ6qAyJqLW2/m1lvE39V4ABfwPFSkPLmOCYG9GMKBb07LlAWdjs3EZgRV8gJVHtyIiIyBSTmy0fA3u/U+63mQC8NA+wc4TWMcExV2I18r6rAVsnIHQrsKIPkJKgdlRERGRKyc2miUDAj8r9rtOBZyZqcrxNTpjgmLOqbZQkRxQFvLgDWN5bWfWViIgsm14P/DteWe9QeO5/QNNhsCRMcMydmGHV9w9lqp8YHb+8FwsCEhFZsrQ04J/3gcCflTpqz88C/F6HpWGCowWiMFP/NYB9ceDyHqXysS5F7aiIiEgN2yYDh+cryU33HwHfgbBETHC0omIzoN+a9DE525SlHSxvHVUiIssW8BNwYJay/8IPSo0bC8UER0sq+gOvLAasbIDjy5VpgUREZBlO/Qlsnqjst58MNO4PS8YER2tqdlYWShP2zQAO/aJ2REREZGxieMLaEcp+0+FAq3dh6ZjgaJHI2p/5SNn/dxxw5m+1IyIiImMJP6nUQ9MlA3W6A52/tpip4HlhgqNVohig72AxVxD4cxhwI0jtiIiIqLDFRCjrE4p1Ciu1BF78BbC2UTsqk8AER6tE9t7tO6BGZ0CXBKzoB8TeUjsqIiIqLKnJyqzZmJuAew2g9zKLqFCcX0xwtExk8aIkt/gPX/wCiF8E8QtBRETmb9N44NpBwMEV6L0ccCqldkQmhQmO1jm6KP/hO7gAVwOATRPUjoiIiJ7UkUXAkYVKrZue8wH36mpHZHKY4FgC8R+++AUQvwhHFgBBi9WOiIiIHtfVg8A/Hyj77ScBNTqqHZFJYoJjKWp0Atp9rOxvfB+4ekjtiIiIqKCibwIr+wNpKcqMqVZj1Y7IZDHBsSSt31N+IcQvxh+Dgfh7akdERET5laYD/hwKxN0CPOoC3X/idPA8MMGxJOIXQfxClK4ORN8A/h7D5RyIiMzFnm+B//Yr6w72+hVwKK52RCaNCY6lEb8QLy8EbOyBcxuBwHlqR0RERI9yZR+w+xtl/7mZQOlqakdk8pjgWKJyDYFnpyj7Wz4Cwk6oHREREeUm7q5SsFWfBvj0Axq8onZEZoEJjqXyfwOo0UUp7S3G4yTFqh0RERE9SAwj+HuUUsus9FNA12lqR2Q2mOBY8nicHj8BJbyAu6HKmlVERGRaDs0Fzm8CbByAVxYB9sXUjshsMMGxZM5uSn0cK2sgeBlweq3aERERkUHEGWDrJ8p+py8Bz/pqR2RWmOBYusotlenjwoaxXK+KiMgU6FKAtW8owwjEmoJNhqodkdlhgkPA0+OUbwYJ94D173DqOBGR2vZ+B4SfUNaXen4W692YaoIze/ZsVK5cGY6OjvD390dgYGCu57Zt2xZWVlYPbd26dcs4Z9CgQQ893rlz56J4K9pkaw/0mAtY2ylTx0+sVDsiIiLLdTNYqXkjdJ0OlCirdkRmyegJzsqVKzF27FhMnjwZR48eRcOGDdGpUyfcupVzV8iaNWsQFhaWsZ06dQo2NjZ45ZXs0+JEQpP1vOXLlxv7rWibZz2gbfpCnP+MA6JuqB0REZHlSU0C1o4A0lKVyvP1eqodkdkyeoIzY8YMDBs2DIMHD0adOnUwd+5cODs7Y+FCsQrqw9zc3ODp6Zmxbd26VZ7/YILj4OCQ7bxSpbhM/BNr+Q5Q3hdIigLWvcmuKiKiorZrKnA7BHB2B7rNYNeUqSY4ycnJCAoKQocOHTJ/oLW1vB8QEJCv11iwYAF69+6NYsWyT43btWsXPDw8ULNmTYwcORJ3797N9TWSkpIQHR2dbaMc2NgCPeYo0xEvbgeOLlE7IiIiy3H9CLD/e2X/+ZlAMXe1IzJrRk1w7ty5A51Oh7Jls/cfivvh4eGPfL4YqyO6qIYOHfpQ99TSpUuxfft2fPPNN9i9eze6dOkif1ZOpk6dCldX14zN29v7Cd+ZhpWpCbRPn5a4ZRIQHaZ2RERE2peanL4+YBpQ/1Wg9vNqR2T2THoWlWi9qV+/Ppo2bZrtuGjReeGFF+RjPXr0wIYNG3D48GHZqpOTiRMnIioqKmO7du1aEb0DM9VsZHpXVTQLABIRFYUD36d3TZUGuqSvOUWmm+C4u7vLAcIRERHZjov7YtxMXuLi4rBixQoMGTLkkT+natWq8meFhobm+LgYr+Pi4pJtozxY2wDPfw9Y2QAh64Bz/6odERGRdt29COxOnzXV+WulCCuZdoJjb28PX19f2ZVkkJaWJu83b948z+euXr1ajp3p16/fI3/O9evX5RiccuXKFUrcJGZV1QdajFH2N77PtaqIiIxBTObY8A6gSwKqtQPqcyFNs+miElPE582bhyVLliAkJEQOCBatM2JWlTBgwADZhZRT95TofipdunS247Gxsfjggw9w8OBBXLlyRSZL3bt3R/Xq1eX0cypEbSYAJSsB0deBnV+qHQ0RkfYcXwFc3gPYOgLdvuOsqUJkCyPr1asXbt++jU8++UQOLPbx8cGmTZsyBh5fvXpVzqzK6ty5c9i3bx+2bNny0OuJLq8TJ07IhCkyMhJeXl7o2LEjpkyZIruiqBDZOwPPzQB+66ks+Ca+WZRvrHZURETaEHcX2Pyhsi/qkLlVVTsiTbHS6y2v2ImYJi5mU4kBxxyPkw9/DgVOrla6rYbtUqaTExHRkxEF/Y4vB8rWA4aLv612akekqc9vk55FRSai01TAsSQQfhI4PE/taIiIzN+VfUpyAytlUgeTm0LHBIcerXgZoMOnyv7Or7jiOBHRk9ClKkviCH6DgQp+akekSUxwKH8aDwDK+Si1cbZ9pnY0RETm68gC4NZpZaXwdpPUjkazmOBQ/mvjiFVtheDflJLiRERUMLG3gR3ps1JFcsOaN0bDBIfyz7sJ4NNX2f/nfVHUSO2IiIjMy/bPlAWNPRsAvoPUjkbTmOBQwYixOA4uwM1jwLFf1Y6GiMh8XA/K/LspWsRFyzgZDRMcKpjiHkq9BsM3kYT7akdERGT6RIu3aPkWGvQGKvqrHZHmMcGhgms6HChTC4i/q8yqIiKivImxizePAvYlgGc5UaMoMMGhghP1GrpMU/YPLwBun1M7IiIi05UUA2yfouy3HQ+UyHuxaSocTHDo8VRtA9TsBuh1wBZOcyQiytW+mUDcLWUphqZvqB2NxWCCQ4/v2c8Ba1vgwmbg0i61oyEiMj1R14GAHzP/Ztraqx2RxWCCQ4/PvTrQZKiyv/ljIE2ndkRERKZFdE2lJgIVWwC1nlM7GovCBIeeTJvxgKMrEHESCP5d7WiIiEzHjaPAiRXKfqcvASsrtSOyKExw6MmIKpxPp6+psuMLIClW7YiIiNSn1wNbPlb2G/QCyjdWOyKLwwSHnlzTYUCpykBsOHDgB7WjISJS39mNwH/7AVtHoP0nakdjkZjg0JOzdQA6pNd12P89EH1T7YiIiNSTmgxsTZ9d2nwM4FpB7YgsEhMcKhx1ugPezYDUBGDXVLWjISJSz9ElwL1LQDEPoNU7akdjsZjgUOEQg+c6pheyOvYbcPu82hERERU9MQ5x9zeZRf0cSqgdkcVigkOFx7tpevG/NGDH52pHQ0RU9A7OAeJuA6WqAI0Hqh2NRWOCQ4Wr/STAyhoIWa+snEtEZCni7irjEIV2HyvL2pBqmOBQ4fKoDTTso+xvm6xMlSQisgT7ZgDJMYBnA6DuS2pHY/GY4FDhazsBsLEHruwFLu5QOxoiIuOLvAYE/qLsd5gMWPPjVW38F6DCV7Ii0GSYsr/tUyAtTe2IiIiMa9fXgC4ZqNwaqNZe7WiICQ4ZTev3APsSQPgJ4PQataMhIjKeW2eB4+lL1XT4lEsymAgmOGQcxUoDLd/KXMJBl6J2RERExrFjijJ7VCymWcFP7WgoHRMcMp5mowBnd+D+ZeD4crWjISIqfDePAWc3KLNH26VXLyaTwASHjMehONDqXWV/97dK+XIiIi3Z+ZVyW/8VwKOW2tFQFkxwyLiaDAGKewJRV4FjS9WOhoio8FwLBC5sAaxsgDbj1Y6GHsAEh4zLzkkZcCzs+Q5ISVQ7IiKiwrHzS+XWpw9Qupra0dADmOCQ8TUeALiUB2JuAkGL1Y6GiOjJXdkPXNoFWNsCT3+gdjSUAyY4ZHx2jsDT7yv7e78DkuPVjoiI6PGJCu2GsTeN+gOlKqsdEamV4MyePRuVK1eGo6Mj/P39ERgYmOu5ixcvhpWVVbZNPC8rvV6PTz75BOXKlYOTkxM6dOiACxcuFME7ocfm008pABh3Czg8X+1oiIge3+XdwH/7lIrthi9vZHkJzsqVKzF27FhMnjwZR48eRcOGDdGpUyfcunUr1+e4uLggLCwsY/vvv/+yPT5t2jTMmjULc+fOxaFDh1CsWDH5momJHN9hsmztMwfh7Z8JJMWoHRER0eO13uxIH3vjOxhwraB2RKRWgjNjxgwMGzYMgwcPRp06dWRS4uzsjIULF+b6HNFq4+npmbGVLVs2W+vNzJkz8fHHH6N79+5o0KABli5dips3b+Kvv/4y9tuhJ9GgN+BWFYi/m7lmCxGROQndDlwPBGwdgdZj1Y6G1EpwkpOTERQUJLuQMn6gtbW8HxAQkOvzYmNjUalSJXh7e8sk5vTp0xmPXb58GeHh4dle09XVVXZ95faaSUlJiI6OzraRCmxsgTYTlP0DPwJJsWpHRERUsNab3V8r+35DgBKeakdEaiU4d+7cgU6ny9YCI4j7IknJSc2aNWXrzt9//43ffvsNaWlpaNGiBa5fvy4fNzyvIK85depUmQQZNpE4kUrq9VRacRLuAUcWqB0NEVH+iVlT1w8rrTeGpWjIZJncLKrmzZtjwIAB8PHxQZs2bbBmzRqUKVMGP//882O/5sSJExEVFZWxXbt2rVBjpgK24rROH5S3fxZnVBGRGbXefKPs+w5i642lJzju7u6wsbFBREREtuPivhhbkx92dnZo1KgRQkND5X3D8wrymg4ODnLgctaNVNTgVaBkJSD+DhC0SO1oiIge7co+4GqAMnOq5dtqR0NqJzj29vbw9fXF9u3bM46JLidxX7TU5Ifo4jp58qScEi5UqVJFJjJZX1OMqRGzqfL7mqQyG7vM6sb7vwdSEtSOiIgob4bWG1m41EvtaMgUuqjEFPF58+ZhyZIlCAkJwciRIxEXFydnVQmiO0p0IRl8/vnn2LJlCy5duiSnlffr109OEx86dGjGDKt33nkHX3zxBdatWyeTH/EaXl5e6NGjh7HfDhWWhn0AV28gNgI4yjWqiMiE/XcAuLIXsLbLXECYTJ6tsX9Ar169cPv2bVmYTwwCFmNrNm3alDFI+OrVq3JmlcH9+/fltHJxbqlSpWQL0IEDB+QUc4Nx48bJJGn48OGIjIxEq1at5Gs+WBCQTJioiyP+UGwcC+z7H9B4oFLxmIjI1Oyeptw26se6N/kQm5SKKevPYFznmihd3AFqsdKLwjIWRnRpidlUYsAxx+OoKDUJ+N5HWaOq23dAE6WVjojIpFYMX/CssubUm0eBUpXUjsikJaboMHjRYQRcugu/SqWwekRz2fOixue3yc2iIgti65DZ3Lv3f0BqstoRERHl3HojutWZ3OQpRZeGMb8flclNcQdbTHquTqEmNwXFBIfUJQbsFS8LRF8HTq5SOxoiokw3g4HQrYCVNasWP4IuTY/3Vh3HtpBbcLC1xvyBfmjoXRJqYoJD6hLjbpqPUfb3zgDSdGpHRESk2DdDua33slKglHIkRrp8/NcprDt+E7bWVpjbzxfNqpaG2pjgkPr8BgOOJYF7F4Ezf6sdDRERcPs8cGadss+ZU3n6etNZLA+8CtEb9b9ePnimlgdMARMcUp9DCcB/RGYrjuWNeyciU7N/pmibAGp2A8pmzuKl7Bbuu4yfd1+S+1NfrI/nG5pOjSAmOGQa/N8A7IoBESeBC1vVjoaILFnkVeDESmWfY29y9e/JMEzZeEbuiynhvZtWhClhgkOmwdkNaPK6sr93OltxiEg9B34A0lKBKm2ACn5qR2OSjly5h3dWBss/1f2aVcTINtVgapjgkOkQg43FOi/XDimVQ4mIilrsrczq6oYlZSibi7djMXTpESSlpqFDbQ98+nxdVaeD54YJDpkOsTqvqBQq7P1O7WiIyBId/AlITQTK+wFVnlY7GpNzOyYJgxYFIjI+BQ0ruGJWn0awtTHNVMI0oyLL1eItwMoGuLgduHlM7WiIyJIkRAKB8zNbb0ywVULtKsXDlh7BtXsJqOjmjAWDmsDZ3ugrPj02JjhkWtyqAPVfzpxRRURUVA7PA5JjAI86QI3OakdjcrVuxv1xAsHXIuHqZIdFg5vAXcV1pvKDCQ6ZnpbvKLch64E7oWpHQ0SWICUBODg3829QlkWgCfhhR2hGIb85fRujWpniMHX8FyTTI2pOPNVJqUER8IPa0RCRJQj+HYi/A7h6A/VeUjsak7LhxE3M2Hpe7n/evR5aVHeHOWCCQ6apVXorTvByICZC7WiISMvEEjFianjGbE47tSMyGcevRco1poTXW1bBa/6mVesmL0xwyDRVbA5UaArokoBDc9SOhoi0TCwRc/8y4OQGNO6vdjQmIzwqUQ4qFtPBn6lZBh91qw1zwgSHTJOYvWBoxTm8EEiMVjsiItIiUalOLssAoOlwwL6Y2hGZzIypN349glsxSahRtricDm5jbV6zypjgkOmq0QVwrwEkRQFBi9SOhoi06NIuIOw4YOukJDgEMWNq4pqTOH49CiWd7TB/QBOUcDS/bjsmOGS6xCwGURdHODgHSE1SOyIi0pr93yu3omuqWGm1ozEJ8/dextpjN2SLzU+vNUbF0s4wR0xwyLQ1eBUoUQ6ICQNOrFI7GiLSkpvBwKWdSnFRMbiYsPv8bUz9N0TuT+pW22xmTOWECQ6ZNlsHoNmozG9aaWlqR0REWmu9EdPCS1WCpbt0OxZjfj+KND3Qy88bA1tUhjljgkOmz3cQ4OAK3L0AnP9X7WiISAvuXwHO/KXst3wbli4mMUXOmIpJTIVvpVL4vIdpLqBZEExwyPQ5ugB+g5X9Az+qHQ0RaYEY16dPA6q1Azzrw5Klpenx7srjuHg7Dp4ujpjTrzEcbG1g7pjgkHnwfwOwtgOuHgCuH1E7GiIyZ/H3gKNLlf0Wb8LSzdx+AdtCImBva42f+/vCo4QjtIAJDpkHFy+g/ivK/oFZakdDRObsyEIgJR4oWw+o+gws2aZT4Zi1/YLcn/pifTT0LgmtYIJD5qPFmMxFOO9dUjsaIjJHotzEoZ8zW2/MfJzJkzgfEYP3VgXL/cEtK6OnbwVoCRMcMh9l6wLVOyj95qL/nIiooES5ibhbQAkvoF5PWKqo+BQMX3oEcck6NK9aGh92Na9lGPKDCQ6ZF0N/+bHflH50IqL8EmUmDItqNhtpsYtq6tL0GLP8KK7cjUf5kk6Y3bcx7Gy0lw5o7x2RtlVpo8x4EP3nRxaoHQ0RmZPQbcCdc4B9CcB3ICzVtM1nsffCHTjaWeOXAb5wK2YPLWKCQ+ZF9Jcblm849AuQkqh2RERkLgwTFERy4+gKS/R38A38vFsZw/jtyw1R10u714EJDpmfui8CLuWVfvSTXL6BiPLh5jHgyl7A2lbpnrJAp25EYfyfJ+T+yLbV8HxDL2gZExwyP6Lf3PAHShT+4/INRPQohiKhdV8CXLU1Wyg/7sQm4Y1fg5CYkoa2Ncvg/Y41oXVFkuDMnj0blStXhqOjI/z9/REYGJjrufPmzUPr1q1RqlQpuXXo0OGh8wcNGiRLSGfdOnfuXATvhExG44FKP7roT7+4Xe1oiMiURV0HTq9V9puPhqVJTk3DqGVHcSMyAVXci+H73o3kSuFaZ/QEZ+XKlRg7diwmT56Mo0ePomHDhujUqRNu3bqV4/m7du1Cnz59sHPnTgQEBMDb2xsdO3bEjRs3sp0nEpqwsLCMbfny5cZ+K2Rqyzc0HqDsB3D5BiLKg6h7o9cBlVsDXj6wJHq9HpPXnUbg5Xso7mCLeQN84epkGbPHjJ7gzJgxA8OGDcPgwYNRp04dzJ07F87Ozli4cGGO5y9btgyjRo2Cj48PatWqhfnz5yMtLQ3bt2f/lu7g4ABPT8+MTbT2kAUu32BlDVzaBUScVjsaIjJFSbFA0BKLbb359eB/WB54Vc7P+KFPI1T3KAFLYdQEJzk5GUFBQbKbKeMHWlvL+6J1Jj/i4+ORkpICNze3h1p6PDw8ULNmTYwcORJ3797N9TWSkpIQHR2dbSMNKFUJqP2Csh/wk9rREJEpCl4GJEUBbtWApzrBkuwPvYPP1p+R+xM618IztTxgSYya4Ny5cwc6nQ5ly5bNdlzcDw8Pz9drjB8/Hl5eXtmSJNE9tXTpUtmq880332D37t3o0qWL/Fk5mTp1KlxdXTM20e1FGtE8ffkGMZsqJkLtaIjIlKTpgIPpX36ajxLfsGEprtyJk+NuRFG/lxqVx/Cnq8LSmPS/9tdff40VK1Zg7dq1coCyQe/evfHCCy+gfv366NGjBzZs2IDDhw/LVp2cTJw4EVFRURnbtWvXivBdkFF5NwEqNAV0ycDh+WpHQ0Sm5Nw/wP0rgFMpoOFrsBQxiSkYuvQIohJS4ONdEl+9VF9OxrE0Rk1w3N3dYWNjg4iI7N+sxX0xbiYv06dPlwnOli1b0KBBgzzPrVq1qvxZoaGhOT4uxuu4uLhk20hDDP3qIsFJSVA7GiIyFQGzlVu/IYC9MyxBqi4No38/htBbsfB0ccQv/X3haGcDS2TUBMfe3h6+vr7ZBggbBgw3b9481+dNmzYNU6ZMwaZNm+Dn5/fIn3P9+nU5BqdcuXKFFjuZkVrPASUrAgn3gOMr1I6GiEzB9SDgagBgbQc0HQZLmTElxtzsOX8bTnY2mDfADx4umb0flsboXVRiiriobbNkyRKEhITIAcFxcXFyVpUwYMAA2YVkIMbUTJo0Sc6yErVzxFgdscXGxsrHxe0HH3yAgwcP4sqVKzJZ6t69O6pXry6nn5MFsrEF/NML/4n+dhb+I6KD6a039V8BSuTdY6AViw9ckbOmRG/U/3r5oH4F7S7DYBIJTq9evWR30yeffCKnfgcHB8uWGcPA46tXr8o6NgZz5syRs69efvll2SJj2MRrCKLL68SJE3IMTo0aNTBkyBDZSrR3717ZFUUWqlE/wMEFuHNeWVCPiCy8sN9fmYOLLcCOsxGYsiFzxlTnepaR1OXFSi/atCyMmCYuZlOJAcccj6Mhmz9Siv5VbQsM+FvtaIhILVsmKQtrVnkaGLgeWhcSFo2X5xxAXLIOvfy88XVP7Q4qLsjnt0nPoiIqkKbDsxT+U77JEJEFFvY7ml7Yr5n2C/uFRSXg9cWHZXLTvGppTOlRT7PJTUExwSFtFf4TA46FQ3PUjoaI1HB8OZAoCvtVBZ7qCC0T08AHLTyMsKhEVCtTDHP7+cLelh/rBrwSpM0p48dXAnF31I6GiIqSmGBwMP3LjZh4oOHCfkmpOrzx6xGci4iBRwkHLHm9KVydLWONqfzS7r8+WSZvf8CrEaBLAo4sUjsaIipKoVuBexcBB1fAR7uF/dLS9Hhv1XEcvKQsoLlocBNUKGUZdX4KggkOaYvoe26WPmvi8DwgNUntiIioqAv7+Q4AHIpDq776JwQbToTB1tpKdkvV9bLs6eC5YYJD2lOnB1CiHBAbAZxeq3Y0RFQUIk4Dl3crEw3EhAON+mXPRczfd1nuf/tKA7R6yl3tkEwWExzSHlt7oMnQzG90llcJgcjyGBbVrP2CUtlcg34/dBVf/XNW7k/oUgsvNqqgdkgmjQkOaZPf64CtIxB+AvjvgNrREJExxd4GTqxW9g1d1Brzd/ANfPTXSbk/sm01jGhTTe2QTB4THNImZzegYe/s3+yISJuOLFQmFpT3BbybQmu2nonA2FXHZWN0/2aVMK5TTbVDMgtMcEi7DN/kzm4E7il91kSkMWIiweH5mVPDNVbkbn/oHYz+/Sh0aXq81Kg8PnuhLgv55RMTHNKuMjWBau3EGrtA4Dy1oyEiYxATCeJuKRML6vaAlgRcvIuhS44gOTUNneqWxbSXG8DamslNfjHBIctoxTn2K5AUo3Y0RFSYRJ+NobCfmFhgY6eplpvBiwORkKLD0zXKYFafRrC14Ud2QfBqkbZVaw+Urg4kRQPBy9WOhogK07VDQFiwMqHAdzC0Ys/523J9qcSUNLStWQa/9PeFg62N2mGZHSY4pG2iVLv/CGX/0FyllDsRaYNhAkH9V4BipaEFu87dwtClR5CUmob2tTzwc39fONoxuXkcTHBI+xr2UUq3ixLuopQ7EZm/yGtAyAZlv9lIaMG2MxEYvjRIjrl5tk5Z/NSvMVtungATHNI+UbK9cX9l39BfT0TmTSzFotcBVZ4GytaFuVt1+Bre+C0IyTplQPHs15jcPCkmOGQZROl2UcL90k7gVoja0RDRk0iOA4KWZE4NN2N6vR6zd4Zi3J8n5FTwno0r4MfXGsPelh/PT4pXkCxDqUpArW6ZY3GIyHwdXwEkRgKlqgA1OsGcVwX/bP0ZfLv5nLwvqhNPf6UB7DhbqlDwKpLlMHzTO74SiL+ndjRE9LhTww/9rOz7vwFYm2c3TmKKDm+tOIbFB67I+5OeqyPXl2IRv8LDBIcsR6UWgGd9IDUBOJrevE1E5uXiDuDOOcC+BODTF+YoIjoRvX45iA0nwmBnY4Xve/tgSKsqaoelOUxwyHKIb0aGVpzA+YAuVe2IiKigDF3MjfoCji4wN8evReKFH/fJW1cnOywe3BTdfcqrHZYmMcEhy1KvJ+DsDkRfB86mTzElIvNw9yJwYYv4tqJMHDDDFcFf/TkAEdFJeMqjONaNaYmW1d3VDkuzmOCQZbFzBPzSK55ysDGReTGMvREDi0tXg7kQdW2++icEb68Izijgt2ZUC1QqXUzt0DSNCQ5ZHr8hgLUtcDUAuBmsdjRElB+JYrmVZZmDi83E1bvxeOXnAPyy55K8P7JtNfwywA8lHLWzbpapYoJDlsdFrDr8YvZvhERk2kRykxwLlKkFVH0G5mDDiZvoNmtvxngbsezC+M61YMMVwYsEExyyTIb1qU79AcTeUjsaIspLmi771HATn0odl5SKiWtOYszvxxCTlAq/SqXwz9ut0amup9qhWRQmOGSZKvgB5f0AXTIQtFjtaIgoLxe2AvcvA46uQINeMPWVwDv+bw+WB16VediYZ6pjxfBmKF/SSe3QLA4THLJchlacw/OB1GS1oyGi3BxKX0Ou8UDA3jQH5kbFp+D91ccxYGEgbkQmoEIpJywb4o/3O9WELSsTq4JXnSxXne5AcU8gNgI485fa0RBRTsTacZd2KWvJNR0GU1xLauOJMHT43278EXRdttoMalEZm995Gi04BVxVtur+eCIV2doDTYYCO79QVhlv8KraERHRgwxjb8RaciUrwpScvB6FKRvOIPCKsvRL1TLFMK1nA/hVdlM7NGKCQxbPdxCw51vg5lHg2mHAu4naERGRgVgzTiysaWKrht+KTpQLZP5x9LpcGsvRzhrDn66GUW2rwdHOPNfG0qIi6aKaPXs2KleuDEdHR/j7+yMwMDDP81evXo1atWrJ8+vXr49//vnnoSbBTz75BOXKlYOTkxM6dOiACxcuGPldkCYVLwPUf1nZZ+E/ItNy7Fdl7TixhpxYS05ld2OT8PW/Z9F2+i6sDlKSmx4+XtjxXluMfbYGkxtLS3BWrlyJsWPHYvLkyTh69CgaNmyITp064datnKfmHjhwAH369MGQIUNw7Ngx9OjRQ26nTp3KOGfatGmYNWsW5s6di0OHDqFYsWLyNRMTE439dkiLDCXfxTic6JtqR0NEglgrLnBe5oQAFaeG345JkpWIW32zE3N3X0R8sg4+3iVlNeKZvRvBizOkTJKVXjSHGJFosWnSpAl+/PFHeT8tLQ3e3t548803MWHChIfO79WrF+Li4rBhQ+Y6Qc2aNYOPj49MaES4Xl5eeO+99/D+++/Lx6OiolC2bFksXrwYvXv3fmRM0dHRcHV1lc9zcTG/xdrICBZ2ViobP/0B0O5jtaMhojPrgFX9AefSwLtnlGVWitjF27FYeuAKVh65hsSUNHmsQQVXvNXuKbSv7QErE6/Ho0UF+fw2agtOcnIygoKCZBdSxg+0tpb3AwICcnyOOJ71fEG0zhjOv3z5MsLDw7OdI96sSKRye82kpCR5UbJuxnA/Lhn95h/C4fQBZ2SGU8aPLAJS2BJIZDKDi30HF2lyo0vTY+uZCPRfcAjtv9uNJQH/yeRGtNgsGtwEf49uiQ51yjK5sfRBxnfu3IFOp5OtK1mJ+2fPns3xOSJ5yel8cdzwuOFYbuc8aOrUqfjss89gbD/uDMW+0Dtye65BOUzsWrtIijuJVq0UnR5JqTp5m6JLk4u7idvUNL3cF7+0Yj9Vp+zr9Hp5myZvIW/F66Tpxeul38/y+lkZfrHF/1tbWUFUHbdKv5X3rZVbW2truS9uRRkIG2tr2FpbyTLldjbK43a21rCT99P3baxgJ59XxH88aj0HuFRQVhk/9SfQqG/R/nwiyhR2Avhvn7JmXJMhRdZa8/exG/jz6A1Zx0YQf4ba1y4rp323qFaaSY2ZsYhZVBMnTpTjgAxEC47oJitsYhE10Te74vBVbDgRJr8FvNGmGka0qQp7G2vEJqUiJjFV3opS3sqtTu7HJafK58YnK8cSknVISFG2xJTM+2Il2sQHbkUCozUiEbK3tVY2G+XWQW42cLDL3BezFwy3YoCfYXOSmzWc7LPct7eBc/p9Z3tbuS+P2dnA1sYWaDoU2PapMtjY5zWTLwdPpFmBP2fWqnLxMupsKPG3+q/gGzhxPSrjeElnO/RuUhF9/SvC283ZaD+fzDjBcXd3h42NDSIiIrIdF/c9PXNek0Mcz+t8w604JmZRZT1HjNPJiYODg9yMzb24A6a+VB/9mlXEZ+vPIPDyPczafgGzd4bK1pKiYmghES0holVE3re2ktU0DS0oWTdDK4zYl/+TLTJKK4yBYdfQmCNuRRuPPkuLj2wNkveV1qGMLf1+qi79Nk20LimtSSnprUoPXh7Z2iQTPl2RXDORMHnZV8Im2MMh/AQ+/P4XXC3ug2IONihmb4tiDrZwdrBB8fT94g7pt45i3wbFHezkuSXSb1m5lOgxxd0BTqzO3nVcSMTfptM3o7Hj7C1sD4nA8SxJjfj79/RT7ujRqLxcM4ozosyfURMce3t7+Pr6Yvv27XImlGGQsbg/ZsyYHJ/TvHlz+fg777yTcWzr1q3yuFClShWZ5IhzDAmNaJERs6lGjjSNOgl1vVyxcngz/HsqHF9uDMlo7hRES0PxLB+QxeQHpg2c5X5my4LSumCb3vJgLW8dROvEA60WSsuGTUYLh2jtKPLunUIgEh/ZtSaSnlTDrR7JusxWqiTDlktLVlJGi1daZutXlpYwQyuYSJoS0lvLDImVeP7lVAessW2JPrY70fLuH/g9vMJjvx/x7yWSnxIOtighkiC5b6fcZhy3y3wsfd9FbmLfTv77skmcLI5YG06XBHg1Aio8eV2qa/ficfDSXRy8dA8HLt5BWFT2MXZibM2LjcqjW4Ny8ksqaYfRu6hE19DAgQPh5+eHpk2bYubMmXKW1ODBg+XjAwYMQPny5eU4GeHtt99GmzZt8N1336Fbt25YsWIFjhw5gl9++UU+Lv7gi+Tniy++wFNPPSUTnkmTJsmZVYYkyhSIOLvWL4eOdcoiPDoxoxVAJCL0MKU1Sek+Kiri25xIbAxdgyLp0YWXBtbuRBfbIPzc2QN3bD2ULsQk3QNdi5nHsm6G7kJDUiWmlz4u0drm4qQkPkryY9i3y9g3PC7uuzil32acZ8uWJDIvuhTg8ILMwn4FTPDFl50zYdE4cS1Sdjkdunwv2xdMw5ePVk+5o30tD7Sr5QEPl6KfnUUaSXDEtO/bt2/LwnxiELBoddm0aVPGIOGrV6/KmVUGLVq0wO+//46PP/4YH374oUxi/vrrL9SrVy/jnHHjxskkafjw4YiMjESrVq3ka4rCgKZGfMBUKMU+XFMkklDDmB23YvbKwbLNgGOtYX1lLzrFrwee/bxArykGesvEJzEVMUkpypir9HFXMSIJEscTleOZt6mIznJMnCtalkQ33b24ZLk9LtEqKBMiJ0NilD0pevCxB88RHwZsRaIiE7IOiLkJFPMA6vbI88uJaIk5FxGD8+ExOB8Ri5CwaJyPiJG/Nw9+URBTu5tVLQ1/sVVxY/eThTB6HRxTxDo4lKezG4EVrwGOJYGxIYB90Sao4lcyLlmXkQBFJ6TfJqYg+oH7mY8rjxmeU1hjl0TLmqE1SHSxZbYgpR9zzOx2E92uImlSxiVl3ueYJMq3+c8C1wOBNhOQ0HIcbkYlICwyEdfvx+O/e/H4724crtxRbsXvSE5KF7OXCU2DCiXRuFIp+FUqJVvPyfI+v/mvTvSgGp2BkpWAyP+Ak6uU9aqKkGgxMSQI5Vwf7zXEeCaR6ESlJz+PTpSyJ01iMwwSj4xPkRuQvam/IJSxZ0qCJBKebOPQDPvp49EMxwxj0uQYNcOtna2c+cauXvMjEnfRbXs/PkXWDBMtk/fjk3EnNll25TreCsY71wORAlt03lMVFzdvyvP1RMtMFfdiqOFZAjU8SqCmZ3HUr1ASXq6ObHUkiQkO0YOsbZTlG7Z8BBycCzQeaHZTxkVdIdHtltH19pgfRtEJqQ+1DmXtXhPdaVm71wzjkpTHUjKqv4rbxJQk3IlNKqT3ZyW7z7JN988YmK8Myn+wbICjoWxARqmBzAH74r4YoC8H7tsoCZRhEz9LPGYpH5ri31508zw4oF/+G6YqA/fFrTJ2TYf4pFTEi8H7WcekJSqlLwyJtNwSU+TsydzMsPsVsAHW6ZrhYkpxeUwkuWIZBLFVLu2MSqWLoZK8dUZFt2JMdClPTHCIctKoH7DzK+B2CHB5N1C1LSyJ+DBXkgdbeLo+/tg20ZIkBmQbkiFxm3Wgdmy2fUMdKOWYGPQdl/4BKu6LD1PD+AqloKXyAVpURJIjyy/IpMc6s/RCekkGcZtRisFKKcEgjolyC4ZyDErBS1GKIbM8g7iVqZMozwBD4czsP1spy5BZnkE5aCjPkKVMgyzeKWarZpZnMJRoEOUZlNmKmaUaZDFQnT5j9mKSuNWlZZSDMNZ1LFXMDqWc7eUmkvAqjjHofuqgfE+VuozFpqpNUM7VSXaFWkpiSYWPCQ5RTpxKAj59gMPzlZLxFpbgFBaRCJR0tpdbYRCz1JTEJzXbdP94QxmA9FIASouD0sLwYMuD0pqUtWVCJz/gM0oRpChlCh6sXSXPEcM+iqg2k6kQCYkoQeEoi2Rmbx17sLSFGH8luhdFGQR562gHV6csM/yc7ORzHkpadk4F9KmAtz/8WrRX662SxjDBIcqNKDImEpxz/wL3LgFuVdWOyOIZuo1cne2KrDaToQ6TbPWQtZmUVg5lUwpVipYlcd+wJErmrTiWvUXFsCSKYamUbEUzRa2wHJpPxCFZgNNQiDO9+KbSCqQcU5ZJMbQeKfdtHizwKZZQSe9yy9r6lK07zlAxvKjqaqUmAUcWGKWwH1k2JjhEuXF/CqjeAQjdBgTOAzortZrIMqhRm8kinVoDxN0GSngBtZ9XOxrSEI7QIsqLKDYmHPsNSIpROxoibRFNU2LtN0GsBWdj/JY5shxMcIjyUq0dULo6kBQNBC9XOxoibbl2CAgLBmwdgcZFW46BtI8JDlFeRJVtw7gA8U0zTXsrtxOp5uAc5bb+K0Cx0mpHQxrDBIfoURr2ARxcgXsXlfE4RPTkoq4DIeuVfQ4uJiNggkP0KA7Fgcb9lf1D6d84iejJiIH7eh1QuTXgmbnWIFFhYYJDlB9NhwFW1sDFHcCts2pHQ2TekuOBoMXKfrP0gfxEhYwJDlF+lKoM1Oyq7BtmfRDR4zmxAkiMVH6vxNpvREbABIcovwzfNI+vAOLvqR0NkRlPDf9Z2W/6hrL2G5ERMMEhyq9KLYGy9YHUBODoUrWjITJPl3YCt88C9sWBRn3VjoY0jAkOUX6JevjNRmQOkNQV3UKPRJqbGu7TF3B0VTsa0jAmOEQFUe9lwNkdiL4OnE2f4kpE+XMnFLiwRVk63f8NtaMhjWOCQ1QQdo6A3+vK/kEONiYqkMD0sTc1OgGlq6kdDWkcExyigmoyBLC2A64dBG4cVTsaIvOQEAkcW6bss7AfFQEmOEQFVcITqPdS9vEERJQ3sWBtShxQpjZQta3a0ZAFYIJD9DgM30BPrwWiw9SOhsi0iQH5hqnhYqC+GLBPZGRMcIgeR/nGQMXmQFoKcHi+2tEQmbZzG4Goq4CTG9Cgl9rRkIVggkP0pIX/jiwEUhLUjobIdBm6csUAfTsntaMhC8EEh+hx1ewGuFYEEu4BJ1apHQ2RaRID8a8GKAPzmwxVOxqyIExwiB6XjW1mLQ/xDVWUoCeinFtvxMB8l3JqR0MWhAkO0ZNo3F8pOX87RClBT0SZxAD802uUfa4aTkWMCQ7RkxCl5kXJeSHgJ7WjITIth+cBaalAxRaAVyO1oyELwwSH6EnJbiorIHQrcPu82tEQmYbkeODIImWfrTekAiY4RE9KlJyv2UXZP8TlG4ikk6uUAfglKwK1uqkdDVkgJjhEhaHZKOU2+Hcg/p7a0RCpSwy4NwwuFkUxrW3UjogsEBMcosJQuRVQtj6QmgAELVY7GiJ1XdwO3D4L2JcAGvVXOxqyUEZNcO7du4e+ffvCxcUFJUuWxJAhQxAbG5vn+W+++SZq1qwJJycnVKxYEW+99RaioqKynWdlZfXQtmLFCmO+FaK8idLzzUcr+4G/AKnJakdEpJ6A2cpt4wGAo4va0ZCFMmqCI5Kb06dPY+vWrdiwYQP27NmD4cOH53r+zZs35TZ9+nScOnUKixcvxqZNm2Ri9KBFixYhLCwsY+vRo4cx3wrRo9XrCRT3BGLE1Ni1akdDpI6IM8DFHYCVdWadKCIVWOn1xqlOFhISgjp16uDw4cPw8/OTx0Sy0rVrV1y/fh1eXl75ep3Vq1ejX79+iIuLg62trRK0lRXWrl372ElNdHQ0XF1dZcuQaF0iKjR7pgM7pgCeDYA39nBRQbI8f49WVg6v0wN4dYna0ZDGFOTz22gtOAEBAbJbypDcCB06dIC1tTUOHTqU79cxvAlDcmMwevRouLu7o2nTpli4cCHyytOSkpLkRcm6ERmFWGvH1gkIPwFc2ad2NERFK/ZW5rIlzceoHQ1ZOKMlOOHh4fDw8Mh2TCQpbm5u8rH8uHPnDqZMmfJQt9bnn3+OVatWya6vnj17YtSoUfjhhx9yfZ2pU6fKjM+weXt7P+a7InoEZzfA57Xs4xCILMXh+YAuGajQFPBuonY0ZOEKnOBMmDAhx0G+WbezZ88+cWCilaVbt26ym+vTTz/N9tikSZPQsmVLNGrUCOPHj8e4cePw7bff5vpaEydOlC1Bhu3atWtPHB/RI6eMn/8XuBOqdjRERSMlQUlwBMOAeyIVZe/3yYf33nsPgwYNyvOcqlWrwtPTE7du3cp2PDU1Vc6UEo/lJSYmBp07d0aJEiXkWBs7O7s8z/f395ctPaIrysHB4aHHxbGcjhMZhXt1oEYXJcE5+BPw3Ay1IyIyvuMrgPi76YX9nlM7GqKCJzhlypSR26M0b94ckZGRCAoKgq+vrzy2Y8cOpKWlyYQkr5abTp06yYRk3bp1cHR0fOTPCg4ORqlSpZjEkOkQ32BFgiMK/7X7WOm6ItKqtDQlmRf8RwI2Bf5oITKfMTi1a9eWrTDDhg1DYGAg9u/fjzFjxqB3794ZM6hu3LiBWrVqyccNyU3Hjh3ljKkFCxbI+2K8jth0Op08Z/369Zg/f76cRh4aGoo5c+bgq6++kvVziEyq8F+5hkrhvyML1I6GyLhCtwF3zgMOLkCjfmpHQyQZNc1etmyZTGrat28vZ0+JAcGzZs3KeDwlJQXnzp1DfHy8vH/06NGMGVbVq1fP9lqXL19G5cqVZXfV7Nmz8e6778qZU+K8GTNmyESKyLQK/40B1gwDAucBLd4CbNnCSBoVkD7Jg4X9yBLq4Jgy1sGhIqFLAWY2AGJuAt1n85stadPNYOCXNoCVDfD2caAkZ6mSxuvgEFk8Gzug2Uhl/8APygKERFoT8KNyW+8lJjdkUpjgEBmT70BlwUGx8OCFrWpHQ1S4Iq8Bp9Yo+yzsRyaGCQ6RMTm6KkmOcCBz/BmRJhyaC+h1QJWnAS8ftaMhyoYJDpGxiW4qa1vgyl7g5jG1oyEqHAmRQNBiZb/F22pHQ/QQJjhExuZaQVlpXDiQPl6ByNwdXQIkxwIedYDq7dWOhughTHCIioJhfMLptUDkVbWjIXoyqcnAwbmZ/22LsghEJoYJDlFRKNcAqNpWGa9wcI7a0RA9mVN/KuUPinsC9V9WOxqiHDHBISoqLdKrbQctARLuqx0N0eMR5Q5E2QPB/w0WsCSTxQSHqKhUaw941AVS4oAji9SOhujxhG4Hbp0G7IoBfoPVjoYoV0xwiIqKGKfQ8i1lX3RTpSSqHRFRwe2fqdz6DgKcSqkdDVGumOAQFSUxm8qlAhB3Czi+XO1oiArmRpBS7kCUPWg+Su1oiPLEBIeoqJdvaD46s/Bfmk7tiIjyb1966039V5TyB0QmjAkOUVGTKy6XBO5dAs5uUDsaovy5EwqErFf2W7KwH5k+JjhERc2hONB0eOY3Yi7CSeYgQMyc0gM1OgMetdWOhuiRmOAQqUEkOLaOwM2jypgGIlMWEwEEp48ZY+sNmQkmOERqKF4GaNRP2d//vdrREOXt0BxAlwRUaApUbK52NET5wgSHSC2yxL01ELoNCD+pdjREOUuMBg4vVPZbvcNlGchsMMEhUotbFaDui8o+W3HIVAUtApKiAPcaQI0uakdDlG9McIjUZBjPcGoNcO+y2tEQZSeKUQb8pOy3eAuw5kcGmQ/+10qkpnINgeodlEU4RV0cIlNy/HcgNhxwKQ806KV2NEQFwgSHSG2t31Nuj/0GxISrHQ2RQpeaWdhPtN7Y2qsdEVGBMMEhUlulFsrMFF1y5irNRGo7vQaI/A9wLq0UpyQyM0xwiEypFUesMh5/T+1oyNKlpQF7Zyj7zUYB9s5qR0RUYExwiEyBGIfjWR9IiQMO/ax2NGTpzv8L3A4BHFyAJkPVjobosTDBITIForaIoRXn0FwgKUbtiMhSiaVD9n6n7Ivkxqmk2hERPRYmOESmovYLQOnqQGKk0lVFpIbLu4EbQcpSIqJ7ishMMcEhMhXWNkCrd5X9gB+VGiRERc3QetN4oLKkCJGZYoJDZErqvwq4VABiI4Bjv6odDVmaa4eBy3sAa1ugxZtqR0P0RJjgEJkSUWvEUN1Y1CBJTVY7IrIku79Rbhv2Bkp6qx0N0RNhgkNkakTNkeKeQPR1pZIsUVEQ425CtwJWNpkD3onMGBMcIlNj55jZiiPGQ+hS1I6ILMHub5XbBq8CblXVjobItBOce/fuoW/fvnBxcUHJkiUxZMgQxMbG5vmctm3bwsrKKts2YsSIbOdcvXoV3bp1g7OzMzw8PPDBBx8gNTXVmG+FqGj5DgKKlQEirwInVqodDWld2HGl9o2VNVtvSDOMmuCI5Ob06dPYunUrNmzYgD179mD48OGPfN6wYcMQFhaWsU2bNi3jMZ1OJ5Ob5ORkHDhwAEuWLMHixYvxySefGPOtEBUtUTlWrP8j7JmurAtEZCy70//G1usJuD+ldjREpp3ghISEYNOmTZg/fz78/f3RqlUr/PDDD1ixYgVu3ryZ53NFy4ynp2fGJlqADLZs2YIzZ87gt99+g4+PD7p06YIpU6Zg9uzZMukh0owmQ5R1gO5fBk79oXY0pFXhp4CzG0S1SeDpD9SOhsj0E5yAgADZLeXn55dxrEOHDrC2tsahQ4fyfO6yZcvg7u6OevXqYeLEiYiPj8/2uvXr10fZsmUzjnXq1AnR0dGytYhIM+yLZU7V3fMtkKZTOyLSIvHfllD3RaBMTbWjISo0tjCS8PBwOT4m2w+ztYWbm5t8LDevvfYaKlWqBC8vL5w4cQLjx4/HuXPnsGbNmozXzZrcCIb7ub1uUlKS3AxEMkRkFkSp/P3fA3dDgdNrgfovqx0RacmtEODM38o+W2/I0ltwJkyY8NAg4Ae3s2fPPnZAYoyOaJERrTRiDM/SpUuxdu1aXLx48bFfc+rUqXB1dc3YvL1Z34HMhEMJoPnozHESbMWhQm+90SvLhJSto3Y0ROomOO+9954cX5PXVrVqVTl25tatW9meK2Y6iZlV4rH8EuN3hNDQUHkrnhsREZHtHMP93F5XdHNFRUVlbNeuXSvo2yZST9PhgKMrcOec0opDVBgizgCnlJZxtt6QFhW4i6pMmTJye5TmzZsjMjISQUFB8PX1lcd27NiBtLS0jKQlP4KDg+VtuXLlMl73yy+/lMmToQtMzNISA5Hr1Mn5G4iDg4PciMySSG7EWJwdXwC7pgJ1egA2RutdJksh/lsytN6Ua6B2NETmM8i4du3a6Ny5s5zyHRgYiP3792PMmDHo3bu3HF8j3LhxA7Vq1ZKPC6IbSsyIEknRlStXsG7dOgwYMABPP/00GjRQfgE7duwoE5n+/fvj+PHj2Lx5Mz7++GOMHj2aSQxpl/8IwMlNGYtzcpXa0ZC5CzsBhKxTZk61nah2NETmVwdHzIYSCUz79u3RtWtXOVX8l19+yXg8JSVFDiA2zJKyt7fHtm3bZBIjnie6w3r27In169dnPMfGxkbW1BG3ojWnX79+Mgn6/PPPjflWiNQfi2OobizWC2J1Y3oSO79Sbuu9xLE3pFlWer1eDwsjZlGJwcZiPE7WGjtEJi05Dvi+IRB3G3h+FuA7UO2IyFzXnJrXTqlaPDqQhf1Is5/fXIuKyJzq4rQamzn7JTWz9AFRgVtvGvRickOaxgSHyJz4vQ6UKAdEXQOOLlU7GjI3Vw8BoduUFcPbjFM7GiKjYoJDZG4rjRsWQxQrjackqB0RmZOdXyi3jfpyxXDSPCY4ROam8QDApQIQEwYcWah2NGQuLu8FLu8BrO1Y94YsAhMcInNj65DZvSBacRK59Ag9gphLsu1TZV8MTi9ZUe2IiIyOCQ6ROfLpC5R+Coi/CwT8qHY0ZOrObgRuHAHsnIGnOfaGLAMTHCJzJCoZt5+k7B/4EYi9rXZEZKp0qcD29DphzUYBJbIvVkykVUxwiMyVKLHv1RhIiUtfNJEoBydWKOuYOZUCWr6ldjRERYYJDpG5srICOqSPqxCDje9fUTsiMjUpicBOseYUlNl3Yl0zIgvBBIfInFVtA1R9BkhLySzgRmRweD4QfR1wKQ80GaZ2NERFigkOkbnrMFm5PbEKCD+ldjRkKhKjgL3TlX2xoKaooURkQZjgEJk7r0ZA3RfFXODMwaREB34AEu4D7jWAhn3UjoaoyDHBIdKCdpOU8vsXNisF3ciyRYcBAbMz/9sQs+6ILAwTHCItKF0N8Bus7G/+EEhLUzsiUtOOL4CUeMDbH6j9vNrREKmCCQ6RVohxFg4uQPgJZWowWaawE0DwMmW/45fKbDsiC8QEh0grirlnLsS5fQqQHK92RKTGkgxbPlLGY9XrCXg3UTsiItUwwSHSEv8RgGtFIOYml3CwROfFGKw9gI0D0D59dh2RhWKCQ6QlYirws+nF//bNBGLC1Y6IioouBdiavnxHs5FAqUpqR0SkKiY4RFpT9yWgQhNlCQcx2JQsQ9Bi4M55wLk00Hqs2tEQqY4JDpHWiEGlYnCpcOw3Fv+zlKJ+u6ZmDjbnkgxETHCINKmif2bxv00TlMGnpF27pwHxdwH3moBverkAIgvHBIdIqzp8Btg6Alf2Amf+UjsaMpbb54BDc5X9Tl+xqB9ROiY4RFolBpm2fEfZ3/wxkByndkRU2ETL3L/jgLRUoGZX4KkOakdEZDKY4BBpWat3gJIVlRWl985QOxoqbCHrgEu7lGnhovWGiDIwwSHSMjunzA++A7OAuxfVjogKiyjkuFkU9QPQ8m3ArYraERGZFCY4RFpX6zmg6jOALllZp4q0Yd//gKhrgKs30OpdtaMhMjlMcIgsYdp4l2mAtS1wfpNS7ZbM273LwP7vlf1OXwL2zmpHRGRymOAQWYIyNZTqtsK/44GURLUjoichWuJ0SUDVtkDtF9SOhsgkMcEhshRPjwOKewL3LwP7OODYbIVsAM79o7TIiZY5rhZOlCMmOESWwtEF6Jxe7VbMqBL1U8i8JEYD/3yg7Ld4CyhTU+2IiEwWExwiSyKqGz/VCUhLAda/A6SlqR0RFYRYW0ysFF+qCtBmnNrREFlugnPv3j307dsXLi4uKFmyJIYMGYLY2Nhcz79y5QqsrKxy3FavXp1xXk6Pr1ixwphvhUgbRHdGt+mAnTNw9QBw7Fe1I6L8un4ECPxF2X/uf0oJACJSJ8ERyc3p06exdetWbNiwAXv27MHw4cNzPd/b2xthYWHZts8++wzFixdHly5dsp27aNGibOf16NHDmG+FSDtE4b9n0uunbJ0ExESoHRE9ik60uL2trC3WoDdQ7Rm1IyIyeUZbtCQkJASbNm3C4cOH4efnJ4/98MMP6Nq1K6ZPnw4vL6+HnmNjYwNPT89sx9auXYtXX31VJjlZiRahB88lonzyHwGcXAWEHQc2TwReXqh2RJSXgB+BiFOAk5syLZyI1GvBCQgIkEmIIbkROnToAGtraxw6dChfrxEUFITg4GDZtfWg0aNHw93dHU2bNsXChQuh52rJRPknFmR8fhZgZQ2c+hO4sFXtiCivmje7vlH2O34BFHNXOyIiy27BCQ8Ph4eHR/YfZmsLNzc3+Vh+LFiwALVr10aLFi2yHf/888/Rrl07ODs7Y8uWLRg1apQc2/PWW2/l+DpJSUlyM4iOjn6s90SkKV4+QLNRSuuAGHA86gDg6Kp2VJSVGAS+7k0gNQGo3BrweU3tiIi024IzYcKEXAcCG7azZ88+cWAJCQn4/fffc2y9mTRpElq2bIlGjRph/PjxGDduHL799ttcX2vq1KlwdXXN2MRYHyIC8MyHQKnKymKcXMbB9ByeB1zZqwwKf/571rwhMmaC895778nxNXltVatWleNjbt26le25qampcmZVfsbO/PHHH4iPj8eAAQMeea6/vz+uX7+erZUmq4kTJyIqKipju3btWgHeMZGG2RcDeswR06uAY79xGQdTIhZG3TpZ2X/2c6B0NbUjItJ2F1WZMmXk9ijNmzdHZGSkHEfj6+srj+3YsQNpaWkyIclP99QLL7yQr58lxumUKlUKDg4OOT4ujuf2GJHFq9QCaD5a6apa9xYwKgBwdlM7KsuWpgP+Gql0TVV5GvB7uCWbiFQaZCzGznTu3BnDhg1DYGAg9u/fjzFjxqB3794ZM6hu3LiBWrVqycezCg0NlVPKhw4d+tDrrl+/HvPnz8epU6fkeXPmzMFXX32FN99801hvhUj72n0MuNcAYsOBf1lATnUBs4FrhwD7EkD32YA1a7ISFZRRf2uWLVsmE5j27dvL6eGtWrXCL7+kF6oCkJKSgnPnzsmuqKzErKgKFSqgY8eOD72mnZ0dZs+eLVuIfHx88PPPP2PGjBmYPDm9KZeICk4UjRNdVWJW1cnVwJl1akdkuW6dVSoWC2JKuKhbREQFZqW3wPnVYhaVGGwsxuOIKstElG7bZ8pCnM7uwKiDQPFHdxFTIRf0W9ARuHkUqP4s0Hc1BxYTPebnN9s9iShT2wmAR10g/g7w1wiuVVXUdn6pJDdiuv4Lok4Rkxuix8UEh4gy2ToAPecDto5A6DZl4DEVjdDtwL7/KfuiCKPLw9XeiSj/mOAQUXZl6wCdv1b2t38GXA9SOyLtE+uBrX1D2fd7HajLtfWInhQTHCJ6mO8goE4PIC0V+GMwkBildkTaJboBRXITdxvwqAN0+krtiIg0gQkOET1MjP0QlXPFDJ7I/5SlHCxvPkLROPA9cGknYOsEvLxImdFGRE+MCQ4R5cypJNBzIWBtC5xeAxxdqnZE2nPtMLB9irLfdRrgUUvtiIg0gwkOEeXOuwnQbpKy/88HwI2jakekHbG3gFUDAL0OqNcTaNRf7YiINIUJDhHlrcVbQI0ugC4JWNlP+WCmJ5OarCQ3MTeVCtLPzeSUcKJCxgSHiPImlgl46Weg9FNA9A1g1UClIB09vk0TgKsBgIML0Pt3wJEFR4kKGxMcIno0UXiuz3LlA/nqAeUDmh5P0GLgyAJlBXdRc8j9KbUjItIkJjhElD/ig1h8IIsP5sPzgaAlakdkfq4eAja+n7nAaY1OakdEpFlMcIgo/8QHcruPlP1/3gf+C1A7IvMRdR1Y1R9ISwHqdAdav6d2RESaxgSHiAqm9ftA7RcAXTKwvLey+jXlLf4e8OtLQGyEUsyv+08cVExkZExwiKhgxAfziz8DFZoCiZHAbz2B6JtqR2W6UhKA5X2AO+eAEl7KCuEOxdWOikjzmOAQUcHZOwOvrUyfWXUd+O1lICFS7ahMT5oO+HMocO0g4OAK9PsTcK2gdlREFoEJDhE9Hmc35QO7eFng1mlgRV8gJVHtqEyHWNpCFEc8uwGwsQf6/K4sZEpERYIJDhE9vlKVgL5/APYlgP/2AWuGAbpUtaMyDbunZU4Hf2keULmV2hERWRQmOET0ZMo1AHovA6ztgJB1wJ9DWAhw1zfArvRVwbt8A9TtoXZERBaHCQ4RPbmqbYBevypdMWf+Av4YrCxHYIndUju/ykxuOnwK+L+hdlREFslKrxe/kZYlOjoarq6uiIqKgotL7iXSdTodUlIs/JtoIbO3t4e1KP1P2nR+i7JelVi3qmY34JXFgK09LIL4U7rjC2DvdOX+s1OAlm+pHRWRRX5+C0xwcrhA4pKEh4cjMpKzQgqbSG6qVKkiEx3SqNBtwPLXlCRHLNL56hLA1gGaJv6MbvsU2D9Tud/pK6D5aLWjItIcJjhPeIHCwsJkcuPh4QFnZ2dYsSBXoUhLS8PNmzdhZ2eHihUr8rpq2cUdSu2X1ESgUiul+0rMutIi0RW3cSxw7Fflfpdp7JYiMoEEx9ZYQZgr0S1lSG5Kly6tdjiaU6ZMGZnkpKamykSHNKpaO+C1VcrUcTG7an4HpcBd6WrQlIT7wKoBwOU9gJU10HU60GSI2lEREQcZP8ww5ka03FDhM3RNiUSSLGDg8ZAtgKs3cO8iML898N8BaMa9y8CCjkpyY18c6LOCyQ2RCWGCkwt2nxgHr6uFEYXthm4HvBorrR1LuwPHV8DsXT2oJGx3zgMu5YHXN3FlcCITwwSHiIyrRFlg0MbMBTrXvgH8NRpIioVZLr2w51tgUVcg/i5QzgcYtgPwrK92ZET0ACY4RFQ0a1e9sgRoM16p7Bv8G/Bza+BGEMxG5DVgyfPKVHC9Dqj3MjD4H6CEp9qREVEOmOBQoXY//fXXX2qHQaZK1D965kNg0AalW+feJWUMy94ZSsuIKTu9FpjbEvhvvzLepsdcoOd8wL6Y2pERUS6Y4BBR0RJrMo3cD9TpAaSlAts/U8azXD0Ek3P3olLTZ/UgIDEKKO8LjNgL+PQRGb3a0RFRHpjgaKzOzNSpU2UhPScnJzRs2BB//PGHLFzYoUMHdOrUSe4L9+7dQ4UKFfDJJ59kzGoaMmRIxnNr1qyJ77///qGfsXDhQtStWxcODg4oV64cxowZI49XrlxZ3r744ouyJcdwnyhHTqWUKsfdZysLdd48BizsCPw5FIi6oXZ0QGI0sPUT4KdmwLmNgJUN0Pp94PXNgFtVtaMjonxgHZz8EElBSnzR/1w75wJ9SxTJzW+//Ya5c+fiqaeewp49e9CvXz9Ze2bJkiWoX78+Zs2ahbfffhsjRoxA+fLlMxIckRyJhGf16tWy/s+BAwcwfPhwmcS8+uqr8pw5c+Zg7Nix+Prrr9GlSxdZaGn//v3yscOHD8vaQYsWLULnzp1hY2NjpItCmiH+227UD3iqI7BjCnD0V+DkaiBkg7LEQdM3gGJFXIsqJQE4vhzYORWIu6Ucq9ZeqUzsUatoYyGiJ2K0SsZffvklNm7ciODgYFn7JD/LHohQJk+ejHnz5snzW7ZsKT9UxYe1gWh5ePPNN7F+/XpZ9r9nz56ypaF48eKFUgkxMTERly9fli0Zjo6OysHkOOArLxS5D2/mu48/KSkJbm5u2LZtG5o3b55xfOjQoYiPj8fvv/8uk5cBAwbgnXfewQ8//IBjx45lu7YPEq0zYskK0QokiIRo8ODB+OKLL3I8X7TcrF27Fj165L5yco7Xl0i4GQxsmghcTa+VY+sI1H8FaDYSKFvXuD9btBodng8ELQYS7inHSldXEhuRgLE7isgkmEQl4+TkZLzyyivyw3bBggX5es60adNkC4NobRAfgJMmTZLdKmfOnMn4MOzbt69cSmHr1q2yKJ/4wBUtDeID3JKFhobKRObZZ5996N+hUaNGcl/8e4gERLTAPJg4CrNnz5ZdUFevXkVCQoJ8ro+Pj3zs1q1bsgJx+/bti/BdkUXx8lFmJYnVyPfNBMKCleUPxFa5NeDzGlD9WaB4mcL5eUkxwKXdwKk/gTN/KzOjhJIVgWajAb/XLWehUCINMlqC89lnn8nbxYsX5+t80Xozc+ZMfPzxx+jevbs8tnTpUpQtW1bOzOnduzdCQkKwadMm2R3i5+cnzxEtEV27dsX06dPh5eVlvK4i0ZpS1MTPzafYWKWmiGg1Ey0tWYnxMoJIgIKCgmT30YULF7Kds2LFCrz//vv47rvvZFJaokQJfPvttzh0SBn4KcblEBmdaCmp+6IyAPnaIeDgHCBkPXBlr7KJKeblGyutKlWfAcrUBJxK5u+1k+OBu6HK61zYAlzZD6QplcslsWZWsxFAza6ANbtYicydyYzBEd0WojtEDIY1EM1Q/v7+CAgIkAmOuC1ZsmRGciOI80VXlfggFgNcjfZH18Sng9apU0cmMqL1pU2bNjme895778lr9e+//8qksFu3bmjXrp18TIyladGiBUaNGpVx/sWLFzP2RcIjBg5v374dzzzzTI6vL9aW4hIMVGi/cxWbKZuoP3N0KXB+ExB+QqmdI7ZdU5Vznd2VNa7cqgGOrtlfR4ydE9PRxRadw+DlUpWBpzopY4HKNSia90ZElpXgiORGEC02WYn7hsfErRjImpWtra0ce2I4J7fxKWLL2oenNSIBES0w7777rhww3KpVq4xBwKKf0t3dXXY/iSSxcePG+OCDDzBw4ECcOHECpUqVkt1VosVs8+bNsnvw119/lS1lYt/g008/lYOTxb+BGGQcExMjX1+MiRIMCZAYOyWSLfG6RE+spDfQ7iNliw4DQrcC5zcD1w8DsRFA/B1lEy0+j+JYEijXUFlWQbQCiXE2HF9DpEkFSnAmTJiAb775Js9zRDdSrVqmNdtAzC4ydJlp2ZQpU+SMKfF+L126JFu7RDIzceJE9OrVSyYo4r4grseWLVtkwrJy5Uq88cYbctCxOE8MFu7Tp49szRGtPQYiIRKDhP/3v//JZEokTS+//HLG46J7S8yyEoPERTfZlStXVLkOpGEu5YDGA5TNMI5GtM6IejViQU8xCyorazvArYrSuiNaeZzdVAmbiEx8FtXt27dx9+7dPM+pWrVqxorRhjE4YtbOo2ZRiQ/katWqyQ9Zw8BWQXS3iPtippRogRDdLPfv3894PDU1VQ5AFjOEcuuiyqkFx9vbO/+zqKjQ8PoSEZHJzaISrQNiMwbxgefp6Sm7OAwJjngjYmzNyJEj5X0x+FUkSmKgrK+vrzy2Y8cO2SUjxurkRnSXGAbaEhERkfYZrZKxGOwqauCIWzHwVOyLzTDbRxBdWWLasiC6RURLj6ixsm7dOpw8eVLWbBEzowx1VWrXri2LyA0bNgyBgYFy/Ieo1SIGIBttBhURERGZHaMNMhYVckU9GwNDLZadO3eibdu2cv/cuXOymclg3LhxiIuLk3VtREuNGCgrpoVn7cpYtmyZTGpEPRZDoT9RO4eIiIjI6JWMTVmBKxlToeH1JSKiohiDw8U2iYiISHOY4ORCDFymwmeBDYZERGTJhf5MhZjiLsb2iHWXxIwxcV8MgKbCSW5EqQFxPUXVYyIiImNhgvMAkdyI8SFiQU+R5FDhEslNhQoV5HpYRERExsIEJwei1aZixYqyiCDXVipcouWGyQ0RERkbE5xcGLpR2JVCRERkfjjImIiIiDSHCQ4RERFpDhMcIiIi0hxbS67FIioiEhERkXkwfG7np6aaRSY4MTEx8tbb21vtUIiIiOgxPsfFkg15sci1qESVYlHjpkSJEoVexE9klyJxunbt2iPXyaDHx+tcNHidiwavc9HgdTb/ay1SFpHceHl5ybp1ebHIFhxxUUSxOWMS/6D8BTI+XueiwetcNHidiwavs3lf60e13BhwkDERERFpDhMcIiIi0hwmOIXMwcEBkydPlrdkPLzORYPXuWjwOhcNXmfLutYWOciYiIiItI0tOERERKQ5THCIiIhIc5jgEBERkeYwwSEiIiLNYYLzGGbPno3KlSvD0dER/v7+CAwMzPP81atXo1atWvL8+vXr459//imyWC3lOs+bNw+tW7dGqVKl5NahQ4dH/rvQ4/33bLBixQpZCbxHjx5Gj9ESr3NkZCRGjx6NcuXKyZkoNWrU4N8OI1znmTNnombNmnBycpKVd999910kJiYWWbzmaM+ePXj++edlNWHxN+Cvv/565HN27dqFxo0by/+Wq1evjsWLFxs/UDGLivJvxYoVent7e/3ChQv1p0+f1g8bNkxfsmRJfURERI7n79+/X29jY6OfNm2a/syZM/qPP/5Yb2dnpz958mSRx67l6/zaa6/pZ8+erT927Jg+JCREP2jQIL2rq6v++vXrRR67lq+zweXLl/Xly5fXt27dWt+9e/cii9dSrnNSUpLez89P37VrV/2+ffvk9d61a5c+ODi4yGPX8nVetmyZ3sHBQd6Ka7x582Z9uXLl9O+++26Rx25O/vnnH/1HH32kX7NmjZiFrV+7dm2e51+6dEnv7OysHzt2rPwc/OGHH+Tn4qZNm4waJxOcAmratKl+9OjRGfd1Op3ey8tLP3Xq1BzPf/XVV/XdunXLdszf31//xhtvGD1WS7rOD0pNTdWXKFFCv2TJEiNGaZnXWVzbFi1a6OfPn68fOHAgExwjXOc5c+boq1atqk9OTi7CKC3vOotz27Vrl+2Y+BBu2bKl0WPVCuQjwRk3bpy+bt262Y716tVL36lTJ6PGxi6qAkhOTkZQUJDs/si6rpW4HxAQkONzxPGs5wudOnXK9Xx6vOv8oPj4eKSkpMDNzc2IkVrmdf7888/h4eGBIUOGFFGklned161bh+bNm8suqrJly6JevXr46quvoNPpijBy7V/nFi1ayOcYurEuXbokuwG7du1aZHFbggCVPgctcrHNx3Xnzh35B0b8wclK3D979myOzwkPD8/xfHGcCu86P2j8+PGyf/jBXyp6suu8b98+LFiwAMHBwUUUpWVeZ/FBu2PHDvTt21d+4IaGhmLUqFEyaRfVYalwrvNrr70mn9eqVSu5SnVqaipGjBiBDz/8sIiitgzhuXwOihXHExIS5PgnY2ALDmnO119/LQfArl27Vg40pMIRExOD/v37ywHd7u7uaoejaWlpabKV7JdffoGvry969eqFjz76CHPnzlU7NE0RA19Fy9hPP/2Eo0ePYs2aNdi4cSOmTJmidmhUCNiCUwDij7qNjQ0iIiKyHRf3PT09c3yOOF6Q8+nxrrPB9OnTZYKzbds2NGjQwMiRWtZ1vnjxIq5cuSJnT2T9IBZsbW1x7tw5VKtWrQgi1/5/z2LmlJ2dnXyeQe3ateU3YdEVY29vb/S4LeE6T5o0SSbtQ4cOlffFLNe4uDgMHz5cJpSii4ueXG6fgy4uLkZrvRH4r1cA4o+K+Da1ffv2bH/gxX3RX54TcTzr+cLWrVtzPZ8e7zoL06ZNk9+8Nm3aBD8/vyKK1nKusyh1cPLkSdk9ZdheeOEFPPPMM3JfTLGlwvnvuWXLlrJbypBACufPn5eJD5ObwrvOYqzeg0mMIankMo2FR7XPQaMOYdboNEQxrXDx4sVyutvw4cPlNMTw8HD5eP/+/fUTJkzINk3c1tZWP336dDl9efLkyZwmboTr/PXXX8vpoX/88Yc+LCwsY4uJiVHxXWjvOj+Is6iMc52vXr0qZwGOGTNGf+7cOf2GDRv0Hh4e+i+++ELFd6G96yz+HovrvHz5cjmVecuWLfpq1arJ2a+UO/F3VZTkEJtII2bMmCH3//vvP/m4uMbiWj84TfyDDz6Qn4OipAeniZsoMYe/YsWK8gNVTEs8ePBgxmNt2rSRf/SzWrVqlb5GjRryfDFVbuPGjSpEre3rXKlSJfmL9uAm/oBR4f73nBUTHONd5wMHDsiSEuIDW0wZ//LLL+UUfSq865ySkqL/9NNPZVLj6Oio9/b21o8aNUp///59laI3Dzt37szx763h2opbca0ffI6Pj4/8dxH/PS9atMjocVqJ/zNuGxERERFR0eIYHCIiItIcJjhERESkOUxwiIiISHOY4BAREZHmMMEhIiIizWGCQ0RERJrDBIeIiIg0hwkOERERaQ4THCIiItIcJjhERESkOUxwiIiISHOY4BARERG05v8QX7n9xoXhZQAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "\n",
    "xvec = np.linspace(0,1, nel)\n",
    "aux = gfrho_old.vec.FV().NumPy()\n",
    "if order == 0:\n",
    "    plt.plot(xvec, aux)\n",
    "else:\n",
    "    rhoplot = aux[::order+1]\n",
    "    plt.plot(xvec, rhoplot )\n",
    "    \n",
    "\n",
    "\n",
    "fesaux = L2(mesh, order = 0)\n",
    "gfrho_aux = GridFunction(fesaux)\n",
    "\n",
    "if ictype == 1:\n",
    "    gfrho_aux.Set(rho0)\n",
    "    aux = gfrho_aux.vec.FV().NumPy()\n",
    "    plt.plot(xvec, aux, label = \"exact\")\n",
    "\n",
    "    errL2 = sqrt(Integrate((gfrho - rho0)**2, mesh))\n",
    "    print(f\"err L2 = {errL2}\")\n",
    "else:\n",
    "    filename = \"exactRP.dat\"   # change to your filename\n",
    "    \n",
    "    # Skip the first two header lines, use comma as delimiter\n",
    "    data = np.loadtxt(filename, delimiter=\",\", skiprows=2)\n",
    "    \n",
    "    # Extract columns\n",
    "    x_dat   = data[:, 0]\n",
    "    rho_dat = data[:, 1]\n",
    "    \n",
    "    # ---- Overlay new plot ----\n",
    "    plt.plot(x_dat, rho_dat, '--', label=\"File data\")\n",
    "\n",
    "plt.legend()\n",
    "\n",
    "\n",
    "\n"
   ]
  }
 ],
 "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.13.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}