Simultaneous Smoothing and Untangling of 2D Meshes Based on Explicit Element Geometric Transformation and Element Stitching

: Mesh quality can affect both the accuracy and efﬁciency of numerical solutions. This paper ﬁrst proposes a geometry-based smoothing and untangling method for 2D meshes based on explicit element geometric transformation and element stitching. A new explicit element geometric transformation (EEGT) operation for polygonal elements is ﬁrstly presented. The transformation, if applied iteratively to an arbitrary polygon (even inverted), will improve its regularity and quality. Then a well-designed element stitching scheme is introduced, which is achieved by carefully choosing appropriate element weights to average the temporary nodes obtained by the above individual element transformation. Based on the explicit element geometric transformation and element stitching, a new mesh smoothing and untangling approach for 2D meshes is proposed. The proper choice of averaging weights for element stitching ensures that the elements can be transitioned smoothly and uniformly throughout the calculation domain. Numerical results show that the proposed method is able to produce high-quality meshes with no inverted elements for highly tangled meshes. Besides, the inherent regularity and ﬁne-grained parallelism make it suitable for implementation on Graphic Processor Unit (GPU).


Introduction
Mesh quality has great effects on both the accuracy and efficiency of numerical solutions based on the finite element method and the finite volume method [1][2][3]. Therefore, meshes with high quality are required for various finite element and finite volume based numerical simulations. For a given geometric model, meshes can be generated using various mesh generation methods [4,5]. Among them Delaunay-based methods are quite popular, since they are well suited for automation, of reasonable complexity, and produce meshes with high quality.
For problems with smooth solutions, a uniform mesh is usually preferred over a nonuniform one because it readily lends itself to efficient solution. Even when adaptivity is required, an isotropic mesh can be preferred for the same reason if it can resolve the solution without using an undue number of mesh points. Anisotropic meshes are favored when there is a need for better alignment of the mesh with certain solution directions, such as those arising due to boundary or interior layers and sharp interfaces. The adaptive moving mesh methods are usually used in these situations, where point locations are dynamically relocated in time while adapting to the evolving structures in the solution [6]. Das et al. [7][8][9][10][11][12][13] solved a series of time-dependent partial differential equations (PDEs) using the adaptive moving mesh methods. Their works have rigorous mathematical derivations and elegant proofs. The comparison results with the methods based on priori meshes show that the adaptive solution on posteriori generated meshes will converge to the exact solution with higher accuracy and reliability. The adaptive moving mesh methods can effectively improve the accuracy of numerical simulations, and the mesh optimization methods studied in this paper try to ensure the feasibility and improve the accuracy of numerical simulations.
However, tangled meshes can occur during mesh generation, mesh optimization, large-scale deformation, moving domains, and mesh morphing. Generally speaking, such meshes cannot be directly used for numerical simulations. Inverted elements should be eliminated before carrying out numerical simulations, since they usually produce invalid solutions [14]. There are many methods for eliminating inverted elements, including remeshing [15], local mesh modification [16], topological transformations [17] and mesh smoothing. Remeshing may not be automatic, especially for models with complex geometries. Local mesh modification can be very effective. However, for a tangled mesh that contains a large number of inverted elements, remeshing the whole computational domain may be a better choice. Topological transformations improve the mesh quality by inserting or removing a node or modifying the node-element connection relationship. Edge flip, edge collapse, and node insertion are three popular topological optimization methods to improve the mesh quality. While remeshing, local mesh modification and topological transformations methods are valuable in mesh untangling, we focus on mesh smoothing or optimization via node-movement in this paper.
Mesh smoothing improves the mesh quality by simply relocating or adjusting node positions while preserving the mesh topology. There exists many node movement methods and the majority of them do not specifically address the mesh untangling problem. For example, Laplace smoothing [18] is widely used to improve mesh quality, where node positions are iteratively updated by the arithmetic mean of neighboring node positions. However, it can lead to a tangled mesh from one that initially is untangled. Therefore, a set of methods known as barrier methods are designed to prevent the initial untangled meshes from becoming tangled. However, these methods can not be used when the initial mesh is tangled.
To address the above issue, there are several methods specialized to untangle an initial tangled mesh, and they generally do not consider the element shape. Therefore, a barrier smoothing technique is needed to obtain a valid and high-quality mesh. Freitag and Plassman tried to untangle the input mesh by maximizing the minimum area/volume of elements within each local patch of the mesh [19]. The system is guaranteed to converge, but not necessarily to an untangled mesh. Knupp used a local non-barrier metric that penalized inverted elements [20] to untangle the mesh. Global optimization of the objective function is recommended when untangling larger and more difficult meshes, as all nodes are updated simultaneously. However, the global optimization strategy is very time-consuming in untangling large-scale meshes. As with the other methods, there is no guarantee of obtaining a valid mesh when using this method. Franks and Knupp proposed two new metrics for simultaneous mesh untangling and smoothing using the target matrix paradigm [21]. This method can effectively eliminate inverted elements and do not require that one specify a parameter as in Reference [20]. The above methods [20,21] are adopted for the well-known mesh quality improvement toolkit (Mesquite) [22].
There are also methods developed for simultaneously eliminating inverted elements and improving the mesh quality. Kim et al. [23] proposed a multiobjective optimization method that combined two or more objective functions into a single objective function, and solved it using a nonlinear conjugate gradient method. Escobar et al. [24,25] proposed a method incorporating a modified quality metric that ensured higher penalty for inverted elements. This method is effective but also does not guarantee an untangled output mesh. Gargallo-Peiró et al. [26] proposed a method optimizing triangular and quadrilateral meshes on parameterized surfaces. The method expresses the quality metrics for planar elements in terms of the parametric coordinates of the nodes, and then relocates the nodes on the parametric space to improve the mesh quality and eliminate the inverted elements in the physical space.
Optimization-based algorithms usually untangle and smooth the input meshes at the expense of higher computational and implementational effort. As a way out of the dilemma between quality improvement effect and computational effort, algorithms that do not need to solve numerical optimization problems have been proposed. Kim et al. [27] proposed a novel approach for simultaneous mesh untangling and smoothing using a Pointer network. This method predicts the approximate solutions for free nodes using a pre-trained network, and there is no need to solve complex numerical optimization problems. A geometry-based algorithm for smoothing 2D/3D meshes is proposed by Vartziotis et al. [28,29], which can produce high quality meshes comparable to those obtained by the optimization-based methods within shorter runtimes. Although their method is designed for mesh smoothing, the method still shows its potential for mesh untangling.
The purpose of this paper is to develop a geometry-based mesh untangling and smoothing approach for 2D meshes. A novel explicit element geometric transformation (EEGT) operation is firstly introduced and analyzed, which can transform an arbitrary initial polygonal element (even inverted) into its regular counterpart if applied iteratively. This is the main driving force behind the proposed approach. In the proposed mesh untangling and smoothing approach, all mesh elements are transformed individually and simultaneously in the first step, then new interior node positions are obtained as weighted means of the scaled transformed element nodes. The proper choice of averaging weights in element stitching ensures that the elements can be transitioned smoothly and uniformly throughout the calculation domain, and the explicit element geometric transformation has strong capability to eliminate inverted elements and improve the mesh quality. To the best of our knowledge, the proposed approach is the first geometry-based framework specially designed for simultaneous mesh untangling and smoothing. Besides, it is also very suitable for parallel implementation, especially on GPU.
In the rest of this paper, a new explicit element geometric transformation operation for arbitrary polygonal elements is introduced in Section 2. Basic properties of the transformation are analyzed and a proof for the convergence of the proposed transformation operation is also given. The proof is based on the use of circulant matrices and the decomposition of a polygon into its eigenpolygons. Then the mesh untangling and smoothing approach for 2D meshes based on explicit element geometric transformation and element stitching is presented in Section 3. Typical numerical examples are tested and comparisons with other methods are given in Section 4. Finally, discussions and limitations are presented in Section 5.

