Finite Element Method - 2D Mesh Generator - Metfem2D

Delaunay triangulation algorithm

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

T = {t_j ∈ 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:

∀ p_i ∈ P ∧ p_i ∉ T

||p_i - u|| > ρ²

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 t₁₂, t₁₃, t₂₃ where t_ij = t_j - t_i, t_i = [t^x_i, t^y_i, 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 A_n are found

A = t₁₂ - t₁₃

A_n = n ⋅ A

in order to determine the triangle orientation. If the quantity A_n > 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 (

t^x₁ t^y₁ (t^x₁)² + (t^y₁)²

t^x₂ t^y₂ (t^x₂)² + (t^y₂)²

t^x₃ t^y₃ (t^x₃)² + (t^y₃)²

)

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 p_i is outside or inside the Delaunay zone (see Fig. 13)

D(T)_i = det (

t^x₁ t^y₁ (t^x₁)² + (t^y₁)² 1

t^x₂ t^y₂ (t^x₂)² + (t^y₂)² 1

t^x₃ t^y₃ (t^x₃)² + (t^y₃)² 1

p^x_i p^y_i (p^x_i)² + (p^y_i)² 1

)

for each p_i ∈ P ∧ p_i ∉ T.
If for any point p_i 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 p_i 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 A_n (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).

Finite Element Method - 2D Mesh Generator - Metfem2D

Delaunay triangulation algorithm

Definition of Delaunay zones

References