We present an accurate and rapid solution of Poisson's equation for space-filling, arbitrarily shaped, convex Voronoi polyhedra (VP); the method is O(NVP), where NVP is the number of distinct VP representing the system. In effect, we resolve the long-standing problem of fast but accurate numerical solution of the near-field corrections, contributions to the potential due to near VP—typically those involving multipole-type conditionally convergent sums, or use of fast Fourier transforms. Our method avoids all ill-convergent sums, is simple, accurate, efficient, and works generally, i.e., for periodic solids, molecules, or systems with disorder or imperfections. We demonstrate the practicality of the method by numerical calculations compared to exactly solvable models.