WIXLIB

CSRI Summer Proceedings 20093QUADRATIC TRIANGLE MESH UNTANGLING AND OPTIMIZATION VIA THETARGET-MATRIX PARADIGMPATRICK KNUPP , NICHOLAS VOSHELL†, AND JASON KRAFTCHECK‡Abstract. Meshes containing high-order nodes sometimes contain inverted elements due to the projection ofcoarse elements onto the domain boundary. Usually, such meshes are problematic and fixing the associated problemsrequires either re-meshing or post-processing techniques. Mesh optimization methods based on the Target-matrixparadigm are studied as a post-processing step to determine if the approach gives a viable solution to both untanglingand improving the quality of a quadratic triangle mesh. Several preliminary algorithms are described, including oneinvolving two successive optimizations. For that algorithm, a first optimization using a non-barrier size metric witha properly constructed target-matrix can potentially untangle the mesh. If the first optimization untangles the mesh,a second optimization of that result using a barrier-based shape metric can improve shape while keeping the meshuntangled. Numerical experiments applying the algorithms on planar quadratic triangle meshes demonstrate thatoptimization-based node-movement methods can successfully untangle and improve element shapes in high-ordermeshes. Further study of these algorithms is needed before they can be used with confidence on realistic meshes.1. Introduction. To achieve high accuracy simulations, second-order finite elementsare sometimes used. Unfortunately there are no reliable tools for meshing curved geometries known to the authors; however mesh generation packages usually follow an a posterioriapproach to create boundary-conforming, quadratic element meshes on complex geometriesby first creating elements with straight sides and then curving boundary sides by projectingmid-side and mid-face nodes onto the bounding geometry [11]. The quality of the resultingboundary elements can be poor, particularly if elements are large compared to the curvatureof the associated geometry. This impacts simulation accuracy, efficiency and, if the Jacobianis negative, even the ability to proceed with the calculation.The convergence theory of finite elements requires that mesh elements not be inverted.Most elements are defined based on a mapping from a logical element and the Jacobian matrix, J, of this mapping. The strict definition of an inverted element is that there exists a point(ξ, η) within the logical element such that det(Jξ,η ) 0. If such a point does not exist, thenthe element is said to be non-inverted. Given an arbitrary linear planar triangle or quadrilateral, it is simple to determine from this criterion whether the element is inverted because thedeterminant is linear in the logical coordinates. However, for hexahedra, quadratic, and otherelement types, it is difficult to find robust and efficient numerical algorithms for determiningwhether a given element is valid.The present work does not focus directly on the detection of inverted elements. Since it ispossible to improve mesh quality without knowing beforehand if the mesh contains invertedelements, the question to be considered here is: given a mesh of quadratic elements, can weimprove shape and size quality, such that the result rarely contains inverted elements? It isrecognized that edge flipping and node insertion are powerful techniques in their own right.However, the approach taken here is based on a node movement strategy called the Targetmatrix paradigm [6]. It is hoped that, once it is understood how to effectively optimize highorder element meshes using node movement alone, all these techniques can be combined intoa single algorithm. The capability to optimize the quality of high-order finite element meshesby node-movement was recently added to the Mesquite mesh quality improvement library[1]. The present study uses the new capability to explore how one might improve meshescontaining high-order nodes. For simplicity, attention is confined to planar quadratic triangle SandiaNational Laboratories, pknupp@sandia.govPennsylvania State University, njv116@cse.psu.edu‡ The University of Wisconsin, kraftche@cae.wisc.edu† The

