Finite Element Analysis
Truth is the offspring of silence and unbroken meditation.Isaac Newton

This program is a demonstration of finite element analysis (FEA)—also referred to as the finite element method (FEM)—with a two dimensional argument as input. As a test case, program FEA reads files defining a two dimensional FEM representation (nodes + elements + values) of the surface $f(x, y) = x^2 + y^2$, and provides a file containing the values of the finite element function at the sample nodes.

Start by downloading the archive that contains the Fortran source code and all the support files needed for the project. Compile the source with:

joel@DeWitt:~$gfortran -g -std=f2008 -Wall -Wextra -O2 -fall-intrinsics FEA.f08 on a Unix/Linux system. The archive contains FEM data for$f(x, y)$on a 5×5 grid of nodes, organized into linear triangles. The sample data seeks the values of this function on a 4×4 grid. The files we'll make use of are fem_sq_nodes.txt, which contains the node coordinates; fem_sq_elements.txt, which contains the nodes that make up each element; fem_sq_values.txt, which contains the values defined at each node. Additionally, we'll need the node coordinates where samples are desired, which are contained in the file sample_sq_nodes.txt. (Note that the _nodes.txt files essentially contain arbitrary$x$-$y$coordinates that are of interest to us.) Having sorted all this out, let's now run the demonstration by entering: joel@DeWitt:~$ ./a.out fem_sq sample_sq

This will produce the file sample_sq_values.txt, which contains the values computed at each sample node. This is basically interpolated FEM data. Let's now proceed by examining what each of the (interesting) subroutines do within the program FEA.

## Apply the Finite Element Method to a Surface

### SUBROUTINE basis_mn_t3

SUBROUTINE basis_mn_t3(t, n, p, phi, dphidx, dphidy)
IMPLICIT NONE
INTEGER(kind=4) :: n
REAL(kind=8) :: area, dphidx(3, n), dphidy(3, n), p(2, n), phi(3, n), t(2, 3)
area = t(1, 1) * (t(2, 2) - t(2, 3)) &
+ t(1, 2) * (t(2, 3) - t(2, 1)) &
+ t(1, 3) * (t(2, 1) - t(2, 2))
IF (INT(area) == 0) THEN
WRITE(*, '(a)') " "
WRITE(*, '(a)') "BASIS_MN_T3 - Fatal error!"
WRITE(*, '(a)') "Element has zero area."
STOP 1
END IF
phi(1,1:n) =     (   ( t(1,3) - t(1,2) ) * ( p(2,1:n) - t(2,2) )     &
- ( t(2,3) - t(2,2) ) * ( p(1,1:n) - t(1,2) ) )
dphidx(1,1:n) =    - ( t(2,3) - t(2,2) )
dphidy(1,1:n) =      ( t(1,3) - t(1,2) )
phi(2,1:n) =     (   ( t(1,1) - t(1,3) ) * ( p(2,1:n) - t(2,3) )     &
- ( t(2,1) - t(2,3) ) * ( p(1,1:n) - t(1,3) ) )
dphidx(2,1:n) =    - ( t(2,1) - t(2,3) )
dphidy(2,1:n) =      ( t(1,1) - t(1,3) )
phi(3,1:n) =     (   ( t(1,2) - t(1,1) ) * ( p(2,1:n) - t(2,1) )     &
- ( t(2,2) - t(2,1) ) * ( p(1,1:n) - t(1,1) ) )
dphidx(3,1:n) =    - ( t(2,2) - t(2,1) )
dphidy(3,1:n) =      ( t(1,2) - t(1,1) )
phi(1:3,1:n) = phi(1:3,1:n) / area !normalize
dphidx(1:3,1:n) = dphidx(1:3,1:n) / area
dphidy(1:3,1:n) = dphidy(1:3,1:n) / area
END SUBROUTINE basis_mn_t3


Given the coordinates of the vertices of a triangle, or T3 element, this routine evaluates the basis functions [PDF] associated with each vertex, and their derivatives with respect to $x$ and $y$. It works directly with these coordinates, and does not refer to a reference element. The sides of the triangle do not have to lie along a coordinate axis.

