Damien Clergeaud

François Guillaume

Pascal Guitton

3D Collaborative Interaction for Aerospace industry

3D collaborative interaction raises issues like awareness of the Virtual Environment that have been well known for some time. They have been the subject of research activities that have led to some interesting results described in the literature. Unfortunately, most of these results have not yet been not implemented in the software systems daily used by industry and remains only theoretical concepts. That is why we recently started a project to transfer some of these results in the aerospace industry for the Airbus Group. Beyond this first transfer target, we also intend to measure the real gains (in real industrial conditions) for the users, and then for the company. This second goal is essential because in most publications, user testing is not satisfying (lack of real users, lack of real procedures, non significant number of subjects). In this paper, we describe the first step of this work in progress and more precisely, the basic interaction features we have developed. Currently, we are designing the first user tests.

PRESENTATION OF THE OVERALL CONTEXT

Virtual reality has been entering the industrial world for a long time. It is now used routinely in design and manufacturing processes mainly in large companies. In this project, we consider a specific industrial context: aerospace with Airbus Group. Airbus Group have been using VR for more than ten years to design and develop their aerospace products. A second major step will be to perform more complex experiments simultaneously involving several remote sites. For instance, if a complex problem occurs during the final stages of assembling a launcher, it may be necessary to involve experts who were involved in the initial design and/or the initial assembly. Traditionally these meetings are held over the telephone or using videoconferencing. Airbus would like to replace them with collaborative VR working sessions between multiple remote sites. It is well known that such sessions raise complex issues particularly for interaction between users. The first is the awareness of the Virtual Environment: when a user is immersed, it may be very difficult to fully perceive -and then, understandthe VE. Of course, the more complex the VE is, the more difficult it is to understand. So, improving awareness is still a big challenge for collaborative interaction in VR and several aspects must be considered [START_REF] Nguyen | A Survey of Communication and Awareness in Collaborative Virtual Environments[END_REF]:

-Components of the VE: it is not easy to fully perceive a dense environment like the inside of a launcher (figure 1). There are lots of technical components (pipes, harnesses equipement…) which may induce occlusions that prevent the user from seeing the whole VE at a single glance. Moreover, it is easy to get disoriented, especially, when complex operations which induce forced body movements are required. Therefore, it is necessary to provide users with efficient spatial referencing tools so that they know precisely where they are in the VE [START_REF] Stoakley | Virtual Reality on a WIM: Interactive Worlds in Miniature[END_REF].

-Other immersed Users: if multiple users are simultaneously immersed, it may be difficult for one user to perceive the location and/or the actions of the other users. Unfortunately, some solutions to these problems can be antagonistic. For instance, adding additional viewing windows to address the problem of remote visualization may decrease the feeling of presence and therefore awareness. Consequently, trying to perform efficient collaborative working sessions is much more complex than simply finding three separate solutions to these three issues. It is therefore necessary to establish balance between sometimes opposing constraints.

In order to provide efficient solutions to this global challenge, we first observed the current industrial procedures in VR experiments. Then we started to develop some basic collaborative interaction features based on previous work. These first features have already been tested and validated by Airbus experts. The second step is to integrate them into the operational software systems used by Airbus. The third step is to test these features in real conditions: real users and procedures, with a significant number of users.

DESCRIPTION OF THE CURRENT SYSTEM

We are designing a system which offers interaction tools that improve awareness for immersed users. The purpose of these tools * Damien.Clergeaud@inria.fr, Francois.Guillaume@airbus.com, Pascal.Guitton@inria.fr is to adress issues that often occur during real VR collaborative sessions : Where are you ? What do you see ? Which component are you talking about ? etc. In current sessions, we observed that these exchanges (questions, answers) are almost exclusiveley oral. They generate noise and interference and cause a lot of repetition between the users. The system we are developing is intended to reduce/eliminate these drawbacks by replacing most of the oral exchanges by VR interaction tools. It seems that it is easier to describe a position in the VE or to designate a 3D object by showing the information in the scene through an appropriate metaphor.

Basic Designation Features

In the VE, users are represented by a virtual head and a virtual hand. In order to ease the manipulation and to access a virtual menu, they can switch their virtual hand to a virtual ray.

To show a specific component, users can either use a virtual ray [START_REF] Hrimech | How 3D Interaction Metaphors Affect User Experience in Collaborative Virtual Environment[END_REF] or highlight it. In the first case, the designation is dynamic and doesn't modify the 3D scene. In the second case, the designation is more persistent and users may choose between several designation metaphors:

• Bounding volume: the component is surrounded by a transparent sphere. • Spotlight: the component is highlighted by a spotlight and ambient light is lowered in order to obtain a more noticeable contrast. • Flashing Texture: the texture of the component is visible intermittently. • Outline shader: the component is surrounded by a color line.

It remains visible if it is hidden by other objects.

Remote Visualization

In order to better perceive what other users are seeing, multiple views are often displayed in a splitscreen [START_REF] Dodds | Using teleporting, awareness and multiple views to improve teamwork in collaborative virtual environments[END_REF][START_REF] Nguyen | Guiding Techniques for Collaborative Exploration in Multi-Scale Shared Virtual Environments[END_REF]. This kind of display is very useful in classical 3D expriments but in a VR context, it can significantly reduce, or even break, the feelings of immersion and presence. For this reason, we enable the immersed user to define 3D camera(s) in the global scene and to display the view(s) in additional window(s) in the VE. A camera and a window are defined as 3D objects in the scene and can be manipulated by users like other objects. They can modify all the parameters of the camera. These cameras can be placed before the session by a master operator. To understand their specific task, users have to check all these preset flows.

Based on the same principle (window and camera), immersed users can get the viewpoint of remote users to see what they are currently doing or what they are talking about.

Component Trajectory

In order to explain/understand a global process, it may be useful to show the trajectory of some components during the progress of the task. For instance, showing the movement of a tool from the workbench to its area of use, or the movement of a component that must be introduced/removed in/from the launcher. This is the reason why we are developing a tool to define trajectories in the virtual environment. The trajectory can be defined during the working session by manually placing and orienting a sequence of points or by capturing a 3D movement using a tracking sensor. These two modes can be mixed ; for instance, by first quickly defining a rough trajectory using some individual points and second, by refining parts that need more precision. A trajectory can be visualized by dynamically moving the component linked to it following the 3D path. The trajectory is a real 3D line (gathering not only the localization of each point but also its orientation) enabling movements of components « around » the trajectory. For instance, we can precisely describe the path of a tool so that it can pass through a small window which forces it to modify its orientation. As with the cameras, a master user can define, before the session, the trajectories of some components during a task. The final user has to play each trajectory to understand the global task.

CONCLUSION

After having developed basic interaction features, we are now designing a test protocol in order to study these tools both quantitatively and qualitatively. We intend to measure some real industrial task executions (time, precision) and to ask users for feedback (comfort, acceptability…). Thanks to these tests, we want to analyze the real benefits taking into account the additional cognitive load. These tests aim to first measure the gains for the users and second try to analyze the contribution of these interaction features towards improving global awareness. More precisely, we intend to quantitatively and qualitatively refine the following figure showing how the interaction features contribute to the collaboration.

Figure 2: This table describes the awareness categories to which the interaction features contribute

In the future, we plan to design new and original interaction features for adding 3D annotations (text, images, sounds) and for helping complex component manipulation involving more than a single user.

Components

  When remote immersed users have to cooperate to perform a collaborative task, they have to know what the other operators are viewing in the VE. So displaying avatars or providing remote visualization is very useful [1, 7, 4].-Procedure: in order to perform complex tasks, users must check the progress of the overall procedure and mainly what other users are doing. Moreover immersed users have to communicate with others to synchronize their own actions with the actions of remote users[START_REF] Fleury | A Generic Model for Embedding Users'Physical Workspaces into Multi-Scale Collaborative Virtual Environments[END_REF][START_REF] Heldal | Immersiveness and symmetry in copresent scenarios[END_REF][START_REF] Hrimech | How 3D Interaction Metaphors Affect User Experience in Collaborative Virtual Environment[END_REF].

Figure 1 :

 1 Figure 1: On the left, an outside view of the digital mock-up of a launcher and on the right, an inside top-down view.


Quentin Mérigot

Jean-Marie Mirebeau

Minimal geodesics along volume preserving maps, through semi-discrete optimal transport

Minimal geodesics along volume preserving maps, through semi-discrete optimal transport

Quentin Mérigot, Jean-Marie

Introduction

The motion of an inviscid incompressible uid, moving in a compact domain X ⊆ R d , is described by the Euler equations [Eul65]

∂ t v + (v • ∇)v = -∇p div v = 0, (1)
coupled with the impervious boundary condition v • n = 0 on ∂Ω. Here v denotes the uid velocity, and p the pressure acts as a Lagrange multiplier for the incompressibility constraint.

As observed by Arnold [START_REF] Arnold | Sur la géométrie diérentielle des groupes de Lie de dimension innie et ses applications à l'hydrodynamique des uides parfaits[END_REF], Euler equations (1) yield in Lagrangian coordinates the geodesic equations along the group SDiff of dieomorphisms of X with unit jacobian, equipped with the L 2 metric. Indeed, let s : [0, 1] → SDiff be a smooth time dependent family of dieomorphisms of X, describing the evolution over time of the position of the uid particles. The position of the particle emanating from x at time t will be denoted indierently s t (x) or s(t, x), while the dieomorphism at time t is denoted s(t) or s t . The Eulerian velocity v and acceleration a of the uid particles are given by v(t, s(t, x)) = ∂ t s(t, x) and a(t, s(t, x)) = ∂ tt s(t, x). The jacobian constraint det ∇s = 1 yields div v = 0, whereas the equation of geodesics on SDiff merely states that acceleration is a gradient, a(t, x) = -∇p(t, x), which is equivalent to (1, Left). This formalism leads to two natural problems for Euler equations:

• The Cauchy problem, forward in time: given the initial position and velocity of the uid particles, nd their subsequent positions at all positive times. This amounts to computing the exponential map on the lie group SDiff.