Explicit Element Geometric Transformation
In 2D numerical simulations, computational meshes usually consists of triangular and quadrilateral elements. The proposed mesh untangling and smoothing approach is mainly based on a new regularizing element geometric transformation, which can transform an arbitrary initial polygonal element (even inverted) into its regular counterpart if applied iteratively. The transformation of neighbor elements is stitched together to improve the overall mesh quality. In this section, such a transformation for arbitrary polygonal elements will be presented and analyzed.

Basic Geometric Transformation
The element geometric transformation for smoothing triangular and quadrilateral meshes has been introduced by Sun et al. in Reference [30]. It is based on transforming elements by a specially designed two-step stretching and shrinking operation (SSO). The transformation, if applied iteratively to a triangular or quadrilateral element, will improve its regularity and quality. The purpose of stretching operation is to regularize the element, while the shrinking operation is performed to preserve the element size. The stretching operation proposed by Sun et al. [30] is shown in Figure 1 for a triangular and a quadrilateral element.  The element geometric transformation operation proposed by Sun et al. [30] is mainly used for mesh smoothing, and the regularizing effect for arbitrary triangular and quadrilateral elements remains to be proved in their paper. In this subsection, the element geometric transformation operation will be generalized for polygonal elements with an arbitrary number of nodes, which may even be inverted.
For a given polygonal element with counterclockwise oriented nodes p k ∈ R 2 , k ∈ {0, 1, ..., n − 1}, the stretching operation was chosen as a natural extension of that for triangular and quadrilateral elements in Reference [30], that is, just moving node p k along the outward directed normal of the edge p (k+1)modn p (k+n−1)modn . However, numerical tests show that such a stretching operation derived from triangular and quadrilateral elements is not applicable to polygonal elements with more than five nodes. This is illustrated in Figure 2 for an inverted hexagonal element, which is still inverted after applying the transformation in Reference [30] for one hundred times. Inverted elements usually cannot satisfy actual computational requirements, so the simple stretching operation mentioned above needs to be modified and redefined in order to regularize inverted polygonal elements. In the simple stretching operation, new position of the node p k depends only on its two adjacent nodes p (k+1)modn and p (k+n−1)modn , so the simple stretching operation is a kind of local transformation.
The modified stretching operation is conducted by moving node p k along the direction obtained by combining the outward directed normal of edges p (k+r)modn p (k+n−r)modn , r ∈ {1, 2, ..., ceil( n−1 4.0 )}, where the function ceil(x) returns the smallest integer no less than x. The outward directed normal denoted by n kr of the edge p (k+r)modn p (k+n−r)modn is calculated in the form of components by the formula: It is obvious that the vector n kr is not a unit vector, and its magnitude is the length of the edge p (k+r)modn p (k+n−r)modn . Then the nodes p k of transformed element can be obtained by moving the node p k along the direction defined as: It can be seen that the position of each node in the modified stretching operation is determined by at least half of the other nodes in the element, which guarantees that the global effect of the modified stretching operation. With a user-specified scaling factor λ, the new node positions of the transformed element can be calculated by: The above formula shows that the moving distance of node p k is proportional to the magnitude of its stretching vector n k . Therefore, the proposed stretching operation is adaptive to element size. The size of a transformed element after the stretching operation is usually larger than that of the initial element. This can be avoided by scaling the element after stretching operation with respect to the original centroid. Such an operation ensures that the transformed element with more regular shape is close to the initial element and has a small change in element size.
In general, there are two schemes for element scaling, the mean edge length or mean area preserving scaling. The mean edge length preserving scaling scheme is used in this paper for the sake of simplicity. With this scheme, the scaled transformed element with new nodes p s k is given by: where c and c represent the initial and transformed element centroid, respectively. The scaling factor k is defined as the mean edge length of the initial element divided by that of the enlarged transformed element. The stretching and shrinking operation together make up a complete element geometric transformation process. The final node positions of the scaled transformed element can be obtained easily and explicitly with Equations (1)-(4). The geometric transformation process of the same inverted hexagonal element as in Figure 2 using the modified transformation operation is shown in Figure 3. It can be seen that the element becomes a nearly regular hexagonal element after applying the modified transformation for one hundred times, which intuitively illustrates the regularizing capacity of the modified transformation operation. In fact, the modified element geometric transformation operation can successfully transform an arbitrary polygonal element into its regular counterpart after enough iterations, which will be proved in the following subsection.

