I guess you are assuming a structured grid (as depicted), as it wouldn’t make sense to do this on an unstructured mesh that is not aligned with the coordinate axis.
What you could do is to build a distance matrix for every internal degree of freedom, similar to: How to represent non-uniform first and second type boundary conditions? - #6 by dokken
and find the 8 closest points.
However, as your grid is structured, I guess you know the location of the neighbouring points. If so, I would simply search for them in tabulate_dof_coordinates and extract the index of the matching degree of freedom (accessed in u.x.array).