Laurent Boudin
email: laurent.boudin@upmc.fr

Laurent Desvillettes
email: laurent.desvillettes@cmla.ens-cachan.fr

Céline Grandmont
email: celine.grandmont@inria.fr

Ayman Moussa
email: ayman.moussa@cmla.ens-cachan.fr

AND Éline Grandmont

Cmla

A M : Cmla

Global existence of solutions for the coupled Vlasov and Navier-Stokes equations

come

Introduction

When we consider a spray, that is a dispersed phase of particles (or droplets) moving inside a continuous fluid, we can choose between several models. One of them, first introduced by Williams in the framework of the combustion theory [START_REF] Williams | Combustion theory[END_REF] (see also [START_REF] Caflisch | Dynamic theory of suspensions with Brownian effects[END_REF]), consists in coupling a kinetic equation with fluid mechanics equations. In such models, the fluid is described through macroscopic quantities (density, velocity, for instance), whereas the particles are represented by a probability density function. We assume that the spray is thin [START_REF] O'rourke | Collective drop effects on vaporizing liquid sprays[END_REF], i.e. the kinetic and fluid equations are coupled through a drag term (which depends on both the fluid macroscopic quantities, and the probability density function and its variables), but not through the volume fraction occupied by the fluid, as in the case of the so-called thick sprays.

In this work, we consider the coupling of the three-dimensional incompressible Navier-Stokes equations with a Vlasov-type equation. This system appears in various situations, for instance, when one wants to model the transport, with no diffusive effects, of an aerosol in a Newtonian, viscous, incompressible airflow, inside the human upper airways [START_REF] Gemci | A numerical and experimental study of spray dynamics in a simple throat model[END_REF][START_REF] Baranger | A modeling of biospray for the upper airways[END_REF]. We here aim to prove global existence of weak solutions for the three-dimensional Vlasov-Navier-Stokes equations, assuming periodicity with respect to the space variable.

The same kind of model has been studied in [START_REF] Hamdache | Global existence and large time behaviour of solutions for the Vlasov-Stokes equations[END_REF] where the author considered, in all dimensions, the unsteady Stokes equations coupled with the Vlasov equation in a bounded domain with reflexion boundary conditions. The main differences between the present work and [START_REF] Hamdache | Global existence and large time behaviour of solutions for the Vlasov-Stokes equations[END_REF] are the following ones.

• We do not neglect the convection term and consider the Navier-Stokes equations.

• Our proof only works in dimension smaller or equal than 3 (because of the use of Sobolev inequalities).

• Our proof is constructive and does not rely on a fixed point argument of Schauder type: we build a sequence of approximated solutions that converges towards a solution of the coupled problem.

Note that Goudon, Jabin and Vasseur [START_REF] Goudon | Hydrodynamic limit for the Vlasov-Navier-Stokes equations. I. Light particles regime[END_REF][START_REF] Goudon | Hydrodynamic limit for the Vlasov-Navier-Stokes equations[END_REF] studied the hydrodynamic limit of the Vlasov-Navier-Stokes system, for some particular regimes of the dispersed phase. In [START_REF] Mellet | Global weak solutions for a Vlasov-Fokker-Planck/Navier-Stokes system of equations[END_REF], Mellet and Vasseur proved the existence of global weak solutions to the compressible Vlasov-Fokker-Planck-Navier-Stokes equation (hence including diffusive effects for the particles and compressibility effects for the fluid), for absorption or reflection boundary conditions, and in [START_REF] Mellet | On the barotropic compressible Navier-Stokes equations[END_REF], they investigate some asymptotic regimes of this system.

Local existence of solutions in the case of the compressible Vlasov-Euler equation was proved by Baranger and Desvillettes [START_REF] Baranger | Coupling Euler and Vlasov equations in the context of sprays: the local-in-time, classical solutions[END_REF]. One can find, in that same article, references on previous works about fluid-kinetic coupling, in particular, when the fluid is modelled by one single equation [START_REF] Domelevo | Existence and stability of travelling wave solutions in a kinetic model of two-phase flows[END_REF][START_REF] Goudon | Asymptotic problems for a kinetic model of two-phase flow[END_REF]. The extension of this result to moderately thick sprays (that is, sprays in which the collisions between droplets has been taken into account) has been done by Mathiaud [START_REF] Mathiaud | Étude de systèmes de type gaz-particules[END_REF].

Date: September 30, 2008.

We present our model in the next section, along with our main results. Those are proven in Sec. 3.

Presentation of the model

Let us consider a constant T > 0 and set T 3 = R 3 /Z 3 . We investigate the interaction of a (monodispersed) spray of particles with a Newtonian, viscous, and incompressible fluid. The fluid is commonly described by macroscopic quantities, such as its density, here assumed to be constant (equal to 1), its velocity u(t, x) and pressure p(t, x) at time t ∈ [0, T ] and position x ∈ T 3 . The evolution of (u, p) is here governed by the incompressible Navier-Stokes equations. The particles are described by a density function f in the phase space which solves a Vlasov-like equation. The quantity f (t, x, v) is the density of particles located in x ∈ T 3 at time t which have the velocity v ∈ R 3 . Equations for the fluid and the spray are coupled through a drag force, which depends on the relative velocity of the fluid and the particles, namely uv and on the density function f . More precisely, f , u and p satisfy the following system

∂ t f + v • ∇ x f + ∇ v • [(u -v)f ] = 0 in (0, T ) × T 3 × R 3 , (1)
∂ t u + (u • ∇ x )u + ∇ x p -∆ x u = - R 3 f (t, x, v)(u(t, x) -v) dv + F ext in (0, T ) × T 3 , (2)
∇ x • u = 0 in (0, T ) × T 3 , (3)
where F ext is an external force applied to the fluid that is periodic in space. The system has to be completed by initial conditions for u and f , which are

f (0, x, v) = f in (x, v), x ∈ T 3 , v ∈ R 3 , (4)
u(0, x) = u in (x), x ∈ T 3 , (5)
where f in ≥ 0 and u in are given.

Remark 1. Note that we here make the thin spray assumption. First, there is no interaction (collisions or coalescences) between particles, which implies that (1) is linear with respect to f . Secondly, the drag force is the only coupling phenomenon between the spray and the fluid and appears in both (1) and (2): the volume fraction of the spray remains negligible with respect to the fluid volume fraction, and is not taken into account explicitly.

If we assume that all quantities are smooth enough, system (1)-( 5) satisfies, at least formally, an energy equality. It is obtained in the following way. We first multiply the Navier-Stokes equations (2) by u and integrate over T 3 . Then we multiply the Vlasov equation (1) by |v| 2 /2 and integrate over T 3 × R 3 . Finally, we add both previous equalities and obtain, after integrations by parts,

(6) 1 2 d dt u 2 L 2 (T 3 ) + 1 2 d dt T 3 ×R 3 f (t, x, v)|v| 2 + ∇ x u 2 L 2 (T 3 ) + T 3 ×R 3 f (u -v) 2 = T 3 F ext • u.
In the previous equation, the first two terms represent the total kinetic energy of the coupled system, whereas the next two ones represent the energy dissipation coming from, respectively, the fluid itself and the drag force. Note that, since f in is nonnegative, the density function f and consequently

T 3 ×R 3 f (u -v) 2 also remain positive.
Remark 2. When no external force is applied, the global energy

1 2 u(t) 2 L 2 (T 3 ) + 1 2 T 3 ×R 3 f (t, x, v
)|v| 2 of the system is decreasing. We can also point out that, in the same case, the global momentum of the system is conserved.

In the sequel, we shall use the Assumption 1.

(7) f in ∈ L ∞ (T 3 × R 3 ), R 3 |v| 2 f in dv ∈ L ∞ ((0, T ) × T 3 ), u in ∈ L 2 (T 3 ), F ext ∈ L 2 (0, T ; L 2 (T 3 )).
Under this assumption, the energy equality gives our functional framework. More precisely, we set

H = w ∈ L 2 (T 3 ) | ∇ x • w = 0 , V = w ∈ H 1 (T 3 ) | ∇ x • w = 0 .
Definition 1. We say that (u, f ) is a weak solution of (1)-( 5) if the following conditions are satisfied:

• u ∈ L ∞ (0, T ; H) ∩ L 2 (0, T ; V) ∩ C 0 ([0, T ]; V ′ ), • f (t, x, v) ≥ 0, for any (t, x, v) ∈ (0, T ) × T 3 × R 3 , • f ∈ L ∞ (0, T ; L ∞ (T 3 × R 3 ) ∩ L 1 (T 3 × R 3 )), • f |v| 2 ∈ L ∞ (0, T ; L 1 (T 3 × R 3 )), • for all φ ∈ C 1 ([0, T ] × T 3 × R 3 ) with compact support in v, such that φ(T, •, •) = 0, (8)
- T 0 T 3 ×R 3 f [∂ t φ + v • ∇ x φ + (u -v) • ∇ v φ] dx dv ds = T 3 ×R 3 f in φ(0, •, •) dx dv,
• for all ψ ∈ C 1 ([0, T ] × T 3 ) such that ∇ x • ψ = 0, for a. e. t,

T 3 u(t) • ψ(t) dx - t 0 T 3 u • ∂ t ψ dx ds + t 0 T 3 (u • ∇ x ) u • ψ dx ds + t 0 T 3 ∇ x u : ∇ x ψ dx ds = - t 0 T 3 ×R 3 f (u -v) • ψ dx dv ds + t 0 T 3 F ext • ψ dx ds + T 3 u in • ψ(0, •) dx. ( 9
)
In the following, the moments of f are denoted

m α f (t, x) = R 3 f (t, x, v)|v| α dv, M α f (t) = T 3 ×R 3 f (t, x, v)|v| α dx dv,
for any t ∈ [0, T ], x ∈ T 3 and α ≥ 0. Note that

M α f (t) = T 3 m α f (t, x) dx.
The main result of this work is hereby stated. Theorem 1. Let T > 0. Under Assumption 1, there exists at least one weak solution (f, u) to (1)-( 5) on (0, T ) in the sense of Definition 1. Moreover, this solution satisfies the energy equality [START_REF] Goudon | Asymptotic problems for a kinetic model of two-phase flow[END_REF], and the L ∞ bound

(10) f L ∞ ((0,T )×T 3 ×R 3 ) ≤ e 3T f in L ∞ (T 3 ×R 3 ) .
Proof of Theorem 1. The proof of Thm. 1 is split into three main steps.

In the first one, we build a sequence of approximated solutions where the convection velocity of the Navier-Stokes equations and the fluid velocity appearing in the Vlasov equation are regularized. We also cut off the high particle velocities in the drag force applied by the spray on the flow. In order to prove the existence of a solution to this regularized but still coupled system, we build a sequence of uncoupled problems where we first solve the kinetic part, then the fluid part. We prove that the sequence of solutions of these regularized uncoupled systems converges towards the solution of the regularized system. This is the object of Prop. 2 below.

In a second step, we prove that the approximated solutions satisfy uniform bounds (with respect to the regularization) on a certain time interval allowing to pass to the limit when the parameter of regularization tends to 0. This is the object of Prop. 1 below.

We eventually prove that the local (in time) solution obtained in the previous step can in fact be extended to the whole interval [0, T ], thanks to the energy estimate.

2.1. The regularized system. First we define the regularized system. Let ε > 0 and define θ ε a mollifier such that

θ ε = ε 3 θ(x/ε) with θ ∈ C ∞ (R 3 ) , θ ≥ 0, and R 3 θ = 1. Moreover, we introduce γ ε ∈ C ∞ (R 3 ) whose support is included in the ball B(0, 1/ε), 0 ≤ γ ε ≤ 1, γ ε = 1 on B(0, 1/2ε
) and such that γ ε (v) -→ 1 as ε goes to zero. We then consider the following regularization of the initial problem:

∂ t f ε + v • ∇ x f ε + ∇ v • [(θ ε ⋆ u ε -v)f ε ] = 0 in (0, T ) × T 3 × R 3 , (11)
∂ t u ε + ((θ ε ⋆ u ε ) • ∇ x )u ε + ∇ x p ε -∆ x u ε (12) = - R 3 f ε (t, x, v)(u ε (t, x) -v)γ ε (v) dv + F ε ext in (0, T ) × T 3 , ∇ x • u ε = 0 in (0, T ) × T 3 , (13)
where F ε ext is a space regularization of the force F ext . We complete the system by initial conditions for u ε and f ε :

f ε (0, x, v) = f ε in (x, v), x ∈ T 3 , v ∈ R 3 , (14)
u ε (0, x) = u ε in (x), x ∈ T 3 , (15) where f ε in and u ε in are C ∞ approximations of f in and u in such that (f ε in ) converges to f in strongly in L p (T 3 ×R 3 ), for all p < ∞, and weakly in w * -L ∞ (T 3 × R 3 ), (M 2 f ε in )
is uniformly bounded with respect to ε, and (u ε in ) strongly converges to u in in L 2 (T 3 ). We can also assume, without any loss of generality, that f ε in has a compact support in v.

