Check if a point is inside an arbitrary hexahedron

user_8044060 0
Viewed 178

I am working on a 3D finite element code, where i face the following problem:

If I take an arbitrary point (say x), how do I figure out, which element it is in?

This can be simplified to: How do I check if an arbitrary point (x) lies inside or outside of an (hexahedral) element?

What I already found:

Contrary to the two approaches above, my problem does not assume right angles nor parallel faces.

Problem sketch:

Notation: (again: though the sketch shows a regular shape, our hexahedron is assumed to be of general shape)

  • 8-node hexahedron topology, nodes: 0,..,7

  • axis: r,s,t

                      t
                  |
         4--------|-------------7
        /|        |            /|
       / |        |           / |
      /  |        |          /  |
     /   |        |         /   |
    /    |        |        /    |
   /     |        |       /     |
  5----------------------6      |
  |      |        |      |      |
  |      |        o------|---------s
  |      |       /       |      |
  |      0------/--------|------3
  |     /      /         |     /
  |    /      /          |    /
  |   /      /           |   /
  |  /      /            |  /
  | /      r             | /
  |/                     |/
  1----------------------2

Data that we have available:

  • coordinates of the nodes (vectors P0 to P7)
  • coordinates of the point we want to check (lets say Px)

Additionaly we assume the nodes are ordered as sketched above.

My approach/solution so far:

  1. calculate the surface (outward) normal vectors

    Use cross products, eg. for the r_pos_normal_vec (pointing out of the plane)

    r_pos_normvec = (P2-P1) x (P5-P1)

    and for the r_neg_normal_vec

    r_neg_normvec = (P4-P0) x (P3-P0)

    similarly for the s and t directions

  2. check two opposite corner nodes (I chose node0 and node 6)

    • For node0

      • calculate vector from P0 to Px:

        P0x = Px - P0

      • calculate inner prodcut of P0x and surfaces adjacent to node 0

        <P0x, r_neg_normal_vec>

        <P0x, s_neg_normal_vec>

        <P0x, t_neg_normal_vec>

    • For node1

      • same scheme as for node 0, whereas P1 instead of P0 and the positive counterparts of the normal vectors are used

  3. Iff all 6 (3 from node0 and 3 from node1) inner products result in negative value -> the point is inside the hexahedron.

Question:

I implemented the functionality described above in my code and ran some tests. It seems to work, from the math side I am quite confident.

Please discuss my approach, I am happy for any hints/clues/recommendations/bug fixes ... Is there some way to make this faster? Alternative solutions?

Note:

  • To speed up the algorithm a box check can be done first:
    • Construct a rectangular box around the hexahedron:
    • Get the min and max values of the node coordinates in each direction.
    • If the point to check (x) is outside this (larger) box, it cannot be inside the hexahedron.
mesh computational-geometry geometry
3 Answers

For any convex polyhedron, establish the implicit equations of the faces (f.i. plane by three points), of the form ax+by+cz+d=0.

When you plug the coordinates of a known point inside the volume (such as the center) in the expression ax+by+cz+d, you will get a set of signs. An arbitrary point is inside if it yields the same signs.

Update:

For maximum performance, you can consider also using an axis-aligned bounding box for quick rejection. This only makes sense if many of the points are outside. Make sur to use a shortcut evaluation so that early rejection can happen.

Note that a rejection test such as X<Xmin is nothing but the above sign test against the plane of equation X-Xmin=0.

user_11965490

I personnally prefer your method, however there also is a way to approach the problem if the hexahedral restricted to parallelepiped. So you can transfer the coordinate of P in the frame $(0; e_1; e_2; e_3)$ to $(P_0, P_0P_1,P_0P_3,P_0P_4)$. We call it $(a,b,c)$, then the point is in that parallelepiped if $a,b,c > 0 a \in [0,1], a+b+c \in [0,1]$.

user_174890620

Because you mentioned that you want to be able to handle arbitrary hexahedrons, I think that your process might be improved if you adjust your s, r, and t normals to account for having slightly warped faces. I would do this by making the following change to r normals (and similar for s and t):

r_pos_normvec = (P6-P1) x (P5-P2)

r_neg_normvec = (P7-P0) x (P4-P3)

This would be important for a case where you shifted node 6 towards node 7 (say 0.9xP6) and had a point at 0.95xP6. Without the warping correction, I believe you would erroneously determine the point as inside the hexahedron.

user_143761960
Related