Vladimir Garanzha

Igor Kaporin

Liudmila Kudryavtseva

François Protais

Nicolas Ray

Dmitry Sokolov

ON LOCAL INVERTIBILITY AND QUALITY OF FREE-BOUNDARY DEFORMATIONS

Keywords: mesh untangling, variational method, polyconvex functional, penalty technique

. Namely, we show the way to guarantee absence of k-covering folds, and so, the local invertibility is assured. We demonstrate enhanced stability of suggested untangling technique which has a potential to make untangling a routine operation over meshes.

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. 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 [START_REF] Winslow | Numerical solution of the quasilinear Poisson equation in a nonuniform triangle mesh[END_REF]. Meanwhile, 4 years before that, Crowley suggested similar ideas in a classified Los Alamos Lab research report [START_REF] Crowley | An equipotential zoner on a quadrilateral mesh[END_REF] which was made public many years later.

In 1972 Godunov [START_REF] Godunov | The use of moving meshes in gas-dynamical computations[END_REF] suggested to control mesh deformations using composition with prescribed invertible mappings thus combining adaptation and guaranteed invertibility. In 1988 Jacquotte [START_REF] Jacquotte | A mechanical model for a new grid generation method in computational fluid dynamics[END_REF] 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 [START_REF] Ivanenko | Construction of nondegenerate grids[END_REF] 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 [START_REF] Charakhch'yan | A variational form of the Winslow grid generator[END_REF].

In 1966 Reshetniak [START_REF] Reshetnyak | Bounds on moduli of continuity for certain mappings[END_REF] 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 [START_REF] Liu | On the shape of tetrahedra from bisection[END_REF]. Note that importance of relations between different shape measures and condition number was underlined in [START_REF] 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[END_REF].

In 1976 J. Ball introduced his theory of finite elasticity based on the concept of polyconvex distortion energies [START_REF] Ball | Convexity conditions and existence theorems in nonlinear elasticity[END_REF]. 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 [START_REF] Ball | Global invertibility of Sobolev functions and the interpenetration of matter[END_REF]. 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 [START_REF] Rumpf | A variational approach to optimal meshes[END_REF]. In [START_REF] Garanzha | The barrier method for constructing quasi-isometric grids[END_REF][START_REF] Garanzha | Variational method for untangling and optimization of spatial meshes[END_REF] 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 [START_REF] Knupp | Hexahedral and tetrahedral mesh untangling[END_REF][START_REF] Freitag | Local optimizationbased simplicial mesh untangling and improvement[END_REF][START_REF] Escobar | Simultaneous untangling and smoothing of tetrahedral meshes[END_REF][START_REF] Toulorge | Robust untangling of curvilinear meshes[END_REF]. 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 [START_REF] Danczyk | Finite element analysis over tangled simplicial meshes: Theory and implementation[END_REF].

The main idea of untangling approach proposed in [START_REF] Garanzha | Foldover-free maps in 50 lines of code[END_REF] 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 [START_REF] Du | Lifting Simplices to Find Injectivity[END_REF] as a major challenge for computational graphics and augmented reality problems. Authors of [START_REF] Du | Lifting Simplices to Find Injectivity[END_REF] 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 [START_REF] Du | Lifting Simplices to Find Injectivity[END_REF] also breaks down when randomized initial guesses for their test set are used, while the results of algorithm from [START_REF] Garanzha | Foldover-free maps in 50 lines of code[END_REF] are not sensitive to initial guesses.

