From e3664a1cad14cf35fa864d0d27f069c1c09f7bca Mon Sep 17 00:00:00 2001 From: Matthew Scroggs Date: Thu, 18 Jun 2026 10:56:26 +0100 Subject: [PATCH 1/2] update some tests --- test/test_jit_expression.py | 38 +++++++++++++++++++------------------ test/test_jit_forms.py | 9 ++++----- 2 files changed, 24 insertions(+), 23 deletions(-) diff --git a/test/test_jit_expression.py b/test/test_jit_expression.py index e92b129c4..f10fc5d98 100644 --- a/test/test_jit_expression.py +++ b/test/test_jit_expression.py @@ -299,14 +299,11 @@ def test_facet_expression(compile_args): consts = np.array([], dtype=dtype) entity_index = np.array([0], dtype=np.intc) quad_perm = np.array([0], dtype=np.dtype("uint8")) - tangents = np.array([coords[1] - coords[2], coords[2] - coords[0], coords[0] - coords[1]]) - midpoints = np.array( - [ - coords[1] + (coords[2] - coords[1]) / 2, - coords[0] + (coords[2] - coords[0]) / 2, - coords[1] + (coords[1] - coords[0]) / 2, - ] - ) + + edges = basix.topology(basix.CellType.triangle)[1] + tangents = np.array([coords[j] - coords[i] for i, j in edges]) + midpoints = np.array([(coords[i] + coords[j]) / 2 for i, j in edges]) + for i, (tangent, midpoint) in enumerate(zip(tangents, midpoints)): # normalize tangent tangent /= np.linalg.norm(tangent) @@ -329,7 +326,7 @@ def test_facet_expression(compile_args): assert np.isclose(np.linalg.norm(output), 1) # Check that facet normal is pointing out of the cell - assert np.dot(midpoint - coords[i], output) > 0 + assert max(np.dot(midpoint - c, output) for c in coords) > 0 def test_facet_geometry_expressions(compile_args): @@ -518,11 +515,12 @@ def test_mixed_mesh_expression(compile_args): # Define exact normals face_normal = np.array([0.0, 0.0, 1.0]) - e0 = coords[2] - coords[1] - e1 = coords[0] - coords[2] - e2 = coords[1] - coords[0] - edges = np.array([e0, e1, e2]) + topology = basix.cell.topology(basix.CellType.triangle) + edges = np.array([coords[j] - coords[i] for i, j in topology[1]]) edge_normals = np.cross(edges, face_normal)[:, :2] + for i, orient in enumerate(basix.cell.facet_orientations(basix.CellType.triangle)): + if orient == 1: + edge_normals[i] *= -1 norms = np.linalg.norm(edge_normals, axis=1, keepdims=True) edge_normals_normalized = edge_normals / norms @@ -554,7 +552,6 @@ def test_mixed_mesh_expression(compile_args): @pytest.mark.parametrize("facet_perm", [0, 1], ids=["no_perm", "perm"]) @pytest.mark.parametrize("local_index", [0, 1, 2], ids=["facet0", "facet1", "facet2"]) def test_expression_facet_perm(compile_args, facet_perm, local_index): - c_el = basix.ufl.element("Lagrange", "triangle", 1, shape=(2,)) mesh = ufl.Mesh(c_el) expr = ufl.SpatialCoordinate(mesh) @@ -586,10 +583,10 @@ def test_expression_facet_perm(compile_args, facet_perm, local_index): entity_index = np.array([local_index], dtype=np.intc) # Perm 0, means that global facet is ordered as - # 1----2 + # 0----1 # and quadrature point is not reflected # Perm 1, means that global facet is ordered as - # 2----1 + # 1----0 # and quadrature point should be reflected quad_perm = np.array([facet_perm], dtype=np.uint8) @@ -606,8 +603,13 @@ def test_expression_facet_perm(compile_args, facet_perm, local_index): ffi.NULL, ) + edges = basix.topology(basix.CellType.triangle)[1] + ordered_points = points * (facet_perm == 0) + (1 - points) * (facet_perm == 1) - facet = np.delete(coords, local_index, axis=0)[:, :2] + facet = [coords[i][:2] for i in edges[local_index]] edge = facet[1] - facet[0] + print(local_index, facet) exact_value = facet[0] + ordered_points * edge - np.testing.assert_allclose(output, exact_value) + print(output) + print(exact_value) + np.testing.assert_allclose(output, exact_value, 1e-10) diff --git a/test/test_jit_forms.py b/test/test_jit_forms.py index 5274977f4..1a5395fc3 100644 --- a/test/test_jit_forms.py +++ b/test/test_jit_forms.py @@ -618,7 +618,7 @@ def lagrange_triangle_symbolic(order, corners=((1, 0), (2, 0), (0, 1)), fun=lamb # vertices eval_points = [S(c) for c in corners] # edges - for e in [(1, 2), (0, 2), (0, 1)]: + for e in basix.topology(basix.CellType.triangle)[1]: p0 = corners[e[0]] p1 = corners[e[1]] if order > 3: @@ -722,7 +722,7 @@ def lagrange_tetrahedron_symbolic( # vertices eval_points = [S(c) for c in corners] # edges - for e in [(2, 3), (1, 3), (1, 2), (0, 3), (0, 2), (0, 1)]: + for e in basix.topology(basix.CellType.tetrahedron)[1]: p0 = corners[e[0]] p1 = corners[e[1]] if order > 3: @@ -738,7 +738,7 @@ def lagrange_tetrahedron_symbolic( for i in range(1, order) ] # face - for f in [(1, 2, 3), (0, 2, 3), (0, 1, 3), (0, 1, 2)]: + for f in basix.topology(basix.CellType.tetrahedron)[2]: p0 = corners[f[0]] p1 = corners[f[1]] p2 = corners[f[2]] @@ -1173,8 +1173,7 @@ def test_derivative_domains(compile_args): @pytest.mark.parametrize("dtype", ["float64"]) -@pytest.mark.parametrize("permutation", [[0], [1]]) -def test_mixed_dim_form(compile_args, dtype, permutation): +def test_mixed_dim_form(compile_args, dtype, permutation, local_index): """Test that the local element tensor corresponding to a mixed-dimensional form is correct. The form involves an integral over a facet of the cell. The trial function and a coefficient f are of codim 0. The test function and a coefficient g are of codim 1. We compare against another From 9b8aa5b96a91afb38cbedc8c26c267bf429d8f58 Mon Sep 17 00:00:00 2001 From: Matthew Scroggs Date: Thu, 18 Jun 2026 13:19:18 +0100 Subject: [PATCH 2/2] make test robust to cell renumbering --- test/test_jit_expression.py | 3 --- test/test_jit_forms.py | 31 ++++++++++++++++++++----------- 2 files changed, 20 insertions(+), 14 deletions(-) diff --git a/test/test_jit_expression.py b/test/test_jit_expression.py index f10fc5d98..61213f9e8 100644 --- a/test/test_jit_expression.py +++ b/test/test_jit_expression.py @@ -608,8 +608,5 @@ def test_expression_facet_perm(compile_args, facet_perm, local_index): ordered_points = points * (facet_perm == 0) + (1 - points) * (facet_perm == 1) facet = [coords[i][:2] for i in edges[local_index]] edge = facet[1] - facet[0] - print(local_index, facet) exact_value = facet[0] + ordered_points * edge - print(output) - print(exact_value) np.testing.assert_allclose(output, exact_value, 1e-10) diff --git a/test/test_jit_forms.py b/test/test_jit_forms.py index 1a5395fc3..5d0c97115 100644 --- a/test/test_jit_forms.py +++ b/test/test_jit_forms.py @@ -1173,13 +1173,18 @@ def test_derivative_domains(compile_args): @pytest.mark.parametrize("dtype", ["float64"]) -def test_mixed_dim_form(compile_args, dtype, permutation, local_index): +@pytest.mark.parametrize("permutation", [[0], [1]]) +def test_mixed_dim_form(compile_args, dtype, permutation): """Test that the local element tensor corresponding to a mixed-dimensional form is correct. The form involves an integral over a facet of the cell. The trial function and a coefficient f are of codim 0. The test function and a coefficient g are of codim 1. We compare against another form where the test function and g are codim 0 but have the same trace on the facet. """ + (facet_index,) = [ + i for i, v in enumerate(basix.topology(basix.CellType.triangle)[1]) if v == [1, 2] + ] + def tabulate_tensor(ele_type, V_cell_type, W_cell_type, coeffs): "Helper function to create a form and compute the local element tensor" V_ele = basix.ufl.element(ele_type, V_cell_type, 2) @@ -1212,7 +1217,7 @@ def tabulate_tensor(ele_type, V_cell_type, W_cell_type, coeffs): A = np.zeros((W_ele.dim, V_ele.dim), dtype=dtype) w = np.array(coeffs, dtype=dtype) c = np.array([], dtype=dtype) - facet = np.array([0], dtype=np.intc) + facet = np.array([facet_index], dtype=np.intc) perm = np.array(permutation, dtype=np.uint8) xdtype = dtype_to_scalar_dtype(dtype) @@ -1261,32 +1266,33 @@ def tabulate_tensor(ele_type, V_cell_type, W_cell_type, coeffs): A_ref = A_ref[1:][:] # The orientation of the interval is assumed to be fixed (since it is given uniquely - # by its global orientation in the mesh). Let us focus on local facet 0. It can have - # two possible orientations depending on the cell topology: + # by its global orientation in the mesh). Let us focus on the local facet between vertices + # 1 and 2. It can have two possible orientations depending on the cell topology: # # 2 1 1 1 # | \ \ | \ \ # | \ \ or | \ \ # 0---1 0 0---2 0 # - # In the second case, local facet 0 is flipped relative to the interval. If we look at + # In the second case, the local facet is flipped relative to the interval. If we look at # the degrees of freedom second-order Lagrange on the triangle and first-order Lagrange # on the interval, we have # # 2 1 1 1 # | \ \ | \ \ - # 4 3 \ 5 3 \ + # B C \ A C \ # | \ \ or | \ \ - # 0--5--1 0 0--4--2 0 + # 0--A--1 0 0--B--2 0 # # Since the trial function is defined on the triangle, the second case (where the - # permutation is 1) is thus the same as swapping cols 1 and 2 and cols 4 and 5 of the + # permutation is 1) is thus the same as swapping cols 1 and 2 and cols A and B of the # first case result. # NOTE: Although we choose to fix the orientation of the interval, the same result # can be obtained by fixing the triangle and considering flipping the interval. if permutation[0] == 1: A_ref[:, [1, 2]] = A_ref[:, [2, 1]] - A_ref[:, [4, 5]] = A_ref[:, [5, 4]] + edge_dofs = [3 + i for i in range(3) if i != facet_index] + A_ref[:, edge_dofs] = A_ref[:, edge_dofs[::-1]] assert np.allclose(A, A_ref) @@ -1660,7 +1666,7 @@ def tabulate_tensor(ele_type, V_cell_type, W_cell_type, coeffs, l_i): A_ref = tabulate_tensor(ele_type, V_cell_type, V_cell_type, coeff_data, local_entity_index) # Map from local edge index to (local) DOF indices on that edge - local_index_to_slice = {0: [2, 3], 1: [1, 3], 2: [1, 2], 3: [0, 3], 4: [0, 2], 5: [0, 1]} + local_index_to_slice = basix.ufl.element(ele_type, V_cell_type, 1).entity_closure_dofs[1] A_ref = A_ref[local_index_to_slice[local_entity_index]] @@ -1832,7 +1838,10 @@ def tabulate_tensor(ele_type, V_cell_type, l_i): V_cell_type = "tetrahedron" # Tabulate the tensor for the mixed-dimensional form - b = tabulate_tensor(ele_type, V_cell_type, 0) + (ridge_index,) = [ + i for i, v in enumerate(basix.topology(basix.CellType.tetrahedron)[1]) if v == [2, 3] + ] + b = tabulate_tensor(ele_type, V_cell_type, ridge_index) # Ref first ridge integral length rl = np.sqrt(y_ext**2 + z_ext**2) / y_ext