<figure><img src="image_1.png" /></figure>

HAL
open science

# Mesh Untangling for Problems with Topological Singularities

Vladimir Garanzha, Liudmila Kudryavtseva, François Protais, Dmitry Sokolov

► To cite this version:

Vladimir Garanzha, Liudmila Kudryavtseva, François Protais, Dmitry Sokolov. Mesh Untangling for Problems with Topological Singularities. NUMGRID Conference 2022 - 11th International Conference on Numerical Geometry, Grid Generation and Scientific Computing, Dec 2022, Moscow, Russia. pp.77-102, 10.1007/978-3-031-59652-0_5 . hal-04886999

HAL Id: hal-04886999

https://inria.hal.science/hal-04886999v1

Submitted on 14 Jan 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
**Chapter 1**

**Mesh Untangling for Problems with Topological Singularities**

Vladimir Garanzha¹, Liudmila Kudryavtseva¹, François Protais², and Dmitry Sokolov²

**Abstract** Mesh untangling is still a hot topic in applied mathematics. Tangled or folded meshes appear in many applications involving mappings or deformations. Despite the fact that a large number of mesh untangling strategies was proposed during the last decades, this problem still persists. Recently we have proposed a penalty-based numerical optimization scheme [GKK$^+$21] that provably untangles 2d and 3d meshes with inverted elements by partially solving a finite number of unconstrained minimization problems using efficiency assumption: constant relative reduction of penalty function for fixed penalty parameter. We have identified key engineering counterexamples which violate the efficiency condition and thus block untangling process. When such topological singularities as 2-coverings or general *k*-coverings are present in the mesh or appear during untangling attempt the penalty function reveals sharp parasitic local minima outside admissible set. This problem is particularly vexing when partially constrained mesh deformation problems are considered. We suggest adaptive procedure which adds phantom mesh elements near topological singularities with careful choice of weights and allows to eliminate topological singularities by moving them outside the domain or by mutual annihilation. We demonstrate enhanced stability of suggested untangling technique which has a potential to make untangling a routine operation over meshes.

**1.1 Introduction**

To manipulate a geometric object inside a computer, the most versatile option is to discretize it and represented by a mesh (polygonal surface or a polyhedral mesh). When computing a large deformation (or a map), a mesh may become tangled, i.e. inverted elements can appear. Untangling takes a very important place in mesh generation: it takes a mesh as an input, and moves the vertices to get rid of foldovers.

Dorodnicyn Computing Center FRC CSC RAS · Université de Lorraine, CNRS, Inria, LORIA, F-54000 Nancy, France
Originally related to Arbitrary Lagrangian-Eulerian (ALE) moving mesh approach, the mesh untangling problem considers a simplicial complex with badly oriented elements and optimizes vertices position in a way that is likely to flip misoriented elements.

Invertibility of deformations is one of central themes in mesh generation research. Historically, numerical simulation of hydrodynamic instability of layered structures required sound mathematical foundations behind moving deforming mesh algorithms. In 1966 Winslow introduced mesh generation method based on inverse harmonic maps [Win66]. Meanwhile, 4 years before that, Crowley suggested similar ideas in a classified Los Alamos Lab research report [Cro62] which was made public many years later.

In 1972 Godunov [GP72] suggested to control mesh deformations using composition with prescribed invertible mappings thus combining adaptation and guaranteed invertibility. In 1988 Jacquotte [Jac88] introduced elastic deformations for mesh generation. Meanwhile in 1988 Ivanenko discovered that finite element approximation of the Winslow functional provides an infinite barrier on the boundary of the set of admissible grids [Iva88] which completely suppresses inverted elements in the deforming meshes. He formulated untangling problem as a separate problem of numerical analysis and it took almost a decade to get recognized by the community [CI97].

In 1966 Reshetniak [Res66] introduced the concept of multidimensional mappings with bounded distortion and found precise relations between condition number of Jacobian matrix and shape distortion measure attributed in mesh generation community to Liu and Joe [LJ94]. Note that importance of relations between different shape measures and condition number was underlined in [Knu00].

In 1976 J. Ball introduced his theory of finite elasticity based on the concept of polyconvex distortion energies [Bal76]. He not only proved Weierstrass-style existence theorem for this class of variational problems, but also formulated a theorem on invertibility of elastic deformations for quite general 3D domains [Bal81]. It is important that Ball invertibility theorem is proved for Sobolev mappings and can be applied directly for finite element spaces, i.e. to deformation of meshes, as was pointed out in [Rum96]. In [Gar00, GKU14] it was suggested quasi-isometric hyper-elastic material which unlike known models provides mappings with bounded global distortion (bounded quasi-isometry constant) as minimizers of elastic energy. Invertibility theorem for deformation of this material was established in the 3D case as well.

Since then there is an abundant research on mesh untangling, mentioning just a few [Knu01, FP00, ERM+03, TGRL13]. However the common opinion is that untangling is a very hard problem and algorithms are not robust enough. As a manifestation of frustration over this problem papers on numerical methods on tangled (sic!) meshes start to appear [DS13].

The main idea of untangling approach proposed in [GKK+21] is to regularize the sum of cell shape distortion and cell volume distortion. This penalty-based untangling scheme works very well in practice for fixed boundary problems. In particular, it solves in a black-box manner 100% tests from the huge 2d and 3d set suggested in
# 1 Mesh Untangling for Problems with Topological Singularities

[DAZ$^{+20}$] as a major challenge for computational graphics and augmented reality problems. Authors of [DAZ$^{+20}$] tested all state-of-art untangling algorithms and claimed that all but their algorithm fail the 100% success threshold. Our tests have shown that algorithm from [DAZ$^{+20}$] also breaks down when randomized initial guesses for their test set are used, while the results of algorithm from [GKK$^{+21}$] are not sensitive to initial guesses.

For fixed boundary problems the approach has a sound theoretical basis. Namely, in the aforementioned work [GKK$^{+21}$, §4.2] it was formulated a “finite untangling theorem”, stating that when the set of locally invertible deformations (admissible meshes) is not empty, one can build a sequence of penalty parameters in such a way, that finite number of partial unconditional minimizations make the mesh admissible. Each optimization step is assumed to be a bit idealized one but already not far from the assumptions for engineering optimization algorithms. This theoretical result is valid for *d*-dimensional simplicial meshes and can be generalized to more general non-simplicial elements.

