High-order quadrature on multi-component domains implicitly defined by multivariate polynomials

A high-order quadrature algorithm is presented for computing integrals over curved surfaces and volumes whose geometry is implicitly defined by the level sets of (one or more) multivariate polynomials. The algorithm recasts the implicitly defined geometry as the graph of an implicitly defined, multi-valued height function, and applies a dimension reduction approach needing only one-dimensional quadrature. In particular, we explore the use of Gauss-Legendre and tanh-sinh methods and demonstrate that the quadrature algorithm inherits their high-order convergence rates. Under the action of h-refinement with q fixed, the quadrature schemes yield an order of accuracy of 2q, where q is the one-dimensional node count; numerical experiments demonstrate up to 22nd order. Under the action of q-refinement with the geometry fixed, the convergence is approximately exponential, i.e., doubling q approximately doubles the number of accurate digits of the computed integral. Complex geometry is automatically handled by the algorithm, including, e.g., multi-component domains, tunnels, and junctions arising from multiple polynomial level sets, as well as self-intersections, cusps, and other kinds of singularities. A variety of numerical experiments demonstrates the quadrature algorithm on two- and three-dimensional problems, including: randomly generated geometry involving multiple high-curvature pieces; challenging examples involving high degree singularities such as cusps; adaptation to simplex constraint cells in addition to hyperrectangular constraint cells; and boolean operations to compute integrals on overlapping domains.