Cheap and stable quadrature on polyhedral elements

Alvise Sommariva; Marco Vianello

REVIEW 3 major objections 4 minor 24 references

Cheap and stable quadrature on polyhedral elements

T0 review · 3 major / 4 minor · reviewed 2026-08-09 · deepseek-v4-flash

Challenge this review Re-run · record.json Download PDF Read on arXiv ↗ Review JSON

Pith's one-line read This paper shows that on any polyhedron, integrals of polynomials up to degree $n$ can be computed exactly by a single matrix-vector product, with no matrix factorization or inversion, using hyperinterpolation in a bounding box and…

desk verdict A clean application of hyperinterpolation to polyhedral quadrature that is cheap and stable, but the exactness claim rests on unstated face-moment quadrature details. read the letter →

arxiv 2502.03446 v1 pith:OE6VROTS submitted 2025-02-05 math.NA cs.NA

Alvise Sommariva , Marco Vianello This is my paper · ORCID

classification math.NAcs.NA MSC 65D32

keywords quadraturepolyhedralelementshyperinterpolationChebyshevmomentsdivergencetheoremtetrahedra-freepolytopalfinitenumericalintegration

verification ladder T0 review T1 audit T2 compute T3 formal

The pith

A machine-rendered reading of the paper's core claim, the machinery that carries it, and where it could break.

The reading

The paper presents a new way to integrate polynomials over arbitrary polyhedral elements without decomposing them into tetrahedra. The key idea is to represent the integral over the polyhedron as an integral over an enclosing box, using hyperinterpolation with respect to a weighted measure, and then to compute the required Chebyshev moments of the polyhedron via the divergence theorem on its planar faces. The resulting quadrature weights are obtained from a single matrix-vector product that is independent of the polyhedron except through the moment vector, so no conditioning issues arise. The authors prove that the rule becomes stable as the degree grows, in the sense that the sum of absolute weights converges to the volume of the polyhedron, even if some weights are negative. If correct, the method offers a cheap and robust building block for high-order polytopal finite element computations.

What carries the argument

The load-bearing mechanism is the product-like hyperinterpolation formula of Theorem 2.1, quoted from [17] and [22], which converts the weighted integral of $h f$ over the box into a sum over the box quadrature nodes with weights $u_i \sum_j \varphi_j(P_i) m_j$; combined with the choice $h(P) = \mathbf{1}_\Omega(P)/\sigma(P)$, this turns the polyhedral integral into a box integral and yields the stability limit (7). The second ingredient is the computation of the Chebyshev moments $m_j$ via the divergence theorem: each moment is reduced to a sum of face integrals of a univariate primitive of the Chebyshev basis, computed by product Gauss–Legendre quadrature on polygonal faces [19]. Together they make the final weight computation a single element-independent Vandermonde-like matrix times the moment vector.

What would settle it

For a polyhedron with a known volume (e.g., a cube or a tetrahedron), take a fixed degree $n$ and compute the quadrature weights from (6) using the Chebyshev basis and a tensorial Gauss–Chebyshev box rule; if the sum of the weights applied to a random polynomial of exact degree $n$ differs from the true integral by more than a small multiple of machine epsilon, then the claimed exactness fails. A sharper test: increase the number of facets while keeping volume fixed and monitor the moment computation via face quadrature; if the polynomial exactness degrades as facet count grows, the practical accuracy asserted for the face moments is the limiting step.

Watch

Extended reading notes

Core claim

On any polyhedron $\Omega \subset \mathbb{R}^3$, the paper constructs a quadrature rule that is exact for every polynomial in $\mathbb{P}^3_n$ by taking a bounding box $B \supset \Omega$, choosing an absolutely continuous measure $\mathrm{d}\mu = \sigma(P)\,\mathrm{d}P$ on $B$, and forming weights $w_i = u_i \sum_{j=1}^N \varphi_j(P_i) m_j$, where $(X,u)$ is a positive-weight algebraic quadrature exact on $\mathbb{P}^3_{2n}$, $\{\varphi_j\}$ is an orthonormal polynomial basis for $\mu$, and $m_j = \int_\Omega \varphi_j(P)\,\mathrm{d}P$ are moments computed by the divergence theorem over the oriented planar faces. The paper's central claim is that this rule is exact for all $f \in \mathbb{P}^3_n$, requires no matrix factorization or inversion, and satisfies the stability limit $\lim_{n\to\infty}\sum_i |w_i| = \mathrm{vol}(\Omega)$; consequently the stability ratio $\sum_i |w_i| / \mathrm{vol}(\Omega)$ tends to 1, so negative weights do not destroy stability.

