Abstract
We present an efficient approach for solving various mesh optimization problems. Our approach is based on Newton’s method, which uses both first-order (gradient) and second-order (Hessian) derivatives of the nonlinear objective function. The volume and surface mesh optimization algorithms are developed such that mesh validity and surface constraints are satisfied. We also propose several Hessian modification methods when the Hessian matrix is not positive definite. We demonstrate our approach by comparing our method with nonlinear conjugate gradient and steepest descent methods in terms of both efficiency and mesh quality.
1. Introduction
Mesh quality improvement and mesh untangling are important topics for partial differential equation- (PDE-) based simulations, because elements of poor quality can be detrimental to accuracy and efficiency of the solution, and an inverted (tangled) element can result in the failure to obtain a PDE solution [1, 2]. Mesh quality improvement (also, mesh untangling) problems are often formulated as nonlinear objective functions [2–5]. Improvement of mesh quality is most often done by minimization of the nonlinear objective function over the whole mesh based on the quality of individual element in the mesh. In order to efficiently minimize the nonlinear objective function, various nonlinear solvers have been developed. These solvers include steepest descent, conjugate gradient, inexact Newton, feasible Newton, quasi Newton, and trust-region methods [6–8]. These solvers can be roughly divided into two groups: gradient-based methods and Hessian-based methods. On one hand, gradient-based methods require less computation and memory than Hessian-based methods, since they do not require Hessian computation. Hessian-based methods on the other hand enjoy the use of more information at the cost of more computation.
The performance of gradient and Hessian-based methods for solving volume mesh optimization problems has been investigated in [9]. Freitag Diachin et al. [10] compared inexact Newton and the block coordinate descent method for solving volume mesh optimization problems. The authors report that the block coordinate descent method is more efficient in achieving a suboptimal solution than the inexact Newton method, because the inexact Newton method incurs high setup costs. The inexact Newton method outperforms the block coordinate method when higher accuracy is desired. The block coordinate descent studied in [10] is basically identical to Newton’s method, which uses the Newton direction as a line search direction. However, it focuses on volume mesh quality improvement problems, which do not include any inverted elements and do not include surface mesh optimization problems. Many researchers have used nonlinear conjugate gradient method (NLCG) for solving mesh optimization problems for various applications [1, 2, 5, 11, 12]. However, it turns out that the NLCG method is slow, prone to getting stuck in local minima, and often fails to converge, especially for surface meshes due to surface mesh constraints.
In this paper, we employ Newton’s method for solving various nonlinear mesh optimization problems. We extend our preliminary results in [13] and present a comparison between Newton’s method and gradient-based methods. The use of Newton’s method for solving nonlinear optimization problem is motivated by the observation that if the nonlinear functional is sufficiently smooth, its optimal points are roots of the derivative function. Thus, a nonlinear optimization problem is transformed into a root finding problem. We choose Newton’s method for solving the mesh optimization problem, because it provides both search direction and a step size using both first-order (gradient) and second-order (Hessian) information. Therefore, for many cases, Newton’s method enjoys a full step length. In addition, close to the solution, its convergence is typically quadratic.
There might exist several potential issues when Newton’s method is used for solving mesh optimization problems. First, computational overhead may increase when we need to compute entries in the Hessian matrix of the objective function [6]. However, if we locally optimize one vertex at a time, Hessian matrix is very small (e.g., 2-by-2 matrix for 2D meshes) and the computational overhead is minimal. A second problem occurs when the Hessian matrix of the objective function is not positive definite. In that case, the Newton direction is not reliable. However, if a Hessian modification method is utilized, the Hessian matrix becomes positive definite and the new search direction becomes reliable. We also develop several Hessian modification methods, which can be easily applied to various mesh optimization problems.
2. Mesh Quality Optimization
Objective function for mesh optimization, , is most often formulated by accumulating a quality metric over all the elements of the mesh and independent variables are coordinates of the movable vertices of the mesh. Let be the quality of th vertex and the total number of vertices in the mesh. Then, is formulated as . We use a condition number quality metric for mesh quality improvement. For trivalent mesh corners in 3D, given by edge vectors , , and , the quality of the vertex corner is given by [11] where , , and . The condition number quality metric is scale-invariant and prefers a right angle corner. Here, the objective function, , is a nonlinear function, since the quality metric is nonlinear and each vertex position affects the quality of all its connected elements.
For initially tangled meshes, we employ Escobar et al.’s [4] modification to the quality metric, which allows us to simultaneously untangle and improve mesh quality. The quality metric for 3D meshes [4] is defined as where is a user-defined constant value. If is positive and is zero, the quality metric is the same as the condition number quality metric. Introducing makes the objective function continuous at zero volume. The choice of for different problem sizes is not well-defined in [4]. The value should be as small as possible (close to zero) in order to remain close to the original condition number quality metric, but bigger values of make the function less steep and ensure more robust gradient computations. In order to satisfy these conflicting requirements, we developed an adaptive and smooth function, which satisfies the above requirements using a sigmoid function, where is a reference element volume for scale invariance, which could be chosen as some average element volume around the element being untangled and is a constant. In our experiments we determined that is a good practical choice. Figure 1 shows an example of an adaptive function using a sigmoid function.
In surface meshes embedded in 3D space, the movement of vertices is subject to the additional constraint that the vertices must remain on the original surface. This is best done by moving the vertices in a parametric space that maps to the surface. This method is superior to the method of allowing the optimization to derive the vertices off the surface and then pulling them back to the surface, since errors can be made in the projection step. The parametric space used can be a global parametric space [14–16] or a local parametric space based on each surface element as detailed in [3]. We use a local parametric space, since a global parametric space requires a huge computational cost. A local parametric space for a triangle element and a quadrilateral element are derived using a barycentric mapping and isoparametric mapping, respectively [3]. The optimization of the objective function with respect to a parametric space introduces additional nonlinearities in the objective function making it harder to optimize.
In any case, optimization of all types of unstructured meshes is challenging and numerical methods often experience problems of slow convergence, getting stuck in local minima far from the optimal vertex positions. Typically, meshes that need to be optimized have thousands to millions of vertices, so in principle, the objective function and the independent vector involve all these coordinates. While it is possible to iteratively compute vertex updates to this large vector, it is quite inefficient to do so for two reasons: (1) the Jacobian of a mesh quality based objective function is quite sparse, since a vertex only affects the quality of its connected elements; (2) if the update of a vertex is to be truncated due to an explicit constraint like mesh validity or surface constraint, then the update of all vertex positions gets truncated.
Therefore, it is much more useful to use a Gauss-Seidel type (local) approach to optimize an objective function by considering the movement of one vertex at a time driven by minimization of the local component of the global objective function. It is important to reemphasize that the local objective function must necessarily include all the terms of the global objective function that the vertex position influences. Otherwise, we will be optimizing some other global function. This local approach is more memory efficient, since the local derivative matrices are small but dense, and it also leads to faster convergence because vertices are not severely constrained by the lack of movement in other far away vertices. For these reasons, we employ local mesh optimization, which optimizes only one free vertex at a time. We visit one free vertex at a time and iteratively move to other vertices. More details on local mesh optimization are described in [3, 8].
3. Nonlinear Conjugate Gradient Method Applied to Mesh Optimization
Let be the (free) vertex coordinate at the th iteration. Nonlinear conjugate gradient (NLCG) method determines a step length, , by using a line search technique along a search direction, . The step length, , decreases along . The is updated by computing The search direction, , is given by where , and is a parameter given by for the Fletcher and Reeves NLCG method [6]. At each iteration, the variation in NLCG is from the computation of [2].
For volume meshes, the line search with the Armijo condition [6] is performed such that the updated vertex position decreases and does not result in a tangled mesh. For surface meshes, the line search with the Armijo condition is performed with local parametric coordinates as described in [3]. The updated parametric space should satisfy both the mesh validity and surface constraints, while decreasing . Here, mesh validity means the updated mesh should not be a tangled mesh. The updated vertex position should stay within parametric bounds to satisfy surface constraints. More details on applying NLCG method to surface meshes are described in [3].
4. Steepest Descent Method Applied to Mesh Optimization
The steepest descent method [6, 9] is a line search technique, which chooses the search direction as The vertex coordinate is updated by computing where is a step length. Similar to the nonlinear conjugate gradient method, the is determined by Armijo condition such that the updated vertex position should sufficiently decrease in the objective function, while preventing a tangled mesh.
5. Newton’s Method Applied to Mesh Optimization
We minimize using Newton’s method by finding vertex coordinates, , such that is zero. We use finite difference approximations to compute the gradient () and Hessian matrix () of . In the th iteration, we compute a Newton update, , by solving . The th vertex position is updated as , where is a step-length. Similar to nonlinear conjugate gradient and steepest descent methods, the parameter is determined by the Armijo condition [6] dictating that the update should lead to a sufficient decrease in the objective function. In practice, we start with and decrease it by 20% until the Armijo condition is satisfied. The Newton direction, , reliably produces a decrease of the objective function only if the matrix, , is positive definite.
When is not positive definite, we perform a Hessian modification step such that all eigenvalues of are positive [6]. In order to make (simply, ) positive definite, we consider three different Hessian modification methods. The first Hessian modification method shown in Algorithm 1 (simply, HM1) is to find a scalar such that is a positive definite matrix, where is a positive constant and is an identity matrix [17]. The value increases until becomes positive definite.
| ||||||||||
The second Hessian modification method shown in Algorithm 2 (HM2) is to modify eigenvalues of by changing the signs of the negative eigenvalues in . By simply flipping the sign of the negative eigenvalues, it is ensured that all eigenvalues of are positive and the updated Newton direction, , is a descent direction.
| ||||
The third Hessian modification method (HM3) is a more sophisticated method that uses a Cholesky factorization [6]. The goal of this Hessian modification method is to find a matrix of small norm such that becomes positive definite and is decomposed using the Cholesky factorization such that . Here, matrix should be sufficiently positive definite but should not be too large [6]. The two positive parameters, and , are chosen such that the following inequality holds: where . Using these bounds, diagonal elements, , are computed using the equation where . Here, values are 2 and 3 for 2D and 3D meshes, respectively. It is known that output matrices computed using HM3 have bounded condition numbers [6].
Our volume mesh optimization algorithm using Newton’s method is summarized in Algorithm 4. Here, is a constant, which starts from and decreases until the updated vertex position, , satisfies the mesh validity. This means the updated vertex position should not tangle the mesh.
As described earlier, surface meshes are optimized by performing optimizations with respect to local parametric coordinates (). We use a condition number quality metric for mesh quality improvement [3]. For surface meshes, the updated vertex coordinates should satisfy both surface constraints and the mesh validity. If the updated vertex coordinates are located beyond the parametric bound, we change the parametric space to the neighborhood parametric space to satisfy the surface constraints. Similar to volume meshes, the Hessian modification step is employed when the Hessian matrix, (), is not positive definite. Algorithm 5 describes the surface mesh optimization algorithm using Newton’s method. Similar to volume meshes, is a constant which starts from and decreases until the updated vertex position satisfies Armijo condition. Here, is a constant, which starts from and decreases until the updated parametric space and corresponding satisfy mesh validity.
6. Numerical Experiments
First, we discuss the Hessian modification methods for solving mesh optimization problems and identify the most efficient one. Next, we compare Newton’s method with the nonlinear conjugate gradient (NLCG) and steepest descent methods in terms of mesh quality and computational cost for both volume and surface meshes. We also discuss whether each method is able to untangle all inverted elements, when initial mesh is tangled.
We compare three different Hessian modification (HM) methods discussed in Section 5 for surface meshes. We focus on surface meshes since, for volume meshes, the difference between the three Hessian modification methods is negligible. The value that is needed for the modified Cholesky factorization method (HM3) is chosen such that it is a sufficiently positive constant. We observe that a small value of close to zero in Algorithm 3 is not desirable for test meshes, since it could result in one of the eigenvalues of the modified Hessian matrix to be zero and thus results in an indefinite Hessian matrix. For our test meshes, we observe that value around one works well.
| ||||||||||||||||||
| ||||||||||||||||||||||||||
| ||||||||||||||||||||||||||||||||||||||||||
Figure 2 shows the worst element quality with respect to the number of Newton iterations for the pig surface mesh (see Figure 6(a)). Here, a smaller value indicates a better mesh quality. We observe that the choice of Hessian modification method considerably affects the mesh quality. We also observe that the modified Cholesky factorization method (HM3) results in the smallest worst element quality compared with the other two Hessian modification methods. The other two methods show a similar worst element quality. Table 1 summarizes the worst element quality of surface meshes (see Figure 6) for three Hessian modification methods when Newton’s method converges to the local optimum. We observe that HM3 results in the smallest worst element quality for three meshes out of our four example meshes. When HM1 and HM2 methods are employed, we observe that the modified Hessian matrix sometimes becomes a nearly singular matrix or struggles with a huge condition number. Also, these two methods sometimes result in a Hessian matrix where the second-order information is not well preserved [6]. The modified Hessian matrix using the HM3 method has a bounded condition number, while preserving the second-order information. Therefore, the optimized mesh using the HM3 method results in better output meshes with good mesh qualities than other two Hessian modification methods. For these reasons, we employ HM3 as our Hessian modification method for the rest of the paper.
6.1. Volume Mesh Optimization and Untangling
We consider three volume meshes, which are shown in Figure 3. The first two meshes are polyhedral meshes, which are generated using the polyhedral mesh generation algorithm in [18]. For the cube mesh, the initial mesh is randomly perturbed such that 60% of the elements are inverted. The 3D slope mesh (Figure 3(a)) does not include any inverted elements and the 2-material model (Figure 3(b)) and cube meshes (Figure 3(c)) are tangled meshes which include inverted elements. For the 2-material model mesh, two different colors indicate two different materials. For the polyhedral mesh such as the 2-material model mesh, the definition of mesh validity is not well-defined [18]. Similar to [18], we perform a star-shape test to verify the mesh validity for this polyhedral mesh. Here, the star-shape test is a quite conservative definition of mesh validity such that all decompositions of tetrahedra from a polyhedron should have positive volumes. Figure 4 shows one example of a star-shape test for polyhedral meshes. More details on the star-shape test for polyhedral meshes are described in [18].
(a) 3D slope model
(b) 2-material model (tangled)
(c) Cube model (tangled)
(a)
(b)
(c)
Figure 5 shows both initial and optimized meshes using Newton’s method for a 2-material model. The initial polyhedral mesh includes several inverted elements (Figure 5(a)). Newton’s method is able to simultaneously untangle and improve element qualities as shown in Figures 5(b) and 5(c). For the 2-material model, both NLCG and steepest descent methods are also able to untangle inverted elements, but they take more number of iterations to untangle all inverted elements compared with Newton’s method. For a cube mesh, we observe that two elements fail to be valid and remain tangled after optimization with NLCG and steepest descent methods. However, Newton’s method is able to untangle all inverted elements. The objective function for simultaneous untangling and smoothing is quite smooth; it converges to the local optimum even if the initial guess (initial mesh) is far from optimal.
(a) 2-material model (initial)
(b) 2-material model (optimized)
(c) 2-material model (optimized)
(a) Pig mesh (surface)
(b) Cow mesh (surface)
(c) Igea mesh (surface)
(d) Slope mesh (surface)
Table 2 shows timing results. We observe that, for our three example meshes, Newton’s method converges between about 50% and 85% faster than the NLCG method. NLCG and steepest descent methods show similar convergence time. Newton’s method is fast since, for most cases, it enjoys the full step length. Also, its convergence is quadratic around the optimal point. Table 3 shows worst element qualities computed using the condition number metric. A smaller value indicates a better mesh quality. The worst element quality using NLCG is up to 20% worse than the one using Newton’s method.
6.2. Surface Mesh Optimization
Figure 6 shows four surface meshes to optimize. Table 4 shows timing results in seconds for various surface meshes until they are optimized. Similar to volume meshes, we observe that Newton’s method converges significantly faster than the NLCG and steepest descent methods. Our experimental results show that time to convergence for Newton’s method is up to 4.8 times faster than for the NLCG method. Table 5 shows worst element qualities for various surface meshes. A smaller value indicates a better mesh quality. We observe that the mesh quality difference between NLCG and Newton’s method is huge for surface meshes. For all these examples, Newton’s method outperforms both the NLCG and steepest descent methods. For the cow (surface) mesh, the output mesh using Newton’s method yields the optimized mesh which has 61% better quality than the mesh that was optimized using the NLCG method. This is because Newton’s method is able to approach the local optimum while NLCG and steepest descent fail to do that when surface constraints exist. Table 6 shows surface mesh quality statistics of the cow mesh (see Figure 6(b)). We observe that Newton’s method outperforms both the NLCG and steepest descent methods in terms of both average element quality and worst element quality. We observe similar statistics for other surface meshes in Figure 6.
7. Conclusion
We have proposed an efficient method for solving various mesh optimization problems using Newton’s method. We observed that Newton’s method significantly outperforms both the nonlinear conjugate gradient (NLCG) and steepest descent methods in terms of both the speed and the mesh quality. Obviously, the major benefit of using Newton’s method is fast convergence. We observe that, for our examples, Newton’s method converges up to 4.8 times faster than the NLCG method for solving mesh optimization problems. We highlight that, even for initially tangled meshes, Newton’s method exhibits robust performance to untangle inverted elements and fast convergence. The benefit of using Newton’s method is more critical for surface meshes due to surface constraints. For many surface meshes with surface constraints, Newton’s method is able to approach the local optimum while NLCG and steepest descent methods fail to do so. For our surface mesh examples, the optimized mesh with Newton’s method shows up to 61% better mesh quality when compared with the optimized mesh with the NLCG method in terms of the worst element quality. We used relatively simple finite difference approximations to compute the gradient and Hessian. We plan to apply more accurate finite difference approximations in the future.
Conflict of Interests
The author declares that there is no conflict of interests regarding the publication of this paper.
Acknowledgments
The author would like to thank Rao V. Garimella and Markus Berndt at the Los Alamos National Laboratory for helpful discussions. This research was supported by Basic Science Research Program through the National Research Foundation of Korea (NRF) funded by the Ministry of Science, ICT & Future Planning (NRF-2014R1A1A1036677).