Vladimir Garanzha

Igor Kaporin

Liudmila Kudryavtseva

François Protais
email: francois.protais@inria.fr

Nicolas Ray

Dmitry Sokolov

Dorodnicyn

Foldover-free maps in 50 lines of code

Fig. 1. Our method of constructing injective maps opens a door for a large variety of applications. This figure shows an example of a thick prismatic mesh layer (shown in green) built around a triangulated surface, a very challenging problem for highly curved objects. Thanks to our method, we are able to compute such a layer free of folds and self-intersections.

Mapping a triangulated surface to 2D space (or a tetrahedral mesh to 3D space) is an important problem in geometry processing. In computational physics, untangling plays an important role in mesh generation: it takes a mesh as an input, and moves the vertices to get rid of foldovers. In fact, mesh untangling can be considered as a special case of mapping where the geometry of the object is to be defined in the map space and the geometric domain is not explicit, supposing that each element is regular. In this paper, we propose a mapping method inspired by the untangling problem and compare its performance to the state of the art. The main advantage of our method is that the untangling aims at producing locally injective maps, which is the major challenge of mapping. In practice, our method produces locally injective maps in very difficult settings, both in 2D and 3D. We demonstrate it on a large reference database as well as on more difficult stress tests. For a better reproducibility, we publish the code in Python for a basic evaluation, and in C++ for more advanced applications.

CCS Concepts: • Computing methodologies → Mesh models.

Additional Key Words and Phrases: Parameterization, injective mapping, mesh untangling

INTRODUCTION

Most geometric objects are represented by a triangulated surface or a tetrahedral mesh. The mapping problem consists in generating a 2D or 3D map of these objects. This is a fundamental problem of computer graphics because it is much easier for many applications to work in this map space than to directly manipulate the object itself. To give few examples, texture mapping stores colors of a surface as images in the map space, remeshing uses global maps in 2D [START_REF] Bommes | Integer-Grid Maps for Reliable Quad Meshing[END_REF] and 3D [START_REF] Gregson | All-Hex Mesh Generation via Volumetric PolyCube Deformation[END_REF][START_REF] Nieser | CubeCover -Parameterization of 3D Volumes[END_REF]]. In addition, mapping algorithms can be used to deform enforcing the constraints. This allows them to optimize at the same time the parameterization and the texture packing [START_REF] Jiang | Simplicial Complex Augmentation Framework for Bijective Maps[END_REF], with a possibility to scale to larger meshes [START_REF] Rabinovich | Scalable Locally Injective Mappings[END_REF].

Recover local injectivity.

Local injectivity can also be recovered for a map with few foldovers present. For example, in 2D [START_REF] Lipman | Bounded Distortion Mapping Spaces for Triangular Meshes[END_REF]] and 3D [START_REF] Aigerman | Injective and Bounded Distortion Mappings in 3D[END_REF], the map is projected on a class of bounded distortion maps. Alternating between projection and optimization steps [START_REF] Fu | Computing Inversion-Free Mappings by Simplex Assembly[END_REF][START_REF] Shahar | Large-scale bounded distortion mappings[END_REF][START_REF] Naitsat | Adaptive Block Coordinate Descent for Distortion Optimization[END_REF][START_REF] Su | Practical Foldover-Free Volumetric Mapping Construction[END_REF]] often allows to produce foldover-free maps with low distortion, as compared to our results in §3.2. The numerical methods are however unlikely to succeed for stiff problems.

Recovering local injectivity is also known as mesh untangling. Originally related to Arbitrary Lagrangian-Eulerian moving mesh approach, the mesh untangling problem considers a simplicial complex with misoriented elements and attempts to flip them back by optimizing the position of the vertices. There is an abundant literature on mesh untangling [START_REF] Du | Lifting Simplices to Find Injectivity[END_REF][START_REF] Escobar | Simultaneous untangling and smoothing of tetrahedral meshes[END_REF][START_REF] Freitag | Local optimization-based simplicial mesh untangling and improvement[END_REF][START_REF] Patrick | Hexahedral and tetrahedral mesh untangling[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 [START_REF] Danczyk | Finite element analysis over tangled simplicial meshes: Theory and implementation[END_REF] investigates a finite element method working directly on tangled (sic!) meshes.

Elastic deformations. To recover local injectivity, we propose a method stemming from the computational physics. It is very important to note that there is rich literature on mesh deformation in the community working on grid generation for scientific computation. Numerical simulation of hydrodynamic instability of layered structures requires sound mathematical foundations behind moving deforming mesh algorithms. In the '60s Winslow and Crowley, independently one from another, introduced mesh generation methods based on inverse harmonic maps [START_REF] Crowley | An equipotential zoner on a quadrilateral mesh[END_REF][START_REF] Winslow | Numerical solution of the quasilinear Poisson equation in a nonuniform triangle mesh[END_REF].

Since then, a lot of effort was spent on mesh generation based on elastic deformations [START_REF] Jacquotte | A mechanical model for a new grid generation method in computational fluid dynamics[END_REF]], but mostly for regular grids. In 1988, at the time of domination of finite difference mapped grid generation methods, S. Ivanenko introduced the pioneering concept of barrier variational grid generations methods guaranteeing construction of non-degenerate grids [START_REF] Charakhch | A variational form of the Winslow grid generator[END_REF][START_REF] Ivanenko | Generation of non-degenerate meshes[END_REF]]. To generate deformations with bounded global distortion (bounded quasi-isometry constant), Garanzha proposed to minimize an elastic energy for a hyperelastic material with stiffening suppressing singular deformations [START_REF] Va Garanzha | The barrier method for constructing quasi-isometric grids[END_REF]]. Invertibility theorem for deformation of this material was established in the 3D case as well [START_REF] Garanzha | Variational method for untangling and optimization of spatial meshes[END_REF].

A solid mathematical ground for these methods was laid by J. Ball who introduced in 1976 his theory of finite elasticity based on the concept of polyconvex distortion energies [START_REF] John | 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 to note 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]].

Our contributions. Inspired by these results on untangling and elastic deformations, we propose a simple numerical recipe for fixed connectivity mesh untangling. This method outperforms recent state of the art on locally injective parameterization [START_REF] Du | Lifting Simplices to Find Injectivity[END_REF] in terms of robustness, quality and supported features.

The main idea of the untangling algorithm is to start from an arbitrary mesh and get directly to the optimal mesh past the barrier delimiting the admissible set. Mesh elements are penalized according to negative value of signed volume. While the general concepts used in our algorithm (deformation of hyperelastic material and penalization of folds) are known since the 80s, a transition to a working algorithm is not trivial. We are the first to actually bring guarantees for the mesh untangling problem.

To sum up, our main contribution is the algorithm whose ingredients are carefuly chosen to get well-founded and well-behaved subproblems:

• This algorithm contains the very first penalization strategy with theoretical guarantees of possibility of untangling in a finite number of steps. • In addition to that, we guarantee boundedness of the deformation energy and the fact that positive definite part of the Hessian matrix is spectrally equivalent to finite element stiffness matrix for Laplace operator. This guarantees absence of extremely stiff problems when crossing the barrier. • Last, but not least, is the "in 50 lines of code". The solution we propose is simple, with virtually no tuning required from the user, leading to easiness of adoption of our mapping approach by a potentially wide range of applications. To ease the reproducibility, we publish the code in Python (refer to Listing 1) for a basic evaluation, and a C++ code in the supplemental material [START_REF] Sokolov | Supplemental material for "Foldover-free maps in 50 lines of code[END_REF]] for more advanced applications.

The rest of the paper is organized as follows: we start with presenting our method in § 2, then we evaluate its performance ( § 3.1 and § 3.2) as well as its limitations ( § 3.3). Then we present theoretical guarantees for our resolution scheme: in § 4.2 we prove that our approximation of Hessian matrix is positive definite, and finally in § 4.3 we prove that our choice of the regularization parameter sequence guarantees that a minimization algorithm1 can find a mesh free of inverted elements in a finite number of steps.

PENALTY METHOD FOR MESH UNTANGLING

In this section we present our method of computing a foldover-free map ì 𝑢 : Ω ⊂ R 𝑑 → R 𝑑 , i.e. we map the domain Ω to a parametric domain. This presentation is unified both for 2D and 3D settings, and by 𝑑 we denote the number of dimensions; in our notations we use arrows for all vectors of dimension 𝑑.

The section is organized as follows: in § 2.1 we give a primer on the variational formulation of mapping problem in continuous settings, then we state our problem in § 2.2 as a regularization of this variational formulation, and finally we present our numerical resolution scheme in § 2.3.

Variational formulation for grid generation

Let us denote by ì 𝑢 ( ì 𝑥) a map to a parametric domain: for the flat 2D case we can write as ì 𝑢 (𝑥, 𝑦) = (𝑢 (𝑥, 𝑦), 𝑣 (𝑥, 𝑦)), and for a 3D map ì 𝑢 (𝑥, 𝑦, 𝑧) = (𝑢 (𝑥, 𝑦, 𝑧), 𝑣 (𝑥, 𝑦, 𝑧), 𝑤 (𝑥, 𝑦, 𝑧)).

Consider the following variational problem:

arg min ì 𝑢 ∫ Ω (𝑓 (𝐽 ) + 𝜆𝑔(𝐽 )) 𝑑𝑥, ( 1
)
where 𝐽 is the Jacobian matrix of the mapping ì 𝑢 ( ì 𝑥), and

𝑓 (𝐽 ) := tr 𝐽 ⊤ 𝐽 (det 𝐽 ) 2 𝑑 , det 𝐽 > 0 +∞, det 𝐽 ≤ 0 𝑔(𝐽 ) := det 𝐽 + 1 det 𝐽 , det 𝐽 > 0 +∞,
det 𝐽 ≤ 0 Problem (1) may be subject to some constraints that we do not write explicitly. To give an example, one may pin some points in the map. In this formulation, functions 𝑓 (𝐽 ) and 𝑔(𝐽 ) have concurrent goals, one preserves angles and the other preserves the area, and thus 𝜆 serves as a trade-off parameter.

As a side note, with 𝜆 = 0 and 𝑑 = 2, Prob.

