A Derivative-Free Mesh Optimization Algorithm for Mesh Quality Improvement and Untangling

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.


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][6][7][8][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][6][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.

Problem Formulation
We mathematically formulate mesh quality metrics and objective functions to optimize.We also describe the local mesh optimization method we used.

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 as Let 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,  −1 becomes an identity matrix.In this case,  imr becomes one.
We choose an untangling beta ( ub ) 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,  ub becomes zero.For the above two quality metrics, a smaller value indicates better element quality.

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 max where  is the number of elements.Then, our goal is to minimize (5).The optimization problem is formulated as min max 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.

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).The seven vertices and elements, which are adjacent to the free vertex, V, are optimized at this time.This local mesh optimization is repeated for other vertices on the mesh until the entire mesh is fully optimized.

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.

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, V. Two (virtual) points, V 1 and V 2 , and V. form an initial simplex.Here,  is computed as  = (simplex diameter) * (minimum edge length) , (7) Figure 2: This figure shows an initial simplex (triangle) to begin when a free vertex V is given.A free vertex V and two (virtual) points, V 1 and V 2 , compose an initial simplex.The  is a user-defined parameter defined in (7).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.

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 { 1 ,  2 ,  3 } such that We find the vertex with the worst cost function on the simplex (i.e.,  3 ) 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 as We consider three possibilities for this reflected point, (−1).
(1) 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.
(2) 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.
(3) 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,  1 ,  2 , and  3 , reflected point, (−1), expansion point, (−2), and contraction point, (1/2).We repeatedly perform the downhill simplex algorithm, which is illustrated in Algorithm 1, until all vertices converge to the local optimum.

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.

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

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  = 0.1 * (minimum edge length).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.

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][6][7], is to minimize 11) or to minimize ∑  =1  2  , 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 Figure 5: These figures show the worst element quality and timing results to converge for various initial simplex diameters.A smaller value indicates a better element quality for the inverse mean ratio quality metric.In (a) and (b), the optimization problem does not converge when the initial simplex diameter is five.In (c) and (d), the optimization problem does not converge when the initial simplex diameter is five.In (e) and (f), all inverted elements are eliminated after optimization when the initial simplex diameters are 0.01 and 0.1.Other initial simplex diameters fail to eliminate inverted elements.In (g) and (h), the optimized meshes fail to untangle inverted elements when the initial simplex diameter is two and five.the th iteration and let   be the search direction.Then, the updated vertex position,  +1 , is computed by where   is a step length.The search direction,   , is denoted by where  PR  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 (c) Optimized cylinder mesh using traditional method (zoomed-in) Figure 7: (a) Initial cylinder mesh (zoomed-in).This mesh has a very skinny element (red) with the element quality of 678.9248 enclosed by the red circle.(b) Optimized cylinder mesh (zoomed-in) using the proposed approach.The downhill simplex method is used to improve the worst element quality.The worst quality element has the worst element quality with 3.73.(c) Optimized cylinder mesh (zoomed-in) using the traditional approach.The traditional approach focuses on improving the average element quality using the nonlinear conjugate gradient method.The worst quality element has the worst element quality with 8.23.
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 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)).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.

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 derivativefree 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 1/  and  is greater than 0. The log-barrier method maximizes  with respect to (c) Optimized cylinder mesh using the traditional method (zoomed-in) Figure 8: (a) Initial hydrocephalus mesh (zoomed-in).This mesh has a very skinny inverted element (red) with the element quality of 551.5212 enclosed by the red circle.(b) Optimized cylinder mesh (zoomed-in) using the proposed approach.The downhill simplex method is used to improve the worst element quality.After performing mesh optimization, all inverted elements are eliminated.(c) Optimized cylinder mesh (zoomed-in) using the traditional approach.The traditional approach focuses on improving the average element quality using the nonlinear conjugate gradient method.The optimized mesh does not have any inverted elements.
Cylinder mesh Bar mesh  The inverse mean ratio quality metric is used to measure the mesh quality.A smaller value indicates a better mesh quality.
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.

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

Figure 1 :
Figure 1: This figure shows an example of local mesh optimization.The seven vertices and elements, which are adjacent to the free vertex, V, are optimized at this time.This local mesh optimization is repeated for other vertices on the mesh until the entire mesh is fully optimized.

Figure 3 :
Figure 3: This figure shows one step of the downhill simplex method in 2D.The current simplex has three points,  1 ,  2 , and  3 .The reflected point (−1) is located at the opposite side of the point,  3 , with the worst cost function.Here, (−2) and (1/2) are expansion and contraction points, respectively.
of mesh optimization Number of inverted elements on the mesh (d) Airfoil mesh

Figure 6 :
Figure 6: Note that a smaller value indicates a better element quality for these figures.The inverse mean ratio quality metric is employed to measure the element quality for (a) and (b): (a) worst element quality as a function of number of iterations of mesh optimization on the cylinder mesh.(b) Worst element quality as a function of number of iterations of mesh optimization on the bar mesh.(c) Number of inverted elements on the hydrocephalus mesh as a function of number of iterations of mesh optimization.(d) Number of inverted elements on the airfoil mesh as a function of number of iterations of mesh optimization.
Optimized cylinder mesh using downhill simplex method (zoomed-in) Optimized hydrocephalus using the proposed method (zoomed-in)

Figure 9 :
Figure9: Comparison of the worst element quality on the optimized mesh between the proposed method and the log-barrier method.The inverse mean ratio quality metric is used to measure the mesh quality.A smaller value indicates a better mesh quality.

Table 1 :
The test mesh configurations.

Table 3 :
Connection between the edge length on the mesh and the best initial simplex diameter.

Table 4 :
Summary of the proposed and traditional mesh optimization methods.