diff --git a/.gitignore b/.gitignore
index 6f7521e..30f47bb 100644
--- a/.gitignore
+++ b/.gitignore
@@ -1,5 +1,6 @@
-## Checkpoint files
+## Output files
 *.h5
+*.json
 
 
 ## Meshes
diff --git a/Makefile b/Makefile
index 25efb39..fc8d189 100644
--- a/Makefile
+++ b/Makefile
@@ -1,6 +1,7 @@
 
 tensor-product-stokes.pdf: tensor-product-stokes.tex tensor-product-stokes.bib
 	cd demo/flow-around-cylinder && make cylinder_flow.pdf && cd ../../
+	cd demo/lbb-stability && make results-2d.pdf && cd ../../
 	pdflatex tensor-product-stokes.tex
 	bibtex tensor-product-stokes
 	pdflatex tensor-product-stokes.tex
diff --git a/demo/lbb-stability/Makefile b/demo/lbb-stability/Makefile
index 63074cf..e764fac 100644
--- a/demo/lbb-stability/Makefile
+++ b/demo/lbb-stability/Makefile
@@ -1,10 +1,8 @@
-all: taylor-hood-2d.png taylor-hood-3d.png taylor-hood-3d-2layer.png
 
-taylor-hood-2d.png: lbb-stability.py
-	python3 lbb-stability.py --dimension 2 --output $@
+results-2d.pdf: results-2d.json make_plots.py
+	python make_plots.py --input $< --output $@
 
-taylor-hood-3d.png: lbb-stability.py
-	python3 lbb-stability.py --dimension 3 --output $@
-
-taylor-hood-3d-2layer.png: lbb-stability.py
-	python3 lbb-stability.py --dimension 3 --num-layers 2 --output $@
+results-2d.json: lbb-stability.py
+	python lbb-stability.py --dimension 2 --element mini --output $@
+	python lbb-stability.py --dimension 2 --element taylor-hood --output $@
+	python lbb-stability.py --dimension 2 --element crouzeix-raviart --output $@
diff --git a/demo/lbb-stability/lbb-stability.py b/demo/lbb-stability/lbb-stability.py
index 80b1965..27ab8f3 100644
--- a/demo/lbb-stability/lbb-stability.py
+++ b/demo/lbb-stability/lbb-stability.py
@@ -1,11 +1,12 @@
+import json
 import argparse
+import pathlib
 import firedrake
 from firedrake import assemble, inner, sym, grad, div, dx
 import matplotlib.pyplot as plt
 import numpy as np
 from petsc4py import PETSc
 from slepc4py import SLEPc
-import tqdm
 
 
 def compute_inf_sup_constant(spaces, b, inner_products):
@@ -16,18 +17,18 @@ def compute_inf_sup_constant(spaces, b, inner_products):
     u, p = firedrake.TrialFunctions(Z)
     v, q = firedrake.TestFunctions(Z)
 
-    A = assemble(-(m(u, v) + b(v, p) + b(u, q)), mat_type='aij').M.handle
-    M = assemble(n(p, q), mat_type='aij').M.handle
+    A = assemble(-(m(u, v) + b(v, p) + b(u, q)), mat_type="aij").M.handle
+    M = assemble(n(p, q), mat_type="aij").M.handle
 
     opts = PETSc.Options()
-    opts.setValue('eps_gen_hermitian', None)
-    opts.setValue('eps_target_real', None)
-    opts.setValue('eps_smallest_magnitude', None)
-    opts.setValue('st_type', 'sinvert')
-    opts.setValue('st_ksp_type', 'preonly')
-    opts.setValue('st_pc_type', 'lu')
-    opts.setValue('st_pc_factor_mat_solver_type', 'mumps')
-    opts.setValue('eps_tol', 1e-8)
+    opts.setValue("eps_gen_hermitian", None)
+    opts.setValue("eps_target_real", None)
+    opts.setValue("eps_smallest_magnitude", None)
+    opts.setValue("st_type", "sinvert")
+    opts.setValue("st_ksp_type", "preonly")
+    opts.setValue("st_pc_type", "lu")
+    opts.setValue("st_pc_factor_mat_solver_type", "mumps")
+    opts.setValue("eps_tol", 1e-8)
 
     num_values = 1
     eigensolver = SLEPc.EPS().create(comm=firedrake.COMM_WORLD)
@@ -41,67 +42,87 @@ def compute_inf_sup_constant(spaces, b, inner_products):
     return λ
 
 