Analysis of the Transformation
In this subsection, the coordinate of node p k is represented by p k = x + iy. The reason for doing this is to take advantage of the properties of circulant determinants. Using a n × n matrix A with zero-based indexes i, j ∈ {0, 1, ..., n − 1} and a n × 1 vector p l = (p l 0 , p l 1 , ..., p l n−1 ) ∈ C n representing the polygonal element with n nodes after l cycles of transformation, the modified stretching operation governed by Equation (3) can be written in the form of p l+1 = Ap l , where A is a circulant Hermitian matrix given by: a n−2 · · · . . . . . . a n−1 a n−1 a n−2 · · · a 1 a 0 and the component can be written as: In the case of triangular elements this yields: and in the case of quadrilateral elements, respectively. With this, the basic properties of the transformation presented in the previous subsection can be analyzed for arbitrary polygonal elements. Proof of Theorem 1.
Since A is a circulant Hermitian matrix, the associated eigenvalues are real-valued. Due to the circulant structure, the eigenvalues µ k , k ∈ {0, 1, ..., n − 1} of A can be easily obtained by multiplying the discrete Fourier matrix R defined by R u,v = ω uv = exp( 2πiuv n ) with the transposed first row of A [31], where i is the imaginary unit. Hence, for arbitrary n ≥ 3 the corresponding eigenvalues of A are given by: µ j = a 0 + a n−1 ω j + a n−2 ω 2 j + · · · + a 1 ω n−1 j , j = 0, 1, ..., n − 1.
Using the fact that ω n−m n ) =ω m j and a n−m =ā m , the eigenvalues can be written as: In particular, the above formula shows that µ 0 = 1 does not depend on the choice of scaling factor λ. Moreover, it holds that |µ 1 | > |µ j | is true for each value of j in the set j ∈ {2, 3, ..., n − 1}.
Theorem 2. If the transformation parameter λ > 0, the eigenvalue µ 1 has the largest absolute value among all the eigenvalues of matrix A. Figure 4 by plotting the absolute value of function f (n, j) for different pair (n, j). Thus, it holds that |µ 1 | ≥ |µ j | is true for each value of j in the set j ∈ {2, 3, ..., n − 1}.
Inspired by the proof in Reference [28], the diagonalization of A will be used to explore the essence of the proposed transformation. Thus, the matrix composed of the associated eigenvectors in the columns, given by T = 1/ √ nR, R = (ν 0 , ν 1 , ..., ν n−1 ) is used, where the eigenvector ν j of matrix A is given by ν j = (1, ω j , ω 2 j , ..., ω n−1 j ) T [31]. It holds that A = TDT * , where D is the diagonal matrix composed of the eigenvalues of matrix A and T * is the conjugate transpose of matrix T. With the decomposition of the n × n identity matrix I in the form of I = n−1 ∑ k=0 I k , where I k denotes the n × n matrix with the only nonzero entry (I k ) k,k =1, the transformation matrix after l cycles of the modified stretching operation can be written as: Therefore, the corresponding polygonal element after l cycles of the modified stretching operation can be decomposed in the following manner: The so-called eigenpolygons e k do depend only on the initial polygon element p 0 . Examples of such decompositions are depicted in Figure 5 for a triangular, quadrilateral, pentagonal and hexagonal element, where the element nodes are marked with different colors in order to denote the element orientation.    The first two eigenpolygons and their corresponding eigenvalues play an important role in determining the characters of element geometric transformation. Due to µ 0 = 1, the first eigenpolygon is given by: which means all nodes of the polygonal element are placed at the centroid of the initial polygonal element p 0 during the transformation. The second eigenpolygon represents the counterclockwise oriented regular polygonal element, due to the expression: implies (e 1 ) (k+1)modn = ω 1 (e 1 ) k , k ∈ {0, 1, ..., n − 1}, which says that node p (k+1)modn of the eigenpolygon e 1 can be obtained by rotating the preceding node p k 360 n degrees around the origin of coordinates. It is not difficult to understand that the remaining eigenpolygons can be generated by the expression (e m ) (k+1)modn = ω m (e m ) k , k ∈ {0, 1, ..., n − 1}.
The properties above lead to the following crucial result building the foundation for the simultaneous mesh untangling and smoothing approach to be introduced in the next section. Proof of Theorem 3. The size of a transformed element can differ significantly from that of the initial element after the stretching operation. This could be easily overcome by scaling the transformed element after each transformation step with respect to the centroid. From a theoretical point of view, using the inverse of µ 1 as scaling factor is the best choice. Of course, the mean edge length preserving scaling scheme is used in practice. The sequence of scaled transformed polygonal elements using the inverse of µ 1 as scaling factor: converges to p ∞ s = lim l→∞ p l s = e 0 + e 1 , which represents the counterclockwise oriented regular polygonal element centered at the centroid of the initial polygonal element.
The speed of convergence for a single polygonal element depends on the ratio ρ = max 2≤k≤n−1 {|µ k /µ 1 |} < 1. The eigenvalues depend only on scaling factor λ, so does the ratio ρ.
For triangular elements, it holds that µ 1 = 1 + √ 3λ, µ 2 = 1 − √ 3λ, and it is obvious that ρ achieves a minimum of 0 when λ is equal to √ 3/3. As for quadrilateral elements, it holds that µ 1 = 1 + 2λ, µ 2 = 1, µ 3 = 1 − 2λ, and it is easy to obtain that ρ achieves a minimum of 1/3 when λ is equal to 1. In a word, the optimal scaling factor λ for a single polygonal element can be found through theoretical analysis. Nevertheless, numerical tests are usually needed to find the appropriate scaling factor in practical mesh optimization problems, which may involve elements of different types and sizes.

