Triangle-Based Spectral Element Method

• Triangle-based SEM is a high-order numerical method that employs polynomial spectral expansions on unstructured triangular meshes for solving complex PDEs efficiently.
• It leverages modal sparsity, advanced basis constructions, and hierarchical operator coupling to achieve scalable algorithms with spectral accuracy and reduced computational cost.
• The method offers adaptive hp-refinement and almost-banded solvers, ensuring robust stability, energy/entropy stability, and optimal error convergence for diverse applications.

The triangle-based spectral element method (SEM) encompasses a suite of high-order numerical techniques employing polynomial spectral expansions and direct-solver frameworks on unstructured triangular meshes. These methods exploit the geometric flexibility of triangles to address partial differential equations (PDEs), integral equations, and time-dependent problems on complex domains. Central to their practical competitiveness is modal sparsity, advanced basis construction, efficient operator coupling, and adaptive hp-refinement schemes. The current state-of-the-art incorporates ultraspherical spectral methods, hierarchical Poincaré–Steklov coupling, tensor-product summation-by-parts (SBP) discretizations, adaptive quadrature, and barycentric tensor evaluation, leading to scalable algorithms with spectral accuracy, robust stability, and demonstrably optimal computational complexity (Fortunato et al., 2020).

1. Triangular Spectral Discretization and Polynomial Bases

Spectral discretization on triangles relies on polynomial bases constructed in reference coordinates via collapsed mappings (Duffy, rectangle-triangle, square-triangle) that facilitate analytic integration, stable numerical evaluation, and sparse operator structure. The leading bases include Chebyshev polynomials (ultraspherical methods), Jacobi polynomials (Legendre, Koornwinder), and Dubiner polynomials in collapsed coordinates.

To realize a modal expansion, consider the reference right triangle $T_\text{ref}$ with vertices $(0,0), (1,0), (0,1)$. Using the Duffy map $$(r,s) \mapsto (x,y) = (\frac14(1+r)(1-s), \frac12(1-s))$$ (Fortunato et al., 2020), functions $u(r,s)$ are expanded as

$$u(r,s) \approx \sum_{i=0}^p \sum_{j=0}^p X_{ij} T_i(r) T_j(s)$$

where $T_k$ is the Chebyshev polynomial of degree $k$. Differentiation is efficiently realized using ultraspherical recurrences and inter-basis conversion operators $S_\lambda$, ensuring sparse representation in higher-order derivatives.

Alternative bases, such as the Proriol-Koornwinder-Dubiner (PKD) modal family, exploit orthonormality and tensor-product structure:

$$\phi^{(\alpha)}(\xi) = \psi_1^{(\alpha_1)}(\eta_1) \psi_2^{(\alpha_1,\alpha_2)}(\eta_2)$$

with $\psi_1^{(\alpha_1)}$ and $\psi_2^{(\alpha_1,\alpha_2)}$ constructed from Jacobi polynomials in collapsed coordinates (Montoya et al., 2023, Montoya et al., 2023). These bases enable exact representation of $L^2$ and $H^1$ projection and interpolation error estimates:

$$\|u - \Pi_N u\|_{L^2(T)} \leq C N^{-r} \|u\|_{H^r(T)}$$

for analytic $u$, spectral convergence in $p$ follows (Samson et al., 2012).

2. Sparse Operator Construction and Almost-Banded Systems

The ultraspherical and Jacobi-based approaches yield sparse differentiation, multiplication, and conversion matrices. Local elliptic systems $L \vec{X} = \vec{F}$ are assembled as block-banded with $O(p)$ bandwidth independent of $p$; Dirichlet conditions lead to additional $O(p)$ dense rows. For $p \gtrsim 5$, the interior block can be efficiently inverted via banded solvers and Woodbury updates in $O(p^4)$ time (Fortunato et al., 2020). In integral equations, operator algebra exploiting Clenshaw's algorithm and three-term Jacobi recurrences produces Galerkin matrices with $\mathcal{O}(N^2)$ nonzeros and exponential convergence for analytic data (Gutleb et al., 2019, Olver et al., 2019).

3. Hierarchical Poincaré–Steklov Coupling and Direct Solvers

Global solution strategies employ hierarchical binary-tree merging of local solution and Dirichlet-to-Neumann (Steklov) operators:

• Each triangle $E$ stores a solution operator $S_E$ and Steklov operator $T_E = D_E S_E$.
• Pairwise merging of triangles along shared edges via a Schur complement produces global interface solution operators and recursive Steklov maps.
• Full mesh assembly through $O(\log(1/h))$ levels yields overall complexity $O(p^4/h^3)$ for build and solve on $O(1/h^2)$ triangles (Fortunato et al., 2020).

Efficient direct solvers and adaptivity are enabled by this reuse of precomputed operators, facilitating fast elliptic solves, implicit and semi-implicit time-steppers, and dynamic updates to the right-hand side in $O(N p^2)$ time (as in the ultraSEM MATLAB package).

4. hp-Adaptivity, Stability, and Error Analysis

Modal sparsity ensures that increasing polynomial degree $p$ does not induce substantial numerical overhead, and local refinement in mesh size $h$ captures singularities. For analytic solutions, maximum-norm error bounds

$\|u - u_p\|_{\infty} \leq C \rho^{-p}$

demonstrate spectral accuracy; $L^2$-errors satisfy

$\|u - u_{hp}\|_{H^k} = O(h^{p+1-k})$

for homogeneous shape-regular meshes. Stability is guaranteed by backward-stable ultraspherical QR factorizations and bounded condition of inter-basis conversion operators (Fortunato et al., 2020). With adaptive $hp$-refinement, super-algebraic convergence in total degrees of freedom $N$ is achievable.

5. Tensor-Product SBP, Entropy Stability, and Efficient Evaluation

Discontinuous spectral-element methodologies for hyperbolic PDEs employ tensor-product SBP operators in collapsed coordinates, enforcing the discrete summation-by-parts property. Chain rule for differentiation and weight-adjusted mass matrix inversion produce robust, energy-stable discretizations:

• Tensor-product structure reduces per-element complexity from $O(p^4)$ to $O(p^3)$ by limiting flux evaluations to $O(p^3)$ node pairs (Montoya et al., 2023, Montoya et al., 2023).
• Modal PKD basis supports sum-factorization for fast matrix-vector computations.
• Entropy variables and flux-differencing guarantee entropy stability and robustness, matching multidimensional SBP schemes in accuracy with reduced computational cost.

High-performance point evaluation and post-processing are further expedited by per-element barycentric tensor kernels, yielding $\mathcal{O}(N)$ evaluation cost for function values and derivatives at arbitrary points (Laughton et al., 2021).

6. Implementation Frameworks and Practical Algorithms

Triangle-based SEMs are available in open-source systems (e.g. ultraSEM for MATLAB), supporting:

• Automated construction of reference triangle domains via Duffy/affine mappings,
• Symbolic evaluation of Jacobian and derivative metrics,
• Fast recurrence-based assembly of modal operators,
• Clenshaw–Curtis quadrature and inter-basis conversion for edge transforms,
• Adaptive hierarchical merging for direct solution,
• Efficient update routines for inhomogeneous right-hand sides (Fortunato et al., 2020).

Iterative explicit time-stepping for continuous SEMs on unstructured triangles leverages lower-triangular pseudo-mass matrices, locally invertible, for efficient high-order projection and Runge–Kutta schemes (Appleton et al., 2019). Stabilization (SUPG, CIP, OSS) and time discretization (SSPRK, deferred correction, IMEX-BDF) enable robust solution of hyperbolic and evolving-surface problems (Michel et al., 2022, Zavalani, 30 Dec 2025).

7. Extensions, Current Directions, and Comparative Properties

Triangle-based SEMs extend naturally to tensor-product SBP operators on curved geometries, higher-dimensional domains, and evolving triangulated surfaces. Merging strategies and Dubiner bases transfer to globally coupled methods for reaction-diffusion equations, surface evolution, and complex interface problems (Zavalani, 30 Dec 2025). Comparative studies show that, for $p\ge5$, tensor-product SBP schemes on triangles outperform multidimensional SBP in computational efficiency without loss in accuracy or robustness (Montoya et al., 2023, Montoya et al., 2023).

Key distinguishing features:

• Modal sparsity and almost-banded structure,
• hp-adaptivity with no barrier to raising $p$,
• Energy/entropy-stable discretization,
• Spectral convergence for analytic solutions,
• Fast hierarchical direct solvers,
• Element-wise operator algebra suitable for unstructured meshes.

Triangle-based spectral element methods are thus foundational to high-order simulation technology on arbitrary domains, supporting scientific, engineering, and applied mathematics research (Fortunato et al., 2020).