Another important theoretical result introduced in [GKK$^{+21}$] (Appendix C) is that untangling strategy guarantees that positive definite part of the Hessian matrix for discrete untangling functional is spectrally equivalent to the finite element Laplacian. It means that potential instability which may spoil the behavior of untangling algorithm [GK99a] when crossing the barrier is eliminated.

However algorithm from [GKK$^{+21}$] has some limitations in the case of free boundaries. In particular it does not provide guarantees of local invertibility, which is the contribution of the current work.

## Our contributions

We extend above results to the case of *free* boundaries. We guarantee that the set of locally invertible deformations (recall that inversion-free does not imply local invertibility!) can be attained in a finite number of minimization steps. More precisely, we propose an algorithmic scheme that adds a small number of phantom triangles for a 2D mesh and tetrahedra for a 3D mesh. This approach guarantees absence of *k*-covering traps, i.e. local minima of the deformation energy where the mesh is free of inverted elements, but total angle around a certain vertex is above $2\pi$ for 2D and above $4\pi$ for 3D. We demonstrate enhanced stability of suggested untangling technique which has a potential to make the mesh untangling a routine operation.

While only simplicial meshes are considered, generalization of presented algorithms to non-simplicial elements (quads, hexes, B-spline elements) is straightforward.

The rest of the paper is organized as follows: first, in § 1.2.1 we recall the untangling method proposed in [GKK$^{+21}$], then in § 1.2.2 we discuss the structure of the admissible set, and we show the limitations of the method for free boundary deformations. Next in § 1.3 we propose a way to alleviate the issue: we introduce phantom element technique for prevention of 2-covering traps. This technique guarantees that a deformation free of inversions is locally invertible.
Finally, we test our methods in § 1.6.

## 1.2 Foldover-free deformations

In this section we go over the method proposed in the paper [GKK+21]. First, in § 1.2.1 we recapitulate the main idea and the resolution scheme, and then in § 1.2.2 we show the main limitation: even without inverted elements, the deformation can be non-invertible.

### 1.2.1 Variational formulation for grid generation

Let us go through the algorithm to compute a foldover-free deformation $\vec{x}(\xi) : \Omega \subset \mathbb{R}^d \to \mathbb{R}^d$. This presentation is unified both for 2D and 3D settings, and by $d$ we denote the number of dimensions; in our notations we use arrows for all vectors of dimension $d$. Consider the following variational problem:

$$ \arg\min_{\vec{x}(\xi)} \int_{\Omega} w(\xi) f(J) d\xi, \qquad (1.1) $$

where $J$ is the Jacobian matrix of the mapping $\vec{x}(\xi)$, $1 \ge w(\xi) \ge w_{\max}$ is the weight function and density of deformation energy is defined by [Gar00]

$$ f(J) = (1 - \theta)f_s(J) + \theta f_v(J), \qquad (1.2) $$

where shape distortion is defined as

$$ f_s(J) := \begin{cases} \frac{1}{d} \text{tr} J^\top J, & \det J > 0 \\ \frac{(\det J)^{\frac{2}{d}}}{d}, & \det J \le 0 \end{cases} \qquad (1.3) $$

while volumetric distortion is defined

$$ f_v(J) := \begin{cases} \frac{1}{2} \left( \det J + \frac{1}{\det J} \right), & \det J > 0 \\ +\infty, & \det J \le 0 \end{cases} \qquad (1.4) $$

The weight function is introduced for better control of mesh quality for particular subsets of mesh cells [Gar00] which are known apriori, say for cells in the boundary layer.

Prob. (1.1) may be subject to some constraints that we do not write explicitly. To give an example, one may constrain position of some vertices. In this formulation,
functions $f_s(J)$ and $f_v(J)$ have concurrent goals, one preserves angles and the other preserves the area, and thus $\theta$ serves as a trade-off parameter.

In finite elasticity $f$ (a sum of shape distortion and volume distortion terms) is called “isochoric-volumetric split”. The idea of such a split goes back to 70s [Flo61], [Pen70], [DBVDBZ88] and remains a hot topic in convexity analysis [VGMN21].

