Next: Smoothing Up: Real Space Meshing Previous: Boundary Curves Discretization


Surface Triangulation

Since the mesh generation is constrained directly to the surface, the standard AFT is not applicable without some modifications allowing to reflect the surface curvature. These amendments, concerning especially the ``ideal'' point placement and the intersection check, become distinct from the detailed description of the individual steps of the implemented AFT presented in this section. The other modifications related to the computational complexity of the algorithm are discussed in Subsection Computational Complexity.

The modified AFT is accomplished according to the following scheme. Initially, the front consisting of edges constituting the boundary of the surface (including inner loops, if any) is established. Once the initial front has been set up, the mesh generation continues on the basis of a triangle removal algorithm, according to the following 11 steps, until the front becomes empty:

  1. The first available edge in the front is made active. The edge is oriented from to and the new triangle is assumed to be generated on its left side when viewing against the outer normal of the surface.

  2. The position of the ``ideal'' point forming a tentative triangle (Fig. 26) is computed as

    (67)


    where

    (68)


    (69)


    The stands for the positional vector, represents the surface outer unit normal, and denotes the element size (extracted from the octree). The subscript refers to a particular point in the 3D space. The coefficient corresponds to the height of the unit equilateral triangle. It should be noted that the weighting of normals in Eq. (69) has only a negligible effect on and therefore the following simplification (depicted in Fig. 26) may be employed

    (70)


  3. Since point is generally out of the surface (Fig. 26), its projection to the surface is calculated using the iterative approach described in Subsection Surface Projection.

  4. The position of the ``ideal'' point is corrected taking into account the local surface curvature (Fig. 27) and the local element size variation (Fig. 28). This corrected ``ideal'' point is referred as and its positional vector is given by

    (71)


    where

    (72)


    and for the modified size holds

    (73)


    Note that correction of the ideal point must be applied repeatedly when a relatively coarse mesh is required on a curved surface.
    The correction due to the local element size variation is demonstrated for a planar case and assumed linear element size variation in Figure 28. The linear profile of the element spacing along the direction given by vector is defined by the size at the base of the tentative triangle and by the element size at the ``ideal'' point which distance from edge equals to . The corrected size then corresponds to such a point which distance from edge is equal to the element size defined by the element spacing profile at that distance. This point can be obtained by the repeated placement of the ``ideal'' point each time using the element size given by the element spacing at the last ``ideal'' point (see points in Fig. 28) which converges to the point corresponding to the intersection of the element spacing profile with the line inclined by passing through the origin.

  5. Point is projected to the surface and the element size corresponding to the projected point is extracted from the octree.

  6. The local neighbourhood of point in terms of a set of terminal octants is established (Fig. 29). It is desirable to make the neighbourhood as small as possible in order to reduce the number of mesh entities that have to be investigated and simultaneously to make it large enough to include all mesh entities that should affect the element creation. Thus the neighbourhood construction is always a compromise between certainty and efficiency. In the presented approach, the following terminal octants are included in the neighbourhood:

    1. octants containing points and ,
    2. octant containing point ,
    3. octants adjacent to the terminal octants selected in case (a) across those of their faces, edges, and corners which are closer to point , where is the relevant from or , than , , and , respectively,
    4. octants adjacent to the terminal octant containing point across those of its faces, edges, and corners which are closer to point than , , and , respectively.
    The denotes the size of the terminal octant comprising point . Note that each octant is considered only once, thus no octant duplicity occurs in the neighbourhood. The advantage of the proposed approach lies in the fact that the neighbourhood is kept relatively small and that the octants specified in cases (c) and (d) can be very efficiently determined by the horizontal octree traversal. In order to decrease the probability that an important mesh entity is not included in the neighbourhood, the neighbourhood is artificially enlarged by introduction of a ratio between the element size and the octant size. The upper bound of this ratio has been determined heuristically and it is currently set to (see Eq. (51)). The lower bound is then set to halve the value of the upper bound.

  7. The neighbourhood is searched for the most suitable candidate to form a new triangle . Two auxiliary geometrical loci are constructed around point to resolve acceptance of the tentative node . Firstly, a sphere with the centre in and radius of is constructed. And secondly, an ellipsoid is derived from the sphere by extension of its half-axis parallel to the active edge by a half (Fig. 30). The selection of point is driven by the following rules. Point  is

    1. identical with point if there is no point inside the sphere , no point inside the ellipsoid directly connected with the active edge , and no point inside the sphere circumscribed to the triangle .
    2. an already existing node inside the sphere circumscribed to the triangle forming the largest angle if there is no point inside the sphere and no point inside the ellipsoid directly connected with the active edge .
    3. taken as node if it falls inside the sphere and there is no other node inside the sphere circumscribed to the triangle .
    4. an already existing node inside the sphere circumscribed to the triangle forming the largest angle . Point is an already existing node inside the sphere (the case depicted in Fig. 30).
    5. taken as node if the already existing edge or forms a sufficiently small angle (currently not larger than ) with the active edge and no point is inside the sphere circumscribed to the triangle .
    6. an already existing node inside the sphere circumscribed to the triangle forming the largest angle , where point is the end node of the already existing edge or being at a sufficiently small angle to the active edge .
    7. taken as node if it falls inside the ellipsoid and there is no other node inside the sphere circumscribed to the triangle . Point is the end node of the already existing edge or .

    8. an already existing node inside the sphere circumscribed to the triangle forming the largest angle , where node falls inside the ellipsoid and is directly connected with the active edge .

    In cases (b), (d), (f), and (h), the validity of point must be always verified. Similarly, point has to be verified in cases (c) and (d). It may happen that although a point is satisfying all the criteria in terms of the metric of the real space, it is still very far from point in the parametric space. In other words, the distance between points and measured on the surface may be much larger than the distance in . However, having access to the parametric space of the surface, the verification can be easily performed. Acceptance of node in cases (g) and (h) improves the mesh quality in the element size transition zones and eliminates creation of chunks of edges when the front closes up. This can be demonstrated on a uniform mesh of a circle with and without the application of the ellipsoidal neighbourhood (Fig. 31). To speed up the neighbourhood search only the active half-space given by point and vector is considered and only nodes on the front (including isolated nodes created by discretization of vertices fixed to the surface) are considered. Furthermore, nodes which are too far from point are ignored. Maximization of angle (for any tentative point ) is used as the primary criterion for the selection of point . This approach exhibits two important advantages. Firstly, maximization of angle is equivalent to the local Delaunay property when discretizing a planar surface. Thus when point is selected in such a way that it maximizes the angle then it is guaranteed that there is no other point in the active half-space inside the circle circumscribed to the tentative triangle . The local character of the Delaunay property consists in the fact that even in such a case, point can be inside the circle circumscribed to any other triangle in the vicinity thus violating the Delaunay property in the global sense. And secondly, the evaluation of angle (more precisely of its cosine) is very cheap. Note that angle is subjected to the correction due to the curvature in order to suppress undesirable element creation in the direction of the higher curvature. The correction is given by the following formulas

    (74)


    (75)


    where denotes the corrected angle , is the projection of point to the surface, and stands for the maximum allowable angle formed by the normal vectors at the ends of a tentative edge (currently set to ). The exact evaluation of , which requires the projection of point to the surface, is substituted by the weighted average of normals and 

    (76)


  8. In order to avoid an overlapping of the tentative triangle with existing triangles and isolated edges in the neighbourhood, an intersection check is carried out. This check can be quite efficiently performed in the parametric space of the surface by testing the intersection of individual edges (represented by straight segments). In this case, the mapping of the whole edge, not only of its nodes, should be taken into account to avoid check failure when relatively long edges on a poorly parameterized surface are considered (Fig. 32). Alternatively, the intersection check can be done in the real space. However, a particular intersection of an edge with a triangle occurring in the parametric space may remain undetected in the real space. Therefore the check, reduced to the test for the intersections of generally non-coplanar edges, has to be carefully implemented. The intersection, specified in terms of the shortest segment connecting the two lines given by the edges, must be verified - the end points of the ``intersection'' segment must lie on the edges, their projections to the surface must be sufficiently close to each other (measured on the surface) and the segment must be approximately parallel with the direction of the projection (perpendicular to the surface). Since it is difficult to quantify what is a relatively long edge or a poorly parameterized surface, the intersection check in the real space is considered to be more robust. If edge or is the already existing edge, the check can be entirely omitted. This reduces the number of intersection investigations to about a half. The intersection check is further enhanced by considering only edges on the front. In the case that an intersection is detected, the selected point is rejected and the search of the neighbourhood is repeated without taking this point into account.

  9. If point has been chosen according to the rule (a), a new node is created and registered into the octree. If any of edges and does not exist, a new edge is formed and connected to its end nodes. Then a new triangle is created and connected to its bounding edges , , and . All newly created mesh entities are classified to the surface being discretized.

  10. The front is updated to account for the newly formed triangle . This means that the already existing triangle edges are removed from the front while the newly created ones are appended to the front. To avoid searching in the front for the edges to be removed from it, these edges are marked as not available, which prevents them to be made active (in the 1st step).

  11. If appropriate, an edge swapping is performed in order to improve the mesh quality. The edge swapping is a topological transformation which changes connectivity properties of the mesh by swapping the diagonal in a quadrilateral (generally non-planar) formed by two adjacent triangles. When discretizing a planar surface the edge swapping is accomplished if edge , classified to the surface, is shared by two triangles and , classified to the same surface, and if node falls inside the circle circumscribed to the triangle (Fig. 33). This criterion can be expressed in terms of angles and using the ratio

    (77)


    The edge swapping is then performed if . This approach can be also generalized to curved surfaces. However, if the surface curvature should be taken into account the ratio is modified, using angles and , to

    (78)


    where

    (79)


    (80)


    The curvature correction coefficients and are given by Eq. (75). It should be mentioned, however, that application of the correction due to the surface curvature has only a negligible effect, unless very coarse meshes are considered.



Next: Smoothing Up: Real Space Meshing Previous: Boundary Curves Discretization

Daniel Rypl
2005-12-07