diff --git a/Jenkinsfile b/Jenkinsfile
index 7684a36f23..fc8446c4e4 100644
--- a/Jenkinsfile
+++ b/Jenkinsfile
@@ -26,7 +26,7 @@ pipeline {
         sh 'mkdir tmp'
         dir('tmp') {
           timestamps {
-            sh '../scripts/firedrake-install --disable-ssh --minimal-petsc --slepc --documentation-dependencies --install thetis --install gusto --install icepack --install pyadjoint --no-package-manager || (cat firedrake-install.log && /bin/false)'
+            sh '../scripts/firedrake-install --package-branch ufl mse --package-branch tsfc mse --package-branch fiat mse --package-branch FInAT mse --disable-ssh --minimal-petsc --slepc --documentation-dependencies --install thetis --install gusto --install icepack --install pyadjoint --no-package-manager || (cat firedrake-install.log && /bin/false)'
           }
         }
       }
diff --git a/requirements-git.txt b/requirements-git.txt
index d935e8dcd8..96651ca2f9 100644
--- a/requirements-git.txt
+++ b/requirements-git.txt
@@ -1,5 +1,5 @@
-git+https://github.com/firedrakeproject/ufl.git#egg=ufl
-git+https://github.com/firedrakeproject/fiat.git#egg=fiat
-git+https://github.com/FInAT/FInAT.git#egg=finat
-git+https://github.com/firedrakeproject/tsfc.git#egg=tsfc
+git+https://github.com/firedrakeproject/ufl.git@mse#egg=ufl
+git+https://github.com/firedrakeproject/fiat.git@mse#egg=fiat
+git+https://github.com/FInAT/FInAT.git@mse#egg=finat
+git+https://github.com/firedrakeproject/tsfc.git@mse#egg=tsfc
 git+https://github.com/OP2/PyOP2.git#egg=pyop2
