Fix coordinate element access of codim >=1 assembly

[jorgensd]

Proper fix for what was attempted in: FEniCS/ufl#513.

Summary

Co-dimensional spatial quantities such as spatialcoordinate and jacobians is required for pull backs or gradients when used in the variational formulation. For instance ufl.grad(u_submesh_codim1) * ds(domain=parent_mesh) requires the coordinate dofs of both the parent and submesh, as the gradient becomes the tangential gradient. The coordinate dofs are already packed with the parent cell.

However, the submesh is ordered with a consistent global dof ordering, which doesn't align with the sub-entity-closure of the parent coordinate dofs.
This is resolved by adding new tables internally in FFCx that uses the quadrature permutations to correct the alignment for arbitrary order coordinate elements.

This PR has been generated with the assistance of CLAUDE (Sonnet 5). I've reviewed all the contents, modified them as I seem fit and take responsibility for the modifications.

[garth-wells]

I checked out the branch, ran ruff/mypy and the new tests (all pass), diffed the generated C against main for standard codim-0 forms (identical modulo set-iteration ordering in comments and table names), and compiled a few adversarial forms to probe the new code path.

The core idea looks right to me: gather the submesh's coordinate dofs out of the parent's coordinate_dofs through a per-entity, per-orientation closure-dofs table, and use permute_subentity_closure_inv so it also works for arbitrary-degree coordinate elements. The table tests are genuinely good — deriving independent ground truth from permute_quadrature_* rather than round-tripping basix against itself is what makes them worth having.

Three things I'd like addressed before merge.

1. needs_facet_permutations is not set when the submesh enters only through geometry

is_mixed_dim in ffcx/ir/integral.py is derived only from extract_arguments + extract_coefficients:

coeffs_and_arguments = ufl.algorithms.analysis.extract_arguments(
    expression
) + ufl.algorithms.analysis.extract_coefficients(expression)
domains = map(ufl.domain.extract_unique_domain, coeffs_and_arguments)

But the new closure-table index is emitted for any SpatialCoordinate/Jacobian on a lower-dimensional domain, whether or not a coefficient or argument lives there. Compiling

domain = ufl.Mesh(basix.ufl.element("Lagrange", "quadrilateral", 1, shape=(2,)))
codomain = ufl.Mesh(basix.ufl.element("Lagrange", "interval", 1, shape=(2,)))
V = ufl.FunctionSpace(domain, basix.ufl.element("Lagrange", "quadrilateral", 1))
v = ufl.TestFunction(V)
a = ufl.inner(ufl.SpatialCoordinate(codomain)[0], v) * ufl.ds(domain=domain)

gives

double x_c0 = coordinate_dofs[quadrilateral_facet_closure_dofs[quadrature_permutation[0]][entity_local_index[0]][0] * 3] * FE3_C0_F_Qdc7[0][0][iq][0] + ...;
...
.needs_facet_permutations = false,

DOLFINx then passes perm = 0 unconditionally — assemble_vector_impl.h and assemble_matrix_impl.h both do

std::uint8_t perm = perms.empty() ? 0 : perms(cell, local_entity);

so the gather uses the parent-local vertex ordering on every facet. That is silently wrong wherever the submesh cell's global orientation is reversed relative to the parent's local sub-entity closure — which is exactly the mismatch this PR exists to correct. No crash, just wrong numbers.

Extending the domain scan to geometric terminals fixes it; I verified this flips the flag to true on the form above and leaves the existing tests passing:

geometry = [
    mt.terminal
    for mt in initial_terminals.values()
    if isinstance(mt.terminal, ufl.geometry.GeometricQuantity)
]
domains = map(ufl.domain.extract_unique_domain, coeffs_and_arguments + geometry)

2. The "-" restriction offset uses the submesh's dof count, not the parent's

In _define_coordinate_dofs_lincomb, offset = num_scalar_dofs * dim is computed from the submesh coordinate element. That is correct only while dof_access is indexed in that element's own layout. Once a closure table is in play, coordinate_dofs is the parent's two-cell buffer, so the offset has to be the parent's. Compiling

domain = ufl.Mesh(basix.ufl.element("Lagrange", "triangle", 1, shape=(2,)))
codomain = ufl.Mesh(basix.ufl.element("Lagrange", "interval", 1, shape=(2,)))
u = ufl.Coefficient(ufl.FunctionSpace(codomain, basix.ufl.element("Lagrange", "interval", 1)))
v = ufl.TestFunction(ufl.FunctionSpace(domain, basix.ufl.element("Lagrange", "triangle", 1)))
a = (ufl.Dx(u, 0) * v)("-") * ufl.dS(domain=domain)

gives

double J1_r1_c0 = coordinate_dofs[triangle_facet_closure_dofs[quadrature_permutation[1]][entity_local_index[1]][0] * 3 + 6] * ...;

+ 6 is 2 * 3 (interval scalar dofs). DOLFINx packs cell + then cell - at num_x_dofs_cell * 3, so for a P1 triangle parent the correct offset is 3 * 3 = 9. As written this reads into the + cell's coordinates.

Either take the offset from the parent element when closure_table is not None, or raise NotImplementedError for restricted mixed-dimensional coordinate gathers — but it shouldn't silently emit this.

3. The codim-2 (ridge) path can't work end-to-end yet