Smoothing and Untangling Approach
The proposed mesh smoothing and untangling approach consists of the following two key steps. First, each element of the input mesh is transformed individually and simultaneously towards its corresponding ideal element by using the transformation introduced in Section 2. In this step, each element is transformed separately without caring about the connection relationship with its neighbors. Therefore, a combination step is needed to stitch the disconnected elements. Here, the element stitching operation is performed by updating the position of each free node by a well-designed weighted means of its corresponding temporary node positions in the scaled transformed elements. The element stitching operation is actually using an explicit formula to approximate the partial derivative of the objective function with respect to the current node position to obtain the displacement of the node.
As the proper choice of a quality metric is important for element stitching operation, the element quality metric used in this paper is first introduced.

Shape Metrics
There are many different shape metrics to measure the quality of a polygonal element. These shape metrics measure the difference between a given polygonal element and a regular polygonal element from different aspects. The mean ratio (MR) metric [32] is one of the commonly used metric in mesh optimization. For a valid polygonal element, the quality number lies within the interval (0,1]. The quality number is equal to 1 only if the element is a regular polygon, and it is defined to be 0 if the element is inverted or degenerated. Overall mesh quality is assessed with the aid of the minimal and average mesh quality numbers given by: where q(E j ) represents the quality number of element E j and N e is the total number of elements.