• The boundary value problem: given some observed initial and nal positions of the uid particles, nd their intermediate states. This amounts to computing a minimizing geodesic on SDiff.

This manuscript is devoted to the boundary value problem only. More precisely, consider an inviscid incompressible uid owing during the time interval [0, 1], and a map s * : X → X giving the nal position s * (x) of each uid particle initially at position x ∈ X. In this paper, we introduce a new numerical scheme for nding a minimizing geodesic joining the initial conguration (

) 2
We denoted by S ⊆ L 2 (X, R d ) the space of maps preserving the Lebesgue measure on X, which in dimension d ≥ 3 is the closure of SDiff. More formally, we let Leb denote the Lebesgue measure on X, rescaled so as to be a probability measure, and we let f #µ denote the push-forward of a measure µ by a measurable map f , dened by f #µ(A) := µ(f -1 (A)). Then the space S is the collection of measurable maps s : X → X obeying s# Leb = Leb.

The motivation for this rst relaxation replacing SDiff with the closed subset S of L 2 is that the optimized functional in (2) does not penalize the spatial derivatives of s, which are involved in the unit jacobian constraint dening SDiff. Despite this relaxation, the optimization problem (2) needs not have a minimizer in s ∈ H 1 ([0, 1], S) in dimension d ≥ 3 [Shn94], and minimizing sequences (s n ) n∈N may instead display oscillations reminiscent of an homogeneization phenomenon. A second relaxation is therefore necessary, see [START_REF] Figalli | Variational models for the incompressible Euler equations[END_REF] for a review.

Generalized ows Brenier introduced in [START_REF] Brenier | The least action principle and the related concept of generalized ows for incompressible perfect uids[END_REF] a second relaxation of the minimizing geodesic problem (2) based on generalized ows, which allow particles to split and their paths to cross. This surprising behavior is an unavoidable consequence of the lack of viscosity in Euler's equations, which amounts to an innite Reynolds number. Generalized ows are also relevant in dimension d ∈ {1, 2} if the underlying physical model actually involves a three dimensional domain X × [0, ε] 3-d in which one neglects the uid acceleration in the extra dimensions [Bre08].

Consider the space of continuous paths (of uid particles) Ω := C 0 ([0, 1], X).

A generalized ow, in Brenier's sense [START_REF] Brenier | The least action principle and the related concept of generalized ows for incompressible perfect uids[END_REF], is a probability measure µ over the space Ω of paths. We denote by e t (ω) := ω(t) the evaluation map at time t ∈ [0, 1]. Given a generalized ow µ ∈ Prob(Ω), the pushforward measure e t #µ can be understood as the distribution of particles at time t under the ow. The generalized ow is incompressible if for every time t ∈ [0, 1], e t #µ coincides with the uniform probability measure on the domain X. Rescaling the domain if necessary we assume that X has unit area, so that the uniform probability measure Leb is simply the restriction of the Lebesgue measure to X.

The use of generalized ows turns the highly non-linear incompressibility constraint s(t) ∈ S into a linear constraint e t #µ = Leb. Note the similarity with Kantorovich's linearization of the non-linear mass preservation constraint in Monge's optimal transport problem. This idea leads to a convex relaxation of the minimizing geodesic distance problem (2), linearizing both X v(x) s 0 (x)

s t (x) X t = 0 t = 1 X X Figure 1:
The motion of inviscid incompressible uids admits three formulations, either (I)

Eulerian based on the local speed v : [0, 1] × X → R d , (II) Lagrangian based on dieomorphisms s(t, •) which integrate the speed: ∂ t s(t, x) = v(t, s(t, x)), or (III) relaxed as a superposition of individual particles paths ω ∈ Ω, weighted by a measure µ.

the objective and the constraints, and for which the existence of a minimizer is guaranteed: (3)

As explained above, the rst constraint in (3) expresses the incompressibility of the generalized ow. The second constraint (e 0 , e 1 )#µ = (s * , s * )# Leb rephrases the prescription of the endpoints (s(0) = s * and s(1) = s * ) in the minimizing geodesic problem (2) by imposing a coupling between particle positions at initial and nal times. This second constraint can equivalently be phrased in term of test functions: for every f ∈ C 0 (X × X, R) one has ˆΩ f (ω(0), ω(1)) dµ(ω) = ˆX f (s * (x), s * (x)) dx.

In other words, ω(1) = s * (s -1 * (ω(0))) for µ-almost any particle path ω ∈ Ω. Any classical ow s ∈ H 1 ([0, 1], S) can be represented by a generalized ow µ s , supported on the paths ω x : t → s(t, x), weighted by the Lebesgue measure on x ∈ X. In this case, if F : Ω → R is a bounded and continuous functional (or a lower-semi-continuous non-negative functional such as

A) one has ˆΩ F (ω) dµ s (ω) = ˆX F (ω x ) dx. (4)
In particular,

ˆΩ A(ω) dµ s (ω) = ˆX ˆ1 0 | ṡ(t, x)| 2 dt dx = ˆ1 0 ṡ(t) 2 dt. (5)
In other words, the µ-average of the particles paths energy, is the time average of the uid kinetic energy. Note that in this paper, • stands for the L 2 (X, R n ) norm, and Pressure. The incompressibility constraint in (3) gives rise to a Lagrange multiplier, called the pressure and which generalises the eld p in (1). Pressure turns out to be the unique maximizer

m 0 m T -1 s * . . . π S (m T -1 ) . . . m 1 m T s * π S (m 1 ) S M N Figure 2:
The geodesic distance d 2 (s * , s * ) along the manifold S of volume preserving maps, represented as a blue curve, is estimated (9) as the length of a chain (m i ) T i=0 in the linear subspace M N , represented as a black line, plus penalizations for the boundary values and the distance from the chain elements to S.

to a concave optimization problem dual to (3), see [Bre93]. Note that, in contrast, the primal problem (3) may have several minimizers up to the notable exception [START_REF] Bernot | Generalized solutions for the Euler equations in one and two dimensions[END_REF] of smooth ows

in dimension d = 1. The pressure is a classical function p ∈ L 2 loc ( ]0, 1[, BV(X)), at least when the domain X is a d-dimensional torus [AF07]
. This regularity is sucient to show that any solution s to (2) (resp. µ-almost any path ω, for any solution µ to (3)) satises

∂ tt s(t, x) = -∇p(t, s(t, x)), resp. ω(t) = -∇p(t, ω(t)). (6)
This implies that the support of µ is contained in the space of solutions to a second-order Ordinary Dierential Equation (ODE), a fact that we will use later in our convergence estimates. In other words, uid particles move by inertia, only deected by the force of pressure. Assume that the Hessian of the pressure is suciently small, namely that

∀t ∈ [0, 1], ∀x ∈ X, ∇ 2 p(t, x) ≺ π 2 Id (7)
in the sense of symmetric matrices 1 . Then using the path dynamics equation (6) Brenier [START_REF] Brenier | The least action principle and the related concept of generalized ows for incompressible perfect uids[END_REF] proved that the relaxed problem (3) admits a unique minimizer µ ∈ Prob(Ω), which is deterministic. In other word, the generalized ow µ stems as in (4) from a possibly non-smooth but otherwise classical minimizer s ∈ H 1 ([0, 1], S) of (2). Inequality (7) is sharp, and several families of examples are known for which uniqueness and/or determinism are lost precisely when the threshold (7) is passed. We present in 5 the rst numerical illustration of this phenomenon.

Numerical scheme

We introduce a new discretization for the relaxation (3) of the shortest path formulation (2) of Euler equations (1). Our approach is numerically tractable in dimension 2, and is the rst to illustrate the transition between classical and generalized solutions occurring at the threshold (7) on the pressure regularity.

For that purpose we need to introduce some notation. As before, we denote Leb the Lebesgue measure restricted to the domain X, which by assumption is a probability measure on X. Denote M := L 2 (X, R d ), and recall that S ⊆ M is the collection of maps s ∈ M preserving the Lebesgue measure on the domain X, i.e. s# Leb = Leb. For every N ∈ N, we let P N be a partition of X into N regions of equal area 1/N , with respect to the Lebesgue measure, and with diameter ≤ C P N -1 d . Finally, we let M N ⊆ M be the N -dimensional subspace of functions which are piecewise constant over this partition. Given two measure-preserving maps s * , s * ∈ S, 1 We write A ≺ B, where A and B are two symmetric matrices, if B -A is positive denite. discretization parameters T, N ∈ N, and a penalization factor λ 1, we introduce the functional which to m ∈ M T +1

N associates E(T, N, λ; m) := T 0≤i<T m i+1 -m i 2 + λ m 0 -s * 2 + m T -s * 2 + 1≤i<T d 2 S (m i ) . (8)
We denoted by d 2 S the squared-distance to the space S of measure-preserving maps, dened as usual by d 2 S (m) = inf s∈S ms 2 . Comparing this with (2), we recognize the standard discretization of the length of the discrete path (m 0 , • • • , m T ), as well as an implementation by penalization of the boundary value constraints. The last term corresponds to a penalization of the incompressibility constraints, or more precisely a penalization of the squared distance between the maps m i and the set S of measure preserving maps. The discrete optimization problem we consider is then

E(T, N, λ) := min m∈M T +1 N E(T, N, λ; m) (9)
The problem (9), considered as a function of m ∈ M However it is non-convex (in fact semi-concave) due to the non-convexity of the set S, forbidding us to guarantee that the global minimum of (9) is found by our numerical solver. Nonetheless, our numerical experiments show that quasi-Newton methods give convincing results, see 5.

Distance to S and semi-discrete optimal transport. Before entering the analysis of (9), we want to emphasize that the inner-subproblems, namely the computation of the squared distances d 2 S (m i ), are numerically tractable thanks to two main ingredients: Brenier's polar factorization [START_REF] Brenier | Polar factorization and monotone rearrangement of vector-valued functions[END_REF], and semi-discrete optimal transport. Brenier's polar factorization theorem asserts that the distance between any m ∈ M and the set S of measure-preserving maps equals the cost of transporting the image measure m# Leb back to the Lebesgue measure, namely

