In order to improve the boundary mesh quality while maintaining the essential characteristics of discrete surfaces, a new approach combining optimization-based smoothing and topology optimization is developed. The smoothing objective function is modified, in which two functions denoting boundary and interior quality, respectively, and a weight coefficient controlling boundary quality are taken into account. In addition, the existing smoothing algorithm can improve the mesh quality only by repositioning vertices of the interior mesh. Without destroying boundary conformity, bad elements with all their vertices on the boundary cannot be eliminated. Then, topology optimization is employed, and those elements are converted into other types of elements whose quality can be improved by smoothing. The practical application shows that the worst elements can be eliminated and, with the increase of weight coefficient, the average quality of boundary mesh can also be improved. Results obtained with the combined approach are compared with some common approach. It is clearly shown that it performs better than the existing approach.

1. Introduction

Numerical simulation is an important component in diverse activities such as medical imaging, engineering design, and cinematic special effects. These simulations routinely rely on tetrahedral meshes to model the physical domain of interest. The popularity of tetrahedral meshes stems from their ability to accurately model extremely complex geometries. High resolution tetrahedral meshes are often particularly desirable, as they can improve numerical accuracy greatly at previously infeasible scales. Fortunately, the power of modern computing and data acquisition technologies has enabled the production of tetrahedral models with enormous size. However, it is normally difficult to acquire mesh in which all the elements are suitable for numerical computation. Poorly shaped tetrahedra in a mesh can result in numerical errors and increase the time cost to find a solution [1, 2]. Hence, there is a market for mesh improvement tools which can enhance the quality of tetrahedra for an existing mesh.

Recently, there are two main techniques to improve the mesh quality. The first modifies topology by inserting or deleting nodes as well as changing connectivity of nodes [3, 4]. The operations include local face swapping and element or vertex insertion/deletion. The second, called smoothing method [5, 6], preserves mesh topology by applying appropriate node placement techniques. Mesh improvement techniques have been shown to be effective in improving tetrahedral meshes. However, most of them regard the boundary configuration of a mesh as “untouchable.” For example, vertices on the boundary cannot be smoothed. Moving the vertices on the boundary of a mesh is fraught with peril. In engineering applications, the mesh is considered to be noise-free and any change to the surface may induce some errors in the computational domain. For this, the quality improvement algorithms usually keep boundary vertices locked in place, which limits their ability to improve the quality of tetrahedra.

As to improving boundary quality of a surface mesh, considerable research has been conducted. Garimella et al. [7, 8] proposed an optimization-based vertex repositioning procedure to improve the quality of a surface mesh containing triangles and quadrilaterals. It has been shown that the method is capable of keeping the nodes on the original mesh faces and close to their original locations. Semenova and Savchenko [9] presented two novel techniques to improve the quality of triangle surface meshes while preserving surface characteristics as much as possible. Frey and Borouchaki [10] proposed a suitable method for the construction of an enriched geometric finite element mesh from a given arbitrary surface triangulation. The initial triangulation is optimized with respect to geometry and element shape quality. Actually, the surface mesh quality has an important effect on high quality solid mesh generation. Surface mesh quality improvement algorithms can reduce the difficulty of generating solid meshes and improving mesh quality. However, when it comes to numerical simulation, they cannot improve its accuracy efficiently.

In this study, we are more concerned with optimizing boundary vertices to obtain better quality tetrahedra, rather than surface meshes. However, little research has been done to improve boundary surface quality in solid meshes. Klingner [11] presented a quadric smoothing method for smoothing boundary vertices of surfaces in solid meshes. It permits all surface vertices to move, but they are encouraged to move along the original surface, and discouraged from making noticeable changes to the shape of the domain. By balancing tetrahedron quality against a quadric error measured at each vertex, the method controls the domain shape error. However, it cannot preserve mesh nodes on the original discrete surface and cannot guarantee boundary conformity in a real sense.

The goal of this paper is to improve the quality of boundary tetrahedra. It should be noted that, unlike Klingner’s work, surface vertices will not be moved. Topology optimization technique is adopted for converting elements that have all their vertices on the boundary into other type of elements, and the optimization-based smoothing algorithm is used for improving mesh quality. Section 2 briefly describes the overall combined algorithm. Section 3 introduces the modified optimization-based smoothing algorithm and Section 4 discusses the topology optimization technique. We then present the results of numerical experiments on several test meshes. Finally, Section 5 concludes the paper.

2. The Overall Combined Optimization Algorithm

