Exact Integration for singular Zienkiewicz and Guzman-Neilan Finite Elements with Implementation

Johannes Storn; Lars Diening; Tabea Tscherpel

REVIEW 1 major objections 3 minor 1 cited by

Exact Integration for singular Zienkiewicz and Guzman-Neilan Finite Elements with Implementation

T0 review · 1 major / 3 minor · reviewed 2026-08-12 · deepseek-v4-flash

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

Pith's one-line read A recursive integration rule computes exact integrals of rational functions on triangles, enabling quadrature-free implementations of rational finite elements.

desk verdict Solid, useful paper: exact recursive integration for rational finite elements, with working code; the central claim holds, and the flagged Lemma 3.2 overstatement does not touch the algorithm. read the letter →

arxiv 2411.09485 v2 pith:K2NBILVC submitted 2024-11-14 math.NA cs.NA

Lars Diening , Johannes Storn , Tabea Tscherpel This is my paper · ORCID

classification math.NAcs.NA MSC 65D3265N3076D0776M10

keywords exactquadraturerationalfunctionssingularZienkiewiczelementGuzmán-NeilanbiharmonicequationStokesequationspressurerobustnessMATLABimplementation

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 establishes the first exact integration rule for the class of rational functions $R_\alpha^\beta = \lambda^\alpha/(1-\lambda)^\beta$ on triangles in 2D, the functions that appear in rational finite element bases. It gives a recursive algorithm that returns the exact integral mean $I(\alpha,\beta)$ for every multi-index pair with finite integral, and proves termination. This matters because such rational elements—the singular Zienkiewicz element for biharmonic problems and the lowest-order Guzmán–Neilan element for Stokes flow—achieve $H^2$-conformity or exact divergence-freeness with only 12 degrees of freedom per triangle, but had previously relied on inexact Duffy-transform quadrature. The paper shows numerically that inexact quadrature can stall eigenvalue convergence and destroy pressure robustness, so the exact rule is what makes these elements fully trustworthy.

What carries the argument

The central object is the rational function $R_\alpha^\beta(\lambda)=\lambda^\alpha/(1-\lambda)^\beta$ on the 2D simplex, where $\lambda$ are the barycentric coordinates and $\alpha,\beta\in\mathbb{N}_0^3$. The argument is carried by three recurrence identities (Lemma 3.8) that express $R_\alpha^\beta$ as a linear combination of rational functions with lower total denominator exponent $|\beta|$ or lower numerator exponent component $\alpha_0$, together with two closed-form base cases: the polynomial integral formula (Lemma 3.3) and the single-pole formula for $\beta=(0,0,\beta_2)$ (Lemma 3.4). Algorithm 2 uses the integrability threshold of Lemma 3.2 as its first check, then routes the recursion through the identities and, when $\alpha_0=0$, through the one-dimensional closed-form evaluation of Algorithm 1.

What would settle it

Implement Algorithm 2 and compare its output on all multi-indices with entries up to, say, 3 to a high-precision Duffy-transform quadrature (e.g., 50-point Gauss rules per direction) on the reference triangle; a single mismatch in either the finite/infinite classification or the numerical value beyond rounding error would falsify the recursion or the integrability threshold. More directly, proving that the two terms in identity (3.7) cancel for some $(\alpha,\beta)$ would break the sharpness argument of Lemma 3.2.

Watch

Extended reading notes

Core claim

The central claim is that the integral mean of $R_\alpha^\beta$ over any triangle can be evaluated exactly by a finite recursion. The recursion uses three algebraic identities that rewrite the integrand into sums of rational functions with smaller parameters, and it is organized by an exact regularity condition: $I(\alpha,\beta)<\infty$ if and only if $\lVert\alpha+\beta\rVert_\infty \le |\alpha|+1$ (Lemma 3.2). Theorem 3.10 proves Algorithm 2, which implements this recursion, terminates and returns the exact value in all finite cases. As a consequence, the paper provides simple MATLAB implementations of the singular Zienkiewicz and the lowest-order Guzmán–Neilan elements whose system matrices are computed exactly, and its numerical experiments demonstrate that inexact quadrature causes stagnation of the biharmonic eigenvalue error and a loss of pressure robustness in the Stokes discretization, both of which are remedied by the exact integration.

Load-bearing premise

The paper's recursion relies on the exactness of its integrability test: it assumes that no exceptional rational function with $\lVert\alpha+\beta\rVert_\infty = |\alpha|+1$ is actually integrable because its two singular pieces cancel. If such a cancellation existed, the algorithm would wrongly return infinity for a finite integral.

