High-order Mesh Generation Based on Optimal Affine Combinations of Nodal Positions

In this paper, we propose a new optimization-based method for high-order mesh generation. Our proposed method is based on affine combinations of nodal positions. For each interior node in the high-order straight-sided mesh, a set of weights that relates the node to its neighbors is calculated. After the high-order boundary nodes are curved to conform with the true boundary, the interior nodes are moved based on the generated weights and the new boundary positions. We present results from applying our method to several aerospace engineering geometries in 2D and 3D which demonstrate the capabilities of our method.


Introduction
The advantage of high-order numerical methods for solving partial differential equations is their higher degree of accuracy compared to low-order numerical methods. A major hurdle in the usage of these methods in the presence of complex geometries is the absence of robust high-order mesh generation methods [23]. In other words, these methods need a high-order mesh that accurately captures the features of the geometry to achieve their full potential [1,10].
The typical approach for generating high-order meshes is to transform a coarse linear mesh [2][3][4][5][6]9,11,12,14,16,17,[19][20][21][22]24]. At a high-level, these transformations usually consist of the following three steps: (1) the low-order mesh is enriched with additional nodes; (2) the new nodes that lie along the boundary of the mesh are moved to the true boundary; (3) the interior nodes are moved based on the boundary deformation. The main challenge of these methods arises from step (2). In particular, the curving of the elements along the boundary can result in invalid mesh elements. With that in mind, these high-order mesh generation methods use different approaches in step (3) in an effort to obtain a valid high-order mesh. Methods for transforming the linear mesh usually fall into two groups. The first group of methods transform the mesh based on the solution to a partial differential equation [3,11,14,24]. The second group of methods are based on optimization of an objective function [2,4,5,9,[16][17][18][19][20][21].
In this paper, we describe an optimization-based approach for generating highorder meshes based on affine combinations of nodal positions. The remainder of this paper is organized as follows. In Section 2, we present our new method for high-order mesh generation. In Section 3, we demonstrate the performance of our proposed method on several aerospace engineering geometries. Finally, in Section 4, we offer some concluding remarks and possible directions for our future work.

High-Order Mesh Generation Based on Affine Combinations of Nodal Positions
In this section, we present our optimization-based method for high-order mesh generation. Our proposed method uses affine combinations of nodal positions to determine the movement of the interior nodes after deforming the boundary. Our method consists of three steps. First, for each interior node in the high-order straightsided mesh, an optimization problem is solved to calculate a set of weights that relates the interior node to its neighbors. Second, the boundary nodes are moved to the true boundary. Third, the new positions of the interior nodes are calculated by solving a linear system of equations using the weights and the new boundary positions. In spirit, this method is similar to the weight-based method that we proposed in [19] with two major differences. The first difference is that we propose an affine combination of nodal positions in this work, as opposed to a convex combination. This change allows us to remove the inequality constraint and log-barrier term, leaving only the equality constraints. We also propose an alternative objective function that when combined with the equality constraints allows us to directly solve the optimization problem via a QR factorization. This change results in simplified computational complexity and faster execution time.
To frame our discussion of the method, we introduce the following notation for the 2D formulation of the problem; the 3D formulation is similar. Let the x-and y-coordinates of the i th interior node be represented as (x i , y i ). In addition, define the x-and y-coordinates of the vertices adjacent to node i as {(x j , y j ) : j ∈ N i }, where N i is the set of neighbors of node i. For each interior node i, this information can be represented as the following linear system, where w i j are the weights: There are several potential choices for the local neighboring set based on use of the low-order nodes, high-order nodes, or both. We include only the high-order nodes as neighbors, as only the high-order boundary nodes move during the boundary deformations. Using either the low-order nodes or both the low-and high-order nodes would dampen the effect the boundary deformation has on the interior nodes, which might lead to tangling near the boundary. Including additional nodes as neighbors would also result in a less sparse matrix when solving (7). Another important consideration is that while the weight calculation is based on only the local neighbors, the position of an interior node is indirectly affected by the deformation of all the interior nodes through the solution of (7).
Adding the additional constraint that the weights sum to one results in the following linear system Aw = b for finding an affine combination of the x-and y-coordinates of the vertices adjacent to node i: Based on the set of neighbors, this linear system will be underdetermined (i.e., A = m × n with m < n) in general. If we assume that A has full rank, we can find one particular solution to our problem by requiring that w has the smallest norm of any solution. This results in the following optimization problem: From the Karush-Kuhn-Tucker (KKT) theory [13], we know that the following conditions must hold for a solution (w * , λ * ) to our problem to be optimal: The Lagrangian of our problem is given by: where λ are the Lagrange multipliers.
Using (3)- (5), we can find the following solution pair (w * , λ * ) as follows: Although we have verified that (w * , λ * ) is a stationary point, we cannot yet claim that it is a minimum. To do so, we must investigate ∇ 2 w L(w * , λ * ): From the second-order sufficient conditions, if w * satisfies (3)-(5) and the following condition is satisfied: where C(w * , λ * ) = {z | ∇ w c(w * ) T z = 0} is the critical cone and c(w) = Aw − b, then our solution is a minimum. Since ∇ 2 w L(w * , λ * ) is symmetric positive definite, the inequality in (6) is satisfied for any choice of z. Thus we can conclude that our solution w * is a minimum of (1)- (2). Now that we have established that w * is our solution, we will discuss calculating it via a reduced QR factorization. Suppose that A T = QR, where Q n×m ,R m×m is upper triangular, and Q T Q = I m×m . Substituting in the QR factorization of A T into w * , we get the following: Rearranging this into linear system form, we have: If we let t = Q T w * , then R T t = b and w * = Qt. Thus calculating w * involves a QR decomposition of A T , solving the lower triangular system R T t = b by forward substitution, and calculating the matrix-vector product Qt.
After calculating the weights, a boundary deformation is applied. The final step is to solve for the new locations of the interior nodes [x I ,ŷ I ] by solving the following global linear system: wherex B andŷ B are the new x-and y-coordinates for the boundary nodes, and A I and A B contain the weights for the interior nodes and boundary nodes, respectively. In this global linear system, each row of the weight matrix corresponds to an interior node with nonzero entries for the node's neighbors and zero entries for the remainder of the row. The resulting global weight matrix is very sparse with irregular structure.
In an effort to shift the nonzero entries closer to the diagonal, we apply the sparse reverse Cuthill-McKee ordering provided in Matlab. In Fig. 1, we show the matrix sparsity plots for the natural node ordering and the updated node ordering for the first two examples in Sec. 3. After applying the matrix reordering, the linear system is solved using a sparse LU factorization.

Numerical experiments
In this section, we demonstrate the results from applying our method to generate several high-order meshes. We use Gmsh [7,8,15,21] to generate the initial straightsided high-order meshes. Our method then uses this mesh to calculate the weights (step 1). Next, we curve the boundary nodes (step 2) using Gmsh. The positions of the interior nodes in the resulting curved high-order mesh are then updated (step 3) by our method. For each example, we show the mesh which results from our method (with high-order nodes visible), the mesh element distortion for the curved boundary mesh generated using Gmsh, and the distortion for the mesh resulting from our method. When reporting the mesh distortion, we list the minimum distortion, maximum distortion, average distortion computed over all elements (referred to as Avg1 in figures), and average distortion computed over curved elements (referred to as Avg2 in figures). The ideal distortion value is 1, indicating that the element is straight. We also list the execution times needed for steps 1 and 3 of our method (excluding I/O) in Tab. 1. The code was run using Matlab R2018a, and the wall-clock execution times were measured on a machine with 16GB of RAM and a Ryzen 7 1700 CPU. All mesh visualizations and distortion evaluations were done using Gmsh. Our first example is a third-order mesh of a square region around a NACA0012 airfoil. In Fig.  2(a-b), we show the mesh resulting from our method and a table of the mesh quality values as measured by the distortion metric. In this example, our method increased the minimum distortion from 0.744 to 0.799, while causing only minor changes in the average distortion.
In our second example, we extrude the NACA0012 airfoil and create a third-order mesh of the resulting region. In Fig. 3(a-b), we show the the mesh resulting from Our third and final example is a second-order mesh of an Airbus A319 aircraft. Unlike our previous examples, this geometry resulted in tangled elements after curving the boundary. Although our method still increased the minimum quality, it was not able to untangle the mesh. To address this, we applied the high-order regularization scheme available in Gmsh as a post-processing step. Aside from changing the target Jacobian range to 0.3 to 2 and fixing the boundary nodes, all other parameters were left at their default values. The untangling for the original mesh took 14.14 seconds, while untangling the mesh resulting from our method required only 1.64 seconds. In Fig. 4(a-c), we show the surfaces of the mesh resulting from our method, a view of a cut through the interior volume mesh, and a table of the mesh quality values. In Fig. 4(c), we list the distortion for the mesh after curving the boundary, the distortion for the mesh resulting from our method, and the distortions of both meshes after applying the regularization scheme in Gmsh. Aside from the third test case, all of these examples were relatively straightforward. In each case, our method increased the minimum distortion when compared to only curving the nodes along the boundary. While additional testing is necessary to confirm this, our results for the third example seem to indicate that our method could be used to reduce the severity of the mesh tangling and thus simplify the work for an untangling method during post-processing.

Concluding remarks and future work
We have presented a new optimization-based method for generating high-order meshes. Our examples have shown that the proposed method based on affine combinations of nodal positions tends to improve the quality of the most distorted elements, while causing minor changes to the least distorted elements. While our approach is optimization-based, we have demonstrated that the optimization problem can be solved directly using a QR factorization as opposed to the typical iterative optimization approach. This change results in lessened computational complexity and reduced execution time.
As part of our future work, we will consider other definitions for the set of neighbors of an interior node. We will also investigate other aspects of the linear system including other node reordering schemes and solvers. Finally, we will apply the untangling method that we proposed in [20] after extending it to 3D.

Acknowledgments
The work of the first author was funded in part by the Madison and Lila Self Graduate Fellowship and NSF CCF grant 1717894. The work of the second author was supported in part by NSF grants CCF 1717894 and OAC 1808553. The authors wish to thank the anonymous referee for his or her careful reading of the paper and for the helpful suggestions which strengthened it.