d 2 S (m) = inf s∈S m -s 2 = W 2 2 (m# Leb, Leb), (10)
where W 2 is the Wasserstein distance for the quadratic transport cost. Given any m ∈ M N , the pushforward measure m# Leb is the sum of N Dirac measures of mass 1/N , located at the N values of the piecewise constant map m on the partition P N . Computing (10) therefore amounts to the computation of the Wasserstein distance between a nitely supported measure and the Lebesgue measure, a special instance of a problem called semi-discrete optimal transport. Ecient computational methods for this problem have been proposed [START_REF] Aurenhammer | Minkowski-Type Theorems and Least-Squares Clustering[END_REF][START_REF] Merigot | A Multiscale Approach to Optimal Transport[END_REF][START_REF] Lévy | A numerical algorithm for L 2 semi-discrete optimal transport in 3D[END_REF],

relying on Kantorovich duality and on tools from computational geometry. Let η be the uniform probability measure over a nite set Y ⊆ R d with cardinal N . Kantorovich duality asserts that

W 2 2 (η, Leb) = sup f :X→R,g:Y →R f (x)+g(y)≤ x-y 2 ˆX f (x) d Leb(x) + ˆX g(y) dη(y).
For any xed g : Y → R the largest f : X → R obeying the constraint is given by f (x) := min y∈Y xy 2g(y).

This function is piecewise quadratic over a partition (Lag y (g)) y∈Y of X into convex polyhedra

Lag y (g) := x ∈ X | ∀z ∈ Y, x -y 2 -g(y) ≤ x -z 2 -g(z) .
Eliminating the optimization variable f , Kantorovitch duality reads

W 2 2 (η, Leb) = sup g:Y →R y∈Y ˆLag y (g) ( x -y 2 -g(y)) dx + 1 N y∈Y g(y). (11)
This partition of X is called the Laguerre diagram of (Y, g) in computational geometry and can be computed in near-linear time in R 2 using existing software 2 [cga]. Specically, we use the 2D regular triangulations package in CGAL [Yvi13]. The maximization problem (11) is an unconstrained, concave and twice continuously dierentiable maximization problem, which is eciently solved via Newton or quasi-Newton methods. Semi-discrete optimal transport has become a reliable and ecient building block for PDE discretizations [START_REF] Benamou | Discretization of functionals involving the Monge-Ampère operator[END_REF].

As discussed above, the chosen discretization of the Euler equations intrinsically leads to an optimal transport problem of semi-discrete type -from a nite collection of Diracs playing the role of uid particles to the standard Lebesgue measure -hence our choice of discretization (11).

A variety of alternative approaches have been developed for other instances of computational optimal transport, based for instance on a pressure-less uid model

Main results

Our results on the proposed discretization (8) are split into two parts. A convergence analysis rst shows how to construct discrete generalized ows out of minimizers of E(T, N, λ), supported on nitely many piecewise linear paths and which converge to a solution µ ∈ Prob(Ω) of (3). This analysis requires upper bounds on the discrete energies E(T, N, λ), which are established in a second part and under adequate a-priori assumptions on the continuous solution. Importantly, two sets of assumptions are considered, encompassing the favorable case where a classical solution to Arnold's formulation (2) exists, and the general case where only Brenier's relaxation (3) is well posed.

Convergence analysis. We present two results establishing the convergence of the discretized problem (9) minimizers towards solutions of (3) as the number of timesteps T grows to innity. The rst proposition shows that one can build short chains of incompressible maps out of minimizers of (9). It involves the quantity

E (T, N, λ; m) := (1 + 4T /λ)E(T, N, λ; m). ( 12
)
2 In R 3 the worst-case complexity of the Laguerre diagram is quadratic in the number of points, but this corresponds to degenerate cases, so that semi-discrete optimal can also be solved eciently on R 3 [START_REF] Lévy | A numerical algorithm for L 2 semi-discrete optimal transport in 3D[END_REF].

. Let m = (m i ) T i=0 ∈ M T +1
N and let (s i ) T i=0 be the chain of incompressible maps dened by: s 0 = s * , s T = s * , and s i is a projection of m i onto S for all 1 ≤ i < T . Then

T 0≤i<T s i+1 -s i 2 ≤ E (T, N, λ; m).
Unfortunately, such discrete chains of incompressible maps need not converge as T → ∞ to a continuous path s : [0, 1] → S, due to the lack of compactness of S. Note also that the interpolates

(1 -δ)s i + δs i+1 , 0 ≤ δ ≤ 1, 0 ≤ i ≤ T , need not belong to the non-convex set S.
The next proposition establishes, under condition, convergence towards a solution of (3) in the sense of generalized ows, a framework which enjoys the desired convexity and compactness properties. In order to state this result, we associate to each m = (m i ) T i=0 in M T +1 N a generalized ow µ m ∈ Prob(Ω) uniformly supported on N paths. For every j ∈ {1, . . . , N }, we denote m j i the (constant) value of m i on the j-th region of the partition P N of X and we dene a path ω j by interpolating linearly between these points, that is

∀0 ≤ i ≤ N, ω j (i/t) = m j i .
The construction of these paths is illustrated in Figure 3. Finally µ m ∈ Prob(Ω) is the uniform probability measure over the nite set of paths {ω 1 , . . . , ω N } ⊆ Ω. We present in 2.2 a reformulation of the discrete energy (8) based on the generalized ow µ m ∈ Prob(Ω), which mimics Brenier's linear relaxation (3) of Euler's equations.

Proposition

1.2. Assume that E (T k , N k , λ k ; m k ) → D * and T k , N k , λ k → +∞ as k → ∞.
Then a subsequence of the generalized ows µ k associated to m k weak-* converges to a ow µ ∈ Prob(Ω) satisfying ´Ω A dµ ≤ D * and obeying the constraints of (3).

Corollary 1.3. Under the assumptions of Proposition 1.2, one has D * ≥ d(s * , s * ). If equality holds, then the limit µ is a solution to (3).

Interestingly, the proof of Proposition 1.2, presented 2.2, relies on a reinterpretation of the discretization (8) of Arnold's geodesic formulation (2) of Euler equations, which turns out to be the natural discretization of Brenier's relaxation (3) as well.

m 0 m 1 • • • m T m 0 m 1 m T Figure 4: Illustration of Remark 1.6.
The set S is the red curve and the set S N is composed of the four large dots. The discretized geodesics is displayed as a dotted line. On the left, the penalization parameter λ is too small and the discretized geodesics takes shortcuts, while on the right λ is too large, leading to overtting.

Upper bounds on the discrete energy.

E(T, N, λ) ≤ d(s * , s * ) + O(T h 2 N λ), • (Classical estimate) with h N = N -1
d assuming that the classical geodesic distance (2) equals the relaxed distance (3), and admits a minimizer with regularity s ∈ L ∞ ([0, 1], H 1 (X)).

• (Relaxed estimate) with h N = N -1 2d (resp. N -1 2 √ ln N if d = 1
), if the pressure eld gradient ∇p is Lipschitz on [0, 1] × X, and the boundary data s -1 * , s * are Lipschitz on X.

Remark 1.5 (Optimal choice of the parameter λ). Since the natural quantity occurring in the

convergence analysis is E (T, N, λ), see (12), it is natural to choose λ = h -1 N = N 1 D so that by Theorem 1.4 E (T, N, λ) ≤ d 2 (s * , s * ) + O(T /λ + T h 2 N λ) ≤ d 2 (s * , s * ) + O(T N -1 D ).
From this point, choosing N, T → ∞ in such way that T N -1 D → 0, will fulll the assumptions required by Proposition 1.2 for convergence.

Remark 1.6 (Finite-dimensional analogy for the choice of λ). Consider a nite d-dimensional manifold S embedded in a Hilbert space M , whose minimising geodesics one wishes to approximate. The minimising geodesics between two points s * , s * ∈ S can be approximated by minimisers of the energy

(m 0 , . . . , m T ) ∈ M T +1 → T 0≤i<T m i+1 -m i 2 +λ   m 0 -s * + m T -s * + 1≤i<T d S N (m i )   ,
where S N is a sampling of the manifold S with N points. The choice of the optimal λ in term of convergence speed depends on the dimension of the manifold S and on the number of points N . Figure 4 illustrates that if λ is too small, the discrete geodesics takes shortcuts, while if λ is too large it forms clusters around points in the discretization S N . Note also that, for a xed N and λ → ∞, the optimization problem becomes combinatorial.

Establishing the upper bounds of Theorem 1.4 requires one to construct, for each choice of N, T , a candidate m ∈ M T +1 N for the discrete problems (9) from a minimizer µ ∈ Prob(Ω) of the continuous problem (3). To do so, we approximate the probability measure µ by a measure µ N uniformly supported over N paths. This amounts to solving a quantization problem [START_REF] Gersho | Vector Quantization and Signal Compression[END_REF] in the innite-dimensional space Ω. Fortunately, the optimal speed of convergence of µ N towards µ can be bounded in terms of the box dimension D of the support of µ (see Denition 4.1). This explains why the decay rate h N = N -1 D in Theorem 1.4 depends, through the constant D, on the structure of the problem solution µ.

When the generalized ow µ is induced by a classical solution, as in (4), the particle trajectories are determined by their initial position x ∈ X ⊆ R d , so that µ is supported on a d-dimensional manifold, implying D = d. In the relaxed estimate however, the trajectories particles obey a second order ordinary dierential equation ( 6) and are thus determined by their initial position and velocity (x, v) ∈ X × R d ⊆ R 2d provided Picard-Lindelof/Cauchy-Lipschitz's theorem applies [START_REF] Earl | Theory of ordinary dierential equations[END_REF]. The quantization dimension of the generalized ow is thus 2d in the worst case, but intermediate dimensions d < D < 2d are also common, see 5.

Outline. Propositions 1.1 and 1.2 are proved in 2. Theorem 1.4 is established in 3 (Classical estimate) and in 4 (Relaxed estimate). Numerical experiments are presented 5.

Remark 1.7 (Monge-Ampere gravitation). Some models of reconstruction of the early universe [Bre11] involve actions of a form closely related to our discrete energy functional (9), for the parameter value λ = 2:

ˆ1 0 1 2 ṁ(t) 2 + inf s∈S m(t) -s 2 dt.

Convergence analysis

We prove Proposition 1.1 in 2.1 and Proposition 1.2 in 2.2, on the length of chains of incompressible maps, and the convergence of generalized ows respectively.

Length of a chain of incompressible maps

The announced Proposition 1.1 immediately follows from Lemma 2.1 below. It relies on a the following identity, valid for any elements a, b of a Hilbert space, and any ε > 0:

(

1 + ε) -1 a + b 2 ≤ a 2 + ε -1 b 2 . ( 13
)
Indeed subtracting the LHS to the RHS of (13) we obtain

(1 + ε) -1 ε 1 2 a -ε -1 2 b 2 ≥ 0.
Lemma 2.1. For any T ∈ N * , any penalization λ > 0, and any (m, s) ∈ (M × S) T +1 one has

T 0≤i<T s i+1 -s i 2 ≤ (1 + 4T /λ)   T 0≤i<T m i+1 -m i 2 + λ 0≤i≤T m i -s i 2   . (14)
Proof. Let 0 ≤ i < T . Choosing a := s i+1s i , and b := m i+1m ia, we obtain

(1 + ε) -1 m i+1 -m i 2 ≤ s i+1 -s i 2 + ε -1 (s i -m i ) -(s i+1 -m i+1 ) 2 ≤ s i+1 -s i 2 + 2ε -1 ( s i -m i 2 + s i+1 -m i+1 2 ).
Summing over 0 ≤ i < T yields

(1 + ε) -1 0≤i<T m i+1 -m i 2 ≤ 0≤i<T s i+1 -s i 2 + 2ε -1 0≤i≤T α i s i -m i 2 , with α 0 = α T = 1, α i = 2 otherwise. Choosing ε = 4T /λ concludes the proof.

Convergence towards generalized solutions

)# Leb on X 2 as required in (3). Lemma 2.2. Let m = (m i ) T i=0 ∈ M T +1 N and µ m ∈ Prob(ω) be the induced generalized ow. Then E(T, N, λ; m) equals ˆΩ A(ω)dµ m (ω)+λ ˆX |m 0 (x)-s * (x)| 2 +|m T (x)-s * (x)| 2 dx+ 1≤i<T t=i/T W 2 2 (e t #µ m , Leb) . ( 15
)
In addition,