+def run(*, dimension, num_cells, element, num_layers=None, vdegree=None):
+    mesh = firedrake.UnitSquareMesh(num_cells, num_cells, diagonal="crossed")
+    h = mesh.cell_sizes.dat.data_ro[:].min()
+
+    cg_1 = firedrake.FiniteElement("CG", "triangle", 1)
+    cg_2 = firedrake.FiniteElement("CG", "triangle", 2)
+    b_3 = firedrake.FiniteElement("Bubble", "triangle", 3)
+    if element == "taylor-hood":
+        velocity_element = cg_2
+        pressure_element = cg_1
+    elif element == "mini":
+        velocity_element = cg_1 + b_3
+        pressure_element = cg_1
+    elif element == "crouzeix-raviart":
+        velocity_element = cg_2 + b_3
+        pressure_element = firedrake.FiniteElement("DG", "triangle", 1)
+
+    elif dimension == 3:
+        cg_z = firedrake.FiniteElement("CG", "interval", vdegree)
+        dg_z = firedrake.FiniteElement("DG", "interval", vdegree - 1)
+        velocity_element = firedrake.TensorProductElement(velocity_element, cg_z)
+        pressure_element = firedrake.TensorProductElement(pressure_element, dg_z)
+        mesh = firedrake.ExtrudedMesh(mesh, num_layers)
+
+    V = firedrake.VectorFunctionSpace(mesh, velocity_element)
+    Q = firedrake.FunctionSpace(mesh, pressure_element)
+
+    # This is always just going to be equal to 1 in our case, but in general
+    # you should compute a characteristic length scale for the domain when
+    # defining the H1 norm.
+    d = mesh.topological_dimension()
+    volume = assemble(firedrake.Constant(1) * dx(mesh))
+    l = firedrake.Constant(volume ** (1 / d))
+
+    b = lambda v, q: -q * div(v) * dx
+    m = lambda u, v: (inner(u, v) + l**2 * inner(grad(u), grad(v))) * dx
+    n = lambda p, q: p * q * dx
+    try:
+        λ = compute_inf_sup_constant((V, Q), b, (m, n)).real
+    except PETSc.Error:
+        λ = np.nan
+    return h, λ
+
+
 parser = argparse.ArgumentParser()
-parser.add_argument('--dimension', type=int, choices=[2, 3])
-parser.add_argument('--log-dx-min', type=int, default=2)
-parser.add_argument('--log-dx-max', type=int, default=5)
-parser.add_argument('--num-layers', type=int, default=1)
-parser.add_argument('--zdegree-min', type=int, default=1)
-parser.add_argument('--zdegree-max', type=int, default=6)
-parser.add_argument('--output')
+parser.add_argument("--dimension", type=int, choices=[2, 3])
+elements = ["taylor-hood", "crouzeix-raviart", "mini"]
+parser.add_argument("--element", choices=elements, default="taylor-hood")
+parser.add_argument("--num-cells-min", type=int, default=4)
+parser.add_argument("--num-cells-max", type=int, default=32)
+parser.add_argument("--num-cells-step", type=int, default=1)
+parser.add_argument("--num-layers", type=int, default=1)
+parser.add_argument("--vdegree", type=int, default=1)
+parser.add_argument("--output")
 args = parser.parse_args()
 