We prove that there exists a solution (u ε , f ε ) of ( 11

)-(15) in L 2 (0, T ; H 1 (T 3 )) ∩ L ∞ (0, T ; L 2 (T 3 )) × L ∞ (0, T ; L ∞ (T 3 × R 3 )) such that ∂ t u ε ∈ L 2 (0, T ; L 2 (T 3
)), and satisfying the following variational formulation:

• for any φ ∈ C 1 ([0, T ] × T 3 × R 3
), for a. e. t,

T 3 ×R 3 f ε (t) φ(t, •, •) dx dv - t 0 T 3 ×R 3 f ε [∂ t φ + v • ∇ x φ + (θ ε ⋆ u ε -v) • ∇ v φ] dx dv ds (16) = T 3 ×R 3 f ε in φ(0, •, •) dx dv, • for all ψ ∈ C 1 ([0, T ] × T 3 ) such that ∇ x • ψ = 0, for a. e. t, (17)
t 0 T 3 ∂ t u ε • ψ dx ds + t 0 T 3 ((θ ε ⋆ u ε ) • ∇ x ) u ε ψ dx ds + t 0 T 3 ∇ x u ε : ∇ x ψ dx ds = - t 0 T 3 ×R 3 f ε (u ε -v) γ ε (v) • ψ dx dv ds + t 0 T 3 F ε ext • ψ dx ds.
We also obtain the following proposition, which is crucial to get the asymptotics ε → 0 and obtain the existence of at least one weak solution to the original problem ( 1

)-(3).

Proposition 1. There exists a time T * ∈ (0, T ], only depending on T and both quantities u in 2

L 2 (T 3 ) + |M 2 f in | + F ext 2 L 2 (0,T ;L 2 (T 3 )) and f in L ∞ (T 3 ×R 3 ) , such that u ε 2 L ∞ (0,T * ;L 2 (T 3 )) + ∇ x u ε 2 L 2 (0,T * ;L 2 (T 3 )) + M 2 f ε L ∞ (0,T * ) (18) ≤ C(T, T * , u in 2 L 2 (T 3 ) + |M 2 f in | + F ext 2 L 2 (0,T ;L 2 (T 3 )) , f in L ∞ (T 3 ×R 3 ) ).
Here C(A, B, ...) denotes a constant depending on A, B...

Moreover, we have

(19) f ε L ∞ ((0,T )×T 3 ×R 3 ) ≤ e 3T f in L ∞ (T 3 ×R 3 ) .
Finally, there exists

