## 3 Smoothing and Untangling Approach
### 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_{i}^{k}\) with \(i\in\{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:  
\[\lambda=0.2(1-q) \tag{17}\]  
for all elements after several attempts, where \(q\) is the mean ratio quality number of the corresponding element. A large value of scaling factor \(\lambda\) 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_{i}^{k+1}\) of each node. The elements containing node \(p_{i}\) will produce \(|E_{i}|\) new temporary nodes \(p_{ij}^{k}\) after the simultaneous transformation. The position of node \(p_{i}^{k+1}\) is updated according to the following formula:  
\[p_{i}^{k+1}=p_{i}^{k}+\frac{\sum_{j\in E(i)}w_{j}(p_{ij}^{k}-p_{i}^{k})}{\sum_{ j\in E(i)}w_{j}}, \tag{18}\]  
in which \(E(i)\) denotes the index set of elements containing node \(p_{i}\), \(p_{ij}^{k}\) is the temporary node corresponding to node \(p_{i}^{k}\) 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:  
\[w_{j}=f(d_{j})g(q(E_{j})), \tag{19}\]  
where \(f\) is the distance function formulated by:  
\[f(d_{j})=(\frac{\bar{I}}{d_{j}})^{c}=(\frac{\bar{I}}{\parallel c_{j}-p_{bj} \parallel})^{c} \tag{20}\]  
and \(g\) is the quality function written as:  
\[g(q(E_{j}))=\left\{ \begin{array}{ll}5.0&if\quad q(E_{j})=0\\ (1-q(E_{j}))^{q}&if\quad q(E_{j})>0.\end{array} \right. \tag{21}\]  
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}\). \(\bar{I}\) 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.