The blob is a polyhedron whose faces are called ‘facets’ here.

Names in the code

Globals

Pt is a triple of reals that serves both as a 3D vector and to locate a point.
T is a triple of vertex ids for some facet.
n as in “Each triangular face is divided into n2 initially congruent triangles.”
h[j] holds the head of the chain of slivers around vertex j.
sph.c

This code subdivides each of the 20 faces of an icosahedron into n2 triangles. This results in VN = 10n2+2 vertices whose coordinates are recorded in array q, and SN = 20n2 triangles each terminating its own sliver.
q[j] holds coordinates for vertex j.
Q[j] holds indices for the 3 vertices that bound sliver j.
L[3j], L[3j+1] & L[3j+2] provide link space for sliver j. Chain flows thru L.
m[j] is the mass of sliver j.

For each vertex v we often need to visit each sliver that v bounds. Each chain link is an index into L, or –1 which ends the chain. The head of the chain for vertex k is in h[k]. The link following k is in L[k]. If the chain rooted in h[i] includes 3k+j and 0⩽j<3 then sliver k includes vertex i. The chain provides each sliver index by following the chain rooted in h and linked thru L. This requires a bit of combinatorial topology.

det(p, q, r) returns the determinant of three P vectors.
sph() initializes a Bucky ball and returns n.

gen.c

nv takes an integer identifying a facet and returns the vector which is the normal to that facet.
centr(j, q) returns the location of the center of facet j.
add & sub return the sum and difference of two vectors.
ip returns the inner (dot) product of two vectors (Pt’s).
R det(Pt a, Pt b, Pt c) returns the determinant of the three vectors.
rf() returns a random float distributed uniformly between 0 and 1.
len(Pt x) returns the length of its argument.
sm(R s, Pt b) returns the product of scalar s and vector b.
cp returns the cross product of two vectors.
prP prints an annotated vector.
R moments(R m1[3], R m2[3][3]) returns moment 0 (the mass) and puts moments 1 and 2 into the two arguments.

puff takes a vertex index and returns the vector sum of the impinging faces. This routine illustrates how the information in h and L may be used to find the facets about a vertex. If that information is needed elsewhere we may abstract this logic.

files with a main()

fac.c

There is now code in main of fac.c to compute moments of the blob, m0, m1 & m2. This code can be moved to a routine when we learn when we want to compute these moments.

test.c Early code to test various functions and initially confirm invariants.
P(t) returns the potential at point t. It is good to only about 1 part per thousand(2n2).
rnd() returns a random normal deviate.
(This tests some off axis formulae that I think I don’t need now.)
blob.c

This code is awkward because C does not pass an array of three R’s by value so I must use type Pt. But I can’t get the math right without tensor notation. Pt a2Pt(R m[3]) turns an array into a Pt and void Pt2a(R a[3], Pt p) stores a Pt into the array a.

This version tries to be cache savvy and keep k*6 floats in cache to avoid 20*n4 cache misses each evaluation of surface potentials. The constant k affects cache efficiency but not the results. I chose k as 10 for my Mac.
rm(char *, m0, m1, m2) merely prints moments 0, 1 and 2.
matrot(m) rotates the vertices about the origin by matrix m.
a2Pt and Pt2a convert between structure Pt and an array of 3 reals.
A performance hack?

Tire Balancing
When a new tire is placed on a car’s wheel it must be balanced so that when the wheel with tire spins about the axis defined by the wheel, there are no periodic forces on the axel. We have the same problem after we have shaved off a bit of blob and moved it. We have a new blob which we must balance. The line of code with comment ‘center’ translates the modified blob so the center of mass is at the origin again. But this translation forces us to adjust the m2 by this formula which we do in the block with the comment ‘turn’. This translation and rotation may be viewed as choosing a new coordinate system; an alias for our blob. The computation of this “theoretical’ tensor needs to be turned with matrix r but I think It will again spin about the z axis. we don’t need that tensor, so for now I won’t finish that computation.