Intelligent Information Management, 2013, 5, 191-195
Published Online November 2013 (http://www.scirp.org/journal/iim)
http://dx.doi.org/10.4236/iim.2013.56021
Open Access IIM
Quality Improvement Algorithm for Tetrahedral Mesh
Based on Optimal Delaunay Triangulation
Shuli Sun, Haoran Bao, Minghui Liu, Yuan Yuan
State Key Laboratory for Turbulence and Complex System, Department of Mechanics and Engineering Science,
College of Engineering, Peking University, Beijing, China
Received September 20, 2013; revised October 21, 2013; accepted November 15, 2013
Copyright © 2013 Shuli Sun et al. This is an open access article distributed under the Creative Commons Attribution License, which
permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
ABSTRACT
The concept of optimal Delaunay triangulation (ODT) and the corresponding error-based quality metric are first intro-
duced. 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.
Keywords: Tetrahedral Mesh; Mesh Quality Improvement; Smoothing; Topological Optimization; Optimal Delaunay
Triangulation
1. Introduction
The mesh improvement procedure can be basically divided
into two categories, topological optimization and geo-
metrical optimization (also called smoothing). Topologi-
cal optimization changes the topology of a mesh, i.e.
node-element connectivity relationship, while the geo-
metrical optimization or smoothing improves mesh qual-
ity 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 dif-
ficult to implement, and with the increase of the number
of design variables, the computational cost significantly
increases. Therefore, in many works, smoothing was com-
monly performed by Laplacian smoothing technique, i.e.
shifting each interior free node to the center of the sur-
rounding polygon or polyhedron. Although computation-
ally inexpensive, this method, may produce invalid or
illegal elements occasionally that are unacceptable in
numerical analysis. The topological optimization im-
proves 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 topo-
logical 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 opti-
mization model. The computational cost of such algo-
rithm is as low as that of Laplacian smoothing. Combin-
ing 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 analy-
sis we find that the smoothing formula used [3] for up-
dating vertex positions also suffers the problem of pro-
ducing illegal elements. To this end, this paper proposes

S. L. SUN ET AL.
Open Access IIM
192
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 in-
tegrating 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 intro-
duces the concept of ODT and the corresponding error-
based 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 smooth-
ing 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. Se-
veral examples are given in Section 6 to illustrate per-
formance of the proposed method.
2. 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 Rn, let Ω be a convex
hull of V, T a triangulation of Ω, and fI,T the piecewise
linear finite element interpolation of a given continuous
function f defined on Ω. Define


 
,
1
,
,,
d
p
IT L
p
p
IT
QT fpff
fx fxx







 (1)
Then for the Delaunay triangulation of V, DT,



22
,,min ,,,1
V
TT
QDTx pQTx pp

 (2)
where TV denotes all possible triangulations of Ω by us-
ing the points in V.
This theorem shows that for all the possible traingula-
tions, Delaunay triangulation makes the linear interpola-
tion error

2
,,QTxp take the minimum.
Based on the linear interpolation error defined in
above theorem,

,,QTfp, Chen et al. [2] derived an
error-based mesh quality metric, i.e. the energy metric
named as EODT in [3],

2 2
1
1
,,1 d
1i
N
ODT i
i
EQTx xxx
n


 (3)
where n is dimension (for tetrahedral mesh n = 3) and Ωi
denotes the first ring of vertex xi.
Chen et al. [2] presented that EODT achieves the mini-
mum if all edge lengths of the triangulation are equal; If
reducing the value of EODT, the differences between edge
length and volume of mesh can also be reduced. There-
fore, reducing the value of EODT by topological or geo-
metrical optimization approach can lead to improvement
of the mesh quality.
3. Mesh Smoothing Based on Interpolation
Error
3.1. Smoothing Scheme
When performing mesh quality improvement, the error-
based mesh quality metric defined in the local patch Ωi
for free vertex xi is used. The local quality metric for
tetrahedral patch Ωi, can be defined according to [2] and
derived as
22
1
44
i
ji kj
ki
i
jk i
TxT
xx
ETxx
 








 (4)
where |Tj| is the volume of tetrahedron Tj with respect
to vertex xi, |Ω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 xi is derived
and given by
2
*1
2i
ji kj
ki
iixj ik
TxT
i
xx
xxT xx
 



 



 (5)
where i
x
j
T is the gradient of the volume of the
tetrahedron Tj with respect to xi.
Formula (5) can be used to directly update the location
of free vertex xi 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 mini-
mizing a simple mesh-dependent energy through glo-
bal updates of both vertex positions and connectivity,
which is very effective to improve the quality of poor-
ly tetrahedrons. However, through numerical investi-
gation and analysis we find Formula (5) for updating
position of xi also suffers the deficiency of producing
illegal elements.
3.2 Algorithm Analysis and Testing
Formula (5) can be derived by making the derivative of
i
E

zero. However, the feasible domain of xi
* should be
restricted in a certain region inside Ωi. The position of xi

S. L. SUN ET AL.
Open Access IIM
193
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.
(a)
(b)
Figure 1. Mesh in cubic region before and after smoothing.
(a) Initial mesh; (b) New mesh after smoothing.
4. Modified Mesh Smoothing Scheme
4.1. Modified Smoothing Scheme and
Corresponding Optimization Model
It is unacceptable for producing illegal elements in smooth-
ing process. There are some strategies to handle this
problem. The simplest strategy is to discard the new po-
sition of xi if illegal element is detected. We test this
strategy with the mesh in Figure 1, and no illegal ele-
ment is found to be produced. However, for such a modi-
fication, the effect of quality improvement is not satis-
factory despite of its high efficiency.
The better way may be to directly solve the optimiza-
tion-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:

22
1
min44
i
ji kj
ki
i
i
ijki
TxT
xx
i
ExT xx
xF
 










 (6)
To skip the extra calculation of the feasible domain,
we prefer to realize the equivalent constraint by intro-
ducing the condition of positive volume, i.e. |Tj| > 0.
Thus we get a modified smoothing scheme and corre-
sponding optimization model as follows:
Step 1) Calculate the new position of xi
* by Formula
(5);
Step 2) Check if xi
* causes illegal element. A container
of Ωi can be first constructed to reduce the checking cost.
If xi
* locates outside of this container, there must exist
illegal element; otherwise, further volume calculation for
tetrahedral elements using the new position of xi
* is
needed.
Step 3) If no illegal element is found to be produced,
accept the new position; otherwise, solve the optimiza-
tion problem min i
E