(f, u) ∈ L ∞ ((0, T * ) × T 3 × R 3 ) × (L ∞ (0, T * ; L 2 (T 3 )) ∩ L 2 (0, T * ; H 1 (T 3 )) × C 0 ([0, T * ]; V ′ ) such that, up to a subsequence, (f ε ) weakly converges to f in w * -L ∞ ((0, T * ) × T 3 × R 3
) and (u ε ) weakly converges to u in w-L 2 (0, T * ; w-H 1 (T 3 )) and in w * -L ∞ (0, T * ; w-L 2 (T 3 )). Moreover, (u ε ) strongly converges to u in L 2 (0, T * ; L 2 (T 3 )).

Finally, in order to prove the existence of a solution to the regularized system ( 11)-( 15), we have to use another system, which is defined in the next subsection.

2.2. The regularized and uncoupled system. For a fixed ε > 0, let us consider the following problem, for any n ∈ N, defined by u 0 ε = u ε in and

∂ t f n ε + v • ∇ x f n ε + ∇ v • [(θ ε ⋆ u n ε -v)f n ε ] = 0 in (0, T ) × T 3 × R 3 , ( 20
)
∂ t u n+1 ε + ((θ ε ⋆ u n+1 ε ) • ∇ x )u n+1 ε + ∇ x p n+1 ε -∆ x u n+1 ε (21) = - R 3 f n ε (t, x, v)(u n ε (t, x) -v)γ ε (v) dv + F ε ext in (0, T ) × T 3 , ∇ x • u n+1 ε = 0 in (0, T ) × T 3 , (22)
supplemented with initial conditions ( 14)-(15).

For the sake of simplicity, we shall drop the subscript ε, and denote, for instance, f n ε or u n ε by f n or u n .

Proposition 2. For any n ∈ N, there exists a unique solution (f n , u n ) of ( 20)-( 22), supplemented with initial conditions ( 14)

-(15), in C 1 ([0, T ]; C ∞ (T 3 × R 3 )) × H 1 (0, T ; L 2 (T 3 )) ∩ L 2 (0, T ; H 1 (T 3 )) .
Moreover it satisfies the following estimates:

f n L ∞ ((0,T )×T 3 ×R 3 ) ≤ e 3T f in L ∞ (T 3 ×R 3 ) , (23) u n L ∞ (0,T ;L 2 (T 3 ))∩L 2 (0,T ;H 1 (T 3 )) ≤ C(ε), (24)
∂ t u n L 2 (0,T ;L 2 (T 3 )) ≤ C(ε), (25)
where C(ε) does not depend on n. Furthermore, the sequence

(f n , u n ) strongly converges in L ∞ ((0, T ) × T 3 × R 3 ) × L ∞ (0, T ; L 2 (T 3
)) towards a variational solution of ( 11) -(15) [that is, satisfying (16), (17)].

Proof of the main results

We first study the properties of the sequence (f n , u n ) n∈N defined by (20)-( 22), and derive bounds needed for the asymptotics n → +∞. We can then obtain the existence of (f ε , u ε ) ε>0 solution of the regularized problems ( 11)-(15), that is Proposition 2.

Proof of Proposition 2.

3.1.1. Existence and uniqueness of (f n , u n ). First we focus on the kinetic equation (20). For any n, let us assume that u n is given and belongs to H 1 (0, T ; L 2 (T 3 )) ∩ L 2 (0, T ; H 1 (T 3 )). Since we have regularized the fluid velocity appearing in the kinetic equation, the unique solution f n of (20) can be defined thanks to the method of characteristics. Indeed, g n = e -3t f n satisfies the transport equation:

∂ t g n + v • ∇ x g n + (θ ε ⋆ u n -v) • ∇ v g n = 0.
The solution of the previous equation can be written in terms of characteristics. As a matter of fact, let us consider the solutions of the following system

dx n dt (t) = v n (t), (26) dv n dt (t) = (θ ε ⋆ u n )(t, x n (t)) -v n (t), (27)
with initial conditions x n (0) = x and v n (0) = v, and set χ n (t, x, v) = (x n (t), v n (t)) for any (t, x, v). Then we have (28)

g n (t, x, v) = f ε in (χ n (t, x, v)), ∀ t, x, v. Since u n ∈ H 1 (0, T ; L 2 (T 3 )), θ ⋆ u n ∈ C 0 ([0, T ]; C ∞ (T 3 × R 3
)), and consequently

f n lies in C 1 ([0, T ]; C ∞ (T 3 × R 3 )).
Next we focus on the modified Navier-Stokes equation ( 21).

Assume that (f n ,

u n ) is given in L ∞ ((0, T ) × T 3 × R 3 ) × L ∞ (0, T ; L 2 (T 3
)). Since we cut off the high particle velocities, the drag force

R 3 f n (v -u n ) γ(v) dv at least belongs to L 2 (0, T ; L 2 (T 3 )).
Applying standard results about the Navier-Stokes system (with a mollified convection term), one can prove the existence and uniqueness of

u n+1 ∈ L ∞ (0, T ; L 2 (T 3 )) ∩ L 2 (0, T ; H 1 (T 3 )) solving (21).
Note that the uniqueness relies on the fact that the convection velocity has been regularized. Moreover, since the convection velocity has been regularized and the initial velocity is smooth, one can check that

∂ t u n ∈ L 2 (0, T ; L 2 (T 3 ))
. By a simple induction argument, these properties hold for all n ∈ N.

In the next subsection, we derive uniform estimates in these spaces with respect to n satisfied by the sequence (f n , u n ).

3.1.2. Uniform estimates for the sequence (f n , u n ). Using (28), we immediately obtain (23). Note that this estimate does not depend on the fluid velocity u n .

To obtain the uniform bounds on (u n ) n∈N , let us multiply (21) by u n+1 and integrate over T 3 . Since the fluid is incompressible and the problem periodic in the space variable, we have

T 3 ((θ ε ⋆ u n+1 ) • ∇ x )u n+1 • u n+1 dx = 0.
Moreover, thanks to (23) and the use the cutoff function γ ε ,

T 3 ×R 3 f n (v -u n ) • u n+1 γ ε (v) dx dv ≤ C(ε) 1 + u n 2 L 2 (T 3 ) + u n+1 2 L 2 (T 3 ) .
Consequently, using the periodicity, the incompressibility and the fact that

T 3 F ε ext • u n+1 dx ≤ 1 2 u n+1 2 L 2 (T 3 ) + 1 2 F ε ext 2 L 2 (T 3 ) , we get 1 2 d dt u n+1 2 L 2 (T 3 ) + ∇ x u n+1 2 L 2 (T 3 ) ≤ C(ε) 1 + u n 2 L 2 (T 3 ) + u n+1 2 L 2 (T 3 ) . This implies that (u n ) n∈N is uniformly bounded in L ∞ (0, T ; L 2 (T 3 )) ∩ L 2 (0, T ; H 1 (T 3 )) (see Lemma 3 in Appendix A).
Estimate (25) is obtained by multiplying (21) by ∂ t u n+1 . Indeed, integrating over T 3 , using estimate (24) on (u n ), the fact that R 3 f n (vu n ) γ(v) dv is uniformly bounded in L 2 (0, T ; L 2 (T 3 )) with respect to n and the fact that the convection velocity has been regularized, we deduce the uniform bound of (∂ t u n ) n∈N in L 2 (0, T ; L 2 (T 3 )).

Remark 3. One could also obtain a uniform (in n) bound of (f n ), for instance, in W 1,∞ ((0, T ) × T 3 × R 3 ).

3.1.3. Convergence. We now prove that the sequence (f n ,

u n ) is convergent in L ∞ (((0, T ) × T 3 × R 3 ) × L ∞ (0, T ; L 2 (T 3 )).
To that purpose, we study the differences w n+1 = u n+1u n and f n+1 -f n , for any n ∈ N.

The function w n+1 satisfies the following weak formulation:

T 3 ∂ t w n+1 • ψ + T 3 ((θ ε ⋆ w n+1 ) • ∇ x ) u n+1 • ψ + T 3 ((θ ε ⋆ u n ) • ∇ x ) w n+1 • ψ + T 3 ∇ x w n+1 : ∇ x ψ = T 3 ×R 3 f n (u n -v) γ(v) • ψ - T 3 ×R 3 f n-1 (u n-1 -v) γ(v) • ψ, for any ψ ∈ L 2 (0, T ; H 1 (T 3 )) such that ∇ x • ψ = 0.
The function ψ = w n+1 is then an admissible test function. The only terms which we must take care of are the nonlinear ones: the convection one

T 1 = T 3 ((θ ε ⋆ w n+1 ) • ∇ x ) u n+1 • w n+1 + T 3 ((θ ε ⋆ u n ) • ∇ x ) w n+1 • w n+1 ,
and the drag force one

T 2 = T 3 ×R 3 f n (u n -v) γ(v) • w n+1 - T 3 ×R 3 f n-1 (u n-1 -v) γ(v) • w n+1 .
Thanks to the space periodicity and the fluid incompressibility, we have

T 3 ((θ ⋆ ε w n+1 ) • ∇ x ) u n+1 • w n+1 = - T 3 ((θ ⋆ w n+1 ) • ∇ x ) w n+1 • u n+1 , and T 3 ((θ ⋆ ε u n ) • ∇ x ) w n+1 • w n+1 = 0.
Consequently, using (24), we successively get

|T 1 | ≤ θ ε ⋆ w n+1 L ∞ (T 3 ) u n+1 L 2 (T 3 ) ∇ x w n+1 L 2 (T 3 ) , ≤ C(ε) w n+1 L 2 (T 3 ) ∇ x w n+1 L 2 (T 3 )
. Finally, we have

|T 1 | ≤ C(ε) w n+1 2 L 2 (T 3 ) + 1 2 ∇ x w n+1 2 L 2 (T 3 ) .
We now take care of T 2 :

T 2 = T 3 ×R 3 (f n -f n-1 )(u n -v) γ ε (v) • w n+1 - T 3 ×R 3 f n-1 (u n -u n-1 ) γ ε (v) • w n+1 .
Taking into account the cutoff function and ( 23)-( 24), we get

|T 2 | ≤ C(ε) f n -f n-1 L ∞ ((0,T )×T 3 ×R 3 ) w n+1 L 2 (T 3 ) + w n+1 L 2 (T 3 ) w n L 2 (T 3 ) ≤ C(ε) f n -f n-1 2 L ∞ ((0,T )×T 3 ×R 3 ) + w n+1 2 L 2 (T 3 ) + w n 2 L 2 (T 3 ) . ( 29
)
The quantity f n -f n-1 2 L ∞ ((0,T )×T 3 ×R 3 ) has now to be estimated. We recall that, for all n ∈ N and (t, x, v), f n (t, x, v) = e 3t f ε in (χ n (t, x, v)) and that f ε in has been mollified. Consequently, we have (30)

f n -f n-1 L ∞ ((0,T )×T 3 ×R 3 ) ≤ C(ε, T ) χ n -χ n-1 L ∞ ((0,T )×T 3 ×R 3 )
. Definition (26)-( 27) of the characteristics χ n = (x n , v n ) and estimate (24) imply that, for any t,

(χ n -χ n-1 )(t) L ∞ (T 3 ×R 3 ) ≤ C(ε) t 0 (θ ε ⋆ u n -θ ε ⋆ u n-1 )(s) L ∞ (T 3 ) ds + t 0 (1 + ||(θ ε ⋆ u n )(s)|| W 1,∞ (T 3 ) ) (χ n -χ n-1 )(s) L ∞ (T 3 ×R 3 ) ds .
Thus, for any t, we can write

(χ n -χ n-1 )(t) L ∞ (T 3 ×R 3 ) ≤ C(ε) t 0 w n (s) L 2 (T 3 ) ds + t 0 (χ n -χ n-1 )(s) L ∞ (T 3 ×R 3 ) ds .
Using Gronwall's lemma, we obtain, for any t ≤ T , (31

) (χ n -χ n-1 )(t) L ∞ (T 3 ×R 3 ) ≤ C(ε, T ) t 0 w n (s) L 2 (T 3 ) ds.
Consequently, with (30) and ( 31), (29) becomes

|T 2 | ≤ C(ε, T ) t 0 w n 2 (s) L 2 (T 3 ) ds + w n+1 2 L 2 (T 3 ) + w n 2 L 2 (T 3 ) .
We finally obtain

1 2 d dt w n+1 2 L 2 (T 3 ) + 1 2 ∇ x w n+1 2 L 2 (T 3 ) ≤ C(ε, T ) t 0 w n (s) 2 L 2 (T 3 ) ds + w n+1 2 L 2 (T 3 ) + w n 2 L 2 (T 3 ) .
Therefore, integrating the previous inequality with respect to t, and using

t 0 s 0 w n (r) 2 L 2 (T 3 ) dr ds ≤ T t 0 w n (r) 2 L 2 (T 3 ) dr,
we obtain, for any t, 0 ≤ t ≤ T ,

w n+1 (t) 2 L 2 (T 3 ) ≤ C(ε, T ) t 0 w n (s) 2 L 2 (T 3 ) ds + t 0 w n+1 (s) 2 L 2 (T 3 ) ds .
Lemma 3 then implies that there exists a constant K(ε, T ) > 0, which does not depend on n, such that

w n 2 L ∞ (0,T ;L 2 (T 3 )) ≤ K n T n n! .
Consequently, the sequence (u n ) n∈N converges in L ∞ (0, T ; L 2 (T 3 )). That also implies, thanks to (30) and (31), that (f n ) n∈N converges in L ∞ ((0, T ) × T 3 × R 3 ).

Remark 4. Note that, at this stage, to study the regularized and uncoupled system, we used the high particle velocity cutoff function to control the force applied by the spray on the fluid. Moreover, the space regularization of the fluid velocity in the Vlasov equation enables to apply the method of characteristics and obtain a strong solution (this will also be true at the limit as n → ∞ for any fixed ε > 0). Finally, the space regularization of the convection fluid velocity simplifies the study of the fluid equations by making the convection term easier to deal with.

3.2.

Asymptotics with respect to n. As ε > 0 is still fixed, we can pass to the limit when n goes to infinity. We already know that (u n ε ) n converges to u ε strongly in L ∞ (0, T ; L 2 (T 3 )) and, thanks to (24), weakly in L 2 (0, T ; H 1 (T 3 )), and that (

f n ε ) n strongly converges to f ε in L ∞ ((0, T ) × T 3 × R 3 ). Thus, since f n ε belongs to C 0 ([0, T ] × T 3 × R 3 ), so is f ε .
These convergences enable us to pass to the limit in n in the following weak formulation, for all t, for all

φ ∈ C 1 ([0, T ] × T 3 × R 3 ):
(32)

T 3 ×R 3 f n ε (t) φ(t) dx dv - t 0 T 3 ×R 3 f n ε [∂ t φ + v • ∇ x φ + (θ ε ⋆ u n ε -v) • ∇ v φ] dx dv ds = T 3 ×R 3 f ε in φ(0) dx dv
Remark 5. Note that here the test functions are not necessarily with compact support in v. Indeed, since

(χ n ε ) n converges to χ ε in L ∞ ((0, T ) × T 3 × R 3 ) when n goes to +∞, we check (33) f ε (t, x, v) = e 3t f ε in (χ ε (t, x, v)
). Since we assumed that f ε in has been chosen with a compact support in v, it is also the case of

f ε (t, x, •) for all (t, x) ∈ [0, T ] × T 3 .
Note that here the solution f ε is in fact a strong solution.

Furthermore, for all ψ ∈ C 1 ([0, T ] × T 3 ) such that ∇ x • ψ = 0, for a. e. t,

t 0 T 3 ∂ t u n+1 ε • ψ dx ds + t 0 T 3 (θ ε ⋆ u n+1 ε ) • ∇ x u n+1 ε ψ dx ds + t 0 T 3 ∇ x u n+1 ε : ∇ x ψ dx ds = - t 0 T 3 ×R 3 f n ε (u n ε -v) γ ε (v) • ψ dx dv ds + t 0 T 3 F ε ext • ψ dx ds.
This yields ( 16) and ( 17). The limit function (f ε , u ε ) thus solves the (variational formulation of the) regularized system ( 11)-( 15).

This ends the proof of Proposition 2.

We now derive uniform bounds in ε for (f ε , u ε ). In a first step, these estimates only hold locally in time (but they only depend on T , u in 2

L 2 (T 3 ) + |M 2 f in | + F ext 2 L 2 (0,T ;L 2 (T 3 ))
and f in L ∞ (T 3 ×R 3 ) : this will be of utmost importance in the last part of the proof of Thm. 1). This constitutes the 3.3. Proof of Proposition 1. A similar argument to the one used in the proof of Proposition 2 ensures that f ε at least belongs to

Uniform estimates for the sequence

(f ε , u ε ). The fact that f ε in L ∞ (T 3 ×R 3 ) ≤ f in ∞ (T 3 ×R 3 ) together with equality (33) implies (19). Therefore, (f ε ) is uniformly bounded in L ∞ ((0, T ) × T 3 × R 3 ).
C 1 ([0, T ] × T 3 × R 3 ). That allows us to write (34) d dt M 2 f ε + 2M 2 f ε ≤ 2 T 3 |θ ε ⋆ u ε (t, x)| m 1 f ε (t, x) dx.
To estimate m 1 f ε in terms of second order moments of f ε , we need the following Lemma 1. Let β > 0 and g be a nonnegative function in L ∞ ((0, T ) × T 3 × R 3 ), such that m β g(t, x) < +∞, for a. e. (t, x). The following estimate holds for any α < β:

m α g(t, x) ≤ 4 3 π g(t, x, •) L ∞ (R 3 ) + 1 m β g(t, x) α+3 β+3 , a. e. (t, x).
the same kind of result can be found in [START_REF] Mellet | Global weak solutions for a Vlasov-Fokker-Planck/Navier-Stokes system of equations[END_REF] or [START_REF] Hamdache | Global existence and large time behaviour of solutions for the Vlasov-Stokes equations[END_REF].

Proof of Lemma 1. We write, for a. e. (t, x),

m α g(t, x) = |v|≤R |v| α g + |v|≥R |v| α g ≤ 4 3 πR α+3 g(t, x, •) L ∞ (R 3 ) + 1 R β-α m β g(t, x). By choosing R = m β g(t, x) 1 β+3
, we find the required result. This ends the proof of Lemma 1.

Applying Lemma 1, we obtain, thanks to (19),

0 ≤ m 1 f ε (t, x) ≤ C m 2 f ε (t, x) 4/5 , which yields (35) m 1 f ε (t) L 5/4 (T 3 ) ≤ C M 2 f ε (t) 4/5 ,
where C only depends on T and f in L ∞ (T 3 ×R 3 ) . This estimate, together with (34), gives

d dt M 2 f ε + 2M 2 f ε ≤ θ ε ⋆ u ε L 5 (T 3 ) m 1 f ε L 5/4 (T 3 ) ≤ C u ε L 5 (T 3 ) (M 2 f ε ) 4/5 .
Setting

H(t) = M 2 f ε (t) 1/5
, we obtain

H ′ + H ≤ C u ε L 5 (T 3 ) .
Since M 2 f ε ≥ 0, and assuming that M 2 f ε in ≤ C M 2 f in , we have, for a. e. t,

H(t) ≤ C 1 + t 0 u ε (s) L 5 (T 3 ) ds ,
and thus

(36) M 2 f ε (t) C 1 + t 0 u ε (s) L 5 (T 3 ) ds 5 ,
where C only depends on the initial data.

Next, we take care of the fluid velocity. The function u ε is an admissible test function for the variational formulation (17). Indeed, by density, (17) remains valid for any ψ ∈ L 2 (0, T ; H 1 (T 3 )). After integration by parts, we have

1 2 d dt u ε 2 L 2 (T 3 ) + ∇ x u ε 2 L 2 (T 3 ) = T 3 ×R 3 f ε u ε • (v -u ε ) γ ε (v) dx dv + T 3 F ext • u ε . Since γ ε ≥ 0, that implies 1 2 d dt u ε 2 L 2 (T 3 ) + ∇ x u ε 2 L 2 (T 3 ) ≤ T 3 |u ε | m 1 f ε dx + F ext L 2 (T 3 ) u ε L 2 (T 3 ) ≤ u ε L 5 (T 3 ) m 1 f ε L 5/4 (T 3 ) + F ext L 2 (T 3 ) u ε L 2 (T 3 ) .
Applying Young's inequality and remembering (35) and (36), we get

(37) 1 2 d dt u ε 2 L 2 (T 3 ) + ∇ x u ε 2 L 2 (T 3 ) ≤ C u ε L 5 (T 3 ) 1 + t 0 u ε L 5 (T 3 ) 4 + C F ext 2 L 2 (T 3 ) + C u ε 2 L 2 (T 3 ) .
Next, using the continuous injection of H 1 (T 3 ) in L 5 (T 3 ) (in dimension less or equal than 3), we obtain

1 2 d dt u ε 2 L 2 (T 3 ) + ∇ x u ε 2 L 2 (T 3 ) ≤ 1 2 ∇ x u ε 2 L 2 (T 3 ) + C u ε 2 L 2 (T 3 ) + C + C t 0 ∇ x u ε L 2 (T 3 ) + u ε L 2 (T 3 ) 8 + C F ext 2 L 2 (T 3 ) . ( 38
) Using t 0 s 0 u ε (r) 8 L 2 (T 3 ) dr ds ≤ T t 0 u ε (r) 8 L 2 (T 3 ) dr,
and Hölder's inequality, we integrate (38) to obtain

1 2 u ε 2 L 2 (T 3 ) (t) + 1 4 t 0 ∇ x u ε 2 L 2 (T 3 ) ≤ 1 4 t 0 u ε 2 L 2 (T 3 ) + 1 2 u ε in 2 L 2 (T 3 ) + C + C t 0 s 0 ∇ x u ε 2 L 2 (T 3 ) 4 + C t 0 u ε 8 L 2 (T 3 ) + C t 0 F ext 2 L 2 (T 3 ) .
A nonlinear Gronwall lemma thus implies that there exist a time T * > 0 and a constant C such that, for almost every t ≤ T * , (39)

u ε (t) 2 L 2 (T 3 ) + t 0 ∇ x u ε (s) 2 L 2 (T 3 ) ds ≤ C.
We recall here that

C := C(T, T * , u in 2 L 2 (T 3 ) + |M 2 f in | + F ext 2 L 2 (0,T ;L 2 (T 3 )) , f in L ∞ (T 3 ×R 3 )
). Finally, estimates (36) and (39) imply that (40)

M 2 f ε (t) ≤ C, a. e. t ≤ T * .
3.4. Asymptotics with respect to ε. In order to pass to the limit in the weak formulations ( 16)-( 17), we need to have some compactness properties for the sequence (f ε , u ε ) ε>0 . In fact, proving that (u ε ) ε>0 is compact in L p (0, T ; L 2 (T 3 )) for some p, is enough for the asymptotics ε → 0. From (17), it is easy to check that (∂ t u ε ) ε>0 is uniformly bounded in L 3/2 (0, T * ; V ′ ). As a matter of fact, once again, the only two terms we have to take care of are the convection and the drag force ones.

For the convection term, we use standard arguments:

T 3 ((θ ε ⋆ u ε ) • ∇ x ) u ε • ψ = - T 3 ((θ ε ⋆ u ε ) • ∇ x ) ψ • u ε . But (u ε ) ε>0 is uniformly bounded in L ∞ (0, T * ; L 2 (T 3 ))∩L 2 (0, T * ; H 1 (T 3
)), which is continuously embedded in L 3 (0, T * ; L 4 (T 3 )). Thus

T * 0 T 3 ((θ ε ⋆ u ε ) • ∇ x ) u ε • ψ ≤ C ∇ x ψ L 3 (0,T * ;L 2 (T 3 )) .
Thus the form

ψ → T * 0 T 3 ((θ ε ⋆ u ε ) • ∇ x ) u ε • ψ is bounded in L 3/2 (0, T * ; V ′ ).
The term coming from the drag force can be estimated as follows:

T * 0 T 3 ×R 3 f ε (u ε -v)γ ε (v)ψ ≤ C u ε L 2 (0,T * ;L 6 (T 3 )) ψ L 2 (0,T * ;L 6 (T 3 )) m 0 f ε L ∞ (0,T * ;L 3/2 (T 3 )) + ψ L 2 (0,T * ;L 5 (T 3 )) m 1 f ε L 2 (0,T * ;L 5/4 (T 3 )) .
But, thanks to Lemma 1, we have

m 0 f ε ≤ C(m 2 f ε ) 3/5 , m 1 f ε ≤ C(m 2 f ε ) 4/5 ,
and we proved that

M 2 f ε L ∞ (0,T * ) ≤ C.
Therefore, the sequence (m 0 f ε ) ε>0 and (m 1 f ε ) ε>0 are respectively bounded in L ∞ (0, T * ; L 5/3 (T 3 × R 3 )) and L ∞ (0, T * ; L 5/4 (T 3 × R 3 )).

Moreover, (u ε ) is bounded in L 2 (0, T * ; H 1 (T 3 )) and thus in L 2 (0, T * ; L 6 (T 3 )).

Eventually, we can write

T * 0 T 3 ×R 3 f ε (u ε -v) γ ε (v) ψ ≤ C ψ L 2 (0,T * ;L 6 (T 3 )) ≤ C ψ L 2 (0,T * ;H 1 (T 3 )) ,
and the drag force

R 3 f ε (u ε -v) γ ε (v) dv ε>0 is uniformly bounded in L 2 (0, T * ; V ′ ).
That ensures that (∂ t u ε ) ε>0 is uniformly bounded in L 3/2 (0, T * ; V ′ ).

Thanks to those bounds, there exists (f, u) ∈ L ∞ ((0, T * )×T 3 ×R 3 )×(L ∞ (0, T * ; L 2 (T 3 ))∩L 2 (0, T * ; H 1 (T 3 ))× C 0 ([0, T * ]; V ′ ) such that, up to a subsequence, (f ε ) weakly converges to f in w * -L ∞ ((0, T * ) × T 3 × R 3 ) and (u ε ) weakly converges to u in w-L 2 (0, T * ; w-H 1 (T 3 )) and in w * -L ∞ (0, T * ; w-L 2 (T 3 )). Thanks to Aubin's compactness lemma, since (∂ t u ε ) ε>0 is uniformly bounded in L 3/2 (0, T * ; V ′ ), we can state that (u ε ) strongly converges to u in L 2 (0, T * ; L 2 (T 3 )).

These convergences allow us to pass to the limit in the following weak formulations:

-

T * 0 T 3 ×R 3 f ε [∂ t φ + v • ∇ x φ + (θ ε ⋆ u ε -v) • ∇ v φ] dx dv ds = T 3 ×R 3 f ε in φ(0) dx dv, , for all φ ∈ C 1 ([0, T * ] × T 3 × R 3
) with compact support in v and such that φ(T * ) = 0, and

- T 3 u ε (t) • ψ(t) dx ds + t 0 T 3 ((θ ε ⋆ u ε ) • ∇ x ) u ε • ψ dx ds + t 0 T 3 ∇ x u ε : ∇ x ψ dx ds = - t 0 T 3 ×R 3 f ε (u ε -v) γ ε (v) • ψ dx dv ds + t 0 T 3 F ε ext • ψ dx ds, a. e. t ≤ T * , for all ψ ∈ C 1 ([0, T * ] × T 3 ) such that ∇ x • ψ = 0.
We then obtain a solution satisfying ( 8)-( 9).

This ends the proof of Proposition 1.

3.5. Global existence. We just derived the existence locally in time of a weak solution (f, u). In this section, we prove that this solution is, in fact, global in time, which constitutes the last step of the proof of Thm. 1. We start with the Lemma 2. The solution (f, u) obtained above satisfies

(41) 1 2 M 2 f (t) + 1 2 u(t) 2 L 2 (T 3 ) + t 0 ∇ x u(s) 2 L 2 (T 3 ) ds + t 0 T 3 ×R 3 f |u -v| 2 dx dv ds ≤ (1 + e T + T e T ) 1 2 u in 2 L 2 (T 3 ) + 1 2 |M 2 f in | + 1 2 F ext 2 L 2 (0,T ;L 2 (T 3 )) , a.e. t ≤ T * , and
f L ∞ ((0,T * )×T 3 ×R 3 ) ≤ e 3T f in L ∞ (T 3 ×R 3 ) .
Proof of Lemma 2. We first note that the second inequality is directly obtained from the definition of f , as the weak limit of (f ε ) for the w * -L ∞ ((0, T * ) × T 3 × R 3 ) topology. For the first inequality, we shall study the energy estimate satisfied by (f ε , u ε ) obtained by choosing once again φ = |v| 2 /2 as a test function in (16), and u ε as a test function in (17):

(42) 1 2 M 2 f ε (t) + 1 2 u ε (t) 2 L 2 (T 3 ) + t 0 ∇ x u ε (s) 2 L 2 (T 3 ) ds + t 0 T 3 ×R 3 f ε |u ε -v| 2 dx dv ds = 1 2 M 2 f ε in + 1 2 u ε in 2 L 2 (T 3 ) + t 0 T 3 F ε ext • u ε dx ds + R ε (t), where R ε (t) = t 0 T 3 ×R 3 f ε |u ε | 2 (1 -γ ε (v)) dx dv ds - t 0 T 3 ×R 3 f ε u ε • v(1 -γ ε (v)) dx dv ds + t 0 T 3 ×R 3 f ε (u ε -θ ε ⋆ u ε ) • v dx dv ds. Let R i ε , i = 1, 2, 3
, denote each of the three terms appearing in R ε and h ε = f ε (1 -γ ε ). We first have, for a. e. t,

0 ≤ R 1 ε (t) ≤ m 0 h ε L ∞ (0,T * ;L 3/2 (T 3 )) u ε 2 L 2 (0,T * ;L 6 (T 3 )) .
Lemma 1 implies that

m 0 h ε (t) L 3/2 (T 3 ) ≤ C M 3/2 h ε (t) 3/2 , and
M 3/2 h ε (t) = T 3 ×R 3 |v| 3/2 f ε (1 -γ ε ) ≤ T 3 ×{|v|≥ 1 2 ε} |v| 3/2 f ε ≤ √ 2ε M 2 f ε (t).
Consequently, from the uniform bound obtained on (M 2 f ε ) ε>0 in L ∞ (0, T * ), we deduce that the sequence ( m 0 h ε L ∞ (0,T * ;L 3/2 (T 3 )) ) converges to 0 when ε goes to 0.

Finally, since (u ε ) ε>0 is uniformly bounded in L 2 (0, T * ; H 1 (T 3 )) ֒→ L 2 (0, T * ; L 6 (T 3 )), we have proven that (R 1 ε ) goes to 0 when ε → 0. To treat the second term R 2 ε , we proceed in the same way. For a. e. t, we have 0 ≤ R 2 ε (t) ≤ m 1 h ε L ∞ (0,T * ;L 6/5 (T 3 )) u ε L 2 (0,T * ;L 6 (T 3 )) .

Lemma 1 implies that m 1 h ε L 5/6 (T 3 ) (t) ≤ C M 9/5 h ε (t) 5/6 , and

M 9/5 h ε (t) = T 3 ×R 3 |v| 9/5 f ε (1 -γ ε ) ≤ T 3 ×{|v|≥ 1 2 ε} |v| 9/5 f ε ≤ (2ε) 1/5 M 2 f ε (t).
Consequently, from the uniform bound obtained on (M 2 f ε ) ε>0 in L ∞ (0, T * ), we deduce that the sequence ( m 1 h ε L ∞ (0,T * ;L 5/6 (T 3 )) ) converges to 0 when ε goes to 0.

Finally, since (u ε ) ε>0 is uniformly bounded in L 2 (0, T * ; H 1 (T 3 )) ֒→ L 2 (0, T * ; L 6 (T 3 )), we have proven that (R 2 ε ) goes to 0 when ε → 0.

We now take care of the last term R 3 ε . Let us write, for any t,

R 3 ε (t) = t 0 T 3 ×R 3 f ε (u ε -θ ε ⋆ u ε ) • vγ µ (v) dx dv ds + t 0 T 3 ×R 3 f ε (u ε -θ ε ⋆ u ε ) • v(1 -γ µ (v)) dx dv ds,
where µ > 0.

From arguments similar to those used for the second term, we can check that 0 ≤

t 0 T 3 ×R 3 f ε (u ε -θ ε ⋆ u ε ) • v(1 -γ µ (v)) dx dv ds ≤ C (2µ) 1/5 ,
where C depends neither on ε nor on µ.

Furthermore, since (f ε ) ε>0 is uniformly bounded in L ∞ ((0, T * ) × T 3 × R 3 ), we have, for a. e. t,

0 ≤ t 0 T 3 ×R 3 f ε (u ε -θ ε ⋆ u ε ) • v γ µ (v) dx dv ds ≤ C u ε -θ ε ⋆ u ε L 1 (0,T * ;L 1 (T 3 )) .
Consequently, since u ε strongly converges in L 2 (0, T * ; L 2 (T 3 )), the sequence

t 0 T 3 ×R 3 f ε (u ε -θ ε ⋆ u ε ) • v γ µ (v) dx dv ds ε>0
goes to 0 when ε → 0. That ensures that (R 3 ε ) goes to 0 when ε → 0. Eventually we have proven that (R ε ) converges to 0.

For the force term, we note that, by Young's inequality,

t 0 T 3 F ε ext • u ε ≤ 1 2 t 0 u ε 2 L 2 (T 3 ) + 1 2 t 0 F ε ext 2 L 2 (T 3 ) .
Equality (42) implies that

1 2 M 2 f ε (t) + 1 2 u ε (t) 2 L 2 (T 3 ) + t 0 ∇ x u ε (s) 2 L 2 (T 3 ) ds + t 0 T 3 ×R 3 f ε |u ε -v| 2 dx dv ds ≤ 1 2 M 2 f ε in + 1 2 u ε in 2 L 2 (T 3 ) + 1 2 t 0 u ε (s) 2 L 2 (T 3 ) ds + 1 2 t 0 F ε ext (s) 2 L 2 (T 3 ) ds.
We then apply Gronwall's lemma to get For the remaining terms, the weak convergences of (f ε , u ε ) ε>0 and the strong convergence of (u ε ) ε>0 enable us to obtain (41).

1 2 M 2 f ε (t) + 1 2 u ε (t) 2 L 2 (T 3 ) + t 0 ∇ x u ε (s) 2 L 2 (T 3 ) ds + t 0 T 3 ×R 3 f ε |u ε -v| 2 dx dv ds ≤ (1 + e t ) 1 2 u ε in 2 L 2 (T 3 ) + 1 2 |M 2 f ε in | +
This ends the proof of Lemma 2.

Thanks to this lemma, we can now end the proof of Thm. 1.

3.6. Extension of the solution. Let T > 0 fixed, and consider a set of admissible initial data and external force. Let us choose K such that (1 + e T + T e T ) 1 2 u in

  Consequently, and thanks to Remark 5, |v| 2 is an admissible test function in the variational formulation (16).

1 2 ( 1 + t e t ) t 0 F ε ext 2 L 2 M 2

 1210222 (0,T ;L 2 (T 3 )) ds, and let ε → 0 in the equation. It easy to check thatM 2 f (t) ≤ lim inf µ→0 (f γ µ )(t), and M 2 (f γ µ ) L ∞ (0,T * ) ≤ lim inf ε→0 M 2 (f ε γ µ ) L ∞ (0,T * ) ≤ lim inf ε→0 M 2 f ε L ∞ (0,T * ) .

Acknowledgements

This work has been partially supported by the French national research agency (ANR) grant # 0003205 Transport et transfert via le système pulmonaire humain. L. Boudin wishes to thank the INRIA research center of Paris -Rocquencourt, and the Reo project-team for the temporary researcher position (délégation) he got there.

Appendix A. Gronwall lemma

In this section, we state and prove the Gronwall lemma we used several times in this work. Lemma 3. Let T > 0 and consider a sequence (a n ) n∈N of nonnegative continuous functions on [0, T ]. Assume that (a n ) satisfies, for any n,

where A, B and C are nonnegative constants.

If A = 0, there exists a constant K ≥ 0 such that

If A > 0, there exists a constant K ≥ 0 depending on A, B, C such that

Proof of Lemma 3: If both A = C = 0, we choose

and the induction is immediate to get (44).

If A > 0 and C = 0, we choose

and the induction is once again immediate to get (45).

If C > 0, we can apply one of the two previous cases to the sequence (α n ) defined, for any n, by That ends the proof of Lemma 3.


Nicolas Ray

Bruno Vallet

Wan-Chiu Chiu Li

Bruno Lévy

N-Symmetry Direction Field Design

Keywords: I.3.7 [Computer Graphics]: Three-Dimensional Graphics and Realism-Color, shading, shadowing, and texture, I.3.5 [Computer Graphics]: Computational Geometry and Object Modeling, G.1.6 [Numerical Analysis]: Optimization, J.6 [Computer Aided Engineering] Vector eld design, N-symmetry direction eld, topology

INTRODUCTION

In Computer Graphics, a wide class of applications requires a smooth direction field over a surface to be defined. For instance, such a direction field was used in [START_REF] Hertzmann | Illustrating smooth surfaces[END_REF] to place hatch strokes in NPR applications, in [START_REF] Turk | Texture synthesis on surfaces[END_REF]] to steer the orientation of features in texture synthesis, or in [START_REF] Alliez | Anisotropic polygonal remeshing[END_REF]] to remesh a surface with cells aligned with a pair of direction fields. A closer look reveals that these applications use objects of higher symmetry than simple direction fields, i.e. objects invariant by rotation of π or π/2. For instance, the "natural" result one would expect when designing two smooth orthogonal direction fields on the cube-like object shown on Figure 1(left) connects the two direction fields around the corners of the cube with a π/2 rotation. To the best of our knowledge, this problem was first mentioned in the Computer Graphics field by Hertzmann and Zorin in their non-photorealistic paper [START_REF] Hertzmann | Illustrating smooth surfaces[END_REF], where they intuitively introduce the notion of "cross-field". They mentioned that (1) characterizing singularities (points where a direction cannot be defined) and especially (2) controlling their number, type and placement are key issues and interesting mathematical problems involving an analog of Euler formula.

