What is the value of f(a,b,c), at the probe's position ?

Analytically, this problem is trivial. Not so numerically! A Finite Element context is assumed in the sequel. Besides, it is assumed that the elements are triangles in the 2-D case and tetrahedrons in the 3-D case.

The major problem encountered is the following: inside which of the elements 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' a finite element.

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

For IO(h,k,l) < 0 a point (h,k,l) is

For IO(h,k,l) = 0 a point (h,k,l) is at the element boundary.

It is possible to consider more complicated forms like quadrilaterals or bricks (which certainly gives rise to very nice 2-D and 3-D visualizations of IO functions) but it was promised that we should restrict ourselves to these

Let's take a look at the IO function of a triangle. We have learned 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 , 1 - h - k ) .The maximum of the IO function is reached for h = k = 1 - h - k --> = 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. The contour lines of this function are triangles, where the contour line with height 0 is the original triangle itself. 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:

| x-x1 x3-x1 x4-x1 | | x2-x1 x-x1 x4-x1 | h = | y-y1 y3-y1 y4-y1 | / DET ; k = | y2-y1 y-y1 y4-y1 | / DET ; | z-z1 z3-z1 z4-z1 | | z2-z1 z-z1 z4-z1 | | x2-x1 x3-x1 x-x1 | | x2-x1 x3-x1 x4-x1 | l = | y2-y1 y3-y1 y-y1 | / DET ; DET = | y2-y1 y3-y1 y4-y1 | ; | z2-z1 z3-z1 z-z1 | | z2-z1 z3-z1 z4-z1 | N(1) = h ; N(2) = k ; N(3) = l ; N(4) = 1-h-k-l . IO = min( h , k , l , 1-h-k-l ) .The maximum of the IO function is reached for h = k = l = 1-h-k-l --> = 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 . The isosurfaces of this function are tetrahedra, where the isosurface with height 0 is the original tetrahedron itself. Hope you get the 4-D picture ;-)

In fact we have an *inverse* problem. The connectivity array determines
the vertex points belonging to a given triangle. What we need here instead 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 ('duale
samenhang' in **Dutch**). This is also equivalent with saying that
Voronoi regions of a Delaunay triangulation are to be constructed.

In the 3-D 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 3-D 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 handled. It
may happen that the search path ends before a suitable element has been found.
This will repeatably be the case with non-convex 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 will
*never be less efficient* than this. Most of the time, the overall gain
in efficiency will be considerable. In the 3-D case, say with 1000 elements in
a convex region, it may be possible to localize the probe point within
approximately 10 steps.

Even *real-time* performance should be within the realm of posibilities:
using the current tetrahedron as the one to begin with, the next tetrahedron to
be probed can be found almost *immediately* by the above algorithm.

fc Random.FOR a.out # Outputs: 't' , 'coor2d.dat' , 'coor3d.dat' voronoi -t < t > topol2d.dat tetras < coor2d.dat > topol3d.datAfter this, you may try to put the programs at work. For example by typing:

fc suna57a.FOR -o probe2d probe2d 0.5 0.5 fc suna57b.FOR -o probe3d probe3d 0.5 0.4 0.6Here the numbers '0.5 0.7' respectively '0.5 0.4 0.6' stand for the ( x , y ) or ( x , y , z ) coordinates of a 2-D respectively 3-D (probe-)point.

Visualization of the search process is feasible for the 2-D case. A little post-processing program was developed for this purpose. At first, standard output has to be re-directed to a file named 'found2d.dat'. Then the program must be loaded, together with a PostScript initialization routine, a drawing routine and a routine for minimizing numerical output. At last the whole thing is executed and the result may be enjoyed at the screen and/or spooled to a PostScript printer. (Search path is painted gray, target triangle is painted black. And I'm too lazy to debug the random gray lines ... ;-)

probe2d 0.5 0.5 > found2d.dat fc Letsee.FOR plots.FOR plot.FOR log10.FOR -o Letsee Letsee gs meshpath.ps lpr -Pps meshpath.psO yeah, the programs contain a number of Fortran 'parameter' statements, by which the size of arrays and the maximum length of files is determined. Their values in the delivery, here, are quite arbitrary:

Random.FOR: parameter (LIMIT=20) suna57b.FOR: parameter (LIMIT=25,MANY=75) suna57a.FOR: parameter (LIMIT=100,MANY=187) Letsee.FOR: parameter (NODES=1500,LINES=4500)At last, all Fortran programs in this WebPage, and others, have been zipped into one (16 KB) file for convenience.