The two dimensional input array t(2, 3) contains the double precision coordinates of the vertices of the triangle. It is common to list these points in counterclockwise order. The 4-byte integer n is also an input, and is equal to the number of evaluation points. The next input is the two dimensional, double precision array p(2, n), which contains the points where the basis functions are to be evaluated. On the output side of this subroutine, the array phi(3, n) holds the values of the basis functions at the evaluation points. Moreover, the output arrays dphidx(3, n) and dphidy(3, n) contain the values of the derivatives at the evaluation points. Finally, the local, double precision variable area is equal to (twice) the area of the triangle.

### SUBROUTINE basis_mn_t6

SUBROUTINE basis_mn_t6(t, n, p, phi, dphidx, dphidy)
IMPLICIT NONE
INTEGER(kind=4) :: n
REAL(kind=8) :: dphidx(6, n), dphidy(6, n), gn(n), gx(n), hn(n), hx(n), p(2, n), phi(6, n), t(2, 6)
!basis function 1: phi(x, y) = g(3, 2) * h(6, 4) / normalization
gx(1:n) = ( p(1,1:n) - t(1,2) ) * ( t(2,3)   - t(2,2) ) &
- ( t(1,3)   - t(1,2) ) * ( p(2,1:n) - t(2,2) )
gn(1:n) = ( t(1,1)   - t(1,2) ) * ( t(2,3)   - t(2,2) ) &
- ( t(1,3)   - t(1,2) ) * ( t(2,1)   - t(2,2) )
hx(1:n) = ( p(1,1:n) - t(1,4) ) * ( t(2,6)   - t(2,4) ) &
- ( t(1,6)   - t(1,4) ) * ( p(2,1:n) - t(2,4) )
hn(1:n) = ( t(1,1)   - t(1,4) ) * ( t(2,6)   - t(2,4) ) &
- ( t(1,6)   - t(1,4) ) * ( t(2,1)   - t(2,4) )
phi(1,1:n) =     ( gx(1:n) * hx(1:n) ) / ( gn(1:n) * hn(1:n) )
dphidx(1,1:n) =  (      ( t(2,3) - t(2,2) ) * hx(1:n) &
+ gx(1:n) * ( t(2,6) - t(2,4) ) ) / ( gn(1:n) * hn(1:n) )
dphidy(1,1:n) = -(      ( t(1,3) - t(1,2) ) * hx(1:n) &
+ gx(1:n) * ( t(1,6) - t(1,4) ) ) / ( gn(1:n) * hn(1:n) )
!basis function 2: phi(x, y) = g(3, 1) * h(4, 5) / normalization
gx(1:n) = ( p(1,1:n) - t(1,1) ) * ( t(2,3)   - t(2,1) ) &
- ( t(1,3)   - t(1,1) ) * ( p(2,1:n) - t(2,1) )
gn(1:n) = ( t(1,2)   - t(1,1) ) * ( t(2,3)   - t(2,1) ) &
- ( t(1,3)   - t(1,1) ) * ( t(2,2)   - t(2,1) )
hx(1:n) = ( p(1,1:n) - t(1,5) ) * ( t(2,4)   - t(2,5) ) &
- ( t(1,4)   - t(1,5) ) * ( p(2,1:n) - t(2,5) )
hn(1:n) = ( t(1,2)   - t(1,5) ) * ( t(2,4)   - t(2,5) ) &
- ( t(1,4)   - t(1,5) ) * ( t(2,2)   - t(2,5) )
phi(2,1:n) = ( gx(1:n) * hx(1:n) ) / ( gn(1:n) * hn(1:n) )
dphidx(2,1:n) =  (      ( t(2,3) - t(2,1) ) * hx(1:n) &
+ gx(1:n) * ( t(2,4) - t(2,5) ) ) / ( gn(1:n) * hn(1:n) )
dphidy(2,1:n) = -(      ( t(1,3) - t(1,1) ) * hx(1:n) &
+ gx(1:n) * ( t(1,4) - t(1,5) ) ) / ( gn(1:n) * hn(1:n) )
!basis function 3: phi(x, y) = g(1, 2) * h(5, 6) / normalization
gx(1:n) = ( p(1,1:n) - t(1,2) ) * ( t(2,1)   - t(2,2) ) &
- ( t(1,1)   - t(1,2) ) * ( p(2,1:n) - t(2,2) )
gn(1:n) = ( t(1,3)   - t(1,2) ) * ( t(2,1)   - t(2,2) ) &
- ( t(1,1)   - t(1,2) ) * ( t(2,3)   - t(2,2) )
hx(1:n) = ( p(1,1:n) - t(1,6) ) * ( t(2,5)   - t(2,6) ) &
- ( t(1,5)   - t(1,6) ) * ( p(2,1:n) - t(2,6) )
hn(1:n) = ( t(1,3)   - t(1,6) ) * ( t(2,5)   - t(2,6) ) &
- ( t(1,5)   - t(1,6) ) * ( t(2,3)   - t(2,6) )
phi(3,1:n) = ( gx(1:n) * hx(1:n) ) / ( gn(1:n) * hn(1:n) )
dphidx(3,1:n) =  (      ( t(2,1) - t(2,2) ) * hx(1:n) &
+ gx(1:n) * ( t(2,5) - t(2,6) ) ) / ( gn(1:n) * hn(1:n) )
dphidy(3,1:n) = -(      ( t(1,1) - t(1,2) ) * hx(1:n) &
+ gx(1:n) * ( t(1,5) - t(1,6) ) ) / ( gn(1:n) * hn(1:n) )
!basis function 4: phi(x, y) = g(1, 3) * h(2, 3) / normalization
gx(1:n) = ( p(1,1:n) - t(1,3) ) * ( t(2,1)   - t(2,3) ) &
- ( t(1,1)   - t(1,3) ) * ( p(2,1:n) - t(2,3) )
gn(1:n) = ( t(1,4)   - t(1,3) ) * ( t(2,1)   - t(2,3) ) &
- ( t(1,1)   - t(1,3) ) * ( t(2,4)   - t(2,3) )
hx(1:n) = ( p(1,1:n) - t(1,3) ) * ( t(2,2)   - t(2,3) ) &
- ( t(1,2)   - t(1,3) ) * ( p(2,1:n) - t(2,3) )
hn(1:n) = ( t(1,4)   - t(1,3) ) * ( t(2,2)   - t(2,3) ) &
- ( t(1,2)   - t(1,3) ) * ( t(2,4)   - t(2,3) )
phi(4,1:n) = ( gx(1:n) * hx(1:n) ) / ( gn(1:n) * hn(1:n) )
dphidx(4,1:n) =  (      ( t(2,1) - t(2,3) ) * hx(1:n) &
+ gx(1:n) * ( t(2,2) - t(2,3) ) ) / ( gn(1:n) * hn(1:n) )
dphidy(4,1:n) = -(      ( t(1,1) - t(1,3) ) * hx(1:n) &
+ gx(1:n) * ( t(1,2) - t(1,3) ) ) / ( gn(1:n) * hn(1:n) )
!basis function 5: phi(x, y) = g(2, 1) * h(3, 1) / normalization
gx(1:n) = ( p(1,1:n) - t(1,1) ) * ( t(2,2)   - t(2,1) ) &
- ( t(1,2)   - t(1,1) ) * ( p(2,1:n) - t(2,1) )
gn(1:n) = ( t(1,5)   - t(1,1) ) * ( t(2,2)   - t(2,1) ) &
- ( t(1,2)   - t(1,1) ) * ( t(2,5)   - t(2,1) )
hx(1:n) = ( p(1,1:n) - t(1,1) ) * ( t(2,3)   - t(2,1) ) &
- ( t(1,3)   - t(1,1) ) * ( p(2,1:n) - t(2,1) )
hn(1:n) = ( t(1,5)   - t(1,1) ) * ( t(2,3)   - t(2,1) ) &
- ( t(1,3)   - t(1,1) ) * ( t(2,5)   - t(2,1) )
phi(5,1:n) = ( gx(1:n) * hx(1:n) ) / ( gn(1:n) * hn(1:n) )
dphidx(5,1:n) =  (      ( t(2,2) - t(2,1) ) * hx(1:n) &
+ gx(1:n) * ( t(2,3) - t(2,1) ) ) / ( gn(1:n) * hn(1:n) )
dphidy(5,1:n) = -(      ( t(1,2) - t(1,1) ) * hx(1:n) &
+ gx(1:n) * ( t(1,3) - t(1,1) ) ) / ( gn(1:n) * hn(1:n) )
!basis function 6: phi(x, y) = g(1, 2) * h(3, 2) / normalization
gx(1:n) = ( p(1,1:n) - t(1,2) ) * ( t(2,1)   - t(2,2) ) &
- ( t(1,1)   - t(1,2) ) * ( p(2,1:n) - t(2,2) )
gn(1:n) = ( t(1,6)   - t(1,2) ) * ( t(2,1)   - t(2,2) ) &
- ( t(1,1)   - t(1,2) ) * ( t(2,6)   - t(2,2) )
hx(1:n) = ( p(1,1:n) - t(1,2) ) * ( t(2,3)   - t(2,2) ) &
- ( t(1,3)   - t(1,2) ) * ( p(2,1:n) - t(2,2) )
hn(1:n) = ( t(1,6)   - t(1,2) ) * ( t(2,3)   - t(2,2) ) &
- ( t(1,3)   - t(1,2) ) * ( t(2,6)   - t(2,2) )
phi(6,1:n) = ( gx(1:n) * hx(1:n) ) / ( gn(1:n) * hn(1:n) )
dphidx(6,1:n) =  (      ( t(2,1) - t(2,2) ) * hx(1:n) &
+ gx(1:n) * ( t(2,3) - t(2,2) ) ) / ( gn(1:n) * hn(1:n) )
dphidy(6,1:n) = -(      ( t(1,1) - t(1,2) ) * hx(1:n) &
+ gx(1:n) * ( t(1,3) - t(1,2) ) ) / ( gn(1:n) * hn(1:n) )
END SUBROUTINE basis_mn_t6


