Reject elements with mismatched sub-entity DOF counts

[garth-wells]

Fixes #1067.

The problem

create_custom_element can be given an element whose sub-entities of the same type carry different numbers of DOFs, for example an anisotropic tensor-product element on a hexahedron. Such an element cannot be represented: the DOF transformations are computed from a single representative sub-entity of each type (compute_entity_transformations picks it with find_first_subentity) and then applied to every sub-entity of that type, and a quadrilateral face holding an a × b grid of DOFs with a != b is not invariant under the face rotation, so its transformation matrix is rank deficient.

We accepted these elements and then did one of three things, depending on rounding:

1. threw an opaque LU decomposition failed: N from precompute::prepare_matrix,
2. read out of bounds while building the sub-entity closure permutations, or
3. silently produced transformations that do not satisfy the group relations of the sub-entity symmetry.

The third is the worst: the element is built without any error and is wrong. The issue reporter had several "passing" axis orderings that fall into this category.

A minimal demonstration, with point-evaluation Lagrange only (no Hermite, no discontinuous factors), building Q(dx, dy, dz) on a hexahedron as a custom element:

element	behaviour before this PR
Q(2,2,2), Q(3,3,3)	correct
Q(3,2,1)	LU decomposition failed: 2
Q(3,1,1)	out-of-bounds read (aborts under libc++ hardening)
Q(3,2,2)	builds; face rotation block is [[0.5, 0.5], [0.5, 0.5]], 20 of 24 base transformations violate T^order == I
Q(4,2,1)	builds; face rotation block has an exact zero singular value, an edge reflection block is [[0]] instead of [[1]]

The fix

Check up front that all sub-entities of a given type carry the same number of DOFs, and throw a descriptive error if not. This catches every case above, including the ones that previously succeeded silently.

The guards around the sub-entity closure permutations test one edge of the sub-entity and then index all of them, which is what reads out of bounds. The new check makes that indexing valid, so this adds a comment saying so rather than repeating the test per edge.

Also fixed

While adding the check, it rejected Hermite on a tetrahedron. That turned out to be a separate pre-existing bug: e-hermite.cpp sized x[3] and M[3] by the number of faces rather than the number of cells, so a tetrahedron reported four cell entities instead of one. Every other element uses cell::num_sub_entities(celltype, 3); Hermite now does too. The extra entries were empty and harmless for the existing code paths, but it is not a valid interpolation layout.

Testing

• New regression test in test/test_custom_element.py.
• Full Python suite: 10126 passed, 337 skipped, 23 xfailed.
• C++ unit tests: 22/22 passed.
• Python demos: 7 passed.

🤖 Generated with Claude Code