From 25f835dc690238852d5c348834a59bec2aec70fd Mon Sep 17 00:00:00 2001 From: Emma Rothwell Date: Wed, 24 Jan 2024 17:14:42 +0000 Subject: [PATCH] added mesh coord change test + refactored --- .../test_restricted_function_space.py | 76 ++++++++----------- 1 file changed, 33 insertions(+), 43 deletions(-) diff --git a/tests/regression/test_restricted_function_space.py b/tests/regression/test_restricted_function_space.py index 01cbb57a10..9787e0d2a6 100644 --- a/tests/regression/test_restricted_function_space.py +++ b/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])