(1) presents a variational formulation of an inverse harmonic map problem. Namely, if we write down the Euler-Lagrange equations for Prob. (1) and interchange the dependent and independent variables2 , we obtain the Laplace equation Δì 𝑥 ( ì 𝑢) = ì 0 (not to be confused with omnipresent Δì 𝑢 ( ì 𝑥) = ì 0!). For this case, Prob. ( 1) is often referred to as Winslow's functional, however Winslow himself has never formulated the variational problem, working with inverse Laplace equations. To the best of our knowledge, the first publication of the variational problem is made by [START_REF] Brackbill | Adaptive zoning for singular problems in two dimensions[END_REF].

There is a rich history behind both terms 𝑓 and 𝑔 of the deformation energy (1), and most publications were overlooked by the computer graphics community. The idea of isochoric-volumetric split, where deformation energy is represented as 𝑓 𝑖 𝐽 /(det 𝐽 ) 1/𝑑 + 𝑓 𝑣 (det 𝐽 ) goes back to the 60s [START_REF] Flory | Thermodynamic relations for high elastic materials[END_REF][START_REF] Penn | Volume Changes Accompanying the Extension of Rubber[END_REF]. Function 𝑓 was also introduced in the theory of 𝑑-dimensional quasi-conformal mappings in the 60s, for example, Reshetnyak [START_REF] Yu | Bounds on moduli of continuity for certain mappings[END_REF] established precise relations between 𝑓 and standard matrix condition number. [START_REF] Ivanenko | Generation of non-degenerate meshes[END_REF]] was the first to use 𝑓 as it is for 2D meshing. It was brought to the graphics community by [START_REF] Hormann | MIPS: An Efficient Global Parametrization Method[END_REF]. First 3D implementation of 𝑓 for mesh generation can be attributed to [Knupp 2000a], however earlier 3D implementations of similar split-based hyperelastic energies for deformations of solids date back to 1988 [START_REF] De Borst | Modelling and analysis of rubberlike materials[END_REF]]. It seems that [START_REF] Va Garanzha | The barrier method for constructing quasi-isometric grids[END_REF]] was the first to use 𝑔 for meshing.

Penalty method

Prob. (1) is known for decades, it provides a simple and efficient tool to optimize the quality of the map. Note that the energy (1) is a polyconvex function (refer to App. B-Rem. 3 for a proof) so Euler-Lagrange equations for optimal deformation satisfy the ellipticity conditions, and therefore it is very well suited for a numerical optimization provided that we have an initial guess in the admissible domain min

Ω 𝐽 ( ì 𝑢) > 0.
The problem, however, is that while being theoretically sound, this problem statement does not offer any practical way to get rid of foldovers in a map, because for a map with foldovers the energy is infinite and provides no indications on how to improve the situation, hence we propose to alter a little the problem statement.

Numerical optimization is hard due to non-positive denominators in 𝑓 and 𝑔. For example, for quasiconformal maps in 2D the dilatation coefficient allows [START_REF] Weber | Computing Extremal Quasiconformal Maps[END_REF]] to avoid them. Here, we need a more general strategy of folds penalization. The idea goes back to [START_REF] Ivanenko | Generation of non-degenerate meshes[END_REF]]. In this paper we improve ( § 4.3) the heuristic penalty technique proposed in [START_REF] Va Garanzha | Regularization of the barrier variational method[END_REF]. Namely, we use a regularization function 𝜒 for a positive value of 𝜀 (Fig. 2):

𝜒 (𝐷, 𝜀) := 𝐷 + √ 𝜀 2 + 𝐷 2 2 , (2)
Then we define a regularized version 𝑓 𝜀 , 𝑔 𝜀 of functions 𝑓 and 𝑔:

𝑓 𝜀 (𝐽 ) := tr 𝐽 ⊤ 𝐽 (𝜒 (det 𝐽, 𝜀)) 2 𝑑 , 𝑔 𝜀 (𝐽 ) := det 2 𝐽 + 1 𝜒 (det 𝐽, 𝜀) , (3)
so that Prob (1) is reformulated as

lim 𝜀→0 + arg min ì 𝑢 ∫ Ω (𝑓 𝜀 (𝐽 ) + 𝜆𝑔 𝜀 (𝐽 )) 𝑑𝑥 (4)
Under certain assumptions3 solutions of Prob. (4) are solutions of Prob. (1), however, Prob. (4) does offer a way of getting rid of foldovers if a foldover-free initialization is not available.

In practice, the map ì 𝑢 is piece-wise affine with the Jacobian matrix 𝐽 being piece-wise constant, and can be represented by the coordinates of the vertices in the parametric domain {ì 𝑢 𝑖 } #𝑉 𝑖=1 . Let us denote the vector of all variables as 𝑈 :

= ì 𝑢 ⊤ 1 . . . ì 𝑢 ⊤ #𝑉 ⊤
, then our optimization problem can be discretized as follows:

lim 𝜀→0 + arg min 𝑈 𝐹 (𝑈 , 𝜀), (5)
where 𝐹 (𝑈 , 𝜀)

:= #𝑇 𝑡 =1 (𝑓 𝜀 (𝐽 𝑡 ) + 𝜆𝑔 𝜀 (𝐽 𝑡 )) vol(𝑇 𝑡 ),
#𝑉 is the number of vertices, #𝑇 is the number of simplices, 𝐽 𝑡 is the Jacobian matrix for the simplex 𝑡 and vol(𝑇 𝑡 ) is the volume of the simplex 𝑇 𝑡 in the original domain.

Resolution scheme

To solve Prob. ( 5), we use an iterative descent method. We start from an initial guess 𝑈 0 , and we build a sequence of approximations 𝑈 𝑘+1 := 𝑈 𝑘 + Δ𝑈 𝑘 . For each iteration we need to carefully choose the regularization parameter 𝜀 𝑘 . Starting from 𝜀 0 := 1, we define the sequence as follows:

𝜀 𝑘+1 := 2 𝜇 𝑘 (𝜇 𝑘 -𝐷 𝑘+1 -) if 𝐷 𝑘+1 -< 𝜇 𝑘 0 if 𝐷 𝑘+1 -≥ 𝜇 𝑘 , (6)
where 𝐷 𝑘+1 -:= min

𝑡 ∈1...#𝑇 det 𝐽 𝑘+1 𝑡
is the minimum value of the Jacobian determinant over all cells of the mesh at the iteration 𝑘 + 1, 𝜎 𝑘 := max 1 10 , 1 -

𝐹 (𝑈 𝑘+1 ,𝜀 𝑘 ) 𝐹 (𝑈 𝑘 ,𝜀 𝑘 )
is descent coefficient and

𝜇 𝑘 := (1 -𝜎 𝑘 ) 𝜒 (𝐷 𝑘+1
-, 𝜀 𝑘 ). This formula is justified by Th. 1 ( § 4.3) on finite untangling sequence.

The simplest way to find Δ𝑈 𝑘 is to call a quasi-Newtonian solver such as L-BFGS [START_REF] Liu | On the Limited Memory BFGS Method for Large Scale Optimization[END_REF]. The only thing we need to implement is the computation of the function 𝐹 (𝑈 𝑘 , 𝜀 𝑘 ) and its gradient ∇𝐹 (𝑈 𝑘 , 𝜀 𝑘 ). It should also be possible to adapt solvers specialized in geometry processing problems [START_REF] Shtengel | Geometric Optimization via Composite Majorization[END_REF][START_REF] Smith | Analytic Eigensystems for Isotropic Distortion Energies[END_REF][START_REF] Zhu | Blended Cured Quasi-Newton for Distortion Optimization[END_REF].

Another option is to compute analytically the Hessian matrix instead of estimating it. The problem, however, is that the Hessian matrix 𝜕 2 𝐹 𝜕 𝑈 𝜕 𝑈 ⊤ is not positive definite. In this paper we propose its approximation that ensures the positive definiteness. The modified Hessian matrix 𝐻 + (𝑈 𝑘 , 𝜀 𝑘 ) of the function 𝐹 with respect to 𝑈 at the point 𝑈 𝑘 is built out of 𝑑 × 𝑑 blocks

𝐻 + 𝑖 𝑗 ≈ 𝜕 2 𝐹 𝜕 ì 𝑢 𝑖 𝜕 ì 𝑢 ⊤ 𝑗 (𝑈 𝑘 , 𝜀 𝑘 ).
Here, the matrix 𝐻 + 𝑖 𝑗 is placed on the intersection of 𝑖-th block row and 𝑗-th block column; the ≈ symbol means that we remove all the terms depending on the second derivative of 𝜒 and second derivatives of det 𝐽 to keep 𝐻 + positive definite. Refer to Appendix A for the formulae, and to § 4.2 for the proof of the positive definiteness of 𝐻 + .

A detailed description of the resolution scheme is given in Alg. 1. Refer to List. 1 and Fig. 3 for a complete working example of Python implementation and the corresponding input/output generated by the code. Note that our method is not limited to simplicial meshes only: in this particular example we evaluate the Jacobian matrix for every triangle forming quad corners, what corresponds to the trapezoidal quadrature rule.

RESULTS AND DISCUSSION

In this section we provide an experimental evaluation of the method.

In the field of computer graphics, any claim about map injectivity always faces a simple sanity check (Fig. 4): take a square and swap any two inner points. Our method successfully avoids the deskreject, so we start this section ( § 3.1) by testing our method on the benchmark [START_REF] Du | Lifting Simplices to Find Injectivity[END_REF]], then we continue with further tests we have found relevant ( § 3.2), and finally we discuss the limitations of the approach in § 3.3.

Author version, 2021.

Initialization

Our result [START_REF] Du | Lifting Simplices to Find Injectivity[END_REF]]

min det 𝐽 ≈ 4 • 10 -4 max 𝜎 1 (𝐽 ) 𝜎 2 (𝐽 ) ≈ 69 min det 𝐽 ≈ 4 • 10 -4 max 𝜎 1 (𝐽 ) 𝜎 2 (𝐽 ) ≈ 72 min det 𝐽 ≈ 4 • 10 -4 max 𝜎 1 (𝐽 ) 𝜎 2 (𝐽 ) ≈ 70 min det 𝐽 ≈ 1 • 10 -20 max 𝜎 1 (𝐽 ) 𝜎 2 (𝐽 ) ≈ 24104 min det 𝐽 ≈ 9 • 10 -17 max 𝜎 1 (𝐽 )
𝜎 2 (𝐽 ) ≈ 16500 max injective maps produced by our method, columns (e) and (f): injective maps produced by [START_REF] Du | Lifting Simplices to Find Injectivity[END_REF]]. The method by [START_REF] Du | Lifting Simplices to Find Injectivity[END_REF]] fails to generate injective maps for 135°and 180°inner cube rotations. The color map illustrates the relative volume scaling: green for det 𝐽 ≈ 1, red for inflation, blue for compression.

