# An Efficient Approach for Solving Mesh Optimization Problems Using Newton’s Method

Authors: Jibum Kim
Journal: Mathematical Problems in Engineering
 (2014)

---

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


---

## Body

## 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,F
(
x
), 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 q
i be the quality of ith vertex and N the total number of vertices in the mesh. Then, F
(
x
) is formulated as ∑
i
=
1
N
‍
q
i. We use a condition number quality metric for mesh quality improvement. For trivalent mesh corners in 3D, given by edge vectors e
1, e
2, and e
3, the quality of the vertex corner is given by [11]
(1)
q
i
=
A
·
L
V
,
where A
=
∥
e
1
×
e
2
∥
2
+
∥
e
2
×
e
3
∥
2
+
∥
e
3
×
e
2
∥
2, L
=
∑
j
=
1
3
‍
∥
e
j
∥
2, and V
=
(
e
1
×
e
2
)
·
e
3. The condition number quality metric is scale-invariant and prefers a right angle corner. Here, the objective function, F
(
x
), 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
(2)
q
i
=
(
A
)
(
L
)
V
+
V
2
+
δ
2
,
where δ is a user-defined constant value. If V 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,
(3)
δ
=
1
1
+
c
e
V
/
V
0
,
where V
0 is a reference element volume for scale invariance, which could be chosen as some average element volume around the element being untangled and c is a constant. In our experiments we determined that 10
<
c
<
100 is a good practical choice. Figure 1 shows an example of an adaptive δ function using a sigmoid function.Figure 1
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

Letx
k be the (free) vertex coordinate at the kth iteration. Nonlinear conjugate gradient (NLCG) method determines a step length, λ
k, by using a line search technique along a search direction, p
k. The step length, λ
k, decreases F
(
x
) along p
k. The x
k is updated by computing
(4)
x
k
+
1
=
x
k
+
λ
k
p
k
.
The search direction, p
k, is given by
(5)
p
k
=
-
∇
F
(
x
k
)
+
ϕ
k
p
k
-
1
,
where p
0
=
-
∇
F
(
x
0
), and ϕ
k is a parameter given by
(6)
ϕ
k
=
∇
F
(
x
k
)
T
∇
F
(
x
k
)
∇
F
(
x
k
-
1
)
T
∇
F
(
x
k
-
1
)
,
for the Fletcher and Reeves NLCG method [6]. At each iteration, the variation in NLCG is from the computation of ϕ
k [2].For volume meshes, the line search with theArmijo condition [6] is performed such that the updated vertex position decreases F
(
x
) and does not result in a tangled mesh. For surface meshes, the line search with theArmijo condition is performed with local parametric coordinates as described in [3]. The updated parametric space should satisfy both themesh validity and surface constraints, while decreasing F
(
x
). 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 p
k as
(7)
p
k
=
-
∇
F
(
x
k
)
.
The vertex coordinate is updated by computing
(8)
x
k
+
1
=
x
k
+
λ
k
p
k
,
where λ
k is a step length. Similar to the nonlinear conjugate gradient method, the λ
k is determined byArmijo 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 minimizeF
(
x
) using Newton’s method by finding vertex coordinates, x, such that ∇
F
(
x
) is zero. We use finite difference approximations to compute the gradient (∇
F
(
x
)) and Hessian matrix (H
F
(
x
)) of F
(
x
). In the kth iteration, we compute a Newton update, p
(
x
k
), by solving p
(
x
k
)
=
-
(
H
F
(
x
k
)
)
-
1
∇
F
(
x
k
). The kth vertex position is updated as x
k
+
1

=