Due to inherent defect, the optimization-based smoothing algorithm cannot improve boundary mesh quality through moving boundary nodes without distorting discrete surface. The boundary conformity cannot be guaranteed when surface nodes are permitted to move. For example, the elements in Figures 1(a), 1(b), and 1(c) might be improved by moving a node if the node is an interior one. However, for the element in Figure 1(d), there are no interior nodes to move to improve its quality. It is called “all boundary element” in the paper. The quality improvement for this type of element with all their nodes on the boundary surface is beyond the capability of smoothing. Therefore, there are a lot of poor quality elements on the boundary after performing the optimization-based smoothing algorithm. To address this problem, a combined optimization algorithm based on optimization-based smoothing and topology optimization algorithm is proposed. The algorithm can improve the overall and boundary mesh quality and won’t destroy boundary integrity of the initial mesh.

The overall scheme is presented. Algorithmic details for each of the major steps in this scheme will be presented later in the paper.

Input: initial mesh, poor quality threshold value , iteration number , objective function , and weight coefficient .

Output: high quality mesh.(1)Compute the number of nodes and elements in initial mesh. If the element is a boundary one, mark the element with the boundary;(2)calculate the mesh quality of initial mesh based on the quality measure (see Section 3.1), and store the nodes of poor quality element in , whose quality is less than ;(3)in , select a node randomly and determine whether it is a boundary node. If the node is a boundary one and all its neighboring nodes are on the boundary, go to Step (4). Otherwise, carry out Step (5).(4)Apply Topology optimization algorithm (see Section 4) to convert elements that have all their vertices on the boundary into other type of elements, and go to Step (6).(5)Apply improved optimization-based smoothing algorithm (see Section 3) to optimize mesh. (6)Repeat Steps (3) to (5) until there are no poor quality elements in .

3. The Improved Optimization-Based Smoothing

3.1. Objective Function

Through studying the quality measures and optimization algorithms of tetrahedral meshes, the error function is derived by means of transforming the measures. The element’ distortion is disposed as the error. And the larger the distortion is, the bigger the value of the error. The maximum value of the error function is infinite, and the minimum is zero. The total error function of the mesh is the sum of the elements’ errors, and it is adopted as the objective function of the optimization-based smoothing. The minimum value of objective function is solved for improving the mesh quality.

Generally speaking, a reasonable quality measure for elements should possess the following attributes.(1)It is invariant under translation, rotation, and scaling.(2)Normalization by an optimal value within a range [0,1], where 1 is for an equilateral tetrahedron and 0 is for a degenerate tetrahedron.(3)Ability to detect all possible badly shaped elements.

So an optimal quality measure is employed, and an error function is derived based on it: where is the volume of tetrahedron, is the surface area of a triangular facet, and is the length of any edge .

The reciprocal of (3.1) is the error function of an element: where denotes a tetrahedron.

When the element is a regular tetrahedron, the value of the function is 1. As one element degenerates, the value tends to be infinite. If the element is reverse or the volume is negative, the value of error is also infinite. The total error function of the mesh is defined as , as shown in the following:

is the objective function of the optimization-based smoothing. In order to control the boundary and interior mesh quality simultaneously, the objective function is expressed as follows: where is the set of all elements, is the set of all interior elements, is the set of all boundary elements, and () is the weight coefficient which indicates the proportion of boundary mesh quality with respect to the overall. The same method is executed for boundary and interior mesh when . Through (3.4), the interior and boundary mesh quality is considered simultaneously.

3.2. Optimization Algorithm

An efficient and robust solver for the large system of equations presented by the optimization problem is needed. Better smoothing algorithms are based on numerical optimization [12, 13]. Early algorithms define a smooth objective function that summarizes the quality of a group of elements (e.g., the sum of squares of the qualities of all the tetrahedra adjoining a vertex) and use a numerical optimization algorithm such as steepest descent or Newton’s method to move a vertex to the optimal location. Freitag et al. [14] proposes a more sophisticated nonsmooth optimization algorithm, which makes it possible to optimize the worst tetrahedron in a group, for instance, to maximize the minimum dihedral angle among the tetrahedra that share a specified vertex. A nonsmooth optimization algorithm is needed because the proposed objective function is not a smooth function of the vertex coordinates in this paper. The gradient of this function is discontinuous wherever the identity of the worst tetrahedron in the group changes. Pseudopodia for nonsmooth optimization algorithm are presented in Algorithm 1.

  function SOLUTION ,  ,  ,  , 
     ▸   = objective function
      = initial step size
      desired degree of accuracy
      iteration number
      initial coordinate for vertices X
 get free vertex
  ▸   = Hessian matrix
   Iteration number
  ▸   optimal search direction
 compute using line search algorithm (see Algorithm 2)
  ▸   = optimal step size
 if   or then
  ▸   optimal point coordinate
  ▸   gradient at vertices
   go to step (3)

It is well known that the line search methods play a pivotal role on optimization problems. Locating a local minimum in the optimization problem with no constrains are prepared. All methods have the basic structure in common. In each iteration, a direction is chosen from the current location . The next location, , is the minimum of the function along the line that passes through in the direction . The line search algorithm is conducted in Algorithm 2. In order to find the minimum of a function , we need to bracket it. To bracket a minimum means fining a triple , , so that and . This indicates that the minimum is in the interval []. The interval search algorithm is in Algorithm 3.

  function LINE_SEARCH ( , ( ( ), , ,   ))
     ▸  f( ) = objective function
      = allowable error
      (f(x), , , ) = the function for search region
                  , see below Algorithm 3
(2)  compute = + 0.382( ) and = + 0.618( )
      ▸   = initial tentative point
      Set i = 1
      i = iteration number
  if < then
(4)   return and break
     ▸   = optimal step size
  else if (f( ) < f( )) go to step (7)
(6)          else   go to step (8)
  set = ,   = ,   = ,  f( ) = f( )
   compute = 0.618 + 0.382 and f( )
            i = i+1 go to step (3)

  function SEARCH_INTERVAL (f( ), , , )
     ▸ f(x) = objective function
      = initial coordinate for vertices
      = initial step size
      = coefficient larger than 1
compute = f( ) and Set j = 0
     ▸ j = iteration number
    = + and compute  f( )
 if   then go to step (6)
 else go to step (7)
    = ,   = ,   =
      f( ) = f( ),  j = j + 1
 if j = 0 then = , = go to step (3)
 else a = min , ,  b = max ,
 return and

4. Topology Optimization

Topology optimization algorithm is employed, which converts boundary elements that have all their nodes on the boundary into other types of boundary elements. Furthermore, the topology optimization does not necessarily improve the quality of elements, but those boundary elements must be eliminated. Pseudopodia for topology optimization are presented in Algorithm 4.

Find all the elements which include the node in the all boundary elements, and denotes
    those elements.
  ▸ : Initial all boundary elements
Edge removal can be performed if an edge in the is shared by three tetrahedra.
Edge Swapping can be performed If an edge in the is shared by four or more tetrahedra.
Face Swapping can be performed If a tetrahedron and share a face
end function

4.1. Edge Removal

Proposed by George [15], edge removal is a topological transformation that removes a single edge from the mesh, along with all the tetrahedra that include it. Figure 2 illustrates replacing three tetrahedra with two. The three tetrahedra (ABCD, ABCE and ACDE) share the edge AC, and ABCD tetrahedron is all boundary element. When AC is removed, the three tetrahedra in the region are transformed to two tetrahedra ABDE and BCDE.

4.2. Edge Swapping

Edge swapping is a more complicated procedure [4]. If there are k tetrahedra containing an interior edge e in the mesh, then e is removed and the original k tetrahedra are replaced with 2k−4 tetrahedra. In Figure 3, the edge TB is perpendicular to the page, and the five tetrahedra (01BT, 12BT, 23BT, 34BT, and 40BT) originally incident to edge TB can be replaced by six new tetrahedra: 012T, 021B, 024T, 042B, 234T, and 234B, where 01BT is all boundary element.

4.3. Face Swapping

Face swapping changes the local connectivity of a simplified mesh which converts all boundary elements into other boundary elements. Each interior face in a tetrahedral mesh separates two tetrahedra consisting of a total of five points. A large number of nonoverlapping tetrahedral configurations can be formed with these five points, but only two of them can be reconnected satisfactorily. Figure 4 shows the face swapping configuration of five points. It is evident that either two or three tetrahedra can be used to fill the convex hull of a set of five points. Switching from two to three tetrahedra requires an additional edge interior to the convex hull. So all boundary elements can be converted into other type of boundary elements.

5. Numerical Experiments