The basic structure of this routine is identical to that of SUBROUTINE basis_mn_t3, except that, in addition to the coordinates of the vertices of an input T6 element, the midside nodes are specified as well. Similar to the previous routine, it works directly with these coordinates, and does not refer to a reference element. This routine requires that the midside nodes be "in line" with the vertices; that is, the sides of the triangle must be straight. However, the midside nodes do not actually have to be halfway along the side of the triangle.

The input/output variables and arrays hold the same meanings as those found in SUBROUTINE basis_mn_t3.

### SUBROUTINE fem2d_evaluate

SUBROUTINE fem2d_evaluate(fem_node_num, fem_node_xy, fem_element_order, fem_element_num, fem_element_node, fem_element_neighbor,   &
fem_value_dim, fem_value, sample_node_num, sample_node_xy, sample_value)
IMPLICIT NONE
INTEGER(kind=4) :: fem_element_num, fem_element_order, fem_node_num, fem_value_dim, sample_node_num, edge,                       &
fem_element_neighbor(3, fem_element_num), fem_element_node(fem_element_order, fem_element_num), i, j, t, t_node(fem_element_order)
REAL(kind=8) :: b(fem_element_order), dbdx(fem_element_order), dbdy(fem_element_order), fem_node_xy(2, fem_node_num),            &
fem_value(fem_value_dim, fem_node_num), sample_node_xy(2, sample_node_num), sample_value(fem_value_dim, sample_node_num),        &
t_xy(2, fem_element_order)
REAL(kind=8), DIMENSION(2) :: p_xy
find_and_evaluate: DO j = 1, sample_node_num !for each sample point, find the element triangle that contains it, and evaluate the finite element function
p_xy(1:2) = sample_node_xy(1:2, j)
CALL triangulation_search_delaunay(fem_node_num, fem_node_xy, fem_element_order, fem_element_num, fem_element_node,            &
fem_element_neighbor, p_xy, t, edge) !find the element triangle that contains the point
t_node(1:fem_element_order) = fem_element_node(1:fem_element_order, t) !evaluate the finite element basis functions at the point in the triangle
t_xy(1:2, 1:fem_element_order) = fem_node_xy(1:2, t_node)
IF (fem_element_order == 3) THEN
CALL basis_mn_t3(t_xy, 1, p_xy, b, dbdx, dbdy)
ELSE IF (fem_element_order == 6) THEN
CALL basis_mn_t6(t_xy, 1, p_xy, b, dbdx, dbdy)
END IF
multiply: DO i = 1, fem_value_dim
sample_value(i, j) = DOT_PRODUCT(fem_value(i, t_node(1:fem_element_order)), b(1:fem_element_order)) !multiply by the finite element values to get the sample values
END DO multiply
END DO find_and_evaluate
END SUBROUTINE fem2d_evaluate


