Is there a way to compute the "correct" distributional answer, -4*pi*dirac_delta(r), as the Laplacian of 1/r (in spherical coordinates on Euclidean space)?

edit

Background: I'm putting together an EM notebook for educational purposes, and would like to be able to show demonstrations of Gauss's Law (for example). However, for completely understandable and obvious reasons, integrating (1/r^2*r_hat).div() over the unit ball just gives 0. Computing the flux of (1/r^2*r_hat) through the unit sphere gives the expected result, of course, but I'd hoped to directly demonstrate that it's proportional to the volume integral.