4Quadratic Triangle Mesh Untangling and Optimizationmeshes, but the approach is generalizable to other types of high-order elements. Studiesinvolving other element types are planned for future work.There have been limited studies that attempt to untangle or improve element shapeswithin meshes that include high-order elements (which are elements that have polynomials of degree greater than one defining their shape). Investigations defining valid regionsfor the placement of mid-nodes to ensure untangled elements (and other properties) werecalculated in [13, 14], this could help in detecting untangled elements, or guiding vertexplacement. In [9, 10], a mesh modification technique is employed to improve the quality ofa high order element mesh (as well as enabling refinement and/or coarsening of the mesh).The technique was developed using vertex insertion and removal, flipping of edges and faces,along with other topological mesh modifications. One method for assessing the quality ofquadratic triangle elements (that could be extended to give quality metrics for mesh optimization) was given in [12]. Proposed optimization-based node-movement strategies for improving quadratic meshes appear to be limited to [2, 3]. In [3], a linear quality metric wasmodified to be sensitive to high-order node positions by adding an angle-based penalty term.The quality metric was extended to quadrilaterals in [2] and a node-movement strategy thatautomatically switched between constrained Laplace smoothing and an optimization procedure based on the linear quality metric and the penalty term was applied to improve quadraticmeshes. Limitations of the method include being applicable only to homogeneous, isotropictwo-dimensional meshes.2. Quadratic Elements. Using the standard basis vectors for R2 one can define a logical triangle which serves as the master element for the transformations used by quadratictriangles. The other elements that are of interest in this work are the active element (which isthe element of interest in the mesh) and the target element (the desired element), which theactive element is compared to. Any given quadratic triangle can then be defined by six points,corner vertices (x0 , x1 , and x2 ) and mid-edge nodes (x3 , x4 , and x5 ). A mapping from point(ξ , η) of the logical triangle to point X of the quadratic triangle can then be defined using thefollowing equations:X({xi }, ξ, η) x0 c1 ξ c2 η c3 ξη c4 ξ2 c5 η2(2.1)c1 3x0 x1 4x3(2.2)c2 3x0 x2 4x5(2.3)c3 4(x0 x3 x4 x5 )(2.4)c4 2(x0 x1 2x3 )(2.5)c5 2(x0 x2 2x5 )(2.6)The vertices of the quadratic triangle are enumerated 0, 3, 1, 4, 2, 5 in counter-clockwiseorder (both vertices and mid-side nodes are in counter-clockwise order with mid-side node 3following vertex 0). The Jacobian matrix for this transformation can be defined at each pointwithin the logical triangle using the following equation: J c1 c3 η 2c4 ξ , c2 c3 ξ 2c5 η(2.7)Note that a linear mapping can be considered a special case of this in which the mid-nodesare mapped to the mean of the two corner vertices they are connected to. This causes c3 c4 c5 0, which gives a linear mapping function with the following constant Jacobian:J [ c1 , c2 ](2.8)Given an active quadratic element (i.e., one that belongs to the mesh to be optimized) and asample point (ξk , ηk ) within the logical element, we have Ak Jactive (ξk , ηk ). Similarly, given

P. Knupp, N. Voshell, and J. Kraftcheck5a quadratic target element, Wk Jtarget (ξk , ηk ). Thus for every sample point, we have the pairAk , Wk . When speaking in general about active and target-matrices, the sample point index issuppressed, giving A and W.From this one can define the Jacobian matrix for the mapping from the target triangleto the active triangle (the inverse of logical 7 target composed with logical 7 active) asfollows:T A W 1(2.9)For an isotropic domain, there is no reason for the target element to break symmetry, so theequilateral triangle is a good choice for the target. For theoretical reasons explained in [4],we choose the area of the target equilateral element to be one-half. Then the Jacobian of thelogical to target mapping (up to multiplication by some rotation matrix) is as follows:#"12 1(2.10)W 432 3 03. Mesh Optimization. Mesh optimization, as the term is used here, is a techniquefor modifying a mesh by moving vertices without changing the connectivity. In mesh optimization one defines an objective function to quantitatively compare the quality of alternativemeshes. The optimal mesh has the best objective function score subject to constraints (on theposition of boundary vertices, for example). The best valid objective function value is foundusing numerical optimization techniques. The objective function is usually some combinationof measurements of the quality of individual components (such as elements or vertices) in themesh. Most effort, as of this writing, has been on optimizing meshes composed of linear elements. For these linear meshes, good element quality metrics have been found, particularlyto measure the shape of triangular elements.In Mesquite there are four main options for the movement of each mesh vertex, whichform optimization constraints. The four options are (a) fixed (no movement permitted), (b)free (any movement permitted), (c) constrained to the geometry (where the vertex is constrained to a lower dimensional space, such as the boundary), and (d) slaved (where a midside node is constrained to the midpoint of the neighboring vertices). The first three optionsapply to all vertices, while (d) applies only to mid-side nodes.Some of the work on improving mesh quality has dealt with meshes where elements areinverted. For quadratic element meshes, the issue frequently arises during mesh creation,where curving the boundary of the initial linear mesh to conform to the geometry causessome elements to be poorly shaped or even inverted. To eliminate inverted elements fromthe mesh via optimization, a number of mesh untangling algorithms have been proposed(see [5] and the references therein). Unfortunately, none of these algorithms guarantee thatthe resulting mesh will, in fact, be untangled. Furthermore, by focusing exclusively on theinversion issue, the untangling techniques ignore other important aspects of mesh quality.To eliminate inverted elements and improve shape in quadratic element meshes we take adifferent approach that does not use a direct untangling algorithm. The approach in this workuses local quality metrics from the Target-matrix paradigm.In the Target-matrix paradigm, each element of the mesh has a C 1 mapping from a targetelement to the active element. The Jacobian of this mapping (evaluated at sample point k) isthe matrix, T k . Then the quality at k is measured by local quality metric µk µ(T k ), whichgives a non-negative real number. The metrics measure the degree to which A is ’close’ to W,in terms of properties such as shape, size, and orientation.There are two basic types of local metrics in the Target-matrix paradigm, those having’barriers’ and those which do not. Barrier metrics are used when the initial mesh is not

