pith. machine review for the scientific record. sign in

arxiv: 2604.03097 · v1 · submitted 2026-04-03 · 🧮 math.NA · cs.NA

Recognition: 2 theorem links

· Lean Theorem

A High-Order Fast Direct Solver for Surface PDEs on Triangles

Gentian Zavalani
Authors on Pith no claims yet

Pith reviewed 2026-05-13 19:09 UTC · model grok-4.3

classification 🧮 math.NA cs.NA
keywords hierarchical Poincaré-Steklovtriangular meshesfast direct solversurface PDEsDubiner polynomialsspectral accuracyO(N log N) complexityunstructured meshes
0
0 comments X

The pith

A triangular hierarchical Poincaré-Steklov method produces fast direct solvers for elliptic surface PDEs on unstructured meshes with O(N log N) complexity and spectral accuracy.

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

The paper introduces a triangular version of the hierarchical Poincaré-Steklov method, called THPS, that works with Dubiner polynomial bases on triangles instead of tensor-product bases on quadrilaterals. Local solution operators and Dirichlet-to-Neumann maps are built on each triangle and then merged up a hierarchy of subdomains. This merging step produces a direct solver whose setup cost and per-solve cost both scale as O(N log N) for N triangular elements. Because the operators are precomputed once and reused, the scheme is especially efficient for repeated solves such as those arising in implicit time-stepping of time-dependent surface PDEs. Numerical tests confirm that the method retains the spectral accuracy of the original HPS approach even on unstructured triangular meshes of complex geometries.

Core claim

By constructing local solution operators and Dirichlet-to-Neumann maps on triangles using orthogonal Dubiner polynomial bases and merging these operators hierarchically, the THPS scheme yields a fast direct solver for elliptic PDEs on surfaces that requires O(N log N) work per solve on a mesh of N elements while preserving spectral accuracy.

What carries the argument

The THPS scheme, which builds and hierarchically merges local solution operators and Dirichlet-to-Neumann maps using Dubiner polynomial bases on triangular elements.

If this is right

  • The precomputed operators can be reused across many time steps in implicit integrators for time-dependent surface PDEs.
  • The method supports high-order convergence on meshes that conform to complex surface geometries without requiring tensor-product structure.
  • Setup cost is amortized over repeated solves, making the approach practical for problems requiring many right-hand sides.
  • Accuracy remains spectral even when the underlying triangulation is unstructured.

Where Pith is reading between the lines

These are editorial extensions of the paper, not claims the author makes directly.

  • The same hierarchical merging strategy could be tested on other element shapes such as tetrahedra if analogous orthogonal bases are available.
  • Because the method separates geometry from the solver, it may allow existing high-order surface codes to adopt fast direct solves without changing their mesh generators.
  • For very large N the logarithmic factor may become noticeable in practice, suggesting a possible need for parallel merging algorithms.

Load-bearing premise

The hierarchical merging of local operators constructed on triangles with Dubiner bases preserves both the O(N log N) complexity and the spectral accuracy of the classical quadrilateral HPS method.

What would settle it

Running the solver on a sequence of successively refined unstructured triangular meshes for a smooth elliptic test problem and observing either superlinear growth in runtime beyond O(N log N) or loss of spectral convergence order.

Figures

Figures reproduced from arXiv: 2604.03097 by Gentian Zavalani.

