Region Discretization

The region discretization is also based on the advancing front technique. An initial front is formed by triangles which are bounding the region under consideration and which have been created during the surface discretization. The mesh generation then continues on the basis of a triangle removal algorithm, using similar steps as in the case of the surface discretization, until there are no remaining triangles in the front:

- The first available triangle in the front is made active. The
nodes , , and are oriented in the anticlockwise sense with
respect to normal
pointing inside (the not yet discretized part of)
the region (Fig. 2.24).
- The location of the ``ideal'' point forming a tentative
tetrahedron (Fig. 2.24) is calculated as
(2.88)

where(2.89)

The coefficient corresponds to the height of a unit equilateral tetrahedron. The positions of points , , and on the sides of the active triangle and associated mesh sizes are calculated according to the following formulas (analogous to those used in Section Surface Discretization (Sequential Mesh Generation))(2.90)

(2.91)

(2.92)

- In this step, the construction of the octree covering a part of the
region is continued (see
Section Octree Data Structure
(Sequential Mesh Generation)).
The terminal octant comprising point is localized. If this octant
is at an octree level lower than is the one corresponding to the maximum
element size specification, it is refined to that level. In this case,
all nodes stored in the original terminal octant are once more
registered into the octree. Note that the octant is typically refined
just by one level. This is caused by the character of the
front propagation inside the region and by the application of the one-level
difference rule during the octree building.
- The corrected position of the ``ideal'' point is computed (this point is now
referred as ) taking into account the local mesh size gradation in
the same way as during the surface discretization
(2.93)

where(2.94)

- The local neighbourhood of point in terms of a set of
terminal octants is established. A similar approach, as described in
Section Surface Discretization
(Sequential Mesh Generation), is adopted. This time, the
following terminal octants are included into 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 tetrahedron . Again, two auxiliary geometrical
loci are constructed around point to resolve the acceptance
of tentative node - a sphere with the centre in
and radius of and an ellipsoid
derived from the sphere by extension of its half-axis parallel to
the active triangle by a half. A similar set of rules as in the
case of the surface discretization (Fig. 2.19) is used to select point .
Point is
- identical with point if there is no point inside the sphere , no point inside the ellipsoid forming a triangle sharing an edge with the active triangle , and no point inside the sphere circumscribed to the tetrahedron .
- an already existing node inside the sphere circumscribed to the tetrahedron if there is no point inside the sphere circumscribed to the tetrahedron , no point inside the sphere , and no point inside the ellipsoid forming a triangle sharing an edge with the active triangle .
- taken as node if it falls inside the sphere and there is no other node inside the sphere circumscribed to the tetrahedron .
- an already existing node inside the sphere circumscribed to the tetrahedron if there is no point inside the sphere circumscribed to the tetrahedron . Point is an already existing node inside the sphere .
- taken as node if the already existing triangle , , or forms a sufficiently small angle with the active triangle and no point is inside the sphere circumscribed to the tetrahedron .
- an already existing node inside the sphere circumscribed to the tetrahedron if there is no point inside the sphere circumscribed to the tetrahedron . Point is the node of the already existing triangle , , or being at a sufficiently small angle to the active triangle .
- taken as node if it falls inside the ellipsoid and there is no other node inside the sphere circumscribed to the tetrahedron . Point is the node of the already existing triangle , , or .
- an already existing node inside the sphere circumscribed to the tetrahedron if there is no point inside the sphere circumscribed to the tetrahedron . Node falls inside the ellipsoid and forms a triangle sharing an edge with the active triangle .

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. When searching in the neighbourhood for the most suitable candidate, the local Delaunay property is actually used as the selection criterion. This reduces the probability of a consequent occurrence of an intersection. The application of a criterion based on the maximal solid angle at node does not bring in any advantage and it is computationally significantly more demanding. From the implementation point of view, the Delaunay circumsphere property is evaluated in terms of the distance of the centre of the sphere circumscribed to the tentative tetrahedron from the active triangle. The position of the centre of the sphere circumscribed to the tetrahedron can be expressed as (Fig. 2.25)(2.95)

where the signed distance (positive in direction of ) is given by(2.96)

The position of the centre and radius of the circle circumscribed to the active triangle are (for completeness) given by(2.97)

with(2.98)

