Preserve direct-sum point indices in dual evaluation

[pbrubeck]

Interpolation into a facet-restricted element on a tensor product cell did not sum factorise. The dual basis of such an element is a direct sum over the entities of the cell, and two separate faults stopped TSFC from contracting each component of that sum on its own.

EnrichedElement._dual_evaluation summed each summand over its own points before it returned, and so reported an empty point_indices. The caller could then only contract the whole concatenated evaluation at once, over every summand's points together, which is the contraction that sum factorisation has to split.

gem.optimise.constant_fold_zero tested a table for zeros with numpy.array_equal(node.array, 0). That function compares shapes before it compares values, so it returned False for every table the pass gave it and folded nothing but a scalar zero. The tabulations that vanish identically therefore survived as dense literals. They are common here: each direct-sum component tabulates its interior basis functions at the boundary points of the other components, and those tables are exactly zero. On a hexahedron 30 of the 49 pairs of components contribute nothing, and all 30 were compiled.

Together the two faults cost a power of the degree. A degree-lowering interpolation between facet-restricted Q elements on a hexahedron now costs O(degree^(dim + 1)) again.

Changes

• EnrichedElement._dual_evaluation leaves each summand's point indices free and returns them, instead of contracting them itself. The evaluations stack along the basis index, which is the free index that gem.unconcatenate.unconcatenate splits downstream, so the caller chooses how to contract each block.
• _constant_fold_zero_literal folds on the value of the table rather than on a shape comparison, so a Literal of any shape that holds only zeros becomes a Zero. Indexed and Product already simplify Zero away, so the expressions those tables multiply disappear with them.
• The fold keeps the literal's dtype, so a table of integers does not become a floating point zero. tsfc.fem builds its quadrature orientation permutation maps as uint literals and indexes them.

The fold stays in this deferred pass rather than moving into Literal.__new__. gem.optimise.select_expression needs every expression it is given to share a structural shape, and folding a literal eagerly breaks that for restricted elements on tensor product cells, which is the case at hand. Commit c8cbefa removed the eager fold for that reason and added this pass in its place; test/tsfc/test_tsfc_274.py in Firedrake still covers it.

Tests

• test/finat/test_dual_basis.py checks that the summand point indices stay free.
• test/gem/test_simplify.py checks that a Literal of any shape folds, that a table holding a nonzero entry does not, and that folding a zero table removes the product it appears in. The three shape cases and the product case fail on the previous predicate.
• tests/tsfc/test_dual_evaluation.py::test_dual_argument_is_sum_factorised[cell1-facet] in Firedrake passes again. Its kernel loses 72% of its flops at degree 16 (530974 → 149476), and the measured cost now approaches degree^(dim + 1) from below, at 2.924, 2.977 and 2.990 over the degree pairs 4→8, 8→16 and 16→32, as the cases that already passed do.
• Firedrake's tests/tsfc (395 tests), this repository's test/gem and test/finat (403 tests), and Firedrake's interpolation regression tests (247 tests) all pass. Interpolation between facet-restricted spaces on an extruded hexahedral mesh reproduces a polynomial that both spaces hold, and satisfies <I^T c, u> == <c, I u> to machine precision, at degrees 3 to 6.

AI was used to draft this change: the first commit with Codex, the second with Claude Code (Opus 5).