I think that the goal is to find a shape of the given volume and angular momentum that minimizes the total energy. I would guess that this means that E = (the gravitational potential minus the potential energy due to centrifugal force) is uniform on the surface. For now I consider this a hypothesis and concentrate on moving surface mass from points of high E to points of low E, maintaining the angular momentum.
Here is a calculation plan. The data are a list, q, of vertex coordinates, and a topology in the form of another list, Q, of triangles each of whose vertices is from the first list, and an angular momentum, L, about the Z axis. The coordinates evolve with time while the topology remains constant. For large L the topology differs from fro small L—two blobs! I do not now have a plan for an automatic transition. Time here is not physics time but a path thru shape space towards configurations of lower energy. The energy loss comes from eliminating internal shear forces which we do not compute. We may or may not ever compute the total energy. Given the data it is easy to compute the tensor of inertia and its eigenvectors. We then do an orthogonal rotation of the vertices so that the angular velocity points in the same direction as the angular momentum. This means that the principle eigenvector is along the Z axis. This describes a body of the right volume and right angular momentum spinning stably, but perhaps with internal shear forces supporting a varying E at the surface. We compute E at each surface point and move mass by adjusting vertex coordinates radially from points of high E to points of low E. We stop when E is nearly uniform. The shape should then be our answer.
It is possible that the only solution is two blobs where E on the surface of each blob is uniform, but differs between the two blobs. This will take a change in topology. Our calculation must be sensitive to this possibility.
Calculation loop:
What happens when the potential gradient is a vector nearly tangent to the surface? If we are near a solution then the gradient is very small there and near the neck of the blob. Things get messy there.
One plan to avoid the code discovering that it can put two blobs in the same place is to maintain a plane which is static in the rotating frame and separates the two blobs. This is easy and cheap to maintain. It probably aids debugging. I conjecture that this eliminates no solutions.
When we have two blobs we must not move stuff from one to the other.
I think I will move each vertex some small fixed fraction of (the distance along the vector sum of the surrounding facet normals)/(length of sum of normals)*(the excess of the potential there over the average surface potential). This should very nearly conserve volume.