Skip to content

Make tests robust to reference cell renumbering #843

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
mscroggs wants to merge 2 commits into
base: main
Choose a base branch
from
Open

Make tests robust to reference cell renumbering #843

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
e3664a1
update some tests
mscroggs Jun 18, 2026
9b8aa5b
make test robust to cell renumbering
mscroggs Jun 18, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
35 changes: 17 additions & 18 deletions test/test_jit_expression.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand All @@ -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):
Expand Down Expand Up @@ -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

Expand Down Expand Up @@ -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)
Expand Down Expand Up @@ -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)

Expand All @@ -606,8 +603,10 @@ 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]
exact_value = facet[0] + ordered_points * edge
np.testing.assert_allclose(output, exact_value)
np.testing.assert_allclose(output, exact_value, 1e-10)
34 changes: 21 additions & 13 deletions test/test_jit_forms.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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:
Expand Down Expand Up @@ -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:
Expand All @@ -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]]
Expand Down Expand Up @@ -1181,6 +1181,10 @@ def test_mixed_dim_form(compile_args, dtype, permutation):
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)
Expand Down Expand Up @@ -1213,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)
Expand Down Expand Up @@ -1262,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)

Expand Down Expand Up @@ -1661,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]]

Expand Down Expand Up @@ -1833,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
Expand Down
Loading