  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. 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 presented.
