GEOS
GEOS copied to clipboard
Robust point location for degenerate configurations
This PR adds a new method to determine if a point lies inside or outsilde a (convex) polyhedral cell. Compared to the current method, the new implementation is slower, but handles gracefully the degenerate conditions of points lying on mesh vertices, edges and faces. It does so by:
- Triangulating each face of a cell in a consistent manner, by always choosing the lowest global-id vertex as root;
- Ordering each triangle again using global ids, so that the edges and vertices are encountered in a consistent order across MPI ranks;
- Computing the winding number (an integer) subtended by each triangle;
- If the point lies on a vertex, edge or face, this will be found consistently by all ranks, whatever the cell from which the face/edge/vertex is seen. This is true because the order in which the vertices are traversed (and thus the result of the calculation) is exactly the same no matter the originating cell or MPI rank;
- Finally, if a degenerate condition is found, the {face,edge,vertex}-to-element maps of the host element are used to choose consistently a single representative element. The adjacent element with the lowest global ID is chosen, ensuring general consistency across MPI ranks.
The only ambiguity left is when the point lies on a face/edge/vertex that is a boundary for an MPI rank, but an interior face/edge/vertex for another rank. In this case, both ranks will claim the point. However, at least one of them will consider the point to belong to a ghost element. This is a desirable outcome, since the desired behaviour in this case can be different according to the specific application (e.g., sources vs. receivers for the seismic wave case). This outcome leaves the choice to the caller method on how to handle this case, using ghost rank information.
The new point location algorithm is used for wave propagation kernels, fixing issue 3104 and issue 2888.