I fixed a problem in my previous 3D dual contouring script today.
For a given mesh grid constructed by
numpy.meshgrid, for example:
A coordinate of
(x0 + step * x, y0 + step * y, z0 + step * z),
whose integer grid coordinate is
(x, y, z),
needs to be indexed from the mesh grid using
s[y, x, z].
Now the 3D heart function can be successfully contoured:
The source code of this script can be found here. It still has some topology issues to be resolved.