Local Element Operations for Curved Simplex Meshes

We introduce the distortion measure for linear elements (2). This measure – known as the inverse mean-ratio shape measure – was first introduced in [33, 34]. This measure is widely used in the meshing community as a shape distortion metric due to its many desirable properties such as the standard invariance to translation, rotation, and scaling, but is also special as a geometric shape measure that also enjoys many of the desirable properties of an algebraic mesh quality metric [35]. In particular, this measure is shown to be convex [36] and well-suited for optimization.

As an aside, we note the metric we are about to introduce is geometric in nature, as are most metrics in use today, as opposed to being PDE or solution based. In the course of using mesh generation and adaptation in the numerical solution of PDEs, we point out that using a metric that is geometric in nature may not perfectly align with the ultimate goal of solving the PDE since the metric employed was not based on the PDE solution. This is addressed in [37] where the author discusses the notion of “bridging the gap between a priori quality metrics and solution-dependent metrics”.

Consider a linear element $E^{I}$ having the desired shape and size (the ideal element) and a corresponding element in physical space $E^{P}$. From these we can define a unique affine mapping $\phi_{E}:E^{I}\rightarrow E^{P}$. It is frequently convenient to introduce the notion of the master element $E^{M}$ and define two additional mappings $\phi_{I},\phi_{M}$ (Figure 3a) that map from the master element to the ideal and physical element, respectively. We can define $\phi_{E}$ in terms of the these mappings as follows

The element shape distortion metric $\eta$ is defined as

where $d$ is the spatial dimension, $||\cdot||_{F}$ is the Frobenius norm, and $\sigma=\det(D\phi_{E})$. Note that for linear elements, since $\phi_{E}$ is affine, its Jacobian $D\phi_{E}$ is constant. This distortion measure quantifies the deviation of the shape of the physical element $E^{P}$ with respect to the ideal element $E^{I}$. This distortion measure is 1 for the ideal element and tends to infinity as the physical element degenerates. It is often convenient to define a corresponding quality measure $q$ as

For high-order elements, we can no longer apply (2) directly because the mapping $\phi_{E}$ is no longer affine, so $D\phi_{E}$ is no longer constant. Roca et al. describes an extension of the linear distortion measure to high order triangular [20] and tetrahedral [21] elements. The idea is that since the linear measure quantifies the local deviation between the ideal and physical elements, we can obtain a high-order element-wise distortion $\hat{\eta}$ measure by integrating the linear measure on the physical high order element

where $|t|$ is the area of the physical element and $r=2$ in this work. Explicit expressions for the isoparametric mappings for the high-order case (Figure 3b) in terms of shape functions are given in [20].

The distortion measure $\eta$ needs to be modified to eliminate the singularity so an optimization procedure can recover from an invalid initial configuration. We use the regularization procedure introduced independently by [31, 32] around the same time, but follow the notation [32] which replaces $\sigma$ in the denominator of (2) by

where $\delta$ is a positive element-wise parameter. This regularized Jacobian $\sigma_{\delta}(\sigma)$ is a monotonically increasing function of $\sigma$ such that $\sigma_{\delta}(0)=\delta$, which tends to 0 when $\sigma$ tends to $-\infty$, which allows us to overcome the vertical asymptote at $\sigma=0$ in the original measure. This parameter $\delta$ is only set for a non-zero value if there exists an invalid element in the mesh under consideration, otherwise it is set to zero for a mesh with all valid elements. Using the regularized Jacobian, we can now modify (2) to obtain a linear distortion measure capable of simultaneous smoothing and untangling

and use this to define a high-order regularized measure $\hat{\eta}_{\delta}$ analogous to (4). Figures 5 and 6 of [19] provide a useful illustration of how this regularized distortion metric compares to the original for a simple example.

What remains is the issue of how to pick this positive elementwise regularization parameter $\delta$. We certainly require the original ($\hat{\eta}$) and regularized ($\hat{\eta}_{\delta}$) distortion measures to have nearby minima, so $\delta$ needs to be sufficiently small. On the other hand, $\delta$ has to be large enough to ensure that $4\delta^{2}$ is significant compared to $\sigma^{2}$ in the expression for $\sigma_{\delta}$. When this idea was originally introduced in [32], this parameter was chosen in a fairly ad hoc manner by testing a few values for a given initial tangled mesh. In [19], a heuristic was developed to choose a constant value of $\delta$ for each element. Their setting deals with high-order mesh generation by curving the boundaries of a well-shaped straight-sided mesh and minimizing the regularized distortion metric. For them, in the mapping $\phi_{E}$, each physical element $E^{P}$ in the curved mesh has the ideal element $E^{I}$ as the corresponding element in the original straight-sided mesh. They propose to choose this parameter $\delta$ solely based off information from the straight-sided ideal element. Defining $\sigma^{*}=-\det\phi_{I}$ and imposing