diff --git a/tests/regression/test_dualmse.py b/tests/regression/test_dualmse.py
new file mode 100644
index 0000000000..333757c6bb
--- /dev/null
+++ b/tests/regression/test_dualmse.py
@@ -0,0 +1,365 @@
+import pytest
+import numpy as np
+from firedrake import *
+
+
+def get_complex(meshtype, hdegree, vdegree=None):
+
+    if meshtype == 'extruded':
+        vdegree = vdegree or hdegree
+        h1_horiz = FiniteElement("EGL", "interval", hdegree)
+        l2_horiz = FiniteElement("EGL-Edge L2", "interval", hdegree-1)
+        h1_vert = FiniteElement("EGL", "interval", vdegree)
+        l2_vert = FiniteElement("EGL-Edge L2", "interval", vdegree-1)
+        hcurlelem = HCurl(TensorProductElement(h1_horiz, l2_vert)) + HCurl(TensorProductElement(l2_horiz, h1_vert))
+        hdivelem = HDiv(TensorProductElement(h1_horiz, l2_vert)) + HDiv(TensorProductElement(l2_horiz, h1_vert))
+        l2elem = TensorProductElement(l2_horiz, l2_vert)
+        h1elem = TensorProductElement(h1_horiz, h1_vert)
+    else:
+        l2elem = FiniteElement('DQ L2', quadrilateral, hdegree-1, variant='dualmse')
+        hdivelem = FiniteElement('RTCF', quadrilateral, hdegree, variant='dualmse')
+        hcurlelem = FiniteElement('RTCE', quadrilateral, hdegree, variant='dualmse')
+        h1elem = FiniteElement('Q', quadrilateral, hdegree, variant='dualmse')
+
+    return h1elem, hcurlelem, hdivelem, l2elem
+
+
+def get_mesh(r, meshtype, meshd=None):
+
+    meshd = meshd or 1
+    if meshtype == 'spherical':
+        mesh = UnitCubedSphereMesh(refinement_level=r, degree=meshd)
+        x = SpatialCoordinate(mesh)
+        mesh.init_cell_orientations(x)
+    elif meshtype == 'planar':
+        mesh = UnitSquareMesh(2**r, 2**r, quadrilateral=True)
+    elif meshtype == 'extruded':
+        basemesh = UnitIntervalMesh(2**r)
+        mesh = ExtrudedMesh(basemesh, 2**r)
+    return mesh
+
+
+def helmholtz(r, meshtype, hdegree, vdegree=None, meshd=None):
+
+    mesh = get_mesh(r, meshtype, meshd=meshd)
+    h1elem, _, _, _ = get_complex(meshtype, hdegree, vdegree=vdegree)
+
+    V = FunctionSpace(mesh, h1elem)
+
+    # Define variational problem
+    lmbda = 1
+    u = TrialFunction(V)
+    v = TestFunction(V)
+    f = Function(V)
+    x = SpatialCoordinate(mesh)
+    if meshtype in ['planar', 'extruded']:
+        f.project((1+8*pi*pi)*cos(x[0]*pi*2)*cos(x[1]*pi*2))
+    elif meshtype == 'spherical':
+        f.project(x[0]*x[1]*x[2])
+    a = (dot(grad(v), grad(u)) + lmbda * v * u) * dx
+    L = f * v * dx
+
+    # Compute solution
+    assemble(a)
+    assemble(L)
+    sol = Function(V)
+    solve(a == L, sol, solver_parameters={'ksp_type': 'cg'})
+
+    # Analytical solution
+    if meshtype in ['planar', 'extruded']:
+        f.project(cos(x[0]*pi*2)*cos(x[1]*pi*2))
+    elif meshtype == 'spherical':
+        f.project(x[0]*x[1]*x[2]/13.0)
+    return sqrt(assemble(dot(sol - f, sol - f) * dx))
+
+# helmholtz on plane
+
+
+@pytest.mark.parametrize(('degree', 'threshold'),
+                         [(2, 2.9),
+                          (3, 3.9)])
+def test_firedrake_helmholtz_dualmse(degree, threshold):
+    diff = np.array([helmholtz(r, 'planar', degree) for r in range(3, 6)])
+    print("l2 error norms:", diff)
+    conv = np.log2(diff[:-1] / diff[1:])
+    print("convergence order:", conv)
+    assert (np.array(conv) > threshold).all()
+
+
+# helmholtz on sphere
+
+@pytest.mark.parametrize(('degree', 'md', 'threshold'),
+                         [(2, 1, 1.9),
+                          (2, 2, 2.9),
+                          (3, 1, 1.87),
+                          (3, 2, 2.9),
+                          (3, 3, 3.9)])
+def test_firedrake_helmholtz_dualmse_sphere(degree, md, threshold):
+    diff = np.array([helmholtz(r, 'spherical', degree, meshd=md) for r in range(2, 5)])
+    print("l2 error norms:", diff)
+    conv = np.log2(diff[:-1] / diff[1:])
+    print("convergence order:", conv)
+    assert (np.array(conv) > threshold).all()
+
+
+# helmholtz on extruded
+
+@pytest.mark.parametrize(('hdegree', 'vdegree', 'threshold'),
+                         [(2, 2, 2.9),
+                          (2, 3, 2.9),
+                          (3, 2, 2.9),
+                          (3, 3, 3.9)])
+def test_firedrake_helmholtz_dualmse_extruded(hdegree, vdegree, threshold):
+    diff = np.array([helmholtz(i, 'extruded', hdegree, vdegree=vdegree) for i in range(3, 6)])
+    print("l2 error norms:", diff)
+    conv = np.log2(diff[:-1] / diff[1:])
+    print("convergence order:", conv)
+    assert (np.array(conv) > threshold).all()
+
+
+def helmholtz_mixed(r, meshtype, hdegree, vdegree=None, meshd=None, useaction=False):
+
+    mesh = get_mesh(r, meshtype, meshd=meshd)
+    _, _, hdivelem, l2elem = get_complex(meshtype, hdegree, vdegree=vdegree)
+
+    V1 = FunctionSpace(mesh, hdivelem, name="V")
+    V2 = FunctionSpace(mesh, l2elem, name="P")
+    W = V1 * V2
+
+    # Define variational problem
+    lmbda = 1
+    u, p = TrialFunctions(W)
+    v, q = TestFunctions(W)
+    f = Function(V2)
+
+    x = SpatialCoordinate(mesh)
+    if meshtype in ['planar', 'extruded']:
+        f.project((1+8*pi*pi)*sin(x[0]*pi*2)*sin(x[1]*pi*2))
+    elif meshtype == 'spherical':
+        f.project(x[0]*x[1]*x[2])
+    a = (p*q - q*div(u) + lmbda*inner(v, u) + div(v)*p) * dx
+    L = f*q*dx
+
+    # Compute solution
+    sol = Function(W)
+
+    if useaction:
+        system = action(a, sol) - L == 0
+    else:
+        system = a == L
+
+    # Block system is:
+    # V Ct
+    # Ch P
+    # Eliminate V by forming a schur complement
+    solve(system, sol, solver_parameters={'pc_type': 'fieldsplit',
+                                          'pc_fieldsplit_type': 'schur',
+                                          'ksp_type': 'cg',
+                                          'pc_fieldsplit_schur_fact_type': 'FULL',
+                                          'fieldsplit_V_ksp_type': 'cg',
+                                          'fieldsplit_P_ksp_type': 'cg'})
+
+    if meshtype in ['planar', 'extruded']:
+        f.project(sin(x[0]*pi*2)*sin(x[1]*pi*2))
+        return sqrt(assemble(dot(sol[2] - f, sol[2] - f) * dx))
+    elif meshtype == 'spherical':
+        _, u = sol.split()
+        f.project(x[0]*x[1]*x[2]/13.0)
+        return errornorm(f, u, degree_rise=0)
+
+
+# helmholtz mixed on plane
+
+@pytest.mark.parametrize(('degree', 'action', 'threshold'),
+                         [(2, False, 2.9),
+                          (3, False, 3.9),
+                          (2, True, 2.9)])
+def test_firedrake_helmholtz_mixed_dualmse(degree, action, threshold):
+    diff = np.array([helmholtz_mixed(r, 'planar', degree, useaction=action) for r in range(3, 6)])
+    print("l2 error norms:", diff)
+    conv = np.log2(diff[:-1] / diff[1:])
+    print("convergence order:", conv)
+    assert (np.array(conv) > threshold).all()
+
+
+# helmholtz mixed on sphere
+
+@pytest.mark.parametrize(('degree', 'md', 'action', 'threshold'),
+                         [(2, 1, False, 1.9),
+                          (2, 2, False, 2.9),
+                          (3, 1, False, 1.87),
+                          (3, 2, False, 2.9),
+                          (3, 3, False, 3.9),
+                          (2, 2, True, 2.9)])
+def test_firedrake_helmholtz_mixed_dualmse_sphere(degree, md, action, threshold):
+    diff = np.array([helmholtz_mixed(r, 'spherical', degree, meshd=md, useaction=action) for r in range(2, 5)])
+    print("l2 error norms:", diff)
+    conv = np.log2(diff[:-1] / diff[1:])
+    print("convergence order:", conv)
+    assert (np.array(conv) > threshold).all()
+
+
+# helmholtz mixed on extruded
+
+@pytest.mark.parametrize(('hdegree', 'vdegree', 'threshold', 'action'),
+                         [(2, 2, 2.9, False),
+                          (2, 3, 2.9, False),
+                          (3, 2, 2.9, False),
+                          (3, 3, 3.9, False),
+                          (2, 2, 2.9, True)])
+def test_firedrake_helmholtz_mixed_dualmse_extruded(hdegree, vdegree, threshold, action):
+    diff = np.array([helmholtz_mixed(i, 'extruded', hdegree, vdegree=vdegree, useaction=action) for i in range(3, 6)])
+    print("l2 error norms:", diff)
+    conv = np.log2(diff[:-1] / diff[1:])
+    print("convergence order:", conv)
+    assert (np.array(conv) > threshold).all()
+
+
+# facet integrals- interior and exterior
+# tested using an interior penalty formulation for the Helmholtz using DG/DQ L2 elements
+
+
+def laplacian_IP(r, degree, meshd, meshtype):
+
+    dIF = dS
+    dEF = ds
+    if meshtype == 'extruded':
+        dIF = dS_v + dS_h
+        dEF = ds_t + ds_b + ds_v
+
+    mesh = get_mesh(r, meshtype, meshd=meshd)
+    _, _, _, l2elem = get_complex(meshtype, degree+1, vdegree=degree+1)
+
+    V = FunctionSpace(mesh, l2elem)
+
+    # Define variational problem
+    u = TrialFunction(V)
+    v = TestFunction(V)
+    f = Function(V)
+
+    x = SpatialCoordinate(mesh)
+    if meshtype in ['planar', 'extruded']:
+        f.project((1+8*pi*pi)*sin(x[0]*pi*2)*sin(x[1]*pi*2))
+    elif meshtype == 'spherical':
+        f.project(x[0]*x[1]*x[2])
+
+    FA = FacetArea(mesh)
+    CV = CellVolume(mesh)
+    ddx = CV/FA
+    ddx_avg = (ddx('+') + ddx('-'))/2.
+    alpha = Constant(4. * degree * (degree + 1.))
+    gamma = Constant(8. * degree * (degree + 1.))
+    penalty_int = alpha / ddx_avg
+    penalty_ext = gamma / ddx
+
+    n = FacetNormal(mesh)
+    aV = (inner(grad(u), grad(v)) + inner(u, v)) * dx  # volume term
+    aIF = (inner(jump(u, n), jump(v, n)) * penalty_int - inner(avg(grad(u)), jump(v, n)) - inner(avg(grad(v)), jump(u, n))) * dIF  # interior facet term
+    aEF = (inner(u, v) * penalty_ext - inner(grad(u), v*n) - inner(grad(v), u*n)) * dEF  # exterior facet term
+    a = aV + aEF + aIF
+    L = f*v*dx
+
+    # Compute solution
+    sol = Function(V)
+
+    solve(a == L, sol, solver_parameters={'pc_type': 'ilu',
+                                          'ksp_type': 'lgmres'})
+
+    # Analytical solution
+    if meshtype in ['planar', 'extruded']:
+        f.project(sin(x[0]*pi*2)*sin(x[1]*pi*2))
+    elif meshtype == 'spherical':
+        f.project(x[0]*x[1]*x[2]/13.0)
+    return sqrt(assemble(dot(sol - f, sol - f) * dx))
+
+
+@pytest.mark.parametrize(('degree', 'meshd', 'meshtype', 'threshold'),
+                         [(2, 1, 'planar', 2.9),
+                          (3, 1, 'planar', 3.8),
+                          (2, 1, 'extruded', 2.9),
+                          (3, 1, 'extruded', 3.8),
+                          (2, 2, 'spherical', 2.9),
+                          (3, 3, 'spherical', 3.9)])
+def test_firedrake_laplacian_IP_dualmse(degree, meshd, meshtype, threshold):
+    diff = np.array([laplacian_IP(i, degree, meshd, meshtype) for i in range(2, 5)])
+    print("l2 error norms:", diff)
+    conv = np.log2(diff[:-1] / diff[1:])
+    print("convergence order:", conv)
+    assert (np.array(conv) > threshold).all()
+
+
+def vector_laplace(r, meshtype, hdegree, vdegree=None):
+    vdegree = vdegree or hdegree
+
+    mesh = get_mesh(r, meshtype)
+    h1elem, hcurlelem, _, _ = get_complex(meshtype, hdegree, vdegree=vdegree)
+    h1elem_high, _, _, _ = get_complex(meshtype, hdegree+1, vdegree=vdegree+1)
+
+    # spaces for calculation
+    V0 = FunctionSpace(mesh, h1elem)
+    V1 = FunctionSpace(mesh, hcurlelem)
+    V = V0*V1
+
+    # spaces to store 'analytic' functions
+    W0 = FunctionSpace(mesh, h1elem_high)
+    W1 = VectorFunctionSpace(mesh, h1elem_high)
+
+    # constants
+    k = 1.0
+    l = 2.0
+
+    xs = SpatialCoordinate(mesh)
+    f_expr = as_vector([pi*pi*(k*k + l*l)*sin(k*pi*xs[0])*cos(l*pi*xs[1]), pi*pi*(k*k + l*l)*cos(k*pi*xs[0])*sin(l*pi*xs[1])])
+    exact_s_expr = -(k+l)*pi*cos(k*pi*xs[0])*cos(l*pi*xs[1])
+    exact_u_expr = as_vector([sin(k*pi*xs[0])*cos(l*pi*xs[1]), cos(k*pi*xs[0])*sin(l*pi*xs[1])])
+
+    f = Function(W1).project(f_expr)
+    exact_s = Function(W0).project(exact_s_expr)
+    exact_u = Function(W1).project(exact_u_expr)
+
+    sigma, u = TrialFunctions(V)
+    tau, v = TestFunctions(V)
+    a = (sigma*tau - dot(u, grad(tau)) + dot(grad(sigma), v) + dot(curl(u), curl(v)))*dx
+    L = dot(f, v)*dx
+
+    out = Function(V)
+
+    # preconditioner for H1 x H(curl)
+    aP = (dot(grad(sigma), grad(tau)) + sigma*tau + dot(curl(u), curl(v)) + dot(u, v))*dx
+
+    solve(a == L, out, Jp=aP,
+          solver_parameters={'pc_type': 'fieldsplit',
+                             'pc_fieldsplit_type': 'additive',
+                             'fieldsplit_0_pc_type': 'lu',
+                             'fieldsplit_1_pc_type': 'lu',
+                             'ksp_monitor': None})
+
+    out_s, out_u = out.split()
+
+    return (sqrt(assemble(dot(out_u - exact_u, out_u - exact_u)*dx)),
+            sqrt(assemble((out_s - exact_s)*(out_s - exact_s)*dx)))
+
+
+@pytest.mark.parametrize(('degree', 'threshold'),
+                         [(2, 1.9),
+                          (3, 2.9),
+                          (2, 1.9)])
+def test_firedrake_helmholtz_vector_dualmse(degree, threshold):
+    diff = np.array([vector_laplace(r, 'planar', degree) for r in range(2, 5)])
+    print("l2 error norms:", diff)
+    conv = np.log2(diff[:-1] / diff[1:])
+    print("convergence order:", conv)
+    assert (np.array(conv) > threshold).all()
+
+
+@pytest.mark.parametrize(('hdegree', 'vdegree', 'threshold'),
+                         [(2, 2, 1.9),
+                          (2, 3, 1.9),
+                          (3, 2, 1.9),
+                          (3, 3, 2.9)])
+def test_firedrake_helmholtz_vector_dualmse_extruded(hdegree, vdegree, threshold):
+    diff = np.array([vector_laplace(i, 'extruded', hdegree, vdegree=vdegree) for i in range(2, 5)])
+    print("l2 error norms:", diff)
+    conv = np.log2(diff[:-1, :] / diff[1:, :])
+    print("convergence order:", conv)
+    assert (np.array(conv) > threshold).all()
diff --git a/tests/regression/test_l2pullbacks.py b/tests/regression/test_l2pullbacks.py
new file mode 100644
index 0000000000..3c92d45bc5
--- /dev/null
+++ b/tests/regression/test_l2pullbacks.py
@@ -0,0 +1,242 @@
+import pytest
+import numpy as np
+from firedrake import *
+
+
+def get_complex(family, hdegree, vdegree=None):
+
+    if family == 'RT':
+        l2elem = FiniteElement('DG L2', triangle, hdegree-1)
+        hdivelem = FiniteElement('RT', triangle, hdegree)
+    elif family == 'BDM':
+        l2elem = FiniteElement('DG L2', triangle, hdegree-1)
+        hdivelem = FiniteElement('BDM', triangle, hdegree)
+    elif family == 'BDFM':
+        l2elem = FiniteElement('DG L2', triangle, hdegree-1)
+        hdivelem = FiniteElement('BDFM', triangle, hdegree)
+    elif family == 'RTCF':
+        l2elem = FiniteElement('DQ L2', quadrilateral, hdegree-1)
+        hdivelem = FiniteElement('RTCF', quadrilateral, hdegree)
+    elif family == 'ext':
+        vdegree = vdegree or hdegree
+        h1_horiz = FiniteElement("CG", "interval", hdegree)
+        l2_horiz = FiniteElement("DG L2", "interval", hdegree-1)
+        h1_vert = FiniteElement("CG", "interval", vdegree)
+        l2_vert = FiniteElement("DG L2", "interval", vdegree-1)
+        hdivelem = HDiv(TensorProductElement(h1_horiz, l2_vert)) + HDiv(TensorProductElement(l2_horiz, h1_vert))
+        l2elem = TensorProductElement(l2_horiz, l2_vert)
+    return hdivelem, l2elem
+
+
+def get_mesh(r, meshtype, meshd=None):
+
+    meshd = meshd or 1
+    if meshtype == 'spherical-quad':
+        mesh = UnitCubedSphereMesh(refinement_level=r, degree=meshd)
+        x = SpatialCoordinate(mesh)
+        mesh.init_cell_orientations(x)
+    elif meshtype == 'spherical-tri':
+        mesh = UnitIcosahedralSphereMesh(refinement_level=r, degree=meshd)
+        x = SpatialCoordinate(mesh)
+        mesh.init_cell_orientations(x)
+    elif meshtype == 'planar-quad':
+        mesh = UnitSquareMesh(2**r, 2**r, quadrilateral=True)
+    elif meshtype == 'planar-tri':
+        mesh = UnitSquareMesh(2**r, 2**r, quadrilateral=False)
+    elif meshtype == 'extruded':
+        basemesh = UnitIntervalMesh(2**r)
+        mesh = ExtrudedMesh(basemesh, 2**r)
+    return mesh
+
+
+def helmholtz_mixed(r, meshtype, family, hdegree, vdegree=None, meshd=None, useaction=False):
+
+    mesh = get_mesh(r, meshtype, meshd=meshd)
+    hdivelem, l2elem = get_complex(family, hdegree, vdegree=vdegree)
+
+    V1 = FunctionSpace(mesh, hdivelem, name="V")
+    V2 = FunctionSpace(mesh, l2elem, name="P")
+    W = V1 * V2
+
+    # Define variational problem
+    lmbda = 1
+    u, p = TrialFunctions(W)
+    v, q = TestFunctions(W)
+    f = Function(V2)
+
+    x = SpatialCoordinate(mesh)
+    if meshtype in ['planar-quad', 'planar-tri', 'extruded']:
+        f.project((1+8*pi*pi)*sin(x[0]*pi*2)*sin(x[1]*pi*2))
+    elif meshtype in ['spherical-quad', 'spherical-tri']:
+        f.project(x[0]*x[1]*x[2])
+    a = (p*q - q*div(u) + lmbda*inner(v, u) + div(v)*p) * dx
+    L = f*q*dx
+
+    # Compute solution
+    sol = Function(W)
+
+    if useaction:
+        system = action(a, sol) - L == 0
+    else:
+        system = a == L
+
+    # Block system is:
+    # V Ct
+    # Ch P
+    # Eliminate V by forming a schur complement
+    solve(system, sol, solver_parameters={'pc_type': 'fieldsplit',
+                                          'pc_fieldsplit_type': 'schur',
+                                          'ksp_type': 'cg',
+                                          'pc_fieldsplit_schur_fact_type': 'FULL',
+                                          'fieldsplit_V_ksp_type': 'cg',
+                                          'fieldsplit_P_ksp_type': 'cg'})
+
+    # Analytical solution
+
+    if meshtype in ['planar-quad', 'planar-tri', 'extruded']:
+        f.project(sin(x[0]*pi*2)*sin(x[1]*pi*2))
+        return sqrt(assemble(dot(sol[2] - f, sol[2] - f) * dx))
+    elif meshtype in ['spherical-quad', 'spherical-tri']:
+        _, u = sol.split()
+        f.project(x[0]*x[1]*x[2]/13.0)
+        return errornorm(f, u, degree_rise=0)
+
+
+# helmholtz mixed on plane
+
+@pytest.mark.parametrize(('family', 'degree', 'celltype', 'action', 'threshold'),
+                         [('RT', 1, 'tri', False, 1.9),
+                          ('RT', 2, 'tri', False, 2.9),
+                          ('BDM', 1, 'tri', False, 1.87),
+                          ('BDM', 2, 'tri', False, 2.9),
+                          ('BDFM', 2, 'tri', False, 2.9),
+                          ('RTCF', 1, 'quad', False, 1.9),
+                          ('RTCF', 2, 'quad', False, 2.9),
+                          ('RT', 2, 'tri', True, 2.9),
+                          ('BDM', 2, 'tri', True, 2.9)])
+def test_firedrake_helmholtz_mixed_l2pullbacks(family, degree, celltype, action, threshold):
+    diff = np.array([helmholtz_mixed(r, 'planar-' + celltype, family, degree, useaction=action) for r in range(3, 6)])
+    print("l2 error norms:", diff)
+    conv = np.log2(diff[:-1] / diff[1:])
+    print("convergence order:", conv)
+    assert (np.array(conv) > threshold).all()
+
+
+# helmholtz mixed on sphere
+
+@pytest.mark.parametrize(('family', 'degree', 'celltype', 'md', 'action', 'threshold'),
+                         [('RT', 1, 'tri', 1, False, 1.9),
+                          ('RT', 2, 'tri', 1, False, 1.9),
+                          ('RT', 2, 'tri', 2, False, 2.9),
+                          ('BDM', 1, 'tri', 1, False, 1.88),
+                          ('BDM', 2, 'tri', 1, False, 1.9),
+                          ('BDM', 2, 'tri', 2, False, 2.9),
+                          ('BDFM', 2, 'tri', 1, False, 1.9),
+                          ('BDFM', 2, 'tri', 2, False, 2.9),
+                          ('RTCF', 1, 'quad', 1, False, 1.67),
+                          ('RTCF', 2, 'quad', 1, False, 1.9),
+                          ('RTCF', 2, 'quad', 2, False, 2.9),
+                          ('RT', 2, 'tri', 1, True, 1.9),
+                          ('BDM', 2, 'tri', 1, True, 1.9)])
+def test_firedrake_helmholtz_mixed_l2pullbacks_sphere(family, degree, celltype, md, action, threshold):
+    diff = np.array([helmholtz_mixed(r, 'spherical-' + celltype, family, degree, meshd=md, useaction=action) for r in range(2, 5)])
+    print("l2 error norms:", diff)
+    conv = np.log2(diff[:-1] / diff[1:])
+    print("convergence order:", conv)
+    assert (np.array(conv) > threshold).all()
+
+
+# helmholtz mixed on extruded
+
+@pytest.mark.parametrize(('hdegree', 'vdegree', 'threshold', 'action'),
+                         [(1, 1, 1.9, False),
+                          (2, 1, 1.9, False),
+                          (1, 2, 1.9, False),
+                          (2, 2, 2.9, False),
+                          (1, 1, 1.9, True)])
+def test_firedrake_helmholtz_mixed_l2pullbacks_extruded(hdegree, vdegree, threshold, action):
+    diff = np.array([helmholtz_mixed(i, 'extruded', 'ext', hdegree, vdegree=vdegree, useaction=action) for i in range(3, 6)])
+    print("l2 error norms:", diff)
+    conv = np.log2(diff[:-1] / diff[1:])
+    print("convergence order:", conv)
+    assert (np.array(conv) > threshold).all()
+
+# facet integrals- interior and exterior
+# tested using an interior penalty formulation for the Helmholtz using DG/DQ L2 elements
+
+
+def laplacian_IP(r, degree, meshd, meshtype):
+
+    dIF = dS
+    dEF = ds
+    if meshtype == 'extruded':
+        dIF = dS_v + dS_h
+        dEF = ds_t + ds_b + ds_v
+        hdivelem, l2elem = get_complex('ext', degree+1, vdegree=degree+1)
+    elif meshtype in ['spherical-quad', 'planar-quad']:
+        hdivelem, l2elem = get_complex('RTCF', degree+1)
+    elif meshtype in ['spherical-tri', 'planar-tri']:
+        hdivelem, l2elem = get_complex('RT', degree+1)
+
+    mesh = get_mesh(r, meshtype, meshd=meshd)
+
+    V = FunctionSpace(mesh, l2elem)
+
+    # Define variational problem
+    u = TrialFunction(V)
+    v = TestFunction(V)
+    f = Function(V)
+
+    x = SpatialCoordinate(mesh)
+    if meshtype in ['planar-quad', 'planar-tri', 'extruded']:
+        f.project((1+8*pi*pi)*sin(x[0]*pi*2)*sin(x[1]*pi*2))
+    elif meshtype in ['spherical-quad', 'spherical-tri']:
+        f.project(x[0]*x[1]*x[2])
+
+    FA = FacetArea(mesh)
+    CV = CellVolume(mesh)
+    ddx = CV/FA
+    ddx_avg = (ddx('+') + ddx('-'))/2.
+    alpha = Constant(4. * degree * (degree + 1.))
+    gamma = Constant(8. * degree * (degree + 1.))
+    penalty_int = alpha / ddx_avg
+    penalty_ext = gamma / ddx
+
+    n = FacetNormal(mesh)
+    aV = (inner(grad(u), grad(v)) + inner(u, v)) * dx  # volume term
+    aIF = (inner(jump(u, n), jump(v, n)) * penalty_int - inner(avg(grad(u)), jump(v, n)) - inner(avg(grad(v)), jump(u, n))) * dIF  # interior facet term
+    aEF = (inner(u, v) * penalty_ext - inner(grad(u), v*n) - inner(grad(v), u*n)) * dEF  # exterior facet term
+    a = aV + aEF + aIF
+    L = f*v*dx
+
+    # Compute solution
+    sol = Function(V)
+
+    solve(a == L, sol, solver_parameters={'pc_type': 'ilu',
+                                          'ksp_type': 'lgmres'})
+
+    # Analytical solution
+    if meshtype in ['planar-quad', 'planar-tri', 'extruded']:
+        f.project(sin(x[0]*pi*2)*sin(x[1]*pi*2))
+    elif meshtype in ['spherical-quad', 'spherical-tri']:
+        f.project(x[0]*x[1]*x[2]/13.0)
+    return sqrt(assemble(dot(sol - f, sol - f) * dx))
+
+
+@pytest.mark.parametrize(('degree', 'meshd', 'meshtype', 'threshold'),
+                         [(1, 1, 'planar-quad', 1.8),
+                          (2, 1, 'planar-quad', 2.9),
+                          (1, 1, 'planar-tri', 1.5),
+                          (2, 1, 'planar-tri', 2.9),
+                          (1, 1, 'extruded', 1.8),
+                          (2, 1, 'extruded', 2.9),
+                          (2, 2, 'spherical-quad', 2.9),
+                          (3, 3, 'spherical-quad', 3.9),
+                          (2, 2, 'spherical-tri', 2.9),
+                          (3, 3, 'spherical-tri', 3.9)])
+def test_firedrake_laplacian_IP_l2pullbacks(degree, meshd, meshtype, threshold):
+    diff = np.array([laplacian_IP(i, degree, meshd, meshtype) for i in range(2, 5)])
+    print("l2 error norms:", diff)
+    conv = np.log2(diff[:-1] / diff[1:])
+    print("convergence order:", conv)
+    assert (np.array(conv) > threshold).all()
diff --git a/tests/regression/test_mse.py b/tests/regression/test_mse.py
new file mode 100644
index 0000000000..2776686094
--- /dev/null
+++ b/tests/regression/test_mse.py
@@ -0,0 +1,291 @@
+import pytest
+import numpy as np
+from firedrake import *
+
+
+def get_complex(meshtype, hdegree, vdegree=None):
+
+    if meshtype == 'extruded':
+        vdegree = vdegree or hdegree
+        h1_horiz = FiniteElement("GLL", "interval", hdegree)
+        l2_horiz = FiniteElement("GLL-Edge L2", "interval", hdegree-1)
+        h1_vert = FiniteElement("GLL", "interval", vdegree)
+        l2_vert = FiniteElement("GLL-Edge L2", "interval", vdegree-1)
+        hcurlelem = HCurl(TensorProductElement(h1_horiz, l2_vert)) + HCurl(TensorProductElement(l2_horiz, h1_vert))
+        hdivelem = HDiv(TensorProductElement(h1_horiz, l2_vert)) + HDiv(TensorProductElement(l2_horiz, h1_vert))
+        l2elem = TensorProductElement(l2_horiz, l2_vert)
+        h1elem = TensorProductElement(h1_horiz, h1_vert)
+    else:
+        l2elem = FiniteElement('DQ L2', quadrilateral, hdegree-1, variant='mse')
+        hdivelem = FiniteElement('RTCF', quadrilateral, hdegree, variant='mse')
+        hcurlelem = FiniteElement('RTCE', quadrilateral, hdegree, variant='mse')
+        h1elem = FiniteElement('Q', quadrilateral, hdegree, variant='mse')
+
+    return h1elem, hcurlelem, hdivelem, l2elem
+
+
+def get_mesh(r, meshtype, meshd=None):
+
+    meshd = meshd or 1
+    if meshtype == 'spherical':
+        mesh = UnitCubedSphereMesh(refinement_level=r, degree=meshd)
+        x = SpatialCoordinate(mesh)
+        mesh.init_cell_orientations(x)
+    elif meshtype == 'planar':
+        mesh = UnitSquareMesh(2**r, 2**r, quadrilateral=True)
+    elif meshtype == 'extruded':
+        basemesh = UnitIntervalMesh(2**r)
+        mesh = ExtrudedMesh(basemesh, 2**r)
+    return mesh
+
+
+def helmholtz_mixed(r, meshtype, hdegree, vdegree=None, meshd=None, useaction=False):
+
+    mesh = get_mesh(r, meshtype, meshd=meshd)
+    _, _, hdivelem, l2elem = get_complex(meshtype, hdegree, vdegree=vdegree)
+
+    V1 = FunctionSpace(mesh, hdivelem, name="V")
+    V2 = FunctionSpace(mesh, l2elem, name="P")
+    W = V1 * V2
+
+    # Define variational problem
+    lmbda = 1
+    u, p = TrialFunctions(W)
+    v, q = TestFunctions(W)
+    f = Function(V2)
+
+    x = SpatialCoordinate(mesh)
+    if meshtype in ['planar', 'extruded']:
+        f.project((1+8*pi*pi)*sin(x[0]*pi*2)*sin(x[1]*pi*2))
+    elif meshtype == 'spherical':
+        f.project(x[0]*x[1]*x[2])
+    a = (p*q - q*div(u) + lmbda*inner(v, u) + div(v)*p) * dx
+    L = f*q*dx
+
+    # Compute solution
+    sol = Function(W)
+
+    if useaction:
+        system = action(a, sol) - L == 0
+    else:
+        system = a == L
+
+    # Block system is:
+    # V Ct
+    # Ch P
+    # Eliminate V by forming a schur complement
+    solve(system, sol, solver_parameters={'pc_type': 'fieldsplit',
+                                          'pc_fieldsplit_type': 'schur',
+                                          'ksp_type': 'cg',
+                                          'pc_fieldsplit_schur_fact_type': 'FULL',
+                                          'fieldsplit_V_ksp_type': 'cg',
+                                          'fieldsplit_P_ksp_type': 'cg'})
+
+    if meshtype in ['planar', 'extruded']:
+        f.project(sin(x[0]*pi*2)*sin(x[1]*pi*2))
+        return sqrt(assemble(dot(sol[2] - f, sol[2] - f) * dx))
+    elif meshtype == 'spherical':
+        _, u = sol.split()
+        f.project(x[0]*x[1]*x[2]/13.0)
+        return errornorm(f, u, degree_rise=0)
+
+
+# helmholtz mixed on plane
+
+@pytest.mark.parametrize(('degree', 'action', 'threshold'),
+                         [(1, False, 1.9),
+                          (2, False, 2.9),
+                          (3, False, 3.9),
+                          (2, True, 2.9)])
+def test_firedrake_helmholtz_mixed_mse(degree, action, threshold):
+    diff = np.array([helmholtz_mixed(r, 'planar', degree, useaction=action) for r in range(3, 6)])
+    print("l2 error norms:", diff)
+    conv = np.log2(diff[:-1] / diff[1:])
+    print("convergence order:", conv)
+    assert (np.array(conv) > threshold).all()
+
+
+# helmholtz mixed on sphere
+
+@pytest.mark.parametrize(('degree', 'md', 'action', 'threshold'),
+                         [(1, 1, False, 1.67),
+                          (2, 1, False, 1.9),
+                          (2, 2, False, 2.9),
+                          (3, 1, False, 1.87),
+                          (3, 2, False, 2.9),
+                          (3, 3, False, 3.9),
+                          (2, 2, True, 2.9)])
+def test_firedrake_helmholtz_mixed_mse_sphere(degree, md, action, threshold):
+    diff = np.array([helmholtz_mixed(r, 'spherical', degree, meshd=md, useaction=action) for r in range(2, 5)])
+    print("l2 error norms:", diff)
+    conv = np.log2(diff[:-1] / diff[1:])
+    print("convergence order:", conv)
+    assert (np.array(conv) > threshold).all()
+
+
+# helmholtz mixed on extruded
+
+@pytest.mark.parametrize(('hdegree', 'vdegree', 'threshold', 'action'),
+                         [(1, 1, 1.9, False),
+                          (2, 1, 1.9, False),
+                          (1, 2, 1.9, False),
+                          (2, 2, 2.9, False),
+                          (1, 1, 1.9, True)])
+def test_firedrake_helmholtz_mixed_mse_extruded(hdegree, vdegree, threshold, action):
+    diff = np.array([helmholtz_mixed(i, 'extruded', hdegree, vdegree=vdegree, useaction=action) for i in range(3, 6)])
+    print("l2 error norms:", diff)
+    conv = np.log2(diff[:-1] / diff[1:])
+    print("convergence order:", conv)
+    assert (np.array(conv) > threshold).all()
+
+
+# facet integrals- interior and exterior
+# tested using an interior penalty formulation for the Helmholtz using DG/DQ L2 elements
+
+
+def laplacian_IP(r, degree, meshd, meshtype):
+
+    dIF = dS
+    dEF = ds
+    if meshtype == 'extruded':
+        dIF = dS_v + dS_h
+        dEF = ds_t + ds_b + ds_v
+
+    mesh = get_mesh(r, meshtype, meshd=meshd)
+    _, _, _, l2elem = get_complex(meshtype, degree+1, vdegree=degree+1)
+
+    V = FunctionSpace(mesh, l2elem)
+
+    # Define variational problem
+    u = TrialFunction(V)
+    v = TestFunction(V)
+    f = Function(V)
+
+    x = SpatialCoordinate(mesh)
+    if meshtype in ['planar', 'extruded']:
+        f.project((1+8*pi*pi)*sin(x[0]*pi*2)*sin(x[1]*pi*2))
+    elif meshtype == 'spherical':
+        f.project(x[0]*x[1]*x[2])
+
+    FA = FacetArea(mesh)
+    CV = CellVolume(mesh)
+    ddx = CV/FA
+    ddx_avg = (ddx('+') + ddx('-'))/2.
+    alpha = Constant(4. * degree * (degree + 1.))
+    gamma = Constant(8. * degree * (degree + 1.))
+    penalty_int = alpha / ddx_avg
+    penalty_ext = gamma / ddx
+
+    n = FacetNormal(mesh)
+    aV = (inner(grad(u), grad(v)) + inner(u, v)) * dx  # volume term
+    aIF = (inner(jump(u, n), jump(v, n)) * penalty_int - inner(avg(grad(u)), jump(v, n)) - inner(avg(grad(v)), jump(u, n))) * dIF  # interior facet term
+    aEF = (inner(u, v) * penalty_ext - inner(grad(u), v*n) - inner(grad(v), u*n)) * dEF  # exterior facet term
+    a = aV + aEF + aIF
+    L = f*v*dx
+
+    # Compute solution
+    sol = Function(V)
+
+    solve(a == L, sol, solver_parameters={'pc_type': 'ilu',
+                                          'ksp_type': 'lgmres'})
+
+    # Analytical solution
+    if meshtype in ['planar', 'extruded']:
+        f.project(sin(x[0]*pi*2)*sin(x[1]*pi*2))
+    elif meshtype == 'spherical':
+        f.project(x[0]*x[1]*x[2]/13.0)
+    return sqrt(assemble(dot(sol - f, sol - f) * dx))
+
+
+@pytest.mark.parametrize(('degree', 'meshd', 'meshtype', 'threshold'),
+                         [(1, 1, 'planar', 1.8),
+                          (2, 1, 'planar', 2.9),
+                          (1, 1, 'extruded', 1.8),
+                          (2, 1, 'extruded', 2.9),
+                          (2, 2, 'spherical', 2.9),
+                          (3, 3, 'spherical', 3.9)])
+def test_firedrake_laplacian_IP_mse(degree, meshd, meshtype, threshold):
+    diff = np.array([laplacian_IP(i, degree, meshd, meshtype) for i in range(2, 5)])
+    print("l2 error norms:", diff)
+    conv = np.log2(diff[:-1] / diff[1:])
+    print("convergence order:", conv)
+    assert (np.array(conv) > threshold).all()
+
+
+def vector_laplace(r, meshtype, hdegree, vdegree=None):
+    vdegree = vdegree or hdegree
+
+    mesh = get_mesh(r, meshtype)
+    h1elem, hcurlelem, _, _ = get_complex(meshtype, hdegree, vdegree=vdegree)
+    h1elem_high, _, _, _ = get_complex(meshtype, hdegree+1, vdegree=vdegree+1)
+
+    # spaces for calculation
+    V0 = FunctionSpace(mesh, h1elem)
+    V1 = FunctionSpace(mesh, hcurlelem)
+    V = V0*V1
+
+    # spaces to store 'analytic' functions
+    W0 = FunctionSpace(mesh, h1elem_high)
+    W1 = VectorFunctionSpace(mesh, h1elem_high)
+
+    # constants
+    k = 1.0
+    l = 2.0
+
+    xs = SpatialCoordinate(mesh)
+    f_expr = as_vector([pi*pi*(k*k + l*l)*sin(k*pi*xs[0])*cos(l*pi*xs[1]), pi*pi*(k*k + l*l)*cos(k*pi*xs[0])*sin(l*pi*xs[1])])
+    exact_s_expr = -(k+l)*pi*cos(k*pi*xs[0])*cos(l*pi*xs[1])
+    exact_u_expr = as_vector([sin(k*pi*xs[0])*cos(l*pi*xs[1]), cos(k*pi*xs[0])*sin(l*pi*xs[1])])
+
+    f = Function(W1).project(f_expr)
+    exact_s = Function(W0).project(exact_s_expr)
+    exact_u = Function(W1).project(exact_u_expr)
+
+    sigma, u = TrialFunctions(V)
+    tau, v = TestFunctions(V)
+    a = (sigma*tau - dot(u, grad(tau)) + dot(grad(sigma), v) + dot(curl(u), curl(v)))*dx
+    L = dot(f, v)*dx
+
+    out = Function(V)
+
+    # preconditioner for H1 x H(curl)
+    aP = (dot(grad(sigma), grad(tau)) + sigma*tau + dot(curl(u), curl(v)) + dot(u, v))*dx
+
+    solve(a == L, out, Jp=aP,
+          solver_parameters={'pc_type': 'fieldsplit',
+                             'pc_fieldsplit_type': 'additive',
+                             'fieldsplit_0_pc_type': 'lu',
+                             'fieldsplit_1_pc_type': 'lu',
+                             'ksp_monitor': None})
+
+    out_s, out_u = out.split()
+
+    return (sqrt(assemble(dot(out_u - exact_u, out_u - exact_u)*dx)),
+            sqrt(assemble((out_s - exact_s)*(out_s - exact_s)*dx)))
+
+
+@pytest.mark.parametrize(('degree', 'threshold'),
+                         [(1, 0.9),
+                          (2, 1.9),
+                          (3, 2.9),
+                          (2, 0.9)])
+def test_firedrake_helmholtz_vector_mse(degree, threshold):
+    diff = np.array([vector_laplace(r, 'planar', degree) for r in range(2, 5)])
+    print("l2 error norms:", diff)
+    conv = np.log2(diff[:-1] / diff[1:])
+    print("convergence order:", conv)
+    assert (np.array(conv) > threshold).all()
+
+
+@pytest.mark.parametrize(('hdegree', 'vdegree', 'threshold'),
+                         [(1, 1, 0.9),
+                          (2, 1, 0.9),
+                          (1, 2, 0.9),
+                          (2, 2, 1.9),
+                          (1, 1, 0.9)])
+def test_firedrake_helmholtz_vector_mse_extruded(hdegree, vdegree, threshold):
+    diff = np.array([vector_laplace(i, 'extruded', hdegree, vdegree=vdegree) for i in range(2, 5)])
+    print("l2 error norms:", diff)
+    conv = np.log2(diff[:-1, :] / diff[1:, :])
+    print("convergence order:", conv)
+    assert (np.array(conv) > threshold).all()