Figure 1
Figure 1. Figure 1: Construction of a surface parametrization over the reference simplex [PITH_FULL_IMAGE:figures/full_fig_p006_1.png] view at source ↗
Figure 2
Figure 2. Figure 2: Interface coupling of two surface elements via continuity of the solution and its [PITH_FULL_IMAGE:figures/full_fig_p010_2.png] view at source ↗
Figure 3
Figure 3. Figure 3: A high-order triangular patch mesh is used to discretize the Laplace–Beltrami [PITH_FULL_IMAGE:figures/full_fig_p014_3.png] view at source ↗
Figure 4
Figure 4. Figure 4: A high-order triangulated sphere mesh is used to approximate a spherical [PITH_FULL_IMAGE:figures/full_fig_p015_4.png] view at source ↗
Figure 5
Figure 5. Figure 5: (Left) Computed solution at t = 1 using time step ∆t = 10−3 . (Right) L∞-error versus polynomial degree for different mesh sizes h. The simulations were performed using the implicit–explicit backward differentiation formula (IMEX–BDF4) scheme. Numerical study of spatial pattern formation in Turing systems. Having established the method’s robustness for this simpler case, we now turn to a more complex and b… view at source ↗
Figure 6
Figure 6. Figure 6: Turing model solution u1 on a Swiss cheese surface. Top: r1 = 3.5, r2 = 0 at t = 0,200,600. Bottom: r1 = 0.02, r2 = 0.2 at t = 0,20,200. Simulated with IMEX-BDF4. Next, we simulate the Turing model on two representative surfaces to investigate the influence of geometry on pattern formation. The first geometry is an asymmetric torus,5 followed by a triangulated Stanford Bunny surface. In both cases, identic… view at source ↗
Figure 7
Figure 7. Figure 7: Turing model solution u1 on the asymmetric torus and Stanford Bunny at times t = 0, 20, and 200 using IMEX-BDF4 with r1 = 0.02, r2 = 0.2. In [PITH_FULL_IMAGE:figures/full_fig_p018_7.png] view at source ↗
Figure 8
Figure 8. Figure 8: Pattern formation process of the reaction–diffusion model on the surface of [PITH_FULL_IMAGE:figures/full_fig_p019_8.png] view at source ↗
Figure 9
Figure 9. Figure 9: Patterns obtained from the coupled system [PITH_FULL_IMAGE:figures/full_fig_p020_9.png] view at source ↗
Figure 10
Figure 10. Figure 10: Patterns obtained from the coupled system [PITH_FULL_IMAGE:figures/full_fig_p020_10.png] view at source ↗
read the original abstract

We develop a triangular formulation of the hierarchical Poincar\'e-Steklov (HPS) method for elliptic partial differential equations on surfaces, allowing high-order discretizations on unstructured meshes and complex geometries. Classical HPS formulations rely on high-order quadrilateral meshes and tensor-product spectral discretizations, which enable efficient algorithms but restrict applicability to structured geometries. To overcome this restriction, we introduce a triangle-based hierarchical Poincar\'e-Steklov scheme (THPS) built on orthogonal Dubiner polynomial bases. As in the classical HPS framework, local solution operators and Dirichlet-to-Neumann maps are constructed and merged hierarchically, yielding a fast direct solver with $O(N \log N)$ complexity for repeated solves on meshes with $N$ elements. The reuse of precomputed operators makes the method particularly effective for implicit time-stepping of surface PDEs. Numerical experiments demonstrate that the proposed method retains spectral accuracy and achieves high-order convergence for a range of static and time-dependent test problems.

Editorial analysis

A structured set of objections, weighed in public.

Desk editor's note, referee report, simulated authors' rebuttal, and a circularity audit. Tearing a paper down is the easy half of reading it; the pith above is the substance, this is the friction.

Referee Report

2 major / 2 minor

Summary. The manuscript develops a triangular hierarchical Poincaré-Steklov (THPS) method for high-order discretization of elliptic surface PDEs on unstructured triangular meshes. Local solution operators and Dirichlet-to-Neumann maps are built using orthogonal Dubiner polynomial bases and merged hierarchically to produce a fast direct solver claimed to have O(N log N) complexity for repeated solves on meshes with N elements, while preserving spectral accuracy. The approach is positioned as an extension of classical quadrilateral HPS methods to complex geometries, with demonstrated utility for implicit time-stepping of time-dependent surface PDEs.

Significance. If the O(N log N) complexity and spectral accuracy claims hold under the non-tensor-product Dubiner discretization, the work would meaningfully extend fast direct solvers to unstructured triangular surface meshes, enabling efficient high-order solutions on complex geometries where quadrilateral HPS methods are inapplicable. The reuse of precomputed operators for time-dependent problems is a practical strength.

