diff --git a/tests/regression/test_restricted_function_space.py b/tests/regression/test_restricted_function_space.py
index 869380b4e0..c9de153eb0 100644
--- a/tests/regression/test_restricted_function_space.py
+++ b/tests/regression/test_restricted_function_space.py
@@ -3,7 +3,7 @@
 from firedrake import *
 
 
-def compare_function_space_assembly(mesh, function_space, restricted_function_space, bcs, res_bcs):
+def compare_function_space_assembly(mesh, function_space, restricted_function_space, bcs, res_bcs=[]):
     x, y = SpatialCoordinate(mesh)
     u = TrialFunction(function_space)
     v = TestFunction(function_space)
@@ -14,7 +14,7 @@ def compare_function_space_assembly(mesh, function_space, restricted_function_sp
     restricted_form = u2 * v2 * dx
 
     normal_fs_matrix = assemble(original_form, bcs=bcs)
-    restricted_fs_matrix = assemble(restricted_form)
+    restricted_fs_matrix = assemble(restricted_form, bcs=res_bcs)
 
     identity = np.identity(normal_fs_matrix.M.nrows)
     delete_rows = []
@@ -32,8 +32,6 @@ def compare_function_space_assembly(mesh, function_space, restricted_function_sp
 
     assert (restricted_fs_matrix.M.nrows == np.shape(normal_fs_matrix_reduced)[0])
     assert (restricted_fs_matrix.M.ncols == np.shape(normal_fs_matrix_reduced)[1])
-    print(normal_fs_matrix_reduced)
-    print(restricted_fs_matrix.M.values)
     assert (np.array_equal(normal_fs_matrix_reduced, restricted_fs_matrix.M.values))
 
 
@@ -44,7 +42,6 @@ def test_restricted_function_space_1_1_square(j):
     bc = DirichletBC(V, 0, 2)
     V_res = RestrictedFunctionSpace(V, name="Restricted", boundary_set=[2])
     res_bc = DirichletBC(V_res, 0, 2)
-
     compare_function_space_assembly(mesh, V, V_res, [bc], [res_bc])
 
 
@@ -53,13 +50,13 @@ def test_restricted_function_space_j_j_square(j):
     mesh = UnitSquareMesh(j, j)
     V = FunctionSpace(mesh, "CG", 1)
     bc = DirichletBC(V, 0, 2)
-    V_res = RestrictedFunctionSpace(V, name="Restricted", bcs=[bc])
+    V_res = RestrictedFunctionSpace(V, name="Restricted", boundary_set=[2])
 
     compare_function_space_assembly(mesh, V, V_res, [bc])
 
 
 def test_poisson_homogeneous_bcs():
-    mesh = UnitSquareMesh(3, 3)
+    mesh = UnitSquareMesh(1, 1)
     V = FunctionSpace(mesh, "CG", 2)
     f = Function(V)
     x, y = SpatialCoordinate(mesh)
@@ -73,7 +70,7 @@ def test_poisson_homogeneous_bcs():
     original_form = inner(grad(u), grad(v)) * dx
 
     bc = DirichletBC(V, 0, 1)
-    V_res = RestrictedFunctionSpace(V, name="Restricted", bcs=[bc])
+    V_res = RestrictedFunctionSpace(V, name="Restricted", boundary_set=[1])
     f2 = Function(V_res)
 
     f2.interpolate(16 * pi**2 * (y-1)**2 * y**2 - 2 * (y-1)**2 - 8 * (y-1)*y
@@ -87,11 +84,12 @@ def test_poisson_homogeneous_bcs():
     L_res = inner(f2, v2) * dx
 
     u = Function(V)
+    bc2 = DirichletBC(V_res, Constant(0), 1)
     u2 = Function(V_res)
+
     solve(original_form == L, u, bcs=[bc])
-    solve(restricted_form == L_res, u2)
+    solve(restricted_form == L_res, u2, bcs=[bc2])
 
-    # might run into problems if other non-boundary nodes evaluate at 0?
     u_data_remove_zeros = u.dat.data[u.dat.data != 0]  # correspond to boundary
     u2_data_remove_zeros = u2.dat.data[u2.dat.data != 0]
 
@@ -99,22 +97,56 @@ def test_poisson_homogeneous_bcs():
 
 
 def test_poisson_inhomogeneous_bcs():
-    mesh = UnitSquareMesh(3, 3)
-    V = FunctionSpace(mesh, "CG", 2) 
+    mesh = UnitSquareMesh(1, 2)
+    V = FunctionSpace(mesh, "CG", 2)
     x, y = SpatialCoordinate(mesh)
-
-    bc = DirichletBC(V, 0, 1)
-    bc2 = DirichletBC(V, 1, 2)
-    V_res = RestrictedFunctionSpace(V, name="Restricted", bcs=[bc, bc2])
-
+    V_res = RestrictedFunctionSpace(V, boundary_set=[1, 2])
+    bc = DirichletBC(V_res, 0, 1)
+    bc2 = DirichletBC(V_res, 1, 2)
     u2 = TrialFunction(V_res)
     v2 = TestFunction(V_res)
     restricted_form = inner(grad(u2), grad(v2)) * dx
     u = Function(V_res)
     rhs = Constant(0) * v2 * dx
+
     solve(restricted_form == rhs, u, bcs=[bc, bc2])
-    return u
-    #assert errornorm(SpatialCoordinate(mesh)[0], u) < 1.e-12
+
+    assert errornorm(SpatialCoordinate(mesh)[0], u) < 1.e-12
+
+
+@pytest.mark.parametrize("j", ["2", "sin(x) * y", "x**2"])
+def test_poisson_inhomogeneous_bcs_2(j):
+    mesh = UnitSquareMesh(1, 1)
+    V = FunctionSpace(mesh, "CG", 2)
+
+    x, y = SpatialCoordinate(mesh)
+
+    bc3 = DirichletBC(V, 0, 1)
+    bc4 = DirichletBC(V, Function(V).interpolate(eval(j)), 2)
+    u = TrialFunction(V)
+    v = TestFunction(V)
+
+    V_res = RestrictedFunctionSpace(V, name="Restricted", boundary_set=[1, 2])
+
+    bc = DirichletBC(V_res, 0, 1)
+    bc2 = DirichletBC(V_res, Function(V_res).interpolate(eval(j)), 2)
+
+    u2 = TrialFunction(V_res)
+    v2 = TestFunction(V_res)
+
+    original_form = inner(grad(u), grad(v)) * dx
+    restricted_form = inner(grad(u2), grad(v2)) * dx
+
+    u = Function(V)
+    u2 = Function(V_res)
+
+    rhs = Constant(0) * v * dx
+    rhs2 = Constant(0) * v2 * dx
+
+    solve(original_form == rhs, u, bcs=[bc3, bc4])
+    solve(restricted_form == rhs2, u2, bcs=[bc, bc2])
+
+    assert np.linalg.norm(np.sort(u2.dat.data) - np.sort(u.dat.data)) < 1.e-12
 
 
 @pytest.mark.parametrize("j", [1, 2, 5])
@@ -125,6 +157,6 @@ def test_restricted_function_space_coord_change(j):
     new_mesh = Mesh(Function(V).interpolate(as_vector([x, y])))
     new_V = FunctionSpace(new_mesh, "CG", j)
     bc = DirichletBC(new_V, 0, 1)
-    new_V_restricted = RestrictedFunctionSpace(new_V, name="Restricted", bcs=[bc])
+    new_V_restricted = RestrictedFunctionSpace(new_V, name="Restricted", boundary_set=[1])
 
     compare_function_space_assembly(new_mesh, new_V, new_V_restricted, [bc])