Consider the sliver which is a tetrahedron with three short edges which define a triangle on the surface, and three long edges that meet at the origin. This sounds as if the origin must be inside the blob but that is not necessary. We approximate this tetrahedron with a tapered stick, or sliver below because the tetrahedron’s potential is too messy to integrate.

Consider a sliver from (0, 0, 0) to (1, 0, 0) with density x² at (x, 0, 0). Its mass is ⅓. We call it the “unit sliver”. The gravitational potential at point (x, y, 0) is p(x, y²) =
∫[λ=0: λ=1] λ²/√((x−λ)² + y²) dλ = (3x(f−g) + f + (2x² − y²)(log((f−(x−1))/(g−x))))/2 where:
f = √((x−1)² + y²)
g = √(x² + y²)

f is the distance from the pointy end of the unit sliver and g from the heavy end at the surface. The y² in p(x, y²) is merely a hack to avoid an unnecessary square root.

Turns out that we don’t need even this simpler formula but it is nice to have around.

This code corroborates these formulae. (slp3 x y z) is the potential of the unit sliver at point (x, y, z).

That is the unit sliver but we need to compute the potential field for slivers rotated about the origin of various thickness and lengths. The volume V of the native sliver is easy to compute and V=m since we presume that the density is 1. Also available is the vector A from the origin to the centroid of the small face of the tetrahedron for which the sliver is a surrogate. We need a computable function P(A, B) so that 3mP(A, B) is the potential at point B, produced by the sliver with mass m whose heavy end is at A.
1. P((1, 0, 0), (x, y, 0)) = p(x, y²)
2. P((1, 0, 0), (x, a, b)) = P((1, 0, 0), (x, c, d)) if a² + b² = c² + d²
3. P(UA, UB) = P(A, B) where U is any orthogonal matrix.
4. P(sA, sB) = s²P(A, B) for any real s.

Equation 2 merely expresses the symmetry of the configuration. Equation 3 follows from freedom to choose any orthonormal coordinate system.
Geometric intuition tells me that
P(A, B) = (A∙A)p((A∙B)/(A∙A), ((A×B)/(A∙A))∙((A×B)/(A∙A)))
but see this which relies on this framework. os is a random 3 by 3 orthogonal matrix, mxm is matrix multiply, (ip A B) is inner product A∙B, (cp A B) is Gibb’s cross product A×B, (sm x A) is scalar multiplication of scalar x by vector A, P and p are as here. We confirm our defining equations.

(older obsolete work)

Data Plan

Slivers and vertices are each numbered with small integers which serve as array indexes. Data for the vertices are 3 floating coordinates, and head of a linked list of indexes for the vertex’s neighboring slivers. Data for the slivers are an array of the three indexes of the bounding vertices, and redundantly, the sliver’s mass. Also the three coordinates of the centroid of the surface triangle are kept.

One array, L, whose length is three times the number of slivers, holds three entries per sliver. Each entry is an index into L or −1. This is where the linked lists live. Each sliver has 3 defining vertices and must provide link space for those.

The intense computation is computing the potential produced by the collective slivers at each vertex. As we seek the minimal energy configuration we adjust two or more vertex locations at a time and recompute vertex potentials everywhere. Probably it is good to keep the potentials and compute the adjusted potentials everywhere by computing the potential change due to a vertex change. (Maybe a SB*VN array of potential contributions of sliver j to vertex k??)

See this too.
Generate a spherical blob from triangles.
potential function with two vector arguments, corroboration of same
Physics