Element Stitching Operation
In the local element transformation step, each element is transformed separately without considering the connection relationship with its neighbors. Therefore, global mesh optimization is accomplished by averaging the temporary node positions obtained by the local element transformation. With the element stitching operation, the mesh quality should be improved gradually, and the mesh elements should be transitioned smoothly and uniformly throughout the whole computational domain.
In this subsection, global indexing will be used to describe the mesh smoothing and untangling approach. Here, the nodes in the k th iteration step of the approach are represented by p k i with i ∈ {1, 2, ..., N}, where N is the total number of mesh nodes. E j represents the j th element in the mesh, and E(i) represents the index set of all elements containing node p i .
At the first step, all elements are transformed simultaneously using the Equations (1)-(4). The scaling factor is set to: for all elements after several attempts, where q is the mean ratio quality number of the corresponding element. A large value of scaling factor λ may increase the risk of inducing inverted elements, so the coefficient 0.2 is used to prevent the element from transforming much too rapidly. Here, an element with poorer quality will undergo a greater degree of transformation in order to improve the convergence speed.
The next step is computing the new position p k+1 i of each node. The elements containing node p i will produce |E i | new temporary nodes p k ij after the simultaneous transformation. The position of node p k+1 i is updated according to the following formula: in which E(i) denotes the index set of elements containing node p i , p k ij is the temporary node corresponding to node p k i in the scaled transformed element E j and w j is the weight of element E j . A proper choice for element weights should put more emphasis on low-quality elements, especially on inverted elements. In addition, elements closer to the boundary should have greater weights, so that the boundary constraints can be used to gradually eliminate inverted elements, and the mesh elements should be transitioned smoothly and uniformly throughout the whole computational domain. In order to satisfy these requirements, the following weight function is introduced: where f is the distance function formulated by: and g is the quality function written as: The symbol d j stands for the Euclidean distance between the centroid c j of element E j and the associated boundary node p bj of element E j . Boundary node that has the shortest distance to the centroid c j is selected as the associated node p bj .l is the average of the mean edge length of all elements in the mesh, and c is a power exponent which is usually set to a value of c = 2.0 in order to smooth the distance-decay effect. In this way, all mesh elements are roughly sorted by their distance to the boundaries. As the proposed approach is iterative, this makes the mesh elements transition smoothly and uniformly throughout the whole computational domain. To be honest, it is very time-consuming to update the distance function of all elements at each iteration step. In practice, the distance function value of each element can be updated every few iterations.
According to the expression of quality function, low-quality elements especially inverted elements account for larger weights in determining the final node positions. This is designed to accelerate the convergence speed of the approach and to improve the ability of mesh untangling. Generally speaking, increasing the value of exponent η will lead to an improved worst element quality as well as a slightly lower average element quality. Therefore, η = 2.0 have been used as a default value to obtain a balanced optimization result.
Note that the default values for η = 2.0 and c = 2.0 usually work well. However, in some cases, they still need to be carefully selected in order to produce a better result.