major comments (2)
  1. [§4.1–4.3] §4.1–4.3: The hierarchical merging procedure for dense DtN maps constructed from Dubiner bases (Eqs. (8)–(11)) is described, but no explicit operation-count analysis is given showing that the per-merge cost remains O(1) with respect to polynomial degree p. Without tensor-product separability, the Schur-complement updates appear to be dense matrix operations whose cost scales with the number of edge degrees of freedom (~p), raising the possibility that total build cost is O(N p^3 log N) rather than the claimed O(N log N) independent of p.
  2. [§5.2] §5.2, Table 2: Reported wall-clock timings for the build phase on successively refined meshes show good scaling with N at fixed p=4, but no data or analysis is provided for increasing p at fixed N, leaving the independence of complexity from polynomial degree unverified.
minor comments (2)
  1. [§3.1] The notation for the reference-element Dubiner basis and the precise definition of the edge restriction operators could be clarified with an explicit diagram or table of degrees of freedom.
  2. [Abstract, §4.3] A few typographical inconsistencies appear in the complexity statements between the abstract and §4.3 (e.g., “O(N log N)” vs. “O(N log N) per solve after O(N log N) setup”).

Simulated Author's Rebuttal

2 responses · 0 unresolved

We thank the referee for the careful reading and constructive feedback on our manuscript. The comments highlight the need for a more explicit complexity analysis and additional numerical verification of the scaling with polynomial degree p. We will revise the manuscript accordingly to address these points while preserving the core claims, which apply to the solve phase for fixed p.

read point-by-point responses

  1. Referee: [§4.1–4.3] §4.1–4.3: The hierarchical merging procedure for dense DtN maps constructed from Dubiner bases (Eqs. (8)–(11)) is described, but no explicit operation-count analysis is given showing that the per-merge cost remains O(1) with respect to polynomial degree p. Without tensor-product separability, the Schur-complement updates appear to be dense matrix operations whose cost scales with the number of edge degrees of freedom (~p), raising the possibility that total build cost is O(N p^3 log N) rather than the claimed O(N log N) independent of p.

    Authors: We agree that an explicit operation-count analysis would improve clarity. The manuscript's O(N log N) claim refers to the solve phase (application of the precomputed hierarchical operators) for repeated solves with fixed p, as is standard when p is treated as a constant. For the build phase, each merge performs dense Schur-complement operations on O(p)-sized blocks, incurring O(p^3) cost per merge. Across O(log N) levels and O(N) merges, the build cost is O(N p^3 log N). We will add a dedicated paragraph in Section 4 detailing this count, distinguishing build from solve phases, and noting that the solve complexity remains O(N log N) independent of p only when p is fixed. revision: yes

  2. Referee: [§5.2] §5.2, Table 2: Reported wall-clock timings for the build phase on successively refined meshes show good scaling with N at fixed p=4, but no data or analysis is provided for increasing p at fixed N, leaving the independence of complexity from polynomial degree unverified.

    Authors: We concur that timings varying p at fixed N would provide useful verification. We will augment Section 5.2 with new experiments and a supplementary table reporting build and solve wall-clock times for p = 2,4,6,8 on a fixed mesh (N ≈ 10^4). These will illustrate the expected O(p^3 log N) build scaling and O(log N) solve scaling (p fixed), confirming the analysis to be added in Section 4. revision: yes

Circularity Check

0 steps flagged

No circularity: THPS is a new construction with independent complexity claim

full rationale

The paper introduces THPS as an explicit new scheme on triangles using Dubiner bases, with local operators and hierarchical merging defined directly in the text rather than reduced to prior HPS results by definition or self-citation. The O(N log N) claim is asserted for the merged operators without any equation or step that equates the output complexity to fitted inputs or renames a known result; the derivation chain remains self-contained as a construction on unstructured meshes. No load-bearing step collapses to a self-citation chain or ansatz smuggled from the authors' prior work.

Axiom & Free-Parameter Ledger

0 free parameters · 2 axioms · 0 invented entities

The central claim rests on standard properties of Dubiner polynomials and the hierarchical operator-merging framework; no free parameters or new entities are introduced in the abstract.

axioms (2)
  • standard math Orthogonal Dubiner polynomial bases provide high-order approximation on triangles
    Invoked for the discretization on triangular elements.
  • domain assumption Hierarchical merging of local operators yields O(N log N) complexity
    Assumed to carry over from classical HPS to the triangular case.