x
k
+
α
p
(
x
k
), where α is a step-length. Similar to nonlinear conjugate gradient and steepest descent methods, the parameter α is determined by theArmijo condition [6] dictating that the update should lead to a sufficient decrease in the objective function. In practice, we start with α
=
1 and decrease it by 20% until theArmijo condition is satisfied. The Newton direction, p
k, reliably produces a decrease of the objective function only if the matrix, H
F
(
x
k
), is positive definite.WhenH
F
(
x
k
) is not positive definite, we perform a Hessian modification step such that all eigenvalues of H
F
(
x
k
) are positive [6]. In order to make H
F
(
x
k
) (simply, H) 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 H
+
β
I is a positive definite matrix, where β is a positive constant and I is an identity matrix [17]. The β value increases until H becomes positive definite.Algorithm 1:HM1: Adding a multiple of the identity on the diagonal.
Chooseβ
>
0
While
H be not positive definite do
H
←
H
+
β
I
Increaseβ
end whileThe second Hessian modification method shown in Algorithm2 (HM2) is to modify eigenvalues of H by changing the signs of the negative eigenvalues in H. By simply flipping the sign of the negative eigenvalues, it is ensured that all eigenvalues of H are positive and the updated Newton direction, p
k, is a descent direction.Algorithm 2:HM2: Changing the signs of negative eigenvalues.
Find a spectral decomposition ofH
=
S
D
S
t, where D is a diagonal matrix.
H
←
S
|
D
|
S
tThe 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 E of small norm such that H
+
E becomes positive definite and is decomposed using the Cholesky factorization such that H
+
E
=
L
D
L
T. Here, matrix D 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:
(9)
D
j
j
≥
η
,
|
m
i
j
|
≤
γ
,
where m
i
j
=
l
i
j
D
j
j. Using these bounds, diagonal elements, D
j
j, are computed using the equation
(10)
D
j
j
=
max
⁡
(
|
C
j
j
|
,
(
θ
j
γ
)
2
,
η
)
,
where θ
j
=
max
⁡
j
<
i
≤
n
⁡
|
C
i
j
|. Here, n 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 Algorithm4. Here, γ is a constant, which starts from γ
=
1 and decreases until the updated vertex position, x
k, 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 (s
k). 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 themesh 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, (H
F
(
s
k
)), 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 α
=
1 and decreases until the updated vertex position satisfiesArmijo condition. Here, γ is a constant, which starts from γ
=
1 and decreases until the updated parametric space s
k and corresponding x
k satisfymesh 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 Section5 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.Algorithm 3:HM 3: Modified Cholesky factorization.
C
← diag(H) (C is a diagonal matrix whose diagonal entries are diagonal entries of H)
E
←
0
D
←
0
for
j
=
1
,
2
,
…
,
n
do
θ
j
←
0
D
j
j
=
max
⁡
(
|
C
j
j
|
,
(
θ
j
γ
)
2
,
η
), with θ
j
=
max
j
<
i
≤
n
⁡
|
C
i
j
|
E
j
j
←
D
j
j
-
C
j
j
H
j
j
←
H
j
j
+
E
j
j
end forAlgorithm 4:Volume mesh optimization and untangling using Newton’s method.
whilenot convergeddo
Compute the gradient(
∇
F
(
x
k
)
) and Hessian (
H
←
H
F
(
x
k
)
)
if
H is not positive definitethen
Perform Hessian modification
end if
Newton direction:p
(
x
k
)
=
-
H
-
1
∇
F
(
x
k
)
Update the vertex position:x
k
←
x
k
+
α
p
(
x
k
) such that α satisfies the A
r
m
i
j
o c
o
n
d
i
t
i
o
n
if
x
k satisfies mesh validitythen
Update the vertex position x
k
else
Backtracking Line Search: set x
k
←
x
k
+
γ
p
(
x
k
),
end if
end whileAlgorithm 5:Surface mesh optimization using Newton’s method.
Map the vertex coordinatex
k to a parametric space s
k [3]
whilenot convergeddo
Compute the gradient(
∇
F
(
s
k
)
) and Hessian (
H
←
H
F
(
s
k
)
) of q
(
s
k
)
if
H is not positive definitethen
Perform Hessian modification
end if
Newton direction:p
(
s
k
)
=
-
H
-
1
∇
F
(
s
k
)
Update the vertex position:s
k
←
s
k
+
α
p
(
s
k
) such that α satisfies the A
r
m
i
j
o  c
o
n
d
i
t
i
o
n
if
s
k satisfies both the mesh validity and surface constraintsthen
Update the vertex position s
k and map s
k to x
k
else
if
s
k does not satisfy the mesh validitythen
Backtracking Line Search: set s
k
←
s
k
+
γ
p
(
s
k
),
such that γ satisfies mesh validity
end if
if
s
k does not satisfy the surface constraintsthen
Change the parametric space to the neighborhood parametric space
end if
Update the vertex positions
k and map s
k to x
k
end if
end whileFigure2 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.Table 1
Worst element quality after convergence for three Hessian modification methods on surface meshes (Figure6). A smaller value indicates a better mesh quality.
Mesh (number of elements)
HM 1
HM 2
HM 3
Pig (surface, 3K)
10.46
9.85
4.35
Cow (surface, 3K)
8.74
15.22
13.28
Igea (surface, 40K)
7.46
7.19
4.20
Slope (surface, 3K)
1.70
1.70
1.70Figure 2
Convergence of three different Hessian modification (HM) methods on the Pig mesh (Figure6(a)).


