{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "fb29e0c0",
   "metadata": {},
   "source": [
    "# L2 — measuring the Brezzi constants (`L2_constants.ipynb`)\n",
    "\n",
    "**Setting.** Companion notebook to Lecture 2. Demo 1 (cells 4–5) measures\n",
    "the constants $\\beta_h$, $\\alpha_h$ of the one-dimensional mixed-Laplacian\n",
    "suite $P_1$–$P_1$ / $P_2$–$P_0$ / $P_1$–$P_0$ of Lecture 1. Demo 2\n",
    "(cells 7–8) examines the triangular $P_1$–$P_1$ Stokes pair: exact\n",
    "spurious pressure modes and the essential inf-sup constant.\n",
    "\n",
    "**Conventions.** In one dimension the $V$-norm is\n",
    "$\\|\\tau\\|_{\\mathrm{div}}^2 = \\|\\tau\\|_0^2 + \\|\\tau'\\|_0^2$ and the\n",
    "$Q$-norm is $\\|\\cdot\\|_0$; for Stokes the velocity norm is\n",
    "$\\|\\nabla\\cdot\\|_0$ and the pressure norm $\\|\\cdot\\|_0$.\n",
    "\n",
    "**Verification.** The notebook is self-contained. Every number printed\n",
    "here is additionally reproduced, in pure NumPy with exact quadrature and\n",
    "independently of Firedrake, by the scripts `verify_beta_1d.py` and\n",
    "`verify_p1p1_tri.py` from the course bundle.\n",
    "\n",
    "**Cell tags.** `[SETUP]`, `[LECTURE]` (shown during the lecture),\n",
    "`[ADD-BACK 1]` (optional lecture material), `[SELF-STUDY S1/S2]`\n",
    "(post-lecture; Problem 4 of exercise sheet 2 builds on these cells)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "043ad630",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Cell 1 [SETUP] -- Firedrake availability (FEM-on-Colab install if needed).\n",
    "try:\n",
    "    import firedrake\n",
    "    _firedrake_source = \"local installation found\"\n",
    "except ImportError:\n",
    "    import os, sys\n",
    "    if \"google.colab\" in sys.modules or \"COLAB_RELEASE_TAG\" in os.environ:\n",
    "        get_ipython().system('wget \"https://fem-on-colab.github.io/releases/firedrake-install-release-real.sh\" -O /tmp/firedrake-install.sh')\n",
    "        get_ipython().system('bash /tmp/firedrake-install.sh')\n",
    "        import firedrake\n",
    "        _firedrake_source = \"installed via FEM-on-Colab\"\n",
    "    else:\n",
    "        raise RuntimeError(\"Firedrake is required to run this notebook.\")\n",
    "def _firedrake_version():\n",
    "    v = getattr(firedrake, \"__version__\", None)\n",
    "    if v:\n",
    "        return v\n",
    "    try:\n",
    "        from importlib.metadata import version\n",
    "        return version(\"firedrake\")\n",
    "    except Exception:\n",
    "        return \"(version not reported)\"\n",
    "\n",
    "print(\"Firedrake\", _firedrake_version(), \"--\", _firedrake_source)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "bf9e64cd",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Cell 2 [LECTURE] -- 1D assembly and the two constants (dense linear algebra).\n",
    "import warnings\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from firedrake import *\n",
    "\n",
    "# petsc4py/Firedrake internals assign to ndarray.shape during matrix\n",
    "# extraction; NumPy >= 2.5 reports this as a DeprecationWarning. The warning\n",
    "# is harmless for these computations and is silenced by message.\n",
    "warnings.filterwarnings(\n",
    "    \"ignore\", message=\"Setting the shape on a NumPy array has been deprecated\")\n",
    "\n",
    "def dense(form):\n",
    "    \"Assemble a bilinear form and return it as a dense NumPy array.\"\n",
    "    M = assemble(form, mat_type=\"aij\").M.handle\n",
    "    m, n = M.getSize()\n",
    "    return np.array(M.getValues(range(m), range(n)))\n",
    "\n",
    "SPACES = {\"P1\": (\"CG\", 1), \"P2\": (\"CG\", 2), \"P0\": (\"DG\", 0)}\n",
    "\n",
    "def matrices_1d(n, vsp, qsp):\n",
    "    mesh = UnitIntervalMesh(n)\n",
    "    V = FunctionSpace(mesh, *SPACES[vsp])\n",
    "    Q = FunctionSpace(mesh, *SPACES[qsp])\n",
    "    s, t = TrialFunction(V), TestFunction(V)\n",
    "    Mv = dense(inner(s, t) * dx)\n",
    "    Kv = dense(inner(s.dx(0), t.dx(0)) * dx)\n",
    "    u, q = TrialFunction(V), TestFunction(Q)\n",
    "    B = dense(inner(u.dx(0), q) * dx)      # rows: Q dofs, columns: V dofs\n",
    "    r, w = TrialFunction(Q), TestFunction(Q)\n",
    "    Mq = dense(inner(r, w) * dx)\n",
    "    return mesh, Mv, Kv, B, Mq\n",
    "\n",
    "def gen_eigs(S, M):\n",
    "    \"Eigenvalues of the symmetric pencil S x = lambda M x (M SPD).\"\n",
    "    L = np.linalg.cholesky(M); Li = np.linalg.inv(L)\n",
    "    Sh = Li @ S @ Li.T\n",
    "    return np.linalg.eigvalsh(0.5 * (Sh + Sh.T))\n",
    "\n",
    "def constants_1d(n, vsp, qsp):\n",
    "    _, Mv, Kv, B, Mq = matrices_1d(n, vsp, qsp)\n",
    "    Gv = Mv + Kv                            # Gram matrix of ||.||_div\n",
    "    lam = gen_eigs(B @ np.linalg.solve(Gv, B.T), Mq)\n",
    "    beta = np.sqrt(max(lam[0], 0.0))\n",
    "    _, s, Vt = np.linalg.svd(B)             # ker B = coefficient vectors of K_h\n",
    "    rk = int(np.sum(s > 1e-10 * s[0]))\n",
    "    Z = Vt[rk:].T\n",
    "    mu = gen_eigs(Z.T @ Mv @ Z, Z.T @ Gv @ Z)\n",
    "    return beta, mu[0], Z.shape[1]\n",
    "\n",
    "b, a, dk = constants_1d(4, \"P1\", \"P0\")\n",
    "print(\"assembly helpers ready: matrices_1d, gen_eigs, constants_1d\")\n",
    "print(f\"smoke test, P1-P0 at n=4: beta_h = {b:.6f}, alpha_h = {a:.6f}, \"\n",
    "      f\"dim K_h = {dk}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "244a7520",
   "metadata": {},
   "source": [
    "## Demo 1 — names for the measured quantities (slides 23–24)\n",
    "\n",
    "The weighted eigenvalue problems of Lemma 2 (matrix form, slide 9) applied\n",
    "to the Lecture 1 suite. Expected exact entries: $\\beta_h = 0$ and\n",
    "$\\alpha_h = 1$ for $P_1$–$P_1$; $\\alpha_h = 1$ for $P_1$–$P_0$;\n",
    "$\\alpha_h/h^2 \\uparrow 1/60$ for $P_2$–$P_0$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fa5032a0",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Cell 4 [LECTURE] -- Demo 1a: the constants table, with assertions.\n",
    "PAIRS = [(\"P1\", \"P1\"), (\"P2\", \"P0\"), (\"P1\", \"P0\")]\n",
    "NS = [8, 16, 32, 64]\n",
    "results = {}\n",
    "print(f\"{'pair':7s}{'n':>4s}  {'beta_h':>13s}  {'alpha_h':>13s}\"\n",
    "      f\"  {'alpha_h/h^2':>12s}  {'dim K_h':>8s}\")\n",
    "for vsp, qsp in PAIRS:\n",
    "    for n in NS:\n",
    "        b, a, dk = constants_1d(n, vsp, qsp)\n",
    "        results[(vsp, qsp, n)] = (b, a, dk)\n",
    "        ratio = f\"{a*n*n:12.6f}\" if (vsp, qsp) == (\"P2\", \"P0\") else f\"{'--':>12s}\"\n",
    "        print(f\"{vsp}-{qsp}  {n:4d}  {b:13.6e}  {a:13.6e}  {ratio}  {dk:8d}\")\n",
    "    print()\n",
    "for n in NS:\n",
    "    b, a, dk = results[(\"P1\", \"P1\", n)]\n",
    "    assert b < 1e-7 and abs(a - 1) < 1e-10 and dk == 1\n",
    "    b, a, dk = results[(\"P1\", \"P0\", n)]\n",
    "    assert abs(a - 1) < 1e-10 and dk == 1 and b > (2/3)**0.5 - 1e-12\n",
    "    b2, a2, dk2 = results[(\"P2\", \"P0\", n)]\n",
    "    assert dk2 == n + 1 and b2 >= b - 1e-12         # monotonicity in V_h\n",
    "ratios = [results[(\"P2\", \"P0\", n)][1] * n * n for n in NS]\n",
    "assert all(abs(r - 1/60) < 1e-4 for r in ratios)\n",
    "assert all(ratios[i] <= ratios[i+1] + 1e-14 for i in range(len(NS)-1))\n",
    "print(\"assertions passed (values match verify_beta_1d.py)\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "85abd8a5",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Cell 5 [LECTURE] -- Demo 1b: constants under refinement; blow-up factors.\n",
    "beta_exact = np.pi / np.sqrt(1 + np.pi**2)\n",
    "fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(9, 3.2))\n",
    "for (vsp, qsp), style in zip(PAIRS, [\"x-\", \"s--\", \"o-\"]):\n",
    "    bs = [results[(vsp, qsp, n)][0] for n in NS]\n",
    "    As = [results[(vsp, qsp, n)][1] for n in NS]\n",
    "    ax1.semilogx(NS, bs, style, label=f\"{vsp}-{qsp}\")\n",
    "    ax2.loglog(NS, As, style)\n",
    "ax1.axhline(beta_exact, ls=\":\", c=\"gray\")\n",
    "ax1.set(xlabel=\"n\", ylabel=\"beta_h\", ylim=(-0.05, 1.05)); ax1.legend()\n",
    "ax2.loglog(NS, [1/(60*n*n) for n in NS], \":\", c=\"gray\")\n",
    "ax2.set(xlabel=\"n\", ylabel=\"alpha_h\")\n",
    "plt.tight_layout(); plt.show()\n",
    "print(f\"exact continuous constant  pi/sqrt(1+pi^2) = {beta_exact:.6f}\")\n",
    "print(\"blow-up factor (1+1/alpha_h)(1+1/beta_h):\")\n",
    "for vsp, qsp in [(\"P2\", \"P0\"), (\"P1\", \"P0\")]:\n",
    "    for n in [8, 64]:\n",
    "        b, a, _ = results[(vsp, qsp, n)]\n",
    "        print(f\"  {vsp}-{qsp}, n={n:3d}: {(1+1/a)*(1+1/b): .3e}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ebb3c754",
   "metadata": {},
   "source": [
    "## Demo 2 — the triangular $P_1$–$P_1$ pair (slides 26–27)\n",
    "\n",
    "Stokes setting, velocity seminorm. The pair possesses **exact** spurious\n",
    "pressure modes: eight on the structured mesh (either diagonal, every\n",
    "$N \\ge 8$) and five on randomly perturbed meshes of the same connectivity\n",
    "(`verify_p1p1_tri.py`). After removal of the exact modes, the essential\n",
    "constant $\\widetilde\\beta_h$ decays approximately like $O(h)$: the pair\n",
    "fails through both failure mechanisms identified in Lecture 1."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0a4a7bb9",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Cell 7 [LECTURE] -- Demo 2a: counting the kernel; one mode.\n",
    "def stokes_p1p1(N, diagonal=\"left\"):\n",
    "    mesh = UnitSquareMesh(N, N, diagonal=diagonal)\n",
    "    V = VectorFunctionSpace(mesh, \"CG\", 1)\n",
    "    Q = FunctionSpace(mesh, \"CG\", 1)\n",
    "    u, v = TrialFunction(V), TestFunction(V)\n",
    "    K = dense(inner(grad(u), grad(v)) * dx)\n",
    "    p, q = TrialFunction(V), TestFunction(Q)\n",
    "    B = dense(-div(p) * q * dx)             # rows: Q dofs, columns: V dofs\n",
    "    r, w = TrialFunction(Q), TestFunction(Q)\n",
    "    Mq = dense(inner(r, w) * dx)\n",
    "    # interior velocity dofs, located by node coordinates (build-stable;\n",
    "    # same technique as L1A cell 16)\n",
    "    X = mesh.coordinates.dat.data_ro\n",
    "    interior_nodes = np.where((X[:, 0] > 1e-12) & (X[:, 0] < 1 - 1e-12)\n",
    "                              & (X[:, 1] > 1e-12) & (X[:, 1] < 1 - 1e-12))[0]\n",
    "    idx = np.sort(np.concatenate([2*interior_nodes, 2*interior_nodes + 1]))\n",
    "    Ki, Bi = K[np.ix_(idx, idx)], B[:, idx]\n",
    "    assert idx.size == 2*(N-1)**2\n",
    "    # the constant pressure must lie in ker B^T exactly (all dofs interior):\n",
    "    assert np.linalg.norm(np.ones(Mq.shape[0]) @ Bi) < 1e-9\n",
    "    return mesh, Ki, Bi, Mq, X\n",
    "\n",
    "def tri_spectrum(N, diagonal=\"left\"):\n",
    "    mesh, Ki, Bi, Mq, X = stokes_p1p1(N, diagonal)\n",
    "    lam = gen_eigs(Bi @ np.linalg.solve(Ki, Bi.T), Mq)\n",
    "    nz = int(np.sum(np.abs(lam) < 1e-10 * lam[-1]))\n",
    "    return mesh, Mq, X, lam, nz\n",
    "\n",
    "mesh, Mq, X, lam, nz = tri_spectrum(16)\n",
    "print(f\"N=16: numerically zero eigenvalues: {nz}   (expected: 8)\")\n",
    "assert nz == 8\n",
    "# one nonconstant mode, for the figure of slide 25\n",
    "_, _, Bi16, Mq16, _ = stokes_p1p1(16)\n",
    "_, s, Vt = np.linalg.svd(Bi16.T, full_matrices=True)\n",
    "Z = Vt[int(np.sum(s > 1e-10 * s[0])):].T\n",
    "one = np.ones(Mq16.shape[0])\n",
    "Zc = Z - np.outer(one, (one @ Mq16 @ Z) / (one @ Mq16 @ one))\n",
    "Qm, R = np.linalg.qr(Zc)\n",
    "mode = Qm[:, int(np.argmax(np.abs(np.diag(R)) > 1e-8))]\n",
    "import matplotlib.tri as mtri\n",
    "cells16 = mesh.coordinates.cell_node_map().values\n",
    "T = mtri.Triangulation(X[:, 0], X[:, 1], cells16)\n",
    "plt.figure(figsize=(4, 3.6))\n",
    "plt.tricontourf(T, mode / np.abs(mode).max(), levels=21, cmap=\"RdBu_r\")\n",
    "plt.gca().set_aspect(\"equal\"); plt.xticks([]); plt.yticks([]); plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "144bb1aa",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Cell 8 [LECTURE] -- Demo 2b: the essential constant under refinement.\n",
    "print(\" N     zero modes   beta_tilde\")\n",
    "bt = {}\n",
    "for N in [8, 12, 16, 24, 32, 40]:\n",
    "    _, _, _, lam, nz = tri_spectrum(N)\n",
    "    assert nz == 8, (N, nz)\n",
    "    bt[N] = np.sqrt(lam[nz])\n",
    "    print(f\"{N:3d}       {nz}        {bt[N]: .6e}\")\n",
    "Ns = sorted(bt)\n",
    "rate = np.polyfit(np.log([1/N for N in Ns[-3:]]),\n",
    "                  np.log([bt[N] for N in Ns[-3:]]), 1)[0]\n",
    "print(f\"fitted decay rate, N=16-32: {rate:.3f}; \"\n",
    "      f\"verify_p1p1_tri.py obtains 0.97 on N=24-40\")\n",
    "print(\"mode-count fragility (verify_p1p1_tri.py): 8 structured, \"\n",
    "      \"5 randomly perturbed, 7 alternating pattern (perturbed)\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3fa29058",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Cell 9 [ADD-BACK 1] -- quadrilateral Q1-P0 beside triangular P1-P1.\n",
    "# Both pairs are failure species (i); the kernel dimensions differ: 2 vs 8.\n",
    "def quad_q1p0(N):\n",
    "    mesh = UnitSquareMesh(N, N, quadrilateral=True)\n",
    "    V = VectorFunctionSpace(mesh, \"CG\", 1)\n",
    "    Q = FunctionSpace(mesh, \"DG\", 0)\n",
    "    u, v = TrialFunction(V), TestFunction(V)\n",
    "    K = dense(inner(grad(u), grad(v)) * dx)\n",
    "    p, q = TrialFunction(V), TestFunction(Q)\n",
    "    B = dense(-div(p) * q * dx)\n",
    "    r, w = TrialFunction(Q), TestFunction(Q)\n",
    "    Mq = dense(inner(r, w) * dx)\n",
    "    X = mesh.coordinates.dat.data_ro\n",
    "    nodes = np.where((X[:, 0] > 1e-12) & (X[:, 0] < 1-1e-12)\n",
    "                     & (X[:, 1] > 1e-12) & (X[:, 1] < 1-1e-12))[0]\n",
    "    idx = np.sort(np.concatenate([2*nodes, 2*nodes + 1]))\n",
    "    lam = gen_eigs(B[:, idx] @ np.linalg.solve(K[np.ix_(idx, idx)],\n",
    "                                               B[:, idx].T), Mq)\n",
    "    return lam\n",
    "\n",
    "lam_q = quad_q1p0(8)\n",
    "_, _, _, lam_t, _ = tri_spectrum(8)\n",
    "print(\"smallest four weighted eigenvalues, N = 8\")\n",
    "print(\"  Q1-P0 (quadrilaterals):\",\n",
    "      np.array2string(lam_q[:4], precision=3, max_line_width=120))\n",
    "print(\"  P1-P1 (triangles)     :\",\n",
    "      np.array2string(lam_t[:9], precision=3, max_line_width=120))\n",
    "nzq = int(np.sum(np.abs(lam_q) < 1e-10 * lam_q[-1]))\n",
    "print(f\"exact kernel dimensions: {nzq} (constant + checkerboard, cf. L1) \"\n",
    "      f\"vs 8\")\n",
    "assert nzq == 2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5569b9f2",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Cell 10 [SELF-STUDY S1 / Problem 4(c*)] -- worst-case data for P2-P0.\n",
    "# Exact discrete solution (w_h, 0) for F = a(., w_h), G = 0, w_h in K_h.\n",
    "# Expected compensated ratios: C_h*h^{3/2} -> 3.33 (single bubble),\n",
    "# C_h*h^2 -> 5.847 (alternating bubbles); cf. exercise sheet,\n",
    "# Problem 4(c*).\n",
    "print(\" n    single: C*h^1.5    alternating: C*h^2\")\n",
    "for n in [8, 16, 32, 64, 128]:\n",
    "    mesh, Mv, Kv, B, Mq = matrices_1d(n, \"P2\", \"P0\")\n",
    "    Gv = Mv + Kv\n",
    "    Vsp = FunctionSpace(mesh, \"CG\", 2)\n",
    "    xs = Function(Vsp).interpolate(SpatialCoordinate(mesh)[0]).dat.data_ro\n",
    "    # midpoint (bubble) dofs of CG2, located by coordinates\n",
    "    mids = np.where(np.abs((xs * n) % 1 - 0.5) < 1e-8)[0]\n",
    "    assert mids.size == n\n",
    "    order = np.argsort(xs[mids]); mids = mids[order]\n",
    "    out = []\n",
    "    for kind in [\"single\", \"alt\"]:\n",
    "        w = np.zeros(Vsp.dim())\n",
    "        if kind == \"single\":\n",
    "            w[mids[n // 2]] = 1.0\n",
    "        else:\n",
    "            w[mids] = (-1.0) ** np.arange(n)\n",
    "        assert np.linalg.norm(B @ w) < 1e-12          # w in K_h\n",
    "        f = Mv @ w                                    # F = a(., w)\n",
    "        dual = np.sqrt(f @ np.linalg.solve(Gv, f))\n",
    "        C = np.sqrt(w @ Gv @ w) / dual\n",
    "        out.append(C)\n",
    "    print(f\"{n:4d}   {out[0]/n**1.5: .5f}            {out[1]/n**2: .5f}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0dc7f7c0",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Cell 11 [SELF-STUDY S2 / Problem 4(d**)] -- the exact modes under mesh\n",
    "# perturbation. Self-contained; the script verify_p1p1_tri.py (course\n",
    "# bundle) reproduces the same counts in pure NumPy and additionally covers\n",
    "# the alternating-diagonal pattern.\n",
    "def p1p1_kernel_count(N, perturb=0.0, seed=0):\n",
    "    mesh = UnitSquareMesh(N, N)\n",
    "    X = mesh.coordinates.dat.data\n",
    "    interior = np.where((X[:, 0] > 1e-12) & (X[:, 0] < 1 - 1e-12)\n",
    "                        & (X[:, 1] > 1e-12) & (X[:, 1] < 1 - 1e-12))[0]\n",
    "    if perturb > 0:\n",
    "        rng = np.random.default_rng(seed)\n",
    "        X[interior] += (perturb / N) * rng.uniform(-1, 1, (interior.size, 2))\n",
    "    V = VectorFunctionSpace(mesh, \"CG\", 1)\n",
    "    Q = FunctionSpace(mesh, \"CG\", 1)\n",
    "    u, v = TrialFunction(V), TestFunction(V)\n",
    "    K = dense(inner(grad(u), grad(v)) * dx)\n",
    "    p, q = TrialFunction(V), TestFunction(Q)\n",
    "    B = dense(-div(p) * q * dx)\n",
    "    r, w = TrialFunction(Q), TestFunction(Q)\n",
    "    Mq = dense(inner(r, w) * dx)\n",
    "    idx = np.sort(np.concatenate([2 * interior, 2 * interior + 1]))\n",
    "    lam = gen_eigs(B[:, idx] @ np.linalg.solve(K[np.ix_(idx, idx)],\n",
    "                                               B[:, idx].T), Mq)\n",
    "    return int(np.sum(np.abs(lam) < 1e-10 * lam[-1]))\n",
    "\n",
    "for label, kw, expected in [\n",
    "        (\"structured\", dict(), 8),\n",
    "        (\"perturbed (0.2h, seed 1)\", dict(perturb=0.2, seed=1), 5),\n",
    "        (\"perturbed (0.2h, seed 2)\", dict(perturb=0.2, seed=2), 5)]:\n",
    "    cnt = p1p1_kernel_count(8, **kw)\n",
    "    print(f\"N=8, {label:26s}: dim ker B^T = {cnt}   (expected: {expected})\")\n",
    "    assert cnt == expected\n",
    "print(\"three of the structured modes are symmetry artifacts and vanish\")\n",
    "print(\"under a generic perturbation; four modes and the constant are\")\n",
    "print(\"topological. Alternating-diagonal pattern, perturbed: 7 modes\")\n",
    "print(\"(verify_p1p1_tri.py).\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6b654ada",
   "metadata": {},
   "source": [
    "## Summary and pointers\n",
    "\n",
    "Measured in this notebook: the two Brezzi constants of the one-dimensional\n",
    "suite (cells 4–5), with the exact entries $\\beta_h(P_1\\text{–}P_1) = 0$\n",
    "and $\\alpha_h(P_1\\text{–}P_1) = \\alpha_h(P_1\\text{–}P_0) = 1$, the sharp\n",
    "ratio $\\alpha_h/h^2 \\to 1/60$, and the limit\n",
    "$\\beta_h \\to \\pi/\\sqrt{1+\\pi^2}$; and the two-layer failure of the\n",
    "triangular $P_1$–$P_1$ pair (cells 7–8), with eight exact spurious modes\n",
    "and an essential constant decaying like $O(h)$.\n",
    "\n",
    "Exercise sheet 2 builds on these measurements. Problem 1 proves the matrix\n",
    "characterizations used in cell 2 and the monotonicity observed in cell 4;\n",
    "Problem 2 proves the uniform bound behind the flat curves of cell 5;\n",
    "Problem 4 extends cell 10 (worst-case data) and cell 11 (the exact modes;\n",
    "the classification question is open).\n",
    "\n",
    "The pure-NumPy scripts `verify_beta_1d.py` and `verify_p1p1_tri.py`\n",
    "reproduce every printed number independently of Firedrake; they are the\n",
    "reference implementations for this notebook."
   ]
  }
 ],
 "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
}