pith-pipeline@v0.9.0 · 5464 in / 1169 out tokens · 35534 ms · 2026-05-13T19:09:02.899307+00:00 · methodology

discussion (0)

Sign in with ORCID, Apple, or X to comment. Anyone can read and Pith papers without signing in.

Lean theorems connected to this paper

Citations machine-checked in the Pith Canon. Every link opens the source theorem in the public Lean library.

What do these tags mean?
matches
The paper's claim is directly supported by a theorem in the formal canon.
supports
The theorem supports part of the paper's argument, but the paper may add assumptions or extra steps.
extends
The paper goes beyond the formal theorem; the theorem is a base layer rather than the whole result.
uses
The paper appears to rely on the theorem as machinery.
contradicts
The paper's claim conflicts with a theorem or certificate in the canon.
unclear
Pith found a possible connection, but the passage is too broad, indirect, or ambiguous to say the theorem truly supports the claim.

Reference graph

Works this paper leans on

36 extracted references · 36 canonical work pages · 1 internal anchor

  1. [1]

    A High-Order Fast Direct Solver for Surface PDEs on Triangles

    Introduction.Partial differential equations (PDEs) on surfaces arise in a wide range of applications, including fluid dynamics, materials science, and biological pattern formation. Their numerical solution requires discretization techniques that accurately capture both the underlying geometry and the associated surface differ- ential operators. Standard a...

    work page internal anchor Pith review Pith/arXiv arXiv 2026

  2. [2]

    Dubiner polynomials on the reference triangle.As in quadrilateral based spectral element [9], triangular spectral element discretizations require an orthogonal polynomial basis together with a suitable set of interpolation nodes. On quadrilateral reference elements, such as the square[−1,1] 2, these ingredients arise naturally from tensor-product construc...

    work page

  3. [3]

    IfΓis not a closed surface,(3.1)may also be subject to boundary conditions, e.g.,u(x) =h(x)forx∈∂Γand some functionh

    High-order parametric surface approximation.We consider a general elliptic surface PDE onΓ, (3.1)L Γu(x) =f(x),x∈Γ, where f(x) is a smooth function onΓ and LΓ is a variable-coefficient linear second- order elliptic surface operator. IfΓis not a closed surface,(3.1)may also be subject to boundary conditions, e.g.,u(x) =h(x)forx∈∂Γand some functionh. HPS ON...

    work page

  4. [4]

    The discretization is derived directly from the strong formulation of the problem

    High-order spectral collocation on a single element.We now describe a spectral collocation method for discretizing(3.1) on a single surface elementEk ⊂Γ. The discretization is derived directly from the strong formulation of the problem. Functions defined onEk ⊂Γ are represented at the same interpolation nodes used for the geometry. Let u be a function def...

    work page

  5. [5]

    DtNb1b1 E1 0 0 DtN b2b2 E2 # +

    Domain decomposition methods.The spectral collocation method described in the previous section converges very quickly as the number of pointsn increases provided the solutionuis smooth. The matrixLEk that arises from the discretization has some structure and contains many zeros, but it is still considerably denser than the matrices produced by finite diff...

    work page

  6. [6]

    Computational asymptotic complexity.For an order-n discretization using triangular spectral elements, the total number of degrees of freedom associated with the mesh{Ek}N k=1 scales asN (n+1)(n+2) 2 .As in the quadrilateral based HPS method [9, 21],thecomputationalcomplexityconsistsofthreemainstagesandisgovernedbythe number of local degrees of freedom. On...

    work page

  7. [7]

    Time-dependentequations.Whilethetriangle-basedHPSschemeisdesigned to solve stationary problems modeled by linear elliptic PDEs, it can also be useful for accelerating time-dependent problems modeled by parabolic PDEs. We start by considering the reaction-diffusion systems, which are widely regarded as key mecha- nisms for pattern formation in a variety of...

    work page

  8. [8]

    activator

    Numerical examples.In this section, we assess the accuracy, convergence, and computational performance of the proposed triangular HPS method on a set of representative surface PDE problems. We begin with the Laplace–Beltrami equation on the sphere, where analytical solutions are available, allowing for a detailed study of spatial convergence. We then cons...

    work page 2000

  9. [9]

    +v 2(1−r 2v1) +q 1u1 +q 2u1v2 +q 3u1v2 2 ∂v 2 ∂t =δ v2 ∆Γv2 +β ′v2 1+ α ′r1 β ′ v1v2 +v 1(γ ′ +r 2v2)−q 2u2v1 −q 3u2 2v1 ∂u 1 ∂t =δ u1∆Γu1 +αu 1(1−r 1v2

    work page

  10. [10]

    Thissetupallowsustoinvestigatehowpatternpropertieschangeduetotheinteraction betweenthetwocoupledsystems

    +u 2(1−r 2u1) ∂u 2 ∂t =δ u2∆Γu2 +βu 2 1+ αr 1 β u1u2 +u 1(γ+r 2u2). Thissetupallowsustoinvestigatehowpatternpropertieschangeduetotheinteraction betweenthetwocoupledsystems. Forthesimulations,weusethefollowingparameter values: •For the(v 1,v 2)system: α ′ =0.398,β ′ =−0.41,γ ′ =−α ′,δ v2 =5×10 −3,δ v1 =0.122δ v2. •For the(u 1,u 2)system: α=0.899,β=−0.91,γ=...

    work page

  11. [11]

    Beyond uni-directional hierarchical wrinkle formation

    Appendix. Tangential differential calculus on surfaces.We review the basic concepts of differential geometry relevant to our setting. We start by describing the parametric representation, then introduce the first fundamental formg, the tangential operators (gradient∇Γ, divergencedivΓ, and Laplace–Beltrami∆Γ). LetΓbe a smooth surface embedded inRd+1 locall...

    work page

  12. [12]

    U. M. Ascher, S. J. Ruuth, and B. T. Wetton,Implicit-explicit methods for time- dependent partial differential equations, SIAM Journal on Numerical Analysis, 32 (1995), pp. 797–823

    work page 1995

  13. [13]

    Barrio, C

    R. Barrio, C. V area, J. Aragón, and P . Maini,A two-dimensional numerical study of spatial pattern formation in interacting Turing systems, Bulletin of mathematical biology, 61 (1999), pp. 483–505

    work page 1999

  14. [14]

    W. L. Briggs, V . E. Henson, and S. F. McCormick,Amultigridtutorial,SIAM,2000

    work page 2000

  15. [15]

    Bunow , J.-P

    B. Bunow , J.-P . Kernevez, G. Joly , and D. Thomas,Pattern formation by reaction- diffusion instabilities: Application to morphogenesis in Drosophila, Journal of theoreti- cal biology, 84 (1980), pp. 629–649

    work page 1980

  16. [16]

    T. A. Da vis,Algorithm 832:UMFP ACKv4.3—an unsymmetric-pattern multifrontal method,ACMTransactionsonMathematicalSoftware(TOMS),30(2004),pp.196– 199

    work page 2004

  17. [17]

    Dubiner,Spectral methods on triangles and other domains, Journal of Scientific Computing, 6 (1991), pp

    M. Dubiner,Spectral methods on triangles and other domains, Journal of Scientific Computing, 6 (1991), pp. 345–390

    work page 1991

  18. [18]

    Dziuk and C

    G. Dziuk and C. M. Elliott,FiniteelementmethodsforsurfacePDEs,ActaNumerica, 22 (2013), pp. 289–396

    work page 2013

  19. [19]

    Fornberg,A Practical Guide to Pseudospectral Methods, Cambridge University Press, 1998

    B. Fornberg,A Practical Guide to Pseudospectral Methods, Cambridge University Press, 1998

    work page 1998

  20. [20]

    Fortunato,A high-order fast direct solver for surface pdes, SIAM Journal on Scientific Computing, 46 (2024), pp

    D. Fortunato,A high-order fast direct solver for surface pdes, SIAM Journal on Scientific Computing, 46 (2024), pp. A2582–A2606

    work page 2024

  21. [21]

    Fortunato and A

    D. Fortunato and A. Townsend,Fast poisson solvers for spectral methods, IMA Journal of Numerical Analysis, 40 (2020), pp. 1994–2018

    work page 2020

  22. [22]

    Gillman and P .-G

    A. Gillman and P .-G. Martinsson,A direct solver withO(n)complexity for variable coefficient elliptic pdes discretized via a high-order composite spectral collocation method, SIAM Journal on Scientific Computing, 36 (2014), pp. A2023–A2046

    work page 2014

  23. [23]

    Isaac,Recursive, parameter-free, explicitly defined interpolation nodes for simplices, SIAM Journal on Scientific Computing, 42 (2020), pp

    T. Isaac,Recursive, parameter-free, explicitly defined interpolation nodes for simplices, SIAM Journal on Scientific Computing, 42 (2020), pp. A4046–A4062

    work page 2020

  24. [24]

    Jeong, Y

    D. Jeong, Y . Li, Y . Choi, M. Yoo, D. Kang, J. P ark, J. Choi, and J. Kim,Numerical simulation of the zebra pattern formation on a three-dimensional model, Physica A: Statistical Mechanics and its Applications, 475 (2017), pp. 106–116

    work page 2017

  25. [25]

    Jonsson, M

    T. Jonsson, M. G. Larson, and K. Larsson,Cutfiniteelementmethodsforellipticprob- lems on multipatch parametric surfaces, Computer Methods in Applied Mechanics and Engineering, 324 (2017), pp. 366–394

    work page 2017

  26. [26]

    Karniadakis and S

    G. Karniadakis and S. Sherwin,Spectral/hp element methods for computational fluid dynamics, American Chemical Society, 2013

    work page 2013

  27. [27]

    Koornwinder,Two-variable analogues of the classical orthogonal polynomials, in Theory and application of special functions, Elsevier, 1975, pp

    T. Koornwinder,Two-variable analogues of the classical orthogonal polynomials, in Theory and application of special functions, Elsevier, 1975, pp. 435–495

    work page 1975

  28. [28]

    R. J. LeVeque,Finite difference methods for ordinary and partial differential equations: steady-state and time-dependent problems, SIAM, 2007

    work page 2007

  29. [29]

    R. W. Lewis, P . Nithiarasu, and K. N. Seetharamu,Fundamentals of the finite element method for heat and fluid flow, John Wiley & Sons, 2004

    work page 2004

  30. [30]

    P . K. Maini, K. J. P ainter, and H. N. P . Chau,Spatial pattern formation in chemical and biological systems, Journal of the Chemical Society, Faraday Transactions, 93 (1997), pp. 3601–3610

    work page 1997

  31. [31]

    Martinsson,A direct solver for variable coefficient elliptic pdes discretized via a compositespectralcollocationmethod,JournalofComputationalPhysics,242(2013), pp

    P .-G. Martinsson,A direct solver for variable coefficient elliptic pdes discretized via a compositespectralcollocationmethod,JournalofComputationalPhysics,242(2013), pp. 460–479. HPS ON TRIANGULAR SURFACES23 [21]P .-G. Martinsson,Fast direct solvers for elliptic PDEs, SIAM, 2019

    work page 2013

  32. [32]

    T. P . A. Mathew,Domain decomposition methods for the numerical solution of partial differential equations, Springer, 2008

    work page 2008

  33. [33]

    Proriol,Surunefamilledepolynomesádeuxvariablesorthogonauxdansuntriangle, Comptes rendus hebdomadaires des seances de l academie des sciences, 245 (1957), pp

    J. Proriol,Surunefamilledepolynomesádeuxvariablesorthogonauxdansuntriangle, Comptes rendus hebdomadaires des seances de l academie des sciences, 245 (1957), pp. 2459–2461

    work page 1957

  34. [34]

    Quarteroni and A

    A. Quarteroni and A. V alli,Numericalapproximationofpartialdifferentialequations, vol. 23, Springer Science & Business Media, 2008

    work page 2008

  35. [35]

    J. C. Strikwerda,Finite difference schemes and partial differential equations, SIAM, 2004

    work page 2004

  36. [36]

    A. M. Turing,Thechemicalbasisofmorphogenesis,Bulletinofmathematicalbiology, 52 (1990), pp. 153–197. [27]R. S. V arga,Matrix iterative methods, Prentice Hall Incorporated, 1962

    work page 1990