Next: Patch Discretization Up: Mesh Generation Previous: Curve Discretization


Surface Discretization

In the present implementation of surface discretization, the mesh generation is constrained directly to the surface. The parametric space of the surface is not explicitly used for the mesh generation but it is used to perform some local operations. A widely known advancing front algorithm with some modifications (respecting the surface curvature) is employed. Firstly, an initial front consisting of edges constituting the boundary of the surface (including inner loops) is established. Once the initial front has been set up the mesh generation continues on the basis of an edge 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 a new triangle is assumed to be generated on the left side of it when viewing against the outer normal of the surface (Fig. 2.14).

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

    (2.64)


    where

    (2.65)


    (2.66)


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

    (2.67)


  3. Since point is generally not constrained to the surface (Fig. 2.14), its projection to the surface is calculated. The projection itself (Fig. 2.15) is accomplished in an iterative manner. Any nearby point (typically an already existing node) on the surface is taken as the first approximation of point . The point is firstly projected to the tangential plane given by point and surface outer unit normal expressed as

    (2.68)


    where gradients and are calculated according to

    (2.69)


    The position of point projected to the tangential plane may be written as

    (2.70)


    The increments of curvilinear coordinate and used to get better approximation of point are given by the relation

    (2.71)


    Note that only two scalar equations from the above vectorial equation can be used to work out and . The selection is done according to the normal vector in such a way that the equation which corresponds to the absolutely maximal component of is not considered. In other words, Eq. (2.71) is solved in that Cartesian plane which is most closely approaching the tangential plane . The increments and are subjected to a normalization the aim of which is to avoid running and out of the range of the parametric space and to prevent a diverging oscillation of the position of point on poorly parameterized surfaces. The approximation of point is then improved by incrementing its curvilinear coordinates by normalized and . The whole process is repeated until the desired accuracy in the approximation of point is achieved. Note that only two iterations are typically used for the projection on a planar well parameterized surface. The complete algorithm of the point projection is summarized in Table 2.3.

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

    (2.72)


    where

    (2.73)


    The modified mesh size , reflecting the local mesh size variation (see Fig. 2.17), may be written as

    (2.74)


    Note that correction must be applied repeatedly when a relatively coarse mesh is required on a curved surface.

  5. Point is projected to the surface using the same technique as in the 3rd step. The mesh 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. 2.18). 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 present 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 octant size. The upper bound of this ratio has been determined heuristically and it is currently set to . The lower bound is then set to half 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 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. 2.19). 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. 2.19).
    5. taken as node if the already existing edge or forms a sufficiently small angle (currently set to ) 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 (geodesic distance) 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 mesh 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. 2.20). 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

    (2.75)


    (2.76)


    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 a computationally expensive projection of point to the surface, is substituted by the weighted average of normals and 

    (2.77)


  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 the edges (represented by straight segments) of the triangles. Caution should be taken when relatively long edges on a poorly parameterized surface are considered (Fig. 2.21). In this case, the mapping of the whole edge, not only of its nodes, should be taken into account. 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. When a planar surface is being discretized the intersection check is accomplished in the surface plane. If the edge or is the 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 generated 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). When discretizing a planar surface the edge swapping is accomplished if edge , classified to the surface, is shared by two triangles and , also classified to the surface, and if node falls inside the circle circumscribed to the triangle (Fig. 2.22). This criterion can be expressed in terms of angles and using the ratio

    (2.78)


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

    (2.79)


    where

    (2.80)


    (2.81)


    The curvature correction coefficients and are given by Eq. (2.76). 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. In planar triangulations, the diagonal edge swapping contributes to the recovery of the Delaunay property which is generally violated in two cases (in the 7th step)

    1. the accepted node falls into the circle circumscribed to a triangle on the opposite side of the active edge,
    2. the accepted node falls into the circle circumscribed to a triangle in the neighbourhood except the triangle on the opposite side of the active edge.

    The first case is resolved by the diagonal edge swapping. The second case is reduced to the first one because the front is implemented as a non-priority queue (first-in - first-out). This is schematically demonstrated in Figure 2.23 where the edge numbers denote edge ordering in the front and the circles correspond to the triangles with the violated Delaunay property. In this way, a 2D boundary constrained Delaunay triangulation can be constructed.

Although a large effort is devoted to the creation of well-shaped elements already during the advancing front procedure, some badly-shaped triangles may be comprised in the grid. Therefore, a Laplacian smoothing technique is used to optimize the shape of the elements. Since the mapping between the parametric and real spaces of the surface is generally not affine, the smoothing is accomplished in the real space using the following iterative relation

(2.82)


where is the node being smoothed, are nodes connected to node by an edge, and is a relaxation parameter given by

(2.83)


The relaxation parameter can be either kept constant, by appropriately modifying weight for each node with respect to the particular number of surrounding nodes, or varying with , while keeping the constant. In the current implementation, the first alternative is used for the smoothing of planar surfaces by setting to and the second alternative is applied for the smoothing of non-planar surfaces by setting to . The optimized grid is usually obtained after only a few cycles of the smoothing (typically up to five). Note that only nodes which are classified to a given surface are subjected to the smoothing. Contrary to the smoothing carried out in 2D, the repositioning of a node is likely to shift the node out of the surface. Therefore the projection is employed to satisfy the surface constraint. The same projection technique as the one used in the advancing front procedure (the 3rd step) is used. Since the original position of the node being smoothed may be used as the initial approximation, better convergence rates are achieved. The stability of the projection can be further improved (especially when very coarse meshes are required) by increasing the weight resulting in a reduction of the relaxation parameter . Note that the smoothing of a planar triangulation can potentionally spoil the Delaunay property. Its recovery may be accomplished by a combination of the smoothing with the diagonal edge swapping. The Laplacian smoothing can be optionally extended by weighting based on the required element size, nodal connectivity, or both. In this case, the first part of Eq. (2.82) can be rewritten as

(2.84)


with equal to for the element size weighting, or to for the connectivity weighting, or to for the combined weighting. The mesh size based weight and connectivity based weight are given by

(2.85)


(2.86)


where and denote the actual and ``optimal'' valency (number of connected nodes on the surface) of node and and are conveniently selected coefficients currently set to and . The optimal valency can be expressed as

(2.87)


where represents the angle (in degrees) in the tangential plane at node which is filled in by the surface and is the number of gaps in that filling. Thus is equal to for surface internal nodes and to for surface boundary nodes on a  continuous boundary curve. In the presented work, a simplified approach has been adopted, in which is equal to for nodes classified to the surface or to otherwise. It should be stated that the computational cost of the otherwise very cheap Laplacian smoothing is considerably increased when the mesh size based weighting is performed and that the application of the weighting does not have a major impact on the mesh quality. Also note that in order to save memory, by storing the nodal position only once, the values of and in Eqs (2.82) and (2.84) are replaced by and , respectively if the -th iteration of the position of node is already available. The deficiency of this replacement, residing in the fact that the potential symmetry of a mesh is disturbed, can be ignored because the advancing front technique typically does not tend to produce symmetrical meshes.



Next: Patch Discretization Up: Mesh Generation Previous: Curve Discretization

Daniel Rypl
2005-12-07