-
-if args.dimension == 2:
-    stability_constants = []
-    for N in tqdm.trange(2**args.log_dx_min, 2**args.log_dx_max + 1):
-        mesh = firedrake.UnitSquareMesh(N, N)
-        d = mesh.topological_dimension()
-        volume = assemble(firedrake.Constant(1) * dx(mesh))
-        l = firedrake.Constant(volume ** (1 / d))
-
-        V = firedrake.VectorFunctionSpace(mesh, 'CG', 2)
-        Q = firedrake.FunctionSpace(mesh, 'CG', 1)
-        try:
-            b = lambda v, q: -q * div(v) * dx
-            m = lambda u, v: (inner(u, v) + l**2 * inner(grad(u), grad(v))) * dx
-            n = lambda p, q: p * q * dx
-            λ = compute_inf_sup_constant((V, Q), b, (m, n))
-        except PETSc.Error:
-            λ = np.nan
-        stability_constants.append(λ)
-
-    fig, axes = plt.subplots()
-    axes.plot(abs(np.array(stability_constants)))
-    fig.savefig(args.output)
-
-elif args.dimension == 3:
-    kmin, kmax = args.zdegree_min, args.zdegree_max + 1
-    fig, axes = plt.subplots()
-    for k in tqdm.trange(kmin, kmax):
-        stability_constants = []
-        for N in tqdm.trange(2**args.log_dx_min, 2**args.log_dx_max + 1):
-            mesh2d = firedrake.UnitSquareMesh(N, N, diagonal='crossed')
-            mesh = firedrake.ExtrudedMesh(mesh2d, args.num_layers)
-            d = mesh.topological_dimension()
-            volume = assemble(firedrake.Constant(1) * dx(mesh))
-            l = firedrake.Constant(volume ** (1 / d))
-
-            V = firedrake.VectorFunctionSpace(mesh, 'CG', 2, vfamily='GLL', vdegree=k + 1)
-            Q = firedrake.FunctionSpace(mesh, 'CG', 1, vfamily='GL', vdegree=k)
-            try:
-                b = lambda v, q: -q * div(v) * dx
-                m = lambda u, v: (inner(u, v) + l**2 * inner(grad(u), grad(v))) * dx
-                n = lambda p, q: p * q * dx
-                λ = compute_inf_sup_constant((V, Q), b, (m, n))
-            except PETSc.Error:
-                λ = np.nan
-            stability_constants.append(λ)
-
-        axes.plot(abs(np.array(stability_constants)), label=str(k))
-
-    axes.legend()
-    fig.savefig(args.output)
-
-else:
-    raise ValueError("Dimension must be 2 or 3!")
+stability_constants = []
+nmin, nmax, nstep = args.num_cells_min, args.num_cells_max, args.num_cells_step
+for num_cells in range(nmin, nmax + 1, nstep):
+    h, λ = run(
+        dimension=args.dimension,
+        num_cells=num_cells,
+        element=args.element,
+        num_layers=args.num_layers,
+        vdegree=args.vdegree
+    )
+    stability_constants.append((h, λ))
+    print(".", flush=True, end="")
+
+data = []
+if pathlib.Path(args.output).is_file():
+    with open(args.output, "r") as output_file:
+        data = json.load(output_file)
+
+result = {
+    "dimension": args.dimension,
+    "element": args.element,
+    "results": stability_constants,
+}
+if args.dimension == 3:
+    result.update({"num_layers": args.num_layers, "vdegree": args.vdegree})
+data.append(result)
+with open(args.output, "w") as output_file:
+    json.dump(data, output_file)
diff --git a/demo/lbb-stability/make_plots.py b/demo/lbb-stability/make_plots.py
new file mode 100644
index 0000000..4087fc0
--- /dev/null
+++ b/demo/lbb-stability/make_plots.py
@@ -0,0 +1,25 @@
+import json
+import argparse
+import numpy as np
+import matplotlib.pyplot as plt
+
+parser = argparse.ArgumentParser()
+parser.add_argument("--input")
+parser.add_argument("--output")
+args = parser.parse_args()
+
+with open(args.input, "r") as input_file:
+    data = json.load(input_file)
+
+fig, ax = plt.subplots()
+ax.set_title("Empirical inf-sup constants")
+ax.set_xscale("log", base=2)
+ax.set_ylim((0.0, 1.0))
+ax.set_xlabel("Mesh spacing")
+ax.set_ylabel("inf-sup constant")
+for entry in data:
+    element = entry["element"]
+    hs, λs = zip(*entry["results"])
+    ax.scatter(hs, λs, label=element)
+ax.legend()
+fig.savefig(args.output, bbox_inches="tight")
diff --git a/tensor-product-stokes.tex b/tensor-product-stokes.tex
index 0c60a4a..1ffee16 100644
--- a/tensor-product-stokes.tex
+++ b/tensor-product-stokes.tex
@@ -272,9 +272,16 @@ \section{Demonstrations}
 
 \subsection{Empirical confirmation of inf-sup stability}
 
-Given a mesh of the spatial domain and a pair of discrete spaces $V^h$, $Q^h$, we can check empirically whether they satisfy the inf-sup condition by solving a generalized eigenvalue problem \citep{rognes2012automated}.
-This empirical check is not a substitute for a formal proof of inf-sup stability, but rather serves as a ``smoke test'' that our numerical implementation and our proofs are not obviously wrong.
+\begin{figure}
+    \begin{center}
+        \includegraphics[width=0.6\linewidth]{demo/lbb-stability/results-2d.pdf}
+    \end{center}
+    \caption{Empirical inf-sup constants for the three 2D element families on the unit square as the mesh is refined.
+    These elements are known to stable for the 2D Stokes problem.}
+    \label{fig:empirical-lbb-2d}
+\end{figure}
 
+Given a mesh of the spatial domain and a pair of discrete spaces $V^h$, $Q^h$, we can check empirically whether they satisfy the inf-sup condition by solving a generalized eigenvalue problem \citep{rognes2012automated}.
 Let $M$ and $N$ be the Riesz representers for the canonical mapping from $V \to V^*$ and $Q \to Q^*$ respectively.
 Consider the following eigenproblem: find velocity-pressure pairs $u$, $p$ and scalars $\lambda$ such that
 \begin{align}
@@ -283,6 +290,12 @@ \subsection{Empirical confirmation of inf-sup stability}
 \end{align}
 The eigenvalues of this problem are all real and positive, and the square root of the smallest eigenvalue is equal to the inf-sup constant for $V^h$, $Q^h$ \citep{malkus1981eigenproblems}.
 We can form this eigenproblem in Firedrake and then call out to the SLEPc package to solve it \citep{hernandez2005slepc}.
+This empirical check is not a substitute for a formal proof of inf-sup stability, but rather serves as a ``smoke test'' that our numerical implementation and our proofs are not obviously wrong.
+
+We performed the empirical inf-sup stability check using the unit square as our footprint domain, with a mesh spacing starting at 1/4 and going down to 1/32.
+Figure \ref{fig:empirical-lbb-2d} shows the results for 2D elements which we know to be inf-sup stable.
+\textcolor{red}{Do the 3D elements.}
+
 
 \subsection{Lid-driven cavity flow}