I work with a grid-based code, this means that all of my quantities are defined on a mesh. I need to compute, for every point of the mesh the divergence of the velocity field. All I have is, for every cell of my mesh, the values of the 3-d velocity in his 26 neighbors. I call neighbors the cells with center (i,j,k) with i,j,k=-1,0,1 and (0,0,0) being the center of the cell. What's the best discretization of the divergence operator under these assumptions?