Editorial extensions

If this is right

The singular Zienkiewicz element can be assembled with exactly computed biharmonic stiffness and mass matrices during an offline phase, eliminating quadrature error from the discretization.
The Guzmán–Neilan Stokes element preserves its exact divergence-free constraint when system matrices are formed with the exact rule, so discrete solutions remain exactly divergence-free.
Inexact quadrature is shown to cause stagnation of eigenvalue error on uniform meshes and a reduced convergence order on graded meshes; exact integration removes this barrier.
The provided MATLAB template demonstrates a general pipeline for implementing rational finite elements, so other rational elements can be studied without developing new quadrature.

Reading between the lines

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

The same recursion identities may extend to rational functions in higher dimensions or with denominator factors that are products of general affine functions, since the algebraic decompositions in Lemma 3.8 are dimension-agnostic in form.
If the numerical stagnation observed with inexact quadrature is generic, then any rational finite element method intended for eigenvalue or long-time structure-preserving simulations should use exact integration from the outset.
The sharp integrability threshold is a reusable analytical tool: it provides a simple criterion for designing rational bubbles with prescribed Sobolev regularity on simplices.
Memoizing intermediate values of $I(\alpha,\beta)$ and exploiting symmetries, as the performance remark suggests, should make the exact rule fast enough for production-scale use.

Share X Bluesky LinkedIn Reddit HN

Signed reviews

No signed human review yet.

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

Here's my read on Diening–Storn–Tscherpel. The core contribution is real: Algorithm 2 recursively integrates rational functions R_α^β on triangles exactly, terminating via a clean integrability criterion (3.10), and the authors back it with MATLAB code and numerical experiments. I'm fairly convinced it is the first exact quadrature for this class; the Duffy-transform route in Schneier's paper is inexact, and the tensor-product NURBS approaches do not cover these non-tensor rationals. The recursion identities in Lemma 3.8 are simple but effective, the termination argument in Theorem 3.10 is sound, and the derivation is from first principles—no fitting parameters anywhere.

The one substantive gap is Lemma 3.2. The sharpness proof for m ≥ 1 leans on a claim about 'linearly independent singularities' that is not demonstrated, and as stated the lemma is actually false for m ≥ 1: take α = (1,0,0), β = 0; then R = λ_1 is in W^{1,2} even though |α| − ∥α+β∥∞ = 0 < 1. The stress-test note is right that this overstatement does not damage the paper, because Algorithm 2 and Proposition 3.9 only invoke the m = 0 case (3.10), which is correct and can be verified directly by checking the exponent near each vertex. The authors should tighten Lemma 3.2 to m = 0 or add the direct proof; it is a minor patch, not a load-bearing flaw.

The numerical experiments are honest and useful: they show stagnation of eigenvalue error with inexact quadrature and destruction of pressure robustness in the Guzmán–Neilan Stokes discretization. This is a practical finding that justifies the exact rule. The code is available via Code Ocean, though I would have liked a commit hash; that is a nit. Citation patterns look appropriate—Schneier as the inexact baseline, Guzmán–Neilan and Ciarlet for element properties, standard quadrature literature.

This paper will be most valuable to anyone implementing rational finite elements or working on structure-preserving discretizations; it removes a real practical barrier. I would send it to a serious referee. The central contribution is formally grounded, reproducible, and the numerical claims are independently checkable. My recommendation: accept after a minor revision that fixes the scope of Lemma 3.2.

Referee Report

1 major / 3 minor

Summary. The paper develops an exact quadrature rule for rational functions of the form λ^α/(1−λ)^β on triangles. The central algorithmic contribution is a recursion (Algorithm 2) that computes the integral mean I(α,β) for all multi-indices with finite integral, with termination guaranteed by a sharp integrability criterion. The method is applied to the rational bubble functions appearing in the singular Zienkiewicz and Guzmán–Neilan finite elements, and MATLAB implementations are provided. Numerical experiments demonstrate that inexact quadrature causes stagnation in biharmonic eigenvalue computations and destroys pressure robustness in the Guzmán–Neilan Stokes discretization.

Significance. If the central claim holds, the paper fills a genuine gap in the literature: no exact quadrature for multivariate non-tensor-product rational functions on triangles was previously available. The derivation is parameter-free and the recursion identities are proven in detail, including a termination argument. The code is publicly available and the numerical experiments are reproducible and directly support the qualitative claims about the effects of inexact quadrature. The main theorem (Algorithm 2) is sound; the only mathematical defect is an overstatement in Lemma 3.2 that does not affect the main result but must be corrected.

