We introduce a sparse and very high order hp-finite element method for the weak form of the Helmholtz equation. The domain may be a disk, an annulus, or a cylinder. The cells of the mesh are an innermost disk (omitted if the domain is an annulus) and concentric annuli.
We demonstrate the effectiveness of this method on PDEs with radial direction discontinuities in the coefficients and data. The discretization matrix is always symmetric and positive-definite in the positive-definite Helmholtz regime. Moreover, the Fourier modes decouple, reducing a two-dimensional PDE solve to a series of one-dimensional solves that may be computed in parallel, scaling with linear complexity. In the positive-definite case, we utilize the ADI method of Fortunato and Townsend to apply the method to a 3D cylinder with a quasi-optimal complexity solve.