Note that function $f$ is not convex, but polyconvex. The notion of polyconvexity is a generalization of the notion of convexity for functions defined on spaces of matrices. A function $\phi(J) : \mathbb{R}^{d \times d} \to \mathbb{R} \cup +\infty$ is said to be polyconvex [Bal76] if there exists a convex function $\Phi(\#J)$, such that $\phi(J) = \Phi(\#J)$, where $\#J$ denotes the set of all minors of $J$.

Any polyconvex function $\phi$ is rank-one convex [Bal76], i.e.

$$ \phi((1-\theta)J + \theta(J + \delta J) \leq (1-\theta)\phi(J) + \theta\phi(J + \delta J), $$

where rank $\delta J = 1$, and thus satisfies the Hadamard-Legendre conditions (ellipticity conditions for the Euler-Lagrange equation of variational Prob.(1.1))

$$ \sum_{i,k,j,m=1}^{d} \frac{\partial^2 \phi}{\partial(J)_{ij} \partial(J)_{km}} p_i p_k q_j q_m \geq 0 $$

for arbitrary vectors $\vec{p}, \vec{q} \in \mathbb{R}^d$.

Since density of energy (1.2) is a polyconvex function satisfying the ellipticity conditions, it is therefore very well suited for a numerical optimization provided that we have an initial guess in the admissible domain $\min_{\Omega} \det J\{\vec{x}\}(\vec{\xi}) > 0$.

Note that if an initial guess is not admissible, then functional (1.1) is not defined. In this case a regularized version can be introduced [GK99b]: we can avoid non-positive denominators in $f_s$ and $f_v$ using a regularization function $\chi$ for a positive value of $\varepsilon$:

$$ \chi(D, \varepsilon) := \frac{D + \sqrt{\varepsilon^2 + D^2}}{2} \qquad (1.5) $$

Then one can define a regularized version $f_\varepsilon$ of functions $f$:

$$ f_{\varepsilon}(J) := (1 - \theta) \frac{\frac{1}{d} \mathrm{tr} J^{\top} J}{(\chi(\det J, \varepsilon))^{\frac{2}{d}}} + \theta \frac{1}{2} \frac{1 + \det^2 J}{\chi(\det J, \varepsilon)}, \qquad (1.6) $$

so that Prob. (1.1) is reformulated as

$$ \lim_{\varepsilon \to 0^+} \arg \min_{\vec{x}(\vec{\xi}) \in \Omega} \int w(\xi) f_{\varepsilon}(J) d\xi \qquad (1.7) $$

The Prob. (1.7) does offer a way of getting rid of foldovers if a foldover-free initialization is not available.

In practice, the map $\vec{x}$ is piecewise affine with the Jacobian matrix $J$ being piecewise constant, and can be represented by the coordinates of the vertices in the computational domain $\{\vec{x}_i\}_{i=1}^V$. Let us denote the vector of all variables as
$X := (\vec{x}_1^\top \dots \vec{x}_{\#V}^\top)^\top$, then optimization Prob. (1.1) has following expression in discrete settings:

$$ \lim_{\epsilon \to 0^+} \arg \min_{X} F(X, \epsilon), \qquad (1.8) $$

$$ \text{where} \quad F(X, \epsilon) := \sum_{k=1}^{\#T} w_k f_{\epsilon}(J_k) \text{vol}(T_k), $$

\#V is the number of vertices, \#T is the number of simplices, $J_k$ is the Jacobian matrix for the k-th simplex and vol($T_k$) is the signed volume of the simplex $T_k$ in the parametric domain. Note that we use term "parametric mesh" for a manifold glued from target elements. For surface flattening problem the parametric domain coincides with surface triangulation. For simplicial elements it is convenient to assign constant weight $w_k$ to each mesh cell.

We denote by $P_k$ the k-th simplex of the computational mesh. Since it is an affine image of $T_k$

$$ P_k = \vec{x}(T_k), $$

we get

$$ \det J_k = \text{vol} P_k / \text{vol} T_k. $$

We assume that vol $T_k > 0$ for all parametric tetrahedra.

To solve Prob. (1.8), we use an iterative descent method, as suggested in [GKK$^{+}$21]. Starting from an initial guess $X^0$, they build a sequence of approximations $X^{k+1} := X^k + \Delta X^k$, carefully choosing the regularization parameter $\epsilon^k$ for each iteration $k$.

### 1.2.2 Discussion: set of admissible deformations

Function $F(X, 0)$ has impenetrable infinite barrier on the boundary of the set of meshes with positive cell volumes

$$ \frac{\text{vol}(P_k)}{\text{vol}(T_k)} > 0, \quad k = 1, \dots, \#T \qquad (1.9) $$

which is finite-dimensional approximation of the set

$$ \det J > 0. \qquad (1.10) $$

This set has a quite complicated structure. For k-th simplex vol($P_k$) is a polylinear function of coordinates of its vertices, hence each term in (1.9) defines a non-convex set. One can hardly expect that intersection of the sets in (1.9) would result in a convex domain. Moreover, Ciarlet [CG82] has proved that barrier property and convexity of the density of deformation energy are incompatible. From his statement it essentially follows that when finite element approximation $S(X)$ of the deformation
energy is bounded in $\mathcal{A}$ and

$$S(X) \to +\infty \text{ when } X \in \mathcal{A}, X \to \partial \mathcal{A},$$

then it cannot be convex function of X. Fortunately, barrier distortion measures can be polyconvex, as shown by J. Ball [Bal76]. For instance, it means that domain $\mathcal{A}$ consists of connected components where each two points can be connected. Consider vector $X+sY \in \mathbb{R}^{\#Vd}$, $0 \le s \le 1$, $X \in \mathcal{A}$, $X+Y \in \mathcal{A}$, and assume that for k-th mesh cell the Jacobian matrix for mesh defined by $X+sY$ is written as

$$J_k + sB_k, \text{ rank } B_k = 1. \qquad (1.11)$$

Since $-\det J$ is rank one convex function of $J$, we get

$$-\det(J_k + sB_k) \le -(1-s)\det J_k - s\det(J_k + B_k),$$

meaning that such a deformation for any $s$ remains inside admissible set. Below we will consider subsets of $\mathcal{A}$ where any pair of points be connected using sequence of rank-one segments.

The admissible set $\mathcal{A}$ may have a quite complicated structure. Even in the case of fixed boundary vertices, it may contain disjoint subsets. An example of disjoint sets is shown on Figure 1.1. For a square with square hole we show 3 admissible meshes with the same boundary conditions and the same connectivity but with different winding numbers. It is not possible to deform one mesh onto another with fixed boundary vertices.

Figure 1.1(d) shows three components of the barrier function and a global penalty function. Secondary subsets create kind of a “rabbit holes” with small attraction domains meaning that if we do not need “principal subset” with large attraction domain, one would need to use global optimization algorithms in order to get into rabbit holes.

For the case of mesh deformations with free boundaries the structure of $\mathcal{A}$ becomes even more complex. The number of disjoint sets may sharply increase. Moreover, parasitic disjoint sets with small energy may appear which present real traps for untangling algorithms. Consider flattening of the surface vertex star consisting of 6 regular triangles, see Figure 1.2(a). While standard flattening would be regular 12-gon, shown in Figure 1.2(b), minimal deformation energy is provided by “2-covering”, shown in Fig. Figure 1.2(c).

Homeomorphism between these two solutions does not exist, so they definitely belong to disjoint subsets of $\mathcal{A}$. The second solution defines flattening which is not locally invertible since one can not find small open neighborhood of the vertex star which can be one-to-one mapped on the open neighborhood of the flat projection.
<figure><img src="image_3.png" /><figcaption>Fig. 1.1: Admissible set with disjoint components: (a) mesh with minimal distortion, (b), (c) untangled meshes with the same boundary conditions but with opposite winding numbers, (d) illustrative barrier mesh functional (red graph) with disjoint subsets and “rabbit holes” matching these subsets on the black graph of penalized mesh functional.</figcaption></figure>

**1.3 Make foldover-free deformations invertible**

In order to explain the idea of the mesh protection algorithm which enforces local invertibility for 2D and 3D mesh deformations, consider simple 2D vertex stars shown in Figure 1.3. We denote by *k*-star the vertex star with triangle valence equal to *k*.

3-stars and 4-stars do not admit 2-coverings or *k*-coverings with *k* > 1 since total angle around vertex larger or equal to 4π would require triangle angles above π, making mesh inadmissible. We may add to the star phantom triangles spanning its central vertex and a pair of outer vertices of the star. For 5-star we may add to the star
<figure><img src="image_4.png" /><figcaption>Fig. 1.2: 2-covering trap free-boundary flattening. (a): the surface to flatten is made of 12 equilateral triangles. (b): an invertible flattening corresponding to a local minimum of the elastic energy. (c): global minimum of the elastic energy. This flattening is free of inverted elements, however is not invertible in the vicinity of the center vertex.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>Fig. 1.3: The concept of protection of 2D vertex stars via additional phantom triangles: 3- and 4- stars does not admit 2-covering. For 5-star one additional (overlapping) triangle effectively prevents double covering, since it basically reduces the configuration to a 4-star. For 6-star 2 triangles suffice to prevent any coverings.</figcaption></figure>
single overlapping triangle, while for 6-star two additional overlapping triangles are enough. Hence the number of unknowns for variational problem is fixed, while the definition (1.9) of admissible set $\mathcal{A}$ is augmented by additional inequalities which cut off disjoint subsets related to this mesh vertex.

In 2D we start by creating flat vertex stars in parametric domain. We presume that all triangles in parametric domain are known. Otherwise it is natural to assume that all parametric cells are unit triangles. We can build flat stars by simple angle-based flattening.

Then recursive aggregation algorithm is applied, which creates the sequence of vertex stars independently around each vertex by aggregating adjacent triangles until basic 3- or 4-stars are created. On each step we identify the adjacent pair of triangles forming quadrilateral (in general non-convex one) which provides the best quality for covering by new phantom triangle $T_1^{(1)}$ as shown in Figure 1.4(a), (b).

<figure><img src="image_6.png" /><figcaption>Fig. 1.4: Recursive aggregation of adjacent triangles to create phantom protective triangles. (a) Initial 8-stencil, (b) 1 level blue phantom triangle is added, (c) three 2 level phantom triangles are added resulting in a 4-star.</figcaption></figure>

Several non-overlapping pairs form 1st level of phantom triangles and reduce the valence of the star. Then we repeat the aggregation getting the second level by creating phantom triangles $T_i^{(2)}$, $i = 1, 3$, see Figure 1.4(c) and continue this procedure until basic star is obtained.

Figure 1.5 explains how aggregation algorithm treats regular boundary vertices. The case when boundary star is an approximate half of the full internal star is considered as a regular case. We apply aggregation algorithm until basic boundary 2-star is created. In particular case shown in Figure 1.5, one first level triangle $T_1^{(1)}$ and two second level triangles $T_i^{(2)}$, $i = 1, 2$ are added.

The irregular boundary case is shown in Figure 1.6 (left). Here aggregation could not create topologically correct boundary star without inverted triangles. The solution is simple: to add single outer phantom triangle thus closing the star and making it an internal one. After that standard aggregation algorithm can be applied.

In the 3D case there is no need to perform cumbersome analysis of vertex stars. Consider the set of 3-sided polyhedral cones originating from the inner mesh vertex $p$. Faces of the cone correspond to three faces of the tetrahedron adjacent to $p$.
<figure><img src="image_7.png" /><figcaption>Fig. 1.5: Regular boundary star: recursive adding of phantom triangles results in a 2-star.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>Fig. 1.6: Irregular boundary star: configuration of the boundary corresponds to sharp incoming corner. Single outer phantom triangle transforms boundary star into internal star.</figcaption></figure>

Intersection of i-th cone with unit sphere defines the spherical triangle $T_i$ with angles $\alpha_i, \beta_i, \gamma_i$. These angles coincide with three dihedral angles of the cone/tet. The number of spherical triangles $n_t$ is equal to the number of tetrahedra adjacent to $p$, while the number $n_v$ of vertices of spherical triangulation is equal the number of mesh edges originating from $p$. Obvious relation $n_t = 2n_v - 4$ holds. Since all tets has positive algebraic volume, the area of each spherical triangle is positive

$$ \text{area}(T_i) = \alpha_i + \beta_i + \gamma_i - \pi > 0. $$

Then

$$ \sum_{i=1}^{n_t} \text{area}(T_i) = \sum_{k=1}^{n_v} \Delta \alpha_k - n_t \pi $$

where $\Delta \alpha_k$ is the sum of angles around k-th vertex. If we impose the condition

$$ \Delta \alpha_k = 2\pi, $$

which means that no 2-coverings are allowed for tetrahedral mesh edges, we get
$$ \sum_{i=1}^{n_t} \text{area}(T_i) = 2\pi n_v - n_t \pi = 4\pi, $$

meaning that 2-covering for vertex star is forbidden.

Hence we consider quasi-2D algorithm which prohibits 2-coverings for tetrahedral fans around each edge of the tet mesh. Figure 1.7 shows sample 2-covering for tetrahedra around an edge.

<figure><img src="image_9.png" /><figcaption>Fig. 1.7: Left: a regular fan of 6 tetrahedra sharing the common edge MN. Right: a double covering around the edge MN.</figcaption></figure>

Figure 1.8 illustrates aggregation algorithm in 3D. The set of tetrahedra around an edge defines 2D triangle arrangement around vertex in the plane orthogonal to the edge. Logics of 2D aggregation algorithm can be used to create phantom tets around the edge keeping the edge fan geometrically and topologically correct. The only difference is the 3D quality criteria for creation of best tets on each level should be used.

For boundary vertices we may encounter regular and irregular stars. For regular star (see Figure 1.9(left)) there exists internal edge $e_{int}$ which is not far from orthogonal to the outer boundary of the star. Then one can apply quasi-2D aggregation around all boundary edges and create new edge star using $e_{int}$ and boundary triangles $T_i$ creating tets as $\text{conv}(e_{int}, T_i)$. Edge-based aggregation should be applied to this edge as well. As a result one can guarantee that total spherical angle around vertex is smaller then $4\pi$. For irregular case shown in Figure 1.9(right) good internal edge for boundary vertex $p$ does not exist hence one should add new phantom vertex $p^*$ outside thus creating new edge and closing the boundary star by adding all tets wither the vertices $p, p^*, p_i, p_{i+1}$ for all boundary vertices $p_i$ of the boundary star. As a result star of the vertex $p$ can be treated as an internal star.
<figure><img src="image_10.png" /><figcaption>Fig. 1.8: Adding phantom tetrahedra around an edge is equivalent to a 2D problem.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>Fig. 1.9: (left) For regular boundary star in 3D standard edge-based aggregation algorithm is applied, (right) for irregular boundary vertex star outer vertex is added creating internal star configuration.</figcaption></figure>

In the “irregular” case the vertex $p$ becomes the internal one and total spherical angle analysis can be applied directly guaranteeing absence of 2-covering. In the regular case on the outer boundary fragment adjacent to $p$ each face satisfies the visibility condition from any point lying on the internal edge $e_{int}$, and total dihedral angle around this edge is equal to $2\pi$. Hence edge $e_{int}$ admits straight prolongation outside the computation domain in such a way that one can complement tetrahedral half-star with outer complementary half-star creating full tet star satisfying dihedral edge constraint for all internal edges.

In principle one can add prismatic layer of phantom cells around all free boundary faces. It would guarantee absence of 2-covering, but the preprocessing becomes quite expensive and considerably increases the number of phantom tets. The main difference between 2D and 3D protection procedure is that in 3D the size of the variational problem increases.
<figure><img src="image_12.png" /><figcaption>Fig. 1.10: (left) Correct phantom triangle for star glued from target triangles, (right) orientation of the same phantom triangle is forced to be wrong due to presence of fixed vertices (dark blue circles) in the computational domain.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>Fig. 1.11: (left) Admissible phantom triangle in the computational domain may force bad quality triangles, (right) conservative solution: split constrained star into sectors and create phantom triangles independently inside each sector prohibiting creation of separating edges.</figcaption></figure>

Generally recursive aggregation serves to reduce the admissible set $\mathcal{A}$ by eliminating its parasitic components. However potentially it may not just make an untangling problem more stiff, but also over-reduce admissible set making untangling problem unresolvable. While it is very hard to make constructive statements about the structure of the admissible set, our claim supported by numerical tests is that for free boundary problems hierarchical aggregation procedure does not create unsolvable untangling problem. Moreover, since all phantom elements can be eliminated after untangling procedure, the quality of resulting mesh is not affected. As a side note, using the set of overlapping simplicial elements for the vertex set is very natural for
variational mesh optimization since non-simplicial elements (quad elements being the simplest example) are modelled in such a way.

However when constrained problem which combines free boundaries with fixed vertices is considered, aggregation potentially may result in empty admissible set. Consider 2d vertex star glued from target triangles shown in Figure 1.10 (left). New phantom triangle has correct orientation. In the presence of the fixed vertices, shown in dark blue, the boundary of the same star in the computational domain has concave fragment which forces inverse orientation of the phantom triangle and empty admissible set A.

Note that the main source of problems is the creation of separating edges, i.e. the edges with two fixed vertices, which actually split our computational domain by creating cuts with prescribed boundary conditions. Figure 1.11 (left) shows that even in the case when phantom triangle is not forced to be misoriented, its separating edge may cut off fragment of the computational domain and make untangling problem too stiff.

We suggest simple conservative aggregation rule in the presence of fixed vertices. Note that fixed vertices shown in dark blue essentially split 2d vertex star in sectors as shown in Figure 1.11 (right).

One have to treat each sector independently, applying aggregation until two triangles cover the full sector. This algorithm does not create separating edges and avoids locking phenomenon shown above. We cannot guarantee that particular aggregation algorithm would always create good admissible set, however one can easily build a full set of phantom triangles for any existing admissible mesh.

## 1.4 Finite untagling theorem and its limitations

Let us restate the main result of [GKK$^{+}$21, Theorem 1], which shows that untangled mesh can be attained in a finite number of approximate unconditional minimizations step.

**Theorem 1.** Let us suppose that the feasible set of untangled meshes (1.9) is not empty. We also suppose that for solving $X^{k+1} \leftarrow \arg\min_{X} F(X, \varepsilon^k)$ we have a minimization algorithm satisfying the following efficiency conditions for some $0 < \sigma < 1$: For each iteration k,

* either the essential descent condition holds

$$F(X^{k+1}, \varepsilon^k) \le (1 - \sigma)F(X^k, \varepsilon^k), \quad (1.12)$$

* or the vector $X^k$ satisfies the quasi-minimality condition:

$$\min_{X} F(X, \varepsilon^k) > (1 - \sigma)F(X^k, \varepsilon^k). \quad (1.13)$$
*Then the feasible set (1.9) is reachable by solving a finite number of minimization problems in X with $\varepsilon^k$ fixed for each problem.*

This theorem provides rigorous foundations for untangling process. However, efficiency condition essentially implies ability of the black box unconditional minimization tool to solve global optimization problem for nonconvex function. This assumption is not practical. The most unpleasant phenomenon is not conconvexity itself, but rather appearance of local minima of $F(X)$ outside feasible set. Numerical experience suggests that for fixed boundary conditions such a local minima even if exist, do not influence the trajectory of solution during untangling. While in the presence of free boundaries or when closed manifold unfolding problem is considered, appearance of the locked mesh state with 2-covering leads to strong local minimum.

**1.5 Practical algorithm for elimination of topological singularities and corresponding local minima of penalty functional**

Adding phantom elements to stars of all mesh vertices in 2d and to stars of all mesh edges in 3d evidently guarantees that meshes with topological singularities such as 2-covering and $k$-covering are not admissible, hence penalty functional would eventually attain arbitrary large values on these meshes with $\varepsilon \to +0$. The resulting discrete problem should process about 2-3 times larger number of cells and generally is stiffer compared to original problem. Unfortunately it does not guarantee successfull untangling when 2-covering is present in the initial mesh.

Empirical evidence suggests that in this case penalty functional has parasitic local minimum outside admissible set and with parameter $\varepsilon$ getting smaller this minimum instead of collapse gets deeper and sharper and moves closer to the boundary of the admissible set from outside. Speculative pattern of the functional is shown in Fig. 1.12(b)-(d). The practical consequence is that mesh is stuck in the position when 1-2 relatively large mesh elements with extremely small signed area/volume are created near topological singularity. Untangling fails! The finite untangling theorem is no longer applicable since efficiency condition for minimization tool-box is violated.

We solve this problem using special adaptation procedure which adds weighted phantom elements. The idea is that we add phantom elements only to mesh vertices with topological singularities and assign increasing sequence of weights to phantom cells: $w_k = 1$ for initial mesh cells, and value $w^j$ is assigned to $j$-th level of phantom cells, where $w > 1$ is a constant. Hence the coarsest phantom cells have largest weights. In 2d we suggest to use $w = 10$, while for 3d the weight $w = 100$ is assigned. After adding phantoms mesh is moved into new optimal position. Observed behaviour is that topological singularity is moved to the nearby vertex. We add weighted phantom elements to the new vertex, while keeping added phantoms. We have found that this simple procedure to be quite efficient since singularities quickly move to the boundary and then disappear.
<figure><img src="image_14.png" /><figcaption>Fig. 1.12: Local minimum for penalty functional $F(X, \varepsilon)$. (a) discrete barrier func-
tional $F(X, 0)$, (b) intermediate penalty function with local minimum $F(X, \varepsilon_1)$, (c)
intermediate penalty function with sharper local minimum $F(X, \varepsilon_2)$, $\varepsilon_2 < \varepsilon_1$, (d) final
penalty function for minimal attainable value $\varepsilon^*$. (e) red line: adaptive assignment of
non-uniform phantom cell weights moves and lowers the local minimum, (f) adaptive
assignment of non-uniform phantom cell weights eventually eliminates parasitic
local minimum.</figcaption></figure>

Fig. 1.13 show setting for basic 2d test: surface mesh is glued from 12 equilateral triangles. 10 phantom triangles are added to counter singularity.

Fig. 1.14 (left) shows initial flat mesh with double covering, while deadlocked
solution which corresponds to unity weights and local minimum of penalty func-
tional is shown in the center. Adaptive weight distribution allows to eliminate local
minimum and attain correct mesh.

Fig. 1.15 illustrates the same algorithm in 3d. Initial mesh contains 12 tetrahedra
around single edges which defined total dihedral angle equal to 4π. 6 second order
tetrahedra with weight equal to 100 and 2 third order tets with weight equal to 10⁴
are added. We use tet shrinking to help understanding of the degenerate initial mesh.
Using unity weights leads to deadlocked solution, while adaptive weights allow to
untangle mesh.
<figure><img src="image_15.png" /><figcaption>Fig. 1.13: Surface mesh consisting of 12 triangle with 10 phantom triangles.</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>Fig. 1.14: Left: Initial flattened mesh with double covering and inverted phantom triangles, center: deadlocked solution: local minimum of penalty functional is attained outside admissible set, right: assignment of largest weights to coarsest phantom triangles eliminates local minimum, correct mesh is attained.</figcaption></figure>

Figures 1.17 and 1.18 illustrated the behaviour of algorithm for annihilation of the singularity. Texture mapping is used to help understand

We can interpret the gradients of the parameterization as two vector fields. Double coverings arise from index -1 singularities of this field (refer to Fig. 1.16). Poincaré-Hopf theorem states that the sum of indices is constant, so an index 1 singularity must be placed somewhere to compensate for the index -1. With free boundary mapping (as in our example in Fig. 1.16), the index 1 is placed outside of the domain, simply adding a loop to the boundary. For the global parameterization case, there is no free boundary, so the index 1 must appear at a vertex. Recall however, that our maps have positive Jacobian determinant over all elements, and thus it is impossible to place singularity 1 at a regular vertex, since it is a pole singularity that would force the map to have degenerate elements.

Fig. 1.17 illustrates simple free-boundary case: regular hexagon is covered by regular triangles. This ideal mesh is tangled by creation of 2-covering: single internal vertex has total angle equal to $4\pi$ (-1 singularity) which is compensated by total
<figure><img src="image_17.png" /><figcaption>Fig. 1.15: Upper row: initial guess includes 12 initial tets plus 6 second order phantom ones and 2 third order phantoms tets. Some of the phantom tets geometrically coincide. Two negative tets are present. Shrinking is used to help understanding the layout. Lower row: (left) - deadlocked degenerate solution, (center-right) good solution for weighted functional, shrinking shows 12 good tets.</figcaption></figure>

<figure><img src="image_18.png" /><figcaption>Fig. 1.16: A double covering example. **Left:** a surface to flatten is made of 12 equilateral triangles. **Middle:** the surface can be mapped to a plane without inverted elements, producing a double covering. The map is not invertible in one point. **Right:** gradients of the parameterization form two orthogonal vector fields (shown in red and blue), both of them have a -1 singularity at the center.</figcaption></figure>

rotation angle $4\pi$ for exterior boundary of the domain (+1 singularity). We show initial mesh and few intermediate meshes as well as mark all added phantom triangles on the ideal parametric mesh. One can clearly see that vertex singularity moves to the boundary by almost shortest trajectory and eventually annihilates via contact with boundary singularity.

Fig. 1.18 shows more elaborated case of 2-covering for fixed boundaries. Note that in practice it is really hard to encounter such a situation: 2-coverings in this case
<figure><img src="image_19.png" /><figcaption>Fig. 1.17: Recovering local injectivity for foldover-free deformations. Here the boundary of the domain is free to move. Upper row: $\xi$ domain, bottom row: $\bar{x}$ domain. Left: initial deformation with two-covering present in the center. Middle and right: we progressively force the two covering outside of the domain by adding virtual triangles (in green) on the problematic 1-rings.</figcaption></figure>

are not stable and are not observed in practice. In order to create this test case we have to fix certain internal vertices in a positions which make 2-covering stable thus creating parasitic local minimum for penalty functional. Note that -1 singularity (vertex with total angle equal to $4\pi$) should be compensated by the +1 singularity (mesh vertex with total angle equal to zero). Figure shows how adaptive weighting allows to merge these singularities.

The speculative explanation of the observed behaviour of suggested algorithm is presented in Fig. 1.12(e)-(f). It seems that adaptive addition of phantom cells makes parasitic local minimum move around boundary of admissible set lowering down in the process and eventually leading to penalty function without parasitic local minima which allows to recover efficiency property for the unconstrained minimization toolbox.

## 1.6 Results and discussion

We start numerical experiments with a simple sanity check: two arbitrary vertices are swapped, and the resulting mesh is untangled. Only two vertices are locked, the rest of the mesh is free to move. Figure 1.19(a) provides the rest shape, the vertices to be swapped are highlighted in red and blue, respectively. Figure 1.19(b) shows that untangling may result in a mesh with double coverings. Untangling of a protected mesh produces the correct result, namely, a rigid transformation of the input mesh (Figure 1.19(c)).
1 Mesh Untangling for Problems with Topological Singularities

<figure><img src="image_20.png" /><figcaption>Fig. 1.18: Recovering local injectivity for deformations with fixed boundary. Note that under properly constrained boundary, it is impossible to obtain a 2-covering (-1 singularity) without inverted elements, since we have to compensate it with a +1 singularity. Left: initial deformation with two-coverings and some inverted elements present. Red dots show the locked vertices. By adding virtual triangles (in green) to the problematic 1-rings, we can merge -1 and +1 singularities, thus recovering local injectivity.</figcaption></figure>

<figure><img src="image_21.png" /><figcaption>Fig. 1.19: Free boundary deformation sanity check: exchanging two vertices (shown in red and blue) in the mesh and then untangling. Only two vertices are locked, the rest of the mesh is free to move. (a): input mesh, (b): foldover-free but not invertible deformation, (c): invertible deformation.</figcaption></figure>

Figure 1.20 shows a little bit more difficult sanity check. The idea is to undo a metal sheet forming, i.e., untangle a mechanical piece model which is almost flat intrinsically. In this test all vertices are free to move, we perform an orthogonal projection of the model on the *Oxy* plane and then untangle the mesh. If the untangling
<figure><img src="image_22.png" /><figcaption>Fig. 1.20: Inverse metal sheet forming. (a): input model, (b): unprotected untangled mesh contains 2-covering traps, (c): protected mesh defines locally invertible flattening.</figcaption></figure>

<figure><img src="image_23.png" /><figcaption>Fig. 1.21: Flattenings of a very rough surface. (a): input surface to flatten, (b): unprotected untangled mesh contains 2-covering traps, (c): protected mesh defines locally invertible flattening.</figcaption></figure>

is executed without adding phantom triangles, the mesh gets caught into 2-covering traps (Figure 1.20(b)). Finally, Figure 1.20(c) shows the result of untangling of the protected mesh. It is easy to see that the inverse metal sheet forming has succeeded.

Next we have performed basically the same test, but on a much more challenging model. Figure 1.21(a) shows the input mesh to flatten, the mesh is highly curved and presents very bad quality elements. Refer to Figures 1.21(a–b) for the untangling with and without protecting elements.
# 1 Mesh Untangling for Problems with Topological Singularities

Then we show that the same problem exists in 3D as well. Figure 1.22 is a stress test for the free boundary untangling, similar to rotated cube-inside-cube stress test considered in [GKU14, GKK$^{+21}$]. We have created an isotropic tetrahedral mesh of a cuboid with two cubic cavities inside. Then we have rotated the boundary of the cavities by 135 degrees around the vertical axis, thus producing inverted tetrahedra. We have constrained the boundary of the cavities, and left the rest of the mesh free to move. As demonstrated in Figure 1.22, without protection one can encounter 2-covering traps in 3D as well, while protection allows to obtain locally invertible deformations. In this particular case the protected deformation is globally invertible.

<figure><img src="image_24.png" /><figcaption>Fig. 1.22: Upper row: two cubic cavities inside elastic cuboid are rotated by 135 degrees. Outer boundary is free. Middle row: untangled mesh with 2-covering stars. Lower row: protected mesh is globally invertible.</figcaption></figure>

Finally we show that it is possible to go beyond simple deformations. Figures 1.23 and 1.24 show a quad-remeshing application. The idea is to define a deformation of the input surface such that if the final quad mesh (the result) undergoes this deformation, it matches a unit, axis aligned grid. The direct application of this idea computes this deformation, applies its inverse to the unit grid, and obtains a quad mesh. In practice, it is better to introduce more degrees of freedom by considering
global parameterizations instead of a deformation. In this case, parameterizations
have some discontinuities that make it possible to represent a much larger family of
quad meshes: the deformed grid can be cut and glued to itself in a non-trivial way.

Running QuadCover [KNP07] to do so often results in a local loss of injectiv-
ity and, as illustrated in Figure 1.23–left, one may need to call a backup solution
such as QEx [EBCK13] to extract a quad mesh despite the tangled parameterization.
Free-boundary variational smoother overcomes this obstacle, producing valid pa-
rameterizations (Figure 1.23–right). We need, however, protect the mesh to avoid
k-coverings, as illustrated in Figure 1.24. Note that in variational method the quad
mesh is represented as a set of overlapping triangles, since in order to guarantee
convexity of quad cells in barrier method we split each quad into four overlapping
triangles corresponding to its corner [Iva88].

**1.7 Conclusion**

We formulate a set of variational problems potentially covering the complete tech-
nological chain for construction of optimal mappings and deformations with fixed
as well as free boundaries. We start with the continuation problem with respect to
parameter ε, this minimization allows us to compute optimal in the average defor-
mations. We illustrate performance of our algorithm with challenging 2D and 3D
numerical tests.

**References**

[Bal76] John M Ball. Convexity conditions and existence theorems in nonlinear elasticity. *Archive for rational mechanics and Analysis*, 63(4):337–403, 1976.

[Bal81] J. M. Ball. Global invertibility of sobolev functions and the interpenetration of matter. *Proceedings of the Royal Society of Edinburgh: Section A Mathematics*, 88(3-4):315–328, 1981.

[CG82] P.G. Ciarlet and G. Geymonat. Sur les lois de comportement en elasticite non-lineaire compressible. *C.R. Acad.Sci. Paris Ser.II*, 295:423 – 426, 1982.

[CI97] AA Charakhch'yan and SA Ivanenko. A variational form of the winslow grid generator. *Journal of Computational Physics*, 136(2):385–398, 1997.

[Cro62] WP Crowley. An equipotential zoner on a quadrilateral mesh. *Memo, Lawrence Livermore National Lab*, 5, 1962.

[DAZ<sup>+</sup>20] Xingyi Du, Noam Aigerman, Qingnan Zhou, Shahar Z. Kovalsky, Yajie Yan, Danny M. Kaufman, and Tao Ju. Lifting simplices to find injectivity. *ACM Trans. Graph.*, 39(4), July 2020.

[DBVDBZ88] R De Borst, PAJ Van Den Bogert, and J Zeilmaker. Modelling and analysis of rubberlike materials. *HERON*, 33 (1), 1988, 1988.

[DS13] Josh Danczyk and Krishnan Suresh. Finite element analysis over tangled simplicial meshes: Theory and implementation. *Finite Elements in Analysis and Design*, 70-71:57–67, 2013.

[EBCK13] Hans-Christian Ebke, David Bommes, Marcel Campen, and Leif Kobbelt. Qex: Robust quad mesh extraction. *ACM Trans. Graph.*, 32(6), November 2013.
1 Mesh Untangling for Problems with Topological Singularities

[ERM<sup>+</sup> 03] José María Escobar, Eduardo Rodríguez, Rafael Montenegro, Gustavo Montero, and José María González-Yuste. Simultaneous untangling and smoothing of tetrahedral meshes. *Computer Methods in Applied Mechanics and Engineering*, 192(25):2775-2787, 2003.

[Flo61] P. J. Flory. Thermodynamic relations for high elastic materials. *Trans. Faraday Soc.*, 57:829-838, 1961.

[FP00] Lori A Freitag and Paul Plassmann. Local optimization-based simplicial mesh untangling and improvement. *International Journal for Numerical Methods in Engineering*, 49(1-2):109-125, 2000.

[Gar00] VA Garanzha. The barrier method for constructing quasi-isometric grids. *Computational Mathematics and Mathematical Physics*, 40:1617-1637, 2000.

[GK99a] VA Garanzha and IE Kaporin. Regularization of the barrier variational method. *Computational mathematics and mathematical physics*, 39(9):1426-1440, 1999.

[GK99b] V.A. Garanzha and I.E. Kaporin. Regularization of the barrier variational method of grid generation. *Comput. Math. Math. Phys.*, 39(9):1426-1440, 1999.

[GKK<sup>+</sup> 21] Vladimir Garanzha, Igor Kaporin, Liudmila Kudryavtseva, François Protais, Nicolas Ray, and Dmitry Sokolov. Foldover-free maps in 50 lines of code. *ACM Transactions on Graphics*, 40(4), 2021.

[GKU14] V.A. Garanzha, L.N. Kudryavtseva, and S.V. Utyuzhnikov. Variational method for untangling and optimization of spatial meshes. *Journal of Computational and Applied Mathematics*, 269:24-41, 2014.

[GP72] S.K. Godunov and G.P. Prokopov. The use of moving meshes in gas-dynamical computations. *USSR Computational Mathematics and Mathematical Physics*, 12(2):182-195, 1972.

[Iva88] SA Ivanenko. Construction of nondegenerate grids. *Zh. Vychisl. Mat. Mat. Fiz*, 28(10):1498, 1988.

[Jac88] Olivier-P Jacquotte. A mechanical model for a new grid generation method in computational fluid dynamics. *Computer methods in applied mechanics and engineering*, 66(3):323-338, 1988.

[KNP07] Felix Kälberer, Matthias Nieser, and Konrad Polthier. Quadcover - surface parameterization using branched coverings. *Computer Graphics Forum*, 26(3):375-384, 2007.

[Knu00] Patrick M Knupp. Achieving finite element mesh quality via optimization of the jacobian matrix norm and associated quantities. part ii—a framework for volume mesh optimization and the condition number of the jacobian matrix. *International Journal for numerical methods in engineering*, 48(8):1165-1185, 2000.

[Knu01] Patrick M Knupp. Hexahedral and tetrahedral mesh untangling. *Engineering with Computers*, 17(3):261-268, 2001.

[LJ94] Anwei Liu and Barry Joe. On the shape of tetrahedra from bisection. *Mathematics of computation*, 63(207):141-154, 1994.

[Pen70] Robert W. Penn. Volume changes accompanying the extension of rubber. *Transactions of the Society of Rheology*, 14(4):509-517, 1970.

[Res66] Yu G Reshetnyak. Bounds on moduli of continuity for certain mappings. *Siberian Mathematical Journal*, 7(5):879-886, 1966.

[Rum96] Martin Rumpf. A variational approach to optimal meshes. *Numerische Mathematik*, 72(4):523-540, 1996.

[TGRL13] Thomas Toulorge, Christophe Geuzaine, Jean-François Remacle, and Jonathan Lambrechts. Robust untangling of curvilinear meshes. *Journal of Computational Physics*, 254:8-26, 2013.

[VGMN21] Jendrik Voss, Ionel-Dumitrel Ghiba, Robert J Martin, and Patrizio Neff. Sharp rank-one convexity conditions in planar isotropic elasticity for the additive volumetric-isochoric split. *Journal of Elasticity*, 143(2):301-335, 2021.

[Win66] Alan M Winslow. Numerical solution of the quasilinear poisson equation in a nonuniform triangle mesh. *Journal of computational physics*, 1(2):149-172, 1966.
<figure><img src="image_25.png" /><figcaption>Fig. 1.23: Classic global parameterization (QuadCover) integrates a frame field (upper left): the quad mesh (lower-left) needs a backup solution for a loss of local injectivity. The variational smoother (right) does not have this problem.</figcaption></figure>

<figure><img src="image_26.png" /><figcaption>Fig. 1.24: The variational smoother can generate a double covering (upper left). It produces an extra singular vertex in the quad remesh (lower left). We can solve this issue using our free boundary version.</figcaption></figure>