for some given tolerance $\tau$ implies

They argue that $\tau$ should be small compared to $\sigma$ and select $\tau=\alpha|\sigma^{*}|$, giving the final value for $\delta$ as

The value $\alpha=10^{-3}$ is observed to work well in practice and accomplish the tradeoff required on the value of $\delta$ (Figure 4). Our setting deals with high-order mesh smoothing, so we choose the ideal element $E^{I}$ to always be the equilateral triangle/tetrahedron. For the unregularized distortion measure, the actual size of $E^{I}$ is irrelevant since the measure is invariant to scaling. This is no longer true for the regularized distortion measure. So in order to maintain the assumption from their setting that each $E^{I}$ and $E^{P}$ are roughly of the same size, we need to scale the ideal element to have the same volume as each corresponding physical element and choose $\delta$ accordingly.

We need to aggregate the distortion measures for individual elements into a single quantity for a group of elements for optimization purposes. What we really want is to optimize for the worst element, since the distortion of the worst element in a mesh is often what drives the solution quality in the context of numerical simulation. But since minimizing the maximum distortion is difficult optimization-wise, we consider the normalized sum of squares of the elementwise distortion of each element $e_{i}$ in a group of elements $M$

where $\hat{\eta}_{\delta}$ is the regularized high-order distortion measure obtained by combining (4) and (6)

For smoothing over the entire mesh, we largely follow the localized approach given in Appendix A of [19], which considers a “Gauss-Seidel” like method which iteratively updates each node by a Newton step with backtracking line search until either the objective function or the amount of node movement satisfies some stopping criteria. Instead of looping over nodes, we loop over patches (Figure 5), which seems to help mitigate some of the “back and forth” the optimization procedure can empirically observed to become stuck in, which is also observed in [18] for a different high-order mesh smoothing procedure. We note this solver offers no guarantees for especially tangled or pathological meshes due to the issues with the regularization method which many authors have pointed out fails quite often in practice such as if the boundary of the mesh is also tangled [38], or the possible failure of the high-order Gauss quadrature rule to detect inversion at the corners of our high-order elements.

Once a group of elements is identified to perform edge/face swaps on, we first perform the straight-sided initial guess. This is done by freezing the boundary, discarding all the interior nodes and replacing them in the flipped configuration by the straight-sided information. This initial flip is extremely prone to tangling even when starting out with fairly well-shaped elements, which is more severe for the 3D case. By the high-order distortion measure (4), many high-order elements which seem valid by visual inspection with all interior nodes contained within the boundary may very well be invalid. Finally, we apply the simultaneous untangling and smoothing procedure described in the previous section to arrive at the curved flip (Figures 9 and 9). We decide to accept the resulting flip if it results in a lower aggregate distortion $\hat{\eta}_{agg}$ from the original configuration. The normalization by the number of elements in the group in the definition is to allow comparison of the original and flipped configuration, since in 3D there are many topological changes which in general may alter the number of elements.

Once an edge is identified to collapse, we first do the straight-sided initial guess. This is done by identifying the patch of elements that contains either endpoint of the edge being collapsed, freezing its boundary, discarding all the interior nodes, and replacing them by the straight-sided edge collapse to the midpoint of the original curved edge. Then we apply the simultaneous untangling and smoothing procedure to arrive the final curved collapse (Figures 9 and 9). If the edge identified for collapse is on the boundary of the domain, information will inevitably be lost but can be mitigated in the high-order case by performing projection based interpolation on the patch boundary by standard techniques encountered in $hp$-FEM methods [39].

Depending on the setting, the edges to be split could either be directly identified or implied by choosing elements. Either way, once an element is identified to be split and the straight sided guess is performed according to one of the templates in 2D (Figure 13) or 3D (Figure 13), its neighbors must also be split by a straight-sided guess with the corresponding template to prevent hanging nodes. Then, we apply simultaneous untangling and smoothing on this group of elements to arrive at the final curved split (Figures 13 and 13). There are many more possible subdivision templates in 3D depending on the number of marked edges (Figure 11 of [9]), and that for any of these templates to determine how to split the edge(s), the amount of “smoothing” required after the edge splits applied to the straight sided guess can differ quite greatly and alter the overall efficiency of the method, but this is outside of the scope of this work.