W 2 2 ((e 0 , e 1 )#µ m , (s * , s * )# Leb) ≤ ˆX |m 0 (x) -s * (x)| 2 + |m T (x) -s * (x)| 2 dx. (16)
Proof. Mimicking (5) at the discrete level, we obtain

ˆΩ A(ω)dµ(ω) = 1 N 1≤j≤N ˆ1 0 | ωj (t)| 2 dt = T N 0≤i<T 1≤j≤N |m j i+1 -m j i | 2 = T 0≤i<T m i+1 -m i 2 ,
On the other hand, Brenier's polar factorization theorem (10

) yields W 2 2 (e t #µ m , Leb) = d 2 S (m),
which gives (24).

The Wasserstein distance from (e 0 , e 1 )#µ m = (m 0 , m 1 )# Leb to (s * , s * )# Leb appears on the LHS of (16). By denition, it is indeed bounded by the cost, appearing on the RHS, of the explicit transport plan π = ((m 0 , m T ), (s

* , s * ))# Leb ∈ Prob(X 2 × X 2 ), which sends (m 0 (x), m T (x)) → (s * (x), s * (x)) for all x ∈ X.
Our next step is to estimate the average violation of incompressibility by the discrete solution.

Lemma 2.3. Let µ m be the generalized ow associated to some m ∈ M T +1

