diff --git a/python/demo/demo_navier_stokes.py b/python/demo/demo_navier_stokes.py
new file mode 100644
index 00000000000..facaf7f816f
--- /dev/null
+++ b/python/demo/demo_navier_stokes.py
@@ -0,0 +1,431 @@
+# ---
+# jupyter:
+#   jupytext:
+#     text_representation:
+#       extension: .py
+#       format_name: light
+#       format_version: '1.5'
+#       jupytext_version: 1.13.6
+# ---
+
+# # Divergence conforming discontinuous Galerkin method for the Navier--Stokes equations
+#
+# This demo ({download}`demo_navier_stokes.py`) illustrates how to
+# implement a divergence conforming discontinuous Galerkin method for
+# the Navier-Stokes equations in FEniCSx. The method conserves mass
+# exactly and uses upwinding. The formulation is based on a combination
+# of "A fully divergence-free finite element method for
+# magnetohydrodynamic equations" by Hiptmair et al., "A Note on
+# Discontinuous Galerkin Divergence-free Solutions of the Navier-Stokes
+# Equations" by Cockburn et al, and "On the Divergence Constraint in
+# Mixed Finite Element Methods for Incompressible Flows" by John et al.
+#
+#
+# ## Governing equations
+#
+# We consider the incompressible Navier-Stokes equations in a domain
+# $\Omega \subset \mathbb{R}^d$, $d \in \{2, 3\}$, and time interval
+# $(0, \infty)$, given by
+#
+# $$
+# \partial_t u - \nu \Delta u + (u \cdot \nabla)u + \nabla p &= f \textnormal{ in } \Omega_t, \\
+# \nabla \cdot u &= 0 \textnormal{ in } \Omega_t,
+# $$
+#
+# where $u: \Omega_t \to \mathbb{R}^d$ is the velocity field,
+# $p: \Omega_t \to \mathbb{R}$ is the pressure field,
+# $f: \Omega_t \to \mathbb{R}^d$ is a prescribed force, $\nu \in \mathbb{R}^+$
+# is the kinematic viscosity, and $\Omega_t := \Omega \times (0, \infty)$.
+#
+# The problem is supplemented with the initial condition
+#
+# $$
+#     u(x, 0) = u_0(x) \textnormal{ in } \Omega
+# $$
+#
+# and boundary condition
+#
+# $$
+#     u = u_D \textnormal{ on } \partial \Omega \times (0, \infty),
+# $$
+#
+# where $u_0: \Omega \to \mathbb{R}^d$ is a prescribed initial velocity field
+# which satisfies the divergence free condition. The pressure field is only
+# determined up to a constant, so we seek the unique pressure field satisfying
+#
+# $$
+#     \int_\Omega p = 0.
+# $$
+#
+#
+#
+# ## Discrete problem
+#
+# We begin by introducing the function spaces
+#
+# $$
+#     V_h^g &:= \left\{v \in H(\textnormal{div}; \Omega);
+#     v|_K \in V_h(K) \; \forall K \in \mathcal{T}, v \cdot n = g \cdot n
+#     \textnormal{ on } \partial \Omega \right\} \\
+#     Q_h &:= \left\{q \in L^2_0(\Omega);
+#     q|_K \in Q_h(K) \; \forall K \in \mathcal{T} \right\}.
+# $$
+#
+# The local spaces $V_h(K)$ and $Q_h(K)$ should satisfy
+#
+# $$
+#     \nabla \cdot V_h(K) \subseteq Q_h(K),
+# $$
+# in order for mass to be conserved exactly. Suitable choices on
+# affine simplex cells include
+#
+# $$
+#     V_h(K) := \mathbb{RT}_k(K) \textnormal{ and }
+#     Q_h(K) := \mathbb{P}_k(K),
+# $$
+#
+# or
+#
+# $$
+#     V_h(K) := \mathbb{BDM}_k(K) \textnormal{ and }
+#     Q_h(K) := \mathbb{P}_{k-1}(K).
+# $$
+#
+# Let two cells $K^+$ and $K^-$ share a facet $F$. The trace of a
+# piecewise smooth vector valued function $\phi$ on F taken approaching
+# from inside $K^+$ (resp. $K^-$) is denoted $\phi^{+}$ (resp.
+# $\phi^-$). We now introduce the average
+# $\renewcommand{\avg}[1]{\left\{\!\!\left\{#1\right\}\!\!\right\}}$
+#
+# $$
+#     \avg{\phi} = \frac{1}{2} \left(\phi^+ + \phi^-\right)
+#
+# $$
+#
+# and jump $\renewcommand{\jump}[1]{[\![ #1 ]\!]}$
+#
+# $$
+#     \jump{\phi} = \phi^+ \otimes n^+ + \phi^- \otimes n^-,
+# $$
+#
+# operators, where $n$ denotes the outward unit normal to $\partial K$.
+# Finally, let the upwind flux of $\phi$ with respect to a vector field
+# $\psi$ be defined as
+#
+# $$
+#     \hat{\phi}^\psi :=
+#     \begin{cases}
+#         \lim_{\epsilon \downarrow 0} \phi(x - \epsilon \psi(x)), \;
+#         x \in \partial K \setminus \Gamma^\psi, \\
+#         0, \qquad \qquad \qquad \qquad x \in \partial K \cap \Gamma^\psi,
+#     \end{cases}
+# $$
+#
+# where $\Gamma^\psi = \left\{x \in \Gamma; \; \psi(x) \cdot n(x) < 0\right\}$.
+#
+# The semi-discrete version problem (in dimensionless form) is: find
+# $(u_h, p_h) \in V_h^{u_D} \times Q_h$ such that
+#
+# $$
+#     \int_\Omega \partial_t u_h \cdot v + a_h(u_h, v_h) + c_h(u_h; u_h, v_h)
+#     + b_h(v_h, p_h) &= \int_\Omega f \cdot v_h + L_{a_h}(v_h) + L_{c_h}(v_h)
+#      \quad \forall v_h \in V_h^0, \\
+#     b_h(u_h, q_h) &= 0 \quad \forall q_h \in Q_h,
+# $$
+#
+# where
+# $\renewcommand{\sumK}[0]{\sum_{K \in \mathcal{T}_h}}$
+# $\renewcommand{\sumF}[0]{\sum_{F \in \mathcal{F}_h}}$
+#
+# $$
+#     a_h(u, v) &= Re^{-1} \left(\sumK \int_K \nabla u : \nabla v
+#     - \sumF \int_F \avg{\nabla u} : \jump{v}
+#     - \sumF \int_F \avg{\nabla v} : \jump{u} \\
+#     + \sumF \int_F \frac{\alpha}{h_K} \jump{u} : \jump{v}\right), \\
+#     c_h(w; u, v) &= - \sumK \int_K u \cdot \nabla \cdot (v \otimes w)
+#     + \sumK \int_{\partial_K} w \cdot n \hat{u}^{w} \cdot v, \\
+#     L_{a_h}(v_h) &= Re^{-1} \left(- \int_{\partial \Omega} u_D \otimes n :
+#     \nabla_h v_h + \frac{\alpha}{h} u_D \otimes n : v_h \otimes n \right), \\
+#     L_{c_h}(v_h) &= - \int_{\partial \Omega} u_D \cdot n \hat{u}_D \cdot v_h, \\
+#     b_h(v, q) &= - \int_K \nabla \cdot v q.
+# $$
+#
+#
+# ## Implementation
+#
+# We begin by importing the required modules and functions
+
+# +
+import numpy as np
+from mpi4py import MPI
+from petsc4py import PETSc
+from ufl import (CellDiameter, FacetNormal, TestFunction, TrialFunction, avg,
+                 conditional, div, dot, dS, ds, dx, grad, gt, inner, outer)
+
+from dolfinx import fem, io, mesh
+
+if np.issubdtype(PETSc.ScalarType, np.complexfloating):
+    print("Demo should only be executed with DOLFINx real mode")
+    exit(0)
+# -
+
+
+# We also define some helper functions that will be used later
+
+
+# +
+def norm_L2(comm, v):
+    """Compute the L2(Ω)-norm of v"""
+    return np.sqrt(comm.allreduce(
+        fem.assemble_scalar(fem.form(inner(v, v) * dx)), op=MPI.SUM))
+
+
+def domain_average(msh, v):
+    """Compute the average of a function over the domain"""
+    vol = msh.comm.allreduce(
+        fem.assemble_scalar(fem.form(fem.Constant(msh, 1.0) * dx)), op=MPI.SUM)
+    return (1 / vol) * msh.comm.allreduce(fem.assemble_scalar(fem.form(v * dx)), op=MPI.SUM)
+
+
+def u_e_expr(x):
+    """Expression for the exact velocity solution to Kovasznay flow"""
+    return np.vstack((1 - np.exp(
+        (Re / 2 - np.sqrt(Re**2 / 4 + 4 * np.pi**2)) * x[0]) * np.cos(2 * np.pi * x[1]),
+        (Re / 2 - np.sqrt(Re**2 / 4 + 4 * np.pi**2)) / (2 * np.pi) * np.exp(
+            (Re / 2 - np.sqrt(Re**2 / 4 + 4 * np.pi**2)) * x[0])
+        * np.sin(2 * np.pi * x[1])))
+
+
+def p_e_expr(x):
+    """Expression for the exact pressure solution to Kovasznay flow"""
+    return (1 / 2) * (1 - np.exp(2 * (Re / 2 - np.sqrt(Re**2 / 4 + 4 * np.pi**2)) * x[0]))
+
+
+def f_expr(x):
+    """Expression for the applied force"""
+    return np.vstack((np.zeros_like(x[0]),
+                      np.zeros_like(x[0])))
+
+
+def boundary_marker(x):
+    return np.isclose(x[0], 0.0) | np.isclose(x[0], 1.0) | np.isclose(x[1], 0.0) | np.isclose(x[1], 1.0)
+# -
+
+# We define some simulation parameters
+
+
+n = 16
+num_time_steps = 25
+t_end = 10
+Re = 25  # Reynolds Number
+k = 1  # Polynomial degree
+
+# Next, we create a mesh and the required functions spaces over it.
+# Since the velocity uses an $H(\textnormal{div})$-conforming function
+# space, we also create a vector valued discontinuous Lagrange space to
+# interpolate into for artifact free visualisation.
+
+# +
+msh = mesh.create_unit_square(MPI.COMM_WORLD, n, n)
+
+# Function spaces for the velocity and for the pressure
+V = fem.FunctionSpace(msh, ("Raviart-Thomas", k + 1))
+Q = fem.FunctionSpace(msh, ("Discontinuous Lagrange", k))
+
+# Funcion space for visualising the velocity field
+W = fem.VectorFunctionSpace(msh, ("Discontinuous Lagrange", k + 1))
+
+# Define trial and test functions
+
+u, v = TrialFunction(V), TestFunction(V)
+p, q = TrialFunction(Q), TestFunction(Q)
+
+delta_t = fem.Constant(msh, t_end / num_time_steps)
+alpha = fem.Constant(msh, 6.0 * k**2)
+
+h = CellDiameter(msh)
+n = FacetNormal(msh)
+
+
+def jump(phi, n):
+    return outer(phi("+"), n("+")) + outer(phi("-"), n("-"))
+# -
+
+# We solve the Stokes problem for the initial condition, omitting the
+# convective term:
+
+
+# +
+a_00 = (1 / Re) * (inner(grad(u), grad(v)) * dx
+                   - inner(avg(grad(u)), jump(v, n)) * dS
+                   - inner(jump(u, n), avg(grad(v))) * dS
+                   + alpha / avg(h) * inner(jump(u, n), jump(v, n)) * dS
+                   - inner(grad(u), outer(v, n)) * ds
+                   - inner(outer(u, n), grad(v)) * ds
+                   + alpha / h * inner(outer(u, n), outer(v, n)) * ds)
+a_01 = - inner(p, div(v)) * dx
+a_10 = - inner(div(u), q) * dx
+
+a = fem.form([[a_00, a_01],
+              [a_10, None]])
+
+f = fem.Function(W)
+u_D = fem.Function(V)
+u_D.interpolate(u_e_expr)
+L_0 = inner(f, v) * dx + (1 / Re) * (- inner(outer(u_D, n), grad(v)) * ds
+                                     + (alpha / h) * inner(outer(u_D, n), outer(v, n)) * ds)
+L_1 = inner(fem.Constant(msh, 0.0), q) * dx
+L = fem.form([L_0, L_1])
+
+# Boundary conditions
+boundary_facets = mesh.locate_entities_boundary(msh, msh.topology.dim - 1, boundary_marker)
+boundary_vel_dofs = fem.locate_dofs_topological(V, msh.topology.dim - 1, boundary_facets)
+bc_u = fem.dirichletbc(u_D, boundary_vel_dofs)
+bcs = [bc_u]
+
+# Assemble Stokes problem
+A = fem.petsc.assemble_matrix_block(a, bcs=bcs)
+A.assemble()
+b = fem.petsc.assemble_vector_block(L, a, bcs=bcs)
+
+# Create and configure solver
+ksp = PETSc.KSP().create(msh.comm)
+ksp.setOperators(A)
+ksp.setType("preonly")
+ksp.getPC().setType("lu")
+ksp.getPC().setFactorSolverType("mumps")
+opts = PETSc.Options()
+opts["mat_mumps_icntl_14"] = 80  # Increase MUMPS working memory
+opts["mat_mumps_icntl_24"] = 1  # Option to support solving a singular matrix (pressure nullspace)
+opts["mat_mumps_icntl_25"] = 0  # Option to support solving a singular matrix (pressure nullspace)
+opts["ksp_error_if_not_converged"] = 1
+ksp.setFromOptions()
+
+# Solve Stokes for initial condition
+x = A.createVecRight()
+try:
+    ksp.solve(b, x)
+except PETSc.Error as e:
+    if e.ierr == 92:
+        print("The required PETSc solver/preconditioner is not available. Exiting.")
+        print(e)
+        exit(0)
+    else:
+        raise e
+
+
+# Split the solution
+u_h = fem.Function(V)
+p_h = fem.Function(Q)
+p_h.name = "p"
+offset = V.dofmap.index_map.size_local * V.dofmap.index_map_bs
+u_h.x.array[:offset] = x.array_r[:offset]
+u_h.x.scatter_forward()
+p_h.x.array[:(len(x.array_r) - offset)] = x.array_r[offset:]
+p_h.x.scatter_forward()
+# Subtract the average of the pressure since it is only determined
+# up to a constant
+p_h.x.array[:] -= domain_average(msh, p_h)
+
+u_vis = fem.Function(W)
+u_vis.name = "u"
+u_vis.interpolate(u_h)
+
+# Write initial condition to file
+t = 0.0
+try:
+    u_file = io.VTXWriter(msh.comm, "u.bp", [u_vis._cpp_object])
+    p_file = io.VTXWriter(msh.comm, "p.bp", [p_h._cpp_object])
+    u_file.write(t)
+    p_file.write(t)
+except AttributeError:
+    print("File output requires ADIOS2.")
+
+# Create function to store solution and previous time step
+u_n = fem.Function(V)
+u_n.x.array[:] = u_h.x.array
+# -
+
+# Now we add the time stepping and convective terms
+
+# +
+lmbda = conditional(gt(dot(u_n, n), 0), 1, 0)
+u_uw = lmbda("+") * u("+") + lmbda("-") * u("-")
+a_00 += inner(u / delta_t, v) * dx - \
+    inner(u, div(outer(v, u_n))) * dx + \
+    inner((dot(u_n, n))("+") * u_uw, v("+")) * dS + \
+    inner((dot(u_n, n))("-") * u_uw, v("-")) * dS + \
+    inner(dot(u_n, n) * lmbda * u, v) * ds
+a = fem.form([[a_00, a_01],
+              [a_10, None]])
+
+L_0 += inner(u_n / delta_t, v) * dx - inner(dot(u_n, n) * (1 - lmbda) * u_D, v) * ds
+L = fem.form([L_0, L_1])
+
+# Time stepping loop
+for n in range(num_time_steps):
+    t += delta_t.value
+
+    A.zeroEntries()
+    fem.petsc.assemble_matrix_block(A, a, bcs=bcs)  # type: ignore
+    A.assemble()
+
+    with b.localForm() as b_loc:
+        b_loc.set(0)
+    fem.petsc.assemble_vector_block(b, L, a, bcs=bcs)  # type: ignore
+
+    # Compute solution
+    ksp.solve(b, x)
+
+    u_h.x.array[:offset] = x.array_r[:offset]
+    u_h.x.scatter_forward()
+    p_h.x.array[:(len(x.array_r) - offset)] = x.array_r[offset:]
+    p_h.x.scatter_forward()
+    p_h.x.array[:] -= domain_average(msh, p_h)
+
+    u_vis.interpolate(u_h)
+
+    # Write to file
+    try:
+        u_file.write(t)
+        p_file.write(t)
+    except NameError:
+        pass
+
+    # Update u_n
+    u_n.x.array[:] = u_h.x.array
+
+try:
+    u_file.close()
+    p_file.close()
+except NameError:
+    pass
+# -
+
+# Now we compare the computed solution to the exact solution
+
+# +
+# Function spaces for exact velocity and pressure
+V_e = fem.VectorFunctionSpace(msh, ("Lagrange", k + 3))
+Q_e = fem.FunctionSpace(msh, ("Lagrange", k + 2))
+
+u_e = fem.Function(V_e)
+u_e.interpolate(u_e_expr)
+
+p_e = fem.Function(Q_e)
+p_e.interpolate(p_e_expr)
+
+# Compute errors
+e_u = norm_L2(msh.comm, u_h - u_e)
+e_div_u = norm_L2(msh.comm, div(u_h))
+
+# This scheme conserves mass exactly, so check this
+assert np.isclose(e_div_u, 0.0)
+p_e_avg = domain_average(msh, p_e)
+e_p = norm_L2(msh.comm, p_h - (p_e - p_e_avg))
+
+if msh.comm.rank == 0:
+    print(f"e_u = {e_u}")
+    print(f"e_div_u = {e_div_u}")
+    print(f"e_p = {e_p}")
+# -
diff --git a/python/doc/source/demos.rst b/python/doc/source/demos.rst
index 18fb1824084..9b2db3f5402 100644
--- a/python/doc/source/demos.rst
+++ b/python/doc/source/demos.rst
@@ -51,4 +51,5 @@ All demos
    demos/demo_pml.md
    demos/demo_half_loaded_waveguide.md
    demos/demo_axis.md
