Hello everyone,

I’m new to fenicsx, and when I try out some functionalities recently, something confuses me. I have read the tutorial of fenicsx and found that u.eval(points, cells) outputs the values of given points. According to the doc, the cells[i] should contain the points[i]. If my understanding is right, the output will be zero if cells[i] doesn’t contain the points[i] because the basis functions of cells[i] will be zeros outside it.

However, things don’t go as expected. I found that eval(points, cells) will output some values whatever cell indexes are assigned. I wonder how eval(cells, points) actually works. Here are some simple reproduction examples:

```
nx, ny = 2, 2
domain = mesh.create_rectangle(MPI.COMM_WORLD, [np.array([0., 0.]), np.array([1., 1.])], [nx, ny], mesh.CellType.triangle)
V = fem.functionspace(domain, ("Lagrange", 2))
u = fem.Function(V)
u.interpolate(lambda x: x[0]**2+x[1]**2)
points = np.array([[0.1, 0.1, 0.], [0.6, 0.8, 0.0]])
u.eval(points, [0, 0])
u.eval(points, [2, 4])
u.eval(points, [7, 7])
```

These "eval"s will all output the correct answer [0.02, 1.0].

And when I change the Lagrange order from 2 to 1, that is, the target function cannot be interpolated perfectly, the outputs will change as the assigned cell indexes change. I have tried to read the C++ code but is still confusing.

An accompanied question is, the computing_colliding_cells function will usually output more than 1 cell indexes for each points. If any cell index input to eval() can give arise to a result, how can we assure that the correct cells are selected? I have read the C++ code of computing_colliding_cells and found that this function simply pushes the cell indexes whose "square_distance"s are smaller than a hard-coded tolerance (1e-20) to the output array. If the dimension of the geometry is small itself, this function doesn’t have the ability to perfectly distinguish among candidate cells.

I would appreciate it a lot if someone can explain what happens behind these functions.