A solution to issue (1) appears in [START_REF] Li | Representing higher-order singularities in vector fields on piecewise linear surfaces[END_REF]] which proposes a data structure for direction fields based on controlling its behavior between discretization points through integers called "period jump" which characterized singularities. The authors then noticed that allowing period jumps to be multiple of 1/2 or 1/4 allowed to handle the direction symmetries, thus singularities of fractionary indices (see Figure 1). They used this to build an interpolation scheme for such symmetric fields, allowing to visualize arbitrary singularities, but examples had to be hand made or converted from an existing vector field (in analytic form or in a classical data structure).

This paper provides a solution to the much more complex issue (2). We start by formalizing N-symmetry directions in the continuous setting as the sets of directions invariant by 2π/N rotations, such that 1,2 and 4-symmetry direction fields correspond to classical direction fields, line fields and cross fields. We then provide a comprehensive study of their topology, and link it with the data structure and notion of period jumps from [START_REF] Li | Representing higher-order singularities in vector fields on piecewise linear surfaces[END_REF]]. This

• 3 allows finally to derive an efficient algorithm to generate a smooth N-symmetry direction field that interpolates a user-defined set of singularities (number, type and placement), as well as a set of geometric constraints on directions.

Before entering the heart of the matter, we give an overview of the previous work in direction field processing.

