sci.math.numanalysis
SUNA57, Efficient 2D & 3D Point Probes
========================================
The shape functions of triangles and tetrahedra can be used for solving quite
other problems than systems of discretized partial differential equations, by
means of a Finite Element Method. The problem we have at hand in this article
can be stated as follows:
Given a point ("probe") with coordinates (a,b,c). Given a (discretized) field
where function values f(x,y,z) are known (at nodal point positions). What is
the value of f(a,b,c), at the probe ?
Analytically, this problem is trivial. Not so numerically. We assume a Finite
Element context in the sequel. And we assume that the elements are triangles
in the 2D case, and tetrahedrons in the 3D case.
The major problem encountered is the following: inside which element the point
probe has to be located? The brute force approximation to this is: simply scan
all elements, and look if the probe point is "inside" or "outside" an element.
Anyway, we first must have a mathematical precise criterion for effectuating
the common sense notions "inside" and "outside". It is here that the element
shape functions turn out to be of much help.
Criterion

 Let the shape functions of a (say linear) element be called N(i) . Then an
 Inside/Outside function is associated with the element, which is defined as:

 IO(h,k,l) = minimum of N(i)(h,k,l) ; here (h,k,l) are local coordinates.

 For IO(h,k,l) > 0 a point (h,k,l) is _inside_ the element;
 For IO(h,k,l) < 0 a point (h,k,l) is _outside_ the element;
 For IO(h,k,l) = 0 a point (h,k,l) is at the element boundary.
It is possible to consider complicated cases like quadrilaterals or bricks (and
it certainly gives rise to nice visualizations of IO functions) but I promised
that we should restrict ourselves to those linear elements.
Let's take a look at the IO function of a triangle. From SUNA02 we learn that
it's shape functions are given by:
 x  x1 x3  x1   x2  x1 x  x1 
h =   / DET ; k =   / DET ;
 y  y1 y3  y1   y2  y1 y  y1 
 x2  x1 x3  x1 
DET =   .
 y2  y1 y3  y1 
N(1) = h ; N(2) = k ; N(3) = 1  h  k ; IO = min(h,k,1hk) .
The maximum of the IO function is reached for h = k = 1hk > = 1/3 ,
hence at the midpoint (barycenter) of the triangle. If we draw straight lines
from the the midpoint towards the vertices, and further, then the whole plane
is subdvided into three regions, one where IO = h , one where IO = k and one
where IO = 1  h  k . The IO function is zero at the triangle sides, positive
inside and negative outside. It's shaped like a mountain with top 1/3 at the
midpoint and three sharply edged slopes downhill. Hope you get the picture ...
Let's take a look at the IO function of a tetrahedron now. It's shape functions
are given by:
 xx1 x3x1 x4x1   x2x1 xx1 x4x1 
h =  yy1 y3y1 y4y1  / DET ; k =  y2y1 yy1 y4y1  / DET ;
 zz1 z3z1 z4z1   z2z1 zz1 z4z1 
 x2x1 x3x1 xx1   x2x1 x3x1 x4x1 
l =  y2y1 y3y1 yy1  / DET ; DET =  y2y1 y3y1 y4y1  ;
 z2z1 z3z1 zz1   z2z1 z3z1 z4z1 
N(1) = h ; N(2) = k ; N(3) = l ; N(4) = 1  h  k  l .
IO = min( h, k, l, 1hk) .
The maximum of the IO function is reached for h = k = l = 1hkl > = 1/4 ,
hence at the midpoint (barycenter) of the tetrahedron. If we nw draw straight
lines from the the midpoint towards the vertices, and further, then the whole
space is subdvided into four regions, one where IO = h , one where IO = k , one
with IO = l , and one where IO = 1  h  k  l . Hope you get the 4D picture.
Back to the 2D case. Instead of the brute force approach, that is simply test
if IO(a,b) >= 0 for each triangle, a much more sophisticated search procedure
now comes into mind. Determine the slope where the probe point is "sliding" on.
If for example IO = h , which is the shape function belonging to vertex (1),
then the probe must be searched in the direction (2)(3). In general, if the IO
function equals the shape function belonging to vertex (k), then searching must
proceed in the direction of the other two vertices. Stated more explicitly: the
next triangle to be investigated is the one that has these two other vertices
in common with the present triangle. This in turn means that we should be able
to lookup quickly vertex points in the connectivity array, preferrably without
having to transverse the whole array. In fact we have an inverse problem. The
connectivity array determines the vertex points belonging to a given triangle.
What we need here is to know which triangles belong to a given vertex point. A
general solution to this problem consists in constructing the inverse topology,
also called "dual" connectivity, or construction of the Voronoi regions. If we
keep a counter in the first element of the dual array, then the accompanying
Fortran program fragment reads as follows:
do ne=1,number_of_triangles
do node=1,3 ! node=1,4 for tetrahedrons
i=topology(node,ne)
dual(i,1)=dual(i,1)+1
counter=dual(i,1)
dual(i,counter+1)=ne
enddo ; enddo
In the 3D case. the situation is completely analogous. If the IO function is
equal to the shape function belonging to vertex (k), then the next tetrahedron
to be investigated is the one that has the other three vertices in common with
the present tetrahedron. A more elegant procedure for actually finding the next
element has been implemented in the 3D program. Three "downhill" vertices of
the neighbouring tetrahedron are known. Now the dual connectivity is scanned
for all tetrahedra which contain these vertices. A scratch array called MERGE
is defined and emptied. Every time an element NE is encountered, it's counter
in MERGE(NE) is incremented. At last, there will be two tetrahedra in MERGE
which have a count equal to 3 : the current element, and the neighbouring one,
to be searched next. The neighbouring element is then selected as "current" and
the search is continued, until the IO function for the probe point has attained
a positive value.
An additional array DONE keeps track of the elements which have been done. It
may happen that the search path ends before a suitable element has been found.
This will repeatably be the case with nonconvex area's / volumes. In that case
the search procedure will be restarted with an element which has not been DONE
already. The search procedure is restarted every time an element is encountered
which has been DONE : it has no sense to continue a deterministic process, the
outcome of which is known. This may happen until all elements are exhausted. At
last, the probe point may be outside the mesh. So our efficient probe algorithm
may turn out to be equivalent to the brute force method, in the end. But it is
_never less_ efficient than that. Most of the time, the gain in efficiency will
be considerable. In the 3D case, say with 1000 elements, in a convex region,
it may be possible to localize the probe point within approximately 10 steps.
Three Fortran programs were developed, meant as an implementation of the above.
 suna57a.for: Efficient 2D probe. Works with data "coor2d.dat", "topol2d.dat"
 suna57b.for: Efficient 3D probe. Works with data "coor3d.dat", "topol3d.dat"
 suna57a.for: 3D probe for bricks. Works with output from "suna56" programs.

* Han de Bruijn; Applications&Graphics  "A little bit of Physics * No
* TUD Computing Centre; P.O. Box 354  would be NO idleness in * Oil
* 2600 AJ Delft; The Netherlands.  Mathematics" (HdB). * for
* Email: Han.deBruijn@RC.TUDelft.NL  Fax: +31 15 78 37 87 * Blood