Over the past, a plenty of approaches for the generation of triangular meshes have been developed from which three basic strategies - the quadtree based approach, the advancing front technique, and the Delaunay triangulation - have proved particularly vital in terms of the efficiency, robustness, and reliability. The basic concepts of these techniques are presented in this section.
The quadtree based approach [5,9] uses a spatial decomposition procedure representing the domain as a set of disjoint cells, called quadrants, stored in a hierarchic tree. The domain to be meshed is placed in the root quadrant, usually a square, that entirely encloses it. The structure of the quadtree is then defined in an recursive manner by subdividing the individual quadrants (starting with the root) into four child quadrants until all quadrants are at the desired level, dictated by the geometry of the domain and element size specification. Since the quadtree structure may be highly irregular, a one-level difference, ensuring that the quadrants sharing an edge do not differ in the size by a factor greater than two, is enforced. This guarantees a gradual variation of the final element density. After that, each quadrant is classified to be either inside, or outside, or boundary with respect to the domain being discretized. No further attention is paid to the outside quadrants. The inside quadrants are meshed using templates, predefined sets of elements filling entirely the quadrant. Since the one-level difference rule essentially limits the number of midside nodes to one per a quadrant edge, the overall number of templates is reduced to 16 (if rotations are taken into account). The discretization of boundary quadrants is usually more difficult because of requirements on a consistent and unambiguous way of treating the domain boundary. In the early version , only a limited number of specific locations in the quadrant was reserved for the intersection of the domain boundary with the quadrant. This allowed to use templates also for the meshing of boundary quadrants but a limit on geometrical complexity of the domain being discretized was imposed. In the more recent implementations , where the use of real intersections virtually eliminates the possibility of the application of templates to the boundary quadrants, the triangle removal concept is used to mesh the boundary quadrants. Simultaneously, a small feature elimination process (on the mesh level) is introduced to remove small elements arising by cutting off small portions of the boundary quadrants by the domain boundary. Although the quadtree based approach proved to be very efficient and robust, it suffers from the limited flexibility of element density control due to the virtual predetermination of the element spacing by the size of the root quadrant and its children.
A very attractive method for the generation of triangular meshes is the advancing front technique [11,45,47,57,60,67,75] which is filling the empty (not yet meshed) space by introducing one new element at a time. The key data structure of this approach is the so called front that is defined as an interface between the meshed and unmeshed parts of the domain. Initially, the front consists of segments forming the boundary (outer as well as inner) of the domain. The generation then continues on the basis of a triangle removal algorithm. It consists in the selection of an appropriate segment from the front and in an attempt to complete it to a proper triangular element either by picking a suitable existing node in the local neighbourhood of the segment or by creating a new node. The front is then updated to account for the newly formed element. The original segment together with already existing segments used to form the new element are removed from the front and the newly created segments are added to the front. This process is repeated until there are no remaining segments in the front. The algorithm respects very well the desired element size because of its high quality runtime point placement strategy, which is in contrast with its early version  using nodes distributed in a global process prior to the discretization. Since this method is based on a hierarchical strategy (the new elements are created from mesh entities of the next lower dimension), the integrity of the boundary is guaranteed. The main difficulty in this approach appears to be the need to track the front as it develops and folds over itself, which requires an intensive computational effort to avoid undesirable intersections. Therefore the spatial localization, intersection checking, and front management are considered to be general bottlenecks of the AFT, increasing its computational complexity. Various data structures and strategies [14,26,66,60] have been proposed to alleviate this negative aspect. It should be mentioned that the AFT can be successfully employed also for the generation of anisotropic meshes (meshes with prestretched and preoriented elements) by simply using the metric given by the transformation between the isotropic and anisotropic spaces. However, extremely high stretches of order over thousand, typically used in the computational fluid dynamics (CFD), are causing numerical difficulties due to truncation and rounding errors. Note that the AFT is also suitable for the discretization of 3D surfaces, as will be discussed in next sections.
The most widely used approach for the generation of triangular meshes is based on the Delaunay triangulation [1,18,21,35,44,52,58,62,64,65] which is a method modifying an existing grid by introducing new nodes. In this technique, the discretization of the domain boundary and an initial mesh (possibly consisting of just a single element) satisfying the Delaunay property and surrounding the whole domain have to be available. Each node of the boundary discretization is then inserted into the mesh. Since the node insertion typically violates the Delaunay property, a local mesh transformation is carried out after each node insertion to recover it. This transformation is usually performed on the basis of the Bowyer's  or Watson's  algorithm using the ideas of the Voronoi diagram. All elements, the circumscribed circle of which contains the node being inserted, are removed from the discretization and the created star-shaped cavity is remeshed taking into account the inserted node. Further nodes are then inserted into the interior of the domain. These nodes are created during the runtime [44,65] or are gathered from a point distribution generated before the actual meshing [1,31]. The former approach provides usually much better node distribution with respect to the desired element size variation. Again, the local transformation process must be accomplished after each internal node is inserted. Note that the insertion of internal nodes can be also governed by alternative (usually combined) approaches not necessarily yielding the DT [55,62]. Finally, the original boundary is recovered (by overriding the Delaunay property) to ensure the compatibility between the model and the mesh. This computationally intensive phase may be skipped if a constrained version of the DT is adopted in which the connections near the boundary are restricted to ensure that the boundary discretization remains intact all the time . The popularity of the DT lies especially in a solid mathematical background and also in the fact that a valid mesh is maintained throughout the grid generation process. Thus the grid represented by a simple data structure serves as an efficient tool for the searching and spatial localization. Similarly as the AFT, the DT can be also employed for the generation of anisotropic meshes using the appropriate metric reflecting the desired stretching and orientation of elements . On the other hand, the use of the Delaunay technique for the triangulation of 3D surfaces is not as straightforward, which is given by the fact that the formulation of the Delaunay property in 2.5D (2D elements in 3D space) is not available. This is usually resolved by the application of the Delaunay property in a plane (either in the tangent plane constructed at inserted node or more commonly in the plane of the element in which the inserted node falls ) onto which elements from the local neighbourhood are projected or using a quasi-Delaunay property based on an empty circumsphere criterion in the 3D space .