# Finite Element Method - 2D Mesh Generator - Metfem2D

## Delaunay triangulation algorithm

Let's remind briefly the main points of the Delaunay triangulation method together with their numerical implementation using Octave and Matlab software. Let P = {pi, i = 1, 2, … ,N} be a set of points in two – dimensional Euclidean plane ℜ2. They are called forming points of mesh. Let us define the triangle T as a set of three mesh points

 T = {tj ∈ P, j = 1,2,3}.

### Definition of Delaunay zones

Using the Delaunay criterion one can generate triangulation where no four points from the set of forming points P are co – circular:

 ∀ pi ∈ P ∧ pi ∉ T ||pi - u|| > ρ2

where u is the center of the T triangle and ρ is its radius. The proposed algorithm consists of the following steps:

• The triangle's bars are given by the following vectors t12, t13, t23 where tij = tj - ti, ti = [txi, tyi, 0] for each i ≠ j and i, j ∈ {1, 2, 3}.
• The cross product of each triangle bars defines a plane. The pseudovector A together with its projection on the normal to the plane n – direction An are found
 A = t12 - t13 An = n ⋅ A
in order to determine the triangle orientation. If the quantity An > 0 the triangle orientation is clockwise, otherwise is counterclockwise.
• The determinant of the square matrix D(T) is built on the basis of the set of triangle's nodes given by the equation

D(T) = det (  tx1 ty1 (tx1)2 + (ty1)2 tx2 ty2 (tx2)2 + (ty2)2 tx3 ty3 (tx3)2 + (ty3)2
)

Figure shows the main idea of the Delaunay criterion. a) Two triangles (with nodes abc and acd, respectively) are not Delaunay triangles, b) after exchange of the edge ac to the edge bd two new triangles abd and bcd replace the old ones. They are both of the Delaunay type. Circles represent the Delaunay zones.

next the following determinant is calculated in order to find out whether a mesh point pi is outside or inside the Delaunay zone (see Fig. 13)

D(T)i = det (  tx1 ty1 (tx1)2 + (ty1)2 1 tx2 ty2 (tx2)2 + (ty2)2 1 tx3 ty3 (tx3)2 + (ty3)2 1 pxi pyi (pxi)2 + (pyi)2 1
)
for each pi ∈ P ∧ pi ∉ T.
• If for any point pi its D(T)i < 0 the triangle T is not the Delaunay triangle (see Fig. 13a). Then one need to find other triangles in the closest neighbourhood of the triangle T corresponding to the number of pi inside the Delaunay zone and recursively exchange the bars between T and each of them (see Fig. 13b).
• Finally, checking whether the new two triangles are the Delaunay triangles takes its place. If so, new ones are accepted unless the change is canceled.

The algorithm ends up with the new triangular mesh Ωnew.
If you wish you can have an insight into the program delaunayTakeTechEase.m that implements the above – presented algorithm using the Octave and Matlab software (it is a part of a free Octave and Matlab course). One can use also the appropriate Octave and Matlab library function.
Help. In order to find an orientation of a triangle T one can check the sign of An (see equation). If it is greater than 0 the triangle orientation is clockwise unless counterclockwise. In the latter case, to ensure the clockwise orientation one can once flip up and down matrix in equation then the triangle orientation turns into the opposite one. Obviously, this flipping results in the change of the sign of the matrix determinant D(T) → -D(T).

### References

1. ^ B. Delaunay, Sur la sphre vide, Izvestia Akademii Nauk SSSR, Otdelenie Matematicheskikh i Estestvennykh Nauk, 7, pp. 793-800, 1934
2. ^The source code of Octave is freely distributed GNU project, for more info please go to the following web page http://www.gnu.org/software/octave/.