The gravitational field of a cube” gives a formula for the gravitational potential in and about a cube of constant density. This code corroborates the formula. In the code u(x,y,z) = u(x,z,y) = u(z,x,y).

V(x, y, z) is the potential both inside and outside the cube [−1, 1]3.

The following seems intuitively true, but it is false:

If u were a potential field then at each point Σi2u/∂xi2 = 4π(density). The numerical test for this fails badly.
It just now occurs to me that this result greatly reduces the cost of computing the potential field of a constant density peculiar shape. Use oct trees to decompose the shape. Perhaps that was the stimulus for the result. Indeed rereading the introduction makes that clear but without mentioning oct trees.

This may allow computing the shape of an accreting incompressible fluid. We know that the answer is a sphere but it would be good to get a 2nd opinion. With that we could examine the shape of a rotating blob of incompressible fluid—the stuff of Roche limits. I want pictures! See this.