with volume constraints
0
j
T to obtain the new position of xi
*.
In order to obtain better numerical stability, the local
quality metric in Equation (4) can be expressed as

2
1
4
i
ji kj
ki
jik
TxT
xx
i
xETxx
 







 (7)
by coordinate translation. Then the constrained optimiza-
tion model in Step 3) can be defined by

2
1
min 4
s. t. 0
i
ji kj
ki
ijik
TxT
xx
j
ExT xx
T
 








 (8)

S. L. SUN ET AL.
Open Access IIM
194
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 |Tj| > 0 can be converted
to the following maximum constraint:

2
1
max 0
4kj
ki
jii k
xT
xx
Txx x








 (9)
Then construct the L penalty function
 
2
1
,max0,
4
i
kj
ki
ii jiik
xT
xx
xExTx xx  




 




(10)
and replace the non-differentiable term in Equation (10)
by its aggregate function. Thus the differentiable objec-
tive function is obtained as:
 

2
,
ln 1exp4
i
kj
ki
ii
jii k
jxT
xx
xEx
qTxx x
q










 







(11)
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 dif-
ferentiable objective function in Equation (11) approaches
to the solution of the original constrained model (8).
4.2. 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 algo-
rithms 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 op-
timization 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 opti-
mization by BFGS algorithm.
The procedure of chaos search algorithm is simply de-
scribed as follows: first map the design variable xi to
chaotic variable, then to the objective function i E

