  Well-conditioned boundary integral methods for the solution of elliptic
boundary value problems (BVPs) are powerful tools for static and dynamic
physical simulations. When there are many close-to-touching boundaries (eg, in
complex fluids) or when the solution is needed in the bulk, nearly-singular
integrals must be evaluated at many targets. We show that precomputing a linear
map from surface density to an effective source representation renders this
task highly efficient, in the common case where each object is "simple", ie,
its smooth boundary needs only moderately many nodes. We present a
kernel-independent method needing only an upsampled smooth surface quadrature,
and one dense factorization, for each distinct shape. No (near-)singular
quadrature rules are needed. The resulting effective sources are drop-in
compatible with fast algorithms, with no local corrections nor bookkeeping. Our
extensive numerical tests include 2D FMM-based Helmholtz and Stokes BVPs with
up to 1000 objects (281000 unknowns), and a 3D Laplace BVP with 10 ellipsoids
separated by $1/30$ of a diameter. We include a rigorous analysis for analytic
data in 2D and 3D.