This routine samples an FEM function on a T3 or T6 triangulation. Note that the sample values returned are the true values of the underlying finite element function. In other words, they are not produced by constructing some other function that interpolates the data at the finite element nodes. Instead, each sampling node is located within one of the associated finite element triangles, and the finite element function is developed and evaluated there. This routine attempts to reproduce the finite element function corresponding to a set of nodal data.

This routine accepts a number of 4-byte integers and integer arrays as input: fem_node_num = number of nodes; fem_element_order = order of the elements, either 3 or 6; fem_element_num = number of triangles; fem_element_node(fem_element_order, fem_element_num) = nodes that make up each triangle; fem_element_neighbor(3, fem_element_num) = index of the neighboring triangle on each side, or $-1$ if there is no neighbor; fem_value_dim = "dimension" of the values; sample_node_num = number of sample nodes. This routine also accepts the following double precision arrays: fem_node_xy(2, fem_node_num) = coordinates of the nodes; fem_value(fem_value_dim, fem_node_num) = finite element coefficient values at each node; sample_node_xy(2, sample_node_num) = sample nodes. On the output side of this subroutine, the array sample_value(fem_value_dim, sample_node_num) holds the sampled values.

Operationally, this routine performs a number of interesting tasks. The DO loop labeled find_and_evaluate examines each sample point, finding the element triangle that contains it, and evaluating the finite element functions at that point. The subroutine triangulation_search_delaunay is used to find the element triangle that contains the current sample point, j, within the loop. Once the corresponding triangle is found, the finite element basis functions are evaluated at that point within the triangle and stored in t_node. Finally, since the goal of this routine is to return the sample values, we'll need to multiply the vector $\mathbf{b}$ (returned from either basis_mn_t3 or basis_mn_t6) by the finite element values. This is performed in the nested loop multiply with the intrinsic function DOT_PRODUCT.