### 6.1. Volume Mesh Optimization and Untangling

We consider three volume meshes, which are shown in Figure3. 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 ofmesh validity is not well-defined [18]. Similar to [18], we perform astar-shape test to verify themesh validity for this polyhedral mesh. Here, thestar-shape test is a quite conservative definition ofmesh validity such that all decompositions of tetrahedra from a polyhedron should have positive volumes. Figure 4 shows one example of astar-shape test for polyhedral meshes. More details on thestar-shape test for polyhedral meshes are described in [18].The element types of the three initial volume meshes are (a) 3D slope: polyhedral, (b) 2-material: polyhedral, and (c) cube: hexahedral.
(a)
3D slope model
(b)
2-material model (tangled)
(c)
Cube model (tangled)Star-shapetest for polyhedral meshes. (a) General polyhedron. (b) Triangulation of face and face center. (c) One tetrahedron in the symmetric decomposition of polyhedron which is used to checkmesh validity. Here, v1 and v2 are two vertices on the face. F and R are face and region centers, respectively.
(a)
(b)
(c)Figure5 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.The element types of the two initial volume meshes. (a) Cut of initial polyhedral mesh showing invalid elements. Here, red elements indicate inverted elements. (b) Cut of optimized polyhedral mesh with no inverted elements. (c) Full polyhedral mesh with no inverted elements.
(a)
2-material model (initial)
(b)
2-material model (optimized)
(c)
2-material model (optimized)The element types of the four surface meshes are (a) Pig: mixed elements, (b) Cow: hexahedral elements, (c) Igea: tetrahedral elements, and (d) Slope: polyhedral elements.
(a)
Pig mesh (surface)
(b)
Cow mesh (surface)
(c)
Igea mesh (surface)
(d)
Slope mesh (surface)Table2 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.Table 2
Time (sec) to optimize various volume meshes.
Mesh (number of elements)
Newton
NLCG
Steepest descent
3D slope (volume, 3K)
442
1,256
1,248
2-material (volume, 2K)
902
1,760
1,756
Cube (volume, 3K)
82
411
409Table 3
Worst element quality computed by the condition number metric for volume meshes.
Mesh (number of elements)
Initial
Optimized (Newton)
Optimized (NLCG)
Optimized (steepest descent)
3D slope (volume, 3K)
1,026.80
2.17
2.36
2.36
2-material (volume 2K)
191.12
1.08
1.29
1.31
Cube (volume, 3K)
1,026.00
1.09
1.17
1.17


### 6.2. Surface Mesh Optimization

Figure6 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.Table 4
Time (sec) to optimize various surface meshes.
Mesh (number of elements)
Newton
NLCG
Steepest descent
Pig (surface, 3K)
61
226
218
Cow (surface, 3K)
90
250
241
Igea (surface, 40K)
297
1,422
1,398
Slope (surface, 3K)
504
1,810
1,802Table 5
Worst element quality computed by the condition number metric for surface meshes. A smaller value indicates a better mesh quality.
Mesh (number of elements)
Initial
Optimized (Newton)
Optimized (NLCG)
Optimized (steepest descent)
Pig (surface, 3K)
45.07
4.35
11.34
11.34
Cow (surface 3K)
18.26
12.68
15.26
16.22
Igea (surface, 40K)
19.08
3.68
3.91
3.76
Slope (surface, 3K)
48.30
1.75
2.06
1.88Table 6
Surface mesh quality statistics of cow mesh (Figure6(b)).
Surface element quality
Initial
Optimized (Newton)
Optimized (NLCG)
Optimized (steepest descent)
1.0–1.5
2,143
2,930
2,922
2,926
1.5–2.0
628
146
144
150
2.0–3.0
216
54
53
51
3.0–4.0
79
9
10
12
4.0–5.0
54
3
3
2
5.0–7.5
24
11
9
10
7.5–10.0
6
4
6
6
10.0–15.0
6
1
0
0
>15.0
2
0
1
1


