Function simplex_duffy_gauss_legendre 

pub fn simplex_duffy_gauss_legendre(
    dimension: usize,
    points_per_axis: usize,
) -> Result<QuadratureRule, Diagnostic>

Construct a positive tensor-product Duffy rule on a unit simplex.

Cube coordinates t_i on [0, 1]^d map to simplex coordinates

x_i = t_i * product_{j < i}(1 - t_j).

The Jacobian is product_i (1 - t_i)^(d - i - 1). With n Gauss–Legendre points per cube axis, the rule integrates every total-degree polynomial through 2n - d exactly. Dimension is an explicit part of the rule, so triangles and tetrahedra use the same construction without pretending that their exactness is identical.

Errors

Returns EQ0804 for dimension zero, insufficient or unsupported axis exactness, point-count overflow, or allocation failure.