+   demos/demo_navier_stokes.md
 
diff --git a/python/dolfinx/fem/petsc.py b/python/dolfinx/fem/petsc.py
index 902bd6e64b7..e83e25fbf9d 100644
--- a/python/dolfinx/fem/petsc.py
+++ b/python/dolfinx/fem/petsc.py
@@ -248,24 +248,13 @@ def _assemble_vector_nest_vec(b: PETSc.Vec, L: typing.List[form_types], constant
 
 # FIXME: Revise this interface
 @functools.singledispatch
-def assemble_vector_block(L: typing.Any,
+def assemble_vector_block(L: typing.List[form_types],
                           a: typing.List[typing.List[form_types]],
                           bcs: typing.List[DirichletBCMetaClass] = [],
                           x0: typing.Optional[PETSc.Vec] = None,
                           scale: float = 1.0,
                           constants_L=None, coeffs_L=None,
                           constants_a=None, coeffs_a=None) -> PETSc.Vec:
-    return _assemble_vector_block_form(L, a, bcs, x0, scale, constants_L, coeffs_L, constants_a, coeffs_a)
-
-
-@assemble_vector_block.register(list)
-def _assemble_vector_block_form(L: typing.List[form_types],
-                                a: typing.List[typing.List[form_types]],
-                                bcs: typing.List[DirichletBCMetaClass] = [],
-                                x0: typing.Optional[PETSc.Vec] = None,
-                                scale: float = 1.0,
-                                constants_L=None, coeffs_L=None,
-                                constants_a=None, coeffs_a=None) -> PETSc.Vec:
     """Assemble linear forms into a monolithic vector. The vector is not
     finalised, i.e. ghost values are not accumulated.
 
@@ -372,18 +361,10 @@ def _assemble_matrix_mat(A: PETSc.Mat, a: form_types, bcs: typing.List[Dirichlet
 
 # FIXME: Revise this interface
 @functools.singledispatch
-def assemble_matrix_nest(a: typing.Any,
+def assemble_matrix_nest(a: typing.List[typing.List[form_types]],
                          bcs: typing.List[DirichletBCMetaClass] = [], mat_types=[],
                          diagonal: float = 1.0,
                          constants=None, coeffs=None) -> PETSc.Mat:
-    return _assemble_matrix_nest_form(a, bcs, diagonal, constants, coeffs)
-
-
-@assemble_matrix_nest.register(list)
-def _assemble_matrix_nest_form(a: typing.List[typing.List[form_types]],
-                               bcs: typing.List[DirichletBCMetaClass] = [], mat_types=[],
-                               diagonal: float = 1.0,
-                               constants=None, coeffs=None) -> PETSc.Mat:
     """Assemble bilinear forms into matrix"""
     A = _cpp.fem.petsc.create_matrix_nest(a, mat_types)
     _assemble_matrix_nest_mat(A, a, bcs, diagonal, constants, coeffs)
@@ -417,18 +398,10 @@ def _assemble_matrix_nest_mat(A: PETSc.Mat, a: typing.List[typing.List[form_type
 
 # FIXME: Revise this interface
 @functools.singledispatch
-def assemble_matrix_block(a: typing.Any,
+def assemble_matrix_block(a: typing.List[typing.List[form_types]],
                           bcs: typing.List[DirichletBCMetaClass] = [],
                           diagonal: float = 1.0,
-                          constants=None, coeffs=None) -> PETSc.Mat:
-    return _assemble_matrix_block_form(a, bcs, diagonal, constants, coeffs)
-
-
-@assemble_matrix_block.register(list)
-def _assemble_matrix_block_form(a: typing.List[typing.List[form_types]],
-                                bcs: typing.List[DirichletBCMetaClass] = [],
-                                diagonal: float = 1.0,
-                                constants=None, coeffs=None) -> PETSc.Mat:
+                          constants=None, coeffs=None) -> PETSc.Mat:  # type: ignore
     """Assemble bilinear forms into matrix"""
     A = _cpp.fem.petsc.create_matrix_block(a)
     return _assemble_matrix_block_mat(A, a, bcs, diagonal, constants, coeffs)