### SUBROUTINE triangulation_order3_neighbor_triangles

SUBROUTINE triangulation_order3_neighbor_triangles(triangle_num, triangle_node, triangle_neighbor)
IMPLICIT NONE
INTEGER(kind=4) :: triangle_num, col(4, 3 * triangle_num), i, icol, j, k, side1, side2, triangle_neighbor(3, triangle_num), tri, &
tri1, tri2
INTEGER(kind=4), PARAMETER :: triangle_order = 3
INTEGER(kind=4) :: triangle_node(triangle_order, triangle_num)
step_1_order3: DO tri = 1, triangle_num !Step 1
i = triangle_node(1, tri)
j = triangle_node(2, tri)
k = triangle_node(3, tri)
IF (i < j) THEN
col(1:4, 3 * (tri - 1) + 1) = [i, j, 1, tri]
ELSE
col(1:4, 3 * (tri - 1) + 1) = [j, i, 1, tri]
END IF
IF (j < k) THEN
col(1:4, 3 * (tri - 1) + 2) = [j, k, 2, tri]
ELSE
col(1:4, 3 * (tri - 1) + 2) = [k, j, 2, tri]
END IF
IF (k < i) THEN
col(1:4, 3 * (tri - 1) + 3) = [k, i, 3, tri]
ELSE
col(1:4, 3 * (tri - 1) + 3) = [i, k, 3, tri]
END IF
END DO step_1_order3
CALL i4col_sort_a(4, 3 * triangle_num, col) !Step 2
triangle_neighbor(1:3, 1:triangle_num) = -1 !Step 3
icol = 1
finalize_order3: DO
IF (3 * triangle_num ≤ icol) EXIT
IF (col(1, icol) /= col(1, icol + 1) .OR. col(2, icol) /= col(2, icol + 1)) THEN
icol = icol + 1
CYCLE
END IF
side1 = col(3, icol)
tri1 = col(4, icol)
side2 = col(3, icol + 1)
tri2 = col(4, icol + 1)
triangle_neighbor(side1, tri1) = tri2
triangle_neighbor(side2, tri2) = tri1
icol = icol + 2
END DO finalize_order3
END SUBROUTINE triangulation_order3_neighbor_triangles


This subroutine determines triangle neighbors. A triangulation of a set of nodes can be completely described by the coordinates of the nodes, and the list of nodes that make up each triangle. However, in some cases, it is necessary to know triangle adjacency information; that is, which triangle, if any, is adjacent to a given triangle on a particular side. This routine creates a data structure recording this information. The primary amount of work occurs in sorting a list of 3 * triangle_num data items. This routine stores the edge data in the array col.

