Thanks for the feedback!
I was able to implement the Bowyer-Watson algorithm for the Delaunay triangulation of 3D points into a tetrahedral mesh:
I also started writing a paper to document all of the math and stuff. There's quite a bit more to it than meets the eye. The paper is very incomplete and a WIP, but it will improve in time.
4 is all you need. In the real world where you can derive a continuous function for mapping position->gravitational vector, you only need 1 directional vector at any point in space. The infinite number of other vectors don't "influence" you, no matter how far or close they are to you.
Remember, this is a vector field, meaning that it is the solution to the differential equations created by multiple bodies of mass. Since it's not possible to have an infinite number of vectors, I have to create a discrete vector field and interpolate (read the introduction to my paper, I think it explains it pretty well).
I too could approximate the vector field a planet's gravitational field would create by spacing the gravity vectors at inverse proportional distances from the planet. It might look something like this: http://i.imgur.com/2fOVpxa.png
Using normals isn't going to cut it for what we have planned. We have lots of cases where gravity won't be collinear with the face's normal vector. That's why I've created a system where you can arbitrarily define gravity, no matter what the mesh looks like.