Here is a scheme that is vastly simpler than any so far. Divide the surface of the blob into triangular facets. To compute the potential at point X divide the blob into units composed of the convex hulls of a facet and X. This is a sliver. The contributed potential of each sliver is easy to compute:
(∫[0:1]x2/x)/(∫[0:1]x2) = (½)/(⅓) = 3/2
P = (3/2)∫(x/|x|)∙da where da is the vector area of an element of the surface and x is the vector from X to the area centroid.
The potential is 3/2 times what it would be if the same total mass were concentrated all at the centroid of the facet.

Some dot products will be negative. Still (∫x∙da)/3 is the volume of the blob and ∫da = 0. da is the vector normal to the facet whose length is proportional to the area. We also compute moments 0, 1 and 2 of the blob, m0, m1 & m2.

Some code

(2017 Feb 1)
Moving the vertices around is like herding cats. They find too many interesting directions to go. There are just now two problems to solve: