  This article presents a high-order accurate numerical method for the
evaluation of singular volume integral operators, with attention focused on
operators associated with the Poisson and Helmholtz equations in two
dimensions. Following the ideas of the density interpolation method for
boundary integral operators, the proposed methodology leverages Green's third
identity and a local polynomial interpolant of the density function to recast
the volume potential as a sum of single- and double-layer potentials and a
volume integral with a regularized (bounded or smoother) integrand. The layer
potentials can be accurately and efficiently evaluated everywhere in the plane
by means of existing methods (e.g. the density interpolation method), while the
regularized volume integral can be accurately evaluated by applying elementary
quadrature rules. We describe the method both for domains meshed by mapped
quadrilaterals and triangles, introducing for each case (i) well-conditioned
methods for the production of certain requisite source polynomial interpolants
and (ii) efficient translation formulae for polynomial particular solutions.
Compared to straightforwardly computing corrections for every singular and
nearly-singular volume target, the method significantly reduces the amount of
required specialized quadrature by pushing all singular and near-singular
corrections to near-singular layer-potential evaluations at target points in a
small neighborhood of the domain boundary. Error estimates for the
regularization and quadrature approximations are provided. The method is
compatible with well-established fast algorithms, being both efficient not only
in the online phase but also to set-up. Numerical examples demonstrate the
high-order accuracy and efficiency of the proposed methodology; applications to
inhomogeneous scattering are mentioned.
