A high-order fast direct solver for surface PDEs

We introduce a fast direct solver for variable-coefficient elliptic partial differential equations on surfaces based on the hierarchical Poincar\'e-Steklov method. The method takes as input an unstructured, high-order quadrilateral mesh of a surface and discretizes surface differential operators on each element using a high-order spectral collocation scheme. Elemental solution operators and Dirichlet-to-Neumann maps tangent to the surface are precomputed and merged in a pairwise fashion to yield a hierarchy of solution operators that may be applied in $\mathcal{O}(N \log N)$ operations for a mesh with $N$ elements. The resulting fast direct solver may be used to accelerate high-order implicit time-stepping schemes, as the precomputed operators can be reused for fast elliptic solves on surfaces. On a standard laptop, precomputation for a 12th-order surface mesh with over 1 million degrees of freedom takes 17 seconds, while subsequent solves take only 0.25 seconds. We apply the method to a range of problems on both smooth surfaces and surfaces with sharp corners and edges, including the static Laplace-Beltrami problem, the Hodge decomposition of a tangential vector field, and some time-dependent nonlinear reaction-diffusion systems.