Untangling and Smoothing Pipeline
A detailed description of the proposed simultaneous untangling and smoothing approach is given in this subsection. Algorithm 1 provides the pseudo-code of this EEGT-based simultaneous untangling and smoothing approach for isotropic 2D meshes.
The element transformation and node positions updating steps are performed iteratively until the current mesh does not contain inverted elements. In step 3 new node positions are not immediately updated, since this would influence the computation of subsequent nodes. Therefore, new node positions computation and applying the new coordinates are separated, ensuring that the computation does not depend on the numbering scheme of nodes and elements, so as to obtain reproducible results. This is important for parallelized implementation of the proposed approach.
After the above mesh untangling process (step 2-4), the mesh contains no inverted elements. However, quality of the mesh obtained by the untangling process may be low, and further optimization is needed. As the approach is geometry-driven, the steps element transformation and node positions updating may lead to the generation of inverted elements. Therefore, it has to be checked for every element with its new nodes during the mesh smoothing process (step 5). All mesh elements, which are inverted, get a note, that their nodes have to be reset. After all elements are checked, all nodes p k+1 i , belonging to one or more inverted elements, are reset to their old coordinates, that is, p k+1 i = p old i . These two steps are repeated iteratively until no invalid elements are left. The described method of handling inverted elements during the mesh smoothing process is simple and parallelizable.
where w j is the weight of element E j ; (b) Store old node positions p old i = p k i and update all inner nodes. 4. Determine the validity of the current mesh Let Count = Count + 1, and update the index set J inverted : (a) If J inverted = NULL, let Valid = 1, and go to step 5 to further improve mesh quality; (b) If J inverted = NULL and Count is less than the maximum number of iterations in the untangling process, go to step 2 for the next iteration; (c) Otherwise, the proposed algorithm fails to untangle the mesh.

Mesh smoothing:
Keep the boundary nodes fixed, and carry out the simultaneous element transformation, node positions updating and revert nodes of inverted elements to their previous position, until the average mesh quality improvement is below a given tolerance or maximal number of iterations in the smoothing process is reached.
The two main steps of the proposed approach, including element geometric transformation and node positions updating can be easily implemented on GPU due to their fine-grained parallelism. Therefore, a GPU-based implementation is also completed by the authors.

Examples
In this section, three tangled meshes will be optimized to evaluate the performance of the proposed approach, and resulting meshes are compared with those obtained by the default untangle wrapper provided in the mesh quality improvement toolkit Mesquite [22]. Mesquite provides three default untanglers based on different metrics, namely Untangle-Beta untangler, Size untangler and Shape-Size untangler. As the three default untanglers in Mesquite are mainly used to eliminate inverted elements, the quality of the mesh after untangling is generally poor, so the Shape Improve Wrapper in Mesquite is used to further improve mesh quality. At last, examples for testing the performance of GPU implementation are given. In this section, the mean ratio (MR) quality metric is also used to compare the quality of meshes obtained by different methods, and the quality number of inverted elements is set to 0.
For the proposed approach, mesh smoothing process (step 5) is terminated if the average mesh quality improvement of two successive meshes is less than 0.001. The described untangling and smoothing approach is implemented with the C++ language and all the tests are executed in an Intel Xeon E5-2637 CPU at 3.50 GHz with 32 GB RAM.