N . Then ˆ1 0 W 2 2 (e t #µ m , Leb) dt ≤ 1 4T 2 E (T, N, λ; m).
Proof. Let t ∈ [0, 1], which we write t = (i + α)/T with 0 ≤ i < T , 0 ≤ α ≤ 1. Then using (13) we obtain for any ε > 0,

W 2 2 (e t #µ m , Leb) = inf s∈S (1-α)m i +αm i+1 -s 2 ≤ (1+ε) α(m i+1 -m i ) 2 + ε -1 inf s∈S m i -s 2
Integrating over t ∈ [0, 1], and using that either α ≤ 1/2 or 1α ≤ 1/2, we obtain

ˆ1 0 W 2 2 (e t #µ, Leb)dt ≤ 1 + ε T 0≤i<T 1 4 m i+1 -m i 2 + ε -1 inf s∈S m i -s 2
which is precisely the announced estimate when ε = 4T /λ.

F (µ k ) ≤ 1 λ k E(T k , N k , λ k ; m k ) → 0 as k → ∞, hence F (µ * ) = 0 which implies
(e 0 , e 1 )#µ * = (s * , s * )# Leb .

Similarly, the functional G : µ → ´1 0 W 2 2 (e t #µ, Leb) dt is weak-* lower semi-continuous on Prob(Ω), as follows from Fatou's lemma and the continuity of µ → W 2 2 (e t #µ, Leb) for any t ∈ [0, 1]. By Lemma 2.3 one has G(µ k ) → 0 as k → ∞, hence G(µ * ) = 0 and therefore the limit generalizered ow obeys the incompressibility constraints ∀t ∈ [0, 1], e t #µ = Leb .

Classical estimate

We establish in this section the rst part of our main result, Theorem 1.4. By assumption, we consider a minimizer s of the shortest path problem (2), and assume that it has regularity s ∈ L ∞ ([0, 1], H 1 (X)). Let also T, N, λ be the discrete problem parameters.

Dene s i := s(i/T ) for all 0 ≤ i ≤ T , and note that s 0 = s * and s T = s * . Let also m i := P N (s i ), for all 0 ≤ i ≤ T , where P N : M → M N denotes the orthogonal projection. This construction is illustrated on Figure 5. The proof, short and elementary, proceeds by estimating the distance from m i to s i and to m i+1 , and then summing over the index 0 ≤ i ≤ T . Distance from m i to the set of incompressible maps. Denote by h N := N -1 d the discretization scale, and recall that each region of the partition P N of X has area 1/N and diameter ≤ C P h N . Let s P i denote the mean of s i on the region P of the partition P N , for all 0 ≤ i ≤ T .

s 0 = s * s i = s(i/T ) m i := π M N (s i ) S M N
Then

s i -m i 2 = P ∈P N ˆP |s i (x) -s P i | 2 dx ≤ C sb (C P h N ) 2 P ∈P N ˆP |∇s i (x)| 2 dx = Ch 2 N ∇s i 2 , (17)
where the Sobolev inequality constant C sb only depends on the dimension, and

C := C sb C 2 P .
Distance from m i to m i+1 . The map P N is 1-Lipschitz, as it is the orthogonal projection onto the linear subspace M N . Hence for any 0

≤ i < T m i -m i+1 2 ≤ s i -s i+1 2 ≤ 1 T ˆi+1 T i T ṡ(t) 2 dt. (18)
Summation and conclusion. Summing (17) and (18) over 0 ≤ i ≤ T we obtain

E(T, N, λ; m) ≤ T 0≤i<T m i+1 -m i 2 + λ 0≤i≤T m i -s i 2 ≤ 0≤i<T ˆi+1 T i T ṡ(t) 2 dt + λ 0≤i≤T Ch 2 N ∇s i 2 ≤ d 2 (s * , s * ) + C T h 2 N λ, where C = C s L ∞ ([0,1],H 1 (X))
, which concludes the proof.

Relaxed estimate

We prove the relaxed estimate in Theorem 1.4 using a quantization of the generalized ow minimizing the relaxed geodesic distance (3), in other words we approximate the solution µ ∈ Prob(Ω) with probability measures (µ N ) N ≥1 equidistributed on a family of N paths. This quantization is a counterpart of the partition P N of the domain (X, Leb) used for the classical estimate 3, which amounts to quantize the initial positions of the uid particles. The proof is split into four parts. In 4.1, we introduce the framework of measure quantization. In 4.2, we specialize it to the case of Brenier's generalized ows. In 4.3 we perform a permutation of the obtained quantization, so as to match the boundary conditions. The announced energy estimate is nally proved in 4.4.

Box and quantization dimensions

In sampling theory, a natural objective is to approximate a given measure with a nite sum of Dirac masses in the Wasserstein distance, or to approximate a given set with a nite point set in the Haussdor distance, see Figure 6. The quantization dimension of a measure, and the box dimension of a set, are dened as the optimal convergence rates of these approximation procedures.

We denote by δ x the Dirac probability measure concentrated at a point x. The Haussdor distance between two closed non-empty subsets A, B of a metric space is dened by Denition 4.1. Let H be a metric space, let µ be a probability measure on H, and let Γ ⊆ H be closed and non-empty. For all N ≥ 1 denote

h N (µ) := inf ω∈H N W 2 µ, 1 N 1≤i≤N δ ω i , r N (Γ) := inf ω∈H N Hauss Γ, {ω i } N i=1 .
The quantization dimension of µ, and the (upper) box dimension of Γ, are dened by

D quant (µ) := lim sup N →∞ ln N -ln h N (µ) , D box (Γ) := lim sup N →∞ ln N -ln r N (Γ)
.

The box dimension D box of a set also called the Minkowski dimension is a variant of the Haussdor dimension in which the set of interest if covered with balls of identical radius.

Indeed, one easily shows that r N (Γ) is the minimal radius allowing to cover Γ with N balls r N (Γ) := inf

ω∈H N min r ≥ 0; Γ ⊆ 1≤i≤N B(ω i , r) . (19)
Box and Haussdor dimension coincide for compact manifolds, but dier in general. For instance, all countable sets have Haussdor dimension zero, whereas one can check that

D box ([0, 1] ∩ Q) d = d, D box 1 n ; n ∈ N * = 1 2 .
The quantization dimension of a measure is bounded above by the box dimension of its support, as shows the following elementary result of quantization theory. We refer the reader to [START_REF] Gersho | Vector Quantization and Signal Compression[END_REF] for more details on this rich subject.

Proposition 4.2. Let H be a metric space, and let µ ∈ Prob(H) be supported on a set Γ. Then D quant (µ) ≤ max{2, D box (Γ)}. More precisely for any D > 0, one has as N → ∞

r N (Γ) = O(N -1 D ) ⇒ h N (µ) = O      N -1 D if D > 2, N -1 2 √ ln N if D = 2, N -1 2 if D < 2. (20)
Proof. Let N ∈ N be xed. For each 1 ≤ i ≤ N let M i ⊆ H be a set of i points such that Γ ⊆ ∪ ω∈M i B(ω, 2r i ), with r i := r i (Γ). We construct a sequence of points ω i ∈ H, and an increasing sequence of measures ρ i supported on Γ and of mass i/N , inductively starting with i = N and nishing with i = 1. Initialization: ρ N := µ.

Induction: for each 1 ≤ i ≤ N , we construct ω i and ρ i-

1 in terms of ρ i . Let indeed ω i ∈ M i be such that B i := B(ω i , 2r i ) satises ρ i (B i ) ≥ 1/N . Such a point exists since |ρ i | = i/N , #(M i ) = i, and supp(ρ i ) ⊆ Γ. Then let ρ i-1 := ρ i -1 N ρ i (B i ) ρ i , so that ρ i -ρ i-1 is a non-negative measure of mass 1 N supported on B i . One has h N (µ) 2 ≤ W 2 2 µ, 1 N 1≤i≤N δ ω i ≤ 1≤i≤N W 2 2 ρ i -ρ i-1 , 1 N δ ω i ≤ 1 N 1≤i≤N (2r i ) 2 .
The comparison (20) of the decay rates of h N (µ) and r N (Γ) immediately follows. Finally the comparison of the dimensions follows from (20).

Quantization dimension of a generalized ow

We specialize in this subsection the measure µ, support Γ and embedding metric space H to which we apply the results of 4.1. Let µ ∈ Prob(Ω) be a generalized ow minimizing the relaxed geodesic distance (3). This measure is concentrated on the set Γ of paths obeying Newton's second law of motion

Γ := {ω ∈ C 2 ([0, 1], X); ∀t ∈ [0, 1], ω(t) = -∇p(t, ω(t))}, (21)
where the pressure gradient ∇p : [0, 1] × X → R d is assumed to have Lipschitz regularity, following the assumptions of Theorem 1.4. We regard Γ as embedded in the Hilbert space H := H 1 ([0, 1], R d ), which plays a natural role in the problem of interest (3) and is equipped with the norm

ω 2 H := ˆ1 0 |ω| 2 + | ω| 2 .
Note that H continuously embeds in C 0 ([0, 1], R d ), hence the evaluation maps e t : H → R d : ω → ω(t) are continuous with a common Lipschitz constant denoted C e .

Notation 4.3 (Problem constants). The following quantities are referred to as the problem constants: the dimension d, the constant C P related to the partitions, the Lipschitz regularity constants C * and C * of the boundary conditions s -1 * and s * , the maximum norm and the Lipschitz regularity regularity constant of the pressure gradient ∇p. Given two expressions A, B, we write A B i A ≤ CB for some constant C that only depends on the problem constants.

Using a general existence and stability result for ODEs, Lemma 4.4, we show that Γ is in bijection with a compact subset of R 2d , in Lemma 4.5. Lemma 4.4. Let F ∈ C 0 ([0, 1] × R n , R n ) be Lipscshitz in the second variable, uniformly with respect to the rst variable: ∃L, ∀t ∈

[0, 1], ∀y 1 , y 2 ∈ R n , |F (t, y 1 )-F (t, y 2 )| ≤ L|y 1 -y 2 |. Then the map R n → L 2 ([0, 1], R n ), associating to each y 0 ∈ R n the solution to the ODE y (t) = F (t, y(t))
with initial condition y(0) = y 0 , is bi-Lipschitz onto its image with regularity constant e L .

Proof. The Picard-Lindelof/Cauchy-Lipschitz Theorem [START_REF] Earl | Theory of ordinary dierential equations[END_REF] guarantees the existence of a unique solution to the considered ODE. If y 1 , y 2 are two such solutions, then for any t ∈ [0, 1]

| ẏ1 (t) -ẏ2 (t)| = |F (t, y 1 (t)) -F (t, y 2 (t))| ≤ L|y 1 (t) -y 2 (t)| Thus for any t ∈ [0, 1] e -Lt |y 1 (0) -y 2 (0)| ≤ |y 1 (t) -y 2 (t)| ≤ e Lt |y 1 (0) -y 2 (0)|, hence by integration over t ∈ [0, 1], as announced e -L |y 1 (0) -y 2 (0)| ≤ y 1 -y 2 L 2 ([0,1]) ≤ e L |y 1 (0) -y 2 (0)|.
The image of the generalized ow µ by the map of initial position and momentum, used in the next lemma, is often called a minimal measure [START_REF] Bernot | Generalized solutions for the Euler equations in one and two dimensions[END_REF].

Lemma 4.5. The map E 0 : Γ → X × R d : ω → (ω(0), ω(0)) is bijective and bi-Lipschitz onto its image, which is a compact set. The Lipschitz régularity constants, and the image diameter, can be bounded in terms of the problem constants, see Notation 4.3.

Proof. The function F (t, x, v) := (v, -∇p(t, x)) is Lipschitz on [0, 1] × X × R d , with a regularity constant determined by that of ∇p. It can be extended to [0, 1] × R 2d . Applying Lemma 4.4 to F in dimension n = 2d, we obtain as announced that E 0 is bi-Lipschitz onto its image. Note that the norm (ω, ω) L 2 ([0,1]) appearing in this estimate is precisely ω H .

By construction, ω(t) and ω(t) are uniformly bounded for all ω ∈ Γ, and all t ∈ [0], respectively by the diameter of X and the maximum norm of ∇p. This implies an upper bound on ω(0), observing for instance that

|ω(0) -ω(1)| = ˆ1 0 ω(t) dt ≥ | ω(0)| - ˆ1 0 | ω(t) -ω(0)|dt ≥ |ω(0)| - ˆ1 0 |ω(t)|dt.
As announced, the image of Γ by E 0 is bounded, hence compact since it is clearly closed.

Since there is no ambiguity, we denote in the rest of this section h N := h N (µ). The information gathered on the set Γ supporting µ yields upper and lower bounds on the decay rate of (h N ), obtained respectively in Corollary 4.6 and Proposition 4.7.

Corollary 4.6. One has h

N N -1 2d (resp. h N N -1 2 √ ln N if d = 1) for all N ≥ 2.
Proof. By Lemma 4.5, the set Γ is in bi-Lipschitz bijection with a compact set K ⊆ R 2d . Hence r N (Γ) r N (K) N -1 2d , and the upper estimate follows from (20).

The quantization scale h N is also bounded below, and is minimal for classical solutions.

Proposition 4.7. One has h N N -1 d for all N ≥ 1. If the generalized ow µ in fact represents a classical solution s to Euler's equations, and ∇ ṡ is bounded on [0, 1] × X, then this lower estimate is sharp:

h N N -1 d .
Proof. Since X is a domain of unit area, there exists c > 0 depending only on the dimension d and such that W 2 (Leb, ν N ) ≥ cN -1 d for any measure ν N supported at N points of R d . The rst point follows: for any measure µ N supported at N points of H

cN -1 d ≤ W 2 (Leb, e 0 #µ N ) = W 2 (e 0 #µ, e 0 #µ N ) ≤ C e W 2 (µ, µ N ) = C e h N .
Second point: for each x ∈ X, introduce the path ω x : [0, 1] → R d : t → s(t, x). Then Φ : (X, Leb) → (Γ, µ) : x → ω x is measure preserving and Lipschitz, with regularity constant denoted C Φ (which we regard as a problem constant). Let ν N be a discrete probability measure, with one Dirac mass of weight 1/N in each region of the partition P N . Since these regions have diameter

≤ C P N -1 d , we conclude that h N ≤ W 2 (µ, Φ#ν N ) = W 2 (Φ# Leb, Φ#ν N ) ≤ C Φ W 2 (Leb, ν N ) ≤ C Φ C P N -1 d .

Permutation of the quantization of µ and boundary conditions

We have shown 4.2 that the generalized ow, solution to Brenier's relaxation of Euler's equations

(3), could be approximated (quantized) by a nite family of particle paths, with a controlled convergence rate. This is not yet sucient however to construct a good candidate for the discrete energy (8), despite its relaxed reformulation (24). For that purpose one needs to permute the paths, so as to match the boundary conditions, as required by the energy term (16).

We show in Lemma 4.8 that an optimal quantization exists, we suitably permute it in Lemma 4.9 and estimate a rst boundary term, and we estimate the second boundary term in Proposition 4.10. In this subsection and the next, we x the integer N and allow ourselves a slight abuse of notation: elements ω j , P j , ρ j , . . . indexed by 1 ≤ j ≤ N do implicitly depend on N , although that second index ω N j , P N j , ρ N j , . . . is omitted for readability.

Lemma 4.8. The inmum dening h N is attained, see Denition 4.1. As a result there exists

(ω j ) N j=1 ∈ H N and probability measures (ρ j ) N j=1 on Γ such that µ = 1 N 1≤j≤N ρ j h 2 N = 1 N 1≤j≤N ˆΓ ω -ω j 2 H dρ j (ω) (22)
Furthermore, ω j is the barycenter of ρ j for each 1 ≤ j ≤ N .

Proof. Let (ω j ) N j=1 be a candidate quantization, and let π be the transport plan associated to W 2 2 ( 1 N N j=1 δ ω j , µ). Then the measures ρ j : A → N π({ω j } × A), 1 ≤ j ≤ N , are probabilities which average to µ, and the transport cost is the RHS of ( 22). The quantization energy, i.e. the squared Wasserstein distance, is decreased when replacing ω j with the barycenter b j of ρ j , 1 ≤ j ≤ N , by the amount 1 N N j=1 |ω j -b i | 2 . Hence ω j = b j for all 1 ≤ j ≤ N if the quantization is optimal. Note also that the barycenter of ρ j belongs to G := Hull(Γ) by construction.

Since Γ is a compact subset of a Hilbert space, the convex hull closure G is also compact, for the strong topology induced by

• H . The quantization energy (ω j ) N j=1 → W 2 2 ( 1 N N j=1
δ ω j , µ) attains its minimum on G N by compactness, and by the previous argument it is the global minimum on H N . Let µ N denote the equidistributed probability on the set {ω j ; 1 ≤ j ≤ N } of Lemma 4.8. Lemma 4.9. The paths (ω j ) N j=1 of the quantization µ N , and the regions (P j ) N j=1 of the partition P N , can be indexed in such way that

1≤j≤N ˆPj |ω j (0) -s * (x)| 2 dx h 2 N . (23)
Proof. Let (P j ) N j=1 be an arbitrary indexation of P N . Let a j be the barycenter of region P j , and b j the barycenter of the image region s * (b j ). For any 1 ≤ j ≤ N one has by assumption ˆPj

|b j -s * (x)| 2 dx ≤ ˆPj |s * (a j ) -s * (x)| 2 dx ≤ C 2 * ˆPj |a j -x| 2 dx ≤ 1 N (C * C P N -1 d ) 2 ,
where C * denotes the Lipschitz regularity constant of s -1 * . Denoting by ν N the discrete measure equidistributed at the N points (b j ) N j=1 , we obtain

W 2 (ν N , Leb) ≤ 1≤j≤N ˆPj |b j -s * (x)| 2 dx N -1 d . On the other hand W 2 (Leb, e 0 #µ N ) ≤ C e W 2 (µ, µ N ) = C e h N . Thus W 2 (ν N , e 0 #µ N ) h N ,
since W 2 denes a distance, and since h N N -1 d by Proposition 4.7. This optimal transport problem between the discrete measures ν N and e 0 #µ N determines an optimal assignment Γ N → B N . Up to a permutation of the paths (ω j ) N j=1 , we assume that this assignment is represented by the common indexation (ω j ) N j=1 of Γ N and (b j ) N j=1 of B N . Finally

1≤j≤N ˆPj |ω j (0) -s * (x)| 2 dx = 1≤j≤N ˆPj |b j -s * (x)| 2 dx + W 2 2 (ν N , e 0 #µ N ) h 2 N .
Proposition 4.10. Using the indexation of Lemma 4.9, one has

1≤j≤N ˆPj |ω j (1) -s * (x)| 2 dx ≤ Ch 2 N .
Proof. The generalized boundary condition of (3) states that ω(1) = s * * (ω(0)) for µ-almost every ω ∈ Γ, where s * * := s * • s -1 * . Denoting by C * * the Lipschitz regularity constant of s * * , we obtain for any 1 ≤ j ≤ N and any x ∈ X

|ω j (1) -s * (x)| 2 = ˆΓ s * * (ω(0)) -s * (x) dρ j (ω) 2 ≤ ˆΓ |s * * (ω(0)) -s * (x)| 2 dρ j (ω) ≤ (C * * ) 2 ˆΓ |ω(0) -s * (x)| 2 dρ j (ω) ≤ 2(C * * ) 2 ˆΓ |ω(0) -ω j (0)| 2 dρ j (ω) + |ω j (0) -s * (x)| 2 .
By summation and integration we conclude, using Lemma 4.9 in the second line,

1≤j≤N ˆPj |ω j (1) -s * (x)| 2 dx 1≤j≤N 1 N ˆΓ |ω j (0) -ω(0)| 2 dρ j (ω) + 1≤j≤N ˆPj |ω j (0) -s * (x)| 2 dx ≤ C 2 e W 2 2 (µ N , µ) + Ch 2 N h 2 N .

Final estimate

We conclude in this subsection the proof of Theorem 1.4 (Relaxed estimate). Let µ ∈ Prob(Ω) be a generalized ow minimizing (3), with an associated Lipschitz pressure gradient. We x the parameters (T, N ), and consider an optimal quantization µ N supported on N paths (ω j ) N j=1 as in Lemma 4.8. For each 0 ≤ i ≤ T , let m i ∈ N be the piecewise constant map on the partition P N dened by ∀1 ≤ j ≤ N, ∀x ∈ P j , m i (x) := ω j (i/T ).

In the following, we bound each of the terms of the energy E(T, N, λ; m) ≥ E(T, N, λ). By Lemma 4.9 and Proposition 1.2, one has respectively

m 0 -s * 2 h 2 N , m T -s * 2 h 2 N .
Regarding the distance to incompressible maps, one has for any 1 ≤ i < T , with t

:= i/T inf s∈S m i -s = W 2 (Leb, e t #µ N ) = W 2 (e t #µ, e t #µ N ) ≤ C e W 2 (µ, µ N ) = C e h N .
We use the Cauchy-Schwartz inequality to bound the action:

T 0≤i<T m i+1 -m i 2 = 1 N 1≤j≤N T 0≤i<T ω j i + 1 T -ω j i T 2 ≤ 1 N 1≤j≤N ˆ1 0 | ωj (t)| 2 dt = 1 N 1≤j≤N ˆ1 0 ˆΓ ω dρ j (ω) 2 ≤ 1 N 1≤j≤N ˆ1 0 ˆΓ | ω| 2 dρ j (ω) dt = d 2 (s * , s * ).
Concluding, the value E(T, N, λ) is bounded by

T 0≤i<T m i+1 -m i 2 +λ   m 0 -s * 2 + m T -s * 2 + 1≤i<T inf s∈S m i -s 2   ≤ d 2 (s * , s * )+O(T h 2 N λ).
5 Numerical experiments

Minimization algorithm and choice of penalization

We rely on a quasi-Newton method to compute a (local) minimum of the discretized problem (9), as summarized in Algorithm 1. This means that we need to compute the value of the functional

m ∈ M T +1 N → T 0≤i<T m i+1 -m i 2 + λ m 0 -s * 2 + m T -s * 2 + 1≤i<T d 2 S (m i ) . ( 24
)
and its gradient, where d 2 S (m) = inf s∈S ms 2 . The only diculty is to evaluate the squared distance d 2 S to the set of measure-preserving vector elds and its gradient. As explained in the introduction, Brenier's Polar Factorization Theorem implies that for any vector-valued function

m ∈ M, d 2 S (m) = W 2 2 (m# Leb, Leb).
When m belongs to M N , the measure m# Leb is nitely supported, and the computation of the Wasserstein distance can be performed using a solver for semi-discrete optimal transport, such as [START_REF] Aurenhammer | Minkowski-Type Theorems and Least-Squares Clustering[END_REF][START_REF] Merigot | A Multiscale Approach to Optimal Transport[END_REF][START_REF] Fernando De Goes | Blue noise through optimal transport[END_REF][START_REF] Lévy | A numerical algorithm for L 2 semi-discrete optimal transport in 3D[END_REF]. The next proposition gives an explicit formulation for the gradient in term of the optimal transport plan. Note that the expression of the gradient of the semi-discrete Kantorovitch functional with respect to the position of the points also appears in the appendix of [START_REF] Fernando De Goes | Blue noise through optimal transport[END_REF].

Recall that M N is the set of piecewise constant functions on the tessellation P

N := (P j ) 1≤j≤N of X. The diagonal D N in M N is the set of functions m in M N such that m(P j ) = m(P k ) for some j = k. The set M N \ D N is a dense open set in M N .
Proposition 5.1. The functional d 2 S is dierentiable almost everywhere on M N and continuously dierentiable on M N \ D N . The gradient of d 2 S at m ∈ M N \ D N is explicit: with x j = m(P j ),

∇d 2 S (m) P j = 2(x j -bary(T -1 (x j ))) (25)
where T : X → m(X) is the piecewise constant optimal transport map between Leb and the nitely supported measure m# Leb, and bary(S) = ´S xdx/ Leb(S) is the isobarycenter of S

Proof. The functional F := d 2 S -• 2 is concave as an inmum of linear functions:

F(m) = d 2 S (m) -m 2 = inf s∈S m -s 2 -m 2 = inf s∈S -2 m|s + s 2 ,
where •, • denotes the L 2 (X) scalar product. This implies in particular that F and therefore d 2 S is dierentiable almost everywhere on M N . Given m in M N \ D N , dene x j = m(P j ) and let T : X → R d be the optimal transport plan from Leb to m# Leb = 1 N N j=1 δ x j . The transport plan is indeed always representable by a function when the source measure is absolutely continuous with respect to the Lebesgue measure. Let V j = T -1 (x j ) be the partition of X induced by this transport plan. Then

F(m) = W 2 2 (m# Leb, Leb) -m 2 = N j=1 ˆVj x j -x 2 -x j 2 dx = m|G(m) + N j=1 ˆVj x 2 dx
where G(m) ∈ M N is the piecewise constant function on X given by G(m)| P j = -2 bary(V j ).

For any m in M N and x j = m (V j ), one has

F(m ) = W 2 2 (m # Leb, Leb) -m 2 ≤ N j=1 ˆVj x j -x 2 -x j 2 dx = F(m) + m -m|G(m)
This shows that G(m) belongs to the superdierential to F at m. In addition, by the continuity of optimal transport plans, the map m ∈ M N \ D N → G(m) is continuous. To summarize, on the open domain M N \ D N the concave function F possesses a continuous selection of supergradient. This implies that F is of class C 1 on this domain, with ∇F (m) = G, and the result follows.

Remark 5.2 (Computation of d 2 S ). The computation of the squared distance to measure-preserving maps d 2 S and its gradient rely on the variational approach used in the proof of the above propo- sition. We use the CGAL library [cga] to evaluate F and its rst and second derivatives, and a simple damped Newton's algorithm to compute its maximum. The implementation is available at https://github.com/mrgt/PyMongeAmpere.

Construction of the initial solution Since the discrete energy (8) is non-convex, the construction of the inital guess is important. We follow a time-renement strategy already used by Brenier [Bre08] to construct a good initial guess. Assuming that we have already a local minimizer for T k = 2 k +1, we use linear interpolation to construct an initial guess for T k+1 = 2 k+1 +1.

The optimization is then performed from this inital guess, using a quasi-Newton algorithm for the energy (8).

Choice of the penalization parameter The optimal choice of λ in (8) depends on the quantization dimension D = D quant (µ) of the generalized solution µ ∈ Prob(Ω) that one expects to recover: namely λ N = N -1 D , see Remark 1.5. We call D the ow dimension, and regard it as as the intrinsic dimensionality of the problem which determines its computational diculty. For a classical solution, this dimension agrees with the ambient dimension i.e. D = d, while for a non-deterministic solution the quantization dimension can be up to 2d. Intermediate dimensions d < D < 2d are also common [START_REF] Brenier | The least action principle and the related concept of generalized ows for incompressible perfect uids[END_REF]. In our numerical experiments we set λ N = N update (m k i ) 0≤t≤T using the L-BFGS scheme see [START_REF] Richard H Byrd | A limited memory algorithm for bound constrained optimization[END_REF]

Visualization of generalized solution

The main interest of numerical experimentation is to visualize generalized solutions to Euler's equation, or equivalently generalized geodesics between two measure-preserving dieomorphisms s * , s * in S.

Gradient of the pressure

Consider a minimizer of the discretized energy (8). Given i ∈ {1, . . . , T -1}, m i minimizes over

M N the functional m → T ( m -m i-1 2 + m i+1 -m 2 ) + λd 2 S (m). This gives T 2 (m i-1 -2m i + m i+1 ) = T λ∇d 2 S (m i ). ( 26
)
This equation is a discretized counterpart of the rule that the acceleration of a geodesic on an embedded manifold, is normal to that manifold (here S plays the role of the manifold, embedded in M, which is internally approximated by the linear space M N ). The second order dierence T 2 (m i-1 -2m i + m i+1 ) approximates a second derivative in time. Comparing (26) to (6), we see that the right hand-side of (26) can be used as an estimation of (minus) the gradient of the pressure.

Geometric data analysis

As in the proof of Theorem 1.4, the discrete minimizer of (8) can converted to a collection of N piecewise-linear curves {ω 1 , . . . , ω N } = Γ N . We recall that the domain X is partitioned into N subdomains (P j ) 1≤j≤N with equal area and we let ω j (i/T ) ∈ R d be the point corresponding to the restriction of m i to the subdomain P j , for each 0 ≤ i ≤ T . Figure 3

C i = ω ∈ Γ N ; ω -i H = arg min 1≤j≤k ω -i H .
Note that l 1 , •

As in Denition 4.1, denote by r i = r i (Γ N ) is the smallest r ≥ 0 such that Γ N can be covered by i balls of radius r. For 1 i N , the ratio log(i)/ log(1/r i (Γ N )) is expected to approximate log(i)/ log(1/r i (supp(µ))) and thus the desired D box (supp µ).

Lemma 5.3. Let ε i := max γ∈Γ d(γ, {γ 1 , . . . , γ i }), where γ i is dened as in (28). Then,

1 - log(2) log(1/ε i ) log(i) log(1/ε i ) ≤ log(i) log(1/r i ) ≤ log(i) log(1/ε i )
Proof. By construction, r i ≤ ε i . Moreover, the balls centered at the points γ 1 , . . . , γ i and with radius ε i /2 are disjoint, so that r i ≥ ε i /2.

Test cases and numerical results

Our two testcases are constructed from two stationary solutions to Euler's equation in 2D. Let s : R + → S be a classical solution to Euler equation in Lagrangian coordinates (6), starting from the identity map. We solve the discretized version (9) of the minimization problem (2)-(3), with s * = s(0) = Id and s * := s(t max ), where t max > 0. For small values of t max the solution to this boundary problem is simply the original classical ow s, but for larger values a completely dierent generalized ow is obtained. In this case the geodesic s in the space of the measure preserving dieomorphisms is no longer the unique shortest path between its boundary values s * and s * . The rst classical behavior is guaranteed if the pressure hessian satises

∇ 2 p ≺ (π/t max ) 2 Id (29)
uniformly on [0, t max ] × X, see (7) and [START_REF] Brenier | The least action principle and the related concept of generalized ows for incompressible perfect uids[END_REF]. In all the numerical experiments, the number of points is set to N = 10 000 and the number of timesteps is T = 2 4 + 1 = 17.

Rotation of the disk

On the unit disk D = {(x 1 , x 2 ) ∈ R 2 ; x 2 1 + x 2 2 ≤ 1}, the simplest stationary solution to Euler's equation (1) is given by a time-independent pressure eld and speed:

p(x 1 , x 2 ) = 1 2 (x 2 1 + x 2 2 ), v(x 1 , x 2 ) = (-x 2 , x 1 ).
The corresponding Lagrangian ow s(t) is simply the rotation of angle t. The largest eigenvalue of ∇ 2 p is 1 at every point in D. Hence by (29) the ow of rotations is the unique minimizer to both the variational formulation (2) and its relaxation (3) with boundary values s * = s(0) = Id and s * = s(t max ), when t max < π. Uniqueness is lost at the critical time t max = π which corresponds to a rotation of angle π, so that the nal dieomorphism becomes s * = s(π) = -Id.

In this situation, the minimization problem (2) has two classical solutions, namely the clockwise and counterclockwise rotations. The relaxation (3) has uncountably many generalized solutions such as, by linearity, superpositions of these two rotations.

Another explicit example of generalized solution was discovered by Brenier [START_REF] Brenier | The least action principle and the related concept of generalized ows for incompressible perfect uids[END_REF]: given a point x ∈ D and a speed v, denote by ω x,v the curve ω x,v (t) = x cos(t) + v sin(t), t ∈ [0, 1]. Then, Brenier's solution is obtained as the pushforward by the map (

x, v) → ω x,v ∈ Ω of the measure on D × R 2 dened by µ(dx, dv) = 1 π H 2 (dx) ⊗ 1 2π 1 -|x| 2 H 1 {|v|= √ 1-|x| 2 } (dv),
where H k denotes the k-dimensional Hausdor measure. In particular, the quantization dimension of the solution is 3 = 2 + 1. We refer to [START_REF] Bernot | Generalized solutions for the Euler equations in one and two dimensions[END_REF] for more examples of optimal ows, and construct four dimensional one. Let µ r be dened by combining (i) a classical rotation on the annulus D \ D(r), with D(r) = {x ∈ R 2 ; |x| ≤ r} and (ii) Brenier's solution rescaled by a factor r on the disc D(r). Then µ r is an optimal generalized ow of quantization dimension 3, whereas the averaged ow ´1 0 µ r dr is also optimal by linearity, and has quantization dimension 4.

Numerical results

The numerical solutions computed by our algorithm for the critical time t max = π are highly non-deterministic. To see this, we select a small neighborhood around several points in the unit disk D and look at the trajectories emanating from this small neighborhood.

As shown in Figure 7, the trajectories emanating from each neighborhood at initial time tend to ll up the disk at intermediate times before gathering again on a small neighborhood at nal time. In addition, each indivual trajectory looks like an ellipse, as in Brenier's explicit solution.

Second, we estimate the box dimension of the support of the numerical solution (as explained in 5.2). The estimated dimension is slightly above 3.

Beltrami ow on the square

On the unit square S = [-1/2, 1/2] 2 , we consider the Beltrami ow constructed from the timeindependent pressure and speed:

p(x 1 , x 2 ) = 1 2 (sin(πx 1 ) 2 + sin(πx 2 ) 2 ) v(x 1 , x 2 ) = (-cos(πx 1 ) sin(πx 2 ), sin(πx 1 ) cos(πx 2 ))

The maximum eigenvalue of the pressure Hessian ∇ 2 p is π 2 , and [START_REF] Brenier | The least action principle and the related concept of generalized ows for incompressible perfect uids[END_REF] implies that the associated ow is minimizing between s * = s(0) = Id and s * = s(t max ) for t max ≤ 1. Because the square has much less symmetries than the disk, generalized solutions constructed from this ow are not as well understood as those occuring from rotations of the disk.

Numerical results Our numerical results suggest the following observations. First, as shown in Figure 9, the computed solutions with boundary values s * = Id and s * = s(t max ) approximate the classical ow if t max < 1, and are non-deterministic generalized ows if t max > 1. This suggests the sharpness of the bound given by [START_REF] Brenier | The least action principle and the related concept of generalized ows for incompressible perfect uids[END_REF]. Interestingly, even for t max > 1, the numerical solutions seem to remain deterministic in a neighborhood of the boundary of the cube.

This can be seen more clearly in Figure 10, where the particles have been divided into clusters using the k-means algorithm (the clustering algorithm is explained in 5.2.2).

The pressure gradient is estimated as in 5.2.1 and is displayed in Figure 11. These pictures seem to indicate a loss of regularity of the pressure near the initial and nal times. This corroborates the result of [START_REF] Ambrosio | On the regularity of the pressure eld of Brenier's weak solutions to incompressible Euler equations[END_REF] according to which the pressure belongs to L 2 loc ( ]0, T [, BV(X)). Figure 7 suggests that the even for t max = 1.5, the reconstructed solution for the Beltrami ow are more deterministic than the solution to the disk inversion. We estimate the box dimension of the support of the solution using the method explained in 5.2.2. The results are displayed in Figure 8. The estimated dimension is D = 2 for the deterministic solution (t max = 0.9) but it increases as the maximum time (and therefore the amount of non-determinism) increases.

Finally, we note that the estimated dimensions for t max ∈ {1.1, 1.3, 1.5} seem to be strictly between 2 and 3, suggesting a fractal structure for the support of the solution. This would need to be conrmed by a mathematical study. [Sch15] Bernhard Schmitzer. A sparse algorithm for dense optimal transport. In Jean-François Aujol, Mila Nikolova, and Nicolas Papadakis, editors, Scale Space and . . . , pages 629641. Springer International Publishing, 2015.

[Shn94] A I Shnirelman. Generalized uid ows, their approximation and applications.

Geometric and Functional Analysis, 1994.

[Yvi13] Mariette Yvinec. 2D triangulations. In CGAL User and Reference Manual. CGAL Editorial Board, 4.2 edition, 2013.

s

  * = Id to the nal one s * minimize ˆ1 0 ṡ(t) 2 dt, subject to ∀t ∈ [0, 1], s(t) ∈ S, s(0) = s * , s(1) = s * .

  d 2 (s * , s * ) := min µ∈Prob(Ω) ˆΩ A(ω)dµ(ω), where A(ω) = ˆ1 0 | ω(t)| 2 dtand subject to ∀t ∈ [0, 1], e t #µ = Leb, (e 0 , e 1 )#µ = (s * , s * )# Leb .

  [BBG02, PPO13], monotone discretizations of Monge-Ampere equation [BFO14, BCM14], regularization via entropy penalization [BCC + 15, CPR15], or adaptive strategies for the linear programming discretization [OR15, Sch15].

Figure 3 :

 3 Figure 3: (Left) A partition P N cuts the domain X into N region of equal area and roughly isotropic shape. (Right) To a sequence (m i ) T i=0 ∈ M T +1 N one can associate N piecewise linear paths (ω j ) N j=1 , by interpolating the map values at the times {0, 1/T, • • • , T /T } for each region of the partition P N . Proposition 1.1. Let m = (m i ) T

Figure 5 :

 5 Figure 5: The classical estimate in Theorem 1.4 is based on projecting the measure preserving maps (s i ) T i=0 ∈ S T +1 onto the nite dimensional space M N . This procedure is symmetrical to the projection of (m i ) T -1 i=1 ∈ M T -1 N onto S involved in the discrete energy optimization (9), see

Figure 6 :

 6 Figure 6: Left: Illustration of measure quantization. Two absolutely continuous probability measures over the disk (with constant density and a truncated Gaussian density) are approximated by nite collections of Dirac masses (shown as black dots) using the algorithm of De Goes et al [dGBOD12]. Right: Covering of the coast of England with N balls of minimal radius, as in (19). Image by A. Monnerot-Dumaine, source Wikipedia.

  Hauss(A, B) := max sup a∈A d B (a), sup b∈B d A (b) , where d Y (x) := inf y∈Y d(x, y).

  justied a-posteriori by the numerical estimation of the quantization dimension of the computed solution, see Figure 8. Note that the numerical error in (1.5) is governed (for a xed number T of time steps) by the quantity λ -1 + h 2 N λ, and that h N = O(N -1 2d ) under the assumptions of Theorem 1.4. The choice λ N = N 1 α thus yields a convergent scheme whenever α > d, although convergence rates are improved if α is close to the ow dimension D, so that λ N ≈ N 1 D ≈ h -1 N . Algorithm 1 Computation of a (local) minimizer for (24). For k = 1 . . . k max T := 2 k for all even t ∈ {0, . . . , T }, let m k i := m k-1 i/2 initialize discrete geodesic for all odd t ∈ {1, . . . , T }, let m k i := 1 2 (m k i-1 + m k i+1 ) Until convergence do optimize for all i, compute d S (m i ) and ∇d S (m i ) see Rem. 5.2 compute the discrete path energy and its gradient see Eq. (24)

  Software. The software developed for generating the results presented in this article is publicly available at https://github.com/mrgt/EulerSemidiscrete Acknowledgement The authors thank Y. Brenier for constructive discussions and introducing them to the topic of Euler equations of inviscid incompressible uids.

  (a) (x, y) = (-0.7, 0) (b) (x, y) = (-0.35, 0) (c) (x, y) = (0, 0) (d) (x, y) = (0.2, 0) (e) (x, y) = (0.35, 0) (f) (x, y) = (0.5, 0)

Figure 7 :

 7 Figure 7: We select particles whose initial position lie in a small disk, and display their trajectories according to the computed solution to (8). (Top) For the inversion of the unit disk (Bottom)For the Beltrami ow on the square, with t max = 1.5.

Figure 8 :

 8 Figure 8: Estimation of the box counting dimension of the support of the computed see 5.2.2. (Left) for the inversion of the unit disk (Right) Comparison between the estimated box counting dimensions of the solutions to the Beltrami ow on the square, depending on the maximum time.

Figure 9 :

 9 Figure 9: (First row) Beltrami ow in the unit square at various timesteps, a classical solution to Euler's equation. The color of the particles depend on their initial position. (Second to fthrow) Generalized uid ows that are reconstructed by our algorithm, using boundary conditions displayed in the rst and last column. When t max < 1 we recover the classical ow, while for t max ≥ 1 the solution is not classical any more and includes some mixing.

Figure 11 :

 11 Figure 11: Estimated pressure gradient for the Beltrami ow on the square with t max = 1.5.

  , s * of a subset S of a Hilbert space H, internally approximated by subspaces (M N ) N ≥0 . Our rst result Proposition 1.1 uses no additional assumptions. In contrast,

	T +1 N , is an N (T +1)d-dimensional optimization
	problem. It makes sense to consider (8) and (9) whenever one attempts to approximate a shortest
	path joining two points s the numerical implementation strategy, and the other results 1.2, crucially rely on the specic
	properties of the set S of measure preserving maps. The key properties of the functional (8) follow from those of the squared distance function d 2 S ,
	which are established in Proposition 5.1. It is continuously dierentiable on an open dense set.

*

  Our main theorem establishes upper bounds on the discrete energies E(N, T, λ), for adequate values of the parameters, as required by the convergence result Proposition 1.2. It relies on a-priori assumptions on the structure of solutions to the continuous problem. The relaxed part of this result addresses the general case, based on Brenier's linear relaxation (3) of the Euler equations, and only requires on an assumption on the

pressure regularity. The classical part provides a faster convergence rate under the assumption that Arnold's geodesic formulation (2) admits a minimizer s : [0, 1] × X → X, of suitable regularity. At a given time t, the uid velocity ṡ(t, x) is uniquely determined by position x ∈ X in the classical case, whereas in the relaxed case it is a probabilisitic superposition ω(t) of velocities, weighted by the generalized ow solution µ ∈ Prob(Ω) conditional to ω(t) = x; this is in a sense reminiscent of non-classical, quantum mechanics, hence the choice of adjective. Recall that the classical (2) and relaxed (3) distances are automatically equal in dimension d ≥ 3, and that the pressure eld gradient ∇p is uniquely determined by the boundary values s * , s * . Theorem 1.4. Let s * , s * ∈ S, let T, N ∈ N, and λ ≥ 0. The relaxed geodesic distance (3) and the discretized minimum (9) satisfy

  This subsection is devoted to the proof of Proposition 1.2. The discrete optimization problem (9) enforces the proximity between the two couplings (e 0 , e 1 )#µ m and (s * , s *

is inspired by Arnold's formulation (2) of Euler equations as a minimizing geodesic problem over S. Our rst lemma shows that, in the equivalent form (24), it could as well be regarded as a discretization of Brenier's relaxation (3). The terms W 2 2 (e t #µ m , Leb), where m ∈ M T +1 N and µ m is the induced generalized ow, can indeed be regarded as penalizations of the relaxed incompressibility constraint e t #µ = Leb in (3). The other term W 2 2 ((e 0 , e 1 )#µ m , (s * , s * )# Leb)

  The convergence announced in Proposition 1.2, nally results from general compactness and continuity arguments.Proof of Proposition 1.2. Let T k , N k , λ k , m k , µ k and D * be as in the statement of Proposition 1.2. The sublevel sets of the action, such as for any ε > 0

{µ ∈ Prob(Ω); ˆΩ A(ω)dµ(ω) ≤ D * + ε}. are weak-* sequentially compact by [Bre93]. Taking a subsequence if necessary, we assume that (µ k ) converges towards a µ * ∈ Prob(Ω), which thus satises ´Ω A dµ * ≤ D * + ε for any ε > 0, hence ´Ω A dµ * ≤ D * as announced.

The functional F : µ → W 2 2 ((e 0 , e 1 )#µ, (s * , s * )# Leb) is weak-* continuous on Prob(Ω). By (16) one has

  illustrates this construction. We regard Γ N as embedded in the Hilbert space H := H 1 ([0, 1], R 2 ) which plays a natural role in the problem of interest, as in 4, and apply techniques from the eld of geometric data analysis.Clustering In order to better visualize the solution, we use the k-means algorithm to divide the set Γ N into k. A distinct particle color is attached to each cluster, see for instance Figure10. The k-means algorithm consists in nding a local minimizer of the optimal quantization problem xed point algorithm, and to divide Γ N into clusters (C i ) 1≤i≤k with

	min 1 ,... k ∈H	1 N	ω∈Γ N	min 1≤i≤k	ω -i	2 H	(27)
	using a simple

  • • , l k automatically belong to Span(Γ N ), hence to the d(T + 1)-dimensional linear subspace of H consisting of piecewise linear paths with nodes ω(t) ∈ R d at times t = i/T , 0 ≤ i ≤ T . This makes (27) tractable.Let µ be the probability measure supported on the set Γ N ⊆ H, with weight 1/N for each element. The k-means problem amounts to nding the probability measure µ k , supported on at most k points of H, which is the closest to µ w.r.t. the Wasserstein distance: (27) equals This quantization problem is closely related to the one dening h N (µ) in Denition 4.1, up to the dierence that the candidate measure µ k need not be equidistributed on the k points of its support.Box dimension A natural objective is to estimate the quantization dimension D quant (µ) of the generalized ow µ ∈ Prob(Ω) minimizing the relaxed problem (3). The probability measure µ N equidistributed on the set Γ N approximates µ, hence we can expect the set Γ N to also approximate supp(µ). The quantization dimension D quant (µ) is dicult to estimate, but by Proposition 4.2 it admits the simpler upper bound D box (supp(µ)). We estimate the latter by applying the furthest point sampling algorithm to the nite metric space Γ N , which denes an ordering on the elements of Γ N as follows: let γ 1 be an arbitrary point of Γ N and dene by

	min #(supp µ k )≤k	W 2 2 (µ, µ k ).

induction γ i+1 := arg max γ∈Γ N d(γ, {γ 1 , . . . , γ i })

‡ ANR grant NS-LBR, ANR-13-JS01-0003-01