Previous work

Our work explains how to design smooth N-symmetry direction fields on mesh surfaces under both geometric and topological constraints, based on the data structure presented in [START_REF] Li | Representing higher-order singularities in vector fields on piecewise linear surfaces[END_REF]]. This problem has already been tackled for several applicative purposes:

-non-photorealistic rendering: In [START_REF] Hertzmann | Illustrating smooth surfaces[END_REF], a smooth cross field is constructed to place strokes for a non-photorealistic rendering application. The construction method shares some common points with ours: the variables used to represent the directions are angles measured relative to a given arbitrary direction field. This paper follows the possibilities of future work and mathematical investigations suggested at the end of their paper.

-texture synthesis: In lapped textures [START_REF] Praun | Lapped textures[END_REF]] and texture particles [START_REF] Dichler | Texture particles[END_REF], a direction field is used to control the orientations and sizes of texture patches distributed over the object. In [START_REF] Praun | Lapped textures[END_REF]], a smooth direction field is generated through radial basis functions, with geodesic distances computed over the surface. The method presented in [START_REF] Turk | Texture synthesis on surfaces[END_REF]] operates at texel level, by using the direction field to steer the anisotropy of a texture synthesizer. The direction smoothing procedure they used is inspired by [START_REF] Gortler | The lumigraph[END_REF], which is based on a multi-resolution Laplacian smoother. A similar procedure is described in [START_REF] Ohtake | Adaptive smoothing tangential direction fields on polygonal surfaces[END_REF], with the addition of non-linear weights that preserve important direction field discontinuities. In [START_REF] Wei | Texture synthesis over arbitrary manifold surfaces[END_REF], 2 and 4-symmetry direction fields are used to steer synthesizing using 2 and 4-symmetry texture samples. Finally, [START_REF] Zelinka | Jump map-based interactive texture synthesis[END_REF] steer their texture generation method using a direction field defined as the gradient of a fair Morse function (it has the same singular points as the function). Based on the study of the Morse complex of smooth harmonic functions [START_REF] Ni | Fair morse functions for extracting the topological structure of a surface mesh[END_REF], this allows the user to control the number and configuration of singularities.

-remeshing and global parameterization: In [START_REF] Alliez | Anisotropic polygonal remeshing[END_REF]], a method is proposed to generate quadrangles aligned with two orthogonal direction fields, obtained by smoothing an estimation of the curvature tensor. The refinement presented in [START_REF] Marinov | Direct anisotropic quad-dominant remeshing[END_REF] operates without a global parameterization and uses, in a certain sense, an explicit version of Ohtake et al.'s non-linear weights [2001] to preserve important features. More recently, [START_REF] Ray | Periodic global parameterization[END_REF]] presented a quad remeshing method based on a periodic global parameterization following a direction field. The global parameterization method presented in [START_REF] Gu | Global conformal surface parameterization[END_REF]] also follows a direction field defined through a pair of holomorphic functions.

Our paper shares some common points with [START_REF] Zelinka | Jump map-based interactive texture synthesis[END_REF] and [START_REF] Gu | Global conformal surface parameterization[END_REF], in particular, the ability of controlling singularities. The main difference is that in these two methods, the direction field is defined to be the gradient of a scalar field, so it is

•

necessarily curl-free. Moreover, we handle symmetry so we can represent a wider class of singularities, with arbitrary indices.

The works that handle symmetries usually focus on the specific case of "cross-fields"

v 1 , v 2 = R n ( v 1 , π/2),
where v 1 and v 2 can be swapped (see the cube-like object of Figure 1), most of the smoothing algorithms mentioned above use a local relaxation procedure, updating values v 1 (p i ) and v 2 (p i ) vertex by vertex. During these computations, v 1 (p i ) will be influenced either by v 1 (p j ) or v 2 (p j ). The algorithm chooses the direction nearest to v 1 (p i ). As a consequence, singularities may appear without control, and the convergence is slow. In contrast, based on a combined topological analysis of both the surface and the direction field, we derive a global formulation of the problem as a quadratic form minimization, yielding a more efficient optimization procedure (conjugate gradient), which computes directly an optimal solution.

Some more recent papers directly address problems related with direction field processing. For instance, in [START_REF] Polthier | Identifying vector fields singularities using a discrete hodge decomposition[END_REF] and [START_REF] Tong | Discrete multiscale vector field decomposition[END_REF]] a procedure for computing the Hodge decomposition of a direction field is described. This decomposition isolates some features of the field, and makes it possible to filter or to enhance them. In [START_REF] Tricoche | Topology simplification of symmetric, secondorder 2D tensor fields[END_REF]], a method is presented to simplify the topology of symmetric, second order 2D direction fields. Recent work on vector field design include [START_REF] Wang | Edge subdivision schemes and the construction of smooth vector fields[END_REF] where subdivision schemes are used to define bases for discrete differential 0-and 2-forms, and [START_REF] Fisher | Design of tangent vector fields[END_REF]] that uses Whitney forms and the Hodge decomposition.

Both approaches allow to design vector fields under user constraints, but do not provide any guarantee that new singularities will not appear, nor tackle the problem of symmetric fields. Finally, a complete toolkit for interactive vector field design based on moving singularities and pair cancellation is presented in [START_REF] Zhang | Vector field design on surfaces[END_REF], and generalized to tensors in [START_REF] Zhang | Interactive design and visualization of tensor fields on surfaces[END_REF]], then to N-rotational symmetry fields (N-RoSy) in [START_REF] Palacios | Rotational symetry field design on surfaces[END_REF]. As in this paper, [START_REF] Palacios | Rotational symetry field design on surfaces[END_REF]] describes a tool to design a smooth N-symmetry direction field from a user-defined set of singularities. However, they start by generating a base field from the user constraints without controlling the topology, then "repair" the field a posteriori by moving or removing (by pair cancellation) the undesired singularities that appeared. In contrast, we ensure topological constraints in the design process a priori, by construction.

Contributions

-We provide a comprehensive study of the topology of N-symmetry direction fields by introducing the concept of turning number. This generalizes the notion of singularity index (see Sections 2.5 and 2.6);

-We prove an analog of the Poincaré-Hopf theorem based on turning numbers, adapt it to N-symmetry direction fields, and generalize it to manifolds with borders. We also prove a theorem relating turning numbers to direction field topology. With both theorems we can fully explain the links between the topology of a surface and the topology of a direction field defined over it. (see Appendix A);

-We apply this theory to the period jumps based structure presented in [START_REF] Li | Representing higher-order singularities in vector fields on piecewise linear surfaces[END_REF]] for representing direction fields. Thus we explain how period jumps control the topology of the direction field (see Section 3); -We describe an algorithm for constrained N-symmetry direction field design. From a user-defined set of singularities and an optional set of points with fixed directions, our algorithm constructs a smooth direction field that satisfies all the constraints. If the indices of the user-defined singularities sum up to the Euler characteristic (2g -2), the construction is guaranteed to generate no other singularity (see Section 5). To our knowledge, this is the first algorithm that performs provably correct direction field design.

DIRECTION FIELDS ON SURFACES WITH BORDERS

This section presents the fundamental tools for studying direction fields defined over surfaces with borders. more specifically, we focus on their topology. Topology is the study of properties that are invariant by continuous deformations (without cutting or gluing anything), called homotopies. In other words, topology tries to answer the following question: under which conditions are two objects homotopic (i.e. can be continuously transformed one into another) ? For oriented surfaces with borders, the answer is that they need to have the same genus g (number of handles) and number of borders b. In other terms, if S 1 and S 2 are two surfaces with borders, we have S 1 ≡ t S 2 ⇔ g(S 1 ) = g(S 2 ) and b(S 1 ) = b(S 2 ), where ≡ t denotes the homotopic equivalence (t stands for topology). What is even more interesting is the structure of the set of homotopy classes (classes of all objects with same topology). For oriented surfaces with borders again, this set is isomorphic to N 2 , since any pair of non negative integers (g, b) can be associated to the class of all surfaces with genus g and b borders. For this reason, we say that these two integers are the two topological degrees of freedom (TDoF) of surfaces. This section addresses the same questions for Nsymmetry direction fields defined over a 2-manifold, which means that we aim at studying under which conditions two N-symmetry direction fields are homotopic. In other words, we want to exhibit the TDoF of these direction fields. To answer these questions, we will introduce the concept of turning numbers of cycles in direction fields. As we will show, the TDoF of a direction field are the turning numbers of a homology basis of cycles, such that these turning numbers capture all the topological structure of the field. Moreover, turning numbers generalize singularities, and we notice that singularities alone do not control all of the field topology.

