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

 Register Blogs Members List Search Today's Posts Mark Forums Read

 July 14, 2015, 19:46 Need an algorithm that searches the cell that an arbitrarily given point is within #1 New Member   Yidong Join Date: Nov 2011 Posts: 23 Rep Power: 6 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

 July 14, 2015, 22:25 #2 Member   Join Date: Jul 2013 Posts: 49 Rep Power: 5 This might be a good start. http://graphics.cs.ucdavis.edu/~joy/...turedGrids.pdf

 July 15, 2015, 03:23 #3 Senior Member   Filippo Maria Denaro Join Date: Jul 2010 Posts: 2,411 Rep Power: 31 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.

 July 15, 2015, 09:17 #4 Member   Join Date: Jul 2013 Posts: 49 Rep Power: 5 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.

July 15, 2015, 11:02
#5
New Member

Yidong
Join Date: Nov 2011
Posts: 23
Rep Power: 6
Thanks for your concise explanation of the algorithm! I totally understand it now.

Quote:
 Originally Posted by Alex C. 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.

July 15, 2015, 11:25
#6
Senior Member

Filippo Maria Denaro
Join Date: Jul 2010
Posts: 2,411
Rep Power: 31
Quote:
 Originally Posted by Alex C. 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.

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

July 15, 2015, 11:27
#7
Member

Join Date: Jul 2013
Posts: 49
Rep Power: 5
Quote:
 Originally Posted by yidongxia Thanks for your concise explanation of the algorithm! I totally understand it now.
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.

 July 16, 2015, 03:38 #8 Senior Member     Paolo Lampitella Join Date: Mar 2009 Location: Italy Posts: 552 Blog Entries: 14 Rep Power: 18 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').

July 16, 2015, 13:12
#9
New Member

Yidong
Join Date: Nov 2011
Posts: 23
Rep Power: 6
I also came across this method in literature search. But thanks a lot for explaining it so clearly here.

Quote:
 Originally Posted by sbaffini 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').

July 16, 2015, 14:36
#10
New Member

Yidong
Join Date: Nov 2011
Posts: 23
Rep Power: 6
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).

Quote:
 Originally Posted by sbaffini 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').

July 16, 2015, 18:30
#11
Senior Member

Michael Prinkey
Join Date: Mar 2009
Location: Pittsburgh PA
Posts: 257
Rep Power: 15
Quote:
 Originally Posted by yidongxia 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:

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.

July 17, 2015, 05:39
#12
Senior Member

Paolo Lampitella
Join Date: Mar 2009
Location: Italy
Posts: 552
Blog Entries: 14
Rep Power: 18
Quote:
 Originally Posted by yidongxia Is there an example algorithm for this part? (e.g., assuming we are dealing with tetrahedron). Thanks
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.
Attached Files
 findpointcell.txt (1.3 KB, 6 views) ispointincell.txt (2.9 KB, 3 views) segint2d.txt (727 Bytes, 2 views) segtriint3d.txt (728 Bytes, 1 views)

 Thread Tools Display Modes Linear Mode

 Posting Rules You may not post new threads You may not post replies You may not post attachments You may not edit your posts BB code is On Smilies are On [IMG] code is On HTML code is OffTrackbacks are On Pingbacks are On Refbacks are On Forum Rules

 Similar Threads Thread Thread Starter Forum Replies Last Post lakeat OpenFOAM Mesh Utilities 57 February 1, 2015 09:25 chenxizh ParaView 4 October 28, 2013 08:44 shrina OpenFOAM Running, Solving & CFD 10 October 3, 2013 14:38 Purushothama CD-adapco 2 May 31, 2010 21:58 Ben Main CFD Forum 5 June 22, 2007 11:23

All times are GMT -4. The time now is 09:00.