𝜎 1 (𝐽 ) 𝜎 3 (𝐽 ) ≈ 37 min det 𝐽 ≈ 0.03 (a) max 𝜎 1 (𝐽 ) 𝜎 3 (𝐽 ) ≈ 58 min det 𝐽 ≈ 0.02 (b) max 𝜎 1 (𝐽 ) 𝜎 3 (𝐽 ) ≈ 759 min det 𝐽 ≈ 0.003 (c) max 𝜎 1 (𝐽 ) 𝜎 3 (𝐽 ) ≈ 659 min det 𝐽 ≈ 0.002 (d) max 𝜎 1 (𝐽 ) 𝜎 3 (𝐽 ) ≈ 524 min det 𝐽 ≈ 0.002 (e) max 𝜎 1 (𝐽 ) 𝜎 3 (𝐽 ) ≈ 3391 min det 𝐽 ≈ 0.0004 (f)
Author version, 2021.

102:6 • Vladimir Garanzha, Igor Kaporin, Liudmila Kudryavtseva, François Protais, Nicolas Ray, and Dmitry Sokolov

Benchmark database

As we have already mentioned, recently Du et al. have published a method for building locally injective maps for constrained boundary problems. There is an elegant idea behind Total Lifted Content (TLC): the authors propose to minimize the total unsigned area of a mesh to untangle. This energy is regularized by lifting simplices to a higherdimensional space. After the regularization the energy becomes smooth, and every global minimum of the energy is achieved by an injective embedding.

Along with their paper, Du et al. have published a valuable benchmark database. It contains a huge number of 2D and 3D constrained boundary injective mapping challenges. To the best of our knowledge, TLC [START_REF] Du | Lifting Simplices to Find Injectivity[END_REF]] and our method are the only ones passing the benchmark without any fail. For 2D challenges, the benchmark contains 3D triangulated surfaces to flatten, and not flat 2D meshes as we have described in § 2.1. Nevertheless, our method can handle it directly because the mapping is still R 2 → R 2 on each triangle.

A representative example from the database is given in top row of Fig. 5. The challenge is to map the "Lucy" mesh statuette from the Stanford Computer Graphics Laboratory to a P-shaped domain. This mesh has a topology of a disk, and its boundary vertices are uniformly spaced on the P-shape boundary. As an initialization to the problem, Du et al. have computed the corresponding Tutte embedding that obviously contains a foldover (Fig. 5-top left). Then the problem boils down to a mesh untangling with locked boundary.

Mapping quality measure. How to measure quality of a map? Well, it depends on the goal. An identity is an unreachable ideal; traditional competing goals are (as much as possible) angle preserving and area preserving maps. We can measure the extreme values of the failure of a map to be conformal or authalic. Our maps being piece-wise affine, the Jacobian matrix 𝐽 is constant per element. Let us define the largest singular value of 𝐽 as 𝜎 1 (𝐽 ), and the smallest singular value as 𝜎 𝑑 (𝐽 ); then the quality of a mapping can be reduced to extreme values of the stretch (max

𝜎 1 (𝐽 )
𝜎 𝑑 (𝐽 ) ) and the scaling (min det 𝐽 ).

For the "Lucy-to-P" challenge (Fig. 5-top row) our map differs from the TLC result by 12 orders of magnitude in terms of minimum scaling, and by two orders of magnitude in terms of maximum stretch. To visualize this difference in scaling, we have provided the close-ups: Fig. 5-top middle shows a map of the Lucy's torch, whereas the same level of zoom on the result by contains not only the torch, but also both wings, the head and the right arm! Note also that the input "Lucy" mesh is slightly anisotropic; our method allows us to prescribe the element target shape, so the dress pleats are clearly visible in our mapping.

Benchmark database. Our method successfully passes all challenges from the benchmark [START_REF] Du | Lifting Simplices to Find Injectivity[END_REF]]. The benchmark consists of 10743 meshes to untangle in 2D and 904 meshes in 3D under locked boundary constraints. In Fig. 7 we provide quality plots of the resulting locally injective maps. These are loglog scatter plots: each dot corresponds to a quality of the corresponding map reduced to two numbers: the maximum stretch (max

𝜎 1 (𝐽 )
𝜎 𝑑 (𝐽 ) ) and the minimum scaling (min det 𝐽 ). Left column of Fig. 7 shows the 2D dataset (10743 challenges) 3D dataset (904 challenges) worst quality measurements for every 2D problem (top) as well as for every 3D challenge (bottom) of the dataset. Our results are shown in blue, whereas TLC results are shown in red. To illustrate the distribution of the elements' quality, for each injective map we have removed 5% of worst measurements: the right column of Fig. 7 shows the maximum stretch and the minimum scaling for the top 95% of measurements.

Note the dot arrangements forming lines in the plot: these dots correspond to the few sequences of deformation present in the database.

Timings. Fig. 8 provides a loglog scatter plot of our running time vs mesh size for all the challenges from the database [START_REF] Du | Lifting Simplices to Find Injectivity[END_REF]]: for each run, the time varies from a fraction of a second to several minutes for the largest meshes. These times were obtained with a 12 cores i7-6800K CPU @ 3.40 GHz. As in Fig. 8, the vertical lines in the 3D dataset plot correspond to the sequences of deformation in the benchmark.

There are two scatter plots superposed, both represent the same resolution scheme with an exception corresponding to the way we compute 𝜀 𝑘 (Alg. 1-line 3). The green scatter plot corresponds to a conservative update rule (Eq. ( 6)) offering guarantees on untangling Author version, 2021. Each dot corresponds to a challenge from the database (10743 in 2D and 904 in 3D). Blue dots show the running times obtained using a heuristic regularization [Garanzha and Kaporin 1999, Eq. 6.3], green dots correspond to Eq. (6).

in a finite number of steps (refer to Th. 1), whereas the blue scatter plot is obtained using the heuristic update rule [Garanzha and Kaporin 1999, Eq. 6.3]. This formula was chosen empirically, however it performs well in the vast majority of situations. For instance, it allows for all the database [START_REF] Du | Lifting Simplices to Find Injectivity[END_REF]] to pass the injectivity test.

Further testing

Sensitivity to initialization. Our next test is the sensitivity to the initialization. We have generated two other initializations for the "Lucy-to-P" challenge: the one with all interior vertices collapsed onto a single point (Fig. 5-middle row), and with the interior vertices being randomly placed withing a bounding square (Fig. 5-bottom row).

Our method produces virtually the same result on all three initializations, whereas TLC generates very different results for the first two, and fails for the third one. It is interesting to note that TLC is heavily depending on the initialization: it alters very little the input geometry. Our experiments with the source code [START_REF] Du | Lifting Simplices to Find Injectivity[END_REF] show that most challenges from the benchmark fail under a random initialization. Moreover, the sanity check (Fig. 4) fails as well.

Large deformation stress test. For our next test we have generated an isotropic tetrahedral mesh of a cube with a cavity, and we rotated the inner boundary to test the robustness of our method to large deformations. Figure 6 shows the results. Our L-BFGS-based optimization scheme succeeds up to the rotation of 135°, and we had to switch to the Newton method to reach the 180°rotation. TLC method had succeeded on 45°and 90°, and failed for the 135°and 180°. Note that as in the previous test, even when the untangling succeeds, TLC alters very little the input map, thus producing heavily stretched tetrahedra, whereas our method evenly dissipates the stress over all the domain.

Free boundary injective mapping. To the best of our knowledge, our method is the only one passing the constrained-boundary benchmark, but able to produce inversion-free maps with free boundary. Since TLC tries to minimize the overall volume, relaxing the boundary constraints results in degenerate maps. Fig. 9 shows two maps obtained with our method: a 2D shape being compressed and the same shape being bent. The boundary is free to move, we lock the vertices shown in red. Refer to Fig. 10-a for the rest shape. The shape behaves exactly as a human would expect it: upon compression the shape chooses one of the two possible results (Fig. 9-a), and successfully passes the bend test (Fig. 9-b), note the geometrical details that are naturally rotated.

Shape-area trade-off 𝜆. Our final tests illustrate the influence of the parameter 𝜆 in Prob. (5) on the resulting map. We have computed three free boundary maps of the rest shape (Fig. 10-a) being stretched. First we chose 𝜆 = 0, that is, only the shape quality term is taken into account in Prob. ( 5). When we optimize for the angles, the area of the triangles is forced to change, refer to Fig. 10-b for the resulting map. Naturally, an area preserving map (𝜆 = 10 4 ) must deform the elements to satisfy the area constraint (Fig. 10-c). Finally, in Fig. 10-d we show an example with a trade-off between the area and angles preservation.

Map quality: comparison with LSCM. To assess quality of our maps, we computed a discrete conformal map of the "Lucy" model by setting 𝜆 = 0; we compare the result to least squares conformal map (LSCM) [START_REF] Lévy | Least Squares Conformal Maps for Automatic Texture Atlas Generation[END_REF]. LSCM is a very widespread method requiring to solve one linear system with a symmetric positive definite matrix. The idea is to compute a 𝑃 1 finite element approximation of the Cauchy-Riemann conditions over all triangles of the mesh. Numerical results are shown in Figure 11. As before, to compare quality of the maps, we use the condition number of the Jacobian matrix

𝜎 1 (𝐽 )
𝜎 2 (𝐽 ) , where 𝜎 1 and 𝜎 2 stand for the singular values of 𝐽 . In our map maximum condition number is equal to 5.1. The vast majority of mesh elements in LSCM method is mapped with very small shape deformation error, however the number of badly distorted elements (the distortion goes beyond 1000) is considerable and even some inverted elements are present.