The intuitive idea behind turning numbers is that when following a closed cycle, a Nsymmetry direction might do an arbitrary number of N th of turns before coming back to its starting point. Imagine you are traveling on earth along a cycle with a compass giving you the north direction, then you can count the number of turns of the compass while following the cycle. If you turned around a tree (there are no trees at earth poles...) you will get one turn, but if you followed the equator or turned around a pole, you will get zero turns. We can do the same on any surface, and with any direction field defined over it. We call this quantity the "turning number" of the field along the cycle (it depends on both the cycle and the direction field). We will show that turning numbers capture the topology of the direction field, because two direction fields are homotopic iff they have the same turning numbers along the 2g + b -1 cycles of a homology basis (a basis for cycles on a 2 manifold). This shows that homotopy classes of direction fields are isomorphic to Z 2g+b-1 , or in other words that the topology of a direction field is entirely defined by 2g + b -1 integers. These integers are the TDoF (topological degrees of freedom) of direction fields.

In this section, we will first give the definitions of the objects that we will be using throughout this paper: surfaces (Section 2.1), cycles (Section 2.2) and direction fields (Section 2.3). Then we define the curvature for cycles and direction fields. We use these curvatures to formally define turning numbers, and exhibit their fundamental properties. We finally explain how turning numbers relate to field singularities.

Surfaces with borders

In this paper, we call surface (or 2-manifold) S a topological space where each point has a neighborhood homeomorphic to the plane or half plane. The surface S considered here is compact, connected and oriented, so that each point has a unique unit normal vector n.

The surface S has g handles and a boundary ∂ S which is the collection of its b borders, for some (non-negative) integers g called the genus and b called the number of borders of S .

Chains and Cycles on S

A chain γ on S is an oriented 1-manifold embedded in S (γ ⊂ S ). Because γ is oriented, it has a unique tangent vector t γ in each point which is also tangent to S . Using this tangent vector, along with the surface normal n, we can define a unique conormal vector n γ = n × t γ on the chain, which ensures that ( t γ , n γ , n) form a natural local orthonormal basis called the Darboux frame (see Figure 4). Notice that chains are not necessarily connected, so the term "set of chains" would be more appropriate (but more cumbersome). We define the following notions on chains (see Figure 3 left): -The reversal -γ of a chain γ is the chain with opposite orientation: t -γ =t γ . We use the notations γ 0 + γ 1 and γ 0γ 1 in place of γ 0 ∪ γ 1 and γ 0 ∪ -γ 1 as it is more practical to handle unions of cycles with various orientations.

-A cycle on S is a chain without endpoints (∂ γ = 0).

-We call ∂ the boundary operator, such that ∂ S is the subset of points of S with neighborhoods homeomorphic to the half plane. This subset is a cycle for which we can choose an orientation by requiring its conormal to point outwards S (this is why we represent the conormal instead of tangent in our figures).

-A chain γ is called exact if there exist a sub-manifold S of S such that γ = ∂ S. An exact chain is necessarily a cycle (a boundary cannot have endpoints).

-Two cycles γ 0 and γ 1 are homological iff γ 0γ 1 is exact (see Figure 5). We will denote ≡ l the homologic equivalence (for cycles), when ≡ t is the homotopic equivalence (for surfaces and direction fields). More formally, the homology for cycles is defined as the quotient set of closed chains (cycles) over exact chains. In particular the 0 of homology is the class of exact cycles, so γ is exact

⇔ γ ≡ l / 0 -H(S ) = {γ H i } i=1..n is called a homology basis on S if: -Linear independence: ∑ a i γ H i ≡ l / 0 ⇔ a 1 = ... = a n = 0 -Spanning: ∀γ ∃a ∈ Z n such that γ ≡ l ∑ n i=1 a i γ H i (

any cycle on S is homological to a formal sum of the basis cycles).

An important result from manifold topology is that a homology basis on S has 2g + b -1 basis vectors, where g is the genus, and b the number of borders of S .

Direction field

A unit tangent vector u on S is a vector that satisfies || u|| = 1 and u • n = 0. We call N-symmetry direction on S a set of N unit tangent vectors on S invariant by rotation • Fig. 5. γ 2 stays homological to γ 1 because γ 2γ 1 is exact(it is the border of a 2-manifold). Note that homological cycles might have a different number of connected components. of 2π/N around the normal. Hence, based on a unit tangent vector u 0 we can build a N-

symmetry direction d = { u k = R n ( u 0 , 2kπ/N)|k ∈ Z} where R n is the rotation around n.
We call direction field d on S a mapping that associates a N-symmetry direction d(P) to each point P ∈ S . In what follows, we will omit the term N-symmetry for brevity.

We call D N (S ) the set of differentiable direction fields, which are direction fields which curvature (see definition below) is defined almost everywhere. Two direction fields d 0 ,

d 1 ∈ D N (S ) are called homotopic if there exists a continuous function Γ : [0, 1] → D N (S ) such that Γ(0) = d 0 , Γ(1) = d 1 .
In other terms, two direction fields are homotopic (or have the same topology) if they can be continuously transformed one into another. As this is a relation of equivalence, we can define homotopy classes as the sets of all direction fields equivalent to a given one. Based on the intuition of what is a direction field in the plane, it would seem that all direction fields are homotopic, then there would be only one homotopy class. However, if the direction field is defined on a surface of arbitrary topology, we will show that it is no longer the case. Homotopy classes of direction fields can be characterized by the turning numbers of the field along some cycles. They correspond intuitively to the number of times the direction turns in a local Darboux frame while moving along the cycle. We are now going to define the notion of curvature for both cycles and direction fields, which will be used for a rigorous definition of the turning number.

Curvature

The curvature of a chain expresses the angular variation of its tangent along the chain. If we call s the arclength on a chain γ, we can define the curvature of γ using the decomposition:

∂ t γ ∂ s = κ γ n γ + κ S n κ S = ∂ t γ ∂ s • n κ γ = ∂ t γ ∂ s • n γ (1)
-κ S measures the normal curvature of S in direction t γ .

-κ γ measures the curvature of γ in the tangent plane of S .

-∂ t γ /∂ s • t γ = 0 by derivation of t γ • t γ = 1 ( t γ is a unit vector)

• 9

We can similarily define the curvature κ d ( t γ ) for the direction field d = { u k } in direction t γ as:

κ d ( t γ ) = ∂ u k ∂ s • u ⊥ k (2)
where u ⊥ k = n× u k is a unit vector orthogonal to u k in the tangent plane, such that ( u k , u ⊥ k , n) is an orthonormal basis (it is the Darboux frame of the streamlines of u k ). The curvature κ d is the same for all k ∈ N so it can be called the curvature of d in direction t γ . The quantity (κ dκ γ )ds can be interpreted as the infinitesimal angular variation dθ of the direction field in the Darboux frame when moving along the cycle (see Figure 4). We are now going to integrate this variation to define the turning numbers.

Turning number

The turning number of a direction field along a cycle corresponds to the number of rotations of the field along this cycle. Turning numbers are characteristic of homotopy classes of direction fields, hence of their topology.

Let γ PQ be a chain with endpoints P and Q, and θ γ (P) be an angle (in R) between d(P) and t γ (P). Then:

θ γ (Q) = θ γ (P) + γ PQ dθ = θ γ (P) + γ PQ (κ d -κ γ )ds (3) is an angle (in R) between d(Q) and t γ (Q).
Especially, for a cycle γ, we have P = Q, so we will have equality of the angles modulo 2π/N. Therefore, the quantity:

T d (γ) = 1 2π γ (κ d -κ γ )ds (4)
is necessarily an integer multiple of 1/N which we call the turning number of d along γ, and corresponds intuitively to the number of N th of turns d does in the Darboux frame along γ (see Figure 6).

In the case T d (γ) = 0, there exists a continuous scalar function θ defined on γ such that (κ dκ γ )ds = dθ . The scalar θ is defined up to a constant which can be chosen such that θ is an angle between t γ and d. In the case T d (γ) = 0, θ cannot be defined in a way to be continuous (we need 2π/N jumps).

Turning numbers have two fundamental properties (see Appendix A for a proof) which make them useful for studying direction field topology:

THEOREM 2.1 BOUNDARY TURNING NUMBER. Let S be a surface (2-manifold with borders) embedded in R 3 , then:

T d (∂ S) = -χ(S) (5)
where χ(S) = 2 -2g(S)b(S) is the Euler characteristic of S.

Theorem 2.1 (Boundary turning number) is equivalent to the Poincaré Hopf theorem with a proper definition for the index of a singularity, which will be developed in next subsection.

It links the turning numbers of homological cycles as γ 1 ≡ l γ 2 ⇔ ∃S such that γ 1 -γ2 = ∂ S Fig. 6. The turning number associated to a generator defines topology that cannot be captured by singularity indices. The difference of topology of theses direction fields without any singular points is defined by the turning numbers of the generator (black cycle) which are respectively -1, 0 and 1 THEOREM 2.2 TOPOLOGICAL EQUIVALENCE. Two direction fields defined over a surface S are homotopic iff they have the same turning numbers along the cycles of a homology basis H(S ) of S :

d 1 ≡ t d 2 ⇔ ∀γ ∈ H(S ), T d 1 (γ) = T d 2 (γ)
Theorem 2.2 shows that direction fields have homotopy classes isomorphic to Z 2g+b-1 . This comes from the fact that the 2g generators of the surface and b -1 surface borders define a homology basis. Hence a direction field has 2g + b -1 = 1χ(S ) topological degrees of freedom (TDoF) on a genus g surface with b borders. In the case χ(S ) = 2 (S is a topological sphere), this number is negative because it is impossible to build a continuous direction field without introducing at least a singularity (a hole), which brings χ(S ) to 1.

Both theorems are proved in Appendix A. With these two theorems, one can exhibit, understand and control all the TDoF of a direction field.

We are now going to explain how turning numbers generalize the usual notion of singularity index.

Singularities

Let v : R 2 → R 2 be a vector field. It is usually assumed that the zero set of v:

{x ∈ R 2 | v(x) =
0} consists of a finite number of distinct points P i which are called the singularities of v.

The singularities can then be classified by their index:

I d (P i ) = 1 2π ∂ Ω(P i ) dθ (6)
where Ω(P i ) is a small neighborhood of P containing no other singularities (zeros) of the vector field, and θ is the angle formed by the vector field and a reference vector. The singularity index is necessarily an integer that equals 1 for sources, vortices and sinks, -1 for saddles. An index of 0 corresponds to a degenerate singularity, which means the corresponding singularity can be removed without creating another singularity.

• 11

The notion of index applies naturally to direction fields and keeps the same definition (6) as it depends only on the direction of the vector. As singularities are zeros, the direction cannot be defined there, hence we identify singularities of a direction field with "holes" in its domain of definition, For N-symmetry vector fields, indices can still be defined by ( 6), but they are now a multiple of 1/N because θ is defined modulo 2π/N. This 2D definition for direction field singularity index cannot be directly extended to surfaces embedded in R 3 because it lacks a reference vector to define θ . However, it can be shown that in R 2 :

T d (∂ Ω(P)) = I d (P) -1 (7)
As the turning number can be extended to surfaces in R 3 because the reference vector is the tangent to the cycle, this allows us to extend the definition of the index of a singularity on a surface in R 3 to:

I d (P) = 1 + T d (∂ Ω(P)) (8)
With this definition for the index, we can generalize the Poincaré-Hopf theorem to Nsymmetry direction fields: Proof: Notice first that this theorem holds for surfaces without borders, so the number of borders is absent in the Euler characteristic. Let us call P i the point of index I i and S h = S \{Ω(P i )} (surface with borders in place of singularities), on which the field is continuous because it does not contain the singularities. Applying the boundary relation (5) to S h yields:

T (∂ S h ) = T (∂ S ) - b ∑ i=1 T (∂ Ω(P i )) = - b ∑ i=1 (I i -1) = -χ(S h ) = -χ(S ) + b
we have T (∂ S ) = 0 because S has no borders, and we get the Poincaré-Hopf theorem by subtracting b in these equalities 2

The idea of replacing singularities with holes allows singularities to be handled as borders. Moreover borders are topological objects that are very well understood through cycle homology. In what follows, we will use again this idea of replacing singularities with holes. We will then characterize the singularity by the behavior of the field along the border of the hole.

Notice that singularities do not capture all the topological degrees of freedom, since a homology basis also contains generators of the surface, that do not enclose any singularity (see Figure 6). This justifies why we need the concept of turning number to capture the direction field topology, and not only the singularity indices. We will now explain how the turning numbers of a direction field in the data structure from [START_REF] Li | Representing higher-order singularities in vector fields on piecewise linear surfaces[END_REF]] depend only on the period jumps. As we have proved in this section that the topology of a direction field is completely defined by its turning number, this means we will have a total control over the direction field topology through period jumps.

DISCRETE DIRECTION FIELD

We will now briefly recall the period jump based discrete representation for N-symmetry direction fields presented in [START_REF] Li | Representing higher-order singularities in vector fields on piecewise linear surfaces[END_REF], and explain its coherence with the continuous theory presented above.

