Quality Improvement Algorithm for Tetrahedral Mesh Based on Optimal Delaunay Triangulation

The concept of optimal Delaunay triangulation (ODT) and the corresponding error-based quality metric are first introduced. Then one kind of mesh smoothing algorithm for tetrahedral mesh based on the concept of ODT is examined. With regard to its problem of possible producing illegal elements, this paper proposes a modified smoothing scheme with a constrained optimization model for tetrahedral mesh quality improvement. The constrained optimization model is converted to an unconstrained one and then solved by integrating chaos search and BFGS (Broyden-Fletcher-Goldfarb-Shanno) algorithm efficiently. Quality improvement for tetrahedral mesh is finally achieved by alternately applying the presented smoothing scheme and re-triangulation. Some testing examples are given to demonstrate the effectiveness of the proposed approach.


Introduction
The mesh improvement procedure can be basically divided into two categories, topological optimization and geometrical optimization (also called smoothing).Topological optimization changes the topology of a mesh, i.e. node-element connectivity relationship, while the geometrical optimization or smoothing improves mesh quality by simply moving or adjusting the position of nodes without changing the topology of a mesh.Usually, the optimization-based smoothing procedure would lead to better results.However, due to the complexity of models and algorithms, the optimization-based smoothing is difficult to implement, and with the increase of the number of design variables, the computational cost significantly increases.Therefore, in many works, smoothing was commonly performed by Laplacian smoothing technique, i.e. shifting each interior free node to the center of the surrounding polygon or polyhedron.Although computationally inexpensive, this method, may produce invalid or illegal elements occasionally that are unacceptable in numerical analysis.The topological optimization improves local and whole mesh quality by changing the topology of a mesh, i.e. node-element connectivity rela-tionship.The most frequently used operations of topological optimization are so-called basic or elementary flips.Delaunay triangulation can also be proved to be a topological optimization tool to improve the quality of a mesh.
Chen et al. [1] presented the concept of ODT (Optimal Delaunay Triangulation) in 2004 and developed a kind of quality smoothing algorithm [2] for triangular mesh by minimizing the linear interpolation error-based mesh quality metric.The algorithm can calculate the optimal location of each interior node directly according to the location of its neighbor nodes without solving the optimization model.The computational cost of such algorithm is as low as that of Laplacian smoothing.Combining this smoothing algorithm and the concept of ODT, Alliez et al. [3] proposed a Delaunay-based variational approach for tetrahedral mesh by minimizing a simple mesh-dependent energy through global and updated both vertex positions and connectivity, which is very effective to improve the quality of "bad" tetrahedrons.
However, through numerical investigation and analysis we find that the smoothing formula used [3] for updating vertex positions also suffers the problem of producing illegal elements.To this end, this paper proposes a modified smoothing scheme for tetrahedral mesh which avoids the possibility of producing illegal elements while maintaining high efficiency of the original scheme.The corresponding optimization model is then solved by integrating chaos search and BFGS algorithm efficiently.Quality improvement for tetrahedral mesh is finally achieved by alternately applying the presented smoothing scheme and re-triangulation.
An outline of this paper is as follows.Section 2 introduces the concept of ODT and the corresponding errorbased quality metric.Then in Section 3, smoothing scheme through minimizing the interpolation error among all triangulations with the same number of vertices in the local patch is discussed and tested.A modified smoothing scheme which avoids the possibility of producing illegal elements is proposed in Section 4. Section 5 gives the quality improvement cycle by alternately applying the presented smoothing scheme and re-triangulation.Several examples are given in Section 6 to illustrate performance of the proposed method.

ODT and Error-Based Quality Metric
Chen et al. [1] presented the concept of ODT and proved the following theorem: For a given finite point set V in R n , let Ω be a convex hull of V, T a triangulation of Ω, and f I,T the piecewise linear finite element interpolation of a given continuous function f defined on Ω. Define Then for the Delaunay triangulation of V, DT, where T V denotes all possible triangulations of Ω by using the points in V.This theorem shows that for all the possible traingulations, Delaunay triangulation makes the linear interpolation error   2 , , Q T x p take the minimum.Based on the linear interpolation error defined in above theorem,   , , Q T f p , Chen et al. [2] derived an error-based mesh quality metric, i.e. the energy metric named as where n is dimension (for tetrahedral mesh n = 3) and Ω i denotes the first ring of vertex x i .Chen et al. [2] presented that E ODT achieves the minimum if all edge lengths of the triangulation are equal; If reducing the value of E ODT , the differences between edge length and volume of mesh can also be reduced.Therefore, reducing the value of E ODT by topological or geometrical optimization approach can lead to improvement of the mesh quality.