Figure 5 shows the distribution of poor quality elements for an impeller in centrifugal pump before and after optimization. The quality for the worst element is 0.08 for initial mesh, 0.15 for , and 0.16 for based on the quality measure . So the improved optimization-based smoothing algorithm can improve the quality of the worst elements in a centrifugal impeller. Furthermore, through increasing the weight coefficient, the quality of boundary and overall elements is improved.

Figure 6 shows the effect of weight coefficient on optimization time. After increasing the weight coefficient , the optimization time for the overall and boundary mesh ascends dramatically. That is because, to obtain optimal solution of objective function, we need to reduce the value of objective function through adjusting interior nodes’ locations. So, it needs more iteration number (i.e., costing more time) when the value of objective function is larger.

The conclusion can be drawn that the proposed algorithm can dramatically reduce the number of poor quality elements in boundary and interior meshes. With the increase of weight coefficient, the overall and boundary mesh quality rises, while the time consumed increases rapidly. So in consideration of optimizing effect and efficiency, it is recommended the weight coefficient lies in the range of 0.7~0.8.

Figure 7 shows two examples to illustrate the effect of different mesh quality improvement algorithms. In Table 1, the number of nodes and elements in each mesh, and the initial mesh quality measured by are given.

In Table 2, the results of each mesh quality improvement technique in terms of the quality of the overall and boundary mesh after optimization are reported. The time consuming is also given. It can be seen that, for the two algorithms, the quality of overall and boundary elements can be improved. And, the combined algorithm performs better in the quality improvement, while it consumes 1.3 times more time than the improved smoothing algorithm.

To evaluate our mesh improvement algorithm, comparison with Freitag and Ollivier-Gooch's method is made. And two cases of mesh (TIRE and RAND2) come courtesy of Freitag and Ollivier-Gooch [4]. TIRE is a tire incinerator with 2,570 vertices and 11,099 tetrahedral elements. Its initial mesh quality can be found in Table 3. RAND2 are lazy triangulations generated by inserting randomly located vertices into a cube one by one. Each vertex is inserted by splitting one or more tetrahedra into multiple tetrahedra. The random meshes have horrible quality and poor dihedral angles at both extremes. The RAND2 mesh has 5,086 vertices and 25,704 tetrahedral elements; its initial quality can be found in Table 3. Table 3 compares the minimum and maximum dihedral angles reported by Freitag and Ollivier-Gooch to that achieved by our mesh improvement algorithm. It can be seen that the mesh quality is bad before optimization, and there are badly shaped elements whose dihedral angles tend to be 0° or 180°. Dihedral angles are improved to be between 14° and 160° for Freitag and Ollivier-Gooch’s algorithm with both swapping and smoothing, and between 23° and 141° for our proposed algorithm.

The two mesh optimization algorithms have been implemented and tested for RAND2 with the distribution of dihedral angles, as shown in Figure 8. Figure 8(a) shows the initial RAND2 mesh, and Figures 8(b) and 8(c) show the optimized mesh by Freitag algorithm and the proposed algorithm, respectively. It can be seen that they can both successfully eliminate poorly shaped elements from the mesh. Comparing with Freitag algorithm, the proposed algorithm is more successful in eliminating poor dihedral angles at both extremes.

6. Conclusions

A new combined algorithm based on optimization-based smoothing and topology optimization is proposed for boundary mesh quality improvement. Error functions of elements determined based on an inverse of measure compose the objective function of a modified optimization-based smoothing. In the objective function, the function terms considering boundary and interior mesh quality are included, and a weight coefficient controlling boundary mesh quality is added, too. By means of minimizing an error function in local iterations, the modified smoothing algorithm can eliminate poor quality elements and improve the overall and boundary mesh quality as much as possible. Topology optimization is employed, by which those bad elements with all their vertices on the boundary are converted into other types of elements whose quality can be improved by the modified smoothing algorithm. Through comparison to Freitag and Ollivier-Gooch’s method with RAND2 mesh used as example, the combined algorithm is confirmed to be capable of improving the boundary mesh quality considerably without distorting discrete surface. With the increase of weight coefficient, the overall and boundary mesh quality rises, while the improvement time increases rapidly. To strike a balance between optimizing effect and time, the weight coefficient is suggested in the range of 0.7~0.8.


This work was supported by The National Natural Science Foundation of China (no. 50825902, 51079062, 51109095, and 51179075), a Project funded by the Priority Academic Program Development of Jiangsu Higher Education Institutions, Natural science fund in Jiangsu Province (BK2010346, BK2009006), and Postgraduate Innovation Foundation of Jiangsu Province (CXLX11_0576).