# 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)