Consistent indexing between mesh and FE solution

Dear all,

I need some help in order to understand how to extract the map indexing linking the the order of the points of the mesh to the order of the FE solution (assuming polynomials of degree equals to 1).
The following code simply implements a voxelization of a FE solution by averaging the FE solution in multiple regions.
I need to find this information in order to modify the FE solution consistently. I think the critical line is the for loop over the indices. This is working with an increasing index that does not keep track of the mesh ordering. Additionally, I know that mesh.vertices returns the indices. However, to me it looks like they are not integer values as I would expect global indices to be. I am looking exactly for this piece of information.

def is_inside(point, min_x, max_x, min_y, max_y):
  return (point[0] >= min_x and point[0] < max_x and
          point[1] >= min_y and point[1] < max_y)

def voxelize(ngmesh, sol, s1D=5):
  # Find points
  new_sol = np.zeros(len(sol))
  points = np.delete( np.array([list(p) for p in mesh.ngmesh.Points()]), 2, 1 )
  # Define regions
  pt_x = np.linspace(-1.0, 1.0, s1D+1) # x-coordinate of points
  pt_y = np.linspace(-1.0, 1.0, s1D+1) # y-coordinate of points
  num_sensors = s1D * s1D # Number of sensors in the library
  # Voxelize each region
  for x_idx in range(pt_x.shape[0]-1):
    for y_idx in range(pt_y.shape[0]-1):
      nodes_idxs_inside = list()
      for vertex in range(mesh.nv):
        if is_inside(points[vertex], pt_x[x_idx], pt_x[x_idx+1], pt_y[y_idx], pt_x[y_idx+1]):
          nodes_idxs_inside.append(vertex)
      print(nodes_idxs_inside)
      nodes_idxs_inside = np.asarray(nodes_idxs_inside, int)
      nodes_idx_inside_mask = np.zeros(ngmesh.nv, dtype=bool)
      nodes_idx_inside_mask[nodes_idxs_inside] = True
      # Assume that the order of the FE solution is the same of the vertices and polynomials of degree = 1
      new_sol[nodes_idx_inside_mask] = sum(sol[nodes_idx_inside_mask]) / sol[nodes_idx_inside_mask].shape
      print(sol[nodes_idx_inside_mask], new_sol[nodes_idx_inside_mask])
  return new_sol

I hope I made myself clear.
Thanks for your help.

Best,
Francesco