Quadrilateral Mesh Example
We first consider a square quadrilateral mesh with a hole. The initial mesh is shown in Figure 6a. The mesh has 168 nodes and 140 quadrilateral elements, including 23 inverted elements. There are 56 boundary nodes in the mesh, and they are fixed during the optimization process. The three default untanglers in Mesquite and the proposed approach (EEGT) in this paper are used to optimize the initial tangled mesh. Table 1 shows mesh quality statistics and the number of inverted elements of the initial mesh and output meshes, and Figure 6 shows output meshes obtained by different methods. The Untangle-Beta untangler successfully eliminates all inverted elements in the initial mesh. The average quality of the mesh is also high enough, but the worst quality of mesh elements is only 0.0207. There are many low-quality elements near the left inner boundary, which can be seen in Figure 6b. Such a mesh can not be directly used in numerical simulations, so the Shape Improve Wrapper in Mesquite is used to optimize the output mesh. The resulting mesh is shown in Figure 6f. The average quality of the mesh is improved to 0.6312, and the quality of the worst element is improved to 0.2297. It can be seen from Figure 6f that after further optimization, the distribution of mesh elements is more uniform in the calculation domain. In this example, the Size untangler does not eliminate all inverted elements, and the resulting mesh still contains 4 inverted elements. The Size untangler tries to make the size of the mesh elements uniform while untangling the mesh, which can be seen in Figure 6c. Different form the Size untangler, the Shape-Size untangler considers both the element size and shape, so it can perform mesh untangling and improve the mesh quality at the same time. Figure 6d shows that the optimized mesh of the Shape-Size untangler is better than those obtained by the other two untanglers. Nonetheless, the optimized mesh still has some low-quality elements near the left boundary. Figure 6e shows the optimized mesh obtained by using the proposed approach. Although the average mesh quality is lower than that obtained by the optimization-based method, the quality of the worst element is higher. In addition, the distribution of elements in the calculation domain is sufficiently uniform. The proposed approach has obvious advantages over the optimization-based methods in terms of running time. The Size untangler takes 0.09 s to eliminate inverted elements, while the Shape-Size untangler takes 0.239 s. The Shape Improve Wrapper takes extra 0.435 s to smooth the untangled mesh obtained by the Untangle-Beta untangler, so the total running time is 0.525 s for Mesquite to optimize the initial mesh. The proposed approach only takes 0.006 s to untangle and smooth the initial mesh, which shows the high efficiency of this method.

Mesh Deformation Example
In computational fluid dynamics simulation, the computational mesh should conform to the changing flow domain at each time step during the simulation process. However, mesh deformation methods may lead to the generation of inverted elements for large mesh deformations. Mesh smoothing and untangling algorithms are of great importance for improving the deformation capacity and deformation effect of mesh deformation methods under large deformations.
In this example, an unstructured NACA0012 airfoil mesh with 5892 nodes and 11,363 triangular elements is tested. The size of the square domain is 20 × 20, while the chord length of airfoil is 1. The airfoil moves 6 units to the both down and right while it rotates 75 degrees counterclockwise around its aerodynamic center, which is a very typical case for mesh deformation. The deformed mesh obtained using the mesh deformation method proposed in Reference [33] is shown in Figure 7a, and the deformation is accomplished in a single step. Due to the large deformation, the deformed mesh has 459 inverted elements. The three untanglers in Mesquite and the proposed approach are used to optimize the deformed mesh. The optimization results are given in Table 2, and the optimized meshes are shown in Figure 7.
None of the three untanglers in Mesquite can untangle the deformed mesh, and the number of inverted elements in the resulting meshes is even larger than that of the initial deformed mesh. Among them, the Size untangler has the worst optimization effect, and some nodes in the output mesh even locate outside the initial boundary.
The proposed approach successfully untangle the deformed mesh. The average and minimal quality of the optimized mesh are 0.8102 and 0.4157, respectively. Unlike the local optimization-based method, the proposed approach can make full use of the boundary constraints to optimize the positions of interior nodes, and element geometric transformation and element stitching operations can improve the mesh quality iteratively during the optimization process.
Of course, there is still room for further optimization of the mesh obtained by our method. Therefore, the Shape Improve Wrapper in Mesquite is used to optimize the mesh obtained by the proposed approach, and the average quality of the mesh is improved from 0.8102 to 0.9006, while the quality of the worst element is reduced to 0.2963. Optimization-based methods usually produce meshes with higher average quality, but they are also more time-consuming. In this example, the Shape Improve Wrapper takes 108.089 s to further optimize the mesh obtained by the proposed approach, while the proposed approach only takes 29.849 s to untangle and smooth the initial deformed mesh.

