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

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.