A simple technique for unstructured mesh generation via adaptive finite elements

This work describes a concise algorithm for the generation of triangular meshes with the help of standard adaptive finite element methods. We demonstrate that a generic adaptive finite element solver can be repurposed into a triangular mesh generator if a robust mesh smoothing algorithm is applied between the mesh refinement steps. We present an implementation of the mesh generator and demonstrate the resulting meshes via several examples.


Introduction
Many numerical methods for partial differential equations (PDE's), such as the finite element method (FEM) and the finite volume method (FVM), are based on splitting the domain of the solution into primitive shapes such as triangles or tetrahedra. The collection of the primitive shapes, i.e. the computational mesh, is used to define the discretisation, e.g., in the FEM, the shape functions are polynomial in each mesh element, and in the FVM, the discrete fluxes are defined over the edges of the cells.
This article describes a simple approach for the triangulation of twodimensional polygonal domains. The process can be summarised as follows: 1. Find a constrained Delaunay triangulation (CDT) of the polygonal domain using the corner points as input vertices and the edges as constraints.
2. Solve the Laplace equation with the given triangulation and the FEM.
3. Split the triangles with the largest error indicator using adaptive mesh refinement techniques. 4. Apply centroidal patch tesselation (CPT) smoothing to the resulting triangulation.
It is noteworthy that the steps 2, 3 and 5 correspond exactly to what is done in any implementation of the standard adaptive FEM; cf. Verfürth [22] who calls it the adaptive process.
The goal of this work is to demonstrate that if the mesh smoothing algorithm of step 4 is chosen properly, the adaptive process tends to produce reasonable meshes even if the initial mesh is of low quality. Thus, we demonstrate that the adaptive process-together with an implementation of the CDT and a robust mesh smoothing algorithm-can act as a simple triangular mesh generator.

Prior work
Two popular techniques for generating unstructured meshes are those based on the advancing front technique [14] or the Delaunay mesh refinement [8,18,20]. In addition, there exists several less known techniques such as quadtree meshing [24], bubble packing [21], and hybrid techniques combining some of the above [15].
Some existing techniques bear similarity to ours. For example, Bossen-Heckbert [4] start with a CDT and improve it by relocating the nodes. Instead of randomly picking nodes for relocation, we use a finite element error indicator that guides the refinement. Instead of doing local modifications, we split simultaneously all triangles that have their error indicators above a predefined threshold.
Persson-Strang [17] describe another technique based on iterative relocation of the nodes. An initial mesh is given by a structured background mesh which is then relaxed by intepreting the edges as a precompressed truss structure. The structure is forced inside a given domain by expressing the boundary using signed distance functions and interpreting the signed distance as an external load acting on the truss. In contrast to the present approach, the geometry description is implicit, i.e. the boundary is defined as the zero set of a user-given distance function.
We do not expect our technique to surpass the existing techniques in the quality of the resulting meshes or in the computational efficiency. However, the algorithm can be easier to understand for those with a background in the finite element method and, hence, it may be a viable candidate for supplementing adaptive finite element solvers with basic mesh generation capabilities.

Components of the mesh generator
The input to our mesh generator is a sequence of N corner points that form a polygon when connected by the edges We do not allow self-intersecting polygons although the algorithm can be generalised to polygons with holes. The corresponding polygonal domain is denoted by Ω C ⊂ R 2 .

Constrained Delaunay triangulation
A triangulation of Ω C is a collection of nonoverlapping nondegenerate triangles whose union is exactly Ω C . Our initial triangulation T 0 is a constrained Delaunay triangulation (CDT) of the input vertices C with the edges (C 1 , C 2 ), (C 2 , C 3 ), . . ., (C N −1 , C N ), (C N , C 1 ) constrained to be a part of the resulting triangulation and the triangles outside the polygon removed; cf. Chew [7] for the exact definition of a CDT and an algorithm for its construction.
An example initial triangulation of a polygon with a spiral-shaped boundary is given in Figure 1. It is obvious that the CDT is not always a high quality computational mesh due to the presence of arbitrarily small angles. Thus, we seek to improve the initial triangulation by iteratively adding new triangles, and smoothing the mesh. Note that the remaining steps do not assume the use of CDT as an initial triangulation-any triangulation with the prescribed edges will suffice.

Solving the Laplace equation
In order to decide on the placement of the new vertices and triangles, we solve the Laplace equation 1 using the FEM and evalute the corresponding a posteriori error estimator. The triangles that have the highest values of the error estimator are refined, i.e. split into smaller triangles.
The Laplace equation reads: find u : The finite element method is used to numerically solve the weak formulation: find u ∈ V such that where w ∈ V if w| ∂Ω C = 0 and Ω C (∇w) 2 dx < ∞. We denote the kth triangulation of the domain Ω C by T k , k = 0, 1, . . ., and use the piecewise linear polynomial space where P 1 (T ) denotes the set of linear polynomials over T . The finite element method corresponding to the kth iteration reads: find u k h ∈ V k h such that The local a posteriori error estimator reads where A T is the area of the triangle T and h T is the length of its longest edge, w | ∂T \∂Ω C denotes the jump in the values of w over ∂T \ ∂Ω C , and n is a unit normal vector to ∂T . The error estimator η T is evaluated for each triangle after solving (4). Finally, a triangle T ∈ T k is marked for refinement if where 0 < θ < 1 is a parameter controlling the amount of elements to split during each iteration. [22] 3

.3 Red-green-blue refinement
The triangles marked for refinement by the rule (6) are split into four. In order to keep the rest of the triangulation conformal, i.e. to not have nodes in the middle of an edge, the neighboring triangles are split into two or three by the so-called red-green-blue (RGB) refinement; cf. Carstensen [5]. Using RGB refinement to the example of Figure 1 is depicted in Figure 2.

Centroidal patch triangulation smoothing
We use a mesh smoothing approach introduced by Chen-Holst [6] who refer to the algorithm as centroidal patch triangulation (CPT) smoothing. The idea is to repeatedly move the interior vertices to the area-weighted averages of the barycentres of the surrounding triangles. The CPT smoothing is combined with an edge flipping algorithm, also described in Chen-Holst [6], to improve the quality of the resulting triangulation. The mesh smoother is applied to the spiral-shaped domain example in Figure 3.

The mesh generation algorithm
In previous sections we presented an overview of all the components of the mesh generation algorithm. The resulting mesh generator is now summarised in Algorithm 1. The total number of refinements M is a constant to guarantee the termination of the algorithm. Nevertheless, in practice and in our implementation the refinement loop is terminated when a quality criterion is satisfied, e.g., when the average minimum angle of the triangles is above a predefined threshold. The entire mesh generation process for the spiralshaped domain example is given in Figure 4. for k ← 1 to M do 5: T k ← CPT(T k ) 7: return T N Figure 4: The entire mesh generation process for the spiral-shaped domain example from left-to-right, top-to-bottom.

Implementation and examples
We created a prototype of the mesh generator in Python for computational experiments [9]. The implementation relies heavily on the scientific Python ecosystem [23]. It includes source code from pre-existing Python packages tri [3] (CDT implementation, ported from Python 2) and the older MITlicensed versions of optimesh [2] (CPT smoothing) and meshplex [1] (edge flipping). Moreover, it dynamically imports modules from scikit-fem [11] (RGB refinement) and matplotlib [13] (visualization). Some example meshes are given in Figure 5. By default our implementation uses the average triangle quality 2 above 0.9 as a stopping criterion which can lead to individual slit triangles. This is visible especially in the last two examples that have small interior angles on the boundary.

An example application
Mayzel-Steinberg-Varshney [16] report on measurements using a microfluidic device which is analogous to electron transport in two-dimensional conductive materials such as graphene. The governing equation is the Stokes system −∇ · (2µ ε(∇u)) + ∇p = 0, where the unknowns are the velocity u : Ω → R 2 and the pressure p : Ω → R, the parameter µ > 0 is given, and ε(∇u) = 1 2 (∇u+∇u T ) is the rate-of-strain tensor.
We extract the boundary ∂Ω of the computational domain Ω ⊂ R 2 from a photograph of the experiment and input it to the mesh generator; the photograph of the experiment is given in Figure 6a and the resulting computational mesh in Figure 6b. In order to simplify the problem setup we choose the computational domain so that the left boundary corresponds to the vertical streamlines in the micro-particle image velocimetry; see Figure 6c. This allows using the Dirichlet boundary condition u x = 0, u y = −U , U > 0, on the left boundary.
We aim to reproduce the location of the swirl and, therefore, the exact parameter values U = µ = 1 are irrelevant as they scale only the maximum values of the velocity and the pressure fields, i.e. the shape of the solution will not change. The numerical result, calculated and visualised with the help of the finite volume package FiPy [12] and the plotting library matplotlib [13], is given in Figure 6d.

Conclusions
We introduced an algorithm for the generation of triangular meshes for explicit polygonal domains based on the standard adaptive finite element method and centroidal patch triangulation smoothing. We presented a prototype implementation which demonstrates that many of the resulting meshes are reasonable and have an average triangle quality equal to or above 0.9. The majority of the required components are likely to exist in an implementation of the adaptive finite element method. Therefore, the algorithm can be a compelling candidate when an existing adaptive finite element solver is supplemented with basic mesh generation capabilities. Technically the approach extends to three dimensions although in practice the increase in (a) A photograph of the experiment. [16] (b) The corresponding mesh from adaptmesh.

(c)
Experimental streamlines from micro-particle image velocimetry. [16] (d) Computational streamlines from FiPy [12] and matplotlib [13]. computational effort can be significant and the quality of the resulting tetrahedralisations has not been investigated.