By way of example, suppose the input information from triangle_node is:

Triangle Nodes -------- ---------------- 1 3 4 1 2 3 1 2 3 3 2 8 4 2 1 5 5 8 2 13 6 8 13 9 7 3 8 9 8 13 2 5 9 9 13 7 10 7 13 5 11 6 7 5 12 9 7 6 13 10 9 6 14 6 5 12 15 11 6 12 16 10 6 11

The output information in triangle_neighbor would be:

Triangle Neighboring Triangles -------- --------------------- 1 -1 -1 2 2 1 4 3 3 2 5 7 4 2 -1 8 5 3 8 6 6 5 9 7 7 3 6 -1 8 5 4 10 9 6 10 12 10 9 8 11 11 12 10 14 12 9 11 13 13 -1 12 16 14 11 -1 15 15 16 14 -1 16 13 15 -1

The 4-byte integers and integer arrays that go into this subroutine are: triangle_num = number of triangles; triangle_node(3, triangle_num) = nodes that make up each triangle. On the output side, the 4-byte integer array triangle_neighbor(3, triangle_num) contains the three triangles that are direct neighbors of a given triangle. triangle_neighbor(1, i) is the index of the triangle which touches side 1, defined by nodes 2 and 3, and so on; triangle_neighbor(1, i) is negative if there is no neighbor on that side. In this case, that side of the triangle lies on the boundary of the triangulation.

So how does this subroutine work?

Step 1. From the list of nodes for triangle t of the form $(i, j, k)$, construct the three neighbor relations

 $\displaystyle{ \begin{matrix} (i, j, 1, t) \text{ or } (j, i, 1, t) \text{;} \\ (j, k, 2, t) \text{ or } (k, j, 2, t) \text{;} \\ (k, i, 3, t) \text{ or } (i, k, 3, t) \text{,} \end{matrix} }$

where we choose $(i, j, 1, t)$ if $i \lt j$, or else $(j, i, 1, t)$.

Step 2. Perform an ascending dictionary sort of the neighbor relations. We only intend to sort rows 1 and 2; the routine we call here sorts rows 1 through 4, but that won't hurt us. What we need is to find cases where two triangles share an edge. Say they share an edge defined by the nodes i and j. Then there are two columns of col that start out as $(i, j, ?, ?)$. By sorting col, we make sure that these two columns occur consecutively. That will make it easy to notice that the triangles are neighbors.

Step 3. Neighboring triangles show up as consecutive columns, with the first 2 entries being identical. Whenever the routine spots this happening, it will make the appropriate entries in triangle_neighbor.

The subroutine triangulation_order6_neighbor_triangles performs a task similar to this one, except that triangle_order = 6 is declared and the quantities passed to i4col_sort_a are reversed.

### SUBROUTINE triangulation_search_delaunay

