Important note: This answer does not answer the actual question, but it was left undeleted per request. Embarrassingly I confused hexahedral and hexagonal. The question is about sorting points into arbitrary hexahedral cells in 3D while this solution sorts points into regular hexagonal cells in 2D, or irregular ones that correspond to some Voronoi tesselation in any dimension. This method is applicable only if the mesh was generated as a Voronoi tesselation in the first place (which does seem to be an occasionally used approach).
I am not sure what you mean by sort here, but I assume you want to sort the point into hexagonal bins on the plane.
Mathematica is what I know, so I am going to show you how to do it in Mathematica, but the method can be ported to other systems. The idea is that a hexagonal lattice is the dual of a triangular one: it can be generated as the Voronoi diagram of a points in triangular arrangement. A point from the cloud belongs to a given hexagon if it is closer to the centre of that hexagon than to the centre of any other hexagon.
This method will work for meshes of different shapes too, as long as they can be generates as the Voronoi diagram of some point arrangement. (E.g. the hexagons don't need to be regular.)
Let's generate the mesh. This is a triangle lattice:
pts = Join @@ Table[{x, Sqrt[3] y}, {x, 0, 4}, {y, 0, 2}];
points = Join[pts, TranslationTransform[{1/2, Sqrt[3]/2}] /@ pts];
Needs["ComputationalGeometry`"]
PlanarGraphPlot[points, LabelPoints -> False]

Its dual is the hexagonal one we are interested in:
DiagramPlot[points, LabelPoints -> False]

This builds a function nf
which finds the index of the hexagon centre that some cloud point is closest to. It is the key to the method:
nf = Nearest[N[points] -> Range@Length[points]];
Now let's generate a cloud of 1000 random points and sort them with nf
:
cloud = RandomReal[{-1/2, 5}, {1000, 2}];
indices = First /@ nf /@ cloud;
indices
contains the indices of the centres that each cloud point is closest to. This is the information we needed. Now we can make a histogram out of them ...
Histogram[indices]

... or colour each of them ...
Show[
DiagramPlot[points, LabelPoints -> False],
Graphics@MapThread[{ColorData[3][#1], Point[#2]} &, {indices, cloud}],
PlotRange -> All, AspectRatio -> Automatic
]

... or do any sort of fancy visualization we want.
tally = Tally[indices];
ListDensityPlot[Join[points, List /@ Sort[tally][[All, 2]], 2],
InterpolationOrder -> 0,
Epilog -> (Text[#2, points[[#1]]] & @@@ tally),
PlotRange -> {{-.5, 5}, {-.5, 5}}, Mesh -> All,
ColorFunction -> (ColorData["BeachColors"][1 - #] &)]

The key point here was the function that finds the closest point to something (Nearest
). Mathematica has this built in, but there's a chance your system doesn't. If this is the case, please see this question on how to efficiently implement such a function (or just go with the naive linear time implementation if you don't have a huge amount of points to process).