Point inside curved finite element

Question

I like to create interpolation functions for second order finite element meshes. For elements with straight edges all is good, but some of my elements may have curved edges as shown in the figure:

I am looking for references on how to test efficiently and robustly if a query point is inside a curved element. The interpolation itself is not a problem. One general way to do it is to find the mapping of the global element in $x$, $y$-coordinates to its mother element in $r$, $s$-coordinates, then map the query point and see if that mapped query point is within the mother element which has straight edges. This seems quite excessive and I am wondering if there are alternatives that you could point me to.

Update: This is for a general interpolation, to be used in plotting functions and the like. The query points can be 'random'. Even if one has a very good algorithm to find the closest element to a query point one needs to check if it is in the element. If that's not the case then a new potential element needs to be found and tested. This means that at a minimum at least one query needs to be done per query point. The expensive testing can be reduced if information if the element is curved or not can be used. In case the element is not curved then a cheap linear test can be used.

Answer 1

I don't think you can do better in general than mapping to the reference element and testing there. If your mapping is somehow special, you might be able to develop a test that's more efficient than the usual method, but I haven't seen one in practice.

One question I have for you is what are you doing that requires lots of checks for which element a point is in?

Answer 2

If your element is quadratic, you can find an implicit equation of the three quadratic edges of the triangle in function of $x,y$, in the form $E_1(x,y)=0$, $E_2(x,y)=0$, $E_3(x,y)=0$ (more references below about how to do that). The three function $E_1$, $E_2$ and $E_3$ measure an (oriented) distance to the curved edges, and the points inside the quadratic triangle are characterized by $E_1(x,y) \ge 0$ and $E_2(x,y) \ge 0$ and $E_3(x,y) \ge 0$.

How to find the implicit form from the Bézier control mesh is explained in Section 4.4 of [1] (it is used in computer graphics to display quadratic curves on graphic processor). Now if you have a large number of quadratic triangles (a mesh) and want to find which one contains a particular point, the same paper proposes a space decomposition data structure that dramatically speeds-up things (Section 5).

There is another reference on the "implicitization" procedure [2].

Note that this technique is interesting if you have many point queries to do (else the preprocessing phase / transform may be prohibitive)

[1] http://research.microsoft.com/en-us/um/people/hoppe/proj/ravg/

[2] GOLDMAN R., SEDERBERG T., AND ANDERSON D. 1984. Vector elimination: A technique for the implicitization, inversion, and intersection of planar parametric rational polynomial curves. CAGD 1, 327-356.

Answer 3

You could adapt one of the traditional Point-in-Polygon algorithms to suit your needs. In particular, you could generalize the "Ray Casting" method. The idea is this: Cast a ray in any direction from your point. Compute the number of times that ray intersects the boundary of your element (note, this could be >1 per side for your curved edges). If the number of intersections is odd, then the point is inside the element. If the nubmer of intersections is even (or zero), the point is outside the element.

The cound is a bit sensitive to corner cases (literally corners: you have to be careful not to double-count an intersection at an element corner), but as long as you take care of that, it's not hard to write a reasonably efficient implementation of this. For the quadratic elements that you have pictured, you will be solving for the intersection between a line and a parabola, which is very easy to do (solve a quadratic equation). Once you have your 0 or 2 solutions to the equation, you check to see if the intersections occur on the segment of the parabola that comprises your element edge.

I've never done this for curved elements before, but I can't see any reason why the method shouldn't work. Generalizing this to 3D should work, in principle, but computing the intersections would probably be more difficult.

Here's a link to the wikipedia page that I used as a reference when I encountered this problem. https://en.wikipedia.org/wiki/Point_in_polygon#Ray_casting_algorithm