diff --git a/tests/submesh/test_submesh_solve.py b/tests/submesh/test_submesh_solve.py
index 08fd1f6e66..f13944414a 100644
--- a/tests/submesh/test_submesh_solve.py
+++ b/tests/submesh/test_submesh_solve.py
@@ -421,3 +421,44 @@ def test_submesh_solve_mixed_poisson_check_sanity_3d(hexahedral, degree, submesh
     sigma_error, u_error = _mixed_poisson_solve_3d(hexahedral, degree, submesh_region)
     assert sigma_error < 0.07
     assert u_error < 0.003
+
+
+@pytest.mark.parallel(nprocs=4)
+@pytest.mark.parametrize('simplex', [True, False])
+@pytest.mark.parametrize('nref', [1, 3])
+@pytest.mark.parametrize('degree', [2, 4])
+def test_submesh_solve_cell_cell_equation_bc(nref, degree, simplex):
+    mesh = RectangleMesh(3 ** nref, 2 ** nref, 3., 2., quadrilateral=not simplex)
+    x, y = SpatialCoordinate(mesh)
+    label_outer = 101
+    label_inner = 100
+    label_interface = 5  # automatically labeled by Submesh
+    DG0 = FunctionSpace(mesh, "DG", 0)
+    f_outer = Function(DG0).interpolate(conditional(Or(Or(x < 1., x > 2.), y > 1.), 1, 0))
+    f_inner = Function(DG0).interpolate(conditional(And(And(x > 1., x < 2.), y < 1.), 1, 0))
+    mesh = RelabeledMesh(mesh, [f_outer, f_inner], [label_outer, label_inner])
+    x, y = SpatialCoordinate(mesh)
+    mesh_outer = Submesh(mesh, dmcommon.CELL_SETS_LABEL, label_outer, mesh.topological_dimension())
+    x_outer, y_outer = SpatialCoordinate(mesh_outer)
+    mesh_inner = Submesh(mesh, dmcommon.CELL_SETS_LABEL, label_inner, mesh.topological_dimension())
+    x_inner, y_inner = SpatialCoordinate(mesh_inner)
+    V_outer = FunctionSpace(mesh_outer, "CG", degree)
+    V_inner = FunctionSpace(mesh_inner, "CG", degree)
+    V = V_outer * V_inner
+    u = TrialFunction(V)
+    v = TestFunction(V)
+    sol = Function(V)
+    u_outer, u_inner = split(u)
+    v_outer, v_inner = split(v)
+    dx_outer = Measure("dx", domain=mesh_outer)
+    dx_inner = Measure("dx", domain=mesh_inner)
+    ds_outer = Measure("ds", domain=mesh_outer)
+    ds_inner = Measure("ds", domain=mesh_inner)
+    a = inner(grad(u_outer), grad(v_outer)) * dx_outer + \
+        inner(u_inner, v_inner) * dx_inner
+    L = inner(x * y, v_inner) * dx_inner
+    dbc = DirichletBC(V.sub(0), x_outer * y_outer, (1, 2, 3, 4))
+    ebc = EquationBC(inner(u_outer - u_inner, v_outer) * ds_outer(label_interface) == inner(Constant(0.), v_outer) * ds_outer(label_interface), sol, label_interface, V=V.sub(0))
+    solve(a == L, sol, bcs=[dbc, ebc])
+    assert sqrt(assemble(inner(sol[0] - x * y, sol[0] - x * y) * dx_outer)) < 1.e-13
+    assert sqrt(assemble(inner(sol[1] - x * y, sol[1] - x * y) * dx_inner)) < 1.e-13