With octrees and this surprising math it is easy to combine cubes into stranger shapes. For a shape with area B and shape veracity ε, B/ε2 cubes are needed.

Another approach for octrees is the surface integral of slivers. At least with octrees we need not guard against negative volumes. Nor need we deal with surface triangulation degrading.

(2017 Feb 1) Octrees are neat but I think I may have a simpler solution. Let me define z-convex to mean that for any x and y the line parallel to the Z axis with that x and y intersects the blob 0 or 2 times. I conjecture that all the solutions I need fit this rule. This results in a very simple and compact representation; two 2D arrays, top and bot, of integers. The 0 penetration case is marked somehow. It is very easy to generate the surface of this blob as a set of unit squares, each with one of six orientations, with integral coördinates, The ends of the z-lines are the principle candidates for additions or deletions. If we assume that the X-Y profile is convex then we must only be prepared to add new planes for max or min X, ditto Y.

We use the ideas developed in the sliver scheme with square facets in place of triangular facets and scan the blob surface, square by square to compute the potential at a surface point.

I sort of dislike adding these presumptions, but I want to see the pictures! Do we need potentials only at the ends of the Z-lines?

Further conjecture: the solutions are preserved if we reverse the sign of z. This means we need only one array! The solution is preserved if we rotate about the Z-axis or reflect in the XZ coördinate plane. Every factor of two helps. Similarly reversing the sign of ω preserves solutions because we assume no motion in relative to the revolving frame and thus no Coriolis force. Thus we consider only those cubes with positive y.

We need to consider both signs for x in order to work on the 2 blob solution for large ω.

This code suggests that an array of shorts is efficient in the floating context.

Combining the above ideas I need one 2D array, Z, of shorts indexed by integer x and y, x may be negative but y is positive. Z[x, y]≥0.
The blob is {(x, y, z) ⎮ |z| < Z[⌊x⌋+ox, ⌊|y|⌋]. ox is a bias to the x subscript since x must range over the interval [−ox, XM]. We will take x=y=0 as the axis of spin. We must describe negative x’s. Some code.

Lets consider computing potentials at the centers of those cubes, in the blob, with at least 3 facets showing. Presumably those with fewer outside facets are not yet candidates for removing. But where do you put them? If we compute potentials at cube centers for cubes not in the blob, but either:

As a check we should presumably never be in a configuration where there is a cube outside the blob adjacent to more then 3 in the blob.

We need a way to visit each of the places where we must compute the potential.

If the goal is to find the maximum assembly energy configuration holding ΣZ constant, then we should omit from the summation the self energy of each cube in the blob. The symmetry imposed by these conjectures means that our inertia tensor will always be diagonal.

Definition: a ‘vertex of the blob’ is where one of the 3 coordinates takes on a maximum, or minimum in the case of x. Special logic is necessary at vertices.
When x is at either extremum, y and z are both 0;
When y is at its max, z=0 and x is not determined by symmetry;
When z is at its max, y=0 and x is not determined by symmetry.

I think I will try a sparse selection of hills and holes.
Whenever 4Z[i, j] > Z[i+1, j] + Z[i−1, j] + Z[i, j+1] + Z[i, j−1] + 1 then we compute potential at (i+½, j+½. Z[i, j]+½) as hill.
Whenever 4Z[i, j] < Z[i+1, j] + Z[i−1, j] + Z[i, j+1] + Z[i, j−1] − 1 then we compute potential at (i+½, j+½. Z[i, j]−½) as hole.
Note the duality. Other special cases are necessary for the vertices of the blob, at least. I am not confident that this is enough. It is time for the computer to think about it.