Smoothing Scheme
When performing mesh quality improvement, the errorbased mesh quality metric defined in the local patch Ω i for free vertex x i is used.The local quality metric for tetrahedral patch Ω i , can be defined according to [2] and derived as where |T j | is the volume of tetrahedron T j with respect to vertex x i , |Ω i | denotes the volume of Ω i .If reducing the value of this local quality metric, the quality of whole mesh will be improved.Then the smoothing scheme can be established by minimizing the local quality metric in Equation ( 4) by the manner of point by point.According to [2,3], this optimization-based smoothing model could be solved explicitly, so that the computational cost is as low as that of Laplacian smoothing.In [3], the optimal position of x i is derived and given by 2 * 1 2 where is the gradient of the volume of the tetrahedron T j with respect to x i .Formula (5) can be used to directly update the location of free vertex x i in tetrahedral mesh through the positions of its neighbors and thus improves the mesh quality.Combining this smoothing algorithm and the concept of ODT, Alliez et al. [3] proposed a Delaunay-based variational approach for tetrahedral meshing by minimizing a simple mesh-dependent energy through global updates of both vertex positions and connectivity, which is very effective to improve the quality of poorly tetrahedrons.However, through numerical investigation and analysis we find Formula (5) for updating position of x i also suffers the deficiency of producing illegal elements.

Algorithm Analysis and Testing
Formula (5) can be derived by making the derivative of i E  zero.However, the feasible domain of x i * should be restricted in a certain region inside Ω i .The position of x i that makes the derivative of i E  zero in this region not always exists.Under some circumstances, especially for poorly mesh, some vertexes after updating by Formula (5) will locate outside of Ω i .Some numerical tests also demonstrate Formula (5) has possibility to lead to illegal elements.For example, for the Delaunay triangulation in a cubic region shown in Figure 1, perform smoothing approach for the interior nodes by Formula (5).As we see, some new positions of vertexes locate far outside the cubic region.

Modified Smoothing Scheme and Corresponding Optimization Model
It is unacceptable for producing illegal elements in smoothing process.There are some strategies to handle this problem.The simplest strategy is to discard the new position of x i if illegal element is detected.We test this strategy with the mesh in Figure 1, and no illegal element is found to be produced.However, for such a modification, the effect of quality improvement is not satisfactory despite of its high efficiency.
The better way may be to directly solve the optimization-based model by minimizing the local quality metric in Equation ( 4) with constraint of feasible domain, which needs to calculate feasible domain i F  first, and then solve a constrained optimization model: To skip the extra calculation of the feasible domain, we prefer to realize the equivalent constraint by introducing the condition of positive volume, i.e. |T j | > 0. Thus we get a modified smoothing scheme and corresponding optimization model as follows: Step 1) Calculate the new position of x i * by Formula (5); Step 2) Check if x i * causes illegal element.A container of Ω i can be first constructed to reduce the checking cost.If x i * locates outside of this container, there must exist illegal element; otherwise, further volume calculation for tetrahedral elements using the new position of x i * is needed.
Step 3) If no illegal element is found to be produced, accept the new position; otherwise, solve the optimization problem min i E  with volume constraints 0 j T  to obtain the new position of x i * .In order to obtain better numerical stability, the local quality metric in Equation ( 4) can be expressed as by coordinate translation.Then the constrained optimization model in Step 3) can be defined by In order to solve this kind of constrained optimization model efficiently, we can convert it to an unconstrained one with differentiable objective function according to [4].First, the volume constraint |T j | > 0 can be converted to the following maximum constraint: Then construct the L  penalty function (10) and replace the non-differentiable term in Equation ( 10) by its aggregate function.Thus the differentiable objective function is obtained as: When the penalty factor α is greater than some threshold value α 0 and the parameter q approaches to infinity, the solution of the unconstrained model for minimizing the differentiable objective function in Equation (11) approaches to the solution of the original constrained model (8).

