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:

- 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).
- The position of the ``ideal'' point forming a tentative triangle (Fig. 2.14)
is computed as
(2.64)

where(2.65)

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)

- 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

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. - 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. - 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.
- 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:
- octants containing points and ,
- octant containing point ,
- 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,
- 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 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
- 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 .
- 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 .
- taken as node if it falls inside the sphere and there is no other node inside the sphere circumscribed to the triangle .
- 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).
- 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 .
- 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 .
- 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 .
- 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 .

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

- 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.
- 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.
- 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).
- 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)- the accepted node falls into the circle circumscribed to a triangle on the opposite side of the active edge,
- 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

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

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

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.

*Daniel Rypl
2005-12-07*