We first demonstrate how node smoothing combined with edge flips can maintain element quality under rotation. The 2D mesh is a circle-in-square mesh of 116 elements, where the circle is centered at the origin with radius 0.5 inside the square $[-1,1]^{2}$ (Figure 14 left). The analogous 3D mesh is a sphere-in-cube mesh of 3395 elements, where the sphere is centered at the origin with radius 0.5 inside the cube $[-1,1]^{3}$ (Figure 15 left). At each step, we rotate the nodes on the circle/sphere, perform high-order node smoothing while keeping the nodes on the circle/sphere fixed, apply curved flips according to Algorithm 1, and smooth again. Algorithm 1 constructs a priority queue of elements with higher distortion than a floating threshold value and considers all topologically possible flips for each element, accepting the one (if any) that results in the greatest improvement. It performs sweeps over all the elements with a reset floating threshold for each sweep until no further improvement is possible.

Rotation cannot exceed much more than a quarter turn with only node smoothing (Figures 14 and 15 middle). A straight sided mesh cannot rotate this far with smoothing alone; the presence of high order nodes grant us the flexibility to maintain valid (albeit highly stretched) elements for more severe deformations. It is well known in the straight-sided case that such greedy flip algorithms like the one used here converge to an optimal triangulation (in the sense of Delaunay) in the 2D case, but offers no such guarantees in higher dimensions [40]. 2D rotation can continue indefinitely, and we can reasonably attain a good aggregate quality under rotation for the 3D case, but some poorly shaped elements emerge as evidenced by the growing maximum element distortion (Table 1); this issue is no doubt magnified by our use of very coarse meshes. We also report the aggregate mesh quality (9), and the numerical results also illustrate a fundamental issue of mesh smoothing discussed at the beginning of Section 2.3. It turns out even if the aggregate distortion of the mesh is very low, this can obscure the fact that there are elements of very low quality present as well. This can also be remedied with more sophisticated mesh cleanup procedures featuring more specialized element operations such as compound operations [41] and sliver exudation [42], but such issues are not the focus of this work. Also note that here we do not nodes on the boundary or the circle/sphere. Allowing these nodes to slide on the fixed boundaries while smoothing the entire mesh would certainly be an improvement to the procedure.

Though this example does not deal with mesh deformation, we demonstrate how curved collapses can be combined with high-order node smoothing to coarsen a mesh (Figure 16). We begin with a rather fine 2D mesh with 690 elements of the square $[-1,1]\times[-1,1]$ with a circular hole at the origin with radius 0.5 (Figure 16, top left). We apply multiple passes of curved edge collapses to coarsen the mesh until convergence according to a chosen uniform size function (Algorithm 2). In this case we choose the ideal edge length $L_{0}=0.1$ as a parameter far enough from the original edge length just to illustrate the procedure to induce collapses. We also coarsen elements on the interior boundary with projection by interpolation in the $L^{2}$ error minimizing sense to maintain an accurate representation of the circular boundary with fewer elements.

We demonstrate the combination of all three curved element operations with mesh smoothing to maintain both element quality and uniform sizing. The 2D mesh is a circle-in-rectangle mesh of 314 elements, where the circle is centered at $(-1,0)$ and has radius 0.5 inside the rectangle $[-2,2]\times[-1,1]$. The analogous 3D mesh is a sphere-in-prism mesh of 3395 elements, where the sphere is centered at $(-1,0,0)$ with radius 0.5 inside the prism $[-2,2]\times[-1,1]\times[-1,1]$. At each step, we translate the nodes on the circle/sphere, perform high-order node smoothing while keeping the nodes on the circle/sphere fixed, apply all three types of local mesh operations according to Algorithm 4, and smooth again. Algorithm 4 applies sweeps of splits (Algorithm 3) and collapses (Algorithm 2) until convergence, and then applies sweeps of flips (Algorithm 1) until convergence. In terms of the parameters for both Algorithms 2 and 3, the ideal edge length is set to be $L_{0}=\text{avg}(L)$, where $L$ is the list of computed original edge lengths. This choice is made simply because after the translation, we want the final mesh to resemble our original mesh, which was roughly uniform and with the same element sizes. As an aside, one could imagine introducing anisotropy to the mesh through these local element operations through a specified size function which would specify the ideal lengths $L_{0}$ at any point in space, but this is beyond the scope of the work.

Translation cannot make it more than halfway across the prism with only node smoothing (Figures 17 and 18). With the addition of flips, we can make it across the entire domain, but this results in elements of highly variable size. This can be remedied with the addition of collapses and splitting to maintain a uniform size, which allows us to practically maintain the original element sizes and qualities after translation across the entire domain (Table 2).