Solution Algorithm
For consideration of efficiency and precision, the scheme that integrates chaos search and BFGS is employed to solve the optimization model.Chaos search has good global search ability, but its convergence rate is low.BFGS method is one of the most representative algorithms in solving unconstrained optimization problem with good numerical stability.The combination of chaos search and BFGS can make full use of both global search ability of chaos search and fast convergence rate near the optimal solution of BFGS.Based on the above consideration, the solution of optimization problem in Step 3) can be divided into two steps: first obtain the approximated solution by chaos search algorithm as the initial input values; then do optimization by BFGS algorithm.
The procedure of chaos search algorithm is simply described as follows: first map the design variable x i to chaotic variable, then to the objective function i E  in Equation (8) directly search the best solution that satisfies the constraints.The volume of each tetrahedron is calculated before evaluating the objective function; therefore, if the current position does not satisfy the con-straints, directly go to the next searching.
Then using the approximated solution by chaos search algorithm as the initial input values, do optimization to the differentiable objective function in Equation (11) by the BFGS algorithm for optimal solution.
Due to stochastic property of chaos search algorithm, the approximated solutions may be different.In practice calculation, four or five approximated solutions obtained by chaos search are taken as the initial inputs for BFGS algorithm.Then the best solution by BFGS algorithm is chosen as the final optimal solution.

Mesh Quality Improvement
To improve mesh quality, smoothing or topological optimization alone may not achieve ideal results.In general, alternately applying smoothing and topological optimization technique can improve mesh quality significantly [5].Smoothing process changes the location and distribution of nodes, and may provide further space of improving mesh quality for topological process.Therefore, this paper combines these two approaches, alternately applying smoothing and topological optimization to achieve better results.The most frequently used operations of topological optimization are so-called basic or elementary flips.Recently, Liu et al. [6] proposed a new topological optimization operation named small polyhedron reconnection (SPR), to search the optimal topological configura- tion in an enlarged region which usually includes 30 ~ 40 tetrahedrons.Delaunay triangulation can also be proved to be a topological optimization tool to improve the quality of mesh.This paper adopts re-triangulation as topological approach, makes a new Delaunay triangulation after nodes updating by the proposed smoothing scheme.Such an approach can reduce the value of quality metric based on interpolation error further [3] (Among all the triangulations, Delaunay triangulation makes the value take minimum [2]).Meanwhile, reconnection of nodes also provides optimization space for next smoothing step.

Examples and Conclusions
The proposed approach is tested by a number of examples to illustrate its effectiveness.Quality improvement is performed by applying the presented smoothing scheme and re-triangulation for topological optimization alternatively.Only interior nodes are repositioned during smoothing, while the nodes on boundaries keep fixed.Three testing meshes are shown in Figure 2. The first testing mesh consists of 10,926 nodes and 58,615 tetrahedrons initially.The second testing mesh consists of 19,149 nodes and 76,188 tetrahedrons.The third mesh includes 6534 nodes and 25,177 tetrahedrons initially.Table 1 shows the statistics of initial quality and quality after 10  optimization cycles.The radius ratio ρ [5] is adopted as the quality measurement for tetrahedral element, which takes the value of 1 for regular tetrahedron and 0 for degenerated element.
From the optimization results, it can be seen that the quality optimization approach presented is able to improve significantly the average quality of the mesh, even the boundary nodes are fixed.More important, the number of poorly elements decreases remarkably and the worst mesh quality also be improved.
The presented smoothing scheme is efficient and robust.No illegal element is produced.Testing results show that the proposed smoothing approach is suitable for combining with the topological optimization procedure, and effective to eliminate poor elements as well as improve mesh quality.

Figure 1 .
Figure 1.Mesh in cubic region before and after smoothing.(a) Initial mesh; (b) New mesh after smoothing.

Figure 2 .
Figure 2. Initial meshes for the three testing examples.(a) The first example; (b) The second example; (c) The third example.