# Point inside curved finite element

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.

• @nicoguaro, thanks for the corrections. There used to be a spellchecker that underlined bad spellings but it seems to be switched off by default and I can not find a way to switch it on. – user21 Jul 20 '15 at 17:55

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?

• Thanks for your time, for your question please see the update. – user21 Jul 21 '15 at 9:26
• I have accepted this for the reason that using mapping makes the code easier. What I could same with geometric transforms was small with respect to other parts of the process. – user21 Sep 11 '15 at 7:03

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)

[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.

• Thanks, these are the kind of references I am looking for! – user21 Jul 22 '15 at 15:44

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

• Is this more efficient than inverting the map to the reference element and doing the standard test for whether a point is inside a straight-sided triangle? – Bill Barth Jul 20 '15 at 18:47
• For a 2nd-order or higher element like this one, can you invert the mapping easily? Obviously for an affine mapping, you could, but I'm not sure how you could do it in the general case. – Tyler Olsen Jul 20 '15 at 18:49
• I see how you could do the inverse mapping with a Newton-Raphson iteration, or the like, but I'm having trouble finding an analytical way to do it. – Tyler Olsen Jul 20 '15 at 18:53
• it should converge in 1 to 2 iterations, do you really need a closed form? – Bill Barth Jul 20 '15 at 19:08
• Not really. I'd estimate that the Newton iterations would cost about the same as the ray-casting though. I was mostly just curious if there was a closed form out there. – Tyler Olsen Jul 20 '15 at 19:09