How to compute the value of a function at specific points using MPI in a parallel computing scenario

I used dolfinx to solve the Poisson equation and obtained the solution uh using MPI parallel computing. I want to know how to compute the value of the solution at a specific point in the domain. I tried the following function, but the results were all NaN. Are there any examples or tutorials on how to compute the value of the solution at specific points in an MPI environment?
def insert_value(self, u, points):
nr_points = len(points)
u_values = np.ones((nr_points,))
cells =
points_on_proc =
tree = bb_tree(self.mesh, self.mesh.geometry.dim)
cell_candidates = compute_collisions_points(tree,
points) # Find cells whose bounding-box collide with the the points
colliding_cells = compute_colliding_cells(self.mesh, cell_candidates,
points) # Choose one of the cells that contains the point

    indices = np.ones((nr_points,), dtype=int) * (-1)
    # roots = np.ones((nr_points,), dtype=int)
    # Each process creates an array of size 3 with values equal to its rank
    sendbuf = np.full(nr_points, 0, dtype=int)
    u_values_send = np.full(nr_points, 0, dtype=float)
    # Root process initializes the receive buffer to gather all arrays
    if self.myrank == 0:
        recvbuf = np.empty(self.size * nr_points, dtype=int)
        u_values_recv = np.full(nr_points * self.size, 0, dtype=float)
        recvbuf = None
        u_values_recv = None
    count = 0
    for i, point in enumerate(points):
        if len(colliding_cells.links(i)) > 0:
            sendbuf[i] = 1
            indices[i] = count
            count += 1
    points_on_proc = np.array(points_on_proc, dtype=np.float64)
    u_values_temp = u.eval(points_on_proc, cells)
    for i in range(nr_points):
        if sendbuf[i] == 1:
            u_values_send[i] = u_values_temp[indices[i]]
    # Gather the arrays at the root process
    self.comm.Gather(sendbuf, recvbuf, root=0)
    self.comm.Gather(u_values_send, u_values_recv, root=0)
    if self.myrank == 0:
        for i in range(nr_points):
            for j in range(self.size):
                if recvbuf[i + j * nr_points] == 1:
                    u_values[i] = u_values_recv[i + j * nr_points]
    u_values = self.comm.bcast(u_values, root=0)
    return u_values

Please make a minimal reproducible example.
The code you have procduced here is

  1. Not properly formatted.
  2. Is not self-contained.

Have you looked at Implementation — FEniCSx tutorial
which shows how to get the local values on each process. Have you checked that these make sense?