{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "cell-00",
   "metadata": {},
   "source": [
    "# L1B — $Q_1$–$P_0$ on the lid-driven cavity\n",
    "\n",
    "Mixed finite element methods — a crash course, Lecture 1.\n",
    "\n",
    "Stokes flow on the unit square, uniform $N \\times N$ quadrilateral mesh,\n",
    "velocities in continuous $Q_1$, pressures in $P_0$. Regularized lid:\n",
    "$\\boldsymbol u = (16x^2(1-x)^2,\\,0)$ on the top side, $\\boldsymbol u =\n",
    "\\boldsymbol 0$ on the remaining boundary. Signs follow the course blueprint:\n",
    "$b(\\boldsymbol v, q) = -(\\operatorname{div}\\boldsymbol v, q)$.\n",
    "\n",
    "A mesh-independent NumPy verification of the results below (kernel dimension,\n",
    "checkerboard amplitude) is contained in `verify_2d.py` in the slide source\n",
    "bundle.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "cell-01",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Firedrake available.\n"
     ]
    }
   ],
   "source": [
    "# Cell 1 [SETUP] -- Firedrake via FEM on Colab (no effect on a local installation)\n",
    "try:\n",
    "    import firedrake\n",
    "except ImportError:\n",
    "    !wget \"https://fem-on-colab.github.io/releases/firedrake-install-release-real.sh\" -O \"/tmp/firedrake-install.sh\"\n",
    "    !bash \"/tmp/firedrake-install.sh\"\n",
    "    import firedrake\n",
    "print(\"Firedrake available.\")\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "cell-02",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "setup complete\n"
     ]
    }
   ],
   "source": [
    "# Cell 2 [LECTURE] -- imports and the cavity problem\n",
    "from firedrake import *\n",
    "import numpy as np\n",
    "import scipy.sparse as sp\n",
    "import scipy.linalg as sla\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "def cavity(N):\n",
    "    mesh = UnitSquareMesh(N, N, quadrilateral=True)\n",
    "    V = VectorFunctionSpace(mesh, \"Q\", 1)\n",
    "    Q = FunctionSpace(mesh, \"DQ\", 0)\n",
    "    W = V * Q\n",
    "    u, p = TrialFunctions(W)\n",
    "    v, q = TestFunctions(W)\n",
    "    a = (inner(grad(u), grad(v)) - div(v) * p - div(u) * q) * dx\n",
    "    L = inner(Constant((0.0, 0.0)), v) * dx\n",
    "    x, y = SpatialCoordinate(mesh)\n",
    "    lid = as_vector([16 * x**2 * (1 - x)**2, 0.0])\n",
    "    bcs = [DirichletBC(W.sub(0), lid, 4),                    # top: id 4\n",
    "           DirichletBC(W.sub(0), Constant((0.0, 0.0)), (1, 2, 3))]\n",
    "    return mesh, W, a, L, bcs\n",
    "\n",
    "def centers(N):\n",
    "    xc = (np.arange(N) + 0.5) / N\n",
    "    return [(float(a), float(b)) for b in xc for a in xc]   # row-major in y\n",
    "\n",
    "def cb_grid(N):\n",
    "    i, j = np.indices((N, N))\n",
    "    return (-1.0) ** (i + j)\n",
    "\n",
    "def eval_at(f, points):\n",
    "    \"\"\"Evaluate a Function at physical points (input ordering preserved).\n",
    "\n",
    "    Recent Firedrake provides the PointEvaluator class and deprecates\n",
    "    Function.at; older installations fall back to Function.at.\n",
    "    \"\"\"\n",
    "    pts = np.asarray(points, dtype=float)\n",
    "    if pts.ndim == 1:                        # 1D mesh: (n,) -> (n, 1)\n",
    "        pts = pts[:, None]\n",
    "    try:\n",
    "        from firedrake import PointEvaluator\n",
    "    except ImportError:\n",
    "        seq = pts[:, 0].tolist() if pts.shape[1] == 1 else pts.tolist()\n",
    "        return np.asarray(f.at(seq))\n",
    "    pe = PointEvaluator(f.function_space().mesh(), pts)\n",
    "    return np.asarray(pe.evaluate(f))\n",
    "\n",
    "print(\"setup complete\")\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "cell-03",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "ConvergenceError raised; message:\n",
      "\n",
      "Nonlinear solve failed to converge after 0 nonlinear iterations.\n",
      "Reason:\n",
      "   DIVERGED_LINEAR_SOLVE\n"
     ]
    }
   ],
   "source": [
    "# Cell 3 [LECTURE] -- the direct factorization fails, as in one dimension\n",
    "N = 32\n",
    "mesh, W, a, L, bcs = cavity(N)\n",
    "w = Function(W)\n",
    "try:\n",
    "    solve(a == L, w, bcs=bcs, solver_parameters={\n",
    "        \"mat_type\": \"aij\", \"ksp_type\": \"preonly\",\n",
    "        \"pc_type\": \"lu\", \"pc_factor_mat_solver_type\": \"petsc\"})\n",
    "    print(\"solver returned without an error (unexpected)\")\n",
    "except Exception as exc:\n",
    "    print(type(exc).__name__, \"raised; message:\\n\")\n",
    "    print(exc)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "cell-04",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "matrix size: 226  zero singular values: 2  smallest four: [7.28138767e-04 7.28138767e-04 1.14224345e-16 3.88637954e-18]\n",
      "null space vs. span{constant, checkerboard}: residual = 1.3375064939976564e-13\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAApYAAAFECAYAAACd7H6zAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjExLjAsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvlcelbwAAAAlwSFlzAAAPYQAAD2EBqD+naQAAUFJJREFUeJzt3XdUVMf7P/A3ggLSBEFRQRAEC5YYu1iwYCAWNGow1lg+GjXGRE1i+2hMTEw+9phojDEmsZdYo8aCiqIIBrH3RrMAKkUEpDy/Pzjcnyu7uJgLi37fr3M8x52duTN3WB6evXtn1khEBERERERE/1IZQw+AiIiIiF4PTCyJiIiISBVMLImIiIhIFUwsiYiIiEgVTCyJiIiISBVMLImIiIhIFUwsiYiIiEgVTCyJiIiISBVMLFXUsWNHhIeHF2sfXbt2RXBwcLH28ao7d+4cfHx88PTpU0MPhej/hM2bN2PIkCHF2seuXbvQt2/fYu3jdTBt2jQsWrTote23uPq5fv06fHx8kJycrPqx/42DBw/inXfeMfQwioSJpYqCg4Px8OHDYu0jJCQECQkJxdrHv+Hv74+QkJAS6+/KlSvw8fHB48ePlbLk5GQEBwcjNze3xMZRUrSd76twbHq9xcbG4uTJk8Xax927d3HixIli7ePf2LFjB/r371+ifX755ZeYM2eORtn58+dx7dq1Eh1HSfZbXP08fvwYwcHByMrKUv3Y/0Z8fDyOHz9u6GEUiYmhB0BF89dff6FWrVqGHoZOR48eRWJiYon1l5qaiuDgYGRnZ5dYn4ZUnOf7f20u6dXSpUsX1KtXz9DD0OnOnTsICwsr0T4vXrwIS0vLEu2T6EV4xVJlKSkp+O9//4vu3btjyJAhuHjxosbzgwYNgo+PDzp06IABAwZgw4YNBY4RGhqKoUOHokuXLpg4cSLi4uKU57799luNYxZW91lTp07F7NmzC5QPHDgQq1atUh5v3LgRgYGB6NatG6ZOnar1CuymTZvQr18/dOvWDQsXLlQSkb59+yI9PR3Tpk2Dj4+Pcvk+OzsbP/74I7p3746uXbti/vz5yMzM1Dhmx44dERQUhClTpsDPzw+//fab1vN41r179zBy5EgAebcI+Pj44Msvv1Sev3nzJsaNG4euXbti/PjxBRLe9PR0zJ8/H926dUPv3r2xePHiFyZV6enpWLhwIXr06IHAwEBs375d4/nExERMnToVfn5+ePfdd7F582aN5/M/pr98+XKhY1uzZg369u2LgIAAfPfdd8jIyNB5vg8fPoSPjw98fHzg6+uLESNG4J9//ilSvy+aSyJ9HDp0CMOGDYO/vz/mzp2r8fsUGhqqvE7ffvttjB8/HjExMRrtMzIyMH/+fPTs2ROBgYH47bffICIAgFOnTmHhwoV61X1WbGwsfHx8cPPmTY3y48ePo1OnTnjy5AkA4P79+5gyZQr8/f0xcOBA/PnnnwWOdf/+fUybNg1vv/02hg4dioiICAB5n1bNnz8fcXFxyjnmx/arV69i1KhR6Ny5M95//32EhoZqHDP/NoIjR46gX79+6NChA9LT018414sWLcLBgwexa9cupc8rV64oz69atQr9+vVDr169sG7dugLtT58+jZEjR8LPzw8jR47E6dOnX9jnmTNnMGrUKPj7+2P8+PG4f/9+gTpq9KtPP/ni4uLQrVs3LF++HACQmZmJxYsXIyAgAL169cL8+fM1bovK/8j77Nmz+OCDD+Dr66vxM7l16xbGjRuHt99+G5988kmBTwiTkpIwY8YM+Pn5oXfv3li7dq3G8+np6crPo1OnThg+fDiOHTumUaewMZw6dQoDBw5EQEAAvvrqK71eC6WOkGqMjY3FwsJCJkyYINu2bZN+/fpJhQoVJC4uTqkTFhYmhw4dkqCgIFm6dKk4OjrK3LlzlefPnz8vZmZm8uWXX8qePXtk4cKF0rRpU+V5Gxsb2bRpk151n7VixQqpWLGiPH36VCm7cOGCAJBz586JiMiECRPE09NTVqxYITt37pT+/fuLq6urpKamKm3GjRsnFSpUkLlz58r27dvl448/lqlTp4qISGhoqJiZmclXX30lhw4dkmPHjomIyPvvvy9OTk7y22+/yerVq8XNzU169OhRYO5sbGxk2rRpsm/fPrl165aIiPTr108WLlyo9ZzS09Plp59+EgCyc+dOOXTokFy4cEGOHj0qAMTT01N+/vln2b59u7Ro0UJat26ttM3MzJRWrVrJW2+9JZs3b5ZNmzZJ06ZNpVevXjp+unn9NW/eXOrUqSMrV66UjRs3SkBAgOzZs0dERDIyMqRWrVrStm1b2bx5s8yfP1/Kly+vMX59xvbzzz+Lg4OD/PLLL/LXX3/JlClTZMiQITrPNzMzUw4dOiSHDh2SvXv3yowZM8TU1FTCwsL07lfXsYn0sWDBAjEzMxMPDw9ZvXq1rFy5UqpUqSKjRo1S6iQmJiqv0507d8rQoUOlYsWK8uDBA6XOiBEjpEGDBrJ+/XrZsmWLDBs2TObNmyciIsuXLxcXFxe96j6vbt268uWXX2qUDRkyRN566y0REYmNjZVq1arJiBEjZMeOHbJixQpxdnaW2bNnK/Wjo6OlSpUq0rlzZ9mwYYOsXLlSWrRoIXfu3JH4+Hj55JNPpGrVqso5RkdHy82bN8XKykoGDhwo27Ztk4kTJ0rZsmXl8OHDBeauYcOGsnbtWjl06JBkZ2fL0aNHpV27dnL37l2t53T16lVp3769vP3220qfKSkpEhAQIBYWFhIYGCjbt2+X+fPnS7ly5WTr1q1K271794qdnZ3MmjVLdu/eLV9//bVYWlrKkSNHdP6Md+3aJaampjJ69GjluB06dFCeV6tfffoZM2aMiIhcuXJFXFxc5MMPP5Tc3FzJzs6Wjh07io+Pj2zcuFH+/PNPad26tfJzFhGJjIwUAOLk5CTff/+9HDx4UBITE5VyBwcHmT9/vmzevFm8vb2lfv36kp2dLSIi2dnZ8sYbb0jz5s1l06ZNsnjxYrGyspJZs2Ypx8/OzlZ+Hvv27ZNZs2aJubm5HDx48IVjuHr1qpQvX17+85//yPbt22XUqFFSvnx5qVy5ss6fS2nExFJFxsbG0r9/f42yRo0ayccff6yzzYYNG8TJyUl5vHz5cqlbt65GncePHyv/fzaxfFHdZyUnJ4uZmZns3LlTKZsyZYo0bNhQREQuX74sZcqUURI6EZHc3FypX7++LFq0SETyElkAsnfvXp19WlhYaASSc+fOiZGRkRw/flwpO3v2rBgZGUlISIhSZmxsLOPGjSswbnd3dyWIaHPy5EkBII8ePVLK8pOoffv2KWX//POPAJD79++LiMiyZcvEw8NDsrKylDr379+XMmXKyJkzZ7T29eOPP4q1tbUkJCRoPf8ffvhB7OzsNOZj8eLFYmNjI2lpaXqPrX///jJ69GitfWg7X21GjRolAwYMKNKc6HtsouctWLBAAMipU6eUsgMHDoixsbHcvn1bZ7umTZtqvPGqVauW/Prrrxp18l/7zyeWhdV93tdffy21atVSHqenp4uNjY2sWrVKRET+85//SO/evTXa/P3332Jubq4kFUOGDJGGDRtKTk6OUufp06eSmZkpIiJLly4Vd3d3jWO8//774u3trVE2ePBgadGihfJ4wYIFWudp69atAkAjJj8vMDBQhg0bplEWEBAgDRs2lNzcXKWsb9++MnDgQOWxh4eH/PDDDxrtJk2aJL6+vjr7cnNzk7Fjx2qUPTvfavWrTz9jxoyRf/75RxwcHGTGjBnKc2vWrBEnJyfJyMhQypKSksTMzEy50JGf1OX/7PPlly9YsEApe/TokVhZWcnq1atFROS3334TS0tLefjwoVJn5cqVUr58eY2y53366acaF1N0jWHw4MHSsWNHjTI/P79XLrHkPZYqe+uttzQe+/n5aazijoqKwtKlS3Hp0iWkpKQgJSUFcXFxyMjIgJmZGVq0aIHr169j3Lhx6N27N5o3bw4LCwutfRWlrrW1Nbp27Yo1a9aga9euEBGsXbsWY8aMAZC38szExATDhw9XPk4SEdy7dw+XLl0CABw+fBjW1tbo3LmzxrF19QkA4eHhqFChAlq2bKmU1a9fHy4uLggPD4e3t7dS3rp16wLt16xZAzs7O53HL0yLFi2U/1evXh1A3ke+lSpVQlBQEFJSUuDn5wfJe4MFEYGxsTEuXbqEBg0aFDjeoUOH4OPjA3t7e43y/PMPCwtDu3btNOajW7duGDt2LK5cuYJGjRrpNbY2bdpg0qRJcHFxgb+/P+rVq1foHAN5HzP+/vvviIqKQnp6OmJiYlCpUqUizQnRv1GpUiWN13iHDh1gbGyMU6dOwcXFBQCwZcsWbN++HXfu3EFWVhZu376N69evK23atGmDb775Bjk5OejUqRNcXV11vvaLUrdfv36YNm0aIiIi0LhxY/z111/Izs5Gz549AQBBQUEoU6YMOnXqpMSC9PR0pKen4/bt23B3d8ehQ4cwdOhQlCnz/+8gK1u2bKFzEhYWhgEDBmiUdevWDevWrUNOTg6MjY0BAE5OTsoc5WvdujUOHToER0fHQvvQpnnz5jAyMlIeV69eHZGRkQCAmJgYXLt2Db///ju2bt2qnO/du3d1LtyLiorCzZs30adPH43y5+f73/arbz///PMPVq9ejVmzZuHDDz9UyoOCgpCeno4uXboox8936dIltGrVSnms7e8NoPk3vEKFCmjRogXCw8PRv39/hIWFoWXLlrC1tVXqdOvWDUOGDMH58+fRpk0bAEBERAR+/fVX3Lp1C0+ePMGdO3dgampaoK/nxxAeHl5gdwV/f39lDl8VTCxVZmVlVeBxUlISgLzVXc2aNUPHjh0xcOBA2NnZ4fLlyxgzZoySWNarVw/h4eFYsWIFRo0ahejoaIwePRrffvttgb6KUhcABgwYgH79+uHx48c4c+YMoqOj8d577wHIW7hhbW2NadOmFWiXH9jS0tJgY2NTpPlITk4uMCdAXqL7/LYO2m5Cb968eZH6e9azQT8/2OWvFE9NTUX9+vUxderUAu1q166t9XhpaWkFkspnJScno0KFChpl1tbWynP6jm3kyJGoVq0a1q9fj0WLFsHIyAjz5s1DYGCg1n4PHToEf39/fPrpp3j77bdhbW2NtWvXat36qrB+if6N53/PjYyMYGlpqcS/efPm4bvvvsPkyZPx7rvvwsLCAjNmzNC4h2zp0qX4/fffsWPHDkycOBEuLi5Yvnw5mjVrVqC/otR1dXWFt7c31qxZg8aNG2PNmjXo0aOHkrCkpqbi3XffRe/evQu0VTv+WVtb4+nTp0hPT1dinrbYZ29vDx8fnyL1l+/5hNfIyEgj9gHAiBEjULNmzULb5UtLSwOAF57/v+1X334yMjKQmZlZIOFMTU1F7dq1tf4d8/Dw0Hisa9FTYX/Dtf088x/nx/gTJ07Ax8cH48aNw8iRI2FtbY1t27Zh165dBfp6fgyFHf9VwsRSZVevXtV4fOXKFdSoUQNA3hU/ABo3+2pbbNOwYUN8//33APLewbRo0QJ+fn5ag0xR6vr7+8PMzAxbt25FaGgo2rdvj2rVqgEA3NzckJiYiFq1aqFKlSpaz83d3R137tzBo0ePNN6xPevZd/P5x7179y5SU1OVX5DMzEzcvn0bbm5uWo9RFM/3py83NzcEBQWhXbt2Gu+wC+Pu7l7gJuznj/n8diiXL19WniuKrl27omvXrgCAuXPnYvDgwejRo4fW8920aZNyo3e+Zxdk6etl55IIyItlT548Qfny5QHkLWR7+PChEv82bNiATz/9FJ988onSJv8Pdj4TExMMGzYMw4YNQ1ZWFoYOHYrRo0cXWIxW1LoA0L9/f3z55ZeYMmUKdu/erbHwzs3NDUlJSYUmcu7u7jh//rzO57X9/ri5uRX4m3D58mU4ODiospr7ZX5nq1evDhOTvD/9+iau+W3Onz+v9dMctfrVt5/WrVtj1qxZ6NOnD4yMjPD+++8DyJvvsLAwtGnTRrkaXFRXr16Fk5OTxuNevXopx9+2bZtG/fwFU/kxfsuWLejYsSO+++47pc7zbXSpUaOG1hziVcO/JCpbunQp7t27ByBvJe6mTZswePBgAHnvPFJSUhAbGwsg7wrm81cX//rrL40tKxwdHVGmTBmtyU9R6gJAuXLl0KdPH6xcuRIbN27U+IimS5cuqF69OkaNGqW8awSA3bt3K328/fbbcHR0xMcff6yssouPj8fWrVuV+pUqVcKdO3eUx507d0alSpUwc+ZMpezbb7+FiYkJunfvrnWcz+rfv3+hm+Hmf4T7bJ/6GD58OK5fv47Zs2crH5dkZWVh8eLFSElJ0dpm6NChOHv2LJYtW6aURUZGKvv3DR48GOHh4fjrr78A5CXQM2fORIcOHZSPnfWxbNky5TUCQCPR13a+VlZWuHbtmrL/2smTJwusVNTHy84lEZB3FenZnSemT58ONzc35eM+KysrjR0tfvvtN5w9e1bjGP/73/+UK1tly5aFvb29znhWlLoA8O677+LBgwcYPXo0bG1t4evrqzw3atQorF+/XuOqUlJSEhYsWKA8HjFiBFatWqXx5nLbtm2Ij48HkPf7k5CQoLESfsiQIVi1apXyBvP+/ftYtGiRXpvJh4SEwMfHR/l7os3z8VYflpaWGDBgAL766iuN/SCvXbuG1atX62zTr18/fPXVV0psysnJwYoVK1Tttyj9dO3aFZs2bcKoUaPw+++/A8ib7/v372P69OnKldKcnBwsW7ZM+Tm9yOzZs5VdSzZs2IBLly4p+5MOHDgQly9fxvr16wHk7Xgyffp0NG/eHHXr1gWQ9zq/efMmMjIyAABnz57FypUr9ep70KBB+OOPP3Djxg0AwO3bt/Hrr7/q1bZUMcB9na8tY2Nj6du3r9jb20ujRo3E1NRUhg4dqtzMnJOTIwEBAWJlZSVNmzaVChUqyDvvvKOxYOL06dPSokULcXZ2lubNm4u1tbWMHj1aOcazi3deVFebI0eOCAAxNzeXlJQUjecuXLggjRo1EhsbG2natKk4ODhI165dJSYmRqlz+vRpqV27ttjb28ubb74p1atXl6CgIOX5//3vf2JhYSHe3t7Ss2dPERE5fPiwVKlSRWrUqCEeHh5ib28vu3btKjB3+aurn/WixTsiIl27dpVKlSpJ27ZtZebMmcpClfT0dKVOQkKCAJDIyEilbMuWLeLo6CjVqlWTxo0bi52dnYwbN05j5fzz1q5dK3Z2duLq6ipeXl7StGlTiYqKUp5fuHChlC9fXurXry+VKlWSevXqyY0bN5Tn9Rnbpk2bxNXVVerUqSONGjWSChUqaCxSeP584+LixNPTU6pWrSpvvPGGODg4iK+vr7IwS99+tR2bSB8LFiyQ6tWrS/v27aVmzZpSo0YNqVixogQHByt1QkJCpGLFiuLp6Sm1a9cWd3d3adKkicbik9mzZ4uDg4M0atRIPD09xcnJSVnk9/zincLq6tK9e3cBoHWh4OzZs8XS0lI8PDykfv364ujoKIsXL9aoM336dDE3N5e6deuKq6urvPPOO/LkyRMREUlNTZWaNWuKm5ubtGvXTtavXy85OTkyYsQIMTc3l0aNGomVlZW8/fbbGjttLFiwQLy8vAqMR5/FO+Hh4WJtbS2NGjWSdu3ayeXLlzVWTef7/PPPNRaFPH78WAYMGCCmpqbSsGFDcXd3lzp16siBAwd09pWSkiI9e/ZUzqVy5coaq6HV6reo/ezYsUPMzMzk999/FxGR3bt3i5OTkzg6OkqTJk3Ezs5OPvjgA+XnlL9w5vlFmPnlgwYNkkqVKkn9+vXF1NS0wK4kP//8s1haWkq9evWkSpUqUqtWLbl48aLyfGJiotSrV08cHR3lzTffFHt7e/H399dY2KVrDFlZWdKnTx/l3B0cHKR3796v3OIdIxEtG3/RSwkODkb9+vVhamqKy5cvw87OTvkY6FnXrl1DYmIiPD09Ua5cOURERKB169bKxwRA3o3O9+7dQ40aNTTu6wsJCUGtWrXg4ODwwrraiAiOHDkCKysrvPnmm1rr3LhxAwkJCfDw8EDFihW1HuPSpUvIzMyEl5cXypUrp/H8nTt3cPv2bQBQbpbOysrC2bNnISLKHGmbu+cX6oSFhcHOzq7A/THPu3r1Ku7duwd7e3tUq1YNkZGRaNu2rfJRUVZWFo4dO4YmTZpofASVk5ODCxcuICsrC3Xq1FE+xitMZmYmzp07BysrK3h6eha4SpKSkoKLFy/CysoKdevW1Xg+OTlZr7Hl5ubiypUryMzMhKenZ4FxPXu+devWRVZWFi5evIisrCzUq1cPCQkJuH//Ppo0aVKkfrUdm+hFYmNjcffuXTRt2hTR0dG4e/cuGjRoAHNzc416aWlpuHjxIszMzFC3bl1cuXIFxsbGGl/6kJmZiYsXL8LU1BSenp5KXLx79y6ioqI0FqDpqqtLTEwMbty4AS8vL40Ymi89PR3nz59H+fLl4enpqfWew+TkZFy6dAlOTk4aH5kCeb9Tly5dwqNHj+Dm5gZnZ2cAeQvkbt68iapVq8LV1VXn3D0rMTER58+fR4sWLWBmZqbznFJSUnD16lU8fvwYjRs3RlRUFExNTTVi5s2bN5GamoqGDRsW6OPatWtwdHSEq6urXrcFxcXFITY2FrVq1dK4p/z8+fOq9luUfi5evIgHDx6gdevWyn2dFy9eREZGBmrXrq0R3x4/fox//vkH3t7eGj/fZ8tTU1Nx7do1uLm5aX2dPH78GOfPn4eFhQW8vLwK3JKQnZ2NixcvKn8jk5KSEBMTo6wZ0DWGfNevX1fuF01NTcX169c1Fh6VdkwsiYiIiEgVvMeSiIiIiFTxUollZmYmrl+/rnPPK20yMjIQExOj8dVKRET/V4kIbt68iQcPHujdJjs7G7GxscrXEBIRlTZFSixjY2Px+eefw83NDR4eHnovoZ8xYwbs7OzwxhtvwN7eHosXL36ZsRIRvfJSUlIwZ84ceHp6olatWhorqQuzZs0aODo6okGDBrCzs8PYsWO5BykRlTpFSiyDg4Nha2tbpF3gV61ahTlz5uDAgQN48OABVq9ejU8++QT79+8v8mCJiF51Fy5cQHx8PPbs2QMvLy+92pw+fRqDBg3C3Llz8fDhQ/zzzz9Ys2YN5s2bV8yjJSIqmpdevGNkZIRVq1YV+Lqq57Vq1Qru7u4aGza3b98etra22LJly8t0TUT0WnjjjTfQqVMnzJ07t9B6o0ePxtGjR3Hu3Dml7JNPPsH27dtx8+bN4h4mEZHeinXxTm5uLk6dOqXxPdFA3q75ur4dgYiINJ08eVJrHL116xYePnxooFERERVUrF/pmJqaiszMzAJ7Idrb2yMxMVFnu8zMTGXneyAvQX348CEqVqyo99fvEdHrTUSQmpqKqlWrvvZfR5mYmKg1juY/9/z+rwDjKBG9WHHE0WJNLJ/diPlZT58+LfR7PGfPnq3xFYBERLrExMQU2Kj6dVOmTBmtcRSAzljKOEpE+lIzjhZrYmllZQUbG5sC33V67969Qk9g8uTJGD9+vPI4OTkZ1atXR9SpN2Bt+XJfLE9Er5eUxzlwefM0rKysDD2UYufs7Kw1jpYpU0bju+SfpSuOononoEyxhv5CPdqXYLC+AcC2c8FvUilphp4DgPMAcA6A4omjqkeXhw8fIjk5WfkqQx8fH+zduxcTJ05U6uzZswc+Pj46j2FqalrgK/8AwNrSGNZWhguIRFT6vI4f66alpeHu3btwdXWFiYkJfHx88NNPPyEnJ0e5Qrl79240a9ZM59eQ6oqjKGMClCn4NXIlxeAx3IDnns/gcwBwHgDOwTPUjKNF+kA9PT0d169fx/Xr1wEA9+/fx/Xr15GQ8P8z7u+//x6NGjVSHk+dOhVHjhzB9OnTcfLkSYwZMwYxMTGYMGGCSqdARPTqyM3NVeLo06dPkZSUhOvXryM2Nlaps3//fnh4eChlY8aMQW5uLgYPHozw8HDMmzcPmzZtwvTp0w11GkREWhUpsYyMjISfnx/8/Pzg7u6OpUuXws/PDwsWLFDq2NnZwc3NTXnctGlT7Nu3DydPnsT777+PuLg4BAcHo2bNmuqdBRHRKyI5OVmJo0+fPsXhw4fh5+eHsWPHKnUsLS3h7u6OsmXzrqg4ODjg6NGjyMrKwpAhQ7B7925s3boV/v7+hjoNIiKtXnofy5KUkpICGxsbPLrauNRcNiYiw0pJzYatZwSSk5NhbW1t6OGUevlxFK5+Bv0IMOdYvMH6BgBj70oG7R8w/BwAnAeAcwAUTxx9vffoICIiIqISw8SSiIiIiFTBxJKIiIiIVMHEkoiIiIhUwcSSiIiIiFTBxJKIiIiIVMHEkoiIiIhUwcSSiIiIiFTBxJKIiIiIVMHEkoiIiIhUwcSSiIiIiFTBxJKIiIiIVMHEkoiIiIhUwcSSiIiIiFTBxJKIiIiIVMHEkoiIiIhUwcSSiIiIiFTBxJKIiIiIVMHEkoiIiIhUwcSSiIiIiFTBxJKIiIiIVMHEkoiIiIhUwcSSiIiIiFTBxJKIiIiIVMHEkoiIiIhUwcSSiIiIiFRhYugBEBFRyXm0LwHWVoYL/cbelQzWNwDkHIs3aP+A4ecA4DwAnAMAQG6W6ofkFUsiIiIiUgUTSyIiIiJSBRNLIiIiIlIFE0siIiIiUgUTSyIiIiJSBRNLIiIiIlIFE0siIiIiUgUTSyIiIiJSBRNLIiIiIlIFE0siIiIiUgUTSyIiIiJSBRNLIiIiIlIFE0siIiIiUgUTSyIiIiJSRZETy9OnT2PQoEFo164dhg8fjuvXr7+wzYYNG/Duu+/Cx8cH/fr1w65du15qsEREr4O4uDh8+OGH8PHxQd++fRESEvLCNkeOHMHgwYPRvn17vPPOO1ixYgVycnJKYLRERPorUmJ5/vx5tG7dGpaWlpg0aRLS0tLQokULxMbG6mwzf/58DBkyBO3atcMXX3yBBg0aICAgAGvXrv3XgycietUkJyejVatWuH37Nj777DPUqFEDHTp0wNGjR3W2OXDgANq3bw9nZ2fMmDEDXbp0waefforJkyeX4MiJiF7MSERE38qBgYGIi4tT3l3n5OSgVq1a6NatGxYsWKC1Tdu2beHu7o6VK1cqZX5+frC2tsbGjRv16jclJQU2NjZ4dLUxrK1M9B0uEb3GUlKzYesZgeTkZFhbWxt6OHqbPXs25s2bh7i4OJiamgIAevTogeTkZBw6dEhrm48++ghHjx5FZGSkUjZp0iRs2bIFV69e1avf0hJHjb0rGaxvAMg5Fm/Q/gHDzwHAeQA4BwCA3Czg9t+qxtEiXbEMCgpCt27dlMfGxsbo0qULDhw4oLNNkyZNcPr0aTx58gQA8PDhQ1y6dAnNmjV7ySETEb26goKC4OvrqySVABAQEICjR48iMzNTa5smTZogKioKcXFxAICnT5/i5MmTjKNEVOro/bY1LS0NDx48QNWqVTXKq1atiqioKJ3tvvvuO0yYMAFOTk5wcXHBzZs3MWHCBEyYMEFnm8zMTI0Am5KSou8wiYhKtaioKDRs2FCjrGrVqsjJyUFcXBzc3NwKtBk0aBCSk5Ph5eUFV1dXxMXFwc/PDz///LPOfhhHicgQ9L5imZWVBQAa77IBwNzcXHlOm/Xr1+OPP/7AjBkzMG/ePHz++eeYN28e9u3bp7PN7NmzYWNjo/xzdnbWd5hERKVaVlaW1jia/5w2J06cwPTp0zF69GjMnz8fX3/9Nfbu3VtoYsk4SkSGoPcVSysrK5QtWxYPHz7UKH/w4AEqVqyos9348ePx8ccfY9y4cQCADh06IDo6Gp9++ineeustrW0mT56M8ePHK49TUlIYFInotWBnZ6c1jgLQGUv/+9//ol27dvjmm28A5MXR3NxcfPzxxxgxYoSSmD6LcZSIDEHvK5bGxsZo0KABTp48qVEeFhaGRo0aaW2Tk5ODlJQUVKtWTaO8atWqBQLrs0xNTWFtba3xj4jodfDmm29qjaNOTk6wt7fX2ubhw4da42hmZibS0tK0tmEcJSJDKNLinWHDhmHz5s24ePEiACAkJARBQUEYNmyYUueXX36Br68vgLxktFWrVli5ciVSU1MBAImJiVi7di3atGmj1jkQEb0yhg0bhsjISOzcuRMAEB0djd9++00jjh45cgQtWrTA3bt3AeTtrrF161ZER0cDADIyMvDzzz+jdu3aOpNRIiJDKNKeEx988AHOnz+PRo0awdXVFVFRUZg8eTJ69Oih1ImNjdV4N75ixQr0798fTk5OcHV1xY0bN9CyZUssWrRItZMgInpV5Me/wMBAODk5ISYmBj169MCUKVOUOg8fPkRYWJiy+Oarr75CdHQ0atWqBQ8PD8TFxcHJyQnr16831GkQEWlVpH0s8yUkJCA2Nhaurq6wtbXVeC42Nhb37t1DkyZNNMrv3buH+/fvo1q1akV+h11a9l8jotLjVd3HMl9qaiquX7+OypUrF9ht49GjR7hy5QoaNWqksdDn0aNHiI6Ohr29PapWrQojIyO9+ystcdTQ+/Zx78I8nAfOAYBi2cfypaKLg4MDHBwctD7n5OQEJyenAuWOjo5wdHR8me6IiF47VlZWOu9Pt7W1RYsWLbSWP/9mnoioNCnyd4UTEREREWnDxJKIiIiIVMHEkoiIiIhUwcSSiIiIiFTBxJKIiIiIVMHEkoiIiIhUwU0hiYj+D7Ht7ACUKWuw/g29d6DB9w2E4ecA4DwAnAMgfz9gdY/JK5ZEREREpAomlkRERESkCiaWRERERKQKJpZEREREpAomlkRERESkCiaWRERERKQKJpZEREREpAomlkRERESkCiaWRERERKQKJpZEREREpAomlkRERESkCiaWRERERKQKJpZEREREpAomlkRERESkCiaWRERERKQKJpZEREREpAomlkRERESkCiaWRERERKQKJpZEREREpAomlkRERESkCiaWRERERKQKJpZEREREpAomlkRERESkCiaWRERERKQKJpZEREREpAomlkRERESkCiaWRERERKQKJpZEREREpAomlkRERESkCiaWRERERKQKJpZEREREpAomlkRERESkChNDD4CIiErOo30JsLYyXOg39q5ksL4BIOdYvEH7Bww/BwDnAeAcAABys1Q/ZJGjS05ODoKCghAVFQUPDw+0a9cORkZGL2x3//59HDp0CEZGRvD19YWdnd1LDZiI6HVw4sQJnD9/Ho6OjvD19YWpqekL26SmpiIoKAgpKSnw8fFB9erVS2CkRET6K9JH4WlpaWjbti1GjhyJw4cPo1+/fujWrRuysgrPeH/66Se4u7tj9erV+Pvvv9GmTRtERkb+q4ETEb2q3n//fXTt2hUHDx7ExIkT0bhxYyQmJhba5sCBA3Bzc8PcuXNx9OhRdOnSBZs3by6hERMR6adIVyy/++47REVF4ezZs7Czs0NUVBTq1auHX375BaNGjdLaJjg4GKNHj8aff/6Jnj17AgAePHiA+HjDX4ImIippW7ZswZo1a3D69Gl4eXkhLS0NTZo0wdSpU7Fs2TKtbeLi4vDOO+9g/Pjx+OKLLwAAT58+xenTp0tu4EREeijSFcsNGzYgMDBQ+RjbxcUFXbt2xYYNG3S2mTdvHlq3bq0klQBQsWJF1KlT5yWHTET06tqwYQPatm0LLy8vAICFhQXef/99bNy4ESKitc2yZctQrlw5TJ06VSkrV64cmjVrViJjJiLSl95XLLOysnDt2rUCCWGdOnUQFBSks92xY8fw0Ucf4cqVKwgJCUGlSpXQtm1b2NjY6GyTmZmJzMxM5XFKSoq+wyQiKtUuXryIdu3aaZTVqVMHSUlJuHv3LqpWrVqgzbFjx9CmTRskJiZi3759sLCwQMuWLVGtWjWd/TCOEpEh6H3FMi0tDSKCChUqaJTb2trqDFgigocPHyI4OBjdunVDSEgIZs2aBQ8PD5w4cUJnX7Nnz4aNjY3yz9nZWd9hEhGVaqmpqVrjKKA7+Xvw4AHi4uLQpk0bBAUF4ZdffoGHhwdWr16tsx/GUSIyBL0TSzMzMwB5QfFZKSkpKF++vNY2RkZGMDc3x40bNxAZGYmVK1ciLCwMbdu2xejRo3X2NXnyZCQnJyv/YmJi9B0mEVGpZm5urjWOAtAZS8uXL48zZ85g7969+OOPP/D3339jypQpGDlyJNLT07W2YRwlIkMoUmLp7OyMW7duaZTfvHkTHh4eOtt5enrC29sbFhYWSpmvry/Onz+v834iU1NTWFtba/wjInodeHh4aI2j5ubmOj/a9vT0hLu7O9zd3ZUyX19fPHnyBDdv3tTahnGUiAyhSIt3unfvjk2bNin37SQnJ2Pnzp3o3r27Uic8PBzLly9XHvfs2RNnz55FTk6OUhYZGQk3Nze99r8kInqddO/eHQcOHMC9e/cAALm5uVi7di26dOkCY2NjAMCtW7fwww8/KFc2e/bsiejoaDx48EA5TmRkJExMTODi4lLyJ0FEpIOR6LpsqMX9+/fRvHlzODs7o3Pnzti6dSuysrJw/PhxWFlZAQC++OILLFy4EElJSQDyPjpv3bo17Ozs4O/vjwsXLmDTpk34888/4e/vr1e/KSkpsLGxwaOrjQ36jRFEVHqkpGbD1jMCycnJr9TVuKysLHTs2BEJCQl47733cPz4cURERCA0NBQ1a9YEAGzbtg09e/bErVu34OrqCgDo1asXLl++jP79+yMhIQHLly/HzJkzMWHCBL36LS1x1NDfNMJvW8nDeeAcAMj75p3bf6saR4t0xbJy5cqIjIxE7969kZSUhBEjRiAsLExJKgGgWbNmGDFihPLYysoKJ06cQL9+/RAfH48GDRrgwoULeieVRESvk7Jly+LAgQP47LPPkJycjI4dO+LcuXNKUgkAbm5uGDNmjEag37x5M2bOnIlHjx7BwcEBhw8f1jupJCIqKUW6YmkopeWdNhGVHq/qFUtDKS1x1NBXaHiVKg/ngXMAwPBXLImIiIiIdGFiSURERESqYGJJRERERKpgYklEREREqmBiSURERESqYGJJRERERKpgYklEREREqmBiSURERESqYGJJRERERKpgYklEREREqmBiSURERESqYGJJRERERKpgYklEREREqmBiSURERESqYGJJRERERKpgYklEREREqmBiSURERESqYGJJRERERKpgYklEREREqjAx9ACIiKjk2HZ2AMqUNVj/OcfiDdY3ABh7VzJo/4Dh5wDgPACcAwBISc2Grae6x+QVSyIiIiJSBRNLIiIiIlIFE0siIiIiUgUTSyIiIiJSBRNLIiIiIlIFE0siIiIiUgUTSyIiIiJSBRNLIiIiIlIFE0siIiIiUgUTSyIiIiJSBRNLIiIiIlIFE0siIiIiUgUTSyIiIiJSBRNLIiIiIlIFE0siIiIiUgUTSyIiIiJSBRNLIiIiIlIFE0siIiIiUgUTSyIiIiJSBRNLIiIiIlLFSyeW6enpRW4jIkhKSsKTJ09etlsiotdGRkbGS7VLSUlBamqqyqMhIvr3ipxYLl68GJUqVYKVlRWcnJywZs0avdvOnj0btra2GDRoUFG7JSJ6bezevRs1a9aEpaUlbG1tMXPmzCK1rVChAho1alSMIyQiejlFSiy3bt2KiRMn4pdffkF6ejq+/PJLDB48GMeOHXth29DQUPzyyy9o1arVSw+WiOhVd/nyZfTs2RMffPABnjx5gm3btmHOnDlYunTpC9veuXMHI0eORPfu3UtgpERERVekxHLhwoXo0aMHunfvjrJly2Lo0KFo3rw5Fi9eXGi7pKQk9O/fH7/++itsbW3/1YCJiF5lS5cuhYuLCyZOnIhy5cqhXbt2GDx4MBYtWlRou9zcXAwYMACffvop6tWrV0KjJSIqGr0TSxHByZMn0aZNG43ydu3aISwsrNC2w4cPR58+feDj4/NSgyQiel2EhYWhbdu2GmXt2rXDlStXkJSUpLPd119/jXLlymHs2LHFPEIiopdnom/F1NRUpKenw97eXqO8UqVKiI+P19lu6dKluHnzJtauXav3oDIzM5GZmak8TklJ0bstEVFpFh8frzWOAkBCQgIqVKhQoE1ISAiWLFmCU6dOwcjISK9+GEeJyBCKvHgnNzdX43F2drbOQHf16lVMnjwZS5YswZMnT5CUlITs7GxkZWUhKSmpwLHyzZ49GzY2Nso/Z2fnog6TiKjU0hZHAWiNpWlpaejXrx++/fZbmJubIykpCZmZmcjNzUVSUhKysrK09sE4SkSGoHdiaWVlBSsrK9y/f1+jPD4+HlWrVtXa5tq1awAAPz8/uLq6wtXVFUFBQdizZw9cXV0RGxurtd3kyZORnJys/IuJidF3mEREpVq1atW0xlEjIyM4OjoWqP/gwQOkpKRg3LhxShxdvHgxoqKi4Orqin379mnth3GUiAxB74/CjYyM4O3tjaCgIHzyySdK+f79+9G6dWvlcUZGBjIzM2FjY4MuXboUuGeoa9euMDMzw+bNm3X2ZWpqClNT0yKcBhHRq6F169ZYvXo1RES5Qrl//3688cYbsLS0BABkZWUhLS0N1tbWqF69eoE4Om3aNKxfvx7Xr1/X2Q/jKBEZQpE+Cv/888/x999/4/vvv8etW7cwffp0XLp0CePHj1fqfPvtt3BxcVF9oEREr4MxY8YoVyBv3ryJ33//HatXr8bkyZOVOrt27YKtrS2io6MNOFIioqIrUmLp4+ODzZs3Y9WqVWjRogUOHjyIPXv2aGx9YWZmBhsbG53HsLS0hIWFxcuPmIjoFebk5ISDBw/iwoULaNmyJebPn4/ly5ejT58+Sp2yZcvCxsYGZcpoD9Hm5uawtrYuqSETEenNSETE0IN4kZSUFNjY2ODR1cawttL703sieo2lpGbD1jMCycnJTLL0kB9H4eoHlClrsHHkHNO9i0hJMPauZND+AcPPAcB5ADgHQPHEUWZpRET/hzzal2DQN+iG/mNu6D/kgOHnAOA8AJwDAECu9l0l/o0ibzdERERERKQNE0siIiIiUgUTSyIiIiJSBRNLIiIiIlIFE0siIiIiUgUTSyIiIiJSBRNLIiIiIlIFE0siIiIiUgUTSyIiIiJSBRNLIiIiIlIFE0siIiIiUgUTSyIiIiJSBRNLIiIiIlIFE0siIiIiUgUTSyIiIiJSBRNLIiIiIlIFE0siIiIiUgUTSyIiIiJSBRNLIiIiIlIFE0siIiIiUgUTSyIiIiJSBRNLIiIiIlIFE0siIiIiUgUTSyIiIiJSBRNLIiIiIlIFE0siIiIiUgUTSyIiIiJSBRNLIiIiIlIFE0siIiIiUgUTSyIiIiJSBRNLIiIiIlIFE0siIiIiUoWJoQdAREQlx7azA1CmrMH6zzkWb7C+AcDYu5JB+wcMPwcA5wHgHABASmo2bD3VPSavWBIRERGRKphYEhEREZEqmFgSERERkSqYWBIRERGRKphYEhEREZEqmFgSERERkSqYWBIRERGRKoq8j+Xt27exZMkSREVFwcPDAx999BEqVdK9F1R2djY2bdqE4OBgZGVloWnTphgyZAhMTU3/1cCJiF5VycnJWLx4Mc6fPw9HR0eMHDkSderUKbTNvn37sGfPHjx48AB169bFyJEjYWtrW0IjJiLST5GuWN66dQuNGzfGzZs34efnh/DwcDRp0gSJiYk623h7e+Ovv/7CG2+8gRYtWuD7779H27ZtkZmZ+a8HT0T0qnny5Am8vb2xe/dudO7cGcnJyWjSpAkiIyN1thk2bBjmz58PFxcXdOrUCfv374eXlxfu3LlTgiMnInoxIxERfSsPGTIEZ8+excmTJ1GmTBlkZmaiZs2aGDhwIL755hutbe7evYsqVaooj2NjY+Hs7IwtW7agZ8+eevWbkpICGxsbPLraGNZW/LIgIsr/xogIJCcnw9ra2tDD0duiRYswffp0xMTEKOP29fVF2bJlsXv3bq1tno+jmZmZcHNzw/DhwzFz5ky9+s2Po3D14zfvGJih5wDgPACcA6B44miRrlju2bMHPXv2RJkyec1MTU3RrVs37NmzR2ebZ4MhADg4OKBs2bJITk5+ieESEb3a9uzZA19fX40g3qdPHxw4cABZWVla2zwfR01NTWFvb884SkSljt6JZXp6Ou7fvw9nZ2eNcmdnZ9y6dUvvDpcsWQIjIyO0b99eZ53MzEykpKRo/CMieh3cunVLaxzNyspCXFycXsc4dOgQzp49Cz8/P511GEeJyBD0Tizz74ksX768RrmlpSUyMjL0OsbBgwcxadIkzJ07Fy4uLjrrzZ49GzY2Nsq/54MwEdGrKjMzU2scBaBXLL1x4wb69u2LoUOHFppYMo4SkSHonVhaWlrC2NgYjx490ih/8OCBXisTjx49iu7du2PKlCkYO3ZsoXUnT56M5ORk5V9MTIy+wyQiKtVsbGy0xlEAL4ylt2/fRocOHdC2bVssW7as0LqMo0RkCHqvhDExMUHdunVx5swZjfLTp0+jQYMGhbY9duwY3n77bUyYMAEzZsx4YV+mpqbcjoiIXksNGzbUGkcrVaqEypUr62wXFRWF9u3bo2nTpli3bh1MTAoP34yjRGQIRVq8M3DgQGzcuBHR0dEAgHPnzmHv3r0YOHCgUmf9+vXo37+/8jg0NBR+fn4YP3683qsXiYheVwMHDkRoaChCQkIAAImJiVi5cqVGHA0LC0OPHj0QH5+3YjQmJgY+Pj5o3Lgx1q9f/8KkkojIUIoUnT7++GOcOHECDRs2RIMGDRAREYGBAwdqJJKXL1/Grl27lMfdu3cHAJw5cwY9evRQyvv27Yu+ffv+y+ETEb1afH19MWXKFHTu3BlNmjTBpUuXUK9ePXzxxRdKnbt372L79u1YuHAhAGDo0KGIioqCl5cXevfurdRr3rw5Jk+eXMJnQESkW5ESy7Jly+LPP//EhQsXEB0djZo1a8LDw0OjTt++fdG0aVPl8a+//oqcnJwCx6pdu/ZLDpmI6NU2a9Ys/Oc//8HFixfh6OiIRo0aaTzfvHlzbN26VflWs//+978YM2ZMgeM8vw0REZGhvdTnKV5eXvDy8tL6XO3atTWSxm7dur3cyIiIXmMuLi46d8eoUqWKxic8bdu2LaFRERH9O0W6x5KIiIiISBcmlkRERESkCiaWRERERKQKJpZEREREpAomlkRERESkCiaWRERERKQKJpZEREREpAomlkRERESkCiaWRERERKQKJpZEREREpAomlkRERESkCiaWRERERKQKJpZEREREpAoTQw+AiIhKzqN9CbC2MlzoN/auZLC+ASDnWLxB+wcMPwcA5wHgHAAAcrNUPySvWBIRERGRKphYEhEREZEqmFgSERERkSqYWBIRERGRKphYEhEREZEqmFgSERERkSqYWBIRERGRKphYEhEREZEqmFgSERERkSqYWBIRERGRKphYEhEREZEqmFgSERERkSqYWBIRERGRKphYEhEREZEqmFgSERERkSqYWBIRERGRKphYEhEREZEqmFgSERERkSqYWBIRERGRKphYEhEREZEqmFgSERERkSqYWBIRERGRKphYEhEREZEqmFgSERERkSqYWBIRERGRKphYEhEREZEqTF6m0dmzZxEVFQUPDw/Url272NoQEb2ubt26hQsXLsDR0RGNGzeGkZFRsbQhIipJRbpi+fTpU/To0QPt27fHwoUL0axZMwwbNgwiomobIqLX2aRJk1CvXj3Mnz8fXbt2hY+PD1JTU1VvQ0RU0oqUWC5cuBDHjh3DmTNnEBQUhNDQUKxduxarV69WtQ0R0etq3759mDNnDg4cOICDBw/iwoULiIqKwhdffKFqGyIiQyhSYrlq1SoEBgbCyckJAODl5QV/f3+sWrVK1TZERK+rVatWoWXLlmjZsiUAoGLFihg6dOgL42hR2xARGYLe91hmZ2fj0qVLGDt2rEZ5gwYN8NNPP6nWBgAyMzORmZmpPE5OTgYApDzO0Xe4RPSay48Hr9ptNWfPnoW3t7dGWYMGDZCQkIB79+7B0dFRlTalNo7mZhm0+5TUbIP2D8DgcwBwHgDOQV7/eXOgZhzVO7F8/PgxcnJyYGtrq1FesWJFJCUlqdYGAGbPno2ZM2cWKHd587S+wyWi/yMePHgAGxsbQw9Db8nJyVpjIgAkJSVpTRJfpg3jqHa2noYeQenAeeAcPEvNOKp3YmlqagoAePLkiUb548ePYWZmplobAJg8eTLGjx+vPE5KSoKLiwuio6NfqT8guqSkpMDZ2RkxMTGwtrY29HD+tdftfIDX75xet/MB8pKt6tWrw87OztBDKRJTU1OtMRFAobG0qG0YR189r9s58XxKv+KIo3onlubm5nB0dER0dLRGeXR0NNzc3FRrA+QF0fyk9Fk2NjavzQ8TAKytrXk+pdzrdk6v2/kAQJkyr9Z2vG5ublpjYrly5VCtWjXV2jCOvrpet3Pi+ZR+asbRIh3J398fW7ZsQW5uLgAgIyMDO3fuhL+/v1Ln4sWL2LFjR5HaEBH9X+Hv74/9+/cr9zwCwKZNm9CpUyeULVsWAHDnzh1s3rxZuUqpTxsiotKgSInl9OnTERsbi169emH58uXo0qULTExMND5u2bhxIwYNGlSkNkRE/1f85z//gaurK9566y38/PPPGDRoEE6cOIFvvvlGqRMeHo4+ffogPj5e7zZERKVBkRJLV1dXnDp1CrVr18bhw4fRpk0bnDx5UrmJHADq1q2LgICAIrV5EVNTU8yYMUPrxzqvIp5P6fe6ndPrdj7Aq3tO5ubmCAkJQUBAAIKDg+Hg4ICIiAg0bNhQqVOtWjX06tULFhYWerd5kVd1vnR53c4HeP3OiedT+hXHORnJq7ZXBxERERGVSq/WXe9EREREVGoxsSQiIiIiVTCxJCIiIiJV6L2PZXHKzMzEsWPH8PjxYzRr1kzrt0io0aYknTt3DtevX4erqysaNWr0wvpPnjxBREQEnjx5gnr16uncm85QEhMTERoaCjMzM7Ru3Rrm5uZ6tz148CDi4+MREBBQpHbFKTs7G8ePH8ejR4/w5ptvwtnZWa92CQkJCA8Ph52dHZo3b16q9lC8cuUKLl26hKpVq6Jp06YwMjIqtL6I4NSpU4iNjYWdnR2aNWtWqm5Kf/r0Kfbv34/c3Fx069ZNrzZJSUk4duwYjI2N0bp1a1haWhbzKEuP3NxchIeH4969e/Dy8oKHh0extClJt2/fxunTp2Fvb4+WLVvC2Ni40PrZ2dk4ffo07t27B09PT3h6lq6vVnny5AlCQkLw9OlTtGrVqkibUkdGRuLKlSto27YtqlatWoyjLJqIiAhER0fD09MTXl5eerV5/PgxQkNDUaZMGXh7exf6hSkl7c6dOzh58iRsbGzQqlUrlCtX7oVtLl++jGvXrsHCwgKNGjUq8C1ZhiQiCA4Oxr1799CrVy+9tifLzMxESEgI0tLS0Lx5c1SuXLnInRrUtWvXxNXVVWrVqiVt27aV8uXLy4oVK1RvU1Kys7Olf//+UqFCBencubNUrFhRunfvLpmZmTrbzJ07V5ycnKR169bSuXNnMTc3l08//bQER124jRs3ioWFhXh7e4uXl5dUrVpVzpw5o1fbgwcPirm5uQCQmJiYYh6pfuLi4qRu3bpSo0YN6dChg5ibm8vcuXNf2G7mzJlSvnx56dixo3Ts2FHatGkjSUlJJTDiF/vwww/F0tJSfH19pXLlyuLj4yOPHz/WWf/+/fvSsGFDqVatmgQEBEjdunXF0dFRwsLCSnDUuk2bNk2qVasmLi4u4uLiolebvXv3io2NjTRv3lwaNWokFStWlJCQkOIdaCmRnJwsrVq1kqpVq4qvr69YWFjIhAkTVG9Tkr7++mvl96169erSoEEDuX//vs76W7ZskZo1a0rjxo2lS5cuYmNjIz179pSMjIwSHLVu//zzj1SuXFkaNGggrVq1EisrK9m2bZtebW/fvi2VK1cWALJz585iHql+0tPTxc/PTxwcHKRz585ibW0tQ4cOldzc3ELbrV27VmxsbKRFixbStWtXqV+/vly6dKmERl24pUuXSvny5cXHx0dq1qwpNWvWlFu3bumsn5mZKQEBAWJjYyNdu3aVFi1aiKWlpaxZs6bkBl2IFStWiIeHh9SsWVMAyKNHj17Y5sqVK+Li4iK1a9dW8qvff/+9SP0aPLFs166ddO7cWbKzs0Uk7wdbrlw5iYqKUrVNSVm+fLlYWlrK1atXRUQkKipK7OzsCk1cVq5cqfEDP3bsmBgZGcmuXbuKe7gvFB8fL5aWljJnzhwREcnNzZWePXvKG2+88cK2CQkJ4uLiIt98802pSix79eolzZo1U/7gbNq0SYyMjOTs2bM62yxbtkzMzMzkxIkTSllYWJjcvXu32Mf7Ilu3bhUTExM5deqUiOTNu5OTk0yaNElnm88++0ycnZ2V5DM3N1f8/Pykbdu2JTLmF5k/f77cu3dPvvrqK70Sy7S0NHFwcJDPP/9cKRs2bJi4urpKVlZWMY60dBg3bpy4u7srceT48eNiZGQkf//9t6ptSkp4eLgAkD179ohI3s+3QYMGMmDAAJ1t/vzzT42/ATExMWJnZydff/11sY/3RXJzc6V27drSv39/pWzGjBlSoUKFF745zcrKkpYtW8qcOXNKVWL51VdfiaOjo9y5c0dERM6fPy+mpqayevVqnW1CQ0OlTJkyGonK7du35fTp08U+3he5evWqmJiYyKpVq0RE5OnTp9KmTRvx8/PT2Wbt2rVSpkwZjeTz888/lwoVKrwwwS4JK1askKtXr8rOnTv1Tiy9vb3F399fcnJyRERk8eLFYmZmJrGxsXr3a9DEMiYmRgBoJFBPnz6VChUq6EzEXqZNSWrXrp289957GmUjRozQKxF7VpUqVUpFQFy+fLmYmZlpXP06cuSIAJDz58/rbJebmytdunSRWbNmyf79+0tNYpmSkiImJiYF3oG5uLgUmoi5urrK6NGji3t4L6VXr17i6+urUTZlyhRxcnLS2WbcuHHy5ptvapSNGjVKWrRoUSxjfFn6JpZbt24VIyMjjUT/4sWLAkAOHz5cjCMsHezt7QvEi1atWhWaiL1Mm5Iybtw4qVOnjkbZkiVLxMzMTNLT0/U+jr+/v/Tp00ft4RXZyZMnBYBEREQoZQ8ePBATE5NCEzERkcmTJ0vv3r0lISGhVCWWtWrVkk8++USjrHv37oUmYgEBAdKqVaviHtpL+eqrr6RSpUpKQiWS92mdkZGRxMfHa22zcuVKMTc31/hEctmyZWJhYaFc+CoN9E0sb926JQBk7969SllGRoZYWVnJwoUL9e7PoDeInTt3DgBQr149paxs2bKoVauW8pwabUrSuXPnNMYGAPXr18eFCxcgem4ZeuHCBdy7d6/AcQzh3LlzqFGjhrJRM5B3PvnP6bJgwQIkJydj0qRJxT7Gorh8+TKys7MLzG29evV0nk90dDRu376Nzp0748qVK9i+fTsiIyP1/nkWN12vudjYWCQlJWltM3HiROTm5mLYsGH4/fffMW3aNPz999+YN29eCYxYfefOnYO9vb3GvdZ16tRB2bJlS0VcKE53795FYmKi1teArnN/mTYlSddrOiMjA9evX9frGCkpKTh58mSpiaOA5t8tOzs7VKtWrdD5DgoKwurVq7Fs2bJiH2NRZGZm4urVq0V+/Rw5cgSdO3fGnTt3sH37doSGhuLp06fFPVy9nDt3Dl5eXhr3zdevXx8iggsXLmht07dvX3To0AHdu3fHr7/+ijlz5uC7777D0qVLX3g/cGmk7XVqamoKT0/PIsUFgy7eyf/e2+dvYK5YsaLOP4gv06YkJScnax1bVlYW0tLSXriYIC0tDQMGDIC3tze6du1anEPVi7bzqVChAoyNjXXOd0REBL799luEh4eXul+uwl4/N27c0Nom/2v11qxZg1OnTqFu3boIDw9HjRo1sHv3boPfqK3rNQfkLWapUKFCgTa2trZo2bIldu3ahYSEBFy7dg116tSBk5NTSQxZddrmAMg7z9IQF4rT6xpHn19I9Oxr+kVEBMOHD4eZmRk+/PDD4hhikSQnJ8PCwqLAQpDC5js+Ph6DBg3CH3/8ATs7OyQmJpbASPWTmpoKESnS6ycnJwePHj1CREQEVq5cifr16+Py5csQEWzfvl3vhT/F5UVxVJty5cqhVatW+Omnn/Dnn38iMTERtra2pW7RmL7UigsGvWKZvwL18ePHGuWPHz/WuUrsZdqUJFNTU61jA/DC8aWnp6N79+7IysrC1q1bS8WKY23nk5GRgZycHJ3nM3LkSHTo0AEnTpzA+vXrcfjwYQDAjh07EBkZWdxDLtTLvH7yy+/fv49Lly5hx44duHr1Ku7du4fp06cX74D18DKvuU8++QSHDx/GxYsXsWPHDly4cAHW1tYaX8f6KtE2B0DpiQvFiXG0oDFjxuDgwYPYvXt3kVZeFxdTU1Okp6cjNzdXo7yw+Z4yZQqcnZ2RkJCA9evXY+vWrQCAo0eP4uDBg8U+5sK8zOvH2NgYJiYmOHXqFCIjI7Fz505cunQJbm5u+OCDD4p9zC/yMq+5n376Cd999x2OHj2KXbt2ISwsDL169YKfn1+peINWVGrFBYNmLu7u7gDyPmp8VnR0NNzc3FRrU5Lc3d0LjC0qKgrOzs4wMdF9gTgjIwPdu3fHvXv3cPDgQdjb2xf3UPXi7u6O2NhYjYB4+/ZtANA53y1btgQAbNu2Ddu2bcOxY8cAAHv37jX4x2y6Xj9RUVE6z6dGjRooU6YMunfvrmzVYG1tjU6dOiEiIqJ4B6wHXa85CwsLndtEBAcHo0uXLsotDmXKlEGvXr1w+vTpVzIguru7IyEhARkZGUpZYmIinjx5UiriQnFycnJCuXLlivSafpk2JUnXaxrQHXfyjR07Fhs3bkRQUJBy246hubu7Izc3F7GxsUpZdnY27ty5o/N86tSpA1dXVyWO7tmzBwAQGhqKkJCQEhm3LlZWVnBwcCjy68fNzQ3t27dXPuUxMTFBQEBAqY6jgO7XXHBwMJo1awYXFxelrE+fPkhKSsLp06eLbazFRbX8qui3gaonNzdXnJ2dZeLEiUrZiRMnBIDGNiH79++X48ePF6mNoXz++efi4uKirDjOysqSOnXqyAcffKDUuXr1qqxbt05ZNZaeni6+vr5St25duXfvnkHGrcvZs2cFgBw4cEAp++KLL8TOzk65Yfnp06eybt06uXnzptZjlKbFOyIiDRs2lPfff195fO3aNTEyMtLY+uPIkSNy6NAh5XHHjh1l+PDhGsdp3bq1vPvuu8U+3heZM2eO2NraSkpKiojk/Y60atVKevfurdSJioqSdevWKQsfOnbsKF26dNE4zpdffimWlpalYjVjvsIW72zYsEEuX74sInmL+kxMTGT9+vXK8z/88IOYm5vrtRLyVdelSxfp1KmT8vjBgwdiYWEhP/zwg1J28uRJjUWP+rQxlHXr1omJiYnGStTevXtrLC5LSEiQdevWycOHD5Wyjz76SCpWrCiRkZElOdwXSktLEysrK40Fptu2bRMjIyNlBxERkR07dii7OzyvtC3eGTx4sLz55pvKYpe0tDSpXLmyfPHFF0qds2fPypYtW5THEyZMkKZNm2ocZ8yYMeLp6Vkygy5E/t+pc+fOKWUffPCB1KxZU3mckpIi69atUxYJTpgwQVxdXeXp06dKnR07dggAuXbtWskN/gUKW7yzd+9eCQ0NFRGRnJwcqVq1qsZC1pCQEAGgsSPKixh8u6HNmzeLiYmJTJw4URYsWCDOzs4SGBioUad58+YaZfq0MZT8LXY6dOggS5YsEX9/f6lcubJGUrV48WIBoGyD0rNnTylXrpwsWrRI1q1bp/x7dgWhIY0YMUIqV64sc+fOlSlTpkjZsmXl119/VZ5/9OiRAJCVK1dqbV/aEsugoCApW7asjBo1Sr7//nvx9PSUTp06aSRUXbp0kY4dOyqPIyMjxcbGRsaOHSsrVqyQ9957TywsLArdoqikPH78WOrWrSstWrSQJUuWSO/evcXGxkZjb7h169YJACUgHjhwQExMTGTQoEHy66+/ymeffSZmZmby3XffGeo0NAQFBcm6devk3XffFXt7e+V34tndCYyNjZVtsETyVsLb2NjI7Nmz5csvvxRzc3ON519nZ8+eFUtLS+nXr5/8+OOP0rhxY2nYsKHGCuqRI0dKrVq1itTGULKzs6Vt27bi5eUlixcvluHDh0u5cuU0Lh4cPXpUAMjJkydFJG+fWQAyadIkjTj67JtiQ1qyZImYmprK9OnT5X//+59UrFhRxo0bp1HHxcWlQFm+0pZY3rx5U+zt7SUgIECWLFkibdq00di+SkRk6tSpUrFiReVxQkKCuLq6Su/evWXFihXy8ccfi6mpqWzatMkAZ1BQQECA1KhRQxYtWiTjxo0TExMT+euvv5TnL126pLEN1q1bt8TW1lY6duwoy5cvl9mzZ4uDg0OBnWEM5dSpU7Ju3Tr57LPPBICsWLFC1q1bp/GGrXHjxhrbYG3YsEFMTEzk008/lfnz50u1atU0nteHwW/i69WrF44ePYqMjAycOXMGX3zxBdasWaNRx9fXF61atSpSG0Oxt7fHP//8g3bt2uHEiRNo0qQJIiMjNRZFeHp6IjAwULmHskqVKujZsyeOHz+ufOyxbds2nDlzxlCnoeGnn37CvHnzcPHiRTx69Aj79u3DkCFDlOfLlSuHwMBA1KhRQ2t7R0dHBAYGonz58iU15EJ16NAB4eHhMDU1RUREBD766CPs2rVL45tq2rZti/bt2yuP33jjDURERMDCwgJHjx5F7dq1ceXKlVLxUZuFhQVCQ0MREBCAEydOoGbNmoiMjETt2rWVOi4uLggMDFS++ahjx464cOECXFxccOTIEYgI/v77b3z22WeGOg0NISEh2LZtG4yMjNCxY0fld+LJkydKncDAQI1z/Prrr7Fy5UrcuHEDsbGx+PPPPzFx4kRDDL/E1a9fX4kzYWFhCAwMREhIiMZ9UU2bNkWXLl2K1MZQjI2NsXfvXowYMQInT56ElZUVIiIi4O3trdRxcHBAYGCgcg+lmZkZAgMDcevWLY04evToUUOdhoZRo0bhr7/+Qnx8PC5fvowff/wRCxcu1KjTrVs3vPnmm1rbm5qaIjAwsNR8K1uNGjUQGRmJevXqITQ0FG+99RbCw8M1Fgs2aNAA77zzjvLY3t4eERERaNSoEYKDg2FmZobw8HD07t3bAGdQ0ObNmzFp0iRl14/jx49r/M5YW1sjMDAQVapUAQC4urriypUr6NSpE0JDQ3H37l38+OOPpSYfOXfuHLZt24aoqCgEBgZi37592LZtG+7evavUeeutt5Tb1wDg3XffRXBwMJ48eYJz585h1qxZ+OOPP4rUr5FIKdkzhYiIiIheaQa/YklERERErwcmlkRERESkCiaWRERERKQKJpZEREREpAomlkRERESkCiaWRERERKQKJpZEREREpAomlkRERESkCiaWRERERKQKJpZEREREpAomlkRERESkCiaWRERERKSK/we1H24MfJY1FgAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 800x340 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Cell 4 [LECTURE / ADD-BACK 3] -- the kernel, computed (N = 8, dense SVD)\n",
    "N = 8\n",
    "mesh, W, a, L, bcs = cavity(N)\n",
    "A_mat = assemble(a, bcs=bcs, mat_type=\"aij\")\n",
    "ai, aj, av = A_mat.petscmat.getValuesCSR()\n",
    "M = sp.csr_matrix((av, aj, ai)).toarray()\n",
    "sv = sla.svd(M, compute_uv=False)\n",
    "tol = 1e-10 * sv.max()\n",
    "print(\"matrix size:\", M.shape[0], \" zero singular values:\",\n",
    "      int((sv < tol).sum()), \" smallest four:\", sv[-4:])\n",
    "\n",
    "# Pressure parts of the null vectors, evaluated at the cell centers.\n",
    "U_, S_, Vt = sla.svd(M)\n",
    "null = Vt[sv < tol]\n",
    "pts = centers(N)\n",
    "grids = []\n",
    "for vec in null:\n",
    "    wn = Function(W)\n",
    "    with wn.dat.vec_wo as vv:\n",
    "        vv.array[:] = vec\n",
    "    grids.append(eval_at(wn.subfunctions[1], pts).reshape(N, N))\n",
    "\n",
    "# An SVD returns an arbitrary orthonormal basis of the null space, which\n",
    "# generically mixes the constant and the checkerboard; the basis is therefore\n",
    "# rotated onto these two directions for display, and the residual of the span\n",
    "# identity is reported.\n",
    "ones = np.ones((N, N)) / N\n",
    "cb = cb_grid(N) / N\n",
    "G = np.array([g.ravel() for g in grids])\n",
    "c1 = G @ ones.ravel()\n",
    "c2 = G @ cb.ravel()\n",
    "resid = np.linalg.norm(G - np.outer(c1, ones.ravel())\n",
    "                         - np.outer(c2, cb.ravel()))\n",
    "print(\"null space vs. span{constant, checkerboard}: residual =\", resid)\n",
    "\n",
    "fig, ax = plt.subplots(1, 2, figsize=(8, 3.4))\n",
    "for a_, g, t in ((ax[0], ones, \"the constant\"),\n",
    "                 (ax[1], cb, \"the checkerboard\")):\n",
    "    a_.imshow(g, origin=\"lower\", cmap=\"cividis\",\n",
    "              vmin=-np.abs(g).max(), vmax=np.abs(g).max(),\n",
    "              extent=[0, 1, 0, 1], interpolation=\"nearest\")\n",
    "    a_.set_title(\"basis vector: \" + t, fontsize=10)\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cell-05",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Cell 5 [LECTURE] -- the regularized solve: inspect the two fields separately\n",
    "N = 32\n",
    "mesh, W, a, L, bcs = cavity(N)\n",
    "u_, p_ = TrialFunctions(W)\n",
    "v_, q_ = TestFunctions(W)\n",
    "a_eps = a - Constant(1e-8) * p_ * q_ * dx\n",
    "w = Function(W)\n",
    "solve(a_eps == L, w, bcs=bcs, solver_parameters={\n",
    "    \"mat_type\": \"aij\", \"ksp_type\": \"preonly\",\n",
    "    \"pc_type\": \"lu\", \"pc_factor_mat_solver_type\": \"mumps\"})\n",
    "u_h, p_h = w.subfunctions\n",
    "\n",
    "pgrid = eval_at(p_h, centers(N)).reshape(N, N)\n",
    "amp = float((pgrid * cb_grid(N)).mean())\n",
    "print(\"checkerboard amplitude of p_h:\", amp,\n",
    "      \"   max |p_h| =\", float(np.abs(pgrid).max()))\n",
    "\n",
    "xs = np.linspace(0, 1, N + 1)\n",
    "upts = [(float(a), float(b)) for b in xs for a in xs]\n",
    "uv = eval_at(u_h, upts)\n",
    "UX = uv[:, 0].reshape(N + 1, N + 1)\n",
    "UY = uv[:, 1].reshape(N + 1, N + 1)\n",
    "\n",
    "fig, ax = plt.subplots(1, 2, figsize=(9, 3.7))\n",
    "vm = np.abs(pgrid).max()\n",
    "ax[0].imshow(pgrid, origin=\"lower\", cmap=\"cividis\", vmin=-vm, vmax=vm,\n",
    "             extent=[0, 1, 0, 1], interpolation=\"nearest\")\n",
    "ax[0].set_title(\"pressure $p_h$\", fontsize=10)\n",
    "X, Y = np.meshgrid(xs, xs)\n",
    "ax[1].streamplot(X, Y, UX, UY, color=\"tab:blue\", density=1.1, linewidth=0.9)\n",
    "ax[1].set_aspect(\"equal\")\n",
    "ax[1].set_title(\"velocity streamlines\", fontsize=10)\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cell-06",
   "metadata": {},
   "source": [
    "## Summary\n",
    "\n",
    "The dimension count ($2(N-1)^2$ interior velocity unknowns against $N^2$\n",
    "pressures) is satisfied for all $N \\ge 4$, and the system is nevertheless\n",
    "singular: the null space of $B^T$ is spanned by the constant (physical, since\n",
    "$Q = L^2_0$) and by the checkerboard (spurious). The velocity field of the\n",
    "regularized solve appears correct while the pressure is dominated by the\n",
    "spurious mode; in mixed methods, each field must be inspected separately.\n",
    "\n",
    "The kernel proof is Exercise 4(a); the singular-value study as a function of\n",
    "$N$ is Exercise 4(b); the triangular case is examined in Exercise 4(c) and\n",
    "classified in Lecture 2.\n"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Firedrake",
   "language": "python",
   "name": "firedrake"
  },
  "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.12.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}