$$\Phi=\iint_{\partial V}\mathbf{g} \cdot d \mathbf{A}=-4 \pi G M$$
Essentially, why is $\Phi$ independent of the distribution of mass inside the surface $\partial V$, and the shape of surface $\partial V$? That is, I'm looking for a mathematical justification for the characterization of flux as a kind of "flow" of gravitational field through a surface.
The reason I worry about that is this. I have encountered the following proof of Gauss' Law for gravity: $$\iint_{\partial V}\frac{-MG}{r^2}\mathbf{e_r} \cdot \mathbf{e_r}d A=\frac{-MG}{r^2}\iint_{\partial V}d A$$ $$=-MG4 \pi$$ However, two assumptions are made in the above proof that I've not convinced myself hold true for the general case:
The Gaussian surface is a sphere.
It is assumed that all mass is concentrated at the centre of the sphere.
Returning to my original question:
- Why is it true that when the Gaussian surface is deformed (not adding or removing any mass to inside the surface, of course), the amount of flux, $\Phi$, through the surface is equal to the flux in the spherical case?
- Why is it true that the positions of the masses inside the surface doesn't affect $\Phi$?
Edit: removed reference to Poisson's equation, changed much of question to make it clearer.
Further edit: I am happy that $$\iint_{\partial V}\mathbf{F}\cdot d \mathbf{A}=\iiint_{ V}( \nabla \cdot \mathbf{F})dV$$
However, my problem is with the following equation: $$\nabla \cdot \mathbf{g}=-\rho G 4 \pi $$
Which can be used to get $$\iint_{\partial V}\mathbf{g}\cdot d \mathbf{A}=\iiint_{ V}(-\rho G 4 \pi)dV=-M G 4 \pi$$