SUBROUTINE triangulation_search_delaunay(node_num, node_xy, triangle_order, triangle_num, triangle_node, triangle_neighbor, p,     &
triangle_index, edge)
IMPLICIT NONE
INTEGER(kind=4) :: node_num, triangle_num, triangle_order, a, b, c, count, edge, i4_uniform, seed,                               &
triangle_node(triangle_order, triangle_num), triangle_index, triangle_neighbor(3, triangle_num)
INTEGER(kind=4), PARAMETER :: dim_num = 2
REAL(kind=8) :: alpha, beta, det, dxp, dxa, dxb, dyp, dya, dyb, gamma, node_xy(dim_num, node_num), p(dim_num)
count = 0
edge = 0
CALL get_seed(seed)
triangle_index = i4_uniform(1, triangle_num, seed)
delaunay: DO
count = count + 1
IF (triangle_num < count) THEN
WRITE(*, '(a)') " "
WRITE(*, '(a)') "TRIANGULATION_SEARCH_DELAUNAY - Fatal error!"
WRITE(*, '(a)') "The algorithm seems to be cycling."
triangle_index = -1
edge = -1
STOP 1
END IF
a = triangle_node(1, triangle_index) !get the nodes of triangle TRIANGLE_INDEX
b = triangle_node(2, triangle_index)
c = triangle_node(3, triangle_index)
dxa = node_xy(1, a) - node_xy(1, c) !using vertex C as a base, compute the distances to vertices A and B, and the point P
dya = node_xy(2, a) - node_xy(2, c)
dxb = node_xy(1, b) - node_xy(1, c)
dyb = node_xy(2, b) - node_xy(2, c)
dxp = p(1)          - node_xy(1, c)
dyp = p(2)          - node_xy(2, c)
det = dxa * dyb - dya * dxb
alpha = (dxp * dyb - dyp * dxb) / det !compute the barycentric coordinates of the point P with respect to this triangle
beta =  (dxa * dyp - dya * dxp) / det
gamma = 1.0D+00 - alpha - beta
IF (0.0D+00 ≤ alpha .AND. 0.0D+00 ≤ beta .AND. 0.0D+00 ≤ gamma) EXIT !if the barycentric coordinates are all positive, then the point is inside the triangle and we're done
IF (alpha < 0.0D+00 .AND. 0 < triangle_neighbor(2, triangle_index)) THEN !at least one barycentric coordinate is negative
triangle_index = triangle_neighbor(2, triangle_index) !if there is a negative barycentric coordinate for which there exists an opposing triangle neighbor closer to the point, move to that triangle
CYCLE
ELSE IF (beta < 0.0D+00 .AND. 0 < triangle_neighbor(3, triangle_index)) THEN !two coordinates could be negative, in which case we could go for the most negative one, or the most negative one normalized by the actual distance it represents
triangle_index = triangle_neighbor(3, triangle_index)
CYCLE
ELSE IF (gamma < 0.0D+00 .AND. 0 < triangle_neighbor(1, triangle_index)) THEN
triangle_index = triangle_neighbor(1, triangle_index)
CYCLE
END IF
IF (alpha < 0.0D+00) THEN !all negative barycentric coordinates correspond to vertices on opposite sides on the convex hull
edge = -2 !note the edge and exit
EXIT
ELSE IF (beta < 0.0D+00) THEN
edge = -3
EXIT
ELSE IF (gamma < 0.0D+00) THEN
edge = -1
EXIT
END IF
END DO delaunay
END SUBROUTINE triangulation_search_delaunay


This subroutine searches a Delaunay triangulation for a point. The algorithm "walks" from one triangle to its neighboring triangle, and so on, until a triangle is found containing point $\mathcal{P}$, or $\mathcal{P}$ is found to be outside the convex hull.

The algorithm computes the barycentric coordinates of the point with respect to the current triangle. If all three quantities are positive, the point is contained in the triangle. If the $i$th coordinate is negative, then $\mathcal{P}$ lies on the far side of edge 1, which is opposite from vertex 1. This gives a hint as to where to search next.

For a Delaunay triangulation, the search is guaranteed to terminate. For other triangulations, a cycle may occur.

Note the surprising fact that, even for a Delaunay triangulation of a set of nodes, the nearest node to $\mathcal{P}$ need not be one of the vertices of the triangle containing $\mathcal{P}$.

The code can be called for triangulations of any order, but only the first three nodes in each triangle are considered. Thus, if higher order triangles are used, and the extra nodes are intended to give the triangle a polygonal shape, these will have no effect, and the results obtained here might be misleading.

The 4-byte integers and integer arrays that go into this subroutine are: node_num = number of nodes; triangle_order = order of the triangles; triangle_num = number of triangles; triangle_node(triangle_order, triangle_num) = nodes that make up each triangle; triangle_neighbor(3, triangle_num) = triangle neighbor list. Here are the double precision entities going in: node_xy(2, node_num) = coordinates of the nodes; p(2) = coordinates of a point $\mathcal{P}$. On the output side, the 4-byte integer triangle_index is the index of the triangle where the search ended. If a CYCLE occurred, then triangle_index = -1. Another output integer is edge, which indicates the position of the point $\mathcal{P}$ in triangle triangle_index. The values of triangle_index have the following meanings: 0, the interior or boundary of the triangle; -1, outside the convex hull of the triangulation, past edge 1; -2, outside the convex hull of the triangulation, past edge 2; -3, outside the convex hull of the triangulation, past edge 3.

### References

• Burkardt, John. "FEM2D_SAMPLE Evaluate a Finite Element Function of a 2D Argument." Florida State University, Department of Scientific Computing.
• Joe, Barry. "GEOMPACK – a software package for the generation of meshes using geometric algorithms." Advances in Engineering Software 13, 325-331 (1991).