added mesh coord change test + refactored

tests/regression/test_restricted_function_space.py

@@ -3,23 +3,17 @@
 from firedrake import *
 
 
-@pytest.mark.parametrize("j", [1, 2, 5])
-def test_restricted_function_space_1_1_square(j):
-    mesh = UnitSquareMesh(1, 1)
-    V = FunctionSpace(mesh, "CG", j)
+def compare_function_space_assembly(mesh, function_space, restricted_function_space, bcs):
     x, y = SpatialCoordinate(mesh)
-    u = TrialFunction(V)
-    v = TestFunction(V)
+    u = TrialFunction(function_space)
+    v = TestFunction(function_space)
     original_form = u * v * dx
 
-    bc = DirichletBC(V, 0, 2)
-    V_res = RestrictedFunctionSpace(V, name="Restricted", bcs=[bc])
-
-    u2 = TrialFunction(V_res)
-    v2 = TestFunction(V_res)
+    u2 = TrialFunction(restricted_function_space)
+    v2 = TestFunction(restricted_function_space)
     restricted_form = u2 * v2 * dx
-    restricted_fs_matrix = assemble(restricted_form) 
-    normal_fs_matrix = assemble(original_form, bcs=[bc])
+    restricted_fs_matrix = assemble(restricted_form)
+    normal_fs_matrix = assemble(original_form, bcs=bcs)
 
     identity = np.identity(normal_fs_matrix.M.nrows)
     delete_rows = []
@@ -41,44 +35,27 @@ def test_restricted_function_space_1_1_square(j):
 
 
 @pytest.mark.parametrize("j", [1, 2, 5])
-def test_restricted_function_space_j_j_square(j):
-    mesh = UnitSquareMesh(j, j)
-    V = FunctionSpace(mesh, "CG", 3)  # this fails for CG1
-    x, y = SpatialCoordinate(mesh)
-    u = TrialFunction(V)
-    v = TestFunction(V)
-    original_form = u * v * dx
+def test_restricted_function_space_1_1_square(j):
+    mesh = UnitSquareMesh(1, 1)
+    V = FunctionSpace(mesh, "CG", j)
 
     bc = DirichletBC(V, 0, 2)
     V_res = RestrictedFunctionSpace(V, name="Restricted", bcs=[bc])
 
-    u2 = TrialFunction(V_res)
-    v2 = TestFunction(V_res)
-    restricted_form = u2 * v2 * dx
-    restricted_fs_matrix = assemble(restricted_form) 
-    normal_fs_matrix = assemble(original_form, bcs=[bc])
+    compare_function_space_assembly(mesh, V, V_res, [bc])
 
-    identity = np.identity(normal_fs_matrix.M.nrows)
-    delete_rows = []
-    for i in range(normal_fs_matrix.M.nrows):
-        row = normal_fs_matrix.M.values[i, :]
-        comparison = row == identity[i, :]
-        if comparison.all():
-            delete_rows.append(i)
 
-    # Remove all rows/columns associated with boundaries (i.e. identity)
-    normal_fs_matrix_reduced = np.delete(normal_fs_matrix.M.values, delete_rows,
-                                         axis=0)
-    normal_fs_matrix_reduced = np.delete(normal_fs_matrix_reduced, delete_rows,
-                                         axis=1)
+@pytest.mark.parametrize("j", [1, 2, 5])
+def test_restricted_function_space_j_j_square(j):
+    mesh = UnitSquareMesh(j, j)
+    V = FunctionSpace(mesh, "CG", 2)  # this fails for CG1
+    bc = DirichletBC(V, 0, 2)
+    V_res = RestrictedFunctionSpace(V, name="Restricted", bcs=[bc])
 
-    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])
-    assert (np.array_equal(normal_fs_matrix_reduced, restricted_fs_matrix.M.values))
+    compare_function_space_assembly(mesh, V, V_res, [bc])
 
 
 def test_poisson_homogeneous_bcs():
-
     mesh = UnitSquareMesh(3, 3)
     V = FunctionSpace(mesh, "CG", 2)   # fails for CG1
     f = Function(V)
@@ -87,7 +64,7 @@ def test_poisson_homogeneous_bcs():
     u = TrialFunction(V)
     v = TestFunction(V)
 
-    f.interpolate(16 * pi**2 * (y-1)**2 * y**2 - 2 * (y-1)**2 - 8 * (y-1)*y 
+    f.interpolate(16 * pi**2 * (y-1)**2 * y**2 - 2 * (y-1)**2 - 8 * (y-1)*y
                   - 2*y**2)*sin(4*pi*x)
 
     original_form = inner(grad(u), grad(v)) * dx
@@ -112,7 +89,20 @@ def test_poisson_homogeneous_bcs():
     solve(restricted_form == L_res, u2)
 
     # 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 -> different ordering for RFS / FS
+    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]
 
     assert (np.all(np.isclose(u_data_remove_zeros, u2_data_remove_zeros, 1e-16)))
+
+
+@pytest.mark.parametrize("j", [1, 2, 5])
+def test_restricted_function_space_coord_change(j):
+    mesh = UnitSquareMesh(1, 2)
+    x, y = SpatialCoordinate(mesh)
+    V = VectorFunctionSpace(mesh, "CG", 2)
+    new_mesh = Mesh(Function(V).interpolate(as_vector([x, y])))
+    new_V = FunctionSpace(new_mesh, "CG", j)  # fails for 2
+    bc = DirichletBC(new_V, 0, 1)
+    new_V_restricted = RestrictedFunctionSpace(new_V, name="Restricted", bcs=[bc])
+
+    compare_function_space_assembly(new_mesh, new_V, new_V_restricted, [bc])