Map quality: comparison with Simplex Assembly. Simplex Assembly (SA) [START_REF] Fu | Computing Inversion-Free Mappings by Simplex Assembly[END_REF]] is a method to compute inversion-free mappings with bounded distortion on simplicial meshes. The idea is to project each simplex into the inversion-free and distortionbounded space. Having disassembled the mesh, the simplices are then assembled by minimizing the mapping distortion, while keeping the mapping feasible.

SA is a quite robust method, nevertheless having some failure cases over the benchmark database. Fig. 12 provides a quality comparison of SA with our quasi-isometric (𝜆 = 1) map for a freeboundary mapping of the "Lucy" mesh. Note that for the 2D case SA energy is exactly the same as the one we use (for the 3D case the energy comes from [Knupp 2000b]), however SA includes an optimization for the distortion bounds, thus reaching better map quality.

Map quality: comparison with LBD. Large-scale Bounded Distortion Mappings (LBD) [START_REF] Shahar | Large-scale bounded distortion mappings[END_REF]] is another method to compute free-boundary mappings. Given an input map (potentially with folds), LBD looks for an injective map as close as possible to the input map, but satisfying some constraints such as the orientation and distortion bounds. our method. The 3D surface to flatten is a regular simplicial mesh of a rectangular patch that was lifted and noised. LBD has an optimization of the distortion bounds, thus the worst element quality of the map by LBD is better than in our map. Note however, that LBD has a lot of elements near the worst bound, whereas our methond is based on the elasticity theory, and providing better overall quality distribution.

Map quality: comparison with ABCD. We conclude the quality comparison by an example of 3D mesh deformation. Adaptive Block Coordinate Descent for Distortion Optimization (ABCD) [START_REF] Naitsat | Adaptive Block Coordinate Descent for Distortion Optimization[END_REF]]. Fig. 14 provides the quality plots. We took a tetrahedral mesh of a combination wrench, and we imposed positional constraints on the vertices located on both ends of the wrench. The deformation by ABCD is computed using ARAP energy, whereas our deformation is a quasi-isometry with 𝜆 = 1. In this experiment, our deformation has a slightly better, but overall comparable quality.

Limitations

While globally performing very well in practice, our method still presents some limitations. We have two main sources of limitations: theoretical limitations as well as very practical ones related to numerical stability of our resolution scheme.

Overlaps. First of all, an inversion-free map does not imply global injectivity. Fig. 15-a provides an example of an inversion-free map with two cases of non-injectivity when optimizing for a map with free boundaries: the map can present global overlaps as well as the boundary can "wind up" around boundary vertices, i.e. the total angle of triangles incident to a vertex can be superior to 2𝜋. Moreover, while being less frequent, similar situations may occur on interior vertices, as mentioned in [START_REF] Weber | Locally Injective Parametrization with Arbitrary Fixed Boundaries[END_REF]. In both cases, the map of the boundary is self-intersecting, as demonstrated in [START_REF] Aigerman | Injective and Bounded Distortion Mappings in 3D[END_REF] for the 3D case. Typically this situation happens near constraints causing a local compression in the shape.

Let us illustrate this behavior on a very simplistic mesh consisting of a single fan of 12 triangles. All vertices are free to move, the target shape is set to be the unit equilateral triangle for all elements. For this problem Fig. 16-a shows a local minimum, and the Fig. 16-c shows the global minimum respecting perfectly the prescribed total angle of 4𝜋 around the center vertex. Both are inversion-free maps, but only the map in Fig. 16-a is a globally injective one. Depending on the initialization and the resolution scheme chosen, we can converge to either minimum. Note, however, that the center vertex has the winding number 1 in one map and 2 in the other, and thus we can not deform continuously one to the other without inverting some elements. Note also that the configurations like in the Fig. 16-b present inverted elements and thus can not be generated by our method. Fig. 17 illustrates the problem on a larger mesh: we took a regular mesh of a square, swapped two vertices, and untangled the resulting mesh. Only two vertices are locked, the rest of the mesh is free to move. Fig. 17-a provides the rest shape, the vertices to be swapped are highlighted in red and blue, respectively. Depending on the minimization toolbox chosen, untangling may result in a mesh with a double covering present (Fig. 17-b) or produce the correct result, namely, a rigid transformation of the input mesh (Fig. 17-c).

It is possible to avoid all overlaps altogether by embedding our shape to optimize into an outer triangulation, and performing a "bi-material" optimization, as in [START_REF] Jiang | Simplicial Complex Augmentation Framework for Bijective Maps[END_REF]]. In this case, both global overlaps and fold-2-coverings are prohibited by the the fact that the outer material must not have inverted elements (refer to Fig. 15-b). The thick prismatic layer in Fig. 1 was generated by a similar procedure: we have generated a very thin layer of triangular prisms around the dragon, and tetrahedralized the exterior bounded by a cube. After calling the untangling procedure, we have obtained an offset surface with exactly the same mesh connectivity as the original dragon mesh.

While this embedding kind of approach works well for certain applications, for other it may be hard to apply. Another option is to add some extra triangles (resp. tets) in Eq. ( 5), but it is out of the scope of this article. Numerical challenges. Even when the problem is well-posed, a robust implementation may present significant difficulties. As we have said above, the quasi-Newtonian optimization scheme performs well for "simple" problems (it passes all the benchmark database!), but may fail for large deformations, where Newton iterations are necessary. While our modified Hessian matrix is symmetric positive definite, note that for stiff problems the Jacobi preconditioned conjugate gradients can fail and one might need the incomplete Choleski decomposition and beyond.

In practice we have found the method being very robust in 2D settings: we have not encountered a practical test case we were not able to treat with our method. In 3D, however, it can fail to the numerical challenges in very anisotropic and highly twisted meshes.

ANALYSIS

This section presents a rigorous analysis of the penalty method. We start with a discussion on invertibility conditions in the discrete and continous case ( § 4.1). Next, in § 4.2 we prove that the modified Hessian matrix 𝐻 + (𝑈 , 𝜀) is indeed positive definite, and finally in § 4.3 we show the origins of Eq. ( 6) for the regularization parameter sequence {𝜀 𝑘 } ... 𝑘=0 . Namely, we prove that if the problem has a solution, then an idealized minimization method can reach the admissible set min det 𝐽 > 0 in a finite number of steps. An immediate consequence of this theorem is that if the problem has a solution, then for some 𝐾 < ∞ the solution arg min 𝑈 𝐹 (𝑈 , 𝜀 𝐾 ) belongs to the admissible set.

Invertibility of the minimizer

A legitimate question arises: does the minimizer of our problem really give us an invertible mapping? This is not an easy question. On discrete level, the functional in Prob. ( 5) has an infinite barrier on the boundary of the set of admissible meshes. It means that any finite energy minimizer is free of inverted elements. Augmenting this barrier property with proper boundary conditions, say with prescribed boundary homeomorphism, one can prove that overall mesh deformation is a homeomorphism [START_REF] Aigerman | Injective and Bounded Distortion Mappings in 3D[END_REF][START_REF] Faivushevna | Problems of homeomorphism arising in the theory of grid generation[END_REF].

On the other hand, for the continuous settings the situation is much more subtle. Variational Prob. (1), while originating from the elasticity, violates the conditions formulated by Ball for its existence [START_REF] John | Convexity conditions and existence theorems in nonlinear elasticity[END_REF]] and invertibility [START_REF] Ball | Global invertibility of Sobolev functions and the interpenetration of matter[END_REF]] theorems. More precisely, the functional does not respect the growth conditions of the theorems. It does not mean that Prob. ( 1) is ill-posed, however we can say very little beyond that.

Why does it work so well in discrete settings? In fact, for every mesh deformation with finite energy the distortion measure for each cell is bounded from above. Hence, our numerical algorithm acts on the space of quasi-isometric homeomorphisms.

Consider following (continuous) problem related to construction of deformations with prescribed quality [START_REF] Va Garanzha | The barrier method for constructing quasi-isometric grids[END_REF]]:

arg min ì 𝑢 ( ì 𝑥) ∫ Ω 𝛽 (𝐽 ) 1 -𝑡 𝛽 (𝐽 ) 𝑑𝑥, 𝛽 (𝐽 ) = 𝑓 (𝐽 ) + 𝜆𝑔(𝐽 ) 𝑑 + 2𝜆 (7)
This integral can be finite only in the case 𝛽 (𝐽 ) < 1 𝑡 . Note that Prob. ( 1) is a special case of Prob. ( 7) for 𝑡 = 0, here parameter 𝑡 is the lower quality bound of the deformation. It is important to note that the density of the deformation energy is polyconvex and thus variational Prob. ( 7) is well-posed. Moreover, the existence theorem can be proved and under proper boundary conditions admissible deformations for Prob. ( 7) are quasi-isometric homemorphisms [START_REF] Garanzha | Variational method for untangling and optimization of spatial meshes[END_REF] both in 2D and 3D, i.e. the invertibility theorem [START_REF] Ball | Global invertibility of Sobolev functions and the interpenetration of matter[END_REF]] can be applied.

One may consider Prob. ( 5) as a minimization of functional ( 7) with arbitrarily small constant 𝑡 > 0. In its turn, Prob. ( 7) with near-zero value of 𝑡 can be seen as a regularized solution of (1).

As a side note, a transformation similar to ( 7) can be applied to a very broad class of polyconvex isochoric-volumetric split based energies which originally violate conditions of the Ball's invertibility theorem. This transformation suppresses singularities of deformations which is quite natural since singular deformations are beyond the scope of assumptions of elastic theory and should rely on other physical models.

Modified Hessian matrix

