I fixed a problem in my previous 3D dual contouring script today.

For a given mesh grid constructed by numpy.meshgrid, for example:

import numpy as np

x = np.linspace(-4, 4, 5)
y = np.linspace(-4, 4, 5)
z = np.linspace(-4, 4, 5)

X, Y, Z = np.meshgrid(x, y, z)

def f(x, y, z):
    return x + y + z

F = np.vectorize(f)

s = F(X, Y, Z)

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:

3D dual contoured heart function

The source code of this script can be found here. It still has some topology issues to be resolved.