ridge_closure_dofs is indexed with quadrature_permutation[0], but DOLFINx passes an empty perms span for IntegralType::ridge and IntegralType::vertex — only exterior_facet gets facet_perms:

for (auto itg_type : {fem::IntegralType::exterior_facet,
                      fem::IntegralType::vertex, fem::IntegralType::ridge})
{
  md::mdspan<const std::uint8_t, md::dextents<std::size_t, 2>> perms
      = (itg_type == fem::IntegralType::exterior_facet)
            ? facet_perms
            : md::mdspan<const std::uint8_t,
                         md::dextents<std::size_t, 2>>{};

and mesh::Topology exposes no edge-permutation accessor at all (only get_facet_permutations). So row 0 (identity) is always selected and the edge reflection is never applied. test_multiple_mesh_codim2_ridge_gradient_quadrature_permutation passes because it hands quadrature_permutation=1 in by hand, which no DOLFINx caller currently does.

That is a DOLFINx-side gap rather than an FFCx bug, but it's worth saying so explicitly in the PR description — otherwise codim-2 looks supported while quietly producing identity-ordered gathers.

Related: the peak_closure_dofs / entity_type == "vertex" branch in definitions.py has no codegen-level test. I couldn't reach it at all — a point-element form trips the pre-existing assert ttype != "ones" in _define_coordinate_dofs_lincomb first. If it really is unreachable today, I'd rather drop it than ship untested code with a table generator and a special case behind it.

Smaller things

• FFCXBackendAccess.entity_permutation doesn't do what its docstring says. It claims "Used both here and, cross-class, from FFCXBackendDefinitions", but nothing in access.py calls it — the identical three lines still sit inline in table_access, and again in FFCXBackendSymbols.element_table. As it stands the PR adds a third copy rather than removing two. It's purely a symbols operation; I'd put it on FFCXBackendSymbols next to element_table and refactor table_access to use it (definitions already has a self.symbols property).

• ~15 lines duplicated verbatim between integral_generator.generate_geometry_tables and expression_generator.generate_geometry_tables. Both build a closure_table_kinds set that can only ever hold one element, since entity_type is fixed per integral — it's a boolean dressed as a set. A single helper in geometry.py taking the entity type, parent element and integrands and returning the parts to append would collapse both copies.

• _FACET_NPERM/_RIDGE_NPERM are keyed by the parent cell, but the permutation count is a property of the entity cell type: interval → 2, triangle → 6, quadrilateral → 8. Keying on entity_celltype (which _closure_dofs_table already computes) removes both dicts, lets facet_closure_dofs and ridge_closure_dofs collapse into one call, and makes prism/pyramid ridges fall out for free rather than needing a hand-maintained entry each.

• _scalar_basix_element uses (scalar_element,) = set(sub_elements) while _define_coordinate_dofs_lincomb uses coordinate_element._sub_element for the same thing. Worth picking one; sub_elements[0] avoids depending on element hashability.

• write_table(tablename, cellname, coordinate_element=None) — the three new table kinds require the element, so omitting it gives an AttributeError on None rather than a useful message. Either split the signature or validate.

• test_facet_closure_dofs_unsupported_cell_type uses @pytest.mark.xfail(raises=NotImplementedError). With xfail_strict off (the default) this silently XPASSes if the error ever stops being raised. pytest.raises(NotImplementedError) is the right idiom for "this must raise".

• integration_domain_coordinate_element: object | None in CommonExpressionIR — it's always a basix.ufl._ElementBase, and typing it as object gives up the mypy coverage the rest of that NamedTuple has.

• Comment volume. Several comments run three to five times longer than the code they annotate — the 11-line nomenclature note inside _define_coordinate_dofs_lincomb, the 20-line docstring on _reference_vertex_permutations, the four-line justification on _FACET_NPERM. The content is accurate, but it reads as generated narrative; a compression pass would help.

• conftest.py as an import target (from conftest import _reference_vertex_permutations) works, but that isn't what conftest is for — a test/utils.py would be more conventional.

[jorgensd]

@garth-wells I've addressed the various comments made by you/CLAUDE.

1. needs_facet_permutations is not set when the submesh enters only through geometry

Fixed

2. The "-" restriction offset uses the submesh's dof count, not the parent's

Fixed

3. The codim-2 (ridge) path can't work end-to-end yet

For the case regarding ridge support in DOLFINx, I've started this work in FEniCS/dolfinx#3904 and FFCx should move independently of what is done in DOLFINx (to avoid mega-feature PRs).

Furthermore on 3;

Related: the peak_closure_dofs / entity_type == "vertex" branch in definitions.py has no codegen-level test. I couldn't reach it at all — a point-element form trips the pre-existing assert ttype != "ones" in _define_coordinate_dofs_lincomb first. If it really is unreachable today, I'd rather drop it than ship untested code with a table generator and a special case behind it.

This PR now depends on #882 which fixes the accessing of spatial coordinates of codim 0 meshes.

The minor comments

All minor comments have been addressed.

Side note

• Comment volume. Several comments run three to five times longer than the code they annotate — the 11-line nomenclature note inside _define_coordinate_dofs_lincomb, the 20-line docstring on _reference_vertex_permutations, the four-line justification on _FACET_NPERM. The content is accurate, but it reads as generated narrative; a compression pass would help.

I've ensured that all comments in the code now are to the point.
However, I find it amusing that AI states this when a mega comment on the PR is exactly the same.