Recall that in our resolution scheme we use the modified Hessian (𝑑 #𝑉 ) × (𝑑 #𝑉 ) matrix 𝐻 + (𝑈 , 𝜀) of the function 𝐹 (𝑈 , 𝜀) built out of 𝑑 ×𝑑 blocks 𝐻 + 𝑖 𝑗 placed on the intersection of 𝑖-th block row and 𝑗-th block column. It is a common practice to add some regularization terms to the Hessian matrix to make it positive definite, but we propose to modify the finite element (FE) matrix assembly procedure by eliminating some terms potentially leading to an indefinite FE matrix.

To this end, we restrict our attention to a single simplex and we study a function 𝜙 (𝐽 ) of the Jacobian matrix defined as follows:

𝜙 (𝐽 ) := 𝑓 𝜀 (𝐽 ) + 𝜆𝑔 𝜀 (𝐽 ) = tr 𝐽 ⊤ 𝐽 (𝜒 (det 𝐽, 𝜀)) 2 𝑑 + 𝜆 det 2 𝐽 + 1 𝜒 (det 𝐽, 𝜀) (8)
Let us denote by 𝑎 ∈ R 𝑑 2 the (column-wise) flattening of the Jacobian matrix 𝐽 , i.e. the vector composed of the elements of 𝐽 . We decompose the 𝑑 2 × 𝑑 2 Hessian matrix of 𝜙 with respect to the Jacobian matrix entries into two parts:

𝜕 2 𝜙 𝜕 𝑎 𝜕 𝑎 ⊤ = 𝑀 + + 𝑀 ± , where 𝑀 + is a positive definite matrix, and 𝑀 ± can be an indefinite matrix Author version, 2021. that we neglect. The matrix 𝑀 ± contains all terms depending on 𝜒 ′′ and second derivatives of det 𝐽 with respect to elements of the Jacobian matrix 𝐽 . Our map is affine on the simplex of interest, therefore its Jacobian matrix 𝐽 is a linear function of the vertices of the simplex. The idea is to compute a positive definite matrix 𝑀 + (𝐽 ), and use the chain rule to get the Hessian matrix with respect to our variables 𝑈 and assemble the matrix 𝐻 + .

So, we choose some arbitrary point 𝐽 0 and we want to show the way to decompose 𝜕 2 𝜙 𝜕 𝑎 𝜕 𝑎 ⊤ (𝐽 0 ) into a sum of 𝑀 + (𝐽 0 ) and 𝑀 ± (𝐽 0 ) with 𝑀 + (𝐽 0 ) > 0. To do so, first we write down the first order Taylor expansion 𝑞(𝐷) of the function 𝜒 (𝐷, 𝜀) around some point 𝐷 0 = det 𝐽 0 :

𝑞(𝐷) := 𝜒 (𝐷 0 , 𝜀) + 𝜕 𝜒 𝜕 𝐷 (𝐷 0 , 𝜀)(𝐷 -𝐷 0 ).
Next we define a function Φ(𝑎, 𝐷) as follows:

Φ(𝑎, 𝐷) := |𝑎| 2 (𝑞(𝐷)) 2 𝑑 + 𝜆 𝐷 2 + 1 𝑞(𝐷)
.

Note that Φ differs a bit from 𝜙: it has one more argument and 𝜒 is replaced by its linearization in the denominator. While this maneuver might seem obscure, light will be shed very shortly. Φ has a major virtue of being convex! The convexity is easy to prove, refer to Appendix B for a formal proof.

Having built a convex function Φ, it is straightforward to verify that the following decomposition holds:

𝜕 2 𝜙 𝜕 𝑎 𝜕 𝑎 ⊤ (𝐽 0 ) = 𝑀 + (𝐽 0 ) + 𝑀 ± (𝐽 0 ), (9)
where

𝑀 + := 𝐼 𝜕 𝐷 𝜕 𝑎 𝜕 2 Φ 𝜕 𝑎 𝜕 𝑎 ⊤ 𝜕 2 Φ 𝜕 𝑎 𝜕 𝐷 𝜕 2 Φ 𝜕 𝐷 𝜕 𝑎 ⊤ 𝜕 2 Φ 𝜕 𝐷 2 𝐼 𝜕 𝐷 𝜕 𝑎 ⊤ , and
𝑀 ± := 𝜕 Φ 𝜕 𝐷 𝜕 2 𝐷 𝜕 𝑎 𝜕 𝑎 ⊤ - 𝜒 ′′ 𝜒 2 𝑑 𝑓 𝜀 + 𝜆𝑔 𝜀 𝜕 𝐷 𝜕 𝑎 𝜕 𝐷 𝜕 𝑎 ⊤ .
The easiest way to check that the equality (9) holds is to note that at the point 𝐽 0 we have 𝑞 = 𝜒, 𝑞 ′ = 𝜒 ′ , and therefore we have

𝜕 𝜙 (𝑎) 𝜕 𝑎 = 𝜕 Φ(𝑎, 𝐷 (𝑎)) 𝜕 𝑎 + 𝜕 Φ(𝑎, 𝐷 (𝑎)) 𝜕 𝐷 𝜕 𝐷 (𝑎) 𝜕 𝑎 .
To calculate the Hessian 𝜕 2 𝜙 𝜕 𝑎 𝜕 𝑎 ⊤ (𝐽 0 ), it suffices to differentiate this expression one more time and add the terms in 𝜒 ′′ that were zeroed out by the linearization.

To sum up, in our computations, for each simplex we approximate the Hessian matrix

𝜕 2 𝜙
𝜕 𝑎 𝜕 𝑎 ⊤ by the 𝑑 2 ×𝑑 2 matrix 𝑀 + and we neglect the term 𝑀 ± . Thanks to the convexity of Φ, it is trivial to verify that for any choice of 𝐽 0 the matrix 𝑀 + is positive definite. Then we use the chain rule over 𝑀 + to get the Hessian matrix with respect to our variables 𝑈 , and we assemble a (𝑑 #𝑉 ) × (𝑑 #𝑉 ) approximation 𝐻 + of the Hessian matrix for the energy function 𝐹 (𝑈 , 𝜀). Matrix 𝐻 + is positive definite provided that at least 𝑑 mesh vertices are fixed. If less than 𝑑 points are fixed, rigid body transformations are allowed. The energy is invariant w.r.t rigid body transformations, so when constraints allow for such transformations, matrix 𝐻 + becomes positive semi-definite. Note that the leading blocks 𝐻 + 𝑖𝑖 are always positive definite. Refer to Appendix A for further details on the finite element assembly procedure. Appendix C proves that positive definite part of the Hessian matrix is spectrally equivalent to finite element stiffness matrix for Laplace operator.

Choice of 𝜀 𝑘

In this section we provide a strategy for the choice of the regularization parameter 𝜀 𝑘 at each iteration. Namely, we prove that an idealized minimization algorithm reaches the admissible set min det 𝐽 > 0 in a finite number of iterations.

Theorem 1. Let us suppose that the admissible set is not empty, namely there exists a mesh 𝑈 * satisfying 𝐹 (𝑈 * , 0) < +∞. We also suppose that we have a minimization algorithm satisfying one of the following efficiency conditions for some 0 < 𝜎 < 1:

• either the essential descent condition holds

𝐹 (𝑈 𝑘+1 , 𝜀 𝑘 ) ≤ (1 -𝜎)𝐹 (𝑈 𝑘 , 𝜀 𝑘 ), (10)
• or the vector 𝑈 𝑘 satisfies the quasi-minimality condition:

min 𝑈 𝐹 (𝑈 , 𝜀 𝑘 ) > (1 -𝜎)𝐹 (𝑈 𝑘 , 𝜀 𝑘 ). (11)
Then the admissible set is reachable by solving a finite number of minimization problems in 𝑈 with 𝜀 𝑘 fixed for each problem. In other words, under a proper choice of the regularization parameter sequence 𝜀 𝑘 , 𝑘 = 0 . . . 𝐾, we obtain 𝐹 (𝑈 𝐾 , 0) < +∞.

Proof. The main idea is to expose an explicit way to build a decreasing sequence {𝜀 𝑘 } ∞ 𝑘=0 such that the sequence {𝐹 (𝑈 𝑘 , 𝜀 𝑘 )} ∞ 𝑘=0 is bounded from above. Then we can prove by contradiction that the admissible set is reachable in a finite number of steps, since if it is not, 𝐹 (𝑈 𝑘 , 𝜀 𝑘 ) must grow without bounds. First of all, the function 𝐹 (𝑈 , 𝜀) can be rewritten as follows

𝐹 (𝑈 , 𝜀) = 𝑖 𝛼 𝑖 𝜓 𝑖 (𝑈 ) 𝜒 (𝐷 𝑖 , 𝜀) , (12)
where 𝐷 𝑖 = 𝐷 𝑖 (𝑈 ) denotes the Jacobian determinant for 𝑖-th simplex of the mesh (𝐷 𝑖 = det 𝐽 𝑖 ), and 𝛼 𝑖 > 0 are positive, separated from zero weights assigned to each simplex. The functions 𝜓 𝑖 (𝑈 , 𝜀) := 𝜒 (𝐷 𝑖 , 𝜀) 1-2 𝑑 tr 𝐽 ⊤ 𝑖 𝐽 𝑖 + 𝜆(𝐷 2 𝑖 + 1) defined according to Eq. ( 3) are positive and bounded from below as

𝛼 𝑖 𝜓 𝑖 (𝑈 , 𝜀) ≥ 𝜆 min 𝑖 𝛼 𝑖 .
Note also that 𝜓 𝑖 (𝑈 , 𝜀) are increasing functions of 𝜀.

Our goal is to build a decreasing sequence {𝜀 𝑘 } ∞ 𝑘=0 such that the sequence {𝐹 (𝑈 𝑘 , 𝜀 𝑘 )} ∞ 𝑘=0 is bounded from above. We split the construction into two parts: first we suppose that at some iteration 𝑘 the essential condition (10) is satisfied, and then we explore the case (11).

Suppose that the condition (10) holds at iteration 𝑘. In order to guarantee that the function does not increase, it suffices to establish the following inequality:

(1 -𝜎)𝐹 (𝑈 𝑘+1 , 𝜀 𝑘+1 ) ≤ 𝐹 (𝑈 𝑘+1 , 𝜀 𝑘 ). (13)
By noting that𝜓 𝑖 (𝑈 𝑘+1 , 𝜀 𝑘+1 ) ≤ 𝜓 𝑖 (𝑈 𝑘+1 , 𝜀 𝑘 ), Ineq. ( 13) is implied if the following condition holds: ∀𝑖 :

(1 -𝜎) 𝜒 (𝐷 Let us show a constructive way to build 𝜀 𝑘+1 such that Ineq. ( 14) is satisfied. Assume that 𝜀 𝑘+1 is found as the solution to equation

𝜒 (𝐷 𝑘+1 -, 𝜀 𝑘+1 ) = (1 -𝜎) 𝜒 (𝐷 𝑘+1 -, 𝜀 𝑘 ), (15)
where 𝐷 𝑘+1 -:= min 𝑖 𝐷 𝑘+1 𝑖 . Obviously, if 𝐷 𝑘+1 -< (1 -𝜎) 𝜒 (𝐷 𝑘+1 -, 𝜀 𝑘 ), there is a unique positive solution to this equation. Now let us evaluate derivative of the function 𝑠 (𝐷) := 𝜒 (𝐷,𝑡 1 ) 𝜒 (𝐷,𝑡 2 ) , 𝑡 1 < 𝑡 2 . It is straightforward to verify that the derivative 𝑠 ′ (𝐷) > 0:

𝑠 ′ (𝐷) = 𝜒 ′ (𝐷, 𝑡 1 ) 𝜒 (𝐷, 𝑡 1 ) - 𝜒 ′ (𝐷, 𝑡 2 ) 𝜒 (𝐷, 𝑡 2 ) 𝑠 = 1 𝐷 2 + 𝑡 2 1 - 1 𝐷 2 + 𝑡 2 2 𝑠 > 0.
This inequality means that Eq. ( 15) for all 𝐷 𝑘+1 𝑖 ≥ 𝐷 𝑘+1 -induces inequality

𝜒 (𝐷 𝑘+1 𝑖 , 𝜀 𝑘+1 ) 𝜒 (𝐷 𝑘+1 𝑖 , 𝜀 𝑘 ) ≥ 𝜒 (𝐷 𝑘+1 -, 𝜀 𝑘+1 ) 𝜒 (𝐷 𝑘+1 -, 𝜀 𝑘 ) = 1 -𝜎,
which is precisely Ineq. ( 14). Hence, if 𝑈 𝑘+1 is an approximate solution of the minimization problem arg min 𝑈 𝐹 (𝑈 , 𝜀 𝑘 ) with fixed parameter 𝜀 𝑘 , we can use the following update rule for 𝜀 𝑘+1 which is the explicit solution of equation Eq. ( 15):

𝜀 𝑘+1 = 2 𝜇 𝑘 (𝜇 𝑘 -𝐷 𝑘+1 -), (16)
where

𝜇 𝑘 := (1 -𝜎) 𝜒 (𝐷 𝑘+1 -, 𝜀 𝑘 ) (17)
This update rule guarantees that Ineq. ( 13) is satisfied; coupled with the assumption (10) of the theorem, this implies the required nongrowth property of the function values sequence: 𝐹 (𝑈 𝑘+1 , 𝜀 𝑘+1 ) ≤ 𝐹 (𝑈 𝑘 , 𝜀 𝑘 ).

(18)

Consider now the case where condition (11) holds at iteration 𝑘. Note that condition (11) essentially means that our current solution 𝑈 𝑘 is very close to the global minimum of 𝐹 (𝑈 , 𝜀 𝑘 ), and thus Ineq. (10) cannot be satisfied. Nevertheless, we can use the same update rule (16) for computation of 𝜀 𝑘+1 . Indeed, with this choice we have

𝐹 (𝑈 𝑘+1 , 𝜀 𝑘+1 ) ≤ 1 (1 -𝜎) 𝐹 (𝑈 𝑘+1 , 𝜀 𝑘 ) ≤ 1 (1 -𝜎) 𝐹 (𝑈 𝑘 , 𝜀 𝑘 ) < 1 (1 -𝜎) 2 min 𝑈 𝐹 (𝑈 , 𝜀 𝑘 ) < 1 (1 -𝜎) 2 min 𝑈 𝐹 (𝑈 , 0).
Here the last inequality provides a global bound on the function values sequence, and it is based on the observation 𝜕 𝜕𝜀 𝜒 (𝐷, 𝜀) > 0.

To sum up, we have shown a way to build a sequence {𝜀 𝑘 } ∞ 𝑘=0 such that the sequence {𝐹 (𝑈 𝑘 , 𝜀 𝑘 )} ∞ 𝑘=0 is bounded from above. Now let us prove that the update rule ( 16), (17) allows to reach the admissible set in a finite number of steps. To do so, we use a reductio ad absurdum argument.

Suppose that the admissible set is never reached for an infinite decreasing sequence {𝜀 𝑘 } ∞ 𝑘=0 built using the update rule ( 16), (17) i.e. 𝐷 𝑘+1 -< 0 ∀𝑘 ≥ 0.

One can readily see that the following identity can be deduced from ( 16), (17), and (2):

(𝜀 𝑘+1 ) 2 = (1 -𝜎) (𝜀 𝑘 ) 2 -4𝜎 (𝜒 (𝐷 𝑘+1 -, 𝜀 𝑘 )) 2 . ( 19
)
In particular, (19) shows the strict decrease of the sequence 𝜀 𝑘 . Moreover, from ( 19) one obviously has

(𝜀 𝑘 ) 2 ≥ 4𝜎 (𝜒 (𝐷 𝑘+1 -, 𝜀 𝑘 )) 2 ,
and using the latter with (19) gives the following inequality:

(𝜀 𝑘 ) 2 -(𝜀 𝑘+1 ) 2 = 𝜎 (𝜀 𝑘 ) 2 + 4𝜎 (1 -𝜎)( 𝜒 (𝐷 𝑘+1 -, 𝜀 𝑘 )) 2 ≥ 4𝜎 2 (𝜒 (𝐷 𝑘+1 -, 𝜀 𝑘 )) 2 + 4𝜎 (1 -𝜎)( 𝜒 (𝐷 𝑘+1 -, 𝜀 𝑘 )) 2 = 4𝜎 ( 𝜒 (𝐷 𝑘+1 -, 𝜀 𝑘 )) 2 .
Hence for some 𝐾 > 0:

(𝜀 0 ) 2 -(𝜀 𝐾 ) 2 ≥ 4𝜎 𝐾-1 𝑘=0 (𝜒 (𝐷 𝑘+1 -, 𝜀 𝑘 )) 2 ≥ 4𝜎𝐾 min 0≤𝑘<𝐾 (𝜒 (𝐷 𝑘+1 -, 𝜀 𝑘 )) 2 ,
with an immediate consequence that for an arbitrarily large 𝐾 we have the following inequality:

max 0≤𝑘<𝐾 1 𝜒 (𝐷 𝑘+1 -, 𝜀 𝑘 ) ≥ √ 4𝜎𝐾 𝜀 0 .
Since all terms 𝛼 𝑖 𝜓 𝑖 (𝑈 ) in ( 12) are bounded from below, the resulting estimate contradicts the boundedness of 𝐹 (𝑈 𝑘 , 𝜀 𝑘 ), thus concluding our proof. □ Remark 1. An important corollary of Th. 1 is that, provided that the admissible set is not empty, there exists an iteration 𝐾 < ∞ such that the global minimum of the function 𝐹 (𝑈 , 𝜀 𝐾 ) belongs to the admissible set. The proof is rather obvious: suppose we have an idealized minimizer such that 𝑈 𝑘+1 = arg min 𝑈 𝐹 (𝑈 , 𝜀 𝑘 ). This minimizer always satisfies the conditions of Th. 1, therefore it can untangle the mesh in a finite number of steps.

Remark 2. In practice the global estimate 𝜎 is not known in advance, and the optimization routine may be far from the ideal. For each minimization step we compute the local descent coefficient 𝜎 𝑘 : (11), it means that minimization procedure for 𝜀 𝑘 failed and theorem cannot be applied (it does not mean that Alg. 1 will not reach the admissible set!). In numerical experiments we use value 𝜎 = 1 10 .

𝜎 𝑘 := 1 - 𝐹 (𝑈

CONCLUSION

Producing maps without inverted elements is a challenge in geometry processing. Inspired by untangling solutions in computational physics, our solution outperforms the state of the art in terms of robustness. It is easy to use since we provide a simple implementation that is free of commercial product dependency (compiler, library, etc.). Moreover, the energy is estimated independently on Author version, 2021.

each triangle / tetrahedra, making it a good candidate to be adapted to more difficult settings including free boundary and global parameterization. and

A COMPREHENSIVE DESIGN FORMULAE

𝑎 ⊤ = ( ì 𝑎 ⊤ 1 . . . ì 𝑎 ⊤ 𝑑 ), 𝑏 ⊤ = ( ì 𝑏 ⊤ 1 . . . ì 𝑏 ⊤ 𝑑 ).

A.1 Gradient

In order to derive expressions for the gradient and the Hessian matrix of 𝐹 , we write down explicitly the Jacobian matrix 𝐽 for the affine map of a simplex 𝑇 with vertices ì 𝑢 0 , ì 𝑢 1 , . . . , ì 𝑢 𝑑 :

𝐽 = ( ì 𝑎 1 . . . ì 𝑎 𝑑 ) = ( ì 𝑢 1 -ì 𝑢 0 ì 𝑢 2 -ì 𝑢 0 . . . ì 𝑢 𝑑 -ì 𝑢 0 ) 𝑆 -1 = = ( ì 𝑢 0 . . . ì 𝑢 𝑑 )𝑍,
where

𝑆 := ( ì 𝑥 1 -ì 𝑥 0 ì 𝑥 2 -ì 𝑥 0 . . . ì 𝑥 𝑑 -ì 𝑥 0 ), det 𝑆 > 0
is the shape matrix, ì 𝑥 𝑖 are vertices of "ideal" or "target" shape for the image of the simplex 𝑇 , and 𝑍 is a (𝑑 + 1) × 𝑑 matrix defined as

𝑍 := {𝑧 𝑖 𝑗 } := -1 . . . -1 𝐼 𝑆 -1
Since the Jacobian matrix is a linear function of ì 𝑢 𝑖 , we have

𝜕 ì 𝑎 𝑖 𝜕 ì 𝑢 ⊤ 𝑗 = 𝑧 𝑗𝑖 𝐼, 𝑖 = 1, . . . , 𝑑, 𝑗 = 0, . . . , 𝑑.
The additive contribution to gradient of 𝐹 from the simplex 𝑇 can be written using correspondence of local indices 0 -𝑑 and global indices 𝑔 0 -𝑔 𝑑 in the list of vertices:

(∇𝐹 ) 𝑔 𝑗 += det 𝑆 𝑑! 𝑑 𝑖=1 𝜕 ì 𝑎 ⊤ 𝑖 𝜕 ì 𝑢 𝑗 𝜕 𝜙 𝜕 ì 𝑎 𝑖 = = det 𝑆 𝑑! 𝑑 𝑖=1 𝑧 𝑗𝑖 𝜕 𝜙 𝜕 ì 𝑎 𝑖 , 𝑗 = 0, . . . ,

A.2 Hessian

The blocks of the non-negative definite part of Hessian matrix of 𝐹 can be updated using the following general formula

𝐻 + 𝑔 𝑗 𝑔 𝑖 += det 𝑆 𝑑! 𝑚,𝑙 𝜕 ì 𝑎 ⊤ 𝑚 𝜕 ì 𝑢 𝑗 𝑀 + 𝑚𝑙 𝜕 ì 𝑎 𝑙 𝜕 ì 𝑢 ⊤ 𝑖 ,
where 𝑀 + 𝑚𝑙 denotes a 𝑑 × 𝑑 block of 𝑑 2 × 𝑑 2 positive definite matrix 𝑀 + defined in Eq. ( 9). Let us provide an explicit expression for the matrix:

𝑀 + = 𝐼 𝑏 𝜕 2 Φ 𝜕 𝑎 𝜕 𝑎 ⊤ 𝜕 2 Φ 𝜕 𝑎 𝜕 𝐷 𝜕 2 Φ 𝜕 𝐷 𝜕 𝑎 ⊤ 𝜕 2 Φ 𝜕 𝐷 2 𝐼 𝑏 ⊤ , where 𝜕 2 Φ 𝜕 𝑎 𝜕 𝑎 ⊤ = 2 𝜒 2 𝑑 𝐼 𝜕 2 Φ 𝜕 𝐷 2 = 2 𝑑 1 + 2 𝑑 |𝑎| 2 𝜒 ′2 𝜒 2+ 2 𝑑 + 𝜆 2 𝜒 -4𝐷 𝜒 ′ 𝜒 2 + 2(1 + 𝐷 2 ) 𝜒 ′2 𝜒 3 𝜕 2 Φ 𝜕 𝑎 𝜕 𝐷 = - 4 𝑑 𝜒 ′ 𝜒 1+ 2 𝑑 𝑎.
Author version, 2021.

Obviously the leading 𝑑 × 𝑑 blocks of the matrix 𝐻 + are strictly positive definite and can be used to build Newton-type minimization algorithm.

B CONVEXITY OF Φ

Recall that the function Φ is defined as follows ( § 4.2):

Φ(𝑎, 𝐷) := |𝑎| 2 𝑞 2 𝑑 + 𝜆 𝐷 2 + 1 𝑞 ,
where 𝑎 ∈ R 𝑑 2 is the (column-wise) flattening of the Jacobian matrix 𝐽 , i.e. the vector composed of the elements of 𝐽 . Lemma 1. ∇∇ ⊤ Φ > 0 Proof. It is straightforward to see that the (𝑑 2 + 1) × (𝑑 2 + 1) Hessian matrix of Φ can be written in the 2 × 2 block representation:

∇∇ ⊤ Φ = 𝑃 + 𝜆𝑄, where 𝑃 := 2 𝑞 2 𝑑 𝐼 -4 𝑑 𝑞 ′ 𝑎 𝑞 1+ 2 𝑑 -4 𝑑 𝑞 ′ 𝑎 ⊤ 𝑞 1+ 2 𝑑 2 𝑑 (1 + 2 𝑑 ) |𝑎 | 2 𝑞 ′2 𝑞 2+ 2 𝑑 and 𝑄 := 0 0 0 2 𝑞 -4𝐷 𝑞 ′ 𝑞 2 + 2(1 + 𝐷 2 ) 𝑞 ′2 𝑞 3
.

It is trivial to verify that 𝑄 ≥ 0, since 𝑄 22 is a strictly positive quadratic function of argument 𝐷. Since the leading blocks of the matrix 𝑃 are positive definite and the Schur complement

𝑃 22 -𝑃 21 𝑃 -1 11 𝑃 12 = |𝑎| 2 𝑞 2+ 2 𝑑 2 𝑑 1 - 2 𝑑 ≥ 0
is non-negative definite, overall convexity of Φ is established. □ Remark 3. Note that we have just proved the convexity of the function Φ. As an immediate consequence we obtain polyconvexity of the functional (1), because it is a particular case of our functional (4) with 𝜒 = 𝑞.

C SPECTRAL BOUNDS FOR THE HESSIAN

Using assumption that function 𝐹 (𝑈 , 𝜀) is bounded, let us derive (non-tight) spectral bounds for the positive definite part of the Hessian matrix at the point 𝐽 = 𝐽 0 , 𝐷 = 𝐷 0 , 𝑞 = 𝜒 (𝐽 0 , 𝜀). Evidently contribution from each simplex is bounded 𝜙 (𝐽 0 ) < 𝐾, meaning that 1 + 𝐷 2 < 𝐾 𝜆 𝜒 (𝐷, 𝜀), |𝑎| 2 < 𝐾 ( 𝜒 (𝐷, 𝜀)) 2/𝑑 From these inequalities and the fact that 𝜆 𝐾 < 𝑞 < √ 𝐷 2 + 𝜀 2 , we can deduce The determinant can be bounded as

𝐷 2 + 𝜀 2 > 𝜆 2 𝐾 2 , 𝐷 2 < 𝐾 2 𝜆 2 ,
det(𝑃 + 𝜆𝑄) = 2 𝑑 2 𝑞 2𝑑 𝜆 2 𝑞 -4𝐷 𝑞 ′ 𝑞 2 + 2(1 + 𝐷 2 ) 𝑞 ′2 𝑞 3 ≥ 2 𝑑 2 +1 𝑞 2𝑑+1 𝜆 1 + 𝐷 2 ,
where the last inequality is obtained by taking the minimum over (𝑞 ′ /𝑞) > 0 considered as independent variable.

The spectral estimates for matrix 𝑃 + 𝜆𝑄 can be expressed via 𝐾, 𝜆, 𝜀 and are guaranteed to be uniformly bounded from below and from above, provided that 𝜆 > 0 and 𝜀 is bounded from above. Since vector 𝑏 is uniformly bounded from above, we immediately get uniform bounds 𝑘 1 𝐼 < 𝑀 + < 𝑘 2 𝐼, 𝑘 1 < 𝑘 2 where parameters 𝑘 𝑖 = 𝑘 𝑖 (𝐾, 𝜆, 𝜀) > 0 are uniformly bounded from above and from below. It means that Q = [ np.matrix('-1,-1;1,0;0,0;0,1'), np.matrix('-1,0;1,-1;0,1;0,0'), # quadratures for 8 np.matrix('0,0;0,-1;1,1;-1,0'), np.matrix('0,-1;0,0;1,0;-1,1') ] # every quad corner

Fig. 2 .

 2 Fig. 2. Regularization function for the denominator in Eq. (3). When 𝜀 tends to zero, 𝜒 (𝜀, 𝐷) tends to 𝐷 for positive values of 𝐷, and to 0 + for negative values of 𝐷.

ALGORITHM 1 :

 1 Fig.3. The input mesh with foldovers and the untangling produced by the Listings 1 and 2. Left: a quad 2D mesh to untangle. The boundary (in red) is locked, and the black mesh is free to move. Right: fold-free result.

Fig. 4 .

 4 Fig. 4. Injective mapping sanity check: swap any two points inside a square. Left: the input problem, all locked points are shown in red. Right: foldoverfree result obtained with our method.

Fig. 5 .

 5 Fig. 5. Constrained boundary injective mapping challenge proposed by[START_REF] Du | Lifting Simplices to Find Injectivity[END_REF]: the "Lucy" mesh is mapped to a P-shaped domain by constraining the vertices shown in red. Left column: three different initializations for the same problem. Middle column: our method produces the same (up to a numerical precision) result on all three initializations, we chose 𝜆 = 0.01. Right column: total lifted content method[START_REF] Du | Lifting Simplices to Find Injectivity[END_REF]] fails to solve for the randomly initialized interior vertices, and produces very different results on other two initializations. Three miniature images of "Lucy" show in red the portion of the surface visible in the corresponding close-ups.

Fig. 6 .

 6 Fig. 6. Constrained boundary injective mapping stress test. We have generated an isotropic tetrahedral mesh of a cube subtracted from a larger cube, and we rotate the inner cube's boundary to test the robustness. Top row and bottom row correspond to two different slices of the same mesh. Columns (a)-(d):injective maps produced by our method, columns (e) and (f): injective maps produced by[START_REF] Du | Lifting Simplices to Find Injectivity[END_REF]]. The method by[START_REF] Du | Lifting Simplices to Find Injectivity[END_REF]] fails to generate injective maps for 135°and 180°inner cube rotations. The color map illustrates the relative volume scaling: green for det 𝐽 ≈ 1, red for inflation, blue for compression.

Fig. 7 .

 7 Fig. 7. Quality plot of the resulting locally injective maps for every challenge from the database provided by [Du et al. 2020]. Our results are shown in blue, whereas the results by Du et al. are shown in red. Each dot corresponds to a quality of the corresponding map reduced to two numbers: the maximum stretch and the minimum scale. Top row: mapping quality on the 2D dataset. Bottom row: mapping quality on the 3D dataset. Left column shows the absolute maximum stretch and absolute minimum scale, whereas the right column shows the maximum stretch and minimum scale for the top 95% of measurements.

Fig. 8 .

 8 Fig.8. Performance of our method tested on the benchmark[START_REF] Du | Lifting Simplices to Find Injectivity[END_REF]. Each dot corresponds to a challenge from the database (10743 in 2D and 904 in 3D). Blue dots show the running times obtained using a heuristic regularization[Garanzha and Kaporin 1999, Eq. 6.3], green dots correspond to Eq. (6).

Fig. 9 .Fig. 10 .

 910 Fig. 9. Free boundary injective mapping. The vertices shown in red are constrained, all other vertices are free to move. (a): a compression test, (b): a bend test. Refer to Fig 10-a for the rest shape.

Fig. 11 .

 11 Fig. 11. Comparison of two discrete conformal maps for the "Lucy" mesh. (a): Least-squares conformal map, (b): our map obtained with 𝜆 = 0. Top row: flattenings, bottom row: log-log histograms of conformal maps element quality.

Fig. 13 .

 13 Fig. 12. Two quasi-isometric maps for the "Lucy" mesh. (a): Simplex assembly, (b): our map obtained with 𝜆 = 1. Top row: flattenings, middle and bottom rows: log-log histograms of element quality.

Fig. 14 .

 14 Fig. 14. Tetrahedral mesh deformation, locked vertices are shown in red. (a): Rest shape, (b): ABCD, (c): our result with 𝜆 = 1. Top row: flattenings, middle and bottom rows: log-log histograms of element quality.

Fig. 15 .Fig. 16 .

 1516 Fig. 15. Free boundary mapping limitations. (a): this mesh presents two kinds of problems, namely, a global overlap and the mesh wrapped around a boundary vertex. (b): Both problems can be avoided by embedding the mesh into an outer triangulation.

Fig. 17 .

 17 Fig. 17. Loss of invertibility under a free boundary deformation: 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.

Fig. 18 .

 18 Fig. 18. On each simplex the map ì 𝑢 ( ì 𝑥) is affine and is entirely defined by the position of the vertices of the domain simplex { ì 𝑥 𝑖 } and its image { ì 𝑢 𝑖 }.

𝑘 1

 1 D ℎ (𝑈 ) < 1 2 𝑈 ⊤ 𝐻 + 𝑈 < 𝑘 2 D ℎ (𝑈 ),where D ℎ (𝑈 ) is the discrete Dirichlet functional for standard simplicial linear finite elements which approximates Dirichlet functional Thus we have demonstrated stability of the positive definite part of the Hessian matrix near the barrier.Author version, 2021. 102:16 • Vladimir Garanzha, Igor Kaporin, Liudmila Kudryavtseva, François Protais, Nicolas Ray, and Dmitry Sokolov Listing 1. A complete L-BFGS-based quad mesh untangling example, the result is shown in Fig.3-right.

G

  qc, quad): # evaluate the Jacobian matrix at the given quadrature point 11 return np.matrix([[U[quad[0] ], U[quad[1] ], U[quad[2] ], U[quad[3] ]],12 [U[quad[0]+n], U[quad[1]+n], U[quad[2]+n], U[quad[3]+n]]]) * Q[qc]13 14 for iter in range(10): # outer L-BFGS loop 15 mindet = min( [ np.linalg.det( jacobian(mesh.x, qc, quad) ) for quad in mesh.quads for qc in range(4) ] ) 16 eps = np.sqrt(1e-6**2 + .04*min(mindet, 0)**2); # the regularization parameter 𝜀 17 18 def energy(U): # compute the energy and its gradient for the map ì 𝑢 19 F,G = 0, np.zeros(2*n) 20 for quad in mesh.quads: # sum over all quads 21 for qc in range(4): # evaluate the Jacobian matrix for every quad corner 22 J = jacobian(U, qc, quad) 23 det = np.linalg.det(J) 24 chi = det/2 + np.sqrt(eps**2 + det**2)/2 # the penalty function 𝜒 (𝜀, det(𝐽 )) 25 chip = .5 + det/(2*np.sqrt(eps**2 + det**2)) # its derivative 𝜒 ′ (𝜀, det(𝐽 )) 26 f = np.trace(np.transpose(J)*J)/chi # quad corner shape quality 27 qc] * np.transpose(dfdj) # chain rule for the real variables 30 for i,v in enumerate(quad): 31 if (mesh.boundary[v]): continue # the boundary verts are locked 32 fmin_l_bfgs_b(energy, mesh.x)[0] # inner L-BFGS loop 36 print(mesh) # print wavefront .obj file Listing 2. A simplistic quad mesh class, the initialization is shown in Fig. 3self): # generate the test problem: a regular 2d grid with upper half shifted 5 n = 8 6 self.x = [ i/n + int(j>=n//2)*3/5 for j in range(n) for i in range(n) ] + \ 7 [ 2*j/n -int(j>=n//2)*3/5 for j in range(n) for i in range(n) ] # 2D geometry 8 self.quads = [ [i+j*n, i+1+j*n, i+1+(j+1)*n, i+(j+1)*n] for j in range(n-1) for i in range(n-1) ] # connectivity 9 self.boundary = [ i==0 or i==n-1 or j==0 or j==n-1 for j in range(n) for i in range(n) ] # vertex boundary flags range(self.nverts):

  𝐷 𝑖 (𝑈 𝑘+1 ) denotes the Jacobian determinant of simplex 𝑖 at iteration 𝑘 + 1.

	102:12 • Vladimir Garanzha, Igor Kaporin, Liudmila Kudryavtseva, François Protais, Nicolas Ray, and Dmitry Sokolov
	𝑖 where 𝐷 𝑘+1	:=
		𝑘+1 𝑖 , 𝜀 𝑘 ) ≤ 𝜒 (𝐷 𝑘+1 𝑖 , 𝜀 𝑘+1 )	(14)
		Author version, 2021.

  𝑘+1 , 𝜀 𝑘 ) 𝐹 (𝑈 𝑘 , 𝜀 𝑘 ) . When 𝜎 𝑘 ≥ 𝜎 one can use the update rule (16), (17) using the local value 𝜎 𝑘 guaranteeing that Ineq. (18) holds. In the case 𝜎 𝑘 < 𝜎 one should check that condition (11) holds for prescribed 𝜎. If positive, we can assign 𝜎 𝑘 = 𝜎 and use update rule (16), (17). If one cannot assure

  Given a map ì 𝑢, let us denote by ì 𝑎 𝑖 , 𝑖 = 1, 2 (, 3) the tangent basis, i.e. vectors forming the columns of the Jacobian matrix 𝐽 . For ex-

	ample, in 2D we have ì 𝑎 1 := 𝜕 𝑢 𝜕 𝑥	𝜕 𝑣 𝜕 𝑥	⊤	and ì 𝑎 2 := 𝜕 𝑢 𝜕 𝑦	𝜕 𝑣 𝜕 𝑦	⊤	.
	Let us denote by ì 𝑏 𝑖 the dual basis, i.e. vectors chosen in the way
	that ì 𝑎 ⊤ 𝑖 settings the dual basis can be written as ì ì 𝑏 𝑗 = 𝛿 𝑖 𝑗 det 𝐽 for all indices 𝑖, 𝑗. In particular, for the 2D 𝑏 1 := 𝜕 𝑣 𝜕 𝑦 ⊤ -𝜕 𝑢 𝜕 𝑦 and
	ì 𝑏 2 := -𝜕 𝑣 𝜕 𝑥	𝜕 𝑢 𝜕 𝑥	⊤	. In the 3D case ì 𝑏 𝑘 = ì 𝑎 𝑖 × ì 𝑎 𝑗 , where 𝑖, 𝑗, 𝑘 is
	cyclic permutation from 1, 2, 3. It is a handy choice of variables, in
	particular, tr 𝐽 ⊤ 𝐽 = 𝑖 | ì 𝑎 𝑖 | 2 and 𝜕 det 𝐽 𝜕 ì 𝑎 𝑖	= ì 𝑏 𝑖 . For further simplifi-

cation of notations we will use 𝜒 for 𝜒 (𝐷, 𝜀), 𝜒 ′ for 𝜕 𝜒 (𝐷,𝜀)

𝜕 𝐷

  𝑞 < 𝐾 2 𝜆 2 + 𝜀 2 . 𝑑 |𝑎| 𝑑 < 𝑑𝐾 𝑑 𝐾 2 𝜆 2 + 𝜀 2 .The last inequality follows from the fact that 𝑏 consists of columns of the matrix cof 𝐽 . Now we can estimate spectral bounds for the matrix 𝑃 + 𝜆𝑄. Clearly,𝜆 max (𝑃 + 𝜆𝑄) ≤ tr(𝑃 + 𝜆𝑄) 𝐷 2 + 𝜀 2 )were used to obtain the last equality.A lower bound for the minimum eigenvalue 𝜆 1 follows from the simple estimate obtained from the arithmetic-geometric mean inequality written for the eigenvalues 𝜆 2 ≤ . . . ≤ 𝜆 𝑑 2 +1 of 𝑃 + 𝜆𝑄 : 0 < 𝜆 min (𝑃 + 𝜆𝑄) = 𝜆 1 = det(𝑃 + 𝜆𝑄)

			=	2𝑑 2 𝑞 2 𝑑	+	2 𝑑	(1 +	2 𝑑	)	|𝑎| 2 𝑞 ′2 𝑞 2+ 2 𝑑	+ 𝜆	2 𝑞	-4𝐷	𝑞 ′ 𝑞 2 + 2(1 + 𝐷 2 )	𝑞 ′2 𝑞 3
				=	2𝑑 2 𝑞 2 𝑑	+	2 𝑑	(1 +	2 𝑑	)	|𝑎| 2 𝑞 ′2 𝑞 2+ 2 𝑑	+	2𝜆 𝑞	1 𝐷 2 + 𝜀 2 1 +	𝜀 4 4𝑞 2 ,
			where the relations
					𝑞 ′ 𝑞	=	1 𝐷 2 + 𝜀 2 √	,	𝑞 =	1 2	(𝐷 + 𝑑 2 +1 𝑘=2 𝜆 𝑘
				≥	det(𝑃 + 𝜆𝑄) 1 𝑑 2 𝑘=2 𝜆 𝑘 𝑑 2 +1	𝑑 2 >	𝑑 2 tr(𝑃 + 𝜆𝑄)	𝑑 2	det(𝑃 + 𝜆𝑄).
	Hence we immediately obtain
	|𝑎| 2 < 𝐾	𝜆 2 + 𝜀 2 𝐾 2	1/𝑑

, |𝑏 | <

The minimization algorithm is subject to conditions of Th. (1); from our numerical experiments we observe that our minimization algorithm almost always respects the conditions.

Attention, this step assumes that the solution of Prob. (1) is a diffeomorphism.

The fact that the solution of Prob. (4) is a diffeomorphism is sufficient (but not necessary) for the equivalence.

Author version, 2021.

* This work is supported by the Ministry of Science and Higher Education of the Russian Federation, project No 075-15-2020-

Author version, 2021.

Author version, 2021.