  We introduce two new methods to obtain reliable velocity field statistics
from N-body simulations, or indeed from any general density and velocity
fluctuation field sampled by discrete points. These methods, the {\it Voronoi
tessellation method} and {\it Delaunay tessellation method}, are based on the
use of the Voronoi and Delaunay tessellations of the point distribution defined
by the locations at which the velocity field is sampled. In the Voronoi method
the velocity is supposed to be uniform within the Voronoi polyhedra, whereas
the Delaunay method constructs a velocity field by linear interpolation between
the four velocities at the locations defining each Delaunay tetrahedron. The
most important advantage of these methods is that they provide an optimal
estimator for determining the statistics of volume-averaged quantities, as
opposed to the available numerical methods that mainly concern mass-averaged
quantities. As the major share of the related analytical work on velocity field
statistics has focussed on volume-averaged quantities, the availability of
appropriate numerical estimators is of crucial importance for checking the
validity of the analytical perturbation calculations. In addition, it allows us
to study the statistics of the velocity field in the highly nonlinear
clustering regime. Although in most cases both the Voronoi and the Delaunay
method lead to equally good velocity field estimates, the Delaunay method may
be slightly preferable. In particular it performs considerably better at small
radii. Note that it is more CPU-time intensive while its requirement for memory
space is almost a factor 8 lower than the Voronoi method.
