Add warped GLL points on a pyramid

[garth-wells]

Adds the gll_warped lattice on a pyramid, so that gll_warped
Lagrange elements can be created on pyramids. This addresses the first
half of #268; Isaac's points on a pyramid are still not implemented.

The construction is the warp and blend one that was removed in #271
after intermittent CI failures. It is conforming: the points on the
quadrilateral base are those of a hexahedron face, and the points on
each of the four triangular faces are those of a tetrahedron face, both
to machine precision for degrees up to 16. The interior lattice is
exactly the interior subset of the full lattice.

The failure in #269

The warp function is built by interpolation, so it vanishes at the ends
of the interval only up to rounding. On a triangle or a tetrahedron a
blending factor damps that residual at every vertex, but on a pyramid
the blending factors all degenerate to one at the apex, so the four
residuals accumulate and the apex comes out as, for example,
(-2.2e-16, -2.2e-16, 1.0) with 1 - z equal to -4.4e-16, i.e.
outside the cell.

The polyset on a pyramid is expressed in the collapsed coordinates
x / (1 - z) and y / (1 - z), guarded by an exact comparison of z
against one. Tabulating at that perturbed apex instead of the exact one
changes values by O(10^3) at degree 6, and which side of the guard a
point lands on depends on the compiler, its flags and FMA contraction —
which fits a failure that appeared on some CI runs and not others.

Two changes make the lattice exact:

• the warp is shifted by the linear interpolant of its two endpoint
  residuals, so that it vanishes at the endpoints exactly;
• the barycentric coordinates are formed from the lattice indices
  rather than as 1 - x - z, so a coordinate that should vanish is
  exactly zero and the degenerate branches are selected exactly. This
  replaces a fabs(l) < 1e-12 test, which was below float epsilon
  and so was unreliable in the float instantiation.

All five vertices are then bit-exact and every coordinate is in the
cell, for degrees 1 to 16 in both float and double. The remaining
2e-16 overshoot on the sloping faces is the same as, and smaller than,
that of the existing triangle and tetrahedron warped lattices, and the
pyramid polyset is singular only at the apex.

The warp function is also now evaluated once for the whole lattice
rather than once per lattice point, and the blending factor, which
appears six times, is a helper function.

Why it is worth having

Estimated Lebesgue constants for the pyramid Lagrange element:

degree	equispaced	gll_warped
6	25.1	7.0
8	130.8	13.7
10	943.3	30.0
12	7625.7	74.7

The condition number of the orthonormal polyset Vandermonde at the
lattice points falls from 7.7e4 to 8.4e2 at degree 10.

Tests

The cases disabled in #271 are re-enabled:

• test_lattice.py::test_pyramid for (gll, warp), which checks that
  the surface points agree with those on a triangle and a quadrilateral;
• the pyramid in test_continuity.py, in both the triangular facet and
  the quadrilateral facet tests, which check continuity against
  tetrahedra and hexahedra — the property the issue asks for;
• the pyramid in the gll_warped list in utils.py, which pulls it
  into the dof transformation, interpolation and tabulation tests.

That is 63 more tests. The full Python suite and the C++ unit tests
pass.

🤖 Generated with Claude Code