Period jumps based discretization

The first step of the discretization of a direction field d such as described in [START_REF] Li | Representing higher-order singularities in vector fields on piecewise linear surfaces[END_REF]] is to sample it at the center of each facet. This is done by choosing a reference direction d 0 in each triangle (for instance one of the oriented edges of the triangle), and defining d by the angle θ if forms with d 0 .

The discrete curvature of d can be defined as the integral of the continuous curvature along each dual edge e * , which is given by applying (3) to e * :

κ d = κ 0 + θ ( f end ) -θ ( f start ) + 2π p N (9)
where the curvature of the base field κ 0 is defined as the oriented angle in (-π, π] between d 0 ( f start ) and d 0 ( f end ), and p is an integer called period jump that corresponds to the ambiguity in the interpolation from d 0 ( f start ) to d 0 ( f start ) (see Figure 7). Once p is chosen, the ambiguity is solved, and the direction can be interpolated on the whole mesh (except at vertices) by the interpolation proposed in [START_REF] Li | Representing higher-order singularities in vector fields on piecewise linear surfaces[END_REF]]. Hence we can define any direction field d through:

(1) The choice of a reference direction d 0 on each facet (i.e. an oriented edge of the facet),

(2) A scalar θ that represents the angle between d and d 0 on each facet, such that d = R( d 0 , θ ) on dual vertices.

• 13

(3) An integer p that represents the behavior of d between discretization points, such that the discrete curvature of d is

κ d = κ 0 + dθ + 2π p N
Using the definition of the turning number (4) and the discrete curvature (9), we derive an expression for the turning number of a cycle γ in the discrete setting:

T d (γ) = 1 2π ∑ e * ∈γ κ d (e * ) -κ γ (e * ) = 1 2π ∑ e * ∈γ κ 0 (e * ) -κ γ (e * ) + 1 N ∑ e * ∈γ p(e * ) (10)
We used the fact that the terms in θ ( f end )θ ( f start ) all cancel out because each θ once with each sign in the sum since γ is a cycle. As the first term is independent of the geometry θ of the field, this means that the turning numbers of d are defined entirely by the period jumps p. Finally, ∑ γ κ γ is computed as the sum of oriented angles between successive edges along all triangles crossed by γ. To compute the singularity index at a vertex v, we can simply choose γ = ∂ v * (v * is a dual facet), in which case ∑ γ κ γ is simply the sum of angles incident to v, and we get the expression for the index:

I d (v * ) = 1 + 1 2π ∑ e * ∈v * (κ 0 -κ γ ) + ∑ e * ∈v * p(e * ) N = I 0 (v) + ∑ e * ∈v * p(e * ) N (11)
We interpret I 0 (v) as the index of the base field, that is computed by writing it:

2πI 0 (v) = ∑ e * ∈v * κ 0 (e * ) -κ γ (e * ) + 2π = ∑ e * ∈v * κ 0 (e * ) + A d (v) (12)
where A d (v) is the angle defect at v. This discretization will now allow us to derive a simple algorithm to compute the period jumps p such that the topology of the field matches user constraints.

ZIPPING ALGORITHM: CONSTRAINING SINGULARITY INDICES

When designing a direction field, we can choose the period jumps p(e * ) on each edge. This completely defines the topology of the represented direction field. We propose in this section an algorithm that computes the period jumps such that the singularities and their indices are exactly the ones defined by a user. However, the period jumps control all TDoF, including those that do not correspond to singularities, such as turning numbers of generators (see Figure 6). We choose to let these TDoF free, and let our optimization play on these additional TDoF to find a smoother solution. As a result the algorithm will not explicitly compute all period jumps, but find some free period jumps and express all the other period jumps as a combination of the free ones.

Problem setting

The problem we tackle in this section is the following:

Given constrained singularity indices I c (v i ) at each interior vertex v i of the mesh, find the integers p such that:

I d (v i ) = I c (v i ) ∀v i (13)
• Using (11), the above expression can be rewritten as:

∑ e * ∈v * p(e * ) = N(I c (v) -I 0 (v)) = ∆I(v) (14)
in other words, we are trying to find integers p on oriented dual edges that sum up to ∆I(v) around each vertex v. We will now propose a greedy algorithm to compute a solution to (14) in linear time, which works by classifying edges into three sets.

Algorithm

As in the continuous setting, the homology basis H(M) of M consists of 2g generators and b -1 borders (we will discuss later the case b = 0). As the turning numbers of the cycles in H(M) do not correspond to singularities, we will leave the turning numbers of these cycles free. As the direction cannot be defined at vertices, the direction field is defined over M h , which is the mesh M with holes at every interior vertices v i allowing a singularity to occur at v i . Each removed interior vertex v i adds a cycle ∂ v * i around it to the homology basis, except one of them that can be chosen arbitrarily (we call it v last ). Hence M h has a homology basis

H(M h ) = H(M) ∪ v i =v last {∂ v * i }.
The turning numbers along these additional cycles in H(M h ) correspond to singularities located on interior vertices, so we want to control them explicitly. Hence, the basic idea in Zipping is to progressively classify the dual edges in three sets:

(1) The set E * 0 of null edges, whose associated period jumps we can set to 0 without constraining any turning number. Hence we have p(e * 0 ) = 0 (2) The set E * f ree of free edges, whose associated period jumps correspond to the turning number of a cycle homological to a cycle in H(M), hence corresponding to a TDoF but not to a singularity (see Figure 9). We choose the name free because these TDoF will be left free in the subsequent optimization. item The set E * dep of dependent edges (they may depend on free edges), whose associated period jumps constrain the turning numbers of cycles ∂ v * i around a single vertex v i , hence constraining the index of a the singularity at v i .

Zipping is a greedy algorithm that computes this classification iteratively. Zipping also computes an expression of the period jumps of edges in E * dep as a function of the free period jumps of edges in E * f ree , such that the indices satisfy their constrained value. ( 1) Fill E * 0 : While it is possible, add to E * 0 edges e * 0 that do not close any cycle, and set their period jumps to 0: p 0 (e * 0 ) = 0, c(e * 0 ) = 0 This is equivalent to building a spanning tree of the dual graph G * .

(2) Fill E * dep : While it is possible, add to E * dep edges that close dual cells (cycles around a single primal vertex) and compute the corresponding period jumps such that they satisfy ( 14).

(3) Add an edge to E * f ree : any remaining edge necessarily closes a cycle (else it would have been added to E * 0 ) that does not enclose a single vertex (else it would have been added to E * dep ). Hence the corresponding cycle is homological to a border or to a generator of the mesh, which turning numbers we want to keep free, so any such edge

• 15
Algorithm 1 Zipping algorithm (see Figure 8)

Build a recovering tree E * 0 of G * // grow black edges ∀e * 0 ∈ E * 0 set p 0 (e * 0 ) ← 0, c(e * 0 ) ← 0 E * do f ← E * \E * 0 V zip ← V i f ree ← 0 // Number of free variables found while E * do f = / 0 do while V zip = / 0 do Take v z ∈ V zip and remove it from V zip if v z / ∈ ∂ M and ∃! unset edge e * z ∈ ∂ v * z then // zip red (then blue) edges Move e * z from E do f to E dep Set p 0 (e * z ) ← ∆I(v z ) -∑ e * ∈∂ v * z \e * z p 0 (e * ), Set c(e * z ) ← -∑ e * ∈∂ v * z \e * z c(e * ), Add the face opposite to v * z across e * z to V zip end if end while // free green edge i f ree ← i f ree + 1 Take e * f ree ∈ E * do f and move it to E * f ree
Set p 0 (e * f ree ) ← 0, c i (e * f ree ) ← δ i,i f ree Add the 2 faces adjacent to e * f ree to V zip end while can be added to E * f ree . These edges have trivially: p 0 (e * f ree ) = 0. The vector c(e * f ree ) only contains zeros except for a 1 that corresponds to e * f ree . (4) After freeing an edge, it becomes possible again to find edges that close dual cells, such that steps 2 and 3 may be iterated until no edge remain. The period jumps computed in step 2 will depend on the period jumps of the edges that have been freed in a step 3 (see Figure 9).

The complete Zipping algorithm is detailed in Algorithm 1. To handle period jumps dependencies, we store all period jumps under the form:

p(e * ) = p 0 (e * ) + c(e * )p f ree where p 0 is an integer, and c a vector of integers which size equals the number of edges already freed.

Topological issues

If M has no border, we have:

∑ v∈V ∆I(v) = ∑ v∈V ∑ e * ∈v *
p(e * ) = 0 as each edge is taken once in each direction. The base field singularities naturally satisfy the Poincaré-Hopf theorem 2.3 (∑ V I d 0 (v) = 2 -2g). Therefore, the index at the last traversed vertex v last satisfies:

∆I(v last ) = -∑ V \v last ∆I(v) , I c (v last ) = 2 -2g -∑ V \v last I c (v) (15)
Hence the index of v last adapts to ensure the Poincaré-Hopf theorem is satisfied. If the constraints on indices have been chosen such that:

∑ V I c (v) = 2 -2g
then the index I d (v last ) will be 0, which proves that no new singularity will appear. As the position of v * last depends on some choices arbitrarily made by the greedy algorithm, it cannot be easily determined, hence it is highly recommended to run the algorithm with constrained indices that sum up to 2 -2g to avoid the appearance of a random (but necessary) singularity.

If M has borders, they are also handled by the Zipping algorithm, but their turning numbers are not constrained. In fact, we can also leave some vertices unconstrained by simply declaring them as border vertices before Zipping. However, in this case, we cannot guarantee that no undesired singularity appears, but only authorize where they can appear. A simple way to ensure that no singularity appear when M has borders, is to triangulate the holes of M before running our Zipping algorithm. The index of a border on such a field will be equal to the sum of the indices of the border vertices. In this case again, no singularity will appear iff the sum of indices of constrained vertices (including border vertices) sum up to 2 -2g.

• 17

APPLICATION TO DISCRETE DIRECTION FIELD DESIGN

We now show a simple use of our algorithm to design smooth direction fields with strong topological constraints.

Problem formulation

Using the notions above, the problem we want to solve is:

Given: -a mesh M =< V , E , F > -a constrained singularity index I c (v)/N given on each vertex v ∈ V . -a constrained direction d c ( f ) given on each facet of a subset F c ⊂ F .
Interpolate the constraints as smoothly as possible. More formally, minimize: 16) subject to the constraints:

E G * (θ , p) = ||κ d || 2 = ∑ e * ∈E * κ d (e * ) 2 = ∑ e * ∈E * θ ( f end ) -θ ( f start ) + κ 0 (e * ) + 2π p(e * ) N 2 (
I d (v) = I d 0 (v) + ∑ e * ∈∂ v * (e * ) = I c (v) ∀v ∈ V (17) d( f ) = d c ( f ) ∀ f ∈ F c (18)
As explained in the previous section, we build a representation for the direction field that implicitly enforces the topological constraints, by expressing all the period jumps p(e * ) as a function of a limited number of free period jumps. We now describe a design algorithm that uses this structure to reduce the problem to a simple quadratic minimization procedure.

Algorithm

The design algorithm links between user inputs and Zipping inputs by computing the indices of the base field and making the difference with user constrained indices, and uses the output of the Zipping algorithm to ensure that the field topology is constrained during the creation of the direction field. The minimization problem is then solved in two passes:

(1) Minimize E G * with respect to θ and p by assuming that p is continuous

(2) Minimize with respect to θ only, with the p being set to their rounded value of the first pass

We use for both passes a standard formula [START_REF] Levy | Numerical methods for digital geometry processing[END_REF]] to solve the problem of minimizing

([A f , A l ][x f , x l ] t -B)
2 where x f are variables and x l are set:

x f = (A t f A f ) -1 A t f (B -A l x l ) (19)
This method is not guaranteed to find the global minimum with respect to the discrete variables p(e * ), but offers good results in practice that satisfy all the constraints. For the

[θ 1 f , k 1 f ree ] t = ([A f ,C] t [A f ,C]) -1 [A f ,C] t (B -A l θ l ) (10) Pass 2: [θ 2 f ] t = (A t f A f ) -1 A t f (B -[C, A l ][rnd(k 1 f ree ), θ l ] t )
where rnd is the rounding to the nearest integer value. (11) Apply rotations θ 2 f to d 0 to get a direction d on each facet of M.

continuous variables θ ( f ), if at least one directional constraint is set, A f is of maximal rank, therefore, by Gramm's theorem, A t f A f is non-degenerate and our algorithm finds the unique minimum. We also noticed that we could improve the visual aspect of the direction field near the constrained directions by adding a Laplacian smoothing term to the energy (16).

RESULTS

