diff --git a/tests/extrusion/test_bdmc_extr.py b/tests/extrusion/test_bdmc_extr.py
new file mode 100644
index 0000000000..8091eb808a
--- /dev/null
+++ b/tests/extrusion/test_bdmc_extr.py
@@ -0,0 +1,46 @@
+import numpy as np
+from firedrake import *
+import pytest
+
+
+def project_bdmc(size, degree, family):
+    mesh = ExtrudedMesh(UnitIntervalMesh(2 ** size), layers=2 ** size)
+    x = SpatialCoordinate(mesh)
+
+    fs = FunctionSpace(mesh, family, degree)
+
+    f = Function(fs)
+    v = TestFunction(fs)
+
+    expr = as_vector((sin(2 * pi * x[0]) * sin(2 * pi * x[1]),
+                      cos(2 * pi * x[0]) * cos(2 * pi * x[1])))
+
+    solve(inner(f-expr, v) * dx(degree=8) == 0, f)
+
+    return np.sqrt(assemble(inner(f-expr, f-expr) * dx(degree=8)))
+
+
+@pytest.mark.parametrize(('testcase', 'convrate', 'degree'),
+                         [((3, 6), 1.9, 1),
+                          ((3, 6), 2.9, 2),
+                          ((3, 6), 3.9, 3),
+                          ((3, 6), 4.9, 4)])
+def test_bdmcf(testcase, convrate, degree):
+    start, end = testcase
+    l2err = np.zeros(end - start)
+    for ii in [i + start for i in range(len(l2err))]:
+        l2err[ii - start] = project_bdmc(ii, degree, "BDMCF")
+    assert (np.log2(l2err[:-1] / l2err[1:]) > convrate).all()
+
+
+@pytest.mark.parametrize(('testcase', 'convrate', 'degree'),
+                         [((3, 6), 1.9, 1),
+                          ((3, 6), 2.9, 2),
+                          ((3, 6), 3.9, 3),
+                          ((3, 6), 4.9, 4)])
+def test_bdmce(testcase, convrate, degree):
+    start, end = testcase
+    l2err = np.zeros(end - start)
+    for ii in [i + start for i in range(len(l2err))]:
+        l2err[ii - start] = project_bdmc(ii, degree, "BDMCE")
+    assert (np.log2(l2err[:-1] / l2err[1:]) > convrate).all()
diff --git a/tests/regression/test_bdmc.py b/tests/regression/test_bdmc.py
new file mode 100644
index 0000000000..6cd881feb1
--- /dev/null
+++ b/tests/regression/test_bdmc.py
@@ -0,0 +1,46 @@
+import numpy as np
+from firedrake import *
+import pytest
+
+
+def project_bdmc(size, degree, family):
+    mesh = UnitSquareMesh(2 ** size, 2 ** size, quadrilateral=True)
+    x = SpatialCoordinate(mesh)
+
+    fs = FunctionSpace(mesh, family, degree)
+
+    f = Function(fs)
+    v = TestFunction(fs)
+
+    expr = as_vector((sin(2 * pi * x[0]) * sin(2 * pi * x[1]),
+                      cos(2 * pi * x[0]) * cos(2 * pi * x[1])))
+
+    solve(inner(f-expr, v) * dx(degree=8) == 0, f)
+
+    return np.sqrt(assemble(inner(f-expr, f-expr) * dx(degree=8)))
+
+
+@pytest.mark.parametrize(('testcase', 'convrate', 'degree'),
+                         [((3, 6), 1.9, 1),
+                          ((3, 6), 2.9, 2),
+                          ((3, 6), 3.9, 3),
+                          ((3, 6), 4.9, 4)])
+def test_bdmcf(testcase, convrate, degree):
+    start, end = testcase
+    l2err = np.zeros(end - start)
+    for ii in [i + start for i in range(len(l2err))]:
+        l2err[ii - start] = project_bdmc(ii, degree, "BDMCF")
+    assert (np.log2(l2err[:-1] / l2err[1:]) > convrate).all()
+
+
+@pytest.mark.parametrize(('testcase', 'convrate', 'degree'),
+                         [((3, 6), 1.9, 1),
+                          ((3, 6), 2.9, 2),
+                          ((3, 6), 3.9, 3),
+                          ((3, 6), 4.9, 4)])
+def test_bdmce(testcase, convrate, degree):
+    start, end = testcase
+    l2err = np.zeros(end - start)
+    for ii in [i + start for i in range(len(l2err))]:
+        l2err[ii - start] = project_bdmc(ii, degree, "BDMCE")
+    assert (np.log2(l2err[:-1] / l2err[1:]) > convrate).all()
diff --git a/tests/regression/test_bdmc_riesz_map.py b/tests/regression/test_bdmc_riesz_map.py
new file mode 100644
index 0000000000..7b8c17809a
--- /dev/null
+++ b/tests/regression/test_bdmc_riesz_map.py
@@ -0,0 +1,70 @@
+from firedrake import *
+import pytest
+import numpy
+
+
+@pytest.fixture(scope='module', params=["div", "curl"])
+def problem(request):
+    return request.param
+
+
+@pytest.fixture(scope='module', params=[1, 2, 3])
+def degree(request):
+    return request.param
+
+
+sp = {"snes_type": "ksponly",
+      "ksp_type": "preonly",
+      "pc_type": "lu",
+      "pc_factor_mat_solver_type": "mumps",
+      "mat_mumps_icntl_14": 200}
+
+
+def error(N, problem, degree):
+    mesh = UnitSquareMesh(N, N, quadrilateral=True)
+
+    if problem == "div":
+        op = div
+        family = "BDMCF"
+    elif problem == "curl":
+        op = curl
+        family = "BDMCE"
+    else:
+        raise ValueError
+
+    V = FunctionSpace(mesh, family, degree)
+    u = Function(V)
+    v = TestFunction(V)
+
+    (x, y) = SpatialCoordinate(mesh)
+    u_ex = as_vector([sin(2*pi*x) * cos(2*pi*y),
+                      x * (1-x) * y * (1-y)])
+
+    if problem == "div":
+        f = u_ex - grad(div(u_ex))
+    else:
+        f = u_ex + curl(curl(u_ex))
+
+    F = inner(op(u), op(v))*dx + inner(u, v)*dx - inner(f, v)*dx
+    bc = DirichletBC(V, project(u_ex, V, solver_parameters=sp), "on_boundary")
+
+    solve(F == 0, u, bc, solver_parameters=sp)
+
+    err = errornorm(u_ex, u, "L2")
+    return err
+
+
+def test_bdmc_riesz_map(problem, degree):
+
+    errors = []
+    for N in [10, 20, 40]:
+        errors.append(error(N, problem, degree))
+
+    convergence_orders = lambda x: numpy.log2(numpy.array(x)[:-1] / numpy.array(x)[1:])
+    conv = convergence_orders(errors)
+    print("errors: ", errors)
+    print("convergence order: ", conv)
+
+    tol = 0.11
+
+    assert (conv > (degree + 1) - tol).all()