Create warped lattices on pyramids

[garth-wells]

Points on a pyramid could only be equispaced; any other lattice type raised Non-equispaced points on pyramids not supported yet. This adds the warp-and-blend construction of Hesthaven and Warburton, Nodal Discontinuous Galerkin Methods, pp. 175-180, so that a pyramid can carry GLL points.

Construction

Each point is displaced by a blend of the one-dimensional warps along the eight edges of the pyramid: the four that meet at the apex, and the two directions of the base. The warp arguments and blending weights are accumulated for every point first, so the warp function, which solves a dense system, is evaluated once.

The blending factors are chosen so that

• the points on the quadrilateral base are those of create_quad, and
• the points on each triangular face are those of create_tri_warped,

which is what lets a mesh mix pyramids with hexahedra and tetrahedra. Face selecting functions interpolate between adjacent triangular faces, and are evaluated only away from the edges on which they are undefined.

Checked as point sets, at n = 2, 3, 4, 5, 6 and 8: the base agrees with the quadrilateral lattice and all four triangular faces agree with the triangle lattice, to within 4e-16.

Keeping the apex on the cell

The warp is evaluated through a Lagrange basis, so it vanishes at the ends of an edge only up to rounding. On a triangle a blending factor damps the displacement at a vertex and that residual is harmless, but the pyramid blending factors all degenerate to one at the apex, where the residuals would otherwise accumulate and move the apex off the cell by a few epsilon. Since the pyramid polyset is expressed in the collapsed coordinates x / (1 - z) and y / (1 - z), a point above the apex tabulates to nonsense.

warp_function_zero_ends evaluates the warp at the two endpoints alongside the requested points and subtracts the linear interpolant of the two residuals, so the displacement is exactly zero at both ends. The apex comes back as exactly (0, 0, 1) for every n tested, and every point lies inside the cell.

Other simplex methods

Isaac's points and the centroid method are not implemented on a pyramid and still raise, with messages naming the method. The methods all agree for n <= 3, so simplex_method::none is accepted there.

Tests

Pyramids are re-enabled in test_continuity_triangle_facet and test_continuity_quadrilateral_facet, in the lattice test for warped GLL points, and in the gll_warped entries of the shared element list in test/utils.py. The gll_isaac entries move to a separate loop that still excludes pyramids.

Python suite 10184 passed / 341 skipped / 23 xfailed; C++ suite 22/22; ruff and clang-format clean.

🤖 Generated with Claude Code

[garth-wells]

Duplicate of #1098, which has the same four files with an identical diff. Closing this one.