For fixed boundary problems the approach has a sound theoretical basis. Namely, in the aforementioned work [1, §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 [START_REF] Garanzha | Foldover-free maps in 50 lines of code[END_REF] (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 [START_REF] Garanzha | Regularization of the barrier variational method[END_REF] when crossing the barrier is eliminated.

However algorithm from [START_REF] Garanzha | Foldover-free maps in 50 lines of code[END_REF] 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 inversionfree 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π for 2D and above 4π 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 § 2.1 we recall the untangling method proposed in [START_REF] Garanzha | Foldover-free maps in 50 lines of code[END_REF], then in § 2.2 we discuss the structure of the admissible set, and we show the limitations of the method for free boundary deformations. Next in § 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 § 4.

FOLDOVER-FREE DEFORMATIONS

In this section we go over the method proposed in the paper [START_REF] Garanzha | Foldover-free maps in 50 lines of code[END_REF]. First, in § 2.1 we recapitulate the main idea and the resolution scheme, and then in § 2.2 we show the main limitation: even without inverted elements, the deformation can be non-invertible.

Variational formulation for grid generation

Let us go through the algorithm to compute a foldoverfree deformation x(ξ)

: Ω ⊂ R d → R d .
x( ξ) Ω f (J) dξ, (1)
where J is the Jacobian matrix of the mapping x( ξ), and density of deformation energy is defined by [START_REF] Garanzha | The barrier method for constructing quasi-isometric grids[END_REF]

f (J) = (1 -θ)fs(J) + θfv(J), (2)
where shape distortion is defined as

fs(J) :=      1 d tr J J (det J) 2 d , det J > 0 +∞, det J ≤ 0 (3)
while volumetric distortion is defined

fv(J) :=    1 2 det J + 1 det J , det J > 0 +∞, det J ≤ 0 (4)
Prob.

(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 fs(J) and fv(J) have concurrent goals, one preserves angles and the other preserves the area, and thus θ 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 [START_REF] Flory | Thermodynamic relations for high elastic materials[END_REF], [START_REF] Penn | Volume Changes Accompanying the Extension of Rubber[END_REF], [START_REF] Borst | Modelling and analysis of rubberlike materials[END_REF] and remains a hot topic in convexity analysis [START_REF] Voss | Sharp rank-one convexity conditions in planar isotropic elasticity for the additive volumetric-isochoric split[END_REF].

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 φ(J) : R d×d → R ∪ + ∞ is said to be polyconvex [START_REF] Ball | Convexity conditions and existence theorems in nonlinear elasticity[END_REF] if there exists a convex function Φ(#J), such that φ(J) = Φ(#J), where #J denotes the set of all minors of J. Any polyconvex function φ is rank-one convex [START_REF] Ball | Convexity conditions and existence theorems in nonlinear elasticity[END_REF], i.e.

D χ(D, ε) ε = 3 ε = 4 ε = 2 ε = 1
φ((1 -θ)J + θ(J + δJ) ≤ (1 -θ)φ(J) + θφ(J + δJ),
where rank δJ = 1, and thus satisfies the Hadamard-Legendre conditions (ellipticity conditions for the Euler-Lagrange equation of variational Prob.( 1)

) d i,k,j,m=1 ∂ 2 φ ∂(J)ij ∂(J) km pip k qjqm ≥ 0 for arbitrary vectors p, q ∈ R d .
Since density of energy ( 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

Ω det J{ x}( ξ) > 0.
Note that if an initial guess is not admissible, then functional (1) is not defined. In this case a regularized version can be introduced [START_REF] Garanzha | Regularization of the barrier variational method of grid generation[END_REF]: we can avoid nonpositive denominators in fs and fv using a regularization function χ for a positive value of ε (Fig. 1):

χ(D, ε) := D + √ ε 2 + D 2 2 (5)
Then one can define a regularized version fε of functions f :

fε(J) := (1 -θ) 1 d tr J J (χ(det J, ε)) 2 d + θ 1 2 1 + det 2 J χ(det J, ε) , (6)
so that Prob. ( 1) is reformulated as

lim ε→0 + arg min x( ξ) Ω fε(J) dξ (7)
The Prob. ( 7) does offer a way of getting rid of foldovers if a foldover-free initialization is not available.

In practice, the map 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 { xi} #V i=1 . Let us denote the vector of all variables as X := x 1 . . . x #V , then optimization Prob. (1) has following expression in discrete settings:

lim ε→0 + arg min X F (X, ε), (8)
where

F (X, ε) := #T k=1 fε(J k ) 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.

We denote by P k the k-th simplex of the computational mesh. Since it is an affine image of T k

P k = x(T k ), we get det J k = vol P k / vol T k .
We assume that vol T k > 0 for all parametric tetrahedra.

To solve Prob. ( 8), we use an iterative descent method, as suggested in [START_REF] Garanzha | Foldover-free maps in 50 lines of code[END_REF]. Starting from an initial guess X 0 , they build a sequence of approximations X k+1 := X k + ∆X k , carefully choosing the regularization parameter ε k for each iteration k.

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 vol

(P k ) vol(T k ) > 0, k = 1, . . . , #T (9)
which is finite-dimensional approximation of the set

det J > 0. ( 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 (9) defines a nonconvex set. One can hardly expect that intersection of the sets in [START_REF] Liu | On the shape of tetrahedra from bisection[END_REF] would result in a convex domain. Moreover, Ciarlet [START_REF] Ciarlet | Sur les lois de comportement en elasticite non-lineaire compressible[END_REF] 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 A and

S(X) → +∞ when X ∈ A, X → ∂ A,
then it cannot be convex function of X. Fortunately, barrier distortion measures can be polyconvex, as shown by J. Ball [START_REF] Ball | Convexity conditions and existence theorems in nonlinear elasticity[END_REF]. For instance, it means that domain A consists of connected components where each two points can be connected. Consider vector

X + sY ∈ R #V d , 0 ≤ s ≤ 1, X ∈ A, X + Y ∈ 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 , rank B k = 1. (11)
Since -det J is rank one convex function of J, we get

-det(J k + sB k ) ≤ -(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 A where any pair of points be connected using sequence of rank-one segments.

The admissible set 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 2. 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 2(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 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 3 Homeomorphism between these two solutions does not exist, so they definitely belong to disjoint subsets of 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-toone mapped on the open neighborhood of the flat projection.

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 4. 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 kcoverings 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 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 (9) of admissible set 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 Figure 6 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 6, one first level triangle T

(1) 1 and two second level triangles T (2)

i , i = 1, 2 are added.
The irregular boundary case is shown in Figure 7 (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. Intersection of i-th cone with unit sphere defines the spherical triangle Ti with angles αi, βi, γi. These angles coincide with three dihedral angles of the cone/tet. The number of spherical triangles nt is equal to the number of tetrahedra adjacent to p, while the number nv of vertices of spherical Hence we consider quasi-2D algorithm which prohibits 2-coverings for tetrahedral fans around each edge of the tet mesh. Figure 8 shows sample 2-covering for tetrahedra around an edge. 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 10(left)) there exists internal edge eint 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 eint and boundary triangles Ti creating tets as conv(eint, Ti). Edgebased 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π. For irregular case shown in Figure 10(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 * , pi, pi+1 for all boundary vertices pi of the boundary star. As a result star of the vertex p can be treated as an internal star.

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 eint, and total dihedral angle around this edge is equal to 2π . Hence edge eint 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. Generally recursive aggregation serves to reduce the admissible set 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 11 (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 12 (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 12 (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.

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 13(a) provides the rest shape, the vertices to be swapped are highlighted in red and blue, respectively. Figure 13(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 13(c)).

Figure 14 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 is executed without adding phantom triangles, the mesh gets caught into 2-covering traps (Figure 14(b)). Finally, Figure 14(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 15(a) shows the input mesh to flatten, the mesh is highly curved and presents very bad quality elements. Refer to Then we show that the same problem exists in 3D as well. Figure 16 is a stress test for the free boundary untangling, similar to rotated cube-inside-cube stress test considered in [START_REF] Garanzha | Variational method for untangling and optimization of spatial meshes[END_REF][START_REF] Garanzha | Foldover-free maps in 50 lines of code[END_REF]. 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 16, 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.

Finally we show that it is possible to go beyond simple deformations. Figures 17 and18 show a quadremeshing 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 [START_REF] Kälberer | QuadCover -Surface Parameterization using Branched Coverings[END_REF] to do so often results in a local loss of injectivity and, as illustrated in Fig- ure 17-left, one may need to call a backup solution such as QEx [START_REF] Ebke | QEx: Robust Quad Mesh Extraction[END_REF] to extract a quad mesh despite the tangled parameterization. Free-boundary variational smoother overcomes this obstacle, producing valid parameterizations (Figure 17-right). We need, however, protect the mesh to avoid k-coverings, as illustrated in Figure 18. 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 [START_REF] Ivanenko | Construction of nondegenerate grids[END_REF].

CONCLUSION

We formulate a set of variational problems potentially covering the complete technological 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 deformations. We illustrate performance of our algorithm with challenging 2D and 3D numerical tests.

Figure 1 :

 1 Figure 1: Regularization function for the denominator in Eq. (6). When ε tends to zero, χ(ε, D) tends to D for positive values of D, and to 0 + for negative values of D.

  (a). While standard flattening would be regular 12-gon, shown in Figure3(b), minimal deformation energy is provided by "2-covering", shown in Fig. Figure3(c).

Figure 2 :Figure 3 :

 23 Figure 2: 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.

Figure 4 :

 4 Figure 4: 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.

1

 1 as shown in Figure5(a), (b).

Figure 5 :

 5 Figure 5: Recursive aggregation of adjacent triangles to create phantom protective triangles. (a) Initial 8stencil, (b) 1 level blue phantom triangle is added, (c) three 2 level phantom triangles are added resulting in a 4-star.

Figure 6 :

 6 Figure 6: Regular boundary star: recursive adding of phantom triangles results in a 2-star.

Figure 7 :

 7 Figure 7: Irregular boundary star: configuration of the boundary corresponds to sharp incoming corner. Single outer phantom triangle transforms boundary star into internal star.

Figure 8 :

 8 Figure 8: Left: a regular fan of 6 tetrahedra sharing the common edge M N . Right: a double covering around the edge M N .

Figure 9

 9 Figure 9 illustrates aggregation algorithm in 3D. The set of tetrahedra around an edge defines 2D triangle

Figure 9 :

 9 Figure 9: Adding phantom tetrahedra around an edge is equivalent to a 2D problem.

Figure 10 :

 10 Figure 10: (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.

Figure 11 :

 11 Figure 11: (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.

Figure 12 :

 12 Figure 12: (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.

  Figures 15(a-b) for the untangling with and without protecting elements.

Figure 13 :

 13 Figure 13: 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.

Figure 14 :Figure 15 :

 1415 Figure 14: Inverse metal sheet forming. (a): input model, (b): unprotected untangled mesh contains 2-covering traps, (c): protected mesh defines locally invertible flattening.

Figure 16 :

 16 Figure 16: 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 is globally invertible.

Figure 17 :

 17 Figure 17: Classic global parameterization (Quad-Cover) 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.

Figure 18 :

 18 Figure 18: 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.