Abstract
We propose a derivative-free mesh optimization algorithm, which focuses on improving the worst element quality on the mesh. The mesh optimization problem is formulated as a min-max problem and solved by using a downhill simplex (amoeba) method, which computes only a function value without needing a derivative of Hessian of the objective function. Numerical results show that the proposed mesh optimization algorithm outperforms the existing mesh optimization algorithm in terms of improving the worst element quality and eliminating inverted elements on the mesh.
1. Introduction
It is well known that mesh qualities affect both the accuracy and efficiency of partial differential equation (PDE) solutions [1, 2]. Elements with poor mesh qualities (also, inverted elements) often occur during mesh generation [3], mesh deformation [4, 5], and mesh optimization [5]. Researchers have proposed many mesh quality improvements and untangling algorithms to solve these problems [5–9]. Mesh smoothing algorithms improve mesh qualities simply by moving vertex positions, while fixing element connectivities. One of the most popular mesh smoothing algorithms is Laplacian smoothing [10], where new vertex position becomes a geometric center of adjacent vertices. However, Laplacian smoothing does not guarantee generating good element qualities and instead often generates inverted elements [11]. Several variants of Laplacian smoothing, such as smart Laplacian smoothing, have been proposed to overcome the limitations of Laplacian smoothing [12, 13].
Optimization-based mesh quality improvement algorithms [2, 11, 13] are now becoming more popular, because these methods are able to offer high-quality meshes though with high computational costs. These methods formulate the nonlinear objective function over the entire mesh domain and improve mesh qualities by minimizing the objective function, while keeping mesh topologies [5, 14]. On the other hand, topological change methods improve mesh qualities by iteratively performing edge splitting, edge swapping, and edge collapsing [15, 16]. Researchers have pointed out that mesh smoothing is preferred to topological changes for many PDE-based applications, since it is able to improve element qualities, while keeping mesh topologies [5]. For many PDE problems, consistent mesh topology is desirable for PDE solution accuracy [5, 17]. Several optimization-based simultaneous untangling and smoothing algorithms have been proposed [5–7]. These methods are known to be able to simultaneously eliminate inverted elements, while improving mesh qualities.
Most previous mesh quality improvement methods have focused on improving the average element qualities on the mesh. However, it is well known that a few poor quality elements significantly diminish the efficiency and accuracy of PDE solutions [1, 7]. The inverted elements in the mesh can even result in invalid PDE solutions [4, 9]. In this paper, we focus on improving the worst element quality in the mesh and also focus on eliminating inverted elements in the mesh. The optimization problem, which improves the worst quality elements, is formulated as a nonlinear optimization problem, which is nonsmooth, and thus, the computation of derivative or Hessian is not available. For this reason, traditional nonlinear optimization solvers such as conjugate gradient, steepest descent, Newton, quasi-Newton, and trust-region methods are not directly applicable for solving these nonsmooth mesh optimization problems [7, 13]. Several derivative-free mesh optimization algorithms have been proposed [7, 12, 18]. For example, Freitag and Plassmann proposed an active set algorithm for solving mesh optimization problems [12]. This method improves element qualities by maximizing the minimum area of an element and solves the problem using linear programming. However, this method requires the mesh quality to be convex. Also, other derivative free mesh optimization algorithms need many initial parameters to be chosen or are quite complex to use in practice [7, 18].
We propose an efficient derivative-free mesh optimization algorithm for mesh quality improvement and untangling. Our proposed method employs a downhill simplex method for improving the worst element quality and eliminating inverted elements on the mesh. Our method does not require the computation of a derivative or Hessian of the objective function, but it iteratively performs reflection, expansion, and contraction, to reach the local optimal point (vertex location). Our algorithm is simple to apply and does not need many initial parameters to be determined. We will show numerically that our proposed mesh optimization algorithm outperforms traditional mesh optimization algorithms in terms of worst element quality and speed.
2. Problem Formulation
We mathematically formulate mesh quality metrics and objective functions to optimize. We also describe the local mesh optimization method we used.
2.1. Mesh Quality Metric
There are many mesh quality metrics for triangles and tetrahedra in the literature. See [1] for more details. Among them, we choose an inverse mean ratio (IMR) quality metric to improve the element shape [19]. The IMR quality metric defines an ideal (reference) element and measures how similar the actual (physical) element to the reference element is. The reference triangle is often considered to be an equilateral triangle for isotropic PDE problems. Let the coordinates of three vertices in the triangle element be . The incidence matrix, , of this (physical) element is defined asLet the incidence matrix for the reference element be . For the isotropic PDE problem, the ideal element is considered to be an equilateral triangle or tetrahedron. For the equilateral triangle, is defined as The IMR quality metric is defined as where is a Frobenius norm. When the actual element and the ideal element are identical, is the same as and thus, becomes an identity matrix. In this case, becomes one.
We choose an untangling beta () quality metric to eliminate inverted elements on the mesh [8]. Inverted elements on the mesh have negative signed areas (volume). Untangling beta quality metric gives a high cost (penalty) to the element, which has a negative signed area. Let the area of the triangle be . The untangling beta quality metric is defined as where is a user-defined value, which is close to zero. For a valid (noninverted) element, becomes zero. For the above two quality metrics, a smaller value indicates better element quality.
2.2. Objective Function
We focus on improving the worst element quality on the mesh. Let the th element quality on the mesh be . The element with the worst quality is denoted as where is the number of elements. Then, our goal is to minimize (5). The optimization problem is formulated as Note that the mesh optimization problem in (6) is a nonsmooth objective function. Thus, traditional nonlinear optimization solvers, which require the computation of derivative and/or Hessian, cannot be used to solve it.
2.3. Local Mesh Optimization
In the literature, there are essentially two kinds of mesh optimization strategies: local (Gauss-Seidal type) optimization and global (Jacobi type) mesh optimization. For local mesh optimization, an entire mesh is decomposed into multiple patches. A patch is defined as a set of adjacent elements around the free vertex to optimize, as shown in Figure 1. When local mesh optimization is performed, we visit one free vertex at a time and iteratively move and optimize other vertices. For global mesh optimization, all vertices on the mesh move simultaneously. That is, the entire mesh becomes a single path.
We use local mesh optimization to improve the mesh quality, because previous researches have shown the efficiency of local mesh optimization compared to the global mesh optimization [3, 20]. When global mesh optimization is used, the update of all vertex positions often gets truncated to satisfy mesh validity. Also, the objective function and corresponding matrices in the problem become too large, when the mesh size is too big (e.g., if the mesh includes millions of vertices).
3. A Downhill Simplex Method for Mesh Quality Improvement and Untangling
The optimization problem in (6) is not a smooth objective function due to the maximum function. Thus, traditional nonlinear optimization solvers cannot be used to solve the optimization problem, because the computation of a gradient and Hessian value is not available. Thus, we propose to use a downhill simplex (amoeba) method to solve (6), which does not use either function derivatives or a Hessian but only uses function evaluations.
The basic idea of a downhill simplex method is to move a free vertex by repeatedly performing three actions: reflections, expansions, and contractions. More details on these three actions will be described in the following sections.
3.1. Initial Simplex
The downhill simplex method requires an initial simplex to begin. A simplex is a triangle and tetrahedron in 2D and 3D, respectively. For the local mesh optimization, an initial simplex is generated around the free vertex to optimize by choosing two (virtual) points. It is natural to assume that the initial simplex size should consider the edge length of the mesh or mesh size. We use the minimum edge length information of the patch around the free vertex to determine the initial simplex size. We expect that the initial simplex size should be smaller than the minimum edge length of the patch. Otherwise, the initial simplex location is placed beyond the patch where the free vertex places and two chosen (virtual) points could be too far away from the optimal vertex position. On the other hand, if the initial simplex is too small, the convergence of mesh optimization could be very slow.
Figure 2 shows the initial simplex around the free vertex, . Two (virtual) points, and , and . form an initial simplex. Here, is computed as where the simplex diameter is a user-defined parameter. We will discuss the effect of the simplex diameter on the convergence and the final mesh quality in more detail in Section 4.
3.2. Downhill Simplex Method
In a single iteration of the downhill simplex method, it seeks to eliminate the (virtual) point with the worst cost function and replaces it with a new one by repeatedly performing three different actions: reflection, expansion, and contraction. Let the three points of the current simplex be such thatWe find the vertex with the worst cost function on the simplex (i.e., ) and compute the reflected point of this vertex using the centroid of the current simplex (see Figure 3). The centroid of three points is computed asNow, these three points, including the reflected point, , compose a new simplex. Here, is denoted by We consider three possibilities for this reflected point, . If the function value at the reflected point is neither worst nor best in the new simplex, then replace the worst point by the reflected point. If function value at the reflected point is better than the one at the current point, we perform expansion, since this means that the direction of reflection is good. If function value at the reflected point is worse than the one at the current point, we perform contraction, since this means that the triangle is too large.
Figure 3 shows a current simplex with three points, , , and , reflected point, , expansion point, , and contraction point, . We repeatedly perform the downhill simplex algorithm, which is illustrated in Algorithm 1, until all vertices converge to the local optimum.
| ||||||||||||||||||||||||||||||||||||||||||||||||
4. Numerical Results
In this section, we describe our numerical results to show the effectiveness of our proposed algorithm, which improves the worst element quality on the mesh using the downhill simplex method.
4.1. Test Meshes
We consider four test meshes to evaluate our algorithm. Cylinder and bar meshes (see Figures 4(a) and 4(b)) include many poor quality elements, which were produced during mesh deformation. We use these two meshes to improve the worst element quality on the mesh, because these meshes include very skinny (poor quality) elements. Hydrocephalus and airfoil meshes (see Figures 4(c) and 4(d)) include some inverted elements. The hydrocephalus mesh is provided by Park et al. [21]. These inverted elements were produced during a mesh deformation process.
(a) Cylinder mesh
(b) Bar mesh
(c) Hydrocephalus mesh
(d) Airfoil mesh
Properties of four test meshes and mesh quality statistics are shown in Tables 1 and 2, respectively. The inverse mean ratio quality metric is used to measure the overall mesh quality. The value of the inverse mean ratio quality metric lies between 1 and for valid (noninverted) elements. A smaller value indicates a better mesh quality and the ideal element has the value of one. Elements with negative values indicate inverted elements on the mesh. The initial hydrocephalus and airfoil meshes have 13 and 42 inverted elements, respectively.
4.2. Effect of Initial Simplex Diameter on Algorithm Performance
When the downhill simplex method is used, the initial simplex is composed of one free vertex to optimize and two (virtual) points as discussed in Section 3.1. Other than the free vertex to optimize, the location of the two points is determined by the initial simplex diameter and the minimum edge length in the patch. For example, if simplex diameter is 1, the initial simplex size becomes a minimum edge length in the patch where a free vertex belongs. We investigate the effect of the initial simplex diameter on mesh optimization performance. Here, the mesh optimization performance means both the mesh quality after performing mesh optimization and the time to optimize (time to converge).
Figure 5 shows the worst element quality and timing to optimize as a function of various initial simplex diameters for test meshes. For these meshes, we observe that the best initial simplex diameters lie between 0.01 and 0.1. For this initial simplex size, the output meshes have the smallest worst element quality with the fastest convergence rate. If the initial simplex diameter is too small (e.g., simplex diameter = 0.001), we observe that the convergence time is very slow. For this case, the initial simplex size is too small compared with the mesh size and thus, the movement of one step of the downhill simplex method is very small. If the initial simplex diameter is too big (e.g., simplex diameter = 5), the optimization problem fails to converge or fails to untangle inverted elements. For these cases, we expect that the initial simplex size is chosen to be too big compared with the element size or edge lengths on the mesh. Thus, we chose . Table 3 shows the minimum and maximum edge lengths on the mesh and the best initial simplex diameter. For all test meshes, the best simplex diameters were placed between the minimum and the maximum edge lengths on the mesh.
(a) Cylinder mesh
(b) Cylinder mesh
(c) Bar mesh
(d) Bar mesh
(e) Hydrocephalus mesh
(f) Hydrocephalus mesh
(g) Airfoil mesh
(h) Airfoil mesh
4.3. Comparison with Existing Mesh Quality Improvement Methods
We compare our proposed method with a traditional (existing) mesh optimization method that focuses on improving the average element quality on the mesh. The traditional mesh optimization method, which has been used in many previous research articles [2, 5–7], is to minimizeor to minimize , where is number of elements. We used the nonlinear conjugate gradient (NLCG) method to minimize (11), because many previous papers have shown the effectiveness of using the NLCG method to minimize (11). Table 4 summarizes both our proposed and traditional mesh optimization methods.
The NLCG method updates the vertex position using a line search technique. Let be a vector of vertex positions at the th iteration and let be the search direction. Then, the updated vertex position, , is computed bywhere is a step length. The search direction, , is denoted by where is a parameter given by for Polak-Ribière NLCG method. See [5, 14] for more details about the nonlinear conjugate gradient method for solving mesh optimization problems. Note that the NLCG method is not applicable to solve our mesh optimization problem, (6), which improves the worst element quality on the mesh, since NLCG requires the gradient computation of the objective function.
Figures 6(a) and 6(b) show the worst element quality with respect to the number of iterations of mesh optimization. We observe that our proposed method, which improves the worst element quality using the downhill simplex method, outperforms the traditional method, which improves the average element quality on the mesh. Our method improves the worst element quality up to 54.6% compared with the traditional mesh optimization method. Figure 7 shows zoomed-in mesh, where the element with the worst element quality is located. The initial mesh has the element with the worst element quality in the middle (see a skinny red element enclosed by the circle). This element is very skinny and the element quality is poor with 678.9248 when the inverse mean ratio (IMR) quality metric is used to measure the element quality. Note that a smaller value indicates a better element quality for the IMR quality metric and the ideal element has value of one. After performing the downhill simplex method to improve the worst element quality, these skinny triangles are eliminated and the worst element quality becomes 3.73 (see Figure 7(b)). The traditional mesh optimization method using NLCG also improves the worst element quality but shows inferior performance compared with the proposed method. The output mesh quality using the proposed algorithm is 54.6% better than the output mesh using the traditional method for the cylinder mesh. The output mesh using the traditional method shows the worst element quality with 10.25 (see red elements in Figure 7(c)).
(a) Cylinder mesh
(b) Bar mesh
(c) Hydrocephalus mesh
(d) Airfoil mesh
(a) Initial cylinder mesh (zoomed-in)
(b) Optimized cylinder mesh using downhill simplex method (zoomed-in)
(c) Optimized cylinder mesh using traditional method (zoomed-in)
Figures 6(c) and 6(d) show the number of inverted elements on the mesh as a function of number of iterations of mesh optimization. We observe that the proposed mesh optimization method takes fewer iterations to untangle inverted elements compared with the traditional mesh optimization method. For these numerical experiments, the untangling beta quality metric was used to eliminate inverted elements on the mesh. When the initial mesh includes inverted elements, these elements become the worst quality elements on the mesh. For the hydrocephalus mesh, there are several inverted elements, which are very skinny (see Figure 8(a)). These elements are eliminated after performing mesh optimization (see Figures 8(b) and 8(c)). Both output meshes have no inverted elements, but the proposed method results in better worst element qualities than the traditional method.
(a) Initial hydrocephalus mesh (zoomed-in)
(b) Optimized hydrocephalus using the proposed method (zoomed-in)
(c) Optimized cylinder mesh using the traditional method (zoomed-in)
4.4. Comparison with a Log-Barrier Interior Point Method
We compare our derivative-free mesh optimization method with a log-barrier interior point method (simply, log-barrier method) [7]. The log-barrier method is a popular derivative-free mesh optimization method, which focuses on improving the worst element quality on the mesh. The log-barrier method is flexible in that it can be used to untangle inverted elements on the mesh [7]. Different from other mesh quality optimization methods, it uses the log-barrier objective function to optimize meshes, which is defined as where and are auxiliary parameters and is the number of elements. A quantity is smaller than and is greater than 0. The log-barrier method maximizes with respect to until convergence. Readers refer to [7] for more details on the log-barrier method. Figure 9 shows comparison between the proposed method and the log-barrier method with respect to the worst element quality on the mesh. Note that some parameters of the log-barrier method could not be fully optimized. We observe that both methods significantly improve the worst element quality. However, the proposed method improves the worst element quality up to 75.0% better than the log-barrier method.
We also compare our derivative-free mesh optimization method with the log-barrier method for the initially tangled meshes (see Figures 4(c) and 4(d)). For the tangled airfoil mesh shown in Figure 4(c), our method takes 3.1 seconds to untangle inverted elements while the log-barrier method takes 323.2 seconds to untangle inverted elements. For the tangled hydrocephalus mesh shown in Figure 4(d), both methods take approximately 3 seconds to untangle inverted elements.
5. Conclusions
We have proposed a derivative-free mesh optimization algorithm, which improves the worst element quality on the mesh. We have formulated the mesh optimization problem as a nonlinear optimization problem and solved it using a downhill simplex (amoeba) method. The downhill simplex method was able to solve the nonlinear optimization problem by computing just a function value without needing a derivative or Hessian of the objective function. We have compared our proposed method with the existing mesh optimization algorithms, in terms of the worst element quality and the convergence time.
Numerical results show that our proposed algorithm is more effective in improving the worst element quality on the mesh compared with the one using the existing mesh optimization algorithms. Specifically, the worst element quality using the proposed algorithm is up to 75.0% better than the existing mesh optimization methods. Moreover, the proposed algorithm is more efficient and faster in eliminating inverted elements on the mesh compared with the existing mesh optimization algorithms.
In this paper, we have focused on improving only one aspect of the mesh: element shape. We plan to apply our derivative-free mesh optimization algorithm on simultaneously improving more than two aspects of the mesh. We also plan to parallelize the downhill simplex algorithm using an OpenMP library.
Conflict of Interests
The authors declare that there is no conflict of interests regarding the publication of this paper.
Acknowledgments
The authors were supported by Basic Science Research Program through the National Research Foundation of Korea (NRF) funded by the Ministry of Science, ICT and Future Planning (NRF-2014R1A1A1036677). The authors would like to thank Shankar P. Sastry for providing the source code of the log-barrier interior point method.