in
Equation (8) directly search the best solution that satis-
fies the constraints. The volume of each tetrahedron is
calculated before evaluating the objective function; there-
fore, 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.
5. Mesh Quality Improvement
To improve mesh quality, smoothing or topological op-
timization alone may not achieve ideal results. In general,
alternately applying smoothing and topological optimiza-
tion 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 pa-
per combines these two approaches, alternately applying
smoothing and topological optimization to achieve better
results. The most frequently used operations of topologi-
cal optimization are so-called basic or elementary flips.
Recently, Liu et al. [6] proposed a new topological opti-
mization operation named small polyhedron reconnec-
tion (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 qual-
ity of mesh. This paper adopts re-triangulation as topo-
logical 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.
6. Examples and Conclusions
The proposed approach is tested by a number of exam-
ples to illustrate its effectiveness. Quality improvement is
performed by applying the presented smoothing scheme
and re-triangulation for topological optimization alterna-
tively. 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 tetra-
hedrons initially. The second testing mesh consists of
19,149 nodes and 76,188 tetrahedrons. The third mesh in-
cludes 6534 nodes and 25,177 tetrahedrons initially. Table
1 shows the statistics of initial quality and quality after 10

S. L. SUN ET AL.
Open Access IIM
195
(a) (b) (c)
Figure 2. Initial meshes for the three testing examples. (a) The first example; (b) The second example; (c) The third exam-
ple.
Table 1. Statistics of initial mesh quality and the quality after optimization.
Initial mesh After 10 optimization cycles
min Average Number of elements with
quality below 0.01 min Average Number of elements with quality below 0.01
Example 1 5.74E−6 0.795 66 0.062 0.899 0
Example 2 2.92E−7 0.804 11 0.081 0.873 0
Example 3 0.008 0.794 11 0.035 0.865 0
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 de-
generated element.
From the optimization results, it can be seen that the
quality optimization approach presented is able to im-
prove significantly the average quality of the mesh, even
the boundary nodes are fixed. More important, the num-
ber of poorly elements decreases remarkably and the worst
mesh quality also be improved.
The presented smoothing scheme is efficient and ro-
bust. 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 im-
prove mesh quality.
7. Acknowledgements
The authors would like to thank the supports of the Na-
tional Natural Science Foundation of China (11172004,
10772004) and Beijing Municipal Natural Science Foun-
dation (1102020).
REFERENCES
[1] L. Chen and J. Xu, “Optimal Delaunay Triangulations,”
Journal of Computational Mathematics, Vol. 22, No. 2,
2004, pp. 299-308.
[2] L. Chen, “Mesh Smoothing Schemes Based on Optimal
Delaunay Triangulations,” Proceedings of 13th Interna-
tional Meshing Roundtable, Williamsburg 19-22 Sep-
tember 2004, pp. 109-120.
[3] P. Alliez, D. Cohen-Steiner, M. Yvinec and M. Desbrun,
“Variational Tetrahedral Meshing,” ACM Transactions on
Graphics, Vol. 24, No. 3, 2005, pp. 617-625.
http://dx.doi.org/10.1145/1073204.1073238
[4] X. S. Li, “An Efficient Approach to a Class of Non-
Smooth Optimization Problems,” Science in China Series
A, Vol. 24, No. 4, 1994, pp. 371-377.
[5] S. L. Sun and J. F. Liu, “An Efficient Optimization Pro-
cedure for Tetrahedral Meshes by Chaos Search Algo-
rithm,” Journal of Computer Science and Technology,
Vol. 18, No. 6, 2003, pp. 796-803.
http://dx.doi.org/10.1007/BF02945469
[6] J. F. Liu, S. L. Sun and D. C.Wang, “Optimal Tetrahe-
dralization for Small Polyhedron: A New Local Trans-
formation Strategy for 3-d Mesh Generation and Mesh
Improvement,” CMES: Computer Modeling in Engineer-
ing & Sciences, Vol. 14, No. 1, 2006, pp. 31-43.