Our framework allows to create and edit direction fields on a surface via topological and geometric constraints. As illustrated in Figure 10, only few constraints are required to create the desired field. The same picture also demonstrates that our method is not affected by the complex geometry of the hairs of the David, whereas the previous methods only based on relaxation ( [START_REF] Ray | Periodic global parameterization[END_REF]] or [START_REF] Hertzmann | Illustrating smooth surfaces[END_REF]) are trapped by this geometry and generate many singularities in this zone.

We have also tested this algorithm on large models from the Stanford 3D Scanning Repository called the statue and Lucy. These two models are treated in less than a minute using TAUCS for the quadratic minimization, and smooth fields are obtained (see Figure 11). Note that these two models have a complex topology (large genus g and number of borders b). The turning number around each border is used to counter the effect of each handle (a small handle is equivalent to a singular point of index 2).

Our algorithm also deals nicely with important constraints that can be applied on geometry (see the rotation constraints in Figure 12) and on topology (see Figure 6).

Finally, we compare with recent work in direction field design. Figure 13 shows that starting from a user defined set of constraints (singularities and directional constraints), our method generates a direction field with no other singularity (Right). To our knowledge, this is the first algorithm that achieves direction field design with provably correct singularity control. Figure 14 shows a more difficult configuration, successfully handled by our algorithm.

CONCLUSION

We have introduced a generalization of direction field that is general enough to manipulate direction fields with different kind of symmetries, which applications range from texture synthesis to quad remeshing.

As the representation we use separates the topology and the geometry of the direction fields, the singularities can be constrained through the topology by using a simple greedy algorithms, and the geometry can be smoothed by simply minimizing a quadratic form. This has been made possible by a complete analysis of the direction field topology including the characterization of the field behavior along cycles, and an extension of the Poincaré-Hopf theorem to rational indices. A new direction field smoothing algorithm has been presented to demonstrate the benefits of our structure. With our algorithm, one can for the first time provably control the placement of singularities. We are confident this new • Fig. 11. Large models with many borders and high genus efficiently processed by our framework.

Fig. 12. Directional constraints are applied to a direction field. Notice that the direction is given by a rotation of the base field that can be greater than 2π.

way to consider and manipulate direction fields will set a clean and practical basis for a broad range of applications in geometry processing.

• 21

A. TURNING NUMBERS FUNDAMENTAL PROPERTIES

We give here an outline of the proof of two fundamental turning number properties.

A.1 Boundary property

The boundary property states that for any direction field on a 2-manifold S, we have T (∂ S) = -χ(S) which generalizes the Poincaré Hopf theorem to N-symmetry direction fields. We first establish three simple results, then prove the equation by structural induction.

LEMMA A.1. The reversal of a cycle has opposite turning number:

T (-γ) = -T (γ) (20)
Proof: integrating in opposite direction changes ds in -ds 2 LEMMA A.2. If we call A and B two sub-manifolds of S h (see Figure 15), the turning numbers of ∂ A and ∂ B are linked by the equation: ). Because we are on the unit disk, the arclength s is equivalent to the angle on the disk boundary, so we have

T (∂ A) + T (∂ B) = T (∂ (A ∪ B)) + T (∂ (A ∩ B)) (21
• t γ = (cos(s), sin(s)) t and n γ = (sin(s), cos(s)) t , which gives:

κ γ = ∂ t γ ∂ s • n γ = 1
Let x = (x = cos(s), y = sin(s)) t be the position vector on the disk boundary. We have:

κ d ( t γ ) = ∂ d ∂ s • d ⊥ = ∂ d ∂ x • d ⊥ , ∂ d ∂ y • d ⊥ t . t γ = κ d . t γ
one can verify that the rotational ∇ ∧ κ d = 0, so κ d is a gradient, or in other words there exist a scalar field θ such that κ d = ∇θ . Thus its integral over a closed loop is zero, so the turning number evaluates to:

T d (∂ D 2 ) = 1 2π ∂ D 2 (κ d -κ γ )ds = 1 2π ∂ D 2 ∇θ . t γ ds - 1 2π 2π s=0 1 • ds = -12
These lemmas allow us to give a proof of Theorem 2.1:

T d (∂ S) = -χ(S) (22)
Proof: It is well known in topology that any orientable surface (with borders) can be cut along g cycles to obtain a sphere (with borders). Hence we start by proving ( 22) for a sphere with b > 0 borders.

If S is a topological disk (b = 1, χ = 1), its boundary is contractible so its turning number is -1 by Lemma A.3 which proves ( 22) for b = 1.

If S is a topological cylinder (b = 2, χ = 0), its boundary is composed of 2 borders: ∂ S = γ 1 + γ 2 . γ 1 and -γ 2 are homotopic so they have the same turning numbers, so by (20) we have:

T (∂ S) = T (γ 1 ) + T (γ 2 ) = T (γ 1 ) -T (-γ 2 ) = 0 which proves (22) for b = 2.

For higher numbers of borders, we prove the property by induction: assume ( 22) is true ∀b ≤ B, and let S have B + 1 borders γ 1 ...γ B+1 . Let A B and A B+1 be two connected submanifolds which boundary contains γ B and γ B+1 (respectively) and intersect in a topological disk. S\{A B ∪ A B+1 } has at most B borders, so it satisfies (22):

T (∂ (S\{A B ∪ A B+1 })) = B-1 ∑ i=1 T (γ i ) + T (∂ {A B ∪ A B+1 }) = B-1 ∑ i=1 T (γ i ) + T (∂ A B ) + T (∂ A B+1 ) -T (∂ {A B ∩ A B+1 }) =
T (∂ S cut ) = n ∑ i=1 T (γ i ) + g ∑ i=1 T (γ G i ) + T (-γ G i ) = 2 -2g -b
as the second sum is null by (20). This finally proves Theorem 2.1 in the general case 2

A.2 Topological equivalence

We will need a simple lemma to prove the topological equivalence:

LEMMA A.4.

T d 1 (γ) = T d 2 (γ) ∀ cycle γ ⊂ S ⇔ T d 1 (γ) = T d 2 (γ) ∀γ ∈ H(S )
where H(S ) is a homology basis on S . In other words, the turning numbers of all cycles on S depend only on the turning numbers of a homology basis.

Proof: By definition of homology basis, any cycle γ is homological to a cycle ∑ i a i γ B i where γ B i ∈ H(S ) are basis cycles, such that ∑ i a i γ B iγ is a boundary. Hence, by Theorem 2.1, we have:

T d 0 (γ) = ∑ i a i T d 0 (γ B i ) + χ(S ) = ∑ i a i T d 1 (γ B i ) + χ(S ) = T d 1 (γ)
because the turning numbers are equal along the cycles of the homology basis. 2

We can now prove the Theorem 2.2 which states that two direction fields are homotopic if and only if they have the same turning numbers along the cycles of a homology basis of S .

Proof: The direct implication is immediate because turning numbers are preserved by homotopy. For the converse, let d 0 and d 1 be the two fields with the same turning numbers along the cycles of a homotopy basis. By Lemma A.4, d 0 and d 1 have the same turning numbers along any cycle γ:

T d 0 (γ) -T d 1 (γ) = γ (κ γ -κ d 0 ( t γ ))ds - γ (κ γ -κ d 1 ( t γ ))ds = γ (κ d 1 -κ d 0 )( t γ )ds = 0 (κ d 1 -κ d 0 )( t γ
) is a linear application which circulation along any closed curve is zero. This means that there exist a continuous scalar field θ that satisfies:

(κ d 1 -κ d 0 )( t γ ) = ∇θ .

Fig. 1 .

 1 Fig. 1. Examples of singularities with fractionary index. Left: Singularity of index 1/4 of a natural smooth cross field on a cube like object. Right: Singularities of index 1/2 and -1/2 on a fingerprint.

Fig. 2 .

 2 Fig.2. Homotopy for direction fields: the middle field is homotopic to the one on the right but not to the one on the left (they can not be transformed continuously one into another)

•Fig. 3 .

 3 Fig. 3. Left: Cycles on a surface with borders (black areas are holes). Only the cycles shown on Right are exact because they are the boundaries (with conormal pointing outwards) of submanifolds of S (in light color).

Fig. 4 .

 4 Fig.4. A Darboux frame on a chain (γ, green) consists of the tangent ( t γ , blue), conormal ( n γ , green) and normal ( n, red). When moving along an infinitesimal portion ds of γ, the direction turns of an infinitesimal angle dθ in the Darboux Frame.

  THEOREM 2.3 POINCARÉ-HOPF. The sum of singularity indices on a closed surface S equals its Euler characteristic: b ∑ i=1 I i = χ(S ) = 2 -2g

•

  Fig. 7. p solves the ambiguity for equal angles with the reference direction d 0 at endpoints of a dual edge e * .

•Fig. 8 .

 8 Fig. 8. Zipping: 1-Grow black edges (width first search) 2-Zip red edges 3-Free green edge 4-Zip blue edges (blue edges depend on the freed green edge)

• 2

 2 Algorithm Design algorithm (1) Set I c (v) = 0 on V \V c // no singularities outside V c (2) Choose a direction d 0 ( f ) on each facet // base field (3) Compute κ 0 (e * ) as the angles between d 0 ( f end ) and d 0 ( f start ) for each dual edge. (4) Compute the angle defects A d (v) at all vertices. (5) Compute the base field indices I d 0 (v) using (12) (6) Compute ∆I(v) = N(I c (v) -I d 0 (v)) (7) Apply the Zipping Algorithm (8) Build the linear system [A f , A l ,C][θ f , θ l , k f ree ] t = B corresponding to (16). Each line of the system corresponds to an edge e * : [A f , A l ] contains +1 and -1 at the indices corresponding to the θ at the 2 extremities of e * , C contains c(e * ) corresponding to k f ree . B contains the κ d 0 (e * )k 0 (e * ). (9) Pass 1:

Fig. 13 .

 13 Fig.13. Left: result of[Fischer et.al 2007] with a single singularity constraint (green dot); Center: result of[Fischer et.al 2007] with the same singularity constraint and directional constraint (green arrow), additional singularities appear (red dots); Right: our result with the same set of constraints. As predicted by our theorem, no additional singularity appears.

Fig. 14 .

 14 Fig.14. A particularly challenging data set for direction field design algorithms. Our algorithm successfully satisfies the four user-defined singularities (green dots) and does not generate any additional singularity.

)

  Proof: The term in κ d is preserved because integrated on the same set∂ A ∪ ∂ B = ∂ (A ∪ B) ∪ ∂ (A ∩ B).The term in κ γ , is preserved by application of the Gauss-Bonnet formula to the equality A + B = A∪B + A∩B 2 LEMMA A.3. Boundaries of topological disks have a turning number of -1. Proof: Let γ C be the boundary of a topological disk D. By definition, there exists a continuous bijective mapping m : D → D 2 where D 2 is the unit disk in R 2 . Using this application, we can define the image field d = m( d) of the direction field on D 2 . The turning number is preserved by m as it is continuous, so T d (∂ D) = T d (∂ D 2

  i ) + 1 = χ(S\{A B ∪ A B+1 }) = -χ(S) + 1 which proves (22) recursively. We can now get to the most general case where S has genus g and b borders γ 1 ...γ b . By definition of the genus, there exists a family γ G 1 , ...γ G g of single cycles such that S\{γ G 1 , ...γ G g } • 25 is connected. This operation does not change the Euler characteristic of the surface (b ← b + 2g, g ← 0 preserves χ = 2 -2gb). The result is a sphere S cut with b + 2g borders γ 1 , ...γ b , γ G 1 , -γ G 1 , ...γ G g , -γ G g , so we can apply (22):

  t γ As θ is defined up to a constant, we can chose it to match the angle (in [0, 2π[ for instance) between d 0 and d 1 . Thus we have d 1 = R( d 0 , θ ) where R is the rotation around the surface normal. As θ is continuous, we can build a continuous function Γ : [0, 1] → D N (S ) defined by Γ(t) = R( d 0 ,tθ ), such that Γ(0) = d 0 and Γ(1) = d 1 which proves that d 0 and d 1 are homotopic. 2

Table I .

 I Timings obtained on a Pentium IV 1.7Ghz

		Lucy	The statue
	genus	45	13
	#borders	47	9
	#triangles	125000	300000
	time	12s	35s

  Fig. 10. Direction field edited by topological and geometric constraints. Red (resp. blue) dots are constrained singularities with positive (resp. negative) index. Orange arrows are directional constraints.

	• 19

ACM Journal Name, Vol. V, No. N, Month 20YY, Pages 1-25.

ACM Journal Name, Vol. V, No. N, Month 20YY.

Acknowledgments

We want to thank the Associate Editor and the anonymous reviewers for their detailed reviews and numerous comments that helped us improve the quality of this paper. We also thank the Stanford Scanning Repository and Digital Michelangelo Project for the 3D models. This work was supported by Microsoft Research Europe (Geometric Intelligence Grant)