"eval(points, cells)" doesn't work as expected

dolfinx
1

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.

2

For the second question, I have tried out some methods and found that scaling will help. I suggest always to scale before implementing subsequent operations, especially when working with small systems, because there are some hard-coded tolerance in the code. However, I’m still confused with the “eval” function.

3

The code gives the same result always as what happens is that when you send in a cell where the point isn’t located, it does an extrapolation when evaluating (and since your solution function is quadratic and you are using a quadratic space, you get the exact result). This is due to the fact that adding a collision detection algorithm inside the eval would cause a major slow-down, and the user is expected to have determined the correct cell a priori (for instance with compute_colliding_cells).

For instance, if you change your input space to be linear, you will get different results for different cells.

You shouldn’t use cell geometries whose coordinates are of the magnitude 10^(-10) in general, as you are touching close to double machine precision (as for instance when computing a Jacobian, one will take the product of such small numbers, resulting in something with less than machine tolerance).
As you say, scale your problem geometry (and subsequent equations) to avoid this.

The eval function only skips cells whose index is negative, i.e. if you send in

print(u.eval(points, np.array([0, -1],dtype=np.int32)))

the latter is skipped.

1 Like
4

Thank you very much for the detailed reply! I understand how the “eval” function works now.