[garth-wells]

Closes #608.

Adds quadrature::type::witherden_vincent, providing the symmetric quadrature rules of F. D. Witherden and P. E. Vincent, On the identification of symmetric quadrature rules for finite element methods, Comput. Math. Appl. 69 (2015) 1232-1241.

What is included

Cell	Max degree
triangle	20
quadrilateral	21
tetrahedron	10
hexahedron	11
prism	10
pyramid	10

All rules have positive weights and points strictly inside the cell. On an interval the family reduces to Gauss-Legendre, so the type can also be used for facet integrals of 2D cells.

On non-simplex cells the rules need considerably fewer points than the collapsed Gauss-Jacobi rules currently used there:

Cell	Degree	Witherden-Vincent	Gauss-Jacobi
hexahedron	5	14	27
hexahedron	9	58	125
prism	9	60	125
pyramid	7	31	125
tetrahedron	9	59	125

Defaults are unchanged

get_default_rule is untouched, so this PR changes no existing results. Switching the default on hexahedra, prisms and quadrilaterals would be a worthwhile follow-up, but it changes every generated kernel and is better considered separately.

Gauss-Jacobi should probably remain the default on pyramids regardless: the collapsed rule is also exact for the rational pyramid polyset (xy/(1-z) is polynomial in collapsed coordinates), whereas the Witherden-Vincent rules are exact for polynomials only. This is noted in the quadrature::type documentation.

Provenance of the data

The coefficients are the paper's tables, taken from the rule files distributed with PyFR (BSD-3-Clause), and mapped from PyFR's [-1, 1]-based reference cells onto the Basix reference cells. This includes the sheared pyramid map, from PyFR's base [-1, 1]^2 at z = -1 with apex (0, 0, 1) to the Basix base [0, 1]^2 at z = 0 with apex (0, 0, 1).

Each transformed rule was checked against exact rational monomial integrals before the C++ was generated; the worst relative error over all 65 rules, after rounding to double, is 7e-16.

Tests

• test/test_quadrature.py: 89 cases checking that every rule integrates all monomials of total degree at most m exactly, that weights are positive, and that unsupported degrees and cell types raise.
• cpp/test/quadrature.cpp: the same properties plus the weight sum against cell::volume, 8403 assertions.

🤖 Generated with Claude Code