VectorTools::point_value() throw exception: point not locally owned

[yiyangzhang37]

The problem is described in deal.II google group. here.

The return of GridTools::find_active_cell_around_point() is ambiguous if the input is a vertex. Such a test seems to be present in the function VectorTools::point_value(). If working with parallel::shared::Triangulation, then we will get an exception saying that the point is not locally_owned even if we have done the cell->is_locally_owned() test.

Attached a test code.

locally_owned_cell_test.zip

[bangerth]

I don't think there is actually a good solution. You are asking that the algorithm tries to find a locally-owned cell when given a point that happens to be a vertex that is shared between a locally owned and a ghost cell.

But if we apply this reasoning to the two processors that are owning these two cells, then that would mean that the algorithm should find different cells on different processors when given the same point. This seems like a generally bad idea -- given floating point inaccuracies, the only useful solution is to declare that the point is either in one of the two cells, or in the other. To have the two processors return different answers is asking for trouble.

[Rombur]

I think the problem is that GridTools::find_active_cell_around_point() only returns one cell. When you ask for a point on a face, an edge, or a vertex it should return more than one cell since there is more than one cell "around" the point.

[bangerth]

But in floating point arithmetic, a point cannot be "on a face or edge". It could be on a vertex. But in practice, I think the better practice is to use points that are just perturbed a tiny bit to avoid the ambiguity -- just move the point by 1e-15 in a random direction and there will be a unique cell in which it is located.

[bangerth]

By the way, if you want to find all cells adjacent to a particular vertex, then there is already GridTools::find_cells_adjacent_to_vertex().

[Rombur]

But in floating point arithmetic, a point cannot be "on a face or edge".

Well that depends on your mesh. Not everybody is using complex spherical meshes some of us are using simple and nice slab geometry :-).

But in practice, I think the better practice is to use points that are just perturbed a tiny bit to avoid the ambiguity -- just move the point by 1e-15 in a random direction and there will be a unique cell in which it is located.

I agree this is the easiest solution.

By the way, if you want to find all cells adjacent to a particular vertex, then there is already GridTools::find_cells_adjacent_to_vertex().

But @jtceleron is not calling find_active_cell_around_point() directly. He calls VectorTools::point_value() which itself uses find_active_cell_around_point().

[yiyangzhang37]

I think the function find_active_cell_around_point() is okay, it has its own way to determine one vertex belongs to one cell, and is deterministic -- this is fine.
The problem is that it seems in the function VectorTools::point_value() there is a test of something like
if( find_active_cell_around_point( vertex )-> is_locally_owned() )
which is too strong.
It would be better if this test is replaced with something like
if ( GridTools::find_cells_adjacent_to_vertex(vertex) ->is_one_of_the_cells_locally_owned() )

[Rombur]

@jtceleron You cannot simply change the assertion because we need a cell that is locally owned. What you propose requires to get all the cells that are around the vertex and then pick one that is locally owned. That means changing find_active_cell_around_point() to return all the cells. It would be easier if you just slightly move your point so that it is included in one only cell.

[tjhei]

VectorTools::point_value() should also work if the cell is a ghost cell (assuming we have the locally relevant DoFs). I assume we currently do not support this so that exactly one processor will return a result, right? But this is also fragile, because we depending on the cached last cell, so we might still run into a situation where no or both processors will return a result, right? (Assume a point is on the interface between two cells and test as "inside" for both)

just move the point by 1e-15 in a random direction

That is not a satisfying solution for point_value. I would rather be more generous and allow ghost cell evaluation. Then the user can decide what they should do if more than one processor find a result.

[bangerth]

But VectorTools::point_value() doesn't cache a cell.

I understand that from a user perspective, the situation may be awkward, but it is consistent and deterministic. So I don't think that there is a bug or undesirable behavior in at least the original testcase provided in the first post above. If there is a problem in VectorTools::point_value(), can someone construct a testcase that shows this?

[drwells]

This is yet another issue related to find_active_cell_around_point() which is being tracked elsewhere.