major comments (1)

[§3.1, Lemma 3.2] The claimed characterization of W^{m,p}-regularity is false for m ≥ 1. A counterexample is R^(1,0,0)_(0,0,0) = λ_1 on the reference triangle: for m = 1 and p = 2, one has |α| − ∥α+β∥∞ = 0, which violates the stated condition |α| − ∥α+β∥∞ > m − 2/p, yet λ_1 ∈ W^{1,2}(T). The sharpness argument invoking 'linearly independent singularities' is not valid in general, because cancellations under differentiation can occur. Since Algorithm 2 and Proposition 3.9 only use the m = 0 case, namely condition (3.10), this error does not affect the exact integration claim; nevertheless, the lemma should be corrected (for example, restricted to m = 0) or replaced by a valid statement for m ≥ 1.

minor comments (3)

[§3.2.3, Algorithm 2] The sorting instruction reads 'Sort ((αj,βj))j=1,2 such that β0 ≤ β1 ≤ β2', but the indices should range over j = 0, 1, 2; please clarify the notation.
[§4, Eq. (4.2)] The notation (DF)^{−⊤} is inconsistent with the code in Figure 4, line 6, which uses (DF)^{-1}, and with the preceding derivation Dλref = Dλ DF; the formula should read (DF)^{-1}.
[§3.2.3, Theorem 3.10] The termination proof states that in the final recursion 'α0 or |β| is reduced', but the earlier branches using Proposition 3.9(c) also reduce α0. For completeness, the proof should specify a global decreasing measure, such as α0 + |β|, for all branches of Algorithm 2.

Circularity Check

0 steps flagged · score 0.0 of 10

No significant circularity: the exact integration recursion is derived from first principles, the termination argument is internally verifiable, and the load-bearing citations are external or non-essential.

full rationale

The central claim, Algorithm 2's exact computation of I(alpha,beta), is derived from Gamma-function identities (Lemma 3.3, Lemma 3.4), elementary one-dimensional integral recursions (Lemmas 3.5-3.6), and algebraic rational-function identities (Lemma 3.8). Each recursion is proven directly in the paper, and the finiteness check used by Algorithm 2 is the m=0 case of Lemma 3.2, which can be verified independently by local integrability analysis near vertices; it is not imported from a self-citation or fitted to data. The finite element properties of the singular Zienkiewicz and Guzman-Neilan spaces are cited from Ciarlet and Guzman-Neilan, but those citations are not used to justify the quadrature formula itself. The numerical experiments compare exact integration to inexact quadrature and are checked against known exact solutions; the observed stagnation and pressure-robustness loss are presented as numerical evidence, not as deductions from the integration formula. Self-citations such as [DST22] and [ST24] appear only in background comparisons and are not load-bearing. The only notable mathematical issue is that Lemma 3.2 states a regularity criterion for all m >= 1 and justifies sharpness by asserting that the two singular terms in (3.7) cannot cancel; that sharpness assertion is not fully demonstrated and the statement appears overbroad for m >= 1. However, Algorithm 2 and Proposition 3.9 invoke only the m=0 condition (3.10), whose sharpness follows directly from explicit vertex-neighborhood estimates, so this is a correctness concern rather than a circularity. No step reduces by construction to its own inputs, no fitted parameter is renamed as a prediction, and no load-bearing uniqueness theorem is imported from the authors' prior work.

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

The paper introduces no free parameters or invented entities. It relies on standard analytic facts and on two finite element theories from the literature (Zienkiewicz, Guzmán-Neilan), which are used as black boxes for the applications but not for the integration formula itself.

assumptions (5)

domain assumption Z_s(T) has 12 local degrees of freedom and defines a C^1-conforming global space.
Used to set up the biharmonic discretization in Section 4.2; proved in [Cia02, Thm. 6.1.4] and cited, not reproven.
domain assumption The Guzmán-Neilan spaces V,Q form an exact discrete de Rham complex and are inf-sup stable.
Used for the Stokes discretization in Section 4.3; proved in [GN14b] and cited, not reproven.
standard math Gamma function identities: Γ(s)=∫_0^∞ t^{s-1}e^{-t}dt and Γ(s)=(s-1)! for integer s.
Used in the proofs of Lemmas 3.3 and 3.4 to evaluate Dirichlet-type integrals.
standard math Polygamma identities, e.g. ∫_0^1 y^α log(y)/(1-y) dy = -ψ^(1)(α+1).
Used in Lemma 3.6 to evaluate the base case J(α1,α2,1,1); taken from [AS64].
standard math Lebesgue integrability near vertices: integrals of r^s over a 2D triangle are finite iff s > -2.
Used in Lemma 3.2 to derive the sharp threshold (3.10).

how reviews work

0 comments

Cite this review

Pith. "Pith review of Exact Integration for singular Zienkiewicz and Guzman-Neilan Finite Elements with Implementation." pith.science (2026). https://pith.science/paper/K2NBILVC

@misc{pith2026241109485,
  author       = {Pith},
  title        = {Pith review of: Exact Integration for singular Zienkiewicz and Guzman-Neilan Finite Elements with Implementation},
  year         = {2026},
  howpublished = {\url{https://pith.science/paper/K2NBILVC}},
  note         = {Machine review of arXiv:2411.09485}
}

read the original abstract

We develop a recursive integration formula for a class of rational polynomials in 2D. Based on this, we present implementations of finite elements that have rational basis functions. Specifically, we provide simple Matlab implementations of the singular Zienkiewicz and the lowest-order Guzman-Neilan finite element in 2D.

Figures

Figures reproduced from arXiv: 2411.09485 by the authors.

**Figure 1.** Degrees of freedom in the Zienkiewicz element and the reduced singular Zienkiewicz (left) as well as in the singular Zienkiewicz (right) element. Dots denote point evaluation, circles denote evaluation of the gradient, and straight lines denote the evaluation of the normal derivative. λj (vk) = δj,k for all j, k = 0, 1, 2. We denote the vector of barycentric coordinates by λ = (λ0, λ1, λ2) ⊤ ∈ [0, 1]3 . We shall f… view at source ↗

**Figure 2.2.** 16]. For this reason, let us first recall the Zienkiewicz element. It is globally [PITH_FULL_IMAGE:figures/full_fig_p004_2_2.png] view at source ↗

**Figure 2.** Degrees of freedom in the Guzm´an–Neilan elements for the vector-valued velocity (left) and the pressure (right). properties. For functions g : R 2 → R we define curl g := ∂x2 g −∂x1 g = R ∇g with R := 0 1 −1 0 . (2.8) As above let bT := λ0λ1λ2 be the element bubble function, let bfj := λj+1λj+2 be the edge bubble functions, and let Bfj be the rational bubble functions as in (2.5), for j = 0, 1, 2. The local … view at source ↗

Figures from the paper (7 more)

**Figure 3.** Figure 3: Degrees of freedom in the reduced Guzm´an–Neilan elements for the vector-valued velocity (left) and the pressure (right). Dots denote point evaluations and arrows evaluation of the normal component. Due to the representation of the velocity space in (2.10) we have that…

**Figure 4.** Figure 4: MATLAB loop over all T ∈ T to obtain the matrix A and right-hand side b in (4.7). all i = 1, 2, 3 and j = 1, . . . , 12, cf. line 13 in [PITH_FULL_IMAGE:figures/full_fig_p023_4.png]

**Figure 5.** Figure 5: Modification replacing lines 15–17 in [PITH_FULL_IMAGE:figures/full_fig_p024_5.png]

**Figure 6.** Figure 6: Computation of the local contributions to assemble the matrices A, B, and b in (4.19). This means PT contains the entries where both basis functions are affine polynomials, RT contains the entries where both basis functions are curls of a rational functions, and MT c…

**Figure 7.** Figure 7: Modification replacing lines 27–29 in [PITH_FULL_IMAGE:figures/full_fig_p029_7.png]

**Figure 8.** Figure 8: Distance |λh−λh|/λh on the square with uniform refinement (left) and L-shaped domain with graded mesh (right), where λh resulted from computations with the quadrature rule in (5.1) and various n. is obtained by the standard AFEM loop with D¨orfler marking and bulk para…

**Figure 9.** Figure 9: Velocity error ∥∇(u−uh)∥L2(Ω) for the Guzm´an–Neilan FEM with uniformly refined triangulation T with #T = 8192 elements and inexact quadrature in (5.1) with n 1D Gauß points. The dashed line represents the velocity error of the Taylor–Hood FEM. [ACFK02] J. Alberty, C. …

Discussion (0). Continue with ORCID to comment.

Forward citations

Cited by 1 Pith paper

Reviewed papers in the Pith corpus that reference this work. Sorted by Pith novelty score. Full citation record

Scott-Vogelius element and iterated penalty method for inhomogeneous Dirichlet boundary conditions
math.NA 2025-09 conditional novelty 6.0 of 10

For inhomogeneous Dirichlet Stokes data, quasi-optimal and pressure-robust error estimates for the Scott-Vogelius element and convergence of the iterated penalty method are established, contingent on a compatibility c...

Reference graph

Works this paper leans on

4 extracted references · 3 canonical work pages · cited by 1 Pith paper

[36]

Application of Randomized Quadrature Formulas to the Finite Element Method for Elliptic Equations

doi: 10.1090/S0025-5718-2013-02753-6 . [GN18] J. Guzm´ an and M. Neilan. “inf-sup stable finite elements on barycentric refinements producing divergence-free approximations in arbitrary dimen- sions”. SIAM J. Numer. Anal. 56.5 (2018), pp. 2826–2844. doi: 10.1137/ 17M1153467. [HCB05] T. J. R. Hughes, J. A. Cottrell, and Y. Bazilevs. “Isogeometric analy- si...

work page Pith review arXiv 2018
[347]

[Zie71] O

doi: 10.1007/s00211-007-0063-4 . [Zie71] O. C. Zienkiewicz. The finite element method in engineering science . The second, expanded and revised, edition of The finite element method in structural and continuum mechanics . McGraw-Hill, London-New York- D¨ usseldorf, 1971, pp. xiv+521. (L. Diening) Department of Mathematics, Bielefeld University, Postfach 1...

work page doi:10.1007/s00211-007-0063-4
[443]

Collision in a cross-shaped domain—a steady 2d Navier-Stokes example demonstrating the importance of mass conservation in CFD

doi: 10.1515/cmam-2020-0023 . [Lin09] A. Linke. “Collision in a cross-shaped domain—a steady 2d Navier-Stokes example demonstrating the importance of mass conservation in CFD”. REFERENCES 35 Comput. Methods Appl. Mech. Engrg. 198.41-44 (2009), pp. 3278–3286. doi: 10.1016/j.cma.2009.06.016. [LLMS17] P. L. Lederer, A. Linke, C. Merdon, and J. Sch¨ oberl. “D...

work page doi:10.1515/cmam-2020-0023 2009
[1508]

Conforming and divergence-free Stokes ele- ments on general triangular meshes

doi: 10.1093/imanum/drt053. [GN14b] J. Guzm´ an and M. Neilan. “Conforming and divergence-free Stokes ele- ments on general triangular meshes”. Math. Comp. 83.285 (2014), pp. 15–

doi:10.1093/imanum/drt053 2014

Pith tools

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

[1] [36]

Application of Randomized Quadrature Formulas to the Finite Element Method for Elliptic Equations

doi: 10.1090/S0025-5718-2013-02753-6 . [GN18] J. Guzm´ an and M. Neilan. “inf-sup stable finite elements on barycentric refinements producing divergence-free approximations in arbitrary dimen- sions”. SIAM J. Numer. Anal. 56.5 (2018), pp. 2826–2844. doi: 10.1137/ 17M1153467. [HCB05] T. J. R. Hughes, J. A. Cottrell, and Y. Bazilevs. “Isogeometric analy- si...

work page Pith review arXiv 2018

[2] [347]

[Zie71] O

doi: 10.1007/s00211-007-0063-4 . [Zie71] O. C. Zienkiewicz. The finite element method in engineering science . The second, expanded and revised, edition of The finite element method in structural and continuum mechanics . McGraw-Hill, London-New York- D¨ usseldorf, 1971, pp. xiv+521. (L. Diening) Department of Mathematics, Bielefeld University, Postfach 1...

work page doi:10.1007/s00211-007-0063-4

[3] [443]

Collision in a cross-shaped domain—a steady 2d Navier-Stokes example demonstrating the importance of mass conservation in CFD

doi: 10.1515/cmam-2020-0023 . [Lin09] A. Linke. “Collision in a cross-shaped domain—a steady 2d Navier-Stokes example demonstrating the importance of mass conservation in CFD”. REFERENCES 35 Comput. Methods Appl. Mech. Engrg. 198.41-44 (2009), pp. 3278–3286. doi: 10.1016/j.cma.2009.06.016. [LLMS17] P. L. Lederer, A. Linke, C. Merdon, and J. Sch¨ oberl. “D...

work page doi:10.1515/cmam-2020-0023 2009

[4] [1508]

Conforming and divergence-free Stokes ele- ments on general triangular meshes

doi: 10.1093/imanum/drt053. [GN14b] J. Guzm´ an and M. Neilan. “Conforming and divergence-free Stokes ele- ments on general triangular meshes”. Math. Comp. 83.285 (2014), pp. 15–

doi:10.1093/imanum/drt053 2014