{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Rates of Convergence\n",
    "\n",
    "Copyright (C) 2010-2020 Luke Olson<br>\n",
    "Copyright (C) 2020 Andreas Kloeckner\n",
    "\n",
    "<details>\n",
    "<summary>MIT License</summary>\n",
    "Permission is hereby granted, free of charge, to any person obtaining a copy\n",
    "of this software and associated documentation files (the \"Software\"), to deal\n",
    "in the Software without restriction, including without limitation the rights\n",
    "to use, copy, modify, merge, publish, distribute, sublicense, and/or sell\n",
    "copies of the Software, and to permit persons to whom the Software is\n",
    "furnished to do so, subject to the following conditions:\n",
    "\n",
    "The above copyright notice and this permission notice shall be included in\n",
    "all copies or substantial portions of the Software.\n",
    "\n",
    "THE SOFTWARE IS PROVIDED \"AS IS\", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR\n",
    "IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,\n",
    "FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE\n",
    "AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER\n",
    "LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,\n",
    "OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN\n",
    "THE SOFTWARE.\n",
    "</details>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "from firedrake import *\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## A Boundary Value Problem\n",
    "\n",
    "Consider\n",
    "$$u_{*} = \\sin(\\omega \\pi x) \\sin(\\omega \\pi y)$$\n",
    "on the unit square with $\\omega = 2$ as a start."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "omega = 2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Refine the Mesh and Check the Error"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Now solving 4x4...\n",
      "Now solving 8x8...\n",
      "Now solving 16x16...\n",
      "Now solving 32x32...\n",
      "Now solving 64x64...\n",
      "Now solving 128x128...\n"
     ]
    }
   ],
   "source": [
    "errsH0 = []\n",
    "errsH1 = []\n",
    "hs = []\n",
    "\n",
    "for nx in [4, 8, 16, 32, 64, 128]: # , 256, 512]:\n",
    "    print(f\"Now solving {nx}x{nx}...\")\n",
    "    mesh = UnitSquareMesh(nx, nx)\n",
    "\n",
    "    V = FunctionSpace(mesh, \"Lagrange\", 1)\n",
    "    Vexact = FunctionSpace(mesh, \"Lagrange\", 7)\n",
    "\n",
    "    x = SpatialCoordinate(V.mesh())\n",
    "    u_exact = interpolate(sin(omega*pi*x[0])*sin(omega*pi*x[1]), Vexact)\n",
    "    f = 2*pi**2*omega**2*u_exact\n",
    "\n",
    "    # all four sides of the square\n",
    "    bc = DirichletBC(V, 0.0, [1, 2, 3, 4])\n",
    "    u = TrialFunction(V)\n",
    "    v = TestFunction(V)\n",
    "\n",
    "    a = inner(grad(u), grad(v))*dx\n",
    "    L = f*v*dx\n",
    "\n",
    "    u = Function(V)\n",
    "    solve(a == L, u, bc)\n",
    "\n",
    "    EH0 = errornorm(u_exact, u, norm_type='L2')\n",
    "    EH1 = errornorm(u_exact, u, norm_type='H1')\n",
    "    errsH0.append(EH0)\n",
    "    errsH1.append(EH1)\n",
    "    hs.append(1/nx)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[1.63571727 1.8993806  1.97405822 1.9934556  1.99840339]\n",
      "[0.8332833  0.95536396 0.98862547 0.99714187 0.99928454]\n"
     ]
    }
   ],
   "source": [
    "errsH0 = np.array(errsH0)\n",
    "errsH1 = np.array(errsH1)\n",
    "hs = np.array(hs)\n",
    "rH0 = np.log(errsH0[1:] / errsH0[0:-1]) / np.log(hs[1:] / hs[0:-1])\n",
    "rH1 = np.log(errsH1[1:] / errsH1[0:-1]) / np.log(hs[1:] / hs[0:-1])\n",
    "print(rH0)\n",
    "print(rH1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x7fcc02d964d0>"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXwAAAD8CAYAAAB0IB+mAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+j8jraAAAgAElEQVR4nO3deXzU1b3H/9fJvpB9DyFhh5AEBNlslapsKkRBqZKorV20tte63NZW27r8em21t/f+eu1te3v9tV5bMAFEQSO4VK21rSA7WdgNZA8hyWSyJ5OZ8/vjm0x2yD4zyef5ePDQ+c53vnMSMh9Ozvt8z1Faa4QQQox/bo5ugBBCiLEhBV8IISYIKfhCCDFBSMEXQogJQgq+EEJMEFLwhRBigvBwdAP6opRKBVIDAgLunz17tqObI4QQLuXw4cOVWuuInseVM8/DX7x4sT506JCjmyGEEC5FKXVYa72453EZ0hFCiAlCCr4QQkwQTlnwlVKpSqmXzGazo5sihBDjhlMWfK11ltb6gaCgIEc3RQghxg2nLPhCCDEhZe+AXyXDs8HGf7N3jOjlnXJaphBCTDjZOyDrYbA0GY/NRcZjgPl3jshbSA9fCCGcwYc/7Sz2HSxNxvER4pQ9/I4br2bOnOnopgghxMjTGurKoOIEXDwBFSeNHn1fzMUj9rZOWfC11llA1uLFi+93dFuEEGJYmkxQcQoq8ozCfvGEUeibazrPmRQNHt7Q1tL79UFxI9YUpyz4QgjhcixNcOm0UdS7Fve60s5zvIMgMhGSNkJUkvH/kfPAL7T3GD6Apy+sfHrEmigFXwghBsNmher8LsMx7X+q80HbjHPcvSFiNkxbYRT1juIeOBmU6vu6HcHshz81hnGC4oxiP0KBLUjBF0KIvmkNtaWdBb2juF86DdaOoRcFodMhah4k32H01iPnGcfch1Be5985ogW+Jyn4QgjRZOreW6842T7O3uVu/4BYo5c+/UudhT1ijjHs4iKk4AshJg5LE1w61T6+ntdZ2OvKOs/xCTKKefKmzuGYiLnGOLuLk4IvhBh/rG3t4+x53Yt7dT7QviS8u7fRQ59+fXt42jHOHtv/OLuLc8qCL/PwhRCAMXPlciGm1lBb0mM45gRcOtM5zq7cIHSG0VOff2dncQ+dBm7ujvm6HEQ2QBFCOKe+pim6e0PS7eDp0z4ccxJauoyzB07unOoYOc8IU8Nnu9Q4+0jobwMUp+zhCyEmMK2hpgDefaL3UgPWFsjOBJ/g9h77lzuLe+Rc8A1xTJtdhBR8IYTjaA2mC1B2DEqPGf8tO27MmumXgh9eGLfj7KNJCr4QYmxoDabznYW9tL24dywx4OZpDMEk3gqxV8HHL0D9xd7XCYqTYj9EUvCFECNPa2NGTM+ee8e8dncvYxgmaQPEXGUU+Mh5xnoyHbwmjfpSAxONFHwhxPDYbO0996NdCnx2Z5jq7mWMtyfdbhT2mI7i7nX5647BUgMTjRR8IcTA2WxQ/Xn3XnvZcWipNZ539zaKe8odnT33iMQrF/f+jPJSAxONFHwhRN9sNqg612NYJhta64zn3b0hOhlSvtyl554I7p6Obbfo15gVfKWUP/A7oBX4WGv96li9txDiCmxWo7h3DVTLs6G13njewweikmHBXV167nOluLuYYRV8pdTLwHqgQmud3OX4TcCLgDvwB631C8DtwE6tdZZSajsgBV8IR7BZofJs7567pcF43sMHolNgQVpnzz1ijhT3cWC4PfxXgN8Af+44oJRyB34LrAaKgYNKqbeAOCCn/TTrMN9XCDEQNitUnunRc8/pUtx9jeK+8O7Onnv4nKEt7Suc3rD+VrXWnyilpvY4vBQ4p7XOB1BKbQNuwyj+ccAxLrN5ulLqAeABgPj4+OE0T4jx5UrryljbjOLetedengOWRuN5T7/24n5PZ889fLYU9wlkNP6mJwNdd+MtBpYBvwZ+o5RaB2T192Kt9UvAS2CspTMK7RPC9fRcV8ZcBG89BBf+acyA6ei5t7U/7+kH0fNh0Ve69NxnT7jFwkR3Y/ZPu9a6AfjaQM6V1TKF6MJmhfef6r2uTFsLHHkFPP0hZj5cfV+XnvssKe6il9Eo+CXAlC6P49qPDZjWOgvIWrx48f0j2TAhnF7H2jKlR6DkSPvNTMc7Z8v0ouDJIinuYkBGo+AfBGYppaZhFPrNQPoovI8Qrq+21CjqJUeMIl96tHPhMHevztkyua9DU3Xv1wfFSbEXAzbcaZmZwPVAuFKqGHhGa/1HpdRDwHsY0zJf1lrnDfK6MqQjxp/G6s5ee0cPvr7ceE65G8sNzF0PkxdB7KLuyw9MWSrryohhkw1QhBgNLXVGkNrRay85Yqzx3iFsVnthX2gU9+gU8PK7/DWvNEtHiHb9bYDilAW/Sw///rNnzzq6OUJcnqXZmCHTtbhXnsG+d2pQPExuL+yxC41g1SfIoU0W45tLFfwO0sMXTsdqMbbV6xqqVpwAW5vxvH9k55DM5EXGjJlJEY5ts5hwZItDIQarY/Ewe3E/0j7Xvdl43ifI6LF/4eHOIh8YK5tzCKfllAVfQlsx5rSGmsLuPffSY50rQ3r6Gb31xd/oHHsPnS7FXbgUGdIRE1Ndee/pkI1VxnPuXsbKkLELO3vuEXNk+qNwGTKkIyauxur2HvvRziJfV2o8p9yMDTrm3NwZqkYldd9qT4hxwikLvgzpiCvqb4piS71xZ2rXGTOm852vC50BU7/YGapGp4CXv+O+DiHGkAzpCNfTcyExMG5cmhRl3MikbcaxwLjO6ZAdM2Z8gx3TZiHGkAzpCNdms8Kl01ByGN79Ye+FxLTVWHpgxQ86Q9VJkY5pqxBOSgq+cD49FxArOWIM03Rs2tGftha44ckxaaIQrsgpC76M4U8wdRe7FPfD7QuItS8U5u5tLP278J7OGTNbNxpj9z0FxY1tu4VwMU5Z8GV55HGs2dwZpnYU99r21bOVW/sCYus6i3tUUu+9VFc+IwuJCTEETlnwxThhaTLuTLUX9yPGnasdQqdD/DWdxT1m/sBmzHQsGCYLiQkxKFLwxciwtsGlk503MpUcNtac6VhjZlI0TL4aFmzunO/uFzr095t/pxR4IQZJCr4YPK2hOr9LcW8PVTv2U/UJMor6Fx/pnBIZGOvYNgshnLPgS2jrZGpLuxf30iPGWDyAh68xFLP4a53FXdaYEWJIdh8t4Zfvnaa0ponYYF8eXzuHDQsnj9j15cYr0Z19GYIjUNL+37oy4znlDlHzjKGZjuIekQjuTtlvEMKl7D5awo/e/xMq9B2UZw3aEoyuvpmfr/nqoIu+3Hglemtt7FyGoKPnXp3f+XzYTJh6nVHgO5Yh8PR1XHuFGMd+9rdXcYvciXKzAKC8atCRO/nZ3zzYsPAHI/IeUvDHg4FsfWe1wMW8Hht3nDTuUAUInGwEqQvvlWUIhBhDrW1W3s47RWPAG7i1F/sOys1Co38WIAVfQO91ZcxFxuO6cmNpgb427vANMYZk5tzcOTwTEOW4r0GICcTcYia7Iod3zh7kQOkxLracA49a3Pqpxm6eNSP23lLwXd2HP+29roylCf7ylPH/HRt3LPlm53z3kKkSqgoxBpramjhVfYqcSznkVuZyuDybiuYS+/O6NYIYnxS+EHcV75Vk0GA19bpGkNfIrQklBd/VNFZ3n+tuLur/3G/vk407hBgjFpuFc6Zz5FblkleZR05lDp/XfI61fdhUtQXT2jgZt9b5XBWZwpfnX8PN86bh42l8PpfmT+GpfzyDRbfYr+mpvHly+b+OWBul4Duz1gYjVC053FnkTRfan1QQPtvowVsae782aIoxo0YIMeJs2kZhbWG34n6q+hQtVqNYT/IMJIBpuNeupK4mGvfWBK6fOZ31K2JZOTcSf+/epXfd9HUAvHjkRcobyon2j+aRRY/Yj48Ep5yW2WUe/v1nz551dHPGRlsrVOR1n+9+6VTn2u5BU9q33Lu6M1T1Cex7bXhPX0j9tdyJKsQIudhwkdyqXHIrjT95VXnUte937OvhS2JoIrG+s6itiSE7P4CiCj883Ny4blY46+fHsjopikAfzyu8y8jpb1qmUxb8DuN2Hr7NZqwp0zEsU9Ieqrb3DvAN7SzsHfPdL7e2+0Bm6QghBsTcYiavKs9e3HMrc7nUdAkAD+XBrJBZJIUnkRKeQpDbdLLzfdibc5GzFfW4KfjCjHDWz49hbVI0If5eDvkapOA7itbGapBdh2VKj0FLrfG8pz/EXtW99x6cIKGqEGOgI1TtWtwL6wrtz08NnEpyeDLJ4ckkhSUxN3Qul2ptZGWX8vbxMk6U1aIULJkaSur8GG5KjiEiwPH7IcuNV2OlZ6hacgQaKozn3DyN5X5TvtzZe5dQVYgx0WZr41zNuW7F/VzNOXuoGukXSUp4ChtnbSQpLImk8CQCvQIBKDM3sSe7jGeyD3O8yJgmuTA+mKfWz2NdSgzRQT4O+7oGQwr+cAwkVJ25sn1Y5mqj2Hu6xg+GEK5Ma01hXWG34n6q+hTNVuNelECvQJLDk1kRt4KU8BSSw5OJ8Ivodo2KumZ2HbpA1vFSDhUY0yWTJwfy5M1zWTc/hrgQvzH/uoZLCv5ADSRUnbwIrv5a91BVCDHqKhoruhX33Kpce6jq4+5DYlgim2Zvshf3KQFTUH0Mm1Y3tPJubjlZx0v57HwVNg1zogL4/prZrJsfy7TwAezX4MSk4PelI1Tt2LSjZ6jqF2b02hNT2+9UlQ2zhRgNe/L39JqmeO3ka8mryiOvsjNYrWgyhk3dlTuzQmaxdupaksOMsfcZwTPw6O82VsDcZOH9vHKyssv457lKrDbN9Ah/HrpxFqnzY5gVFTBWX+6ok9BWa2N2S9c9VcuO9w5Vu86YkVBViFG3J38Pz376rH0YBkCh0HTWrITABCNUbS/uc0Ln4Otx5QX+6lva+ODERd7OLuWTM5W0Wm1MCfVl/fxY1s+PYV5MYJ+/AbiKiRPaXmmKYkeo2rX33jVUjU7uDFUnX22Mw0uoKsSoa7O18XnN5+RW5pJTmcObn79JW8eOae00mgDPAP7z+v9kXtg8gryDBnz9plYrH52q4O3sUj46VUFLm42YIB++ck0C6xfEsiAuyKWL/ECMWcFXSk0HfgwEaa03jcqb9LWQ2JsPwZn3jLH2ksNQU9DRos5QtWMBsehk8HD8lCohxjutNUV1RfbinleVx8mqk/befIBXQK9i36HeUs81sdcM6H1a2qz87fQl3s4u44OTF2lstRI+yZvNS6awfkEsV8eH4OY2vot8VwMq+Eqpl4H1QIXWOrnL8ZuAFwF34A9a6xf6u4bWOh/4hlJq5/CafBl9LSRmbYHcnZ2h6uKvS6gqxBi71HiJnMoc+12quZW51LYaw6be7t4khhqhasec9/iAeNa+vpayhrJe14r2j77se1msNv5xrpK3j5fxfl45dS1thPh5cttVk0ldEMOyaWG4T6Ai39VAe/ivAL8B/txxQCnlDvwWWA0UAweVUm9hFP/ne7z+61rrimG39krMxf08oeCx3FF/eyEE1LbWkleZZy/sOZU5VDR2hqozg2eyOmG1vbjPCJ6Bp1vvZQceWfRIrzF8H3cfHln0SK9zrTbN/vwq3s4u5Z3ccmoaLQT4eHBTcjTrF8TyhRlheLq7jd4X7SIGVPC11p8opab2OLwUONfec0cptQ24TWv9PMZvA2MvKK7v1SOD4sa+LUJMAM1tzZyqPkVelbGAWF5lHhdqL9ifjw+IZ3HUYntxnxs6d0ChKlx5MTGbTXOowMTb2aXszSmjsr4Vfy93Vs+LYv38WK6bHY63h+RvXQ1nDH8y0LW6FgPL+jtZKRUG/AxYqJR6sv0fhr7OewB4ACA+Pn5wLVr5dN8Lia18enDXEUL00hGqdi3uZ01nadPGWHuEbwTJ4cmkzki1L0UwmFC1LxbzVTSce4K6miYCgn1pnT6bo4Umso6XsTenjPLaZnw83Vg5N4r182O4YW6kfblh0duYhbZa6yrgwQGc9xLwEhjTMgf1Jh2zcWQhMSGGRWtNcV0xuVW59uJ+svokTW1GZyrAM4Ck8CTuS77PPi0yyn9kd03bfbSEJ9/IocliLH1QUtPEv+44jga83N340pwInpw/l1WJUX0uNyx6G853qQSY0uVxXPuxYeuyPPLgXzz/TinwQgxSZVOlsStT+/ruuVW5mFvMgBGqzg2dyx2z7iApPInksGTiA+NxU6M7Jv7CO6fsxb6DBoL9PPnkBzeM6XLD48VwCv5BYJZSahpGod8MpI9Eo7TWWUDW4sWL7x+J6wkhOtW11tkD1Y7NOy42XgSMUHVG8AxWxa+yF/eZITP7DFVHg82m+fu5SrYdKKS8trnPc8yNFin2QzTQaZmZwPVAuFKqGHhGa/1HpdRDwHsYM3Ne1lrnjUSjhtXDF0LYtVhbei3/2zVUnRIwhUVRi0gOSyYlImVQoepIKjc389qhIrYfKqLY1ESovxeTvN2pb7H2Ojc2eOzbN17I0gpCuIi+1pXpuv2d1Wblc/Pn9l57bmVut1A13DfcPt6eEp5CUvjwQ9XhaLPa+NuZS2QeKOKjUxexafjizDA2L4lnTVIU7+SUdxvDB/D1dOf521PYsHCyw9o9msxZWVT86r9oKyvDIyaGyMceJSg1ddDXcakNUCbkFodCXEZf68p4u3tz+8zb8XT3JLcyt1eoOi98nrE6ZFgySeFJRPlFOcXSASU1TWw/WMSOg0WU1zYTPsmbLy+OY/OSKSSEdV+NcvfREn753mlKa5qIDfbl8bVzxnWxL3vqaXRzl7WDfHyI+befDrrou1TB7yA9fCEMq15bZR9n78nLzYu5YXONXntYEsnhySQEJox6qDoYFquND09WsO1gIX87Y2wXuGJWBGlL41mZGCk3RQFnb1xJW2lpr+MesbHM+ujDQV1r4iyeJoSLq2+t7wxVq7qHqj0pFPvv3j9moepgFVY1su1gIa8dLuZSXQvRgT5894aZ3LlkiktuIDJabK2tfRZ7gLay3stLDJVTFnwJbcVE0WJt4XT16W4bd1wwX7AvARw3KY6FEQv5p+Wf9rVnuor2j3a6Yt/SZuUvJy6SeaCQf56rwk3BjXMjSVsaz5dmR+AhvXk7S0kJpu07qNnZ/xJjHjExI/Z+TlnwZVqmGI+sNiv55vxuxf2M6Yx9VcgwnzBSwlO4Zdot9uGZYJ9goO8x/P7WlXGUzy/Vs/1gETsPF1Pd0MrkYF++t3o2X148xWX2fB0L2maj4dN9mDIyqP/4YwAm3XgDXlOnYdq6tdcYfuRjj47YeztlwRfC1WmtKakvIbcql9xLRnE/UXXCHqpO8pxEUlgSX5n3Ffu2e5cLVa+0royjNFusvJtbTsaBQg6cr8bDTbEqMYq0ZfFcOzN8wq5K2Rer2UzNrl3UZG6jtaAA99BQwu6/n5C77sQzNhYAnzmzR2SWTn+cMrSVWTrC1VQ1VdnH2ztuaDK1GBtfe7l5MTd0rn0BsaTwJKYGTnWqUHWwTpfXkXmgkF1HSzA3WUgI82PzknjuuHoykQHSm++q+eRJTBkZmLPeRjc347twISHp6QSsXYObl9eovKfM0hFihNS31nOi6oTRe28fnulYt91NuTE9aLq9154cnsys4Fl4ujvXOPtQNLa28XZ2GdsOFHKksAYvdzfWJkeTtmQKy6eHTaiNRK7E1tpK3XvvYXo1g6Zjx1A+PgSlrickPR2fxMRRf3+ZpSPEELRaW41QtUtxP28+3y1UXRCxgLsT7yY5PJnE0ET8PMfX7JO8UjOZBwp582gpdS1tzIjw5yfrErl9URyh/qPTQ3VVltJSI4R97TWs1dV4JSQQ9eQTBG3ciHug4zdckoIvRDurzcp58/luxf206XS3UDU5PJmbpt1kD1VDfEIc3OrRUd/SxlvHStl2sJDsYjPeHm6sS4lh89J4lkwNcYobuJyFttlo2LcPU0Ym9X/9KwCTbriBkPQ0/K+5BuXmPEN3TlnwZVqmGG1aa0obSrutMXOi6gSNbY0A+Hv620PVjuUIov2jx3Wh01qTXWz05t86Xkpjq5U5UQE8mzqPjQvjCPJz/WGpkWStrcW8axemzG20XrhghLDf/KYRwk52zruBZQxfTAgdoWrXAt8Rqnq6eXYLVZPDkpka5Nqh6mCYmyy8eayEzANFnCyrxdfTndQFRm9+4ZTgcf2P3FAYIWwm5rffRjc14XvVVYSkpxFw002jFsIOlozhi3HlcguJNVgajFC1S3EvbTDuYlQoZgTP4EtTvkRymFHgZ4fMHheh6mBorTlcYCLzQBF7ckpptthIig3kuQ3J3HZVLAGy/HA3Rgj7PqaMDJqOHkX5+BC4fh0haWn4JiU5unkDJj184XL6ugnJw82D+eHzMbeYyTfn20PVyZMm23vtyeHJzAubN+5C1cEwNbTyxtESth0o5GxFPZO8Pbj1qljSlsSTEue4lTOdlT2E3bkTa1UVngnxhKSlEbxxI+5Bzvv9kh6+cHlWm5ULtRd4/sDz3Yo9GPutHqs4xrVx17J26lr7fPdQn1AHtdZ5aK3Zn1/NtoOFvJNbTmubjaumBPOLO1JYPz9WtgfsQWtN4759VGdkUP9Rewh7/fWEpKXh/8UvOFUIO1jyNy2cktaasoaybssQ5FXm2UPVPl+D5rcrfzuGrXRulfUtvH64mO0Hi8ivbCDAx4O0JVPYvDSexBjHTxF0NtbaWsy7d2PKyDRC2JAQwr7xDYLvuguvOOcMYQfLKQu+zNKZeKqbq+13qHZMi6xurgaMUHVOyBxunXEryeHJvHjkRS41Xep1jWj/6LFuttOx2TT//LySbQeKeP9EORarZsnUEP7lhpnckhKDr5e7o5vodJpPnTJC2KwsI4RdsIDYf/8FAWvX4ubt7ejmjSinLPiyeNr41mhpJK8qr1txL6kvAYxQdXrQdK6bfJ191szskNl4uXfOfvBw83D6hcTGWkVtM68dLmbbwUKKqpsI8fPkK9dMZfOSKcyKCnB085yObm2l9v2/GCHskSMob28jhE1Pd6kQdrCcsuCL8cNitXDGdIbcylxyKnPIq8oj35yPTdsAiPWPJTk8mbvm3GUPVf09/S97TWddSGy09dz96XurZxPi70XmgUI+PFWB1aa5ZnoY318zh7VJ0fh4Sm++J0tZGabt26l5rT2EjY8n8oc/JHjjBtyDgx3dvFEns3TEiLFpGxfMF8ityiXnklHcT1WfwmKzABDqE2rfkSk5PJmksCTCfMMc3GrXsPtoSa/9XRWggfBJXtxxdRybl8QzLfzy/1hORB0hrCkzk7oPPwKtjRA2PQ3/L37RpUPY/sgsHTGitNaUN5Qbq0O2B6p5VXk0WBoA8PPwY17YPO5JvIekcKPIx/rHyk08Q/Tv753qVuzBKPahfp58+sRKvDzGX9EaLmtdHeZduzFlZtJ6/jzuwcGEfePrBN+1edyEsIMlBV8MiKnZ1G22TE5ljj1U9XDzYE7IHNZPX2+f8z4taBrubjKkMFyV9S1sP1hEaU1zn8+bGi1S7HtoPn0a06sZ9hDWZ8F8Yn/xgnEn7DgLYQdLCr7opdHSyImqE93Wd+8aqk4Lmsa1k68lOTyZlPCUXqGqGB6tNYcKTGzZV8A7uWVYrBovDzda22y9zo0N9nVAC52Pbm2l9i9/wZSRSdPhw50hbFo6vsnjN4QdLKcs+DItc+xYrBbO1Jyx99pzK3O7haox/jEkhydz55w7SQlPITE0kUlekxzc6vGprtnC7qMlbN1fyOmLdQT4eHD3sgTuWR5PbkltrzF8X093Hl87x4EtdjxLeXlnCFtZaYSwP/gBwbdvnBAh7GBJaDuB2LSNC7UX7MU9r9IIVVttrQCEeIfYx9s7lv+VUHX0nSyrZev+AnYfLaGh1Ury5EDuXZ5A6oJY/Lw6+2Q9Z+k8vnYOGxZOvLForTWN+/djysik7qOPwGZj0pe+RMjd6eM2hB0s2fFqHOtrIbFbpt3CxcaL3bbcy6vKo95SD4Cvhy/zwuYZhT08ieSwZCZPmiyh6hhpabPyTk45W/cXcKjAhLeHG+vnx3LvNQksiAuSv4d25qyszj1eo6LwW7aM5pwcWvPzcQ8OJnjTHQRv3oxXXJyjm+pUpOCPU30tJOam3PB196WhzZgx4+HmweyQ2fZee3J4MtODpkuo6gBF1Y28+lkhOw4VUd3QytQwP+5ZnsCmq+MI9pMcpCtzVhZlTz2Nbu4eWHtMmULEv3yHwJtvnvAhbH9kWuY40mhp5GT1SXIrc/nN0d/0WkjMpm3YsPHk0ieNUDV0Nt7u8sFwFKtN8/HpCrbuL+DjM5dQwOp5UdyzPIEvzgiXvWD7oFtbufjzn/cq9gBYrQRv2DD2jRoHpOA7OYvNwlnT2W6LiH1e87k9VO1Pc1sz6YnpY9RK0ZeOKZUZnxVSUtNEZIA3371xFmlLpxATJLNr+mIpL6dmxw5Mr72G1VTT5zltZWVj3KrxQwq+E7FpGwW1Bd2K+6mqzlA12DuYpPAkVsavJDnMWP43fU86ZQ29PwCykJhjaK05cL6arZ8V8m77lMovzAjjx+sSWT0vCk93CRR70lrT+NlnmF7N6AxhV6ygKScHa3V1r/M9YmIc0MrxQQq+g2itudh4sVtxP1F5gjpLHWCEqomhiaTNTbOv7R43Ka5XmPfIokdkITEnUNdsYdfRErbuL+DMxXoCfDy4Z3kCdy9LYGakTGPti7W+HvPuN407YT//HPegIELv+yohmzfjNWVKn2P4yseHyMcedWCrXZsU/DFibjF323IvtyqXyqZKADyUB7NCZnHztJvt68wMNFSdqAuJOYsTpbVs/cyYUtnYaiVlchC/uCOl15RK0an5zBlMGRmY38pCNzbik5JCzPPPE3jzTbj5+NjPC0pNBeicpRMTQ+Rjj9qPi8Eb01k6SqkNwDogEPij1vr9y53vqrN0Gi2NnKo+1a24F9UV2Z+fFjTNPiSTEp7CnNA5Eqq6kGaLlXdyy9iyr4AjhTV4e7iRuiCWe5cnsGCK3OzTF93aSt0HH2DKyKTx0CGUlxeBt9xCyN3p+KakOLp5486wZ+kopV4G1gMVWuvkLsdvAp9tT0QAACAASURBVF4E3IE/aK1f6O8aWuvdwG6lVAjwH8BlC74rsNgsnDOdsy/9m1uZy7mac/ZQNdo/muSwZO6YdYd9+d8AL1mf3BUVVjXy6oECdhwswtRoYVq4Pz9ZlyhTKi/DcvEiNdt3YHptB9ZLlXjGxRH5+PcJuv12PEJCHN28CWcwv3O+AvwG+HPHAaWUO/BbYDVQDBxUSr2FUfyf7/H6r2utK9r//yftr3MpNm2jsLawW3E/VX2KFmsLAEHeQSSHJXPDlBvsQzPhvuEObrUYDqtN89dTFWzZX8AnZy/hphSrEiO5d/lUvjAjTKZU9sEIYQ9gysig7sMPwWbDf8V1hKan43/ttSh3uf/DUQZc8LXWnyilpvY4vBQ4p7XOB1BKbQNu01o/j/HbQDfKSBxfAN7RWh8ZaqPHysWGi/YhmY67VXuGqh0bdySHJRMX0DtUFa7pUl0L2w8WknmgyD6l8uEbZ5G2NJ7oIJ8rX2ACstbXY36zPYQ91zuEFY433FRpMlDU5XExsOwy538XWAUEKaVmaq1/3/MEpdQDwAMA8fHxw2zewJlbzPYt9zrWmenYN7UjVL1p2k3dQlUPNwnlxhOtNZ+dr2br/gLeyzP2g/3izDB+si6RVTKlsl8tZ89SnZFB7ZtvYWtsxCc5mZif/5zAW27uFsIKxxvTiqW1/jXw6yuc8xLwEhih7WDfo691ZXrOWGlqa7KHqh3FvbCu0P781MCpLItZZi/uc0Lm4OMhP7jjVW2zhV1HjCmVZyvqCfTx4N7lU7l7eTwzImRKZV+0xdIZwh48KCGsixhuwS8Buv6uFtd+bFiGujxyz3VlyhrKePbTZymtLyXYJ9jowbeHqlZtLDMb5RdFcngyG2dttIeqgV6Bw/0ShAvIKzWzdX8hbx4zplTOjwvi3++YT+qCWHy9ZJy5L5aLFdTs2EHNjh20XbokIayLGdS0zPYx/Lc7ZukopTyAM8BKjEJ/EEjXWueNROMGOy1zzc41fd512iHQK9Dea08OM/4b4RcxEk0VLqLZYmVvThlb93dOqbx1QSz3yJTKfmmtaTxw0AhhP/jACGGvu5aQ9HQmXXedhLBOaCSmZWYC1wPhSqli4Bmt9R+VUg8B72HMzHl5JIr9UHv45Q3l/T63d+NeCVUnsIKqBjLaV6k0NVqYLlMqr8ha34D5rTepycyk5ew5I4T96lcJ2XwXXmOYr4mRM5hZOmn9HN8L7B2xFhnXzAKyFi9efP9gXhftH91nDz/GP4YpgTJLYKKx2jQfdUypPHMJdzfF6sQo7r0mgS/MCJN//PvRcvYspsxMzLvfNELYpCRifvYzAtfdIiGsi3PKaSZD7eHLujICoKKumR3tq1SWmpuJCvTm0VWz2LxEplT2R1ss1H34oRHCHjhghLA332yEsPPnO7p5YoSMuw1QBjJLR7i+ntv9fX/NbGKCfdmyv4D3cstps2munRnOPcvjWZkoUyr7Y6mooGbHa0YIW1GB5+TJhKRtJuiOOySEdWGy45UYN3YfLem1obcCNBDk68mmq+O4e1k802VKZZ+01jQePGjsCfvBB9DWhv+K6whJS2PSihUSwo4DLrXj1VCHdMTE8Mv3Tncr9mAU+2A/T/Y9sVKmVPajZwjrFhRE6L33EpK2WULYCcIpC/5QQ1sxvrW0WXk3t5ySmqY+nzc3WqTY96Hl3DlMGZmY33wTW0NDZwh7y824+crOWxOJUxZ8IboqqWni1f0FbD9YRFVDK+5uCqut91BkbLAUrw5GCPsRpowMI4T19CTwlpsJSU/HZ/58maE0QTllwZchHWGzaf5+rpIt+wr46NRFAFYmRnHv8gSq6lr40e7cbsM6vp7uPL52jqOa6zQsFRXUvPYaNdvbQ9jYWCK+968E33EHHqGhjm6ecDCnLPgypDNxmRstvHa4iK37C7hQ1UiYvxffvn4GaUvjiQvxs5+n3FS3WTqPr53DhoWTHdhyx9Fa03ToENUZGdT9pT2Eve46op99lklfkhBWdHLKgi8mnpxiM1v2X+Ct46U0W2wsTgjhsdWzuSk5Gm+P3gVrw8LJE67Am7Oyum33F/7tb0ObBVNGJi1nz+IWGEjoPfcYIWxCgqObK5yQFHzhMM0WK3uyy9iyv4BjRTX4erqzcWEc9yyPJyk2yNHNcyo9N/RuKy2l/KmnAPCZN4+Ynz1H4C23SAgrLksKvhhzRdWNbP2sc6vA6RH+PJM6jzuujiPQx9PRzXNKFf/vr+zFviv38HCmvr5TQlgxIE5Z8CW0HX9sNs3fzlxiy/4C/nq6Ajcl69oMREcI21bW9yqw1qoq+d6JAXPKgi+h7fhhamhlx6Eitn5WQFF1E+GTvPnuDTNJWxZPTJAMP/RFa03T4cOYMjKoff8v0NaG8vZGt7T0OtcjJsYBLRSuyikLvnB9x4pq2LKvgKzsUlrbbCydFsoP1s5lbVI0Xh6yrk1fbA0NmLOyjBD2zBkjhL37bkLSNtOUk9NtDB9A+fgQ+dijDmyxcDVS8MWIabZYeet4KVv3F5BdbMbfy507F8dx7/KpzIkOcHTznFZLfr5xJ+zu3djq6/Gel0jMc/9G4Lp19hDWa+pUgG6zdCIfe5Sg1FQHtnzkWCwWiouLae4jpxD98/HxIS4uDk/PgWVfsniaGLYLlQ28+lkBOw4VY26yMCtyEvdek8DGhZMJkBC2T7qtjbqPPsKUmUnjvv0oT08CbrqJkPQ0fK+6asKNy58/f56AgADCwiTPGSitNVVVVdTV1TFt2rRuz8niaWJEWW2av56q4M/tm4t4uCnWJkVz7zUJLJsWKh/afrRdukTNzp2Ytu+grbwcj9gYIh57jOBNd+ARFubo5jlMc3MzU6dOlZ+bQVBKERYWxqVLlwb8Gqcs+BLaOq+q+ha2Hyri1f2FlNQ0ERXozWOrZrN56RSiAmVzkb5orWk6cgTTqxnU/uUvYLHg/8UvEv3UT5h0/fVyJ2w7KfaDN9jvmVMWfOFctNYcKaxhy74L7M0pp9Vq45rpYfxkXSKr5snmIv0xQti3MWVm0nL6NG4BAYSmpxG8eTPePX4FF2IsSMEX/WpsbeOtY6X8eV8BJ8pqCfD2IH1ZPPcsj2dmpISw/WnJP2/sCbtrlxHCJiYS/W8/JWjdOtz8/K58ASFGiRR80Uv+pXq27C9g5+Fi6prbmBsdwM82JrPhqsn4e8uPTF90Wxt1f/2rsRzxvv3g6Ung2rWEpKfju3DihbDjUUNDA9/5znfw8vLi+uuv5+6773Z0kwZNPr0CgDarjQ9PVbBlXwH/OFeJp7vi5uQY7r0mgcUJIVKw+tFWWdkZwrZPl4x49FGCv7xpQoewrup///d/OXbsGP/zP/9jP5acnMxrr73GoUOH2LRpE6mpqdx1111S8IXruVTXwrYDhWQcKKTM3ExMkA/fXzObO5dMITJAQti+aK1pOnrUCGHff98IYb/wBaJ//CMjhPWQj5WrysnJYdGiRfbHzc3NXLhwgdmzZ7N7925SUlIAcB9i0K61RmuNm5tbn48H+rqhcsqfTJmWObq01hy8YGLL/gLezS3DYtVcNyucZ29NYuXcSDwkhO2TrbGxM4Q9dQq3gABC0jYTsjkN7+kSwo6l3UdLRmU/hOzsbL72ta/ZH+fk5DB79mzc3d2Ji4ujuLiYq666CpvN1ufrt27dyq9//WtaW1tZtmwZv/vd7ygqKmLt2rUsW7aMw4cP87vf/Y4HHnjA/njv3r28/vrrvPzyywB885vf5NFHH+XChQvdXrd3714Shrnstdx4NYE0tLSx62gJW/cXcKq8jgAfD7589RTuXh7PjIhJjm6e02rJP49pWybmXbux1dXhPXcuIelpBK1fLyHsCDl58iSJiYkDOnf30RKefCOn145nz9+eMuyiHxYWxqRJk+xDmPX19axfv55XXnmFhoYGHnroIXx8fLj22mt7DemcPHmSH/zgB7zxxht4enryne98h+XLl7NixQqmT5/Op59+yvLly7lw4UK3x4cPH+a+++5j//79aK1ZtmwZW7duJSQkpNt5/enre+dSN16JkXX2Yh1b9xfw+pES6lvamBcTyAu3p3DrVbH4ecmPQF90Wxv1H3+MKSODhk/3dQlh0/BduFAyjVH0/2TlcaK0tt/njxbW0Grt3sNuslj5wc5sMg8U9vmaebGBPJOadNn3LSoqIiIiglOnTtmPPfTQQ/a7WP39/fm///u/fl//4YcfcvjwYZYsWWK0qamJyMhIVqxYQUJCQrei3fXxP/7xDzZu3Ii/vz8At99+O3//+9+59dZbe71uuOTTPg709evtuvkx/OXERbbsK2BffhVe7m6smx/DPcsTWBQfLAWrH21VVdS8thPT9u1GCBsdTcSjjxC8aRMe4eGObp6AXsX+SscHKicnh6Sk7v8onDhxgttuu21Ar9da89WvfpXnn3++2/ELFy7Yi3mHno/7M9DzBkoKvovr+ettSU0T33/tOE+/mUNts5XJwb784KY53LV4CmGTvB3cWudkhLDHjOWI33uvPYS9hqgfPUnADTdICDvGrtQT/+ILH1FS09Tr+ORgX7Z/65ohv292djbz5s3rdiwvL88e1F7JypUrue2223jssceIjIykurqaurq6K77uuuuu47777uOJJ55Aa82uXbvYsmXLkL6GK5GfZBf3y/dOdxvLBGizaVraNH/4ymJumBuJu5v05vtia2zE/PbbmDK30XLyJG6TJhGyeTMhaZvxnj7d0c0T/Xh87Zw+x/AfXztnWNfNyclh/fr19sfV1dVorYmOjh7Q6+fNm8dzzz3HmjVrsNlseHp68tvf/vaKr1+0aBH33XcfS5cuBYzQduHChVy4cGHIX0t/JLR1YXXNFlKefb/P5xRw/oV1Y9sgF9Fy/jw127ZR88YuI4SdM4eQ9HSC1q/DbYR/hRYDM5jQFkZvlo4rktB2nDtzsY4t+wp440hxv+fEBstuUl1pq9UIYV/NoOHTT40QdvVqQu5Ox3fRIsk0XMyGhZMnbIEfDin4LsJitfFeXjlb9hXw2flqvDzcWD8/hoQwP37/cf6I/3o7XrRVVVGz83VM27fRVtoewj7ysBHCRkQ4unlCjKkxK/hKqUTgESAc+FBr/T9XeIkAys3NZB4oJPNAIRV1LcSF+PLEzXO5c/EUQv29AEgI9Zdfb7vQWtN07BimjEzq3n0XbbHgd81yop54goAbb5QQVkxYA/rJV0q9DKwHKrTWyV2O3wS8CLgDf9Bav9DfNbTWJ4EHlVJuwJ8BKfj90FqzP7+aLfsv8F7eRWxa86XZEbxwTQJfmt07hJ2Iv96as7J6bfcXsGpVewibScsJI4QNvusuQtLTJIQVgoH38F8BfoNRqAFQSrkDvwVWA8XAQaXUWxjF//ker/+61rpCKXUr8G1gdOYcubi6Zgu7jpawZV8BZyvqCfbz5BvXTuPuZfEkhEmY2MGcldVtQ++20lJKn3gSPJ+G5ma8Z88m+tlnCUpdLyGsEF0MqOBrrT9RSk3tcXgpcE5rnQ+glNoG3Ka1fh7jt4G+rvMW8JZSag+QMdRGjzdnLtbx530X2HWkhIZWKymTg/j3TfO5dUEsPp6yG1JPFb/6L3uxt7NaUZ6exG/dgu/VV0sIK0QfhjOYORko6vK4GFjW38lKqeuB2wFvYO9lznsAeAAgPj5+GM1zbn2FsKnzY7n3mgSumhLs6OY5rbbqatpKS/t8Tre04Le410w0IUS7MUuvtNYfAx8P4LyXgJfAmIc/uq0aewMJYUV3Wmuajx+nOiODunfe7fc8j5iYMWyVEK5nOAW/BJjS5XFc+7FhG2/LIw82hBUGW1MTtXv2UJ2RYYSw/v4E33knHjHRVP7mt92GdZSPD5GPPerA1grh/IZT8A8Cs5RS0zAK/WYgfSQapbXOArIWL158/0hcz1EkhB2a1oICTJnbqNm1C5vZjPesWUQ/+wxBqan2ENYzKqrXLJ2g1FQHt1yMZ/n5+fzsZz/DbDazc+dORzdnSAY6LTMTuB4IV0oVA89orf+olHoIeA9jZs7LWuu8kWiUq/fwT5fXsWW/hLCDoa1W6v/2ibEc8T/+AR4eBK5ZTUhaGr6LF/cKYYNSU6XAixF3uS0OExMT+eMf/8imTZsc2MLhGegsnbR+ju/lMgHsULliD78jhP3zvgIOSAg7YG3V1dTsfJ2abduwlJbiERlJ+HcfIvjLX8YzMtLRzRMTzOW2OBwJssVhH1yph19ubibjQCHbJIQdMK01zdnZxnLEe98x7oRdtozIH/6QgBtvQHl6OrqJwtll74APfwrmYgiKg5VPw/w7h3/Zy2xxOBDOvsWh/V8OZ/xz9dVXa2dks9n0P89d0t/eekhPf3KPnvrE2/qrL3+mPzxZrtusNkc3z2lZGxu1aedOnb/xdn1izlx9atHVuuyn/6abz551dNOEg504cWLgJx/frvVzUVo/E9j557ko4/gwhYaG6vj4eJ2QkKATEhJ0WFiY/upXv6q11rqyslJ/61vf0tOnT9c///nP+/wa1q9fr1tbW7XWWn/729/Wf/rTn/T58+e1Ukrv27dPa617PT506JBOTk7W9fX1uq6uTs+bN08fOXKk13n96et7BxzSfdRUp+zhOysJYYemdwg7k+hnniYw9VbcJ8n3TfTwzhNQntP/88UHwdrS/ZilCd58CA7/qe/XRKfAzf2u/AJceYvDsLAwfv/73/f7etnicIicbUinZwg7Py6IX26aT6qEsP3SViv1n3yCKSOThr//HTw8CFi9itD09D5DWCEGrGexv9LxAZItDh1EO0FoKyHs0LSZTNTs3EnNtu1YSkqMEPah9hA2SkJYMQBX6Inzq2QwF/U+HjQFvrZnyG8rWxxOQB0hbOaBQi5JCDsg2h7CZlL7zjvo1lb8li4l8vHHCVh5o4SwYmStfBqyHjaGcTp4+hrHh0G2OHSwsdriUGvNvvwqtu4v6HYn7FfkTtjLsjU3U7tnL6aMDJrz8nDz8yNowwZjT9hZsxzdPOFCBrvF4WjN0nFFLr/F4ViN4dc1W3jjSAlb9hdwTkLYAWstLDRC2DfewGY24zVzBlFPP0XQrbdJCCvGxvw7J2yBHw6nLPijPYYvIezg9RnCrlpFSHoafkuWSAgrhAtwyoI/GiSEHZo2kwnz669jytxmhLARERLCCuGixl3B3320pNv+rvevmEZ1g0VC2EFqys7G9GpGjxD2+wSsXCkhrBAuyikL/lDH8HcfLeHJN3JoslgBKKlp4tm3TgBw/RwJYa/EHsJmZtKcm4ubnx/Bm+4gJC1NQlghxgGnLPhDHcP/5Xun7cW+q6hAb1752tKRat6401pYiGnbdsyvv45VQlghxi2nLPhDVVrT1Ofxitrh3YE3Hmmrlfq//91Yjvjv/wA3NwJWG8sR+y2VEFaI8WhcFfzYYF9K+ij6scG+DmiNc7KHsNu2YykuNkLY73yH4DvvlBBWiHFuXBX8x9fO6TaGD+Dr6c7ja+c4sFXOoSknxwhh9+41QtglS4j8/vckhBViAnHKgj/U0HbDwskA3WbpPL52jv34RGNrbqZ27zvGnbDtIWzQHbcTkpaGzwht6CDERDAetjcEWVphXDBnZXXb3zX0K/fSdukS5p3tIeyMGYSkpxF02224T5rk6OYK0cugl1YYJR9++CGvvPJKv4uXbdq0yekKvssvrSAGzpyVRdlTT6ObmwFoKy2l4oVfgFIErFlDSHq6hLBCDNDx48dZuHCho5sxaoa3QaJwuIr//E97se/KIzKSuBf/C/9lS6XYi3FnT/4e1uxcw/w/zWfNzjXsyR/6sshdHT9+nPLyclasWEF8fDwffPDBiFzXWUjBd1FNOTmUPvkj2sov9vl8W0XFGLdIiLGxJ38Pz376LGUNZWg0ZQ1lPPvpsyNS9I8fP05ERASffPIJL774Iq+++ioAVVVVPPjggxw9erTXBieuRIZ0XIitpaUzhM3JQfn5ofz80I2Nvc71iIlxQAuFGL5fHPgFp6pP9ft89qVsWm2t3Y41W5t5+p9Ps/NM3+Prc0Pn8sOlP7zs+1osFqqqqvje975nfxwcbKyzdaXtDV2FFHwX0FpcjCkzE/Prb2CtqcFrxgyifvITgjbcRv1f/9ptDB9A+fgQ+dijDmyxEKOnZ7G/0vGBOnnyJAsWLMDNzRj4yM7OJjk5eVjXdDZOWfCdbU9bR9A2Gw3/+AemVzOo/+QT407YlSuNELbLuHxQaipAt1k6kY89aj8uhKu5Uk98zc41lDWU9Toe4x/D/930f0N+3+PHj7NgwQL74+zs7AHvZ+sqnLLgO8Oeto5iramh5vU3MG3bhqWoCPfwcMK//aBxJ2w/W6UFpaZKgRcTxiOLHuHZT5+l2dr5W62Puw+PLHpkWNc9fvw4S5YssT/Ozc2VHr4YHU05uZgyM6ndswfd0oLv4quJfOxRAlatQnnJMs5CdFg3fR0ALx55kfKGcqL9o3lk0SP240P1H//xH90e5+fnD+t6zkgKvgPZWlqofecdTBmZNGdnozr2hE1Px2eO3AkrRH/WTV837AI/EUnBd4DW4mJqtm2jZufrRgg7fTpRP/4xQRtuwz0gwNHNE0KMU1Lwx4g9hM3IpP5vfzNC2BtvNPaEXb5cbo4SQow6KfijzFpTQ80bu4wQtrAQ9/Bwwh78FiF33dVvCCuEEKNBCv4oacrNw5SR0RnCXn01EY88TODq1RLCCiEcQgr+CLK1tFD37rtUZ2TQfDwb5etrhLBpm/GZO9fRzRNCTHBS8EdAa3EJNdvbQ1iTCa9p04j60Y8I2rhBQlghhNMY04KvlPIH/gY8q7V+eyzfe6Rpm42Gf/7TCGE//thYjnjljcadsBLCCiGc0IAKvlLqZWA9UKG1Tu5y/CbgRcAd+IPW+oUrXOqHwI4httUpWGtqqNm1G9O2TCwFhbiHhRH2rQeMEFYWLBNCOLGB9vBfAX4D/LnjgFLKHfgtsBooBg4qpd7CKP491w/9OrAAOAH4DK/JjtGU1xHC7kU3N+O7aBER332YwDUSwgox3u3evZs9e/ZQW1vLN77xDdasWePoJg3JgAq+1voTpdTUHoeXAue01vkASqltwG1a6+cxfhvoRil1PeAPzAOalFJ7tda2Ps57AHgAID4+fsBfyGjoCGFNGZk0HT9uhLCpqYTcnS4hrBDjUH9bHG7YsIENGzZgMpn4/ve/P74Lfj8mA0VdHhcDy/o7WWv9YwCl1H1AZV/Fvv28l4CXwNjTdhjtGzJLSQmmbdup2bnTCGGnTiXqR08StGED7oGBjmiSEGIMXGmLw+eee45/+Zd/GcMWjawxn6WjtX7lSuc4YnlkI4T9FFNGhnEnLDDpxhsISUvD/5prUG6yOZgQzsKclTUqS4IfP36cqKgoVqxYwYULF3j55ZdZtWoVWmueeOIJbr75ZhYtWjQCX4FjDKfglwBTujyOaz82bGO5PLLVbKZm1y5Mme0hbGgoYfffT8hdd+IZGzvaby+EGCRzVla3TX/aSkspe+ppgGEX/ePHj3P33XfzySefsGvXLl599VVWrVrFf//3f/PBBx9gNps5d+4cDz744LC/DkcYTsE/CMxSSk3DKPSbgfSRaNRY9PCbT5ygOiOD2rf3GCHswoVEPPRdAtauwU1CWCEcpvznP6flZP9bHDYdP45u7b67lW5upuzHP6Fmx2t9vsY7cS7RP/rRZd/3clscPvzwwzz88MOD+TKc0kCnZWYC1wPhSqli4Bmt9R+VUg8B72HMzHlZa503Eo0arR6+rbW1M4Q9dgzl40NQ6npjOeLExJF8KyHEKOlZ7K90fKBki8N2Wuu0fo7vBfaOaIuGqa+xPb9FizBt32GEsNXVeCUkEPXkEwRt3CghrBBO5ko98bM3rqSttLTXcY/YWBK2/LmPVwyMbHHoIEMd0ulrbK/0h0+AzQZubky64QZC0iWEFcKVRT72aLfPOYDy8SHysUeHdV3Z4tBBhjqkU/Gr/+r2QwCAzYbbpElMf3M3npMnj2ArhRCO0BHMjvQsHdni0EGG2sNvK+u9kz2AraFBir0Q40hQauqITMOcaJxyXENrnaW1fiAoKGhQr/PoZy2b/o4LIcRE4pQFf6giH3sU5dN9qZ6RGNsTQojxwCmHdIZqtMb2hBCjT2sty4oPktaDW33GKQv+cG68krE9IVyPj48PVVVVhIWFSdEfIK01VVVV+PgMfAFiNdh/IcbS4sWL9aFDhxzdDCHEKLNYLBQXF9Pcc5aduCwfHx/i4uLw9PTsdlwpdVhrvbjn+U7ZwxdCTCyenp5MmzbN0c0Y98ZVaCuEEKJ/TlnwlVKpSqmXzGazo5sihBDjhlMW/KHOwxdCCNE/pw5tlVKXgIIeh4OAwXT9B3P+QM4NByoH8f7jwWC/56NprNoyku8z3GsN9fUj/bM/mHPlc+JYCVrriF5HtdYu9Qd4abTOH8i5wCFHfw+c/Xs+Htoyku8z3GsN9fUj/bM/mHPlc+Kcf5xySOcKskbx/MFee6Jwpu/LWLVlJN9nuNca6utH62ffmX4enInTf1+cekjHGSmlDuk+5rcKITrJ58Q5uWIP39FecnQDhHAB8jlxQtLDF0KICUJ6+EIIMUFIwRdCiAlCCr4QQkwQUvBHkFJqg1Lq/1NKbVdKrXF0e4RwRkqp6UqpPyqldjq6LRONFPx2SqmXlVIVSqncHsdvUkqdVkqdU0o9cblraK13a63vBx4E7hrN9grhCCP0OcnXWn9jdFsq+iKzdNoppVYA9cCftdbJ7cfcgTPAaqAYOAikAe7A8z0u8XWtdUX76/4TeFVrfWSMmi/EmBjhz8lOrfWmsWq7kPXw7bTWnyilpvY4vBQ4p7XOB1BKbQNu01o/D6zveQ1lbNXzAvCOFHsxHo3E50Q4jgzpXN5koKjL4+L2Y/35LrAK2KSUenA0GyaEExnU50QpFaaU+j2wi5Ju3AAAALpJREFUUCn15Gg3TnSSHv4I0lr/Gvi1o9shhDPTWldh5FxijEkP//JKgCldHse1HxNCdJLPiYuQgn95B4FZSqlpSikvYDPwloPbJISzkc+Ji5CC304plQnsA+YopYqVUt/QWrcBDwHvASeBHVrrPEe2UwhHks+Ja5NpmUIIMUFID18IISYIKfhCCDFBSMEXQogJQgq+EEJMEFLwhRBigpCCL4QQE4QUfCGEmCCk4AshxAQhBV8IISaI/x/2Z8d0EHSjgAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.loglog(hs, errsH0, \"o-\", label=\"$H^0$ error\")\n",
    "plt.loglog(hs, errsH1, \"o-\", label=\"$H^1$ error\")\n",
    "plt.loglog(hs, hs**1, \"o-\", label=\"$h^1$\")\n",
    "plt.loglog(hs, hs**2, \"o-\", label=\"$h^2$\")\n",
    "plt.legend()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- Now play with $\\omega$ and the element order."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}