Need an algorithm that searches the cell that an arbitrarily given point is within
 

Such an algorithm is common in commercial CFD post-processing software. But I want to implement in my own CFD post-processing code:

In finite volume or finite element mesh:

[1] Specify an arbitrary point in space
[2] Find a cell or element that the given point can be within
[3] So that the solution value at this point can be interpolated by solution data.

Thanks in advance for recommending any literature for such an algorithm

This might be a good start.

http://graphics.cs.ucdavis.edu/~joy/...turedGrids.pdf

Many times ago I developed a simple method on unstructured grid.
It is based on the linear shape function on triangle. Whena point is outside the chosen triangle, you get alway a negative value for the shape function avaluated using it. Only when the point is inside a triangle, you get positive values.
I used for Lagrangian transport, the method then searches at next time step only adjacent triangles to save computational time.

I'll give my idea for a 2D triangular mesh (I don't know if it is documented, and I've never tried it myself)

For a given cell and a given point.
1- Calculate the area of the cell (or have it already calculated in a database)
2- Calculate the area of the triangle formed by every edge of the cell and the given point (this gives 3 area)

If the sum of the 3 area is greater than the area of the cell, the point is outside the domain.

If the point is on an edge from the considered cell, one of the triangle (from step 2) will have no area (3 co-linear points).

If the point is on an edge vertex of the considered cell, two of the triangle (from step 2) will have no area (3 co-linear points).

If the sum of areas is equal to the cell's area and all triangle have positive area, the given point is inside the given cell.

Cycle through each cell to find your answer. This might be inefficient if it is done randomly, but if you have neighbors cell in your mesh database, and have a good first guess, it might be pretty quick.

This idea can be carried over to 3D using volumes and faces.

Thanks for your concise explanation of the algorithm! I totally understand it now.

yes, this is exactly the method of linear shape function, they define nothing else that the normalized areas defined by a point in a triangle ;)

I would not describe it as 'THE' algorithm. It seems to me that it is 'one of the many' solutions. It is probably not the most efficient one, probably far from it in fact, but it has the large advantage to be pretty straight forward to program.

Usually, for convex polyhedral(polygonal) cells in 3D(2D), a much faster test is the triangle(segment)/segment intersection test between:

1) The segment connecting the given point and the cell centroid and

2) A triangular face of the cell (a face segment in 2D). For polyhedral cells, each face of the cell is usually decomposed in triangles before doing the test (this is usually needed in any case).

Looping on the faces of a cell, if no intersection is detected then the point is inside the cell, otherwise you have a pretty good guess for which is the next cell to test (the one from the other side of the intersected face).

You might want to find the first cell to test by first putting the cell centroids in a smart data structure like bins (i.e., a structured uniform grid) or an octree.

You might also want to implement some smart heuristics so that, during the search (when your first guess is wrong), you do not end up on a boundary with no more guesses (besides the one 'the point is outside the domain').

I also came across this method in literature search. But thanks a lot for explaining it so clearly here.

Is there an example algorithm for this part? (e.g., assuming we are dealing with tetrahedron). Thanks


Looping on the faces of a cell, if no intersection is detected then the point is inside the cell, otherwise you have a pretty good guess for which is the next cell to test (the one from the other side of the intersected face).

Here is a note I put together for testing the face intersection:

https://drive.google.com/file/d/0B4F...ew?usp=sharing

This was in the context of looking for particles escaping a cell, but it is the same idea. The test line segment would extend from the current cell centroid to the target point. As others have mentioned, if it intersects the face, move to the cell that is on the other side of that face. And just keep walking across the mesh until you find a cell whose centroid-to-test point line segment doesn't intersect a face. Then you know you are done.

PS. I should note that this algorithm works for general planar faces directly. There is no need to decompose quad/polygon faces into triangles. I haven't seen this algorithm published anywhere, but that is probably only because I haven't looked very hard.

You can find attached some routines i wrote some time ago. They are part of a larger piece of code, so some parts might not look clear (i.e., the data structures), but the logic
is there.

The routines are trivial and with no error handling (mostly because the grid is checked before executing them), so take them just as inspiration (especially the intersection tests).

Also, no attempt is made to recover from an erroneous path leading to the mesh boundary and the wrong answer. For example, if at a convex mesh corner there are cells with large differences in sizes, the algorithm, as it is, will probably fail.