Triangular Mesh Example
The final example is a triangular mesh with two inner boundaries. The initial mesh is shown in Figure 8a. The mesh is generated by randomly moving the positions of interior nodes, and it has 5383 nodes and 10083 triangular elements, 4914 of which are inverted. The Size and Shape-Size untanglers in Mesquite and the proposed approach (EEGT) are used to optimize the initial tangled mesh. The Untangle-Beta untangler fails to output the optimized mesh within 1800s, so its optimization result is not given in this example.
The output meshes are shown in Figure 8. Both the Size and Shape-Size untanglers fail to untangle the initial mesh, and the optimized meshes of these two untanglers contain 1825 and 2636 inverted elements, respectively. All of the inverted elements are successfully removed by the proposed approach. Table 3 shows mesh quality statistics and the number of inverted elements of the initial mesh and output meshes of different methods. It can be seen that the proposed approach is effective in eliminating inverted elements and producing meshes with high quality, even for highly tangled meshes. Table 3. Mesh quality statistics and the number of inverted elements for triangular mesh example. The mean ratio quality metric is used to measure the mesh quality.

Methods
Time (

Implementation on Gpu
In recent years, the programmable Graphic Processor Unit (GPU) has evolved into a highly parallel, multi-threaded processor with tremendous computational horsepower and very high memory bandwidth. More specifically, the GPU is especially well-suited to address problems that can be expressed as data-parallel computations, that is, the same program is executed on many data elements in parallel. Due to fine-grained parallelism, the proposed approach is accelerated by GPU to demonstrate its potential. The GPU used in this test is NVIDIA Quadro K2200, and its detailed parameters are listed in Table 4. Again, mesh smoothing process is terminated if the average quality numbers of two consecutive meshes deviate less than 0.001. Optimization results on CPU and GPU are given in Table 5. More than 10 times speedup is achieved by GPU implementation in the test platform, which shows the great potential of the proposed approach in solving large-scale mesh optimization problems. With the increase of the number of mesh elements.

Conclusions and Future Work
A novel and efficient geometry-based approach for simultaneously untangling and smoothing isotropic 2D meshes has been proposed for the first time in this paper. It has the following features: (1) The proposed approach is based on a new explicit element geometric transformation operation, which can transform an arbitrary polygonal element into its regular counterpart if applied iteratively. It is more reliable and efficient, and its convergence is proved by theoretical analysis.
(2) The proposed approach involves two main steps, namely independent element geometric transformation and node positions updating. It is very suitable for parallel implementation, especially on GPU.
(3) The proposed approach has a strong ability to eliminate inverted elements, and the optimized mesh elements can be transitioned smoothly and uniformly throughout the whole computational domain.
Limitations. The proposed approach can fail on meshes with extreme difference in element size. It may also fail when the initial positions of interior nodes are fairly random. Both situations may lead to the failure of eliminating all inverted elements, because the node positions updating may lead to worse node coordinates. As the element geometric transformation aims at transforming a polygonal element into its regular counterpart, the proposed approach based on the regularizing element transformation cannot be used to optimize anisotropic meshes.
Future Work. Further study will be focused on the extension of the proposed approach to volume meshes and anisotropic meshes. Beyond that, it can be combined with topology modification techniques and local mesh regeneration techniques to form a more reliable method for untangling and smoothing of 2D meshes.
Author Contributions: Conceptualization, S.S. and Z.G.; methodology, S.S. and Z.G.; software, Z.G.; theoretical analysis, Z.G. and M.G.; writing-original draft preparation, Z.G.; writing-review and editing, S.S., Z.G. and M.G. All authors have read and agreed to the published version of the manuscript.
Funding: This research was funded by National Key R&D Program of China grant number 2017YFC0803300 and 2018YFC0809700.