6Quadratic Triangle Mesh Untangling and Optimizationinverted; the barrier guarantees that the optimal mesh is not inverted at the sample points.Non-barrier metrics are used when the initial mesh contains one or more inverted samplepoints. The optimal mesh resulting from a non-barrier metric may or may not contain invertedelements. Since the initial quadratic meshes that we wish to improve are assumed to containelements which are inverted, we use non-barrier metrics to optimize the initial mesh. Ifoptimization with a non-barrier metric succeeds in creating a non-inverted mesh, the resultcan be further optimized using a barrier metric.The first non-barrier metric of interest is the Shape (S) metric given by:µS (T ) kT k2F 2 det(T )(3.1)In [6, 7] it was shown that µS 0 for all T and µS 0 if and only if T is a scaled rotationmatrix. In that case, A has the form A sRW with arbitrary scalar s and arbitrary rotationR. Then the shape of the active matrix is the shape of the target-matrix, which in the presentcase corresponds to the equilateral triangle. This is a non-barrier metric, so it can potentiallyuntangle an initially inverted mesh, as well as improve shape. The barrier form of the shapemetric isµS B (T ) kT k2F2 det(T )(3.2)The second non-barrier metric of interest is the Size (Sz) metric, given by:µS z (T ) (det(T ) 1)2(3.3)The metric obeys µS z 0 for all T and µS z 0 if and only if det(T ) 1, i.e., det(A) det(W).Thus, at the minimum, the local areas of the active and target matrices agree. Because thismetric has no barrier, it can potentially untangle an inverted mesh as well as improve relativesize. It does not, however, encourage the shape of the active element to be close to the shapeof the target element. Thus, its primary use in the present application is to encourage meshuntangling because det(W) 21 .Finally, we consider the Shape Size (SS) metric, and it’s barrier form:qµS S (T ) kT k2F 2 kT k2F 2 det(T ) 2(3.4)µS S B (T ) 1µS S (T )2 det(T )(3.5)It was shown in the previously-cited references that the metric is non-negative and is minimized for A R W so that, at the minimum, both the shape and size of A agree with the shapeand size of the Target-matrix. The metric can also be used to encourage the untangling ofmeshes.Many of the issues involved in measuring sample point quality within an element aredescribed in [8]. That work specifically considers quality metrics of the form investigatedhere, including results for quadratic triangular elements. It also introduces a relevant propertyknown as label invariance. Note that x0 can be assigned to any corner vertex, thus there arethree labellings for the corner vertices and mid-side nodes of a quadratic triangle elementthat conform to the naming conventions. A local quality metric is label invariant if all threelabellings of the active element give the same quality. Different conditions are derived for themetric formulation, sample point selection, and target element selection that guarantee labelinvariance of the local quality. Specifically, if the target element is a straight edged equilateraltriangle and the sample points are selected to respect the symmetry of the element labellings,then the resulting quality metric will be locally label invariant.

