@andre try again :)

@VsevolodA. Of couse one can do that numerically for a few points. If one is looking for the gravitational potential (as function) on a domain $\varOmega \subset \mathbb{R}^d$ and the density $\varrho$ is supported in (nearly) the whole domain, the complexity of FEM (with suitable linear solver, not Mathematica's standard direct solver) may be lower than numerical integration, in particular if one is satisfied with rather coarse, qualitative information.

I'd rather say, the problem of FEM in this context lies in the infinites of the domain. This can be cured by coupled FEM-BEM systems (BEM means boundary element method) where the interior of, say, a large ball of radius $R$ that encompesses all "complicated" events is modelled by finite elements and the interaction of the sphere of radius $R$ with the sphere at infinity is modelled by BEM (this involves only integrals over the spheres, not the domain in between.

All in all, both methods (the integration and FEM) have their value and their applications. No need to be snobbish.

@andre Good points. I am very with you, in particular with point 3): Easy it is not.