## 6.1. Volume Mesh Optimization and Untangling

We consider three volume meshes, which are shown in Figure3. 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 ofmesh validity is not well-defined [18]. Similar to [18], we perform astar-shape test to verify themesh validity for this polyhedral mesh. Here, thestar-shape test is a quite conservative definition ofmesh validity such that all decompositions of tetrahedra from a polyhedron should have positive volumes. Figure 4 shows one example of astar-shape test for polyhedral meshes. More details on thestar-shape test for polyhedral meshes are described in [18].The element types of the three initial volume meshes are (a) 3D slope: polyhedral, (b) 2-material: polyhedral, and (c) cube: hexahedral.
(a)
3D slope model
(b)
2-material model (tangled)
(c)
Cube model (tangled)Star-shapetest for polyhedral meshes. (a) General polyhedron. (b) Triangulation of face and face center. (c) One tetrahedron in the symmetric decomposition of polyhedron which is used to checkmesh validity. Here, v1 and v2 are two vertices on the face. F and R are face and region centers, respectively.
(a)
(b)
(c)Figure5 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.The element types of the two initial volume meshes. (a) Cut of initial polyhedral mesh showing invalid elements. Here, red elements indicate inverted elements. (b) Cut of optimized polyhedral mesh with no inverted elements. (c) Full polyhedral mesh with no inverted elements.
(a)
2-material model (initial)
(b)
2-material model (optimized)
(c)
2-material model (optimized)The element types of the four surface meshes are (a) Pig: mixed elements, (b) Cow: hexahedral elements, (c) Igea: tetrahedral elements, and (d) Slope: polyhedral elements.
(a)
Pig mesh (surface)
(b)
Cow mesh (surface)
(c)
Igea mesh (surface)
(d)
Slope mesh (surface)Table2 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.Table 2
Time (sec) to optimize various volume meshes.
Mesh (number of elements)
Newton
NLCG
Steepest descent
3D slope (volume, 3K)
442
1,256
1,248
2-material (volume, 2K)
902
1,760
1,756
Cube (volume, 3K)
82
411
409Table 3
Worst element quality computed by the condition number metric for volume meshes.
Mesh (number of elements)
Initial
Optimized (Newton)
Optimized (NLCG)
Optimized (steepest descent)
3D slope (volume, 3K)
1,026.80
2.17
2.36
2.36
2-material (volume 2K)
191.12
1.08
1.29
1.31
Cube (volume, 3K)
1,026.00
1.09
1.17
1.17


## 6.2. Surface Mesh Optimization

Figure6 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.Table 4
Time (sec) to optimize various surface meshes.
Mesh (number of elements)
Newton
NLCG
Steepest descent
Pig (surface, 3K)
61
226
218
Cow (surface, 3K)
90
250
241
Igea (surface, 40K)
297
1,422
1,398
Slope (surface, 3K)
504
1,810
1,802Table 5
Worst element quality computed by the condition number metric for surface meshes. A smaller value indicates a better mesh quality.
Mesh (number of elements)
Initial
Optimized (Newton)
Optimized (NLCG)
Optimized (steepest descent)
Pig (surface, 3K)
45.07
4.35
11.34
11.34
Cow (surface 3K)
18.26
12.68
15.26
16.22
Igea (surface, 40K)
19.08
3.68
3.91
3.76
Slope (surface, 3K)
48.30
1.75
2.06
1.88Table 6
Surface mesh quality statistics of cow mesh (Figure6(b)).
Surface element quality
Initial
Optimized (Newton)
Optimized (NLCG)
Optimized (steepest descent)
1.0–1.5
2,143
2,930
2,922
2,926
1.5–2.0
628
146
144
150
2.0–3.0
216
54
53
51
3.0–4.0
79
9
10
12
4.0–5.0
54
3
3
2
5.0–7.5
24
11
9
10
7.5–10.0
6
4
6
6
10.0–15.0
6
1
0
0
>15.0
2
0
1
1


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



---

*Source: 273732-2014-10-21.xml*