Next: Mesh Generator Evaluation Up: Mesh Generation Previous: Shell Discretization

## 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:

1. 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).

2. 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)

3. 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.

4. 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)

5. 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:

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.
Note that a different strategy (compared to the one used for neighbourhood construction in Section Surface Discretization (Sequential Mesh Generation)) is used in case (c) resulting, under some circumstances, in an increase of the neighbourhood.

6. 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

1. 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 .
2. 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 .
3. taken as node if it falls inside the sphere and there is no other node inside the sphere circumscribed to the tetrahedron .
4. 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 .
5. 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 .
6. 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 .
7. 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 .
8. 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.

7. 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.

8. 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.

9. 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

 (2.99)

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.

Next: Mesh Generator Evaluation Up: Mesh Generation Previous: Shell Discretization

Daniel Rypl
2005-12-07