{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Bad Discretizations for 2D Stokes\n",
    "\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>\n",
    "----\n",
    "\n",
    "Follows [Braess](https://doi.org/10.1017/CBO9780511618635), Section III.7.\n",
    "\n",
    "**Of note:** This notebook contains a recipe for how to compute negative Sobolev norms, in one of the folds in part II, below.\n",
    "\n",
    "(Thanks to [Colin Cotter](https://www.imperial.ac.uk/people/colin.cotter) and [Matt Knepley](https://cse.buffalo.edu/~knepley/) for tips!)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import numpy.linalg as la\n",
    "import firedrake.mesh as fd_mesh\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "from firedrake import *"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 44,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x7fa71527ab90>"
      ]
     },
     "execution_count": 44,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAP4AAAD4CAYAAADMz1tMAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+j8jraAAAR1klEQVR4nO3df4xV5Z3H8fengJKNY60w7hIGHLC4YWwt6q0/qoug7QZoA9ngWkjdbhuU1pbuNpjdtWljXU22Wxvd2EirbLbb1rRQ6ma3kzqsSUEloaV1rNYfsJgphXKpKyMqjm1HZtzv/nEv9joMzJm5Z+4Pns8rmeSecx+e5zsXPvf84JzzKCIws7S8rd4FmFntOfhmCXLwzRLk4JslyME3S9DEeg08derUaG9vr9fwZkl4/PHHX4yI1qHr6xb89vZ2uru76zW8WRIk7RtuvXf1zRLk4JslyME3S1DdjvHNmsHAwADFYpH+/v56l3JCkydPpq2tjUmTJmVq7+CbnUCxWKSlpYX29nYk1bucYUUEhw4dolgsMmvWrEx/ZsRdfUnfkHRQ0jPHeV+SviqpR9JTki4cZd1mDau/v58pU6Y0bOgBJDFlypRR7ZVkOcb/JrDoBO8vBuaUf1YDX888ulkTaOTQHzXaGkfc1Y+IbZLaT9BkGfDtKN3fu0PSGZKmRcTzo6rkOLbcd20e3ZiNyZkXrOHV3lPqXcYxTm89p6o/n8cx/nRgf8VysbzumOBLWk1pr4CZM2eO2PE1m+/jPXt2ceTXL+RQZlYCavmMglqPV48xm3e8JXOvp6+vL5e+xuqmf7iVH23dxtQpZ7Llvx9g4sRTeXHP/zJ79p+Muc+antyLiPXAeoBCoZDpb+arvJveDU+Pa12VJFHLh5PUerx6jNnM4+3atYvps+eesE13dzeFQiGX8Ybzqc+s5e8/dwsf/ehHmT57HgcPvkZU+WWUx//jHwBmVCy3ldeZWQ7mz5/PmWeemWufeWzxO4E1kjYClwCH8zq+N2s012y+75h1fX19tPQ+Pqb+Hlj8iWpLGpMRgy9pA7AAmCqpCHwRmAQQEfcCXcASoAf4HfDx8SrWzPKR5az+yhHeD+DTuVVk1sCG20KP9zH+ePC1+mYJcvDNGtzKlSu57LLL2L17N21tbXznO9+quk9fq2/W4DZs2PCW5YMHX6v62gJv8c0S5OCbJcjBN0uQg2+WIAffLEEOvlmCHHyzBrd//34WLlxIR0cH5513HuvXf63qPv3/+GYNbuLEidx5551ceOGF9PX1MW/eBRQKl3DOOdPG3meO9ZnZOJg2bRrTppVC3tLSwpw5f8oLL1R3A6yDbzYK997z42PW9fX10b3j2PVZfHLN+0bVfu/evTzzzFO85z3VPdPWx/hmTeK1115j+fLl3H77P9PS0lJVX97im43CcFvoWtyWOzAwwPLly/nIRz7CBz+4zNfqm53sIoJVq1Yxd+5c1q5dm0ufDr5Zg9u+fTv3338/W7duZd68eVx11ft45JEtVfXpXX2zBnfFFVe85anBvi3XzMbEwTdLkINvliAH3yxBTXFyr9azlZ7s49VjzGYdb/Pmzfz2t78dsV13d3cu42XxjndMp9rfrimC36zzrjXiePUYs5nH27VrF3Pn1nfuvKHyOKvfFME3S1l/fz/z58/n9ddfZ3BwkMWLl3LjjX9TVZ8OvlmDO/XUU9m6dSunnXYaAwMDXHLJ+7j00ss555wPjblPn9wza3CSOO2004DSNfuDgwNVn8PwFt9sFO7c8uFj1vX19fHo4bHdLXfT1d/L1O6NN97goosuoqenh49//AbmzfNtuWYnvQkTJvDkk09SLBb5+c8f57nn/qeq/rzFNxuF4bbQtTyrf8YZZ3DFFfPZtu1hFi9eOOZ+vMU3a3C9vb288sorAPz+97/n0Ue3Mnv2O6vqM1PwJS2StFtSj6Sbh3l/pqSHJT0h6SlJS6qqysze9Pzzz7Nw4ULOP/983vve93LllVdx1VUfqKrPEXf1JU0A1gEfAIrAY5I6I2JnRbMvAJsi4uuSOoAuoL2qyswMgPPPP58nnnjizeVa3ZZ7MdATEXsi4giwEVg2pE0Ap5dfvx34TVVVmdm4yhL86cD+iuVieV2lW4HrJBUpbe0/M1xHklZL6pbU3dvbO4ZyzSwPeZ3cWwl8MyLagCXA/ZKO6Tsi1kdEISIKra2tOQ1tZqOVJfgHgBkVy23ldZVWAZsAIuInwGRgah4Fmln+sgT/MWCOpFmSTgFWAJ1D2vwauBpA0lxKwfe+vFmDGjH4ETEIrAEeAnZROnv/rKTbJC0tN7sJuEHSL4ANwMei1veamllmma7ci4guSiftKtfdUvF6J3B5vqWZWaU33niDQqHA1Kl/zL33/ntVffnKPbMmcffdd4/4UJCsHHyzJlAsFnnwwQe5/vrrc+nPN+mYjcKW+649Zl1fXx+HHx/bbblXf2JTpnaf/exnueOOO6q+Yu8ob/HNGtwPf/hDzjrrLC666KLc+vQW32wUhttCj/dtudu3b6ezs5Ouri76+/t59dVXWbt2DT/4wX+MuU9v8c0a3Je+9CWKxSJ79+5l48aNXH75fO66656q+nTwzRLkXX2zJrJgwQI6OgqeLdfMRs/BN0uQg282gma47WS0NTbFMX6zTrjYqOPVY8xmHe/uu+/m8OHDTJx44qjUc9LMiODQoUNMnjw5cx9NEfxmnXCxEcerx5jNPN7AwADFYpH+/v7jttm3bx9nn312LuNlcejQC/T39/P6kVfeXDd58mTa2toy99EUwTerl0mTJjFr1qwTtuno6KjpF9u99/yYRx59lI3f/9yY+/AxvlmCHHyzBDn4Zgly8M0S5OCbJcjBN0uQg2+WIAffLEEOvlmCHHyzBDn4Zgly8M0S5OCbJcjBN0uQg2+WIAffLEGZgi9pkaTdknok3XycNtdK2inpWUnfzbdMM8vTiE/gkTQBWAd8ACgCj0nqjIidFW3mAJ8DLo+IlyWdNV4Fm1n1smzxLwZ6ImJPRBwBNgLLhrS5AVgXES8DRMTBfMs0szxlCf50YH/FcrG8rtK5wLmStkvaIWnRcB1JWi2pW1J3b2/v2Co2s6rldXJvIjAHWACsBP5V0hlDG0XE+ogoREShtbU1p6HNbLSyBP8AMKNiua28rlIR6IyIgYj4FfAcpS8CM2tAWYL/GDBH0ixJpwArgM4hbf6L0tYeSVMp7frvybFOM8vRiMGPiEFgDfAQsAvYFBHPSrpN0tJys4eAQ5J2Ag8DfxcRh8araDOrTqYJNSKiC+gasu6WitcBrC3/mFmDa4qZdJp13rVGHa8eY3q8/Hz4mn+i2tGaIvjNOu9aI45XjzE9Xr6OTqFVDV+rb5YgB98sQQ6+WYIcfLMEOfhmCXLwzRLk4JslyME3S5CDb5YgB98sQQ6+WYIcfLMEOfhmCXLwzRLk4JslyME3S5CDb5YgB98sQQ6+WYIcfLMEOfhmCXLwzRLk4JslyME3S5CDb5YgB98sQU0xhdbJPA9aPcarx5geLz+eO28cnOzzrtVjTI+XL8+dZ2Zjkin4khZJ2i2pR9LNJ2i3XFJIKuRXopnlbcTgS5oArAMWAx3ASkkdw7RrAf4W+GneRZpZvrJs8S8GeiJiT0QcATYCy4ZpdzvwZaA/x/rMbBxkCf50YH/FcrG87k2SLgRmRMSDJ+pI0mpJ3ZK6e3t7R12smeWj6pN7kt4G3AXcNFLbiFgfEYWIKLS2tlY7tJmNUZbgHwBmVCy3ldcd1QK8C3hE0l7gUqDTJ/jMGleW4D8GzJE0S9IpwAqg8+ibEXE4IqZGRHtEtAM7gKUR0T0uFZtZ1UYMfkQMAmuAh4BdwKaIeFbSbZKWjneBZpa/TFfuRUQX0DVk3S3Habug+rLMbDz5yj2zBDn4Zgly8M0S5OCbJcjBN0uQg2+WIAffLEEOvlmCHHyzBDn4Zgly8M0S5OCbJcjBN0uQg2+WIAffLEEOvlmCmmIKrZN5HrR6jFePMT1efjx33jg42eddq8eYHi9fnjvPzMbEwTdLkINvliAH3yxBDr5Zghx8swQ5+GYJcvDNEuTgmyXIwTdLkINvliAH3yxBmYIvaZGk3ZJ6JN08zPtrJe2U9JSkLZLOzr9UM8vLiMGXNAFYBywGOoCVkjqGNHsCKETE+cADwB15F2pm+cmyxb8Y6ImIPRFxBNgILKtsEBEPR8Tvyos7gLZ8yzSzPGUJ/nRgf8VysbzueFYBm4d7Q9JqSd2Sunt7e7NXaWa5yvXknqTrgALwleHej4j1EVGIiEJra2ueQ5vZKGR5As8BYEbFclt53VtIej/weeDKiHg9n/LMbDxk2eI/BsyRNEvSKcAKoLOygaQLgPuApRFxMP8yzSxPIwY/IgaBNcBDwC5gU0Q8K+k2SUvLzb4CnAZ8X9KTkjqP052ZNYBMD9uMiC6ga8i6Wypevz/nusxsHPnKPbMEOfhmCXLwzRLk4JslyME3S5CDb5agppg772SeALEe49VjTI+XH0+aOQ5O9gkX6zGmx8uXJ800szFx8M0S5OCbJcjBN0uQg2+WIAffLEEOvlmCHHyzBDn4Zgly8M0S5OCbJcjBN0uQg2+WIAffLEEOvlmCHHyzBDn4Zgly8M0S5OCbJcjBN0uQg2+WIAffLEEOvlmCMgVf0iJJuyX1SLp5mPdPlfS98vs/ldSed6Fmlp8Rgy9pArAOWAx0ACsldQxptgp4OSLeCfwL8OW8CzWz/GSZSedioCci9gBI2ggsA3ZWtFkG3Fp+/QBwjyRFDtOL/Nk7n+bOLR+utpvMPvSFs07q8eoxpsfL129OPwxcWVUfGimbkq4BFkXE9eXlvwIuiYg1FW2eKbcplpd/WW7z4pC+VgOrAWbOnHnRvn37RizwL/5xOoODg6P6pcxOdtvXDfDSSy+N2E7S4xFRGLq+pnPnRcR6YD1AoVDItDfwn188MK41mTWl26v741lO7h0AZlQst5XXDdtG0kTg7cCh6kozs/GSJfiPAXMkzZJ0CrAC6BzSphP46/Lra4CteRzfm9n4GHFXPyIGJa0BHgImAN+IiGcl3QZ0R0Qn8G/A/ZJ6gJcofTmYWYPKdIwfEV1A15B1t1S87gf+Mt/SzGy8+Mo9swQ5+GYJcvDNEuTgmyVoxCv3xm1gqRcY+dI9mAq8OGKr+mjk2sD1VaORa4Ps9Z0dEa1DV9Yt+FlJ6h7uksNG0Mi1geurRiPXBtXX5119swQ5+GYJaobgr693ASfQyLWB66tGI9cGVdbX8Mf4Zpa/Ztjim1nOHHyzBDVM8Bv5gZ4ZalsraaekpyRtkXR2rWrLUl9Fu+WSQlLN/psqS22Sri1/fs9K+m6tastSn6SZkh6W9ET573dJDWv7hqSD5SdcDfe+JH21XPtTki7M3HlE1P2H0u2+vwRmA6cAvwA6hrT5FHBv+fUK4HsNVNtC4I/Kr2+sVW1Z6yu3awG2ATuAQqPUBswBngDeUV4+q5E+O0on0W4sv+4A9tawvvnAhcAzx3l/CbAZEHAp8NOsfTfKFv/NB3pGxBHg6AM9Ky0DvlV+/QBwtSQ1Qm0R8XBE/K68uIPSU4pqJctnB6WHNX0Z6G+w2m4A1kXEywARcbDB6gvg9PLrtwO/qVVxEbGN0vMtjmcZ8O0o2QGcIWlalr4bJfjTgf0Vy8XyumHbRMQgcBiY0iC1VVpF6Vu4Vkasr7wLOCMiHqxhXZDtszsXOFfSdkk7JC2qWXXZ6rsVuE5SkdIzKT5Tm9IyGe2/zTfV9GGbJztJ1wEFqn32cY4kvQ24C/hYnUs5nomUdvcXUNpT2ibp3RHxSl2r+oOVwDcj4k5Jl1F60tS7IuL/6l1YNRpli9/ID/TMUhuS3g98HlgaEa/XoK6jRqqvBXgX8IikvZSOBTtrdIIvy2dXBDojYiAifgU8R+mLoBay1LcK2AQQET8BJlO6QaYRZPq3OaxanagY4STGRGAPMIs/nGQ5b0ibT/PWk3ubGqi2CyidJJrTiJ/dkPaPULuTe1k+u0XAt8qvp1LadZ3SQPVtBj5Wfj2X0jG+avj3287xT+59kLee3PtZ5n5r9Qtk+AWXUPq2/yXw+fK62yhtQaH0Tft9oAf4GTC7gWr7EfAC8GT5p7ORPrshbWsW/IyfnSgdiuwEngZWNNJnR+lM/vbyl8KTwJ/XsLYNwPPAAKU9o1XAJ4FPVnx268q1Pz2av1dfsmuWoEY5xjezGnLwzRLk4JslyME3S5CDb5YgB98sQQ6+WYL+H/8gjH52kmbUAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "nx = 7\n",
    "mesh = UnitSquareMesh(nx, nx, quadrilateral=True)\n",
    "\n",
    "triplot(mesh)\n",
    "plt.gca().set_aspect(\"equal\")\n",
    "plt.legend()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Part I: The Checkerboard Instability\n",
    "\n",
    "$\\let\\b=\\boldsymbol$Build a $Q^1$-$Q^0$ mixed space:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 45,
   "metadata": {},
   "outputs": [],
   "source": [
    "V = VectorFunctionSpace(mesh, \"CG\", 1)\n",
    "W = FunctionSpace(mesh, \"DG\", 0)\n",
    "\n",
    "Z = V * W"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 46,
   "metadata": {},
   "outputs": [],
   "source": [
    "x = SpatialCoordinate(mesh)\n",
    "x_w = interpolate(x[0], W)\n",
    "y_w = interpolate(x[1], W)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 47,
   "metadata": {},
   "outputs": [],
   "source": [
    "x_array = x_w.dat.data.copy()\n",
    "y_array = y_w.dat.data.copy()\n",
    "\n",
    "x_idx = (np.round(x_array * nx * 2) - 1)//2\n",
    "y_idx = (np.round(y_array * nx * 2) - 1)//2\n",
    "                  \n",
    "checkerboard = (x_idx+y_idx) % 2 - 0.5\n",
    "\n",
    "q = Function(W, checkerboard)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 48,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.colorbar.Colorbar at 0x7fa6dce34590>"
      ]
     },
     "execution_count": 48,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXIAAAD8CAYAAABq6S8VAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+j8jraAAAaRUlEQVR4nO3dfZBddZ3n8fcn3dPdsqA2iTZZAibUBMeIFg+9gOXiZCFIZLcSa3UsGB3jFppiFFeK0TKUuywV1yoYd7TW2vjQ4zCChTyIu9i1ZGRChGFXBdMzRiChkBAZ6BhBIIIuktDw2T/uCVw7t/uem77d957weVXd4jz8zjnfc7r59MnvngfZJiIiqmtepwuIiIiZSZBHRFRcgjwiouIS5BERFZcgj4iouAR5RETFNQ1ySVdJelzSfVPMl6QvSdoh6R5JJ7e/zIiI7iFppaQHitxbN0Wb90naLmmbpG/VTV8j6cHis6Yt9TS7jlzSO4DfAtfYPqHB/HOBjwPnAqcB/932ae0oLiKi20jqAX4GnA2MA1uA821vr2uzFLgRONP2Hkmvt/24pCOBMWAYMPCPwCm298ykpqZn5LbvBJ6apslqaiFv23cBr5W0cCZFRUR0sVOBHbZ32t4HXE8tB+t9BNiwP6BtP15MPwfYZPupYt4mYOVMC+qd6QqAo4FH68bHi2m7JzeUtBZYC3DYq3pOOfJVfW3YfHv1HzbR6RKmtPfZdvy42q9bj1mOV+u68ZjN6+3hkcd/+4Tt181kPWcsH/Cep14s1Xbbvc9vA56rmzRie6QYbpR5k3shjgeQ9AOgB7jc9vemWPbosvswlTn9qRUHYgRgYOFRPsH/fi43X8ql/+smhoa675cZ4MKTVjDQe3inyzhAtx6zHK/WdeMxW3LmaXzlug/980zXs+epF/nOLQtKtf2jY3c/Z3t4BpvrBZYCy4FFwJ2S3jKD9U2rHVet7AKOqRtfVEyLiDgUlcm8cWDU9vO2f06tT31pyWVb1o4gHwU+WFy9cjrwtO0DulUiIg4RW4ClkpZI6gPOo5aD9W6mdjaOpAXUulp2ArcC75Q0KGkQeGcxbUaa/vtO0nVFQQskjQP/BfgDANtfBTZSu2JlB/As8B9mWlRERLeyPSHpImoB3ANcZXubpPXAmO1RXg7s7cALwKdsPwkg6bPU/hgArLc93cUkpTQNctvnN5lv4GMzLSQioipsb6R2Els/7bK6YQOXFJ/Jy14FXNXOenJnZ0RExSXIIyIqLkEeEVFxCfKIiIpLkEdEVFyCPCKi4hLkEREVlyCPiKi4BHlERMUlyCMiKi5BHhFRcQnyiIiKS5BHRFRcgjwiouIS5BERFZcgj4iouAR5RETFJcgjIiouQR4R0SJJKyU9IGmHpHXTtHuPJEsaLsYXS/qdpK3F56vtqKfpOzsjIuJlknqADcDZwDiwRdKo7e2T2h0BfAK4e9IqHrJ9Yjtryhl5RERrTgV22N5pex9wPbC6QbvPAlcCz812QQnyiIjWHA08Wjc+Xkx7iaSTgWNs39Jg+SWSfiLpHySd0Y6C0rUSEa8Ie144jJueOblk61sWSBqrmzBie6TMkpLmAV8APtRg9m7gWNtPSjoFuFnSm20/U7KwhhLkEREHesL28BTzdgHH1I0vKqbtdwRwAnCHJICjgFFJq2yPAXsBbP+jpIeA44H6PxotS9dKRERrtgBLJS2R1AecB4zun2n7adsLbC+2vRi4C1hle0zS64ovS5F0HLAU2DnTgjp2Rt532ASX3nBTpzY/pXXnLKdfA50uo6GLr93I4Pzu+9vbrccsx6t13XjMfrjpOLiu01W8zPaEpIuAW4Ee4Crb2yStB8Zsj06z+DuA9ZKeB14ELrT91Exr6liQCxga6r6enX4NMNB7eKfLaGhw/rwcsxbkeLWuG49ZX586XcIBbG8ENk6adtkUbZfXDX8H+E676+muP70REdGyBHlERMUlyCMiKi5BHhFRcQnyiIiKS5BHRFRcgjwiouIS5BERFVcqyJs9RF3SsZJuL57odY+kc9tfakRENNI0yOseov4uYBlwvqRlk5r9J+BG2ydRe+7Al9tdaERENFbmjLzMQ9QNvLoYfg3wi/aVGBER0ykT5E0fog5cDnxA0ji15w98vNGKJK2VNCZpbGLCB1FuRERM1q4vO88HvmF7EXAu8M3i4eq/x/aI7WHbw7293fcgnIiIKioT5M0eog5wAXAjgO0fAQPAgnYUGBER0ysT5NM+RL3wCHAWgKQ3UQvyX7Wz0IiIaKxpkNueAPY/RP1+alenbJO0XtKqotlfAB+R9FNqj4D/kO10gkdEzIFST5Bv9hB129uBt7e3tIiIKCN3dkZEVFyCPCKi4hLkEREtKvHYkgsl3Stpq6T/W383vKRLi+UekHROO+pJkEdEtKDkY0u+Zfsttk8E/hL4QrHsMmpX/r0ZWAl8uVjfjCTIIyJa0/SxJbafqRv9F9QeY0LR7nrbe23/HNhRrG9GSl21EhFRdc9MDLDpl39UsvUtCySN1U0YsT1SDDd6bMlpk9cg6WPAJUAfcGbdsndNWnbyI09aliCPiDjQE7aHZ7IC2xuADZL+lNoTYte0pbIG0rUSEdGaMo8tqXc98O6DXLaUBHlERGuaPrZE0tK60X8LPFgMjwLnSeqXtARYCvx4pgWlayUiogW2JyTtf2xJD3DV/seWAGO2R4GLJK0Angf2UHSrFO1uBLYDE8DHbL8w05oS5BERLSrx2JJPTLPs54DPtbOedK1ERFRcx87I9z7bw4UnrejU5qd08bUbGZzfnX/f1p2znH4NdLqMA3TrMcvxal03HrM3nj3Y6RK6Xge7VsRA7+Gd2/wUBufPY2ioO3uc+jWQY9aCHK/WdeMxmzfzGx8Ped15WhAREaUlyCMiKi5BHhFRcQnyiIiKS5BHRFRcgjwiouIS5BERFZcgj4iouAR5RETFJcgjIiouQR4RUXEJ8oiIikuQR0RUXII8IqLiEuQRERWXII+IaJGklZIekLRD0roG898h6Z8kTUh676R5L0jaWnxGJy97MLrz6fYREV1KUg+wATgbGAe2SBq1vb2u2SPAh4BPNljF72yf2M6aEuQREa05FdhheyeApOuB1cBLQW774WLei3NRUII8Il4R9u3r5eHx15VtvkDSWN34iO2RYvho4NG6eePAaS2UMlCsewK4wvbNLSzbUII8IuJAT9genqV1v8H2LknHAd+XdK/th2aywlJfdjbr2C/avE/SdknbJH1rJkVFRHSxXcAxdeOLimml2N5V/HcncAdw0kwLahrkdR377wKWAedLWjapzVLgUuDttt8MXDzTwiIiutQWYKmkJZL6gPOAUlefSBqU1F8MLwDeTl3f+sEqc0b+Use+7X3A/o79eh8BNtjeA2D78ZkWFhHRjWxPABcBtwL3Azfa3iZpvaRVAJL+laRx4E+Ar0naViz+JmBM0k+B26n1kc84yMv0kZfp2D8eQNIPgB7gctvfm7wiSWuBtQB9R5X+0iEioqvY3ghsnDTtsrrhLdS6XCYv90PgLe2up11fdvYCS4Hl1Iq/U9JbbP+6vlHxre8IwMDCo9ymbUdEvKKV6Vop07E/Dozaft72z4GfUQv2iIiYZWWCvEzH/s3Uzsb3d+AfD+xsY50RETGFpkFepmO/mPekpO3UOvA/ZfvJ2So6IiJeVqqPvETHvoFLik9ERMyhPP0wIqLiOnaLft9hE1x6w02d2vyU1p2znH4NdLqMhi6+diOD87vvb2+3HrMcr9Z14zH74abj4LpOV9HdOhbkAoaGuu9RL/0aYKD38E6X0dDg/Hk5Zi3I8WpdNx6zvj51uoSu111/eiMiomUJ8oiIikuQR0RUXII8IqLiEuQRERWXII+IqLgEeURExSXIIyIqLkEeEVFxCfKIiIpLkEdEtEjSSkkPSNohaV2D+f2Sbijm3y1pcd28S4vpD0g6px31JMgjIlogqQfYALwLWAacL2nZpGYXAHts/yHwReDKYtll1F7O82ZgJfDlYn0zkiCPiGjNqcAO2ztt7wOuB1ZParMauLoYvgk4S5KK6dfb3lu8FnNHsb4ZSZBHRBxogaSxus/aunlHA4/WjY8X02jUpnjL2tPA/JLLtqy7nlcZETFLtE/0P9JXtvkTtodns552yhl5RERrdgHH1I0vKqY1bCOpF3gN8GTJZVuWII+IaM0WYKmkJZL6qH15OTqpzSiwphh+L/D94t3Go8B5xVUtS4ClwI9nWlC6ViIiWmB7QtJFwK1AD3CV7W2S1gNjtkeBvwG+KWkH8BS1sKdodyOwHZgAPmb7hZnWlCCPiGiR7Y3AxknTLqsbfg74kymW/RzwuXbWk66ViIiKS5BHRFRcgjwiouIS5BERFZcgj4iouAR5RETFJcgjIiouQR4RUXEJ8oiIikuQR0RUXMdu0d/7bA8XnrSiU5uf0sXXbmRwfnf+fVt3znL6NdDpMg7Qrccsx6t13XjM3nj2YKdL6HodfNaKGOg9vHObn8Lg/HkMDXXnI2j6NZBj1oIcr9Z14zGbN/M3oR3yuvO0ICIiSisV5M3eGF3X7j2SLKkyb9aIiKi6pkFe8o3RSDoC+ARwd7uLjIiIqZU5Iy/zxmiAzwJXAs+1sb6IiGiiTJA3feuzpJOBY2zfMt2KJK3d/1Zq82LLxUZExIFm/GWnpHnAF4C/aNbW9ojtYdvDyvesERFtUSZNm731+QjgBOAOSQ8DpwOj+cIzIl6JJB0paZOkB4v/NrwQXtL3JP1a0v+eNP0bkn4uaWvxObHZNssE+bRvjLb9tO0FthfbXgzcBayyPVZi3RERh5p1wGbbS4HNxXgjnwf+bIp5n7J9YvHZ2myDTYPc9gSw/43R9wM37n9jtKRVzZaPiHiFWQ1cXQxfDby7USPbm4HftGODpW4va/bG6EnTl8+8rIiIjlogqb5XYcT2SMllh2zvLoZ/CQwdxPY/J+kyijN623una9yd9wlHRLRZzz444p9dtvkTtqf8nk/SbcBRDWZ9pn7EtiWV3mjhUmp/APqAEeDTwPrpFkiQR0S0yPaUT/yT9JikhbZ3S1oIPN7iuvefze+V9LfAJ5stk2sAIyLaaxRYUwyvAb7bysJF+CNJ1PrX72u2TII8IqK9rgDOlvQgsKIYR9KwpK/vbyTp/wDfBs6SNC7pnGLWtZLuBe4FFgD/tdkG07USEdFGtp8EzmowfQz4cN34GVMsf2ar28wZeURExSXIIyIqLkEeEVFxCfKIiIpLkEdEVFyCPCKi4hLkEREVlyCPiKi4BHlERMUlyCMiKi5BHhFRcQnyiIiKS5BHRFRcgjwiouI69hjbvsMmuPSGmzq1+SmtO2c5/RrodBkNXXztRgbnd9/f3m49ZjlerevGY/bDTcfBdZ2uort1LMgFDA113+PQ+zXAQO/hnS6jocH583LMWpDj1bpuPGZ9fep0CV2vu/70RkREyxLkEREVlyCPiGgjSUdK2iTpweK/gw3avEHSP0naKmmbpAvr5p0i6V5JOyR9qXgJ87QS5BER7bUO2Gx7KbC5GJ9sN/A22ycCpwHrJP3LYt5XgI8AS4vPymYbTJBHRLTXauDqYvhq4N2TG9jeZ3tvMdpPkcWSFgKvtn2XbQPXNFp+sgR5RMSBFkgaq/usbWHZIdu7i+FfAkONGkk6RtI9wKPAlbZ/ARwNjNc1Gy+mTau7rjOKiJglPc+Z1zy0t3nDmidsD081U9JtwFENZn2mfsS2JbnROmw/Cry16FK5WdJB31iTII+IaJHtFVPNk/SYpIW2dxddJY83WdcvJN0HnAH8AFhUN3sRsKtZPelaiYhor1FgTTG8Bvju5AaSFkl6VTE8CPxr4IGiS+YZSacXV6t8sNHykyXIIyLa6wrgbEkPAiuKcSQNS/p60eZNwN2Sfgr8A/DfbN9bzPso8HVgB/AQ8HfNNpiulYiINrL9JHBWg+ljwIeL4U3AW6dYfgw4oZVt5ow8IqLiEuQRERVXKsglrZT0QHHL6AF3KUm6RNJ2SfdI2izpDe0vNSIiGmka5JJ6gA3Au4BlwPmSlk1q9hNg2PZbgZuAv2x3oRER0ViZM/JTgR22d9reB1xP7RbUl9i+3fazxehd/P51kBERMYvKBPnR1G4h3a/ZLaMXMMXlMpLW7r/ldWKi4c1OERHRorZ+2SnpA8Aw8PlG822P2B62Pdzbm7d+RES0Q5nryHcBx9SNN7xlVNIKas8Z+OO6p3pFRMQsK3NGvgVYKmmJpD7gPGq3oL5E0knA14BVtqd9rkBERLRX0yC3PQFcBNwK3A/caHubpPWSVhXNPg8cDny7eOPF6BSri4iINit1i77tjcDGSdMuqxue8klgERExu3JnZ0RExSXIIyIqLkEeEVFxCfKIiIpLkEdEVFyCPCKi4hLkEREV17FXve19tocLT+q+y88vvnYjg/O78+/bunOW06+BTpdxgG49ZjlerevGY/bGswc7XUJLJB0J3AAsBh4G3md7zxRtXw1sB262fVEx7Q5gIfC7otk7m90x38F3doqB3sM7t/kpDM6fx9BQd77KtF8DOWYtyPFqXTces3nq6XQJrVoHbLZ9RfEinnXAp6do+1ngzgbT31+8u7OU7jwtiIiortXA1cXw1cC7GzWSdAowBPz9TDeYII+IaK8h27uL4V9SC+vfI2ke8FfAJ6dYx98Wz636z5KaPvO7O/99FxHRZnpuH333P9q8Yc0CSfVdGyO2R15al3QbcFSD5T5TP2Lbkhq9ReejwEbb4w1y+v22d0k6AvgO8GfANdMVmyCPiDjQE7aHp5o53YMCJT0maaHt3ZIWAo2+qHwbcIakj1J7cmyfpN/aXmd7V7GN30j6FrXXbU4b5OlaiYhor1FgTTG8Bvju5Aa232/7WNuLqXWvXGN7naReSQsAJP0B8O+A+5ptMEEeEdFeVwBnS3oQWFGMI2lY0tebLNsP3CrpHmArtbex/XWzDaZrJSKijWw/CZzVYPoY8OEG078BfKMY/n/AKa1uM2fkEREVlyCPiKi4BHlERMUlyCMiKi5BHhFRcQnyiIiKS5BHRFRcgjwiouIS5BERFZcgj4iouAR5RETFJcgjIiouQR4RUXEJ8oiIikuQR0RUXII8IqLiEuQRERWXII+IqLgEeURExZUKckkrJT0gaYekdQ3m90u6oZh/t6TF7S40IqIKJB0paZOkB4v/Dk7R7lhJfy/pfknb9+empCVFju4ocrWv2TabBrmkHmAD8C5gGXC+pGWTml0A7LH9h8AXgSubrTci4hC1DthseymwuRhv5Brg87bfBJwKPF5MvxL4YpGne6jl67R6SxR1KrDD9k4ASdcDq4HtdW1WA5cXwzcB/0OSbHuqlR7W/xtef8GPSmx+bm0a76PvMXW6jIYW/cetzFNPp8s4QLcesxyv1nXjMXvhtdvguk5X0ZLVwPJi+GrgDuDT9Q2Kk+Fe25sAbP+2mC7gTOBP65a/HPjKdBssE+RHA4/WjY8Dp03VxvaEpKeB+cATk4pfC6wtRvdec9nW+0psv2oWMGm/22d8dlbb3Czu02xqerwqul/TmuE+dex3rJk3znQFz0z86tbvPfaVBSWbD0gaqxsfsT1Sctkh27uL4V8CQw3aHA/8WtL/BJYAt1E7cx8Efm17omg3Ti1fp1UmyNumOBAjAJLGbA/P5fbnwqG4X4fiPsGhuV+H4j5Bbb9mug7bK9tRC4Ck24CjGsz6zKRtWlKjnole4AzgJOAR4AbgQ8B3D6aeMkG+CzimbnxRMa1Rm3FJvcBrgCcPpqCIiG5ne8VU8yQ9Jmmh7d2SFvJy33e9cWBrXZf1zcDpwFXAayX1FmfljfL2AGWuWtkCLC2+Se0DzgNGJ7UZBdYUw+8Fvj9d/3hExCGsPg/X0Pgsewu1wH5dMX4msL3Izdup5eh0y/+epkFe/FW4CLgVuB+40fY2SeslrSqa/Q0wX9IO4BKm/pa2Xtn+pqo5FPfrUNwnODT361DcJ6jWfl0BnC3pQWBFMY6kYUlfB7D9AvBJYLOkewEBf10s/2ngkiJP51PL12kpJ84REdWWOzsjIiouQR4RUXGzHuSH4u39JfbpkuKW23skbZb0hk7U2apm+1XX7j2SLKnrL3Mrs0+S3lf8vLZJ+tZc13gwSvwOHivpdkk/KX4Pz+1Ena2QdJWkxyU1vL9ENV8q9vkeSSfPdY1dy/asfYAe4CHgOKAP+CmwbFKbjwJfLYbPA26YzZrmaJ/+DXBYMfzn3b5PZferaHcEcCdwFzDc6brb8LNaCvwEGCzGX9/putu0XyPAnxfDy4CHO113if16B3AycN8U888F/o7aF4OnA3d3uuZu+cz2GflLt/fb3gfsv72/3mpqt6FC7fb+s4rbVLtV032yfbvtZ4vRu6hdC9rtyvysAD5L7VkQz81lcQepzD59BNhgew+A7UbX/HabMvtl4NXF8GuAX8xhfQfF9p3AU9M0WQ1c45q7qF2+t3Buqutusx3kjW7vn3y76e/d3g/sv72/W5XZp3oXUDuL6HZN96v4p+wxtm+Zy8JmoMzP6njgeEk/kHSXpLbd/TeLyuzX5cAHJI0DG4GPz01ps6rV//deMeb0Fv1XGkkfAIaBP+50LTMlaR7wBWq3ER9Keql1ryyn9i+nOyW9xfavO1rVzJ0PfMP2X0l6G/BNSSfYfrHThUX7zfYZeSu391OR2/vL7BOSVlB77sIq23vnqLaZaLZfRwAnAHdIephaH+Vol3/hWeZnNQ6M2n7e9s+Bn1EL9m5WZr8uAG4EsP0jYIDaA7WqrNT/e69Esx3kh+Lt/U33SdJJwNeohXgV+lyhyX7Zftr2AtuLbS+m1ve/yvaMH2Y0i8r8/t1M8chRSQuodbXsnMsiD0KZ/XoEOAtA0puoBfmv5rTK9hsFPlhcvXI68LRffsrgK9tsf5tK7Zvmn1H7lv0zxbT11EIAar9g3wZ2AD8Gjuv0N8Bt2KfbgMeArcVntNM1t2O/JrW9gy6/aqXkz0rUuoy2A/cC53W65jbt1zLgB9SuaNkKvLPTNZfYp+uA3cDz1P6ldAFwIXBh3c9qQ7HP91bh92+uPrlFPyKi4nJnZ0RExSXIIyIqLkEeEVFxCfKIiIpLkEdEVFyCPCKi4hLkEREV9/8BHGGFWo1DBBMAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "ax = plt.gca()\n",
    "l = tricontourf(q, axes=ax)\n",
    "triplot(mesh, axes=ax, interior_kw=dict(alpha=0.05))\n",
    "plt.colorbar(l)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Assemble the discrete coefficients representing $\\int (\\nabla \\cdot \\boldsymbol u) q$:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 54,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,\n",
       "        -0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0., -0.,  0.,  0.,\n",
       "         0.,  0.,  0., -0.,  0., -0.,  0.,  0.,  0.,  0.,  0.,  0., -0.,\n",
       "         0.,  0.,  0.,  0.,  0.,  0.,  0., -0.,  0.,  0.,  0., -0.,  0.,\n",
       "         0.,  0.,  0., -0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.],\n",
       "       [ 0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,\n",
       "         0.,  0.,  0.,  0.,  0.,  0.,  0.,  0., -0.,  0.,  0.,  0.,  0.,\n",
       "         0.,  0., -0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0., -0.,\n",
       "         0.,  0.,  0.,  0.,  0.,  0., -0.,  0.,  0.,  0.,  0.,  0.,  0.,\n",
       "        -0.,  0.,  0.,  0., -0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.]])"
      ]
     },
     "execution_count": 54,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "utest, ptest = TestFunctions(Z)\n",
    "bcs = [DirichletBC(Z.sub(0), Constant((0, 0)), (1, 2, 3, 4))]\n",
    "coeffs = assemble(div(utest)*q*dx, bcs=bcs)\n",
    "coeffs.dat.data[0].round(5).T"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Is this bad news?\n",
    "\n",
    "$$b (\\b{v}, q) = \\int_{\\Omega} \\nabla \\cdot \\b{v}q.$$\n",
    "\n",
    "Needed:\n",
    "\n",
    "> There exists a constant $c_2 > 0$ so that (*inf-sup* or *LBB condition*):\n",
    "> $$ \\inf_{\\mu \\in M} \\sup_{v \\in X} \\frac{b (v, \\mu)}{||v||_X ||\\mu||_M} \\geqslant c_2 .$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Part II: Is Removing the Checkerboard Sufficient?\n",
    "\n",
    "Suppose we consider the space that has the checkerboard projected out. Is that better?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[-3.  2.  3. -1. -2. -3. -0.  1.  2.  3.  1.  0. -1. -2. -3. -2. -1. -0.\n",
      "  1.  2.  3.  3.  2.  1.  0. -1. -2. -3. -3. -2. -1. -0.  1.  2.  3.  2.\n",
      "  1.  0. -1. -3. -2. -1. -0.  3.  2.  1. -3. -2.  3.]\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "<matplotlib.colorbar.Colorbar at 0x7fcdfb9e8f50>"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAWwAAAD8CAYAAABTjp5OAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+j8jraAAAZ3ElEQVR4nO3dfZRcdZ3n8fcnXV1dpJMOPWkJSFDgGNGWQYI9CMMOZATWkHWT1VE32cMKY8YsswvHOQIObFxkcXfPKGd84MiIrbIgywAxjtIrUQwQZFXCEIfHJMKEB0kgEB5CcMxA0uS7f9TNnEpXddftdHXfezuf1zl9qHvrV/f3/XVVf7j51X1QRGBmZvk3JesCzMwsHQe2mVlBOLDNzArCgW1mVhAObDOzgnBgm5kVRNPAlnStpG2SHh3meUm6StImSQ9LOqH1ZZqZ5YOkiqS/l/SQpPWS/nuDNp+RtCHJxDslvb0VfafZw74OmD/C82cBc5KfZcA3xl6WmVluvQF8ICLeCxwPzJd00pA2DwB9EXEcsBL4Uis6bhrYEXEP8MoITRYB342qtcDBkg5rRXFmZnmTZN0/JYvtyU8MabMmInYmi2uB2a3ou9SCbRwObK5Z3pKs2zq0oaRlVPfCmXpQ2/t+76ByC7pvrbaD9mRdwrAGd7dlXUJDkc+ykAJlXUQju3JZFQAdld1Zl1DnzcEyTz//u5ci4i1j2c4fzavE9lfS/X2vf2T3euD1mlX9EdG/d0FSG/Ar4B3A1RFx3wibWwr8ePQV12tFYKeWDLgfoHLYoXFsfGQiu0/luBsep+uQStZlNHTL/3gf5c4ZWZdR5/ljX6bU3ZV1GXUOrTxFpWda1mXU2f25diqV/P2+AD5/5VeYNWtCY6Gp21ZcyIUXXfSbsW5n+yt7+P5tPanavuttW1+PiL7hno+IN4HjJR0M/EDSsRFR9z2fpLOBPuC0/Sx7H604SuRZ4Iia5dnJOjOzSS0iXgXW0OB7PklnAMuBhRHxRiv6a0VgDwCfSI4WOQnYERF10yFmZpOBpLcke9ZIOgg4E/j1kDZzgW9SDettreq76b99JN0EzAN6JG0BPk91kp2IuAZYBSwANgE7gT9tVXFmZjl0GHB9Mo89BVgRET+SdAWwLiIGgCuBacD3JAE8ExELx9px08COiCVNng/gv4y1EDOzIoiIh4G5DdZfVvP4jPHo22c6mpkVhAPbzKwgHNhmZgXhwDYzKwgHtplZQTiwzcwKwoFtZlYQDmwzs4JwYJuZFYQD28ysIBzYZmYF4cA2MysIB7aZWUE4sM3MCsKBbWZWEA5sM7OCcGCbmRWEA9vMrCAc2GZmoyDpCElrJG2QtF7Sp0do+weSBiV9tBV9N72no5mZ7WMQuDAi/kHSdOBXklZHxIbaRslNer8I/LRVHXsP28xsFCJia0T8Q/L4t8BG4PAGTS8Avg9sa1XfDmwzs/0k6Uiqd1C/b8j6w4EPA99oZX+eEjGzA8L2N6ey8rUTUra+rUfSupoV/RHRX9tC0jSqe9B/ERGvDdnAV4G/jIg9kva/6CEc2GZm9V6KiL7hnpTUTjWsb4yIv2vQpA+4OQnrHmCBpMGI+OFYinJgm5mNgqop/B1gY0R8uVGbiDiqpv11wI/GGtaQYWCXpw5y6S0rs+p+WBedt5j29s6sy2io/MmnKHVNzbqMOoOPd2RdQkM7fjaNnZX8vZenXnwXnd35/J1d8sF5dKiSdRn7OObM7qxLGOoU4D8Cj0h6MFn3X4G3AUTENePVcWaBLWDWrPzt4Le3d1KpdGVdRkNvdm2n0jMt6zLqlLZ2UOrO3++svbKTcueMrMuo09ndQdch+QrFvTpUoVLK12dsitqyLmEfEfFzqhGWtv25rerbR4mYmRWEA9vMrCAc2GZmBeHANjMrCAe2mVlBOLDNzArCgW1mVhAObDOzgkgV2JLmS3pM0iZJlzR4/m3JBb0fkPSwpAWtL9XM7MDWNLCTi3BfDZwF9AJLJPUOafY5YEVEzAUWA3/T6kLNzA50afawTwQ2RcSTEbELuBlYNKRNAHvPTZ4BPNe6Es3MDNIF9uHA5prlLdTfXeFy4GxJW4BVVO+0UEfSMknrJK0bHIz9KNfM7MDVqi8dlwDXRcRsYAFwg6S6bUdEf0T0RURfqdS6i3qbmR0I0gT2s8ARNcuzk3W1lgIrACLiXqBC9aLdZmbWImkC+35gjqSjJJWpfqk4MKTNM8DpAJLeTTWwX2xloWZmB7qmgR0Rg8D5wO1U7w68IiLWS7pC0sKk2YXApyQ9BNwEnBsRnqQ2M2uhVHcQiIhVVL9MrF13Wc3jDVTvwmBmZuPEZzqamRWEA9vMrCAc2GZmoyTpWknbJD06Qpt5kh6UtF7Sz1rRrwPbzGz0rgPmD/ekpIOpXqJjYUS8B/hYKzp1YJuZjVJE3AO8MkKT/wD8XUQ8k7Tf1op+Ux0lYmZWdK8NVlj9/LtStr6tR9K6mhX9EdE/iu7eCbRLuhuYDnwtIr47itc35MA2M6v3UkT0jeH1JeB9VE8oPAi4V9LaiHh8LEU5sM3MWm8L8HJE/A74naR7gPcCYwpsz2GbmbXercC/klSSNBV4P9UzxcfEe9hmZqMk6SZgHtCTXFb680A7QERcExEbJf0EeBjYA3w7IoY9BDAtB7aZ2ShFxJIUba4Ermxlv54SMTMriMz2sN/Y2cZ5c8/Iqvthze3/OZ3dHVmX0dDAVcezs9KZdRl1BntfzbqEhsof2Eypa2rWZdRZc3Ev7e35ex8BPnvjt+iema/9uF+uPrp6DVDLckpEVErTsut+GJ3dHXQdUsm6jIbaK52UO2dkXUad0vRBSt1dzRtOsFJlKpWe/H3Gdre3U6nk7/cF0D1zCrNm5WumtFz23an2ytf/Ss3MbFgObDOzgnBgm5kVhAPbzKwgHNhmZgXhwDYzKwgHtplZQTiwzcwKwoFtZlYQDmwzs4JwYJuZFYQD28ysIBzYZmYF4cA2MysIB7aZWUE4sM3MRknSfEmPSdok6ZIGz79N0hpJD0h6WNKCVvTrwDYzGwVJbcDVwFlAL7BEUu+QZp8DVkTEXGAx8Det6NuBbWY2OicCmyLiyYjYBdwMLBrSJoC9txWaATzXio7zdS8gM7NxsmtXiae3vCVt8x5J62qW+yOiP3l8OLC55rktwPuHvP5y4KeSLgA6gZbcwNaBbWZW76WI6BvD65cA10XEX0s6GbhB0rERsWcsRaWaEmk2wZ60+bikDZLWS/rbsRRlZpZjzwJH1CzPTtbVWgqsAIiIe4EK0DPWjpsGdpoJdklzgEuBUyLiPcBfjLUwM7Ocuh+YI+koSWWqXyoODGnzDHA6gKR3Uw3sF8facZo97DQT7J8Cro6I7QARsW2shZmZ5VFEDALnA7cDG6keDbJe0hWSFibNLgQ+Jekh4Cbg3IiIsfadZg47zQT7OwEk/QJoAy6PiJ8M3ZCkZcAygPKhqSf/zcxyJSJWAauGrLus5vEG4JRW99uqLx1LwBxgHtX5nHsk/X5EvFrbKPmWtR+gctihY/6/jZnZgSTNlEiaCfYtwEBE7I6Ip4DHqQa4mZm1SJrATjPB/kOqe9dI6qE6RfJkC+s0MzvgNQ3slBPstwMvS9oArAEujoiXx6toM7MDUao57BQT7AF8JvkxM7Nx4GuJmJkVRGanppenDnLpLSuz6n5YF523mPb2zqzLaKj8yacodU3Nuow6g493ZF1CQzt+No2dlfy9l6defBed3fn8nV3ywXl0qJJ1Gfs45szurEvIjcwCW8CsWfm7lEl7eyeVSlfzhhl4s2s7lZ5pWZdRp7S1g1J3/n5n7ZWdlDtnZF1Gnc7uDroOyVco7tWhCpVSvj5jU9SWdQm54SkRM7OCcGCbmRWEA9vMrCAc2GZmBeHANjMrCAe2mVlBOLDNzArCgW1mVhAObDOzgnBgm5kVhAPbzGyUJM2X9JikTZIuafB8h6Rbkufvk3RkK/p1YJuZjYKkNuBq4CygF1giqXdIs6XA9oh4B/AV4Iut6NuBbWY2OicCmyLiyYjYBdwMLBrSZhFwffJ4JXC6JI21Ywe2mVm9Hknran6W1Tx3OLC5ZnlLso5GbZK7du0AZo61qPxd39TMbBxol+h4ppy2+UsR0Tee9ewP72GbmY3Os8ARNcuzk3UN20gqATOAMd/n1oFtZjY69wNzJB0lqQwsBgaGtBkAzkkefxS4K7n37Zh4SsTMbBQiYlDS+cDtQBtwbUSsl3QFsC4iBoDvADdI2gS8QjXUx8yBbWY2ShGxClg1ZN1lNY9fBz7W6n49JWJmVhAObDOzgnBgm5kVhAPbzKwgHNhmZgXhwDYzKwgHtplZQTiwzcwKwoFtZlYQDmwzs4LI7NT0N3a2cd7cM7Lqflhz+39OZ3dH1mU0NHDV8eysdGZdRp3B3lezLqGh8gc2U+qamnUZddZc3Et7e/7eR4DP3vgtumfmaz/ul6uPhpuyriIfMryWiKiUpmXX/TA6uzvoOqSSdRkNtVc6KXfOyLqMOqXpg5S6u7Iuo06pMpVKT/4+Y7vb26lU8vf7AuieOYVZs/J1iaFyecw3apk08vW/UjMzG1aqwG52h+Cadn8iKSTl7k4NZmZF1zSwU94hGEnTgU8D97W6SDMzS7eHneYOwQBfoHor99dbWJ+ZmSXSBHbTOwRLOgE4IiJuG2lDkpbtvQtxsGfUxZqZHcjG/KWjpCnAl4ELm7WNiP6I6IuIPvn7TjOzUUmTms3uEDwdOBa4W9LTwEnAgL94NLMDkaTfk7Ra0j8m/+0ept2XJK2XtFHSVZKaHr+YJrBHvENwROyIiJ6IODIijgTWAgsjYl2q0ZmZTS6XAHdGxBzgzmR5H5L+EDgFOI7qDu8fAKc123DTwI6IQWDvHYI3Aiv23iFY0sLRjMLM7ACwCLg+eXw98O8atAmgApSBDqAdeKHZhlOd0tTsDsFD1s9Ls00zsxzrkVQ7S9AfEf0pXzsrIrYmj58HZg1tEBH3SloDbAUEfD0iNjbbcL7OQTUzGydtu2D6byJt85ciYtjv4STdARza4KnltQsREZLqOpX0DuDdVL8TBFgt6Y8i4v+NVJQD28xslCJi2CvXSXpB0mERsVXSYcC2Bs0+DKyNiH9KXvNj4GRgxMD2sXVmZq01AJyTPD4HuLVBm2eA0ySVJLVT/cKx6ZSIA9vMrLX+CjhT0j8CZyTLSOqT9O2kzUrgCeAR4CHgoYj4v8027CkRM7MWioiXgdMbrF8H/Fny+E3gP412297DNjMrCAe2mVlBOLDNzArCgW1mVhAObDOzgnBgm5kVhAPbzKwgHNhmZgXhwDYzKwgHtplZQTiwzcwKwoFtZlYQDmwzs4JwYJuZFURml1ctTx3k0ltWZtX9sC46bzHt7Z1Zl9FQ+ZNPUeqamnUZdQYf78i6hIZ2/GwaOyv5ey9PvfguOrvz+Tu75IPz6FAl6zL2ccyZ3VmXkBuZBbaAWbPydznu9vZOKpWurMto6M2u7VR6pmVdRp3S1g5K3fn7nbVXdlLunJF1GXU6uzvoOiRfobhXhypUSvn6jE1RW9Yl5IanRMzMCsKBbWZWEA5sM7MWkvQxSesl7ZHUN0K7gyWtlPRrSRslndxs2w5sM7PWehT4CHBPk3ZfA34SEe8C3kuKu6bn71s/M7MCi4iNAJKGbSNpBnAqcG7yml3Armbb9h62mVm9Hknran6WtXj7RwEvAv9b0gOSvi2p6TGo3sM2swNC2+vBjCfeSNv8pYgYaf75DuDQBk8tj4hbU2y/BJwAXBAR90n6GnAJ8N+avcjMzEYhIs4Y4ya2AFsi4r5keSXVwB6Rp0TMzCZYRDwPbJZ0TLLqdGBDs9c5sM3MWkjShyVtAU4GbpN0e7L+rZJW1TS9ALhR0sPA8cD/arZtT4mYmbVQRPwA+EGD9c8BC2qWHwSGnSdvxHvYZmYF4cA2MyuIVIEtab6kxyRtklT3Taakz0jaIOlhSXdKenvrSzUzO7A1DWxJbcDVwFlAL7BEUu+QZg8AfRFxHNXDU77U6kLNzA50afawTwQ2RcSTyemTNwOLahtExJqI2JksrgVmt7ZMMzNLE9iHA5trlrck64azFPhxoyckLdt7qufgYKSv0szMWvulo6SzqR6mcmWj5yOiPyL6IqKvVBr+wihmZlYvzXHYzwJH1CzPTtbtQ9IZwHLgtIhIfcK+mZmlk2YP+35gjqSjJJWBxcBAbQNJc4FvAgsjYlvryzQzs6aBHRGDwPnA7VQvsL0iItZLukLSwqTZlcA04HuSHpQ0MMzmzMxsP6U6NT0iVgGrhqy7rObxWK9cZWZmTfhMRzOzgnBgm5kVhAPbzKwgHNhmZgXhwDYzKwgHtplZQTiwzcwKIrNbhL2xs43z5ubv8O25/T+ns7sj6zIaGrjqeHZWOrMuo85g76tZl9BQ+QObKXVNzbqMOmsu7qW9PX/vI8Bnb/wW3TPztR/3y9VHw01ZV5GepCuBfwvsAp4A/jQiGv6RJJevXgc8GxEfarbtDO/pKCqladl1P4zO7g66DqlkXUZD7ZVOyp0zsi6jTmn6IKXurqzLqFOqTKXSk7/P2O72diqV/P2+ALpnTmHWrHzd6rVcLtyF4lYDl0bEoKQvApcCfzlM209TPYM81QciX/8rNTMruIj4aXJJDxjh/gCSZgP/Bvh22m07sM3Mxs8nGeb+AMBXgc8Ce9JuLF//9jEzGyd6fRfljZubN6zqkbSuZrk/Ivr/ZVvSHcChDV63PCJuTdosBwaBG+tqkT4EbIuIX0mal7YoB7aZWb2XIqJvuCebXfBO0rnAh4DTI6LR7bVOARZKWgBUgC5J/ycizh5pu54SMTNrIUnzqU51LKy51+0+IuLSiJgdEUdSvcfAXc3CGhzYZmat9nVgOrA6uT/ANQCS3ipp1cgvHZmnRMzMWigi3jHM+ueABQ3W3w3cnWbb3sM2MysIB7aZWUE4sM3MCsKBbWZWEA5sM7OCcGCbmRWEA9vMrCAc2GZmBeHANjMrCAe2mVlBOLDNzArCgW1mVhAObDOzgnBgm5kVhAPbzKwgHNhmZgXhwDYzKwgHtplZQTiwzcwKIlVgS5ov6TFJmyRd0uD5Dkm3JM/fJ+nIVhdqZlYEkr4g6eHkBrw/lfTWBm2Ol3SvpPVJ23+fZttNA1tSG3A1cBbQCyyR1Duk2VJge3Lzya8AX0zTuZnZJHRlRBwXEccDPwIua9BmJ/CJiHgPMB/4qqSDm204zV3TTwQ2RcSTAJJuBhYBG2raLAIuTx6vBL4uSRERw210asdvOWTpvSm6n1i/ffmf+eff5nOmqPcPtzGlLX83uj+yczcq5a+u8p7XmbK9Lesy6sRHxJQp+asLYPWWMuUXlHUZ+3h55vezLmFUIuK1msVOoC4HI+LxmsfPSdoGvAV4daRtp/krOxzYXLO8BXj/cG0iYlDSDmAm8FJtI0nLgGXJ4hvfvezBR1P0XzQ9DBn3JDAZxwSTc1yTcExbAI4Z61ZeG3zx9p+88I2elM0rktbVLPdHRH/aviT9T+ATwA7gj5u0PREoA0802+6E7hYlA+4HkLQuIvomsv+JMBnHNRnHBJNzXJNxTFAd11i3ERHzW1ELgKQ7gEMbPLU8Im6NiOXAckmXAucDnx9mO4cBNwDnRMSeZv2mCexngSNqlmcn6xq12SKpBMwAXk6xbTOzwomIM1I2vRFYRYPAltQF3EY15Nem2Viaydr7gTmSjpJUBhYDA0PaDADnJI8/Ctw10vy1mdlkJWlOzeIi4NcN2pSBHwDfjYiVabfddA87mZM+H7gdaAOujYj1kq4A1kXEAPAd4AZJm4BXqIZ6M6nngwpmMo5rMo4JJue4JuOYoFjj+itJxwB7gN8A5wFI6gPOi4g/Az4OnArMlHRu8rpzI+LBkTYs7wibmRVDPo9fMzOzOg5sM7OCGPfAnoyntacY02ckbUhOOb1T0tuzqHO0mo2rpt2fSIpkTi7X0oxJ0seT92u9pL+d6Br3R4rP4NskrZH0QPI5XJBFnaMh6VpJ2yQ1PD9DVVclY35Y0gkTXWPmImLcfqh+SfkEcDTVA8MfAnqHtPnPwDXJ48XALeNZ0wSN6Y+BqcnjP8/7mNKOK2k3HbgHWAv0ZV13C96rOcADQHeyfEjWdbdoXP3AnyePe4Gns647xbhOBU4AHh3m+QXAjwEBJwH3ZV3zRP+M9x72v5zWHhG7gL2ntddaBFyfPF4JnC4pX+fG7qvpmCJiTUTsTBbXUj12Pe/SvFcAX6B6rZjXJ7K4/ZRmTJ8Cro6I7QARsW2Ca9wfacYVQFfyeAbw3ATWt18i4h6qR5kNZxHVw+AiqsctH5yceHLAGO/AbnRa++HDtYmIQaqncs4c57rGIs2Yai2luleQd03HlfwT9IiIuG0iCxuDNO/VO4F3SvqFpLWSWnY23DhKM67LgbMlbaF64sYFE1PauBrt396kk78r9kwiks4G+oDTsq5lrCRNAb4MnJtxKa1WojotMo/qv4TukfT7ETHiRXgKYAlwXUT8taSTqZ4ncWykOP3Z8mu897BHc1o7BTmtPc2YkHQGsBxYGBFvTFBtY9FsXNOBY4G7JT1NdQ5xIOdfPKZ5r7YAAxGxOyKeAh6nGuB5lmZcS4EVABFxL1ChemGoIkv1tzeZjXdgT8bT2puOSdJc4JtUw7oIc6LQZFwRsSMieiLiyIg4kurc/MKIGPNFecZRms/fD6nuXSOph+oUyZMTWeR+SDOuZ4DTASS9m2pgvzihVbbeAPCJ5GiRk4AdEbE166Im1AR887uA6l7LE1QvcgJwBdU/dqh+kL4HbAL+Hjg6629iWzCmO4AXgAeTn4Gsa27FuIa0vZucHyWS8r0S1ameDcAjwOKsa27RuHqBX1A9guRB4F9nXXOKMd0EbAV2U/2Xz1Kqp3WfV/NeXZ2M+ZEifP5a/eNT083MCsJnOpqZFYQD28ysIBzYZmYF4cA2MysIB7aZWUE4sM3MCsKBbWZWEP8fX4xTIOqaq/YAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "ramp_checkers = (x_idx-(nx//2))*(-1)**(x_idx+y_idx)\n",
    "print(ramp_checkers)\n",
    "\n",
    "q_ramp = Function(W, ramp_checkers)\n",
    "\n",
    "ax = plt.gca()\n",
    "l = tricontourf(q_ramp, axes=ax)\n",
    "triplot(mesh, axes=ax, interior_kw=dict(alpha=0.05))\n",
    "plt.colorbar(l)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Check that this new function `q_ramp` is orthogonal to the checkerboard `q`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "-4.5102810375396984e-17"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "#clear\n",
    "assemble(q_ramp*q*dx)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We would like to computationally check whether the inf-sup condition is obeyed.\n",
    "\n",
    "> There exists a constant $c_2 > 0$ so that (*inf-sup* or *LBB condition*):\n",
    "> $$ \\inf_{\\mu \\in M} \\sup_{v \\in X} \\frac{b (v, \\mu)}{\\|v\\|_X \\|\\mu\\|_M} \\geqslant c_2 .$$\n",
    "\n",
    "What are the $X$ and $M$ norms?\n",
    "<details style=\"background-color: #ddd\">\n",
    "<summary>Show answer</summary>\n",
    "$\\|\\cdot\\|_X=\\|\\cdot\\|_{H^1}$ and $\\|\\cdot\\|_M=\\|\\cdot\\|_L^2$.\n",
    "</details>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "How do we check the inf-sup condition?\n",
    "<details style=\"background-color: #ddd\">\n",
    "<summary>Show answer</summary>\n",
    "We're choosing a specific $\\mu$ (=`q`) here, so we need to check that a (mesh-independent\n",
    "$$ \\sup_{v \\in X} \\frac{b (v, \\mu)}{\\|v\\|_X }\n",
    "=\\sup_{v \\in H^1} \\frac{b (v, \\mu)}{\\|v\\|_{H^1} } \\ge c_2  \\|\\mu\\|_{L^2}\n",
    "$$\n",
    "So we should really be computing the $H^{-1}$ norm of the functional $b(\\cdot, \\mu)$.\n",
    "</details>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "How do we evaluate that quantity?\n",
    "<details style=\"background-color: #ddd\">\n",
    "<summary>Show answer</summary>\n",
    "Find the $(H^1)^2$ Riesz representer $\\b u$ of $f(\\b v):=b(\\b v, \\mu)$, i.e. $\\b u$ so that\n",
    "$$(\\b u,\\b v)_{(H^1)^2} = b(\\b v, \\mu) \\qquad(\\b v\\in (H^1_0)^2).$$\n",
    "Then, evaluate $\\|u\\|_{(H^1)^2}$.\n",
    "\n",
    "This works because\n",
    "$$\n",
    "\\|f\\|_{H^{-1}}\n",
    "=\\sup_{v\\in H^1}\\frac{|f(v)|}{\\|v\\|_{H^1}}\n",
    "=\\sup_{v\\in H^1}\\frac{|(u,v)|_{H^1}}{\\|v\\|_{H^1}}\n",
    "=\\|u\\|_{H^1}.\n",
    "$$\n",
    "Equivalently, we may evaluate $\\sqrt{f(u)}=\\sqrt{(u,u)_{H^1}}  =\\|u\\|_{H^1}$.\n",
    "</details>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Write a function with arguments <tt>(V, q)</tt> to evaluate that quantity."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "#clear\n",
    "def hminus1_norm(V, q):\n",
    "    # Evaluates the H^{-1} norm of b, where V is\n",
    "    # the function space for the first argument.\n",
    "    \n",
    "    u = TrialFunction(V)\n",
    "    v = TestFunction(V)\n",
    "    \n",
    "    riesz_rep = Function(V)\n",
    "    solve((inner(grad(u), grad(v)) + inner(u,v))*dx == div(v)*q*dx, riesz_rep)\n",
    "    \n",
    "    return norm(riesz_rep, \"H1\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Now computing nx=9...\n",
      "Now computing nx=19...\n",
      "Now computing nx=39...\n",
      "Now computing nx=79...\n",
      "Now computing nx=159...\n",
      "Now computing nx=319...\n"
     ]
    }
   ],
   "source": [
    "h_values = []\n",
    "h1_norms = []\n",
    "\n",
    "for e in range(6):\n",
    "    nx = 10 * 2**e - 1\n",
    "    print(f\"Now computing nx={nx}...\")\n",
    "    mesh = UnitSquareMesh(nx, nx, quadrilateral=True)\n",
    "    \n",
    "    V = VectorFunctionSpace(mesh, \"CG\", 1)\n",
    "    W = FunctionSpace(mesh, \"DG\", 0)\n",
    "\n",
    "    Z = V * W\n",
    "    \n",
    "    x = SpatialCoordinate(mesh)\n",
    "    x_w = interpolate(x[0], W)\n",
    "    y_w = interpolate(x[1], W)\n",
    "    \n",
    "    x_array = x_w.dat.data.copy()\n",
    "    y_array = y_w.dat.data.copy()\n",
    "\n",
    "    x_idx = (np.round(x_array * nx * 2) - 1)//2\n",
    "    y_idx = (np.round(y_array * nx * 2) - 1)//2\n",
    "    \n",
    "    odd_checkers = (x_idx-(nx//2))*(-1)**(x_idx+y_idx)\n",
    "    q_ramp = Function(W, odd_checkers)\n",
    "    \n",
    "    h_values.append(1/nx)\n",
    "    h1_norms.append(hminus1_norm(V, q_ramp)/norm(q_ramp, \"L2\"))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAbAAAAEnCAYAAADILRbRAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAAPYQAAD2EBqD+naQAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+j8jraAAAgAElEQVR4nO3deXxV1bn/8c+TMIQhCTIGBRVRQGQIqVKH1uJUqVO1giDV1ns7+gPHarXWatHWsUJBuXa4V4tVcBbFFlTEgaolKgQMQURBmcdAEoaEJOf5/bEPNITMOckZ8n2/XrzI2Vl77+fEI9/stddey9wdERGReJMU7QJEREQaQgEmIiJxSQEmIiJxSQEmIiJxSQEmIiJxSQEmIiJxSQEmIiJxSQEmIiJxqVW0C5CGMzMDDgeKol2LiEgEpQIbvJaZNhRg8e1wYF20ixARaQK9gPU1NVCAxbcigLVr15KWlhbtWkREGq2wsJDevXtDHXqWFGAJIC0tTQEmIi2OBnGIiEhcUoCJiEhcUheiiIhEVHnIyV6dz5aiYrqnpjC8T2eSkyzi51GAiYhIxMzN3cjE2XlsLCg+sK1negp3XjiQkYN6RvRc6kIUEZGImJu7kaufXHRQeAFsKijm6icXMTd3Y0TPpwATEZFGKw85E2fnUdWTx/u3TZydR3moxmeT60UBJiIijZa9Ov+QK6+KHNhYUEz26vyInVMBFofMbLyZ5QHZ0a5FRARgS1H14dWQdnWhAItD7j7N3QcCw6Ndi4gIQEqrusVJ99SUiJ1ToxBFRKRR5n+6mdteyq2xjQEZ6cGQ+khRgImISIPsLinj9/9czoyFawDomZbCxsJiDA4azLH/CbA7LxwY0efBFGAiIlJvi9bs4MZncvhy+x4AfvSNPtx8bn/eXrHlkOfAMproOTCrZbkViWFmlgYUFBQUaDJfEWkWpeUhpr65kmlvfU7Ig4eU/zB6KKcd2/VAm8bMxFFYWEh6ejpAursX1tRWV2AiIlInn2/ZxQ3P5PDJ+gIALs48nInfHUR6u9YHtUtOMk7p26XJ61GAiYhIjUIh54kPvuTeOZ9SUhYivV1rfn/JIC4YcnhU61KAiYhItTYVFHPz80tYsHIbAN88risPjhpKRnrkhsM3lAJMRESqNHvJBm6flUvB3lLatkri1+cfz5UnH4VZ5GeWbwgFmIiIHKRgTyl3vJLLyzkbABjSK51Jl2VybPeOUa7sYAowERE54L3Pt3HTc0vYWFBMcpIx/oxjuebMY2mdHHsTNynAYoiZvQSMAN5091FRLkdEWpDi0nLun/spj7/3JQBHd2nPpDGZZB15WHQLq4ECLLZMAR4DfhjtQkSk5chdX8D1z+Tw+ZZdAHz/60fy6/OPp32b2I6I2K6uhXH3t81sRLTrEJGWoTzk/OmdL5j8xmeUhZyuHdvy4KghnDGge7RLq5PY69RsADM7wsyeNLPtZrbXzD4xsxMjePzTzWy2mW0wMzezi6tpN97MvjSzYjNbaGaaLV5EYtJX23dz2Z8/4MHXVlAWckaekMHrN5weN+EFCXAFZmaHAe8BbwHfAbYCxwE7qml/GpDt7qWVtg8Etrv75ip26wAsIejee7Ga444BJgE/BxYC1wOvmVl/d98SbpND1T/zb7v7hlreqohIo7k7z3y4lrtezWPPvnI6tm3FxItO4HtZR8TM8Pi6ivsAA24B1rr7f1XYtrqqhmaWBEwDVprZWHcvD2/vD8wnCKAHKu/n7nOAOeG21dVxI/BXd3883O7nwPnAfwP3hY+TWd83JyISKVuLSvjVi0uZt3wLAMP7dOah0UPp3bl9lCtrmEToQrwI+MjMnjOzLWa22Mx+UlVDdw8B5wHDgCfMLMnM+hKE1yx3PyS86sLM2gBfA+ZVOtc84JSGHLOW82lFZhGpl9eXbWLkH99l3vIttElO4rbzBjDzJyfHbXhBYlyBHQNcTXD1dA9wEjDVzPa5+/TKjd19g5mdCSwAZhAEzLzwMRqqK5AMVO5+3AwMqOtBzGweMBToYGbrgNHu/kHldu4+DZi2fzb6BlctIglvV0kZd8/O45mP1gIwICOVyWMyOb5n/K9gkQgBlgR85O63hV8vNrNBBPeiDgkwAHdfY2ZXAu8Aq4AfeQysK+PuZ0e7BhFJHB9+mc+Nz+awNn8vZvDTbx7Djd/uR9tWydEuLSISIcA2AnmVti0HLq1uBzPrAfwFmE1wxTYZuKYRNWwDyoEelbb3ADY14rgiIvW2ryzE5Hmf8ad3vsAdjujUjkmXDeXrxzT9EifNKREC7D2gf6Vt/YCvqmpsZl2BNwlCbnS47dtmVuLuNzWkAHffZ2YfA2cBs8LnSQq/fqQhxxQRaYgVm4q44Zkc8jYGa0FemtWL3140kNSU1rXsGX8SIcAmA++b2W3As8Bw4KfhPwcJh8ocgnAb4+5lQJ6ZnQPMN7P17j65iv06AsdW2NTHzDKBfHdfE942CZhuZh8RDK64nmD4/eMRep8iItUKhZzH3lvNA6+tYF9ZiMPat+be7w1m5KCe0S6tyVgM3PppNDO7ALiX4Pmv1cAkd/9rNW3PARa4e3Gl7cOAre6+rop9RhA8Z1bZdHe/qkK7CcDNQAaQA1zr7gsb8p7qYv8gjoKCAtLS4v+GrIg0zPqde7np2SV8sGo7AGf078b9o4bQPTX6a3bVV2FhIenp6QDp7l5YU9uECLCWSgEm0rK5O7Ny1nPHy8soKi6jXetkbr/geMYNPzLuHkrerz4BlghdiCIiLc7OPfv49Uu5/OOTjQBk9u7E5DGZ9OnaIcqVNR8FmIhInHnns63c/NwSthSVkJxkXHfWcfy/EX1pFYNrdjUlBZiISJzYu6+ce+cs54kPgkHWx3TrwB/HZDKkV6coVxYdCjARkTiwZO1Obngmh1XbdgNw1alHc8vIAbRrkxgPJTeEAkxEJIaVlYeY9tYXTJ2/kvKQ0yOtLQ+OGsrp/bpFu7SoU4CJiMSoVVt3ccOzS1iydicA5w/pye8vHkSn9m2iXFlsUICJiMQYd+fJhWu45x/L2VtaTmpKK3538SAuGnp43A6PbwoKMBGRGLKlsJhfvrCUt1dsBeDUvl34w+ihHN6pXZQriz0KMBGRGDHnk43c9tIn7NhTSptWSdwycgD/derRJCXpqqsqCjARkSgrLC7lt68s48VF6wE44fA0Jo/JpF+P1ChXFtsUYCIizaA85GSvzmdLUTHdU1MY3qczyUnGv1dt5xfPLmH9zr0kGVw9oi/XndWPNq1a1kPJDaEAExFpYnNzNzJxdh4bC/4zh3hGWlsG9+rEvOWbcYcjO7dn0mVDOfHozlGsNL4owEREmtDc3I1c/eQiKk+bvqmwhE15mwEYe1Jvbr9gIB3b6p/k+tBPS0SkiZSHnImz8w4Jr4oOa9+a318ymGQN1Kg3dbLGITMbb2Z5BAtnikiMyl6df1C3YVV27Ckle3V+M1WUWBRgccjdp7n7QILVp0UkRm0pqjm86ttODqYAExFpIp3at65Tu3hcOTkW6B6YiEgT+GxzEb9/dXmNbQzISA+G1Ev9KcBERCLI3ZmRvYa7ZudRUhYiNaUVRcVlGBw0mGP/kI07LxyoARwNpAATEYmQnXv2cesLnzB32SYATu/XjYdGD+Xjr/IPfQ4sPYU7LxzIyEE9o1Vu3DP3mgZ4SiwzszSgoKCggLS0tGiXI9KiLVy1neufyWFjQTGtk41fnjuAH32jz4F5DKubiUMOVlhYSHp6OkC6uxfW1FZXYCIijVBWHuLh+Z/z8PyVhByO7tKehy/PYnCv9IPaJScZp/TtEqUqE5MCTESkgdbv3Mv1Ty/mwy93AHBpVi8mfvcEzajRTPRTFhFpgDmfbOSWF5ZSWFxGx7bBgpMXDzsi2mW1KAowEZF62LuvnLtezWNm9hoAhvbuxNSxmRzVpUOUK2t5FGAiInW0fGMh185czMotuzCDn3+rLzee04/WyZoTIhoUYCIitXB3/v7vr/jdP5azryxEt9S2TL4sk28c1zXapbVoCjARkRrk797HL59fyrzlwdInZ/Tvxh9GD6VLx7ZRrkwUYCIi1Xj/i23c8EwOmwtLaJOcxK3fGcB/nXY0Znp+KxYowEREKiktDzFl3kqmvf057nBMtw48fPkwTjg8vfadpdkowEREKlibv4frnl7MojU7ARhzYm/uvGgg7dvon8tYo/8iIiJhs5ds4LYXP6GopIzUtq2453uDuXDo4dEuS6qhABORFm/PvjJ++8oynv1oHQBZR3Ziythh9O7cPsqVSU0UYCLSouWuL+DapxezautuzGDCGcdy3VnH0UrPdsU8BZiItEjuzmPvfcn9cz5lX3mIHmltmTwmk1P76tmueKEAiyFm9hIwAnjT3UdFuRyRhLV9Vwk3PbeEt1ZsBeDs43vwwKghdO7QJsqVSX0owGLLFOAx4IfRLkQkUf1r5TZueDaHrUUltGmVxO3nH8+VJx+lZ7vikAIshrj722Y2Itp1iCSi0vIQD73+GX9+9wvc4bjuHZl6+TCO76nFYONVQt2lNLNbzczN7I8RPu7pZjbbzDaEj39xNe3Gm9mXZlZsZgvNbHgk6xCRhvlq+25G/ekD/vROEF7jvn4kr0z4hsIrziXMFZiZnQT8DFhaS7vTgGx3L620fSCw3d03V7FbB2AJQffei9UcdwwwCfg5sBC4HnjNzPq7+5Zwmxyq/pl/29031FS3iDTMrMXruX1WLrtKykhLacUDo4YwclDPaJclEZAQAWZmHYGngJ8At9fQLgmYBqw0s7HuXh7e3h+YTxBAD1Tez93nAHPCbas7/I3AX9398XC7nwPnA/8N3Bc+TmYD3l5V72M8MJ4Eu4IWiaRdJWXc8XIuLy5aD8BJRx/GH8cO44hO7aJcmURKovwDOA34h7vPq6mRu4eA84BhwBNmlmRmfQnCa5a7HxJedWFmbYCvAQfOHz7XPOCUhhyzJu4+zd0HAuqiFKnC0nU7uWDqAl5ctJ4kg+vPPo6ZPzlZ4ZVg4v4KzMzGAlnASXVp7+4bzOxMYAEwgyBg5gFXN6KMrkAyULn7cTMwoK4HMbN5wFCgg5mtA0a7+weNqEukRQmFnP/712oeeO1TSsudw9NT+OPYYQzv0znapUkTiOsAM7PeBEPPz3H34rru5+5rzOxK4B1gFfAjd/cmKrPO3P3saNcgEq+2FBXzi2eXsGDlNgBGnpDBfZcOplN7PduVqOI6wAi67boDiyrcm0oGTjezCUDb/fe5KjKzHsBfgNkEV26TgWsaUcc2oBzoUWl7D2BTI44rInXw9oot3PTcErbt2kfbVkncceFAxg0/Us92Jbh4D7A3gcGVtj0OfArcX014dQ3vtxwYDfQD3jazEne/qSFFuPs+M/sYOAuYFT5PUvj1Iw05pojUbl9ZiAdf+5S/LlgNwICMVKZePox+PVKjXJk0h7gOMHcvAnIrbjOz3QTD4XMrtw+HyhzgK2CMu5cBeWZ2DjDfzNa7++Qq9usIHFthUx8zywTy3X1NeNskYLqZfQRkEwyj70AQqCISYau37ebamYv5ZH0BAD845ShuO+94UlonR7kyaS5xHWD15e4hM7sNWODu+ypsX2JmZwNbq9n1ROCtCq8nhf+eDlwVPsYzZtYNuAvIAHKAkdU8VyYiDeTuvLBoPXe8nMuefeV0at+aBy4dwrdPyIh2adLMLAbGLkgDmVkaUFBQUEBammYUkMRXVFzK7bNyeTkneO7/5GM6M3lMJj3TNTw+URQWFpKeng6Q7u6FNbVtsiswMyt3d13Li0hELF6zg2ufXsza/L0kJxk3nH0cV484luQkDdRoqZqyC1GfKhFptFDI+dO7XzDp9c8oCzlHdGrH1Msz+dpRerarpWvKAFPfpIjUWXnIyV6dz5aiYrqnpjC8T2e27SrhxmdzeO/z7QCcP6Qn91wymPR2raNcrcSCegWYmX0dGAecSjBQYS/BcPQ5wEx3L4h4hSKS8ObmbmTi7Dw2FvxnPoLD2remtNzZVVJGu9bJTLzoBEaf2EvPdskBdQ4wM3sVWEfw8O99BCP2UgiGl38LeN7MHnb3V5qiUBFJTHNzN3L1k4sO6bLZsSdYMOKITu2Y/t/DObZ7x+YvTmJafa7ArnD3nZW27SIYLp4DTDGzThGrTEQSXnnImTg7r8b7DeXu9OnaodlqkvhR59no94eXmb1hZg+Z2VVmlmVmbSu3ERGpi+zV+Qd1G1ZlU0Ex2avzm6kiiScNGcTxCXAMUAqMAc40s8+BZcBSd/9dBOsTkQS2pahuc3DXtZ20LA0JsHPc/cD8g2Z2HnAy8CowJFKFiUji211SVqd23VNTmrgSiUcNWdCyILyCMQDu/k/gQnfPdvf/jVxpIpKo3J3p73/JnS8vq7GdAT3TU7Sel1SpIVdgPwOeNbN3gCUECzbW7dcoEWnxCvaU8ssXlvDasmCa0CFHpB+YkLfiYI79g+XvvHCgZtuQKtX7CszdlxGsofUv4GhgA3BeZMsSkUT08Vc7OG/qAl5btpnWycYdFwzk5Qmn8egVWWSkH9xNmJGewqNXZDFyUM8oVSuxrskm890/F6KZvUXwi1XlX6H2n9jCX//N3Z9okmISlCbzlXgRCjl/fncVf3h9BeUh56gu7Xn48mEM6fWfJ2+qmolDV14tT0xM5rufu5/R1OcQkdgVTAe1hHc/C1YrunDo4dxzySBSUw6eDio5yTilb5dolChxqtEBZmY9CRZ2LKml3RvAUoJh+EuBZbXtIyLx7f3Pt3HdMzlsLSqhbaskJl50AmNO6q3poCQiInEF9negr5m94O431dBOz4+JtBBl5SGmvrmSh9/6HHc4rntHHhmXRf+M1GiXJgmk0QHm7mcDmNmAWprq+TGRFmBTQTHXPr34wOwZY07szW8vOoF2bbQ8oERWvUchmtl0Mzum8nZ3/7SWXfX8mEiCm//pZr4z5V2yV+fToU0yU8Zmcv+oIQovaRINuQKbCTxhZp8Bd7v76jrup+fHRBLUvrIQD772KX9dEPxzcMLhaTwyLkuT8EqTqneAuftcYG64C/BJM8sjCLI1tey3zMxOAi4GBqPnx0QSwtr8PUyYuZgla4O5vK869Wh+dd4A2rbSVZc0rcbcA3sNWAX8AviMYG2wGrn7PuDZ8B+pxMxeAkYAb7r7qCiXI1Krf36ykVteWEpRcRlpKa14cPRQzj0hI9plSQtR7wAzs1eAfkAbYCXwKXB9hOtqqaYAjwE/jHYhIjUpLi3n7lfzeGph0PGSdWQnpl4+jF6HtY9yZdKSNOQK7HfAcncvqkvjCjNxVNsEzcQBgLu/bWYjol2HSE0+37KLCTMW8emm4J+Aq0f05cZz+tE6uSFzg4s0XEPugWXXs32TzsRhZlcDVxPMywjBc2V3ufucCJ7jdOBm4GtAT+ASd59VRbvx4XYZBANVrqnvz0sklj3/8Tp+MyuXvaXldOnQhkljMvlWv27RLktaqDr/yhReifnHZtat0vYkMzvFzP7HzK6KeIW1WwfcShAuJwLzgZfN7ISqGpvZaWbWuortA82sRzXn6EAQSOOrK8LMxgCTgIlAVrj9a2bWvUKbHDPLreLP4XV6pyJRsrukjBufzeGm55awt7ScU/t2Yc5131R4SVTVeTJfM+sA/BgYB3QHdgDtCAZvvAM86u4LK7Qvd/eoDEMys3zgZnf/v0rbk4BFBPfuxrp7eXh7f4L3MMndH6jl2E4VV2BmthD40N0nVDjXWuBhd7+vHrWPACbUZRCHJvOV5pC3oZAJMxexautukgxuOLsf/++MYzXRrjSJJpnM1913EwwymGJmbYAuQLG772hMsZFkZsnAaIIrpg8qf9/dQ+Hh/+8SPMt2JdCH4KptVm3hVcN52xBcAd5b6VzzgFMacsxazjee4GpQNx2kybg7Ty5cw92v5rGvLERGWgpTxmby9WM04a7EhoaMQuwH3AR0BZaa2VR3z494ZfWraTBBYKUAuwiukPKqauvuG8zsTGABMIMgYOYR3EdrqK5AMrC50vbNBA9s10k48IYCHcxsHTDa3asK4mnAtP1XYA2uWqQaBXtLufWFpczJ3QTAmQO684fRQ+ncoU2UKxP5j4aMQnyR4Eosh+Bez1wzu9Xd50e0svpZAWQC6cAoYLqZfauGEFsTvvp6h+BZth95Uy2MVg/755UUiaactTuZMGMR63bspXWyccvIAfzoG300g7zEnIYEWLG7/zX89Yfhh2/fILhyiIrwA9Kfh19+HJ7x4zqC6asOER6s8RdgNsHq0pOBaxpRwjagHKg8CKQHsKkRxxVpNqGQ83//Ws39cz+lLOT07tyOhy/PIrN3p9p3FomCOgeYmU0luOp6y8wmuPsj4W9tB0JNUVwjJAFtq/qGmXUF3gSWE9wv6we8bWYltSwHUy1332dmHwNnAbPC50kKv36kpn1FYsH2XSXc9NwS3loRLDp5/uCe3HvpYNJSDhmwKxIz6nMFNpdg2ZMjgO+Y2Q0EIdCXYEmUqDCze4E5wBoglWCU5Ajg3CraJoXbfgWMcfcyIM/MzgHmm9l6d59cxX4dgWMrbOpjZpkEC3nunwNyEkHX5UdANsHsJB2AxyPyRkWayL9Xbee6pxezuTBYdPKOCwcybviR6jKUmFefUYj/BP65/3V45N1Agq7DQeFtFm7rBDNsNIfuwBMEDxgXEKz2fK67v1G5YXhk4G3AgnC34/7tS8zsbGBrNec4EXirwutJ4b+nA1eFj/FM+Bm5uwgeZM4BRrp75YEdIjGhPOQ8PH8lU99cScihb7cOTPt+FgMy9EiGxIc6PwdW40HMfkxwxXFceNNKYEqFe2UV297i7vc3+qSi58CkwTYXFnP90zl8sGo7AKO+1ou7vnsC7dtEYpF2kYZrkufAqmNmdwE3ENzr2T/k+xTgITPr7e53VNplsJm9DPzA3QsqHKcz8BOFm0jTenvFFn7x7BK2795H+zbJ/O7iQXwvq1e0yxKpt0j8unU18FN3n1lh2ytmthR4GDgowNz9ivD8he+b2fcJuhqvBU5D94tEmkxpeYg/vL6CP7+zCoDje6bxyLhh9O3WMcqViTRMJAKsNfBRFds/ruH4/0swKOIjgiHo44Ef75/aSUQia23+Hq59ejGL1wSLTl558lH8+vzjSWmtRSclfkViKqK/U/UsFj8Fnqq80cx+R7CGWAfgdIIQO49gJgsRibC5uZs4f+oCFq/ZSWpKKx79fhZ3XzxI4SVxL1J3bH9kZt8G/h1+/XXgSIL5BveP2HN3/wXBc2NZFe5/XWBmvyHoUhzt7qsjVJNIi1ZcWs69/1zO9A++AiCzdycevnwYvTtr0UlJDI0ehRhesLIu3N3PrOE45wDT3L1fowpqQTQKUaqzausuJsxYTN7GYBDXz04/hpvO7a9FJyXmNesoxEgtWOnub4Qn2RWRRnhp8TpufymX3fvK6dyhDQ9dNpQz+nevfUeROBNTD324+7po1yASr/bsK+POl5fx3MfB/0YnH9OZKWOH0SMtJcqViTSNSDwHVvk5r4O4+11V7PMbd787/PVv3f23lbeLSN19uqmQCTMW8/mWXSQZXHvWcVxz5nFadFISWiTugS2utKk1wSKRZcAX7p5VxT7Z7j48/PWi/W0qbpfa6R5Yy1IecrJX57OlqJjuqSkM79OZJIOZ2WuZOHsZJWUhuqe2ZcrYYZzSV4tOSnxq7ntgwypvC//D+jfgpWp2szp8LSJhc3M3MnF2HhsLig9s65HWliMOa8eir4Jnu0b078ZDo4fSpWOVCzGIJJwmuQfm7oVmdifBelt/r6pJHb4WEYLwuvrJRYf8z7G5sITNhSUkGdwycgA/+eYxJKnLUFqQphxTmx7+U5VWZrY/PA0OzG4fU4NKRKKtPORMnJ1X4292h3Vow48VXtICRWIQx7WVNxEsbXIlwdpbVZkB3A38iv9cdU0EZlbTXqRFyl6df1C3YVW279pH9up83feSFicSVzw3VHodIlhXazpwbzX7PATca2ZvA73NbD6w0N0fjEA9IgljS1HN4VXfdiKJJBKDOPo0YB8Hbg2vdHwU8JW772psLSKJplsdB2R0T9WzXtLyRKILsR3BcPw94ddHAZcAee7+ek37hkNrWWNrEElE23eV8Od3v6ixjQEZ6cGQepGWJhKDOF4GfgBgZp2AbOAXwMvhdb+qZGbTzeyYCJxfJOFkr87nvKkLeOezbbQKD86oPERj/+s7LxyoB5alRYpEgGUBC8JfjwI2EXQL/oBgocrqzCSYrf4xM6t3N6RIIgqFnEfmr2TsXz5gc2EJfbt14NVrv8GfrsgiI/3gbsKM9BQevSKLkYN6RqlakeiKxCCO9kBR+OtvAy+6e8jM/k0QZFVy97nAXDM7D3jSzPKAu919TQRqEok7W4tKuPHZHBas3AbA97KO4O7vDqJD21YMyEjjnIEZh8zEoSsvackiEWCfAxeb2UvAucDk8PbuQI3TgIS9Bqwi6Hb8DNDdaGlx3v9iG9c9ncPWohJSWidx93cHMfrE3ge1SU4yDZUXqSASAXYXwXNdk4E33f2D8PZvA5XnSTzAzF4B+gFtgJUEqzRfH4F6ROJGech5eP5Kpr65kpDDcd078j/fz+K4HqnRLk0k5kViGP3zZvYvgoeXl1T41ptUPxciwO+A5e5eVEMbkYS1paiY65/O4f0vtgNw2Ym9mHjRINq1SY5yZSLxISJTN7n7JoLBGxW3ZdeyT43fF0lk/1q5jeufWcy2Xfto3yaZ3108iO9l9Yp2WSJxRXMPijSjsvIQU95cySNvfY47DMhI5ZFxWRzbvWO0SxOJOwowkWayubCYa2YuJnt1PgCXDz+SOy8cSEprdRmKNIQCTKQZvPPZVm54Jof83fvo0CaZey8dwkVDD492WSJxTQEm0oTKykM89MZnPPp2MCXUwJ5pTPt+Fn26dohyZSLxTwEm0kQ27NzLtTMX89FXOwC48uSj+PX5x6vLUCRCFGAiTWD+p5u58dkl7NxTSmrbVtx36RDOH6Ipn0QiSQEmEkGl5SEefG0Ff3l3FQCDj0jnkXHDOKqLugxFIk0BJhIh63bs4ZqZi1m8ZicAV516NL86bwBtW6nLUKQpKKxvYn8AAA57SURBVMBEIuD1ZZu4+fmlFOwtJS2lFQ+MGsrIQRnRLkskoSnARBphX1mI++Z8ymPvrQZgaO9OPHL5MHp3bh/lykQSnwIshoRn9B9BMCnyqCiXI7VYm7+HCTMWsWRdAQA//kYffjlyAG1aRWKZPRGpjQIstkwBHgN+GO1CpGZzczdy8/NLKSouI71dax4aPZSzB/aIdlkiLYoCLIa4+9tmNiLadUj1SsrKuecfy5n+wVcAZB3ZiYfHZXFEp3ZRrkyk5Yn7vg4z+5WZfWhmRWa2xcxmmVn/CJ/jdDObbWYbzMzN7OJq2o03sy/NrNjMFprZ8EjWIdH15bbdXPro+wfC62ffOoZnfnaKwkskShLhCuxbwDTgQ4L3cw/wupkNdPfdlRub2WlAtruXVto+ENju7purOEcHgrXOHgNerKoIMxsDTAJ+DiwkWJzzNTPr7+5bwm1yqPpn/m1331CXNyvR8erSDdz6wifsKinjsPatmXRZJmcM6B7tskRatLgPMHcfWfG1mV0FbAG+Brxb6XtJBGG30szGunt5eHt/YD5BAD1QxTnmAHPCbasr5Ubgr+7+eLjdz4Hzgf8G7gsfJ7Mh71Gip7i0nLtfzeOphWsAOOnow5h6+TB6puuqSyTa4r4LsQrp4b/zK3/D3UPAecAw4AkzSzKzvgThNcvdDwmvujCzNgSBOa/SueYBpzTkmLWcb7yZ5QFaFLQJrdq6i0v+532eWrgGMxh/Rl9m/uRkhZdIjIj7K7CKwldYfwTec/fcqtq4+wYzOxNYAMwgCJh5wNWNOHVXIBmo3P24GRhQ14OY2TxgKNDBzNYBo939g8rt3H0aMM3M0oCCBlct1Xo5Zz23vfgJu/eV06VDGyaPyeT0ft2iXZaIVJBQAUbQPTgI+EZNjdx9jZldCbwDrAJ+5O7eDPXVyN3PjnYNLV1xaTm/fWUZT3+4FoCTj+nMlLHD6JGWEuXKRKSyhAkwM3sEuAA43d3X1dK2B/AXYDZwEjAZuKYRp98GlAOVHwTqAWxqxHGlGX2+ZRfjn1rEis1FmME1Zx7HdWcdR3JStfc9RSSK4j7ALBhV8TBwCTDC3VfX0r4r8CawHBgN9APeNrMSd7+pITW4+z4z+xg4C5gVPk9S+PUjDTmmNK8XPl7H7bNy2VtaTteObZkyNpPTju0a7bJEpAZxH2AE3YbjgO8CRWa2fwbVAnffW7FhOFTmAF8BY9y9DMgzs3OA+Wa23t0nVz6BmXUEjq2wqY+ZZQL57r4mvG0SMN3MPiIYXHE9wfD7xyP1RiXy9uwr446Xl/H8x8FF+2nHdmHymEy6p6rLUCTWWQzc+mkUM6vuDfyXu/+tivbnAAvcvbjS9mHA1qq6H8OzY7xVxTmmu/tVFdpNAG4GMoAc4Fp3X1i3d1J/+wdxFBQUkJaW1lSnSVifbS5i/FOLWLllF0kG15/dj/FnHKsuQ5EoKiwsJD09HSDd3Qtrahv3AdaSKcAaxt157qN13PFKLsWlIbqntmXK2GGc0rdLtEsTafHqE2CJ0IUoUme7S8q4fVYuLy1eD8A3j+vK5DGZdO3YNsqViUh9KcCkxVi+sZDxMxaxautukpOMG8/px9Xf6kuSugxF4pICTBKeuzMzey0TZy+jpCxERloKD48bxklHd452aSLSCAowSSjlISd7dT5biorpnprC8T1T+c3Ly5i9JJgr+Yz+3Xjoskw6d2gT5UpFpLEUYJIw5uZuZOLsPDYW/GeAaXKSUR5ykpOMX57bn5988xh1GYokCAWYJIS5uRu5+slFVB5TWx4Kttx4Tj9+9q2+zV+YiDSZRJyNXlqY8pAzcXbeIeG1nwFP/vurA2EmIolBASZxL3t1/kHdhpU5sLGgmOzVh6ywIyJxTAEmcW9LYfXhdVC7orq1E5H4oACTuFZYXMpTC7+qU1vNbyiSWDSIQ+LWJ+sKGD9jEWvy99TYzoCM9BSG99FzXyKJRFdgEnfcnenvf8mlj77Pmvw99DqsHbeOHIARhFVF+1/feeFATdIrkmB0BSZxpbC4lFueX8qc3GCd0HNP6MEDo4aS3q41R3dtf8hzYBnpKdx54UBGDuoZrZJFpIloNvo41tJmo1+6bicTZixmTf4eWicbt513PFedejTBmqaByjNxDO/TWVdeInFEs9FLQtnfZfj7fy6ntNzpdVg7po3LYmjvToe0TU4yLYsi0kIowCSmFewNugznLju0y1BEWjYFmMSspet2Mn7GItbm7622y1BEWi4FmMScyl2GvTu345HLq+4yFJGWSwEmMaVyl+HIEzK4f9QQdRmKyCEUYBIzlqzdyYSZ/+ky/PV5x/NDdRmKSDUUYBJ17s7f3v+Se9RlKCL1oACTqCrYW8ovn1/Ca8s2A+oyFJG6U4BJ1OSs3cmEGYtYtyPoMrz9/IH84JSj1GUoInWiAJNm5+48/t6X3DvnP12G08ZlMaSXugxFpO4UYNKsCvaUcvPzS3g9L+gy/M6gDO67VF2GIlJ/CjBpNhW7DNskJ/Hr849Xl6GINJgCTJpc5S7DIzu3Z9q4LAb3So92aSISxxRg0qSq6jK8f9QQ0lLUZSgijaMAkyaTs3Yn459axPqdQZfh7Rccz5Unq8tQRCJDASYR5+489t6X3KcuQxFpQgowiaiCPaXc9PwS3gh3GZ43OBhlqC5DEYk0BZhEzOI1O5gwY7G6DEWkWSjApNHcnf/712rum/MpZSHnqC5Bl+GgI9RlKCJNRwEmjbJzzz5uem4p85YHXYbnD+7JvZcOVpehiDQ5BZg0WOUuw99ccDxXqMtQRJqJAkzqTV2GIhILFGBSL+oyFJFYoQCTOlu0ZgfXVOwyvHAgV3z9SHUZikhUKMCkVuoyFJFYpACTGgVdhkuYt3wLAOcP6cl93xtMqroMRSTKFGBSrYO6DFslcccFA/m+ugxFJEYowOQQ7s7/LljN/XODLsOju7TnEXUZikiMUYDJQSp3GV4wpCf3qstQRGKQAkwO+PirHVw7U12GIhIfFGAtUHnIyV6dz5aiYrqnpnDiUYfx+PureWDuCnUZikjcUIC1MHNzNzJxdh4bC4oPbGvbKomSshCgLkMRiR8KsBZkbu5Grn5yEV5p+/7wGje8N7+/ZLC6DEUkLiRFuwBpHuUhZ+LsvEPCq6K3VmwlVFMDEZEYogBrIbJX5x/UbViVjQXFZK/Ob6aKREQaRwHWQmwpqjm86ttORCTaFGAtRPfUlIi2ExGJNgVYCzG8T2d6pqdQ3fAMA3qmpzC8T+fmLEtEpMEUYC1EcpJx54UDAQ4Jsf2v77xwIMlJGoEoIvFBAdaCjBzUk0evyCIj/eBuwoz0FB69IouRg3pGqTIRkfozd42bjldmlgYUFBQUkJaWVuf9Ks/EMbxPZ115iUhMKCwsJD09HSDd3QtraqsHmVug5CTjlL5dol2GiEijqAtRRETikgJMRETikroQE0BhYY3dxCIicaM+/55pEEccM7MjgHXRrkNEpAn0cvf1NTVQgMUxC6aNPxwoquMu2cDwJionEsduzDHqu2992telbU1tUgl+0ehF3f9bxZOm/FxF+/yJ+rmORLum/FynAhu8loBSF2IcC//HrfE3lIrMLFTbsNSGisSxG3OM+u5bn/Z1aVtTmwrL0xQ11c8/mprycxXt8yfq5zoS7Zr4c12n42kQR8syLcaP3Zhj1Hff+rSvS9um/NnGumi/d32u698+0u2iQl2IIk1s/wPn1OHBTJF4EQufa12BiTS9EmBi+G+RRBH1z7WuwEREJC7pCkxEROKSAkxEROKSAkxEROKSAkxEROKSAkxEROKSAkwkhphZbzN728zyzGypmY2Odk0ijWVmL5nZDjN7PqLH1TB6kdhhZj2BHu6eY2YZwMdAP3ffHeXSRBrMzEYQzG/4Q3cfFanj6gpMJIa4+0Z3zwl/vQnYBnSOblUijePub9MEE1krwETqwcxON7PZZrbBzNzMLq6izXgz+9LMis1soZk1aCZyM/sakOzuaxtduEg1mvMzHWkKMJH66QAsAcZX9U0zGwNMIphiJyvc9jUz616hTY6Z5Vbx5/AKbToDTwA/bcL3IgLN9JluCroHJtJAZubAJe4+q8K2hcCH7j4h/DoJWAs87O731fG4bYE3gL+6+98jX7lI1ZrqMx3ebwQwQffARGKQmbUBvgbM27/N3UPh16fU8RgG/A2Yr/CSaIvEZ7opKcBEIqcrkAxsrrR9M5BRx2OcBowBLg53y+SY2eAI1ihSH5H4TGNm84DngPPMbJ2ZRST8tCKzSAxx93+hXywlwbj72U1xXP2PIhI524ByoEel7T2ATc1fjkijxfRnWgEmEiHuvo/gweOz9m8L3/A+C/ggWnWJNFSsf6bVhShSD2bWETi2wqY+ZpYJ5Lv7GoLhxtPN7CMgG7ieYJjy481erEgdxPNnWsPoReohPBT4rSq+Nd3drwq3mQDcTHCTOwe41t0XNleNIvURz59pBZiIiMQl3QMTEZG4pAATEZG4pAATEZG4pAATEZG4pAATEZG4pAATEZG4pAATEZG4pAATEZG4pAATEZG4pAATkVqZ2R/MbFbtLUWajwJMROoiE1ga7SJEKlKAiUhdDAWWRLsIkYoUYCJSIzPrRbC0PGb2hpntMbMVZvb1KJcmLZwCTERqkxn+ezxwD8HV2BrgvqhVJIICTERqlwnkA5e5+1vuvhJ4BegW3bKkpVOAiUhtMoGX3X1bhW19gM+jVI8IoAATkdplAv+uYltOFGoROUABJiLVMrNU4BhgcaVvKcAk6hRgIlKToUA58Mn+DWZ2FHAYCjCJMgWYiNQkE1jh7sUVtg0Ddrr7l9EpSSRg7h7tGkREROpNV2AiIhKXFGAiIhKXFGAiIhKXFGAiIhKXFGAiIhKXFGAiIhKXFGAiIhKXFGAiIhKXFGAiIhKXFGAiIhKXFGAiIhKXFGAiIhKX/j+023QaWZke5gAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 400x300 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.figure(figsize=(4,3), dpi=100)\n",
    "plt.loglog(h_values, h1_norms, \"o-\")\n",
    "plt.xlabel(\"$h$\")\n",
    "z = plt.ylabel(r\"$\\sup_{v\\in X}\\frac{b(v,q)}{\\|q\\|}$\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "What does this mean?\n",
    "\n",
    "<details style=\"background-color: #ddd\">\n",
    "<summary>Show answer</summary>\n",
    "Since, apparently,\n",
    "$$\\sup_{v \\in X} \\frac{b (v, \\mu)}{\\|v\\|_X \\|\\mu\\|_M} =O(h)$$\n",
    "as $h\\to 0$, there cannot be a mesh-independent lower bound $c_2$ for this quantity. A discrete inf-sup condition does not hold for this discretization.\n",
    "</details>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Matt Knepley added this follow-up: (on the Firedrake Slack channel):\n",
    "\n",
    "> You can also see the instability by looking at the condition number of the Schur complement of the full system, which grows with N.\n",
    ">\n",
    ">  There's also the ASCOT approach (Automated testing of saddle point stability conditions in the [FEniCS book](https://fenicsproject.org/book/)). [Florian Wechsung] forward-ported that code to Firedrake. Additionally in the case of the Stokes problem, see this nice paper: https://www.waves.kit.edu/downloads/CRC1173_Preprint_2017-15.pdf\n",
    "\n",
    "Note that none of these approaches require manually identifying the problematic functions."
   ]
  },
  {
   "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
}