Add an API for extracting cross-section information

[lucacarniato]

@veenstrajelmer would like to add an API that exports the geometrical information about cross sections. The use case is shown in this figure

For each crossed face the API should return its index AND the coordinates of the two intersections (green dots). The order must be increased along the line direction (in the figure above the order will be [0,1,5,6,7,11]).

Important notes:

• The face indices are provided in the netCDF file and the Mesh2D parameter in mkernel_mesh2d_set must contain the face nodes information
• To limit the search of the potential crossed faces, the m_facesRTree RTree with a search radius could be used.

[Huite]

For what it's worth, I've implemented this for the cell tree: https://deltares.github.io/numba_celltree/api.html#numba_celltree.CellTree2d.intersect_edges

There's another way of intersecting a mesh, which is by finding the intersection with the boundary, "entering" the cell and finding which internal edge it intersects, then entering the next cell through the edge connectivity. This obviously requires a full check with all entries and exists, otherwise you terminate too early when e.g. your mesh is a donut. However, the bigger problem is robustness: your line might intersect a cell vertex (almost) exactly, you end up in the wrong neighboring cell, or on the edge between two cells (and similarly with floating point issues). Searching a tree is more robust.

I don't think a search radius is the way. Your RTree has a bounding box for every node, and you can run a box-line intersection test as you descend the tree. This is a relatively cheap test (e.g. Cohen-Sutherland or Liang-Barsky if bounding boxes are axis aligned), then check the children; if the child is a leaf, you run an actual line-polygon test (e.g. Cyrus-Beck). This is what I do as well here. (The code for an RTree should be simpler than for a CellTree, as a CellTree has these kinda weird overlapping "buckets", whereas an RTree has ordinary bounding boxes.)

[lucacarniato]

@Huite, thank your useful suggestions. After thinking a bit more about the problem, I was indeed considering using the edge node connectivity to limit the number of edges to check for intersections.

As you mention, using a box-line intersection is also another option, which should be supported by the RTree implementation of boost.

[arthurvd]

@Huite : which line intersection support have you got in xugrid? The example on its frontpage (using uda.ugrid.sel) shows slicing through instructured data, which eventually comes from numba_celltree.CellTree2d.intersect_edges(), right?

Does that one have full support for an intersection request as sketched in Luca's above picture?

[Huite]

Hi @arthurvd, yes that's correct. The intersect_edges() method from the CellTree2d will return the indices which can be used to fetch the values of the faces.

Taking a good look, I think there's a mistake in the docstring of intersect_edges of the CellTree2d, it returns the xy coordinates rather than the length of the intersected edges.

I haven't tested the performance very exhaustively, but the method is based on the ray tracing method which is in vtk.