Load-bearing premise

The central claim stands on a theorem from earlier work that is quoted rather than re-proved: it requires the polyhedron's volume measure, divided by the chosen box weight, to be square-integrable on the box, and it requires the box cubature to be exact on degree $2n$ polynomials with positive weights; neither condition is verified in detail for the polyhedra used in the experiments.

Editorial extensions

If this is right

For any polyhedron, the rule integrates all polynomials of total degree at most $n$ exactly, with no tetrahedral subdivision.
The weight vector is produced by one matrix-vector product whose matrix depends only on the reference box and the chosen cubature, so it can be reused for all elements in a polyhedral mesh.
The stability ratio $\sum_i |w_i|/\mathrm{vol}(\Omega)$ tends to 1 as $n$ grows, so the formula is stable despite possible negative weights.
For smooth integrands $f \in C^{k+1}$, the error decays like $O(n^{-k})$ on the bounding box by (8), so the method applies to non-polynomial functions as well.
Computational cost is dominated by moment computation, and with a product Chebyshev basis the matrix-vector product can be accelerated by an FFT-based discrete cosine transform.

Reading between the lines

Editorial extensions of the paper, not claims the author makes directly.

With minimal or near-minimal box cubature formulas, the node count could drop from $(n+1)^3$ to about $N = \dim \mathbb{P}^3_n$, which would make the quadrature itself cheaper; the paper notes such formulas exist but treats node count as secondary.
The stability limit (7) suggests the rule may remain usable for nonconvex and multiply connected polyhedra exactly as for convex ones; the numerical tests include such domains but do not isolate the effect of concavity on face moment accuracy.
A natural stress test is to push the face quadrature of [19] to higher degrees on polyhedra with many small facets, where moment errors could dominate; the paper does not report such a breakdown analysis.

Share X Bluesky LinkedIn Reddit HN

Editorial analysis

A structured set of objections, weighed in public.

Desk editor's note, referee report, and a circularity audit.

Desk Editor's Note

read the letter

Quick take: the paper gives a cheap way to build tetrahedra-free quadrature on polyhedra by applying the authors' own nontensorial Clenshaw-Curtis / hyperinterpolation product formula to h = 1_Omega/sigma, with moments computed via the divergence theorem. The twist is that the weight vector is a single matrix-vector product with a polyhedron-independent matrix, and the sum of absolute weights provably tends to vol(Omega). If you do polytopal FEM, that's a genuinely useful trick.

The derivation from (5) to (6) is clean, and the numerical tests on three quite different polyhedra (nonconvex, 760-facet convex, multiply connected) show errors near machine precision and stability ratios <= 2 decreasing toward 1. The cputimes are sensible and compare favorably in order of magnitude with their own NNLS-based TetraFreeQ, though not head-to-head with [10].

The main soft spot is exactly what the stress-test note says: exactness of the rule for P_n^3 is only as good as the moments m_j, and the paper does not state the exactness degree of the face quadrature used to compute them. If those face rules are not exact for the polynomial integrands, both the exactness and the stability table are corrupted by an uncontrolled delta. This may be a non-issue in practice – near machine-precision results on random polynomial tests suggest the moments are accurate for those shapes – but the paper needs to state the quadrature degree and to test on more distorted or many-faceted polyhedra. The authors cite [19,21] for the face formulas and [21, §2.1] for the moment computation, but a self-contained sentence would save the referee a lot of work.

The stability theorem is imported from Sloan and their own earlier paper, and the limit (3) is nontrivial for discontinuous h. That's acceptable for a research note, but they should at least verify the hypotheses (h in L^2_mu) explicitly for the Chebyshev measure; it's true, but a remark would help.

Minor: the code is preliminary, no commit hash; the comparison with Langlois et al. [10] is qualitative, not benchmarked. Both are fixable.

