I am inspired by this paper (which gives the gravitational potential field of a cube analytically) to derive the potential of a 3D simplex. Decomposing a body into simplices has some advantages over cubes. The plan is to derive the potential here and concurrently verify results with Scheme code here and accompanying notes. That code catches my many errors before they propagate.

In each case the uniform density of the n-simplex is 1 and its mass is 1/n!. r2 = X2 + Y2 + Z2.

The 3D field produced by a 0-simplex at the origin is 1/r. The potential field transforms with Euclidean similarities but not affine.

For the 1-simplex {(x, 0, 0) | 0≤x≤1} is already slightly complex: p1(X, Y, Z) = ∫[x=0:x=1] dx/√((X-x)2 + Y2 + Z2) at the point (X, Y, Z).
∫[x=0:x=1] dx/√((X-x)2 + Y2 + Z2) = ∫[x=0:x=1] dx/√(x2 − 2Xx + r2)
= (log(2√(r2+(−2X)x+x2) + (−2X) + 2x))[x=0:x=1] = log(2√(r2+(−2X)+1) + (−2X) + 2) − log(2√(r2) + (−2X))
= log(2√(r2 − 2X + 1) − 2X + 2) − log(2r − 2X) = log((2√(r2 − 2X + 1) − 2X + 2)/(2r − 2X))
= log((√(r2 − 2X + 1) − X + 1)/(r − X)).
See p1 here.

To explore simplification we consider pd(X, Y, Z) = p(X, Y, Z)+p(−X, Y, Z) = p((X+1)/2, Y/2, Z/2). This is the potential for {(x, 0, 0) | −1≤x≤1} with mass = 2.
See (pq x y z) = (let ((f (lambda (x) (- (sqrt (+ (sq x) (sq y) (sq z))) x)))) (log (/ (f (- x 1)) (f (+ x 1)))))).

We choose the triangle with vertices at (0, 0, 0), (1, 0, 0), (x', y', 0). Any triangle is, for some x' and y', similar to this. The potential at (X, Y, Z) is:
p2(x', y')(X, Y, Z) = ∫[λ=0, λ=1] p1((X−λx')/(1−λ), (Y−λy')/(1−λ), Z/(1−λ))


Celestial Mechanics and Dynamical Astronomy (ISSN 0923-2958), vol. 59, no. 3, p. 253-278
http://www.wseas.us/e-library/conferences/2011/Penang/ICOPOW/ICOPOW-06.pdf