P. Knupp, N. Voshell, and J. Kraftcheck7An important issue concerns how to select the sample point locations within an element.As mentioned in the Introduction, detecting inverted elements requires, in principle, an infinite number of sample points. However, the present goal is not to detect inverted elements,but simply to improve quality and attempt to untangle inverted elements. Thus, we investigatewhether it is sufficient to use a small number of sample points for this task. In particular, wechoose to place sample points at all the corners and mid-side nodes of each triangular elementin the mesh.1 Numerical experiments will reveal whether this collection of sample points issufficient for untangling realistic meshes.Qualities, µk , at different sample points are combined into an objective function to beused in our numerical optimization procedure. To do so, this study uses a power-mean, givingthe following objective function:X1µk(3.6)OF N S ample Point k4. Numerical Examples. Because this is a preliminary investigation, the issue of algorithm efficiency is not considered. Instead, algorithm effectiveness (robustness) is investigated. The numerical optimizations were carried out using Mesquite’s Feasible Newtonsolver, global patch smoothing, and a termination criterion of 400 iterations. The solutionswere well-converged on the small meshes used in this study.The a posteriori approach to quadratic mesh generation tends to localize poor quality (orinverted) elements on a curved boundary (rather than the interior of the mesh), which guidesexperiment selection. The initial experiment involves a small patch of six quadratic elementssharing a vertex, with a fixed boundary (see Figure 4.1), where the two triangles at the bottomof the patch are curved so that one element is initially inverted.The figure depicts an investigation into which metrics and which combinations of free,fixed, and slaved vertices give the best results. The pictures are arranged in normal (left toright and top to bottom) order with the initial configuration first, then the results, with all midside nodes slaved, optimizing for Size, Shape, and Shape Size (respectively). All of thesemeshes are inverted, however optimizing for Size comes closest to untangling, and optimizing for Shape Size comes close to untangling while maintaining better shape. After thesecome the results optimizing for the Shape Size metric, but freeing some mid-side nodes,first freeing only the interior mid-side nodes in the initially inverted element, then freeing allthe internal mid-side nodes. Both of these final two techniques untangled the mesh (with theincreased freedom creating more curved edges but better shape).After modifying the geometry of this mesh a few times it was found that the Size metricuntangled more cases than the Shape Size metric (which untangled more cases than theShape metric). An instructive investigation performed was one in which the leftmost vertex ofthe small patch was moved to the right and the mesh was optimized for different metrics withall internal mid-side nodes free. See figure 4.2 for a depiction of the results. A big challenge tountangling this patch is that it significantly reduces the quality of the two leftmost elements(especially the lower one) and perturbing the leftmost vertex exacerbates this competition(making untangling the element more challenging). Optimizing with the Shape Size metricwas unable to untangle perturbations less than 0.4, whereas optimizing with the Size metricwas able to untangle a perturbation of 0.6 (but not 0.7). In both cases the point at which theproblem with untangling first occurs is the lower left vertex (that starts out inverted). Notethat there exists an untangled mesh even when the leftmost vertex has been moved by 0.7, sothere are cases where optimizing for the Size metric is insufficient to find a valid untangledmesh configuration.1 Notethat, in this approach, a mesh vertex can have multiple associated sample points, as can a mid-side node.

8Quadratic Triangle Mesh Untangling and OptimizationFig. 4.1. Small Mesh Experiment: The initial mesh is given in the upper left figure. Then the Shape, Size,and Shape Size metrics are used (with all mid-nodes slaved) in the upper middle, upper

Zhaofang Wen and S. Scott Collis Sandia National Laboratories January 11, 2010 A Department of Energy National Laboratory SAND2010-3083P Sandia National Laboratories is a multi-program laboratory operated by Sandia Corporation, a wholly owned subsidiary of Lockheed Martin Corporation, for the U.S. Department of Energy's National