Verdict: this is a solid, useful note, not a breakthrough. The idea is simple and correct, and the implementation is genuinely cheap. It deserves a serious referee, who should push for details on the face quadrature and a direct comparison with [10] on the same polyhedra. I'd cite it if I worked in this area.

Referee Report

3 major / 4 minor

Summary. The manuscript proposes a tetrahedra-free quadrature rule for integrating polynomials over arbitrary polyhedral domains in three dimensions. After embedding the polyhedron Ω in a Cartesian bounding box B with measure σ dP, the authors apply a hyperinterpolation product formula from [17,22] to obtain weights w_i = u_i Σ_j φ_j(P_i) m_j, where the moments m_j = ∫_Ω φ_j dP are computed via the divergence theorem as sums of face integrals. The resulting rule is exact on P^3_n, does not require matrix factorization or inversion, and has a theoretical stability ratio ∑_j |w_j|/vol(Ω) tending to 1. Numerical experiments on three polyhedra of different topology report relative errors close to machine precision, cputimes between 10^-3 and about 1 second, and stability ratios below 2.

Significance. If the result holds, it is a genuinely cheap and stable tetrahedra-free quadrature tool for polytopal FEM: the Vandermonde-like matrix is element-independent, only one matrix-by-vector product is needed, and the stability bound is inherited from hyperinterpolation theory. The formula has no fitted or free parameters, and the authors provide a public Matlab implementation [18]. The numerical section is reproducible and includes diverse polyhedra (nonconvex, many-faceted, multiply connected), although the tests are limited in scope and benchmarked only against the authors' earlier method [21].

major comments (3)

[Section 3, step (ii) and Eq. (6)] The exactness of the formula on P^3_n is contingent on computing the moments m_j exactly, but the paper cites [19]/[21] for the face integrals without stating the exactness degree, quadrature order, or behavior of those face rules on distorted or many-faceted faces. If m_j is computed with perturbation δ_j, the quadrature error for f ∈ P^3_n is Σ_j δ_j ∫_B φ_j f dμ, which is generally nonzero. Please specify the exactness properties of the face quadrature, prove that the implemented moments are exact for the Chebyshev basis up to degree n, or derive an explicit moment-error bound, and test on strongly skewed and many-faceted polyhedra.
[Section 3, Figure 2 and Table 2] The numerical evidence tests only random polynomials of the special form (a_k x + b_k y + c_k z + d_k)^n and reports average logarithmic errors. Such integrands can miss moment perturbations, and averages are not informative about worst-case exactness on the whole space P^3_n. Report maximum relative errors over a basis of P^3_n, or over a large set of dense random polynomials, with a high-accuracy reference integral, and report worst-case stability ratios over several random rotations or scalings of the polyhedra.
[Section 2, Theorem 2.1] The stability limit (7) and the error bound (8) are inherited from a cited theorem whose proof is not reproduced. Since h = I_Ω/σ is discontinuous on B, the paper would be stronger if it verified the L^2_μ hypothesis and the convergence of the hyperinterpolation projections in this specific polyhedral setting, or gave a concise proof sketch for the limit (3) under h = I_Ω/σ. This would turn the stability claim from an external citation into a checked hypothesis of the present application.

minor comments (4)

[Throughout] There are several typographical and stylistic errors, including 'constantly growing', 'exahustivity', 'diﬀerently', 'non negligible', and 'Tchakaloﬀ'; these should be corrected.
[Section 3, Table 1] The claimed speed advantage over [21] ('up to one order of magnitude lower') is not documented in this paper; give a side-by-side timing comparison or state explicitly where the comparison is reported.
[Eq. (8)] The notation E_n(f;B) is defined, but the asymptotic statement '∼ 2 vol(Ω) E_n' should be phrased more carefully: for finite n the bound is (vol(Ω) + Σ_i |w_i|) E_n(f;B), and the asymptotic form uses (7) as n → ∞.
[Section 3, step (i)] The text mentions minimal or near-minimal cubature formulas as possible low-cardinality choices, but the implementation only uses tensor Gauss-Chebyshev nodes; if the low-cardinality variants are not tested, they should be described as future work rather than as part of the tested method.

Circularity Check

0 steps flagged · score 0.0 of 10

No significant circularity: formula (6) is a direct application of a proven hyperinterpolation identity with no fitted parameters, and the self-citations are to independent published results.

full rationale

The central construction is not circular. Equation (6) follows from Theorem 2.1, quoted from [17] and [22], by setting h = I_Omega/sigma on a bounding box and using a positive quadrature rule exact on P^3_{2n}. The weights are defined through the moments m_j = ∫_Omega phi_j dP; if these moments are computed exactly, the identity ∫_Omega f dP = Σ w_i f(P_i) is exact for every f in P^3_n by the hyperinterpolation product-like formula. No parameter is fitted to the target integrals, and the stability ratios in Table 2 are compared with the theoretical limit (7), not used to define the rule. The reliance on the authors' prior face quadrature [19] and [21] is a reduction to a lower-dimensional, independently published quadrature problem, not a tautological reuse of the target result; the benchmark reference [21] is a different tessellation-based method, so using it as a comparison does not make the test circular. The paper's dependence on Theorem 2.1 from [17] and [22] is a citation of a proven, parameter-free result whose assumptions do not include the polyhedral quadrature target, so it is real evidence rather than a self-citation chain. Any concern that the implemented face-moment computation may not be certified exact at high degree is a correctness/robustness issue, not a circularity of the derivation.

Assumptions & free parameters 0 free parameters · 6 assumptions · 0 invented entities

The method introduces no new entities and no fitted parameters: nodes and weights come from known Gauss-Chebyshev rules, the basis is the classical product Chebyshev basis, and moments are computed via the divergence theorem. All building blocks are prior mathematical results; the contribution is their assembly and numerical validation. The main burdens are the imported stability theorem and the delegated face-quadrature accuracy.

assumptions (6)

standard math Theorem 2.1: the product-like hyperinterpolation formula (1) and its stability limit (3), quoted from [17] and [22]
Invoked in Section 2 as the engine of the method; equations (6) and (7) are direct consequences. The proof is not reproduced in this paper.
standard math The divergence theorem converts volume moments into sums of integrals over oriented planar faces of the polyhedron (Section 3, step ii)
Used to compute m_j = integral over Omega of phi_j dP via flux integrals; requires the polyhedron boundary to be a closed oriented piecewise planar surface.
standard math Existence of positive-weight quadrature rules on the cube exact on P^3_{2n}, e.g., the tensorial Gauss-Chebyshev rule with (n+1)^3 nodes
Used in step (i) of the implementation; the correctness of the underlying box rule is assumed, not shown.
domain assumption h = I_Omega/sigma lies in L^2_mu(B), i.e., 1/sigma is integrable on Omega
Stated in Section 2 as the condition for writing equation (5); satisfied for Lebesgue and Chebyshev measures on bounded polyhedra.
domain assumption The face integrals in the moment computation are computed accurately by the quadrature rules of [19]/[21]
Section 3 delegates the face quadrature to prior work without an accuracy analysis at high polynomial degree or for many-faceted polyhedra.
standard math Multivariate Jackson theorem for the error bound (8) for non-polynomial integrands (cited to [15])
Used only in Remark 2.1 to justify the extension of the rule to functions well approximated by polynomials on the bounding box.

how reviews work

0 comments

Cite this review

Pith. "Pith review of Cheap and stable quadrature on polyhedral elements." pith.science (2026). https://pith.science/paper/OE6VROTS

@misc{pith2026250203446,
  author       = {Pith},
  title        = {Pith review of: Cheap and stable quadrature on polyhedral elements},
  year         = {2026},
  howpublished = {\url{https://pith.science/paper/OE6VROTS}},
  note         = {Machine review of arXiv:2502.03446}
}

read the original abstract

We discuss a cheap tetrahedra-free approach to the numerical integration of polynomials on polyhedral elements, based on hyperinterpolation in a bounding box and Chebyshev moment computation via the divergence theorem. No conditioning issues arise, since no matrix factorization or inversion is needed. The resulting quadrature formula is theoretically stable even in the presence of some negative weights.

Figures

Figures reproduced from arXiv: 2502.03446 by the authors.

**Figure 1.** Examples of polyhedral domains. Left: Ω1 (nonconvex, 20 facets); Center: Ω2 (convex, 760 facets); Right: Ω3 (multiply connected, 20 facets). 4 6 8 10 12 14 16 18 20 10-14 10-12 10-10 4 6 8 10 12 14 16 18 20 10-15 10-14 10-13 10-12 10-11 4 6 8 10 12 14 16 18 20 10-15 10-14 10-13 10-12 10-11 [PITH_FULL_IMAGE:figures/full_fig_p005_1.png] view at source ↗

**Figure 2.** Relative errors E(gk) of the tetrahedra-free rule over 200 polynomial integrands of the form gk = (akx+bky+ckz+dk) n on the three polyhedra of [PITH_FULL_IMAGE:figures/full_fig_p005_2.png] view at source ↗

Discussion (0). Continue with ORCID to comment.

Reference graph

Works this paper leans on

24 extracted references · 24 canonical work pages

[10]

Langlois, T

C. Langlois, T. van Putten, H. B´ eriot, E. Deckers, Frug al numerical integration scheme for polytopal domains, Eng . Comput. (2024)

work page 2024
[18]

Sommariva, CheapQ, Matlab codes for cheap and stable tetrahedra-free quadrat ure on polyhedral elements, https://github.com/alvisesommariva/CheapQ

A. Sommariva, CheapQ, Matlab codes for cheap and stable tetrahedra-free quadrat ure on polyhedral elements, https://github.com/alvisesommariva/CheapQ

work page
[21]

Sommariva, M

A. Sommariva, M. Vianello, TetraFreeQ: tetrahedra-fr ee quadrature on polyhedral elements, Appl. Numer. Math. 20 0 (2024), 389–398

work page 2024
[19]

Sommariva, M

A. Sommariva, M. Vianello, Product Gauss cubature over polygons based on Green’s integration formula, BIT Numeric al Mathematics 47 (2007), 441–453

work page 2007
[1]

Antonietti, P

P.F. Antonietti, P. Houston, G. Pennesi, Fast numerical integration on polytopic meshes with applications to disco ntinuous Galerkin ﬁnite element methods, J. Sci. Comput. 77 (2018) 33 9–1370

work page 2018
[2]

Boroomand, N

B. Boroomand, N. Niknejadi, Adaptive quadrature/cubat ure rule: Application to polytopes, Comp. Methods Appl. Mec h. Engrg. 403 (2023), 115726

work page 2023
[3]

Cangiani, E.H

A. Cangiani, E.H. Georgoulis, P. Houston, hp-version di scontinuous Galerkin methods on polygonal and polyhedral m eshes, Math Models Methods Appl. Sci. 24 (2014), 2009–2041

work page 2014
[4]

Chin, J.B

E.B. Chin, J.B. Lasserre, N.Sukumar, Numerical integra tion of homogeneous functions on convex and nonconvex polyg ons and polyhedra, Comput. Mech. 56 (2015), 967–981

work page 2015

Show all 24 references

[5]

E.B. Chin, N. Sukumar, An eﬃcient method to integrate pol ynomials over polytopes and curved solids, Comput. Aided Geom. Des. 82 (2020), 101914

2020
[6]

Cools, An Encyclopaedia of Cubature Formulas, J

R. Cools, An Encyclopaedia of Cubature Formulas, J. Comp lexity 19 (2003), 445–453 (https://nines.cs.kuleuven.be/research/ecf/)

2003
[7]

De Marchi, M

S. De Marchi, M. Vianello, Y. Xu, New cubature formulae an d hyperinterpolation in three variables, BIT Numer. Math. 49 (2009), 55–73

2009
[8]

Dunkl and Y

C.F. Dunkl and Y. Xu, Orthogonal Polynomials of Several V ariables, Encyclopedia of Mathematics and its Application s, vol. 81, Cambridge University Press, Cambridge, 2001. 6

2001
[9]

Hubrich, A

S. Hubrich, A. D¨ uster, Numerical integration for nonli near problems of the ﬁnite cell method using an adaptive sche me based on moment ﬁtting, Comput. Math. Appl. 77 (2019), 1983– 1997

2019
[11]

Mason, D.C

J.C. Mason, D.C. Handscombe, Chebyshev Polynomials, C hapman&Hall/CRC, 2002

2002
[12]

The Mathworks, alphaShape, Polygons and polyhedra from points in 2-D and 3-D, https://www.mathworks.com/help/matlab/ref/alphashape.html
[13]

M¨ oller, Kubaturformeln mit minimaler Knotenzah l, Numer

H.M. M¨ oller, Kubaturformeln mit minimaler Knotenzah l, Numer. Math. 25 (1975/76), no. 2, 185–200

1975
[14]

Mousavi, N

S.E. Mousavi, N. Sukumar, Numerical integration of pol ynomials and discontinuous functions on irregular convex p olygons and polyhedrons, Comput. Mech. 47 (2011) 535–554

2011
[15]

Ple´ sniak, Multivariate Jackson Inequality, J

W. Ple´ sniak, Multivariate Jackson Inequality, J. Com put. Appl. Math. 233 (2009), 815–820

2009
[16]

Potts, T

D. Potts, T. Volkmer, Fast and exact reconstruction of a rbitrary multivariate algebraic polynomials in Chebyshev form, 2015 International Conference on Sampling Theory and Appli cations (SampTA), IEEE, pp. 392–396

2015
[17]

Sloan, Interpolation and Hyperinterpolation ove r General regions, J

I.H. Sloan, Interpolation and Hyperinterpolation ove r General regions, J. Approx. Theory 83 (1995) 238–254

1995
[20]

Sommariva, M

A. Sommariva, M. Vianello, Computing approximate Feke te points by QR factorizations of Vandermonde matrices, Com - put. Math. Appl. 57 (2009), 1324–1336

2009
[22]

Sommariva, M

A. Sommariva, M. Vianello, R. Zanovello, Nontensorial Clenshaw-Curtis cubature, Numer. Algorithms 49 (2008), 40 9–427

2008
[23]

Sudhakar, J.P

Y. Sudhakar, J.P. Moitinho De Almeida, W.A. W all, An acc urate, robust, and easy-to-implement method for integrati on over arbitrary polyhedra: application to embedded interfa ce methods, J. Comput. Phys. 273 (2014), 393–415

2014
[24]

Sudhakar, A

Y. Sudhakar, A. Sommariva, M. Vianello, W.A. W all, On th e use of compressed polyhedral quadrature formulas in embedded interface methods, SIAM J. Sci. Comput. 39 (2017) B 571–B587. 7

2017

Pith tools

Reviewed August 9, 2026 · model on record in the stance chip above.

[1] [10]

Langlois, T

C. Langlois, T. van Putten, H. B´ eriot, E. Deckers, Frug al numerical integration scheme for polytopal domains, Eng . Comput. (2024)

work page 2024

[2] [18]

Sommariva, CheapQ, Matlab codes for cheap and stable tetrahedra-free quadrat ure on polyhedral elements, https://github.com/alvisesommariva/CheapQ

A. Sommariva, CheapQ, Matlab codes for cheap and stable tetrahedra-free quadrat ure on polyhedral elements, https://github.com/alvisesommariva/CheapQ

work page

[3] [21]

Sommariva, M

A. Sommariva, M. Vianello, TetraFreeQ: tetrahedra-fr ee quadrature on polyhedral elements, Appl. Numer. Math. 20 0 (2024), 389–398

work page 2024

[4] [19]

Sommariva, M

A. Sommariva, M. Vianello, Product Gauss cubature over polygons based on Green’s integration formula, BIT Numeric al Mathematics 47 (2007), 441–453

work page 2007

[5] [1]

Antonietti, P

P.F. Antonietti, P. Houston, G. Pennesi, Fast numerical integration on polytopic meshes with applications to disco ntinuous Galerkin ﬁnite element methods, J. Sci. Comput. 77 (2018) 33 9–1370

work page 2018

[6] [2]

Boroomand, N

B. Boroomand, N. Niknejadi, Adaptive quadrature/cubat ure rule: Application to polytopes, Comp. Methods Appl. Mec h. Engrg. 403 (2023), 115726

work page 2023

[7] [3]

Cangiani, E.H

A. Cangiani, E.H. Georgoulis, P. Houston, hp-version di scontinuous Galerkin methods on polygonal and polyhedral m eshes, Math Models Methods Appl. Sci. 24 (2014), 2009–2041

work page 2014

[8] [4]

Chin, J.B

E.B. Chin, J.B. Lasserre, N.Sukumar, Numerical integra tion of homogeneous functions on convex and nonconvex polyg ons and polyhedra, Comput. Mech. 56 (2015), 967–981

work page 2015