I'm working on a fluid simulation project where two spheres are simulated using smoothed particle hydrodynamics (SPH). I have generated two sets of particles inside and on the surface of a sphere. You can see the result in this picture:
scatter plot of spheres' particles
I am going to extract the surface of the particles by the marching cubes algorithm. Currently, I am using "scikit-image" implementation of marching cubes in python, and I got this result:
Reconstructed surfaces and particles
In the above picture, two small spheres are scattered plots of my particles, and two large surfaces are the reconstructed surface that I got from the marching cube algorithm.
I have two problems:
1- As you can see, the surfaces are not smooth enough. So, I'm looking for an efficient method to make these surfaces smoother with a more sphere shape. The function that I used for surface reconstruction of my SHPs is:
Phi(x) = - c + sum_{j in N(x)} of (V_j * W(x - x_j, h))
where c is the surface threshold, x is the position where we want to evaluate, sum_{j in N(x)} is the sum over all particles j in the neighborhood of point x, V_j is the volume of particle j, W is the SPH kernel function, x_j is the position of particle j, and h is the kernel smoothing length
2- The surfaces are supposed to exactly lie on particles as they are constructed to cover the outer surfaces of the particles. But as you can see in the picture, two surfaces are located somewhere far from the particles. I think my voxelization is not correct.
The code that I have written for surface reconstruction is:
def surface_reconstruction(rezolution=20, surfacelevel=0, bounds=(-4, 4)):
# First I creat my voxel grid
x = np.linspace(bounds[0], bounds[1], rezolution)
y = np.linspace(bounds[0], bounds[1], rezolution)
z = np.linspace(bounds[0], bounds[1], rezolution)
xx, yy, zz = np.meshgrid(x, y, z)
surface = np.zeros((len(xx), len(yy), len(zz)))
# Then I loop over each voxel's corner and find a value for that point to check whether
# it is inside the mesh or no?
for i in range(len(xx)):
for j in range(len(yy)):
for k in range(len(zz)):
for p1 in env.particles:
# p1.X is the position of particle p1. Here, I compute the distance
# between each particle and voxel coordinate
dX = p1.X - [xx[i, j, k], yy[i, j, k], zz[i, j, k]]
r = np.linalg.norm(dX)
# p1.H is a kernel radius of particle p1
if r <= p1.H:
# p1.vol is particle p1 volume
# Wij is a kernel function of SPH
surface[i, j, k] += -surfacelevel + p1.vol * Wij(r, p1.H)
return surface
Then I apply marching cubes on "surface"
verts, faces, normals, values = skimage.measure.marching_cubes(surface_reconstruction(), 0)
So, I am wondering
1- How can I make the surfaces smoother?
2- Why surfaces are larger than particles scale and far from particles?
