Discretization of the divergence operator

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?

Well I assume you must be dealing with something like an incompressible fluid if the only parameter you need to use is the velocity.

For an incompressible fluid, the divergence of a point measures the net amount of fluid coming out of that point, so unless there is a "source" or "sink" of fluid in that region, then there is no divergence. If you discretize the fluid into cubical regions, where each cube has a given velocity vector, then you can simply calculate the volume of fluid each cube sends through each of its square faces. If a given face separates cubes A and B, then the relevant quantity for that face is (volume A sends through face)+(volume B sends through face). The volume of incompressible fluid passing through the face is equal to the dot product of the velocity vector with the surface's normal vector (oriented outwards from the core of the cube). If, for example, the face separating cubes A and B is parallel to the y-z plane, then the normal vector is [itex]\hat{\mathbf{x}}[/itex] for cube A and [itex]-\hat{\mathbf{x}}[/itex] for cube B, so the relevant quantity for the face is [itex]\mathbf{v_A}\cdot\hat{\mathbf{x}}-\mathbf{v_B}\cdot\hat{\mathbf{x}}[/itex]. So to get a measure of the divergence for an entire cube, you might average the six values for each of the faces.