where stands for the centre of edge , is the length of edge , and denotes the angle at node . Considering only the nodes in the active half-space given by point and vector , then the tetrahedron, formed by the active triangle and a point for which takes the minimum value, satisfies the empty circumsphere property. To speed up the neighbourhood search only the active half-space and only the nodes on the front (including the isolated nodes created by the discretization of vertices or curves inside the region) are considered. The nodes which are too far from point are ignored. If no point has been selected the following actions are taken to resolve the problem. Firstly, a new ``ideal'' point is established, using a reduced element size , and the search of the neighbourhood is repeated. And secondly, if this ``ideal'' point repositioning repeatedly fails, existing tetrahedrons in the local neighbourhood are deleted, updating simultaneously the front, and the created cavity is remeshed. - To avoid an overlapping of the tentative tetrahedron with existing
tetrahedrons and isolated triangles and edges in the
neighbourhood, an intersection check is performed. The intersection
check is computationally much more demanding, with respect to the one
performed during the surface discretization,
for the following reasons. Firstly, it is done fully in 3D and
secondly, a significantly larger number of mesh entities must be
considered as the neighbourhood becomes truly three-dimensional.
Also the number of cases in which the intersection check may be
entirely omitted is considerably reduced. In the present
implementation, the non-existing boundary edges and triangles bounding
the tentative tetrahedron are tested against the intersection with edges
and triangles on the front in the neighbourhood. The intersection
check is further limited only to those edges and triangles from the
neighbourhood which have at least one point in common with the (solid) sphere
circumscribed to the tentative tetrahedron. This criterion effectively
replaces the reduction of the number of mesh entities involved in the
intersection check based on the investigation of their bounding box.
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. As a consequence, the local Delaunay
property is violated resulting in an increased probability of
occurrence of further intersections and in a necessity to repeat the neighbourhood
search again, leading, in the most pathological case, to the repositioning
of the ``ideal'' point or local remeshing.
- 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. Also, if necessary, new triangles , , and are
created and connected to their bounding edges.
Then a new tetrahedron is generated and stored at appropriate
sides of its bounding triangles , , , and . All
newly created mesh entities are classified to the region being discretized.
- The front is updated to account for the newly formed tetrahedron . The already existing triangles bounding the tetrahedron are removed from the front and the newly created ones are appended to the front. Again, to avoid searching in the front for the triangles to be removed from it, these triangles are marked as not available, which prevents them to be made active (in the 1st step).

Contradictory to a 2D (or surface) triangulation, a 3D discretization suffers from a higher percentage of poor quality elements. This is, first of all, typical for a 3D Delaunay triangulation in which high quality triangles may form degenerated tetrahedral elements with volume close (or even equal) to zero. Unfortunately, these elements, called slivers, are likely to be present also in meshes created by the advancing front technique. Since the Delaunay property is used as the governing criterion for the neighbourhood search in the present approach, the probability of a sliver occurrence is even higher. Therefore some postprocessing has to be employed to improve the quality of the mesh. In this case, however, the smoothing process alone is not effective enough. The application of the Laplacian smoothing technique improves the mesh in average but the quality of sliver elements is usually deteriorated even more. It is therefore essential to combine the smoothing with some topologically based technique. In the presented work, two different topological transformations are employed (Fig. 2.26). While the first one removes a flat tetrahedron by swapping a triangle, the second one consists in a nodal reconnection in an octahedron. The actual transformation is accomplished only if its application does not yield a tetrahedron of the worse quality than is the quality of the worst tetrahedron involved in the transformation and if no violation of the topological compatibility of the mesh occurs. The latter condition limits the use of topological transformations near the region boundary. From this point of view, the introduction of ellipsoid in the local neighbourhood search is very favourable because its application tends to produce such element configurations (near the boundary) which can be subjected to the transformation without the lost of the topological compatibility. Note that a node insertion (e.g. in the centroid of the sliver) followed by an appropriate nodal reconnection is another effective way to remove slivers. Nevertheless, also this technique has a limited use near the boundary. The postprocessing then consists of several cycles of the Laplacian smoothing (optionally with a weighting) alternating with the topological transformation process. This technique proves to be effective and yields relatively high quality meshes in only a few cycles (typically about five). The actual implementation of the Laplacian smoothing for the region mesh optimization does not include the relaxation parameter as it was the case in the surface mesh optimization. The position of node being smoothed in the -th iteration can be then expressed as

where are nodes connected to node by an edge and is either equal to (see Eq. (2.85)) for the element size weighting, or to (see Eq. (2.86)) for the connectivity weighting, or to for the combined weighting, or is set to if no weighting is required. Only nodes which are classified to the region are subjected to the smoothing. Similarly as in the smoothing of a surface mesh, the values of and are replaced by and , respectively if the -th iteration of the position of node is already available. Again, a simplified approach for the evaluation of the ``optimal'' valency, required when the nodal connectivity based weighting (Eq. (2.86)) should be applied, has been employed, thus is equal to for nodes classified to the region or to otherwise. Note that the value has been obtained from an assumption that there are twice as many triangles as nodes in a ``closed'' 2D uniform mesh and that nearly equilateral tetrahedrons are required to fill in the full solid angle.

*Daniel Rypl
2005-12-07*