Nicolas Bremard
email: nicolas.bremard@lifl.fr

Laurent Grisoni
email: laurent.grisoni@lifl.fr

Bruno De Araujo
email: bruno.araujo@inria.fr

Interaction events in contactless gestural systems: from motion to interaction

Keywords: H.5.2 [Information Interfaces and Presentation]: User interfaces-Interaction styles, Prototyping parameterized gesture, contactless interaction, numerical art

Contactless interactive systems appear to have an interesting context of application in interactive art, however practical example often involve low interactivity; we think this is due to the fact that this interaction paradigm is still poorly understood, and tools for interaction design are still missing. We present a simple interaction technique (Point-and-Move) that allows to point-and-click without instrument through distance; through this example, we show that the presented concept and framework has nice and flexible properties. Finally, we illustrate the approach on the two artistic installation we participated to, and also provide uncontrolled user experience feedback on the presented interaction technique, from these installations.

INTRODUCTION

An important challenge about contactless gestural interaction is to understand how to take the best from user gesture expressivity: WIMP paradigm is based on as-elementary-aspossible interaction input (control of a pointer position on a screen, and event generation through clicking). By comparison, user hand movements has a much larger expressivity, from which interactive applications only take poor benefit. Such richness can be, on some specific application contexts, a difficulty (e.g. within the multi-touch context, it is known that users do not matter about how many fingers actually touch the surface [START_REF] Wobbrock | User-defined gestures for surface computing[END_REF]; understanding such variability is a key point for multi-touch systems [START_REF] Rekik | Towards many gestures to one command: A user study for tabletops[END_REF]). In the context of gestural interaction, the problem is still very large, since Permission to make digital or hard copies of part or all of this work for personal or classroom use is granted without fee provided that copies are not made or distributed for profit or commercial advantage, and that copies bear this notice and the full citation on the first page. Copyrights for third-party components of this work must be honored. For all other uses, contact the owner/author(s). Copyright is held by the author/owner(s). MOCO'14, June 16-17 2014, Paris, France. ACM 978-1-4503-2814-2/14/06 ...$15.00. http://dx.doi.org/10.1145/2617995.2618027. cases of use are not clearly identified yet, and many open question remains unsolved, the most important ones probably being why interacting with a computer using hands, and how to achieve that. Indeed, hand movements are classically the mean to act on physical world, complement inter-human communication, mimic something . It is a recent idea to think about user hand movements as a mean to interact with an application. The idea is all the more complex as structured tasks are foreseen.

One of the first question that is to be studied is how to handle the non-separable nature of user movements: grasping a glass, or pointing at something, involves one single gesture, and the intrinsic nature of such movement is very much different from that of grasping a mouse and moving it for interaction:from a computer science point of view, interaction gesture would often have to be interpreted in terms of action (what is being done?) and parameters (on which part of the application? In which manner?). In this article we address the question of how to generate interaction events from user gestures. We aim at providing elements that are independent of any standard recognition technique; instead, our goal is to provide a framework for fast prototyping of explicit gestural interaction recognition and design.

Few tools are currently existing for fast prototyping of contactless and full-body interactions. The context of interactive art is an interesting application context to such systems, as several recent movement low-cost acquisition systems open the way to accessible installation (however these low-cost system are often of low space and time precision). We introduce an approach based on simple, basic numerical tools, that can be easily implemented and combined by interaction designer, so that any target system can be prototyped. Our approach is flexible, can be easily adapted to specific input, and does not involve any sophisticated recognition method; as a result it can be easily reproduced and adapted to specific interaction design problems. We illustrate our method within two artistic installations where we propose an interaction technique combining pointing and selection within gestural, distant and non-instrumented, interaction context using low-cost non-intrusive tracking devices such as the Kinect depth sensor or RGB cameras. Finally, we discuss user feedback from two artistic installations.

RELATED WORK

Researches have already been conducted to improve standard gestures recognition processes within contactless interactions. Wilson et al. ([15], [START_REF] Wilson | Recognition and interpretation of parametric gesture[END_REF]) proposed to improve Hiden Markov Model (HMM) by extracting meaningful parameters from the move and modifying HMM's output with it. He demonstrated that meaningful parameters could be computed from moves and be used outside a gesture, now called a parameterized gesture.

Karam and Schraefel [START_REF] Karam | A taxonomy of gestures in human computer interactions[END_REF] proposed a classification of gestures based on four key elements: gesture styles, application domains, input and output technologies. They determined six categories : deictic, manipulative, semaphoric, gesticulation, sign language and multiple. It have been reviewed by Aignier et al. [START_REF] Aigner | Understanding mid-air hand gestures: A study of human preferences in usage of gesture types for hci[END_REF] making the distinction between static and dynamic gestures, admitting that dynamic is a key for gestures classification. While these classifications are good to describe gesture in an abstract way, they do not really describe how gestures can be decomposed into parameters to be recognized and applied on a given interaction scenario.

Several works have focused on how to define gesture within the user interaction space and discuss the idea to interact within a space that user would imagine in front of him/her. Rateau et al. [START_REF] Rateau | Mimetic interaction spaces : Controlling distant displays in pervasive environments[END_REF] create this surface from a gesture, whilst Chan et al. [START_REF] Chan | Touching the void: Direct-touch interaction for intangible displays[END_REF] and Gustafson et al. [START_REF] Gustafson | Imaginary interfaces: Spatial interaction with empty hands and without visual feedback[END_REF] hang it to the nondominant hand. Even if users can create as many surfaces as they need, these interactions are still in two dimension and are restraining the application with a cursor.

A special focus have been done regarding pointing gestures to a surface for emulating 2D cursor. Pierce et al. [START_REF] Pierce | Image plane interaction techniques in 3d immersive environments[END_REF] use it to manipulate virtual objects. Other studies try to recognize a command, such as a click : with non dominant hand [START_REF] Schwaller | Two handed mid-air gestural hci: Point + command[END_REF] or with fingers [START_REF] Vogel | Distant freehand pointing and clicking on very large, high resolution displays[END_REF]. These works create a two dimensional interaction too, and are not enough for creating a high level of interactivity. Ren et al. [START_REF] Ren | 3d selection with freehand gesture[END_REF] propose a three dimensional selection technique in which the hand's direction is used to emulate a command (choose an option). This technique constantly check if the hand is moving and trigger the associated command, which avoid any other gesture to be associated with this hand. Choumane et al. [START_REF] Choumane | Buttonless clicking: Intuitive select and pick-release through gesture analysis[END_REF] suggest to use the dynamic, represented by acceleration, of the gesture to select an object. The technique based on acceleration let the hand free to do other gestures. That proves the use of gesture dynamic has more advantages than direction alone for gestures recognition.

Regarding tools for easing interaction design, some graphical tools have been proposed such as ICON [START_REF] Dragicevic | Input device selection and interaction configuration with icon[END_REF] and Flowstate [START_REF] Appert | FlowStates: prototypage d'applications interactives avec des flots de données et des machines à états[END_REF]. These tools are efficient when used with keyboards and devices with two dimensions such as mice and touch screen; but they are unusable in three dimensions : they have not been designed for it. Our proposal presents a solution to this problem easing the design of gestural based interactions using 3D data gathered from the user.

GESTURAL INTERACTION DESIGN FRAMEWORK

We present here a new tool for interaction design, fully compatible with contactless interactions and parameterized gestures. A parameterized gesture is a gesture from which meaningful values related to the gesture can be extracted; such as a pointed position for a pointing gesture or a direction for a manipulative action. Some of these attributes can be extracted from the gesture itself and gathered information from the tracking device without requiring any information. For example, it can be a given position, a speed, an acceleration, a direction or a duration associated to the ges-ture. Others attributes are related to the environment and the application context which the gesture is subject to. For example, which object (virtual or real) or person someone is pointing at while doing a pointing gesture.

General overview

Our tool is based on three components : an input device abstraction, a user model and set of observers for gesture recognition. We follow a plugin architecture to represent the input devices and forward the input data to the other components using an update method. First we update our user model converting the data gathered from the input device for example the skeleton data of tracked user when using the Kinect device. Our user model represents a full-body user by a set of joints. Each joint has three dimensional position and orientation. The set of joints can be customized , i.e. they can be ignored if the data is missing or not needed by the designer, or added if they are available by a given input device. This information is then forwarded to the observers designed to recognize a set of gestures attributes.

Observers

The Observers are basic methods, called at each input update, that evaluate an observer value. These values are tested within each interaction automata associated to interaction gestures. The following set of observers are modeled by our approach and can be combined depending on the kind of gestures they are designed to recognize.

zero-Order Observer : it detects when a given joint is close to a given state. For example it can represent the positional or angular distance to a recorded position or orientation using use absolute and relative coordinates.

First-Order Observer : these observers recognize if a joint is moving in a specific direction, according to the coordinate system, or if the speed or the angular speed of its orientation are above,below or between specified values.

Second-Order Observer : it recognizes if a joint is moving at, above, below or between specified accelerations.

Time-Switch Observer : these observers are designed to detect predefined changes about a joint. Changes can be in position, speed, acceleration, or orientations for example.

World-Aware Observers: it integrates additional informations related with the environment or the type of interaction.

In most examples, they represent either physical objects or screen-displayed widgets related to a given gesture.

Personalized Observer : it allows to create your own gesture recognizer integrating it seamlessly in our architecture.

Composite Observer : it enables to combine several observers and create a small automata with them.

Gesture example: Point and Select

Providing generic interaction techniques, such as pointand-click, in contactless interaction, is a question to which no current interaction paradigms answer. Indeed, in public context, instrumented interaction is rarely available, as a result no physical button can achieve the click feature. Moreover, most current contactless acquisition devices are still with low resolution and frequency, making data noisy, and body reconstruction usually coarse.

In the practical artistic request described in this article, we had the two different options to explore, resulting in two way to select with different accuracies due to contrast in the gesture dynamic. These pointing techniques, proposed by [START_REF] Jota | A comparison of ray pointing techniques for very large displays[END_REF][START_REF] Choumane | Buttonless clicking: Intuitive select and pick-release through gesture analysis[END_REF], have been implement using a Time-Switch Observer based on the acceleration and a world-aware observer combining the gesture performed by the user and the object pointed by the user in our installation:

The virtual hammer with click-and-point: In this technique, the dynamic of the gesture is done first, then pointing is used to determine where the click has been performed. This technique is efficient with big object users can point easily because they must stop their gesture on it.

The virtual gun using point-and-click : Here, users can point on entities to be selected and trigger a click at the pointing location. This technique is more accurate due to the dynamic done after pointing the position users want.

INTERACTIVE ART USE CASES

Our approach was demonstrated in two artistic installations created in collaboration with artists from Le Fresnoy , a digital art school from Tourcoing, France. . Both installations consist in offering to a broad audience an interactive experience where participants could interact through gestures with an artwork. Such scenario allowed to test our gestural approach in a public situation and enabled artists to take advantage of our contact-less interaction techniques. In addition, we aimed on assessing if our gestural interaction design could fit artist needs in a real artwork installation without constraining their creativity.

Damassama

The Damassama is an interactive artwork created by Leonore Mercier where the artist or the visitor can play music using a set of percussion instruments. Through contact-less gestures, the performer can interact with physical instruments mimicking a maestro conducting a real orchestra. The in- stallation consists of a set of 24 singing bowls placed radially into two rows as depicted by Figure 1. Each bowl produces a different tone, resulting from mechanical hammer's rebound which can be stopped using a muffler. The visitor can take place in the middle of the artwork and its body gestures are gathered using a Kinect camera device. To mimic the performance of a maestro, we devise a set of both pointing and semaphoric gestures such as swipe, tap or arm opening.

The virtual hammer gesture is used as a directional tap gesture from any hand to sing a given bowl. Such gesture observe the acceleration to control the hammer hitting force. While the hand remains pointing to the same location, the pointed bowl is hit repetitively. Adjacent bowls can be hit successively by swiping the hand to its location. Using each hand, the visitor can sing several bowls concurrently. If the pointing gesture is performed with both hands to a singing bowl, the muffler is triggered stopping the bowl tone. If the performer opens their arms all singing bowls are stopped at the same time. Finally, when both hands are getting closer to each other in front of the player, the last notes are replayed as maestro musical partition.

Tempo Scaduto

The second artwork installation is the Tempo Scaduto which was created by Vincent Ciciliato . It consists of an interactive movie where the visitor can control the movie storyline performing gestures in front of the screen as illustrated in Figure 3. The synopsis of the movie tries to revive both historical mafia facts and crime scenes. The visitor can access such information by playing the role of an assassin choosing their targets to evolve between historical facts while following the storyline proposed by the artist. The revival of the scenes is presented using panoramic like views where people appear and disappear. By mimicking a gun shot with their hands, the visitor accesses movie sequences explaining historical facts. Due to the characteristic of the gesture, two cameras in the installation allow the tracking of both visitor body and their hand gestures. First, a Kinect camera is placed in front of the visitor to track their arms direction. A second camera is placed aside the user looking to their hand shape to recognize specific poses as shown in Figure 5. The virtual gun gesture combines observers using both cameras and takes advantage of the affordance provided by the handgun shape pose related to the historical context of this interactive artwork as illustrated in Figure 4.

Discussion

Both interactive art installations were designed favoring the usage of pantomimic gestures familiar to the visitors in the thematic of the artwork, i.e. mimicking a maestro or a handgun child game.While we placed a notice to describe supported gestures on both installations, most of the participants preferred to observe previous visitors before taking part of the proposed interactive experience. Visitors seem to enjoy the installation taking pictures of them self while gesturing. The gestural interaction, without requiring to instrument the visitor to track its body, provides a new dimension to the artwork and to the user experience even if the thematic of the second installation might dissuade them.

Regarding our gestural based approach, it demonstrated to be flexible enough to propose a wide range of gestures from a basic set of observers. The modular approach provided by the observer pattern and its composition allows to easily define new gestures as well as it easy the integration of different approaches . For example, the second installation relies on two distinct recognition techniques and it was done by two separate research team: one specialized in vision based algorithms and the other on 3D tracked skeleton analysis. In addition, it allows to handle the interactive requirement of each installation.

Finally we noticed problems inherent to gesture recognition. First, the visitors were getting closer and closer to the camera when the recognition was not behaving correctly, making their body tracking even harder to achieve. The different granularity requirements of both hand pose detection and body tracking can be an issue as seen by the setup of the second artwork installation.

CONCLUSION AND FUTURE WORK

In this paper, we introduced a tool adapted to contactless interaction to represent parametrized gestures and a new contactless interaction technique for clicking. We validated all of them in art installation from which we collected interesting observations. As future work, a study for confirming or improving our parameter classification is a direct continuation of our reflection to adapt existing taxonomies to parameterized gesture.Our toolkit for easing interaction design implements basics for contactless and full-body interaction. However a graphic interface, similar to ICON and Flowstate, is an interesting next step.

Figure 1 :

 1 Figure 1: Damassama installation setup using 24 singing bowls actuated mechanically. A Kinect camera is placed in the center underneath the first row.

Figure 2 :

 2 Figure 2: The artist mimicking a maestro using its own artwork installation.

Figure 3 :

 3 Figure 3: Installation setup used by the Tempo Scaduto experience.

Figure 4 :

 4 Figure 4: The Tempo Scaduto author demonstrating his work during an exhibit.

Figure 5 :

 5 Figure 5: The three hand states used by the Tempo Scaduto installation.


Sarah Delcourte
email: delcourte@math.univ-lyon1.fr

Pascal Omnes
email: pascal.omnes@cea.fr

A

Keywords: Stokes equations, vorticity-velocity-pressure formulation, boundary conditions, discrete duality, finite volumes, arbitrary meshes, a priori estimates

A discrete duality finite volume discretization of the vorticity-velocity-pressure formulation of the 2D Stokes problem on almost arbitrary two-dimensional grids

Sarah Delcourte, Pascal Omnes

Introduction

Let Ω be a bounded, open, connected but not necessarily simply connected polygon of R 2 , whose boundary is denoted by Γ; we consider the numerical approximation by means of finite volumes of the solution of the Stokes equations: with given functions f , g, σ, p d and ω d and a given real number m ω , find (u, p) satisfying These conditions are written here in the case of simply connected domains, but they will be extended in the core of the article to non-simply connected domains. There are compatibility conditions between the data (g, σ) in (1.3) and (1.4), the data (ω d , σ) in (1.5) and the data (m ω , σ) in (1.6). They will be discussed in subsection 4.1.

-∆u + ∇p = f in Ω, (1.1
As recognized for example by Dubois et. al. [START_REF] Dubois | Vorticity-velocity-pressure and stream functionvorticity formulations for the Stokes problem[END_REF], these non-standard conditions can be treated in a very general and natural way thanks to the vorticity-velocity-pressure formulation of the Stokes problem (for earlier works based on different approaches, we refer to [START_REF] Bègue | Les équations de Stokes et de Navier-Stokes avec des conditions aux limites sur la pression[END_REF] and [START_REF] Girault | Incompressible Finite Element Methods for Navier-Stokes Equations with Nonstandard Boundary Conditions in R 3[END_REF]). Since -∆u = ∇ × ∇ × u -∇∇ • u, using (1.2), and introducing the vorticity ω, Eq. (1.1) may be split as

∇ × ω + ∇p = f + ∇g in Ω , (1.7) ∇ × u = ω in Ω .
(1.8)

The mathematical analysis of system (1.7)-(1.8)-(1.2) with various boundary conditions was led in several references, among which [START_REF] Amara | Vorticity-velocitypressure formulation for Navier-Stokes equations[END_REF][START_REF] Amara | Vorticity-velocity-pressure formulation for the Stokes equations[END_REF][START_REF] Azaïez | Spectral discretization of the vorticity, velocity and pressure formulation of the Navier-Stokes equations[END_REF][START_REF] Dubois | Vorticity-velocity-pressure formulation for the Stokes problem[END_REF][START_REF] Dubois | Vorticity-velocity-pressure and stream functionvorticity formulations for the Stokes problem[END_REF]. Finite element methods for the vorticityvelocity-pressure formulation were analyzed in [START_REF] Amara | Vorticity-velocitypressure formulation for Navier-Stokes equations[END_REF][START_REF] Amara | Vorticity-velocity-pressure formulation for the Stokes equations[END_REF][START_REF] Dubois | First vorticity-velocity-pressure numerical scheme for the Stokes problem[END_REF][START_REF] Abboud | A priori and a posteriori estimates for threedimensional Stokes equations with nonstandard boundary conditions[END_REF]. Spectral methods were considered in [START_REF] Azaïez | Spectral discretization of the vorticity, velocity and pressure formulation of the Navier-Stokes equations[END_REF][START_REF] Bernardi | Spectral discretization of the vorticity, velocity and pressure formulation of the Stokes equations[END_REF] and [START_REF] Proot | Least-Squares spectral elements applied to the Stokes problem[END_REF], where a least-square formulation is used.

In the present work, we shall be interested in a finite volume generalization of the Marker and Cell (MAC) scheme on very general meshes. The MAC scheme was developed initially in [START_REF] Harlow | Numerical calculations of time dependent viscous incompressible flow of fluids with a free surface[END_REF] on staggered rectangular grids and extended to the so-called covolume scheme using Delaunay-Voronoi mesh pairs, as reviewed in [START_REF] Nicolaides | Covolume methods[END_REF]. We note that the orthogonality property of these mesh pairs might be in certain cases a drawback, in particular in the context of adaptive mesh refinement. The standard MAC scheme discretizes (1.1)-(1.2), while the covolume scheme discretizes (1.7)-(1.8)-(1.2). Given a (primal) mesh, the MAC and covolume schemes use as velocity unknowns the normal components of the velocity field with respect to the edges of the control volumes, while the pressure unknowns are located at their circumcenters. In addition, the covolume scheme considers vorticity unknowns at the vertices of the primal mesh.

The generalization of the MAC scheme we propose is a new application of the "discrete duality finite volume" (DDFV) method [START_REF] Delcourte | Discrete duality finite volume method for second order elliptic problems[END_REF]. Originally developed for linear diffusion equations [START_REF] Domelevo | A finite volume method for the Laplace equation on almost arbitrary two-dimensional grids[END_REF][START_REF] Hermeline | A finite volume method for the approximation of diffusion operators on distorted meshes[END_REF][START_REF] Hermeline | Approximation of diffusion operators with discontinuous tensor coefficients on distorted meshes[END_REF], the DDFV method was extended to nonlinear diffusion [START_REF] Andreianov | Discrete duality finite volume schemes for Leray-Lions-type elliptic problems on general 2D meshes[END_REF][START_REF] Boyer | Finite volume method for 2D linear and nonlinear elliptic problems with discontinuities[END_REF][START_REF] Coudière | A 3D discrete duality finite volume method for nonlinear elliptic equations[END_REF], convection-diffusion [START_REF] Coudière | The discrete duality finite volume method for convectiondiffusion problems[END_REF], electro-cardiology [START_REF] Andreianov | Convergence of discrete duality finite volume schemes for the cardiac bidomain model[END_REF][START_REF] Coudière | A 2D/3D Discrete Duality Finite Volume Scheme. Application to ECG simulation[END_REF], drift-diffusion and energy-transport models [START_REF] Chainais-Hillairet | Discrete duality finite volume schemes for two-dimensional driftdiffusion and energy-transport models[END_REF], electro-and magnetostatics [START_REF] Delcourte | A discrete duality finite volume approach to Hodge decomposition and div-curl problems on almost arbitrary two-dimensional meshes[END_REF], electromagnetism [START_REF] Hermeline | A finite volume method for the approximation of Maxwell's equations in two space dimensions on arbitrary meshes[END_REF], and Stokes flows [START_REF] Delcourte | Développement de méthodes de volumes finis pour la mécanique des fluides[END_REF][START_REF] Krell | Stabilized DDFV Schemes for Stokes Problem with Variable Viscosity on General 2D Meshes[END_REF][START_REF] Krell | Finite volume method for general multifluid flows governed by the interface Stokes problem[END_REF][START_REF] Krell | The discrete duality finite volume method for Stokes equations on three-dimensional polyhedral meshes[END_REF]. The advantage of this covolume-like method is that it allows the use of almost arbitrary meshes, including very distorted, degenerating, or highly non-conforming meshes (see the numerical tests in [START_REF] Domelevo | A finite volume method for the Laplace equation on almost arbitrary two-dimensional grids[END_REF]). The name of the method comes from the definition of discrete gradient and divergence operators verifying a discrete Green formula.

In order to get rid of the orthogonality constraints that restrict the use of MAC and covolume schemes to certain families of meshes, the price to pay in the DDFV framework is to discretize both velocity components on the edges of the control volumes, while pressure and vorticity unknowns are associated with the centroids of the primal cells and to their vertices. Then, we integrate both components of (1.7) on cells associated with the edges (the so-called "diamond mesh") and (1.8) and (1.2) on both primal and dual cells associated to the vertices. This process uses discrete versions of the divergence, gradient and curl operators which appear in (1.7), (1.8) and (1.2). These discrete operators are known to satisfy properties analogous to those verified by the continuous operators [START_REF] Delcourte | A discrete duality finite volume approach to Hodge decomposition and div-curl problems on almost arbitrary two-dimensional meshes[END_REF]. With these properties, we show that the solution of the DDFV discretization of the Stokes equations with any of the conditions (1.3) to (1.6) can be reduced to the solution of four discrete Laplace equations involving the pressure, the vorticity and the potentials stemming from the discrete Hodge decomposition of the velocity.

Other approaches with finite volumes were developed in [START_REF] Blanc | A staggered finite volume scheme on general meshes for the generalized Stokes problem in two space dimensions[END_REF][START_REF] Eymard | A staggered finite volume scheme on general meshes for the Navier-Stokes equations in two space dimensions[END_REF][START_REF] Eymard | On a stabilized collocated finite volume scheme for the Stokes problem[END_REF][START_REF] Eymard | Convergence analysis of a colocated finite volume scheme for the incompressible Navier-Stokes equations on general 2 or 3D meshes[END_REF]. In [START_REF] Blanc | A staggered finite volume scheme on general meshes for the generalized Stokes problem in two space dimensions[END_REF][START_REF] Eymard | A staggered finite volume scheme on general meshes for the Navier-Stokes equations in two space dimensions[END_REF], velocity unknowns are located at the cell centers, and pressure unknowns at the vertices. But an orthogonality condition is also needed: the segments connecting the cell centers must be orthogonal to the edges of the mesh, which restricts this scheme to particular families of meshes. The same kind of orthogonality condition is needed in [START_REF] Eymard | On a stabilized collocated finite volume scheme for the Stokes problem[END_REF][START_REF] Eymard | Convergence analysis of a colocated finite volume scheme for the incompressible Navier-Stokes equations on general 2 or 3D meshes[END_REF], where velocity and pressure unknowns are all located at the cell centers. This choice requires some stabilization obtained by penalizing the discrete incompressibility equation. Such a penalization is also needed in the previous approaches that used the DDFV ideas [START_REF] Krell | Stabilized DDFV Schemes for Stokes Problem with Variable Viscosity on General 2D Meshes[END_REF][START_REF] Krell | Finite volume method for general multifluid flows governed by the interface Stokes problem[END_REF][START_REF] Krell | The discrete duality finite volume method for Stokes equations on three-dimensional polyhedral meshes[END_REF] because in these works, velocity unknowns are located at the cell centers and vertices, while pressure is defined on the edges (or faces in 3D) of the cells, and thus the pressure gradient is proved to be injective only on simplices, and there might be spurious pressure modes on more general meshes if no penalization is added. In the approach developed in the present article, such a penalization is not needed: The pressure gradient is computed from two pressure differences, one related to pressure unknowns located at the vertices of an edge, and another related to pressure unknowns located at two neighboring cell centers. This leads to the fact that on any mesh, a vanishing gradient implies that pressure is a constant on the primal and on the dual meshes. Moreover, in our work, the important role played in the Helmholtz-Hodge decomposition of vector fields by internal boundaries in the case of non-simply connected domains is very simply taken into account through the Euler formula that links the number of edges, cells, vertices and internal boundaries (see [START_REF] Delcourte | A discrete duality finite volume approach to Hodge decomposition and div-curl problems on almost arbitrary two-dimensional meshes[END_REF] for more details on this). This wouldn't have been the case if we had chosen the same unknown locations as in [START_REF] Krell | Stabilized DDFV Schemes for Stokes Problem with Variable Viscosity on General 2D Meshes[END_REF][START_REF] Krell | Finite volume method for general multifluid flows governed by the interface Stokes problem[END_REF][START_REF] Krell | The discrete duality finite volume method for Stokes equations on three-dimensional polyhedral meshes[END_REF].

Another important result we obtain here is that there exists a way to discretize the data f such that if g = 0 and if f is the gradient of the pressure (resp. the curl of the vorticity), then our scheme computes exactly (up to the precision of 1D quadrature formulas) the pressure and a vanishing vorticity (resp. the vorticity and a vanishing pressure) at the cell centers and vertices.

As compared to previous work with DDFV schemes, the novelty here is to realize that the vorticity-velocity-pressure formulation of the Stokes equations associated to an appropriate choice of the degrees of freedom reduces the numerical solution of this problem to four discrete DDFV Laplace equations. Another novelty is that in the course of the numerical analysis of the scheme, we are led to evaluate the difference between some functions and their pointwise values on the one hand, and their mean-value projections on the other hand, both on the primal and dual meshes. These kinds of estimations had never been used previously in the DDFV context. Moreover, we had to combine this with some Poincaré inequalities which have recently been proved [START_REF] Le | Discrete Poincaré inequalities for arbitrary meshes in the discrete duality finite volume context[END_REF].

When g = 0 in Eq. (1.2), another feature of this scheme is that it satisfies the notion of "reinforced incompressibility" introduced in a finite volume element context in [START_REF] Heib | Nouvelles discrétisations non structurées pour des écoulements de fluides à incompressibilité renforcée[END_REF] and in a finite element context in [START_REF] Bernardi | More pressure in the finite element discretization of the Stokes Problem[END_REF] to overcome non-perfectly divergence free velocity modes that may appear in unsteady Navier-Stokes simulations performed with Crouzeix-Raviart [START_REF] Crouzeix | Conforming and nonconforming finite element methods for solving the stationary Stokes equations[END_REF] finite elements of lowest order. Indeed, since these elements involve pressure unknowns located at the triangle barycenters only, the incompressibility constraint is satisfied only around these barycenters, but the discrete velocity divergence, when computed around the vertices of the mesh, may not vanish or even be small. A cure to this problem, proposed in [START_REF] Bernardi | More pressure in the finite element discretization of the Stokes Problem[END_REF] and [START_REF] Heib | Nouvelles discrétisations non structurées pour des écoulements de fluides à incompressibilité renforcée[END_REF], is to add pressure unknowns at the vertices of the mesh, which introduces incompressibility constraints around these nodes. Thus, on triangular meshes, the scheme presented here has exactly the same unknowns as (but is however not equivalent to) those involved in [START_REF] Bernardi | More pressure in the finite element discretization of the Stokes Problem[END_REF] and [START_REF] Heib | Nouvelles discrétisations non structurées pour des écoulements de fluides à incompressibilité renforcée[END_REF], and incompressibility conditions are written on each triangle (primal cell) and around each vertex (dual cell) of the mesh. The advantage of the scheme we present here is that it handles meshes that are much more general than triangular.

The convergence analysis shows that pressure, vorticity, their gradients and the velocity field are all first order accurate on general meshes. Moreover, for families of meshes where diamond cells are almost all parallelograms, pressure and vorticity gradients superconverge with the order 1.5. The numerical tests included in this work validate these theoretical results, and even go beyond: We observe that pressure and vorticity converge on general meshes with the order 2, while velocity converges with the order 2 for families of triangular meshes or for which the diamond cells are almost all parallelograms.

The paper is organized as follows: in Section 2, we introduce the notations associated with the primal, dual and diamond meshes. In Section 3, we define the discrete gradient and vector curl on the diamond cells, and the corresponding adjoint discrete divergence and scalar curl on the primal and dual cells and state their properties. In Section 4, we present the DDFV schemes for the steady Stokes problems with the various conditions (1.3)-(1.6) extended to non-simply connected domains. Section 5 is devoted to the convergence analysis of the DDFV scheme with conditions (1.3). Finally, section 6 presents numerical results over unstructured and non-conforming meshes.

Definitions and notations

Let Ω be a bounded connected polygon of R 2 , with boundary denoted by Γ. We suppose that the domain has Q holes. Throughout the paper, we shall assume that Q > 0, but the results also hold for the case Q = 0. Let Γ 0 denote the exterior boundary of Ω and let Γ q , with q ∈ [1, Q], be the interior polygonal boundaries of Ω, so that Γ = Γ 0

q∈[1,Q] Γ q .

Construction of the primal mesh

We consider a first partition of Ω (named primal mesh) composed of elements T i , with i ∈ [1, I], supposed to be convex polygons. With each element T i of the mesh, we associate a node G i located at the barycentre of T i . The area of T i is denoted by |T i |. We shall denote by J the total number of edges of this mesh and by J Γ the number of these edges which are located on the boundary Γ and we associate with each of these boundary edges its midpoint, also denoted by G i with i ∈ [I + 1, I + J Γ ]. By a slight abuse of notations, we shall write i ∈ Γ q if and only if G i ∈ Γ q .

Construction of the dual mesh

We denote by S k , with k ∈ [1, K], the nodes of the polygons of the primal mesh. Around each S k , we construct a polygon P k by joining the points G i associated to the elements of the primal mesh (and possibly to the boundary edges) of which S k is a node, through the midpoints of the edges. The area of P k is denoted by |P k |. The family of cells (P k ) k∈ [1,K] constitutes a second partition of Ω, which we name dual mesh. We shall also write k ∈ Γ q (respectively k ∈ Γ) if and only if S k ∈ Γ q (resp. S k ∈ Γ). Figure 1 displays an example of a non-conforming primal mesh and its associated dual mesh.

Pk Ti Sk

Gi Figure 1: An example of a primal mesh and its associated dual mesh.

Construction of the diamond mesh

With each edge of the primal mesh, denoted by A j (whose length is

|A j |), with j ∈ [1, J],
we associate a quadrilateral named "diamond cell", denoted by D j . When A j is not on the boundary, this cell is obtained by joining the points S k1(j) and S k2(j) , which are the two nodes of A j , with the points G i1(j) and G i2(j) associated with the elements of the primal mesh which share this edge. When A j is on the boundary Γ, the cell D j is obtained by joining the two nodes of A j with the point G i1(j) associated with the only element of the primal mesh of which A j is an edge and to the point G i2(j) associated with A j (i.e. by convention i 2 (j) is element of [I + 1, I + J Γ ] when A j is located on Γ). The cells D j constitute the "diamond-mesh", which is a third partition of Ω. The area of the cell D j is denoted by |D j |. Such cells are displayed on Fig. 2. We shall also write j ∈ Γ q (respectively j ∈ Γ) if and only if A j ⊂ Γ q (resp. j ⊂ Γ).

Definitions of geometrical elements

The following geometrical elements are displayed on Fig. 3. The unit vector normal to A j is denoted by n j and is oriented so that its dot product with G i1(j) G i2(j) is positive. We further denote by A ′ j the segment [G i1(j) G i2(j) ] (whose length is |A ′ j |) and by n ′ j the unit vector normal to A ′ j oriented so that S k1(j) S k2(j) • n ′ j ≥ 0. We also define the unit vectors τ j and τ ′ j such that (n j , τ j ) and (n ′ j , τ ′ j ) are orthogonal positively oriented bases of R 2 . We denote by M j the midpoint of A j and by A ′ j1 (respectively that A j is an edge of T i and for each k ∈ [1, K] the set E(k) of integers j ∈ [1, J] such that S k is a node of A j . We define for each j ∈ [1, J] and each k such that j ∈ E(k) (resp. each i such that j ∈ V (i)) the real-valued number s ′ jk (resp. s ji ) whose value is +1 or -1 whether n ′ j (resp. n j ) points outward or inward P k (resp. T i ). We define n ′ jk := s ′ jk n ′ j (resp. n ji := s ji n j ) and remark that n ′ jk (resp. n ji ) always points outward P k (resp. T i ). In the same way, we set n ′ jk1 := s ′ jk n ′ j1 and n ′ jk2 := s ′ jk n ′ j2 .

A ′ j2 ) the segment [G i1(j) M j ] (resp. [M j G i2(j) ]) and by n ′ j1 (resp. n ′ j2 ) the unit normal vector to A ′ j1 (resp. A ′ j2 ) oriented so that |A ′ j |n ′ j = |A ′ j1 |n ′ j1 + |A ′ j2 |n ′ j2 . We define for each i ∈ [1, I] the set V (i) of integers j ∈ [1, J] such G i 1 2 S G i 1 2 k S k j D G G 1 i S 2 k S k 1 i 2 D j Figure 2: Examples of diamond cells. i 2 G 1 S S k 1 G i 2 k A n n' j D j j j A' j S k 1 A' j G 2 i G i 1 j n' M j n' j1 A' j1 n' j2 A' j2 S k 2

A hypothesis on the mesh regularity

In the sequel of this article, we shall obtain error estimates under the following hypothesis (see Figure 4 for the notations).

Hypothesis 2.1. The angles of the subtriangulation G i M j S k of the diamond-cells D j are greater than an angle θ * which is strictly positive and independent of the mesh:

∃θ * , 0 < θ * < π 2 such that min{α m , β m , γ m } ≥ θ * , m = 1, 2, 3, 4. i G 1 G 2 S M j S k 1 i 2 k 4 β 1 γ α 1 β 1 γ 2 β 2

The discrete operators and their properties

We may approach the gradient ∇ = ∂ ∂x , ∂ ∂y T and the vector curl ∇× = ∂ ∂y , -∂ ∂x T by discrete counterparts on the diamond cells D j (see [START_REF] Domelevo | A finite volume method for the Laplace equation on almost arbitrary two-dimensional grids[END_REF][START_REF] Delcourte | A discrete duality finite volume approach to Hodge decomposition and div-curl problems on almost arbitrary two-dimensional meshes[END_REF]).

Definition 3.1. For any φ = (φ T i , φ P k ) ∈ R I+J Γ × R K , the discrete gradient ∇ D h and vector curl ∇ D
h × are defined by their values over the diamond cells D j (see Fig. 3):

(∇ D h φ) j := 1 2 |D j | φ P k2 -φ P k1 |A ′ j |n ′ j + φ T i2 -φ T i1 |A j |n j , (3.1)
(∇ D h × φ) j := - 1 2 |D j | φ P k2 -φ P k1 |A ′ j |τ ′ j + φ T i2 -φ T i1 |A j |τ j . (3.2)
Note that formulas (3.1)-(3.2) are exact for any affine function φ if we set φ P k := φ(S k ) and

φ T i := φ(G i ), for any (i, k).
Next, for a discrete vector field u, given by its values u j on the cells D j , we define its discrete divergence and curl on the primal and dual cells.

Definition 3.2. For u = (u j ) ∈ (R 2 ) J , the discrete divergence ∇ T,P h • := (∇ T h •, ∇ P h •) and discrete curl ∇ T,P h × := (∇ T h ×, ∇ P h ×) are

defined by their values over the primal cells T i and the dual

T i A j τ ji ji n jk1 n' jk1 τ ' jk2 τ jk2 n' k A' A' j2 j1 P '
Figure 5: Edges and unit vectors for the discrete divergence and curl cells P k (see Fig. 5)

(∇ T h • u) i := 1 |T i | j∈V (i) |A j |u j • n ji , (3.3) (∇ P h • u) k := 1 |P k | j∈E(k) |A ′ j1 |n ′ jk1 + |A ′ j2 |n ′ jk2 • u j + j∈E(k)∩Γ 1 2 |A j |u j • n j . (∇ T h × u) i : = 1 |T i | j∈V (i) |A j |u j • τ ji , (3.4) (∇ P h × u) k : = 1 |P k | j∈E(k) |A ′ j1 |τ ′ jk1 + |A ′ j2 |τ ′ jk2 • u j + j∈E(k)∩Γ 1 2 |A j |u j • τ j .
Remark that if the node S k is not on the boundary Γ (i.e. if k / ∈ Γ), then the set E(k) ∩ Γ is empty. On the contrary, if P k is a boundary dual cell, then the set E(k) ∩ Γ is composed of the two boundary edges which have S k as a vertex. In this case, the quantity

j∈E(k)∩Γ 1 2 |A j |u j • n j (resp. j∈E(k)∩Γ 1 2 |A j |u j • τ j ) is an approximation of ∂P k ∩Γ u • n(ξ) dξ (resp. ∂P k ∩Γ u • τ (ξ) dξ.)
For a given vector field u, it is straightforward to check that formulae

(3.3) (resp. (3.4)) are the exact mean-values of ∇•u (resp. ∇•τ ) over T i and an inner P k if |A j |u j •n ji = Aj u•n ji ds (resp. |A j |u j • τ ji = Aj u • τ ji ds) and |A ′ j1 |n ′ jk1 + |A ′ j2 |n ′ jk2 • u j = A ′ j1 u • n ′ jk1 ds + A ′ j2 u • n ′ jk2 ds (resp. |A ′ j1 |τ ′ jk1 + |A ′ j2 |τ ′ jk2 • u j = A ′ j1 u • τ ′ jk1 ds + A ′ j2 u • τ ′ jk2 ds
). Now, we define discrete scalar products and a trace operator.

Definition 3.3. For (φ, ψ) ∈ R I × R K 2 and (u, v) ∈ (R 2 ) J 2 , we set: (u, v) D := j∈[1,J] |D j | u j • v j , (3.5) (φ, ψ) T,P := 1 2 i∈[1,I] |T i | φ T i ψ T i + k∈[1,K] |P k | φ P k ψ P k . (3.6)
The associated norms will be denoted by || • || D and || • || T,P . We also define discrete scalar products on the boundaries Γ q and Γ: for sets of values (u j , v j ) j∈Γ we set:

(u, v) Γq,h := j∈Γq |A j | u j v j , (u, v) Γ,h := q∈[0,Q] (u, v) Γq,h . (3.7)
In particular, we shall often consider these boundary scalar products for the trace φ of a given φ ∈ R I+J Γ × R K , defined by Definition 3.4. (Trace operator) For any φ ∈ R I+J Γ × R K , and for any j ∈ Γ, we set

φj := 1 4 φ P k1(j) + 2φ T i2(j) + φ P k2(j) . (3.8)
With these scalar products, the discrete operators verify properties which are analogues to properties of the continuous operators [START_REF] Delcourte | A discrete duality finite volume approach to Hodge decomposition and div-curl problems on almost arbitrary two-dimensional meshes[END_REF][START_REF] Domelevo | A finite volume method for the Laplace equation on almost arbitrary two-dimensional grids[END_REF]: Proposition 3.5. The following discrete analogues of the Green formulae hold:

(∇ T,P h • u, φ) T,P = -(u, ∇ D h φ) D + (u • n, φ) Γ,h , (3.9)
(∇ T,P h × u, φ) T,P = (u, ∇ D h × φ) D + (u • τ , φ) Γ,h , (3.10)
for all u ∈ R J 2 and all φ = (φ T , φ

P ) ∈ R I+J Γ × R K .
The following two propositions state discrete equivalents of the continuous properties ∇ • (∇×) = 0, ∇ × ∇ = 0 and ∇ × (∇×) = -∇ • ∇ (in two dimensions). Note however the hypothesis needed on the boundary dual cells: Proposition 3.6. For all φ = (φ T , φ P ) ∈ R I+J Γ × R K , the following equalities hold:

∇ T,P h • (∇ D h × φ) i,k = 0, ∀i ∈ [1, I], ∀k / ∈ Γ, (3.11)
∇ T,P h × (∇ D h φ) i,k = 0, ∀i ∈ [1, I], ∀k / ∈ Γ. (3.12)
Moreover, these formulae still hold on boundary dual cells

P k , k ∈ Γ, if for any q ∈ [0, Q],
there exist two real numbers (c T q , c P q ) such that φ T i = c T q and φ P k = c P q uniformly over Γ q .

Proposition 3.7. For all φ = (φ T , φ P ) ∈ R I+J Γ × R K , the following equality holds:

(∇ T,P h × ∇ D h × φ) i,k = -(∇ T,P h • ∇ D h φ) i,k , ∀i ∈ [1, I], ∀k ∈ [1, K]. (3.13)
Finally, a discrete Helmholtz-Hodge decomposition of vector fields in non simply connected domains holds [START_REF] Delcourte | A discrete duality finite volume approach to Hodge decomposition and div-curl problems on almost arbitrary two-dimensional meshes[END_REF]: Proposition 3.8. Let (u j ) j∈[1,J] be a discrete vector field defined by its values on the diamond cells D j . Then, there exist unique φ = (φ T i , φ P k ) and ψ = (

ψ T i , ψ P k ) both in R I+J Γ × R K and (c T q , c P q ) in R Q × R Q such that: u j = (∇ D h φ) j + (∇ D h × ψ) j , ∀j ∈ [1, J], (3.14
)

with i∈[1,I] |T i |φ T i = k∈[1,K] |P k |φ P k = 0 , (3.15)
ψ T i = 0 , ∀i ∈ Γ 0 , ψ P k = 0 , ∀k ∈ Γ 0 , (3.16) ∀q ∈ [1, Q] , ψ T i = c T q , ∀i ∈ Γ q , ψ P k = c P q , ∀k ∈ Γ q . (3.17)
Moreover, the decomposition (3.14) is orthogonal, in the sense that

(∇ D h φ, ∇ D h × ψ) D = 0, (3.18)
a property that is true for all φ and all ψ verifying (3. 16) and (3.17).

Formulae (3.15) are discrete analogues (respectively stated on the primal mesh and on the dual mesh) of the condition Ω φ = 0 that appears in the continuous Helmholtz-Hodge decomposition, while formulae (3.16) and (3.17) are discrete analogues of the boundary conditions ψ = 0 on Γ 0 and ψ = c q on Γ q . Remark 3.9. We may also write a similar decomposition by changing the conditions (3.15)- (3.17

|T i |ψ T i = k∈[1,K] |P k |ψ P k = 0 , (3.19)
φ T i = 0 , ∀i ∈ Γ 0 , φ P k = 0 , ∀k ∈ Γ 0 , (3.20) ∀q ∈ [1, Q] , φ T i = c T q , ∀i ∈ Γ q , φ P k = c P q , ∀k ∈ Γ q . (3.21)
Finally, we mention discrete analogues of the so-called Poincaré inequalities, which were proved in [START_REF] Le | Discrete Poincaré inequalities for arbitrary meshes in the discrete duality finite volume context[END_REF] (respectively Theorem 3.2 and Theorem 3.6): (3.15), (3.16) and (3.17) hold; then, there exists a constant C(θ * ), depending only on θ * and on Ω such that

Proposition 3.10. Let φ = (φ T i , φ P k ) and ψ = (ψ T i , ψ P k ) both in R I+J Γ × R K and (c T q , c P q ) in R Q × R Q be such that
||φ|| T,P ≤ C(θ * )||∇ D h φ|| D , ||ψ|| T,P ≤ C(θ * )||∇ D h ψ|| D . ( 3
u • n = σ over Γ , ω |Γ0 = ω d ; ω |Γq = ω d + c q , ∀ q ∈ [1, Q], Γq u • τ = k q , ∀ q ∈ [1, Q], Ω p(x)dx = 0 , (4.1)
u • n = σ over Γ , p |Γ0 = p d ; p |Γq = p d + c q , ∀ q ∈ [1, Q], Γq u • τ = k q , ∀ q ∈ [1, Q], Ω ω(x)dx = m ω , (4.2)
u • τ = σ over Γ , ω |Γ0 = ω d ; ω |Γq = ω d + c q , ∀ q ∈ [1, Q], Γq u • n = k q , ∀ q ∈ [1, Q], Ω p(x)dx = 0 , (4.3)
u • τ = σ over Γ , p |Γ0 = p d ; p |Γq = p d + c q , ∀ q ∈ [1, Q], Γq u • n = k q , ∀ q ∈ [1, Q], Ω ω(x)dx = m ω , (4.4)
where the constants (c q ) ∈ R Q have to be determined. Before going into details of the discretization of these equations, we discuss the compatibility conditions of the right-hand sides of these sets of equations.

Compatibility conditions

Consider Eqs. (1.2), (1.7) and (1.8). First, when associated with one of the sets of conditions (4.1) or (4.2), the data has to verify

Ω g(x)dx = Γ σ(ξ)dξ. (4.5)
Then, when associated with the set (4.4), integration of (1.8) over Ω and application of the Green formula and of the boundary conditions yield

Γ σ(ξ)dξ = m ω . (4.6)
Finally, in the case of Eqs. (1.2), (1.7) and (1.8) associated with the set (4.3), the situation is more involved. Indeed, we infer from (1.8) that Γ σ(ξ)dξ = Ω ω(x)dx. However, the quantity Ω ω(x)dx is not a data given by (4.3), but is a result of the computation of ω through the Hodge decomposition of f + ∇g expressed by Eq. (1.7) associated with the conditions over p and ω expressed in the set (4.3). This may be interpreted as an implicit compatibility condition between the boundary conditions ω d and σ in (4.3). Further details on how to handle this will be given in subsection 4.2.

Discretization of the Stokes equations in vorticity-velocity-pressure formulation

In this subsection, we are interested in the approximation of the continuous problem given by (1.2)-(1.7)-(1.8) associated with one of the sets of conditions (4.1) to (4.4). We choose to approach the solution of this problem by a vector (u j ), with j ∈ [1, J], which discretizes the velocity field by values defined over the diamond cells of the mesh, and by scalars (ω T i , ω P k ) and

(p T i , p P k ), with i ∈ [1, I + J Γ ], k ∈ [1, K],
which discretize the vorticity and the pressure fields by values defined over the primal and dual cells, as well as boundary edges of the mesh. The problem will be solved in two steps. In the first, we use the Hodge decomposition of f + ∇g (see prop. 3.8) to solve for p and ω. In the second, we solve a div-curl problem for u.

Step 1: The discrete Hodge decomposition of the data f + ∇g reads:

find p = (p T i , p P k ) i∈[1,I+J Γ ], k∈[1,K] , ω = (ω T i , ω P k ) i∈[1,I+J Γ ], k∈[1,K] and (c T q , c P q ) q∈[1,Q] such that (∇ D h p) j + (∇ D h × ω) j = f D j + (∇g) D j , ∀j ∈ [1, J], (4.7)
together with one of the following sets of conditions

       ω T i = ω d (G i ) , ∀i ∈ Γ 0 , ω T i = ω d (G i ) + c T q , ∀i ∈ Γ q , ∀q , ω P k = ω d (S k ) , ∀k ∈ Γ 0 , ω P k = ω d (S k ) + c P q , ∀k ∈ Γ q , ∀q , i∈[1,I] |T i | p T i = k∈[1,K] |P k | p P k = 0, (4.8)
for a given (up to constants to be determined on each internal boundary) vorticity field ω d on the boundary (see the corresponding equations in (4.1) and (4.3)), or

       p T i = p d (G i ) , ∀i ∈ Γ 0 , p T i = p d (G i ) + c T q , ∀i ∈ Γ q , ∀q , p P k = p d (S k ) , ∀k ∈ Γ 0 , p P k = p d (S k ) + c P q , ∀k ∈ Γ q , ∀q , i∈[1,I] |T i | ω T i = k∈[1,K] |P k | ω P k = m ω , (4.9)
for a given (up to constants to be determined on each internal boundary) pressure field p d on the boundary (see the corresponding equations in (4.2) and (4.4)).

In (4.7), the following definitions have been used

f D j = 1 |D j | Dj f (x) dx , (∇g) D j = 1 |D j | Dj ∇g (x) dx ∀j ∈ [1, J]. ( 4
:= (s D j ) j∈[1,J] with s D j := f D j + (∇g) D j , find (ω T i , ω P k ) i∈[1,I+J Γ ],k∈[1,K] and (c T q , c P q ) q∈[1,Q] such that                -(∇ T h • ∇ D h ω) i = (∇ T h × s) i , ∀i ∈ [1, I], -(∇ P h • ∇ D h ω) k = (∇ P h × s) k , ∀k / ∈ Γ, (∇ D h ω • n, 1) Γq ,h = -(s • τ , 1) Γq,h , ∀q ∈ [1, Q], -k∈Γq |P k | (∇ P h • ∇ D h ω) k = k∈Γq |P k | (∇ P h × s) k , ∀q ∈ [1, Q], ω T i = ω d (G i ), ∀i ∈ Γ 0 , ω P k = ω d (S k ), ∀k ∈ Γ 0 , ω T i = ω d (G i ) + c T q , ∀i ∈ Γ q , ω P k = ω d (S k ) + c P q , ∀k ∈ Γ q , ∀q ∈ [1, Q] , (4.11)
and, once ω has been computed, find

(p T i , p P k ) i∈[1,I+J Γ ],k∈[1,K] such that            (∇ T h • ∇ D h p) i = (∇ T h • s) i , ∀i ∈ [1, I], (∇ P h • ∇ D h p) k = (∇ P h • (s -∇ D h × ω)) k , ∀k ∈ [1, K], (∇ D h p) j • n j = (s j -(∇ D h × ω) j )) • n j , ∀j ∈ Γ, i∈[1,I] |T i | p T i = k∈[1,K] |P k | p P k = 0.
(4.12)

Proof. Applying the discrete vector curl operator to (4.7) on any primal cell and on any inner dual cell yields the first two lines of (4.11), thanks to (3.12) and (3.13).

Next, for a given q ∈ [1, Q], we consider the element ψ ∈ R I+J Γ × R K which has the following values: ψ T i = 1, ∀i ∈ Γ q and ψ T i = 0 everywhere else, and ψ P k = 0 everywhere. Then, we compute the scalar product (3.5) of Eq. (4.7) with ∇ D h × ψ:

(∇ D h × ω, ∇ D h × ψ) D + (∇ D h p, ∇ D h × ψ) D = (s, ∇ D h × ψ) D . (4.13)
Using the orthogonality of ∇ D h p and ∇ D h × ψ (see last line of Prop. 3.8), using the discrete Green formula (3.10) and Eq. (3.13), we infer

-(∇ T,P h • ∇ D h ω, ψ) T,P -(∇ D h × ω • τ , ψ) Γ,h = (∇ T,P h × s, ψ) T,P -(s • τ , ψ) Γ,h . (4.14)
Further, using the fact that ψ vanishes everywhere but on the boundary points G i , the first term in the left-hand side and the first term in the right-hand side of (4.14) vanish. Next, since

∇ D h × ω • τ = -∇ D h ω
• n, and using the definition (3.7) and the fact that the values of ψ on the boundaries imply that ψj = 1/2 for all j ∈ Γ q and ψj = 0 for all j ∈ Γ ′ q for any q ′ = q, Equ. ( 4

.14) implies -(∇ D h ω • n, 1 2 ) Γq,h = (s • τ , 1
2 ) Γq,h and thus the third line of (4.11). Last, we consider the element ψ ∈ R I+J Γ × R K which has the following values: ψ T i = 1, ∀i ∈ Γ q and ψ P k = 1, ∀k ∈ Γ q and vanishing everywhere else. Then, the dot product of Eq. Once ω has been computed, the derivation of (4.12) from (4.7)-(4.8) is obvious. We stress that in the second equation of (4.12), there holds (∇ P h • ∇ D h × ω) k = 0 for all inner dual cells (k / ∈ Γ) but that this property might not necessarily be true for boundary dual cells (k ∈ Γ), see Prop. 3.6.

Step 1 ends with the fact that (4.11) and (4.12) are well-posed: Lemma 4.2. It was shown in [START_REF] Delcourte | A discrete duality finite volume approach to Hodge decomposition and div-curl problems on almost arbitrary two-dimensional meshes[END_REF]Proposition 5.2] that systems of types (4.11) and (4.12) both have a unique solution.

Now, let us describe

Step 2 of the calculations, which consists in solving for the velocity.

Step 2: Once (ω T i , ω P k ) i∈[1,I],k∈ [1,K] has been computed through Step 1, we solve a div-curl problem for u:

given (k q ) q∈[1,Q] , find (u j ) j∈[1,J] such that                  (∇ T,P h • u) i,k = g T,P i,k , ∀i ∈ [1, I], ∀k ∈ [1, K], (∇ T,P h × u) i,k = ω T,P i,k , ∀i ∈ [1, I], ∀k / ∈ Γ, u j • n j = σ j , ∀j ∈ Γ, (u • τ , 1) Γq,h = k q , ∀q ∈ [1, Q], k∈Γq |P k | (∇ P h × u) k = k∈Γq |P k | ω P k , ∀q ∈ [1, Q], (4.15)
in the case of a given normal velocity field on the boundary (Eqs. (4.1) or (4.2)) or

             (∇ T,P h • u) i,k = g T,P i,k , ∀i ∈ [1, I], ∀k / ∈ Γ, (∇ T,P h × u) i,k = ω T,P i,k , ∀i ∈ [1, I], ∀k ∈ [1, K], u j • τ j = σ j , ∀j ∈ Γ, (u • n, 1) Γq,h = k q , ∀q ∈ [1, Q], k∈Γq |P k | (∇ P h • u) k = k∈Γq |P k | g P k , ∀q ∈ [1, Q], (4.16)
in the case of a given tangential velocity field on the boundary (Eqs. (4.3) or (4.4)). In (4.15) and (4.16), we have set

g T i = 1 |T i | Ti g(x) dx ∀i ∈ [1, I] , g P k = 1 |P k | P k g(x) dx ∀k ∈ [1, K], (4.17)
σ j = 1 |A j | Aj σ(ξ) dξ ∀j ∈ Γ. ( 4
|T i |(∇ T h • u) i = k∈[1,K] |P k |(∇ P h • u) k = j∈Γ |A j |u j • n j , (4.19
)

i∈[1,I] |T i |(∇ T h × u) i = k∈[1,K] |P k |(∇ P h × u) k = j∈Γ |A j |u j • τ j . ( 4
|T i |g T i = k∈[1,K] |P k |g P k = j∈Γ |A j |σ j . (4.21)
This relation is true thanks to the definitions (4.17

|T i |g T i = k∈[1,K] |P k |g P k = Ω g(x) dx and j∈Γ |A j |σ j = Γ σ(ξ) dξ .
That the right-hand sides of the above two equalities are identical follows from (4.5). Further, because of (4.20), the right-hand sides in (4.16) must satisfy

i∈[1,I] |T i |ω T i = k∈[1,K] |P k |ω P k = j∈Γ |A j |σ j . (4.22)
For the set of conditions (4.9), the first two terms in (4.22) are equal to m ω and the last term in (4.22) is equal to Γ σ(ξ) dξ. These two quantities are identical thanks to (4.6).

On the other hand, for the set of conditions (4.8), the values of the first two terms in (4.22) are never imposed, but, rather, are results of the computations involved in the first step of our procedure, see Eq. (4.11), so that the compatibility condition (4.22) may not be verified in general. A possible way to overcome this problem is to change ω T i into ω T i + c T and ω P k into ω P k +c P for all i ∈ [1, I +J Γ ] and all k ∈ [1, K] in (4. [START_REF] Hermeline | A finite volume method for the approximation of diffusion operators on distorted meshes[END_REF]), with two constants c T and c P computed so that (4.22) holds. Doing so does not change the value of ∇ D h ω, so that the modified ω still verifies system (4.11), but with modified boundary conditions. This may be interpreted as an implicit compatibility condition between the boundary conditions ω d and σ in (4.3). Now, using the discrete Hodge decomposition of (u j ) j∈[1,J] , each of the problems (4.15) and (4.16) may be split into two independent subproblems involving the potentials. We only detail the resulting systems for problem (4.15) using the Hodge decomposition with boundary conditions (3.15)

(φ T i , φ P k ) i∈[1,I+J Γ ],k∈[1,K] such that            (∇ T h • ∇ D h φ) i = g T i , ∀i ∈ [1, I], (∇ P h • ∇ D h φ) k = g P k , ∀k ∈ [1, K], (∇ D h φ) j • n j = σ j , ∀j ∈ Γ, i∈[1,I] |T i | φ T i = k∈[1,K] |P k | φ P k = 0, (4.23
)

and find (ψ T i , ψ P k ) i∈[1,I+J Γ ],k∈[1,K] and (c T q , c P q ) q∈[1,Q] such that                    -(∇ T h • ∇ D h ψ) i = ω T i , ∀i ∈ [1, I], -(∇ P h • ∇ D h ψ) k = ω P k , ∀k / ∈ Γ, -(∇ D h ψ • n, 1) Γq ,h = k q , ∀q ∈ [1, Q], - k∈Γq |P k | (∇ P h • ∇ D h ψ) k = k∈Γq |P k | ω P k , ∀q ∈ [1, Q], ψ T i = ψ P k = 0, ∀i ∈ Γ 0 , ∀k ∈ Γ 0 , ∀q ∈ [1, Q], ψ T i = c T q , ψ P k = c P q , ∀i ∈ Γ q , ∀k ∈ Γ q . (4.24)
Proof. The proof is given in [25, Proposition 5.1].

Step 2 ends with the fact that (4.23) and (4.24) are well-posed: Once these two subproblems have been solved, the vector u is then reconstructed by u j = (∇ D h φ) j + (∇ D h × ψ) j .

Error estimates

Obtaining error estimates usually relies on regularity assumptions on the solution of the problem.

In order to apply results given in [START_REF] Delcourte | A discrete duality finite volume approach to Hodge decomposition and div-curl problems on almost arbitrary two-dimensional meshes[END_REF][START_REF] Domelevo | A finite volume method for the Laplace equation on almost arbitrary two-dimensional grids[END_REF], we shall assume more regularity on the vorticity and pressure fields (ω, p) given by Proposition 5.1 below, and on the velocity potentials given by Proposition 5.3 below.

Proposition 5.1. Let (f , g, σ, ω d ) belong to L 2 (Ω) 2 × H 1 (Ω) × H 1/2 (Γ) × H 1/2 (Γ)
, and let

(k q ) q∈[1,Q] be
(C q ) q∈[1,Q]
are the exact solution of the Hodge decomposition of f + ∇g:

   ∇ × ω + ∇p = f + ∇g in Ω, ω|Γ0 = ω d ; ω|Γq = ω d + C q ∀ q ∈ [1, Q], Ω p(x)dx = 0, (5.1)
and, once ω has been determined, û ∈ Hdiv(Ω) ∩ Hcurl(Ω) solves the div-curl problem:

       ∇ • û = g in Ω, ∇ × û = ω in Ω, û • n = σ on Γ, Γq û • τ = k q , ∀q ∈ [1, Q].
(5.2) Hypothesis 5.2. We assume that ω and p given by Proposition 5.1 belong to H 2 (Ω).

The velocity field û solving (5.2) may be found by the following Hodge decomposition.

Proposition 5.3. Let (g, ω, σ) belong to H 1 (Ω) × H 2 (Ω) × H 1/2 (Γ), and let (k q ) q∈[1,Q] be a set of given real numbers; let û be the exact solution of problem (5.2). Then, there exist φ and ψ both in H 1 (Ω) and a set of real numbers

(C q ) q∈[1,Q] such that û = ∇ φ + ∇ × ψ,
where φ and ψ are the solutions of

   ∆ φ = ∇ • û = g in Ω, ∇ φ • n = û • n = σ on Γ, Ω φ = 0,

and

(5.3)

     -∆ ψ = ∇ × û = ω in Ω, ψ|Γ0 = 0; ψ|Γq = C q ∀ q ∈ [1, Q], Γq ∇ ψ • n = -k q .
(5.4)

Proof. The Hodge decomposition of û and the determination of φ and ψ through (5.3) and (5.4) are direct consequences of [39, Theorem 3.2 and Corollary 3.1].

Hypothesis 5.4. We suppose that φ and ψ given by Proposition 5.3 belong to H 2 (Ω).

Due to re-entrant corners related to the internal polygonal boundaries Γ q , the H 2 regularity of the potentials is not a consequence of the regularity of the data (g, ω, σ).

Error estimates in the

I + J Γ ], (Πv) T i = v(G i ), ∀k ∈ [1, K], (Πv) P k = v(S k ).
We note that, under Hypothesis 5.2, ω and p belong to H 2 (Ω), which implies they are continuous on Ω, so that considering Πω and Πp makes sense.

Theorem 5.6. Let p and ω be the solutions of the numerical scheme (4.7)- (4.8), and let (p, ω) be the exact solutions of (5.1). Then, if all diamond-cells are convex and under assumptions 2.1 and 5.2, there exists a constant C(θ * ), independent of h, such that

||∇ D h (p -Πp)|| D + ||∇ D h (ω -Πω)|| D ≤ C(θ * ) h ( p 2,Ω + ω 2,Ω ). (5.5)
Proof. Firstly, taking the mean-value of each term of the first line of Equ. (5.1) on a diamond cell D j and using (4.10), we get:

1 |D j | Dj (∇p + ∇ × ω)(x)dx = f D j + (∇g) D j , ∀j ∈ [1, J]. (5.6)
Then, since (5.6) and (4.7) have the same right-hand side, we infer that:

(∇ D h p) j + (∇ D h × ω) j = 1 |D j | Dj (∇p + ∇ × ω)(x)dx, ∀j ∈ [1, J]. ( 5
(∇ D h ε p ) j + (∇ D h × ε ω ) j = 1 |D j | Dj (∇p(x) -(∇ D h Πp) j )dx + 1 |D j | Dj (∇ × ω(x) -(∇ D h × Πω) j )dx, ∀j ∈ [1, J].
(5.9)

According to (4.8) and the second line of (5.1), the error ε ω satisfies on the boundary

(ε ω ) T i = 0 , ∀i ∈ Γ 0 , (ε ω ) T i = c T q -C q , ∀i ∈ Γ q , ∀q , (ε ω ) P k = 0 , ∀k ∈ Γ 0 , (ε ω ) P k = c P q -C q , ∀k ∈ Γ q , ∀q .
(5.10) Therefore, (∇ D h ε p ) and (∇ D h × ε ω ) are orthogonal for the scalar product (•, •) D , as recalled at the end of Proposition 3.8. Consequently, multiplying (5.9) by |D j |(∇ D h ε p ) j and summing for all j ∈ [1, J], we obtain:

||∇ D h ε p || 2 D = j∈[1,J] Dj (∇p(x) -(∇ D h Πp) j ) • (∇ D h ε p ) j dx (5.11) + j∈[1,J] Dj (∇ × ω(x) -(∇ D h × Πω) j ) • (∇ D h ε p ) j dx.
In the very same way, multiplying (5.9) by |D j |(∇ D h × ε ω ) j and summing for all j ∈ [1, J], we obtain:

||∇ D h ε ω || 2 D = j∈[1,J] Dj (∇p(x) -(∇ D h Πp) j ) • (∇ D h × ε ω ) j dx (5.12) + j∈[1,J] Dj (∇ × ω(x) -(∇ D h × Πω) j ) • (∇ D h × ε ω ) j dx.
The right-hand sides of (5.11) and (5.12) may be bounded using the traditional P 1 Lagrange interpolations of p and ω on a submesh of the diamond mesh, obtained by splitting each diamond cell, along any of its diagonals, into two triangles. The details of the calculations may be found in [START_REF] Domelevo | A finite volume method for the Laplace equation on almost arbitrary two-dimensional grids[END_REF], starting with inequality (62) of that reference, and then using Lemma 5.11, where the norm in the right-hand side of inequality (66) has to be replaced by the H 2 (Ω) norm of p and ω.

Error estimates in the L 2 -norm for the pressure and the vorticity

We shall use the discrete Poincaré inequalities recalled in Proposition 3.10 to infer error estimates for the discrete L 2 norm of the errors in the pressure and the vorticity. On the one hand, Proposition 3.10 can be applied directly to the error ε ω because it verifies (5.10), which are exactly conditions (3.16) and (3.17) that Proposition 3.10 requires. Thus, we infer from (5.5) the following theorem Theorem 5.7. There exists a constant C(θ * ) that does not depend on h such that ||ω -Πω|| T,P ≤ C(θ * ) h ( p 2,Ω + ω 2,Ω ).

(5.13)

On the other hand, since (p -Πp) does not in general verify the vanishing mean-value condition (3.15) (because Πp doesn't), we may not apply the discrete Poincaré inequality in a straightforward way. However, defining

c T := i∈[1,I] |T i |(Πp) T i |Ω| , ( Πp) T i := (Πp) T i -c T , ∀i ∈ [1, I + J Γ ],
(5.14)

c P := k∈[1,K] |P k |(Πp) P k |Ω| , ( Πp) P k := (Πp) P k -c P , ∀k ∈ [1, K], (5.15)
we have that Πp verifies the vanishing mean-value condition (3.15) and

∇ D h ( Πp) = ∇ D h (Πp) (5.16)
because ( Πp) and (Πp) only differ by a constant on the primal and dual meshes, and definition (3.1) of the discrete gradient implies that this does not affect its values. Further, ||p -Πp|| T,P ≤ ||p -Πp|| T,P + ||Πp -Πp|| T,P .

(5.17)

Estimating the first term in the right-hand-side of (5.17) may be performed through the discrete Poincaré inequality applied to (p -Πp).

Lemma 5.8. There exists a constant C(θ * ) that does not depend on h such that

||p -Πp|| T,P ≤ C(θ * ) h ( p 2,Ω + ω 2,Ω ). (5.18)
Proof. From Prop. 3.10, it holds that

||p -Πp|| T,P ≤ C(θ * )||∇ D h (p -Πp)|| D = C(θ * )||∇ D h (p -Πp)|| D (5.19)
because of (5.16). Applying (5.5), we get (5.18).

To estimate the second term in the right-hand-side of (5.17), we first recall (5.14) and (5.15). This leads to ||Πp -

Πp|| 2 T,P = |Ω| 2 (c T ) 2 + (c P ) 2 . (5.20)
Lemma 5.9. Let c T be defined by (5.14); then

|c T | ≤ |Ω| -1/2 i ||p -(Πp) T i || 2 L 2 (Ti) 1/2 . (5.21)
Proof. From (5.14), we have that

|Ω|c T = i |T i |(Πp) T i = i Ti
(Πp) T ip (x)dx because i Ti p(x)dx = Ω p(x)dx = 0. Then, a continuous and a discrete Cauchy-Schwarz inequality lead to (5.21) through

|Ω||c T | ≤ i |T i | 1/2 ||p -(Πp) T i || L 2 (Ti) ≤ |Ω| 1/2 i ||p -(Πp) T i || 2 L 2 (Ti) 1/2 .
So, what remains to evaluate is the L 2 norm of (Πp) T ip on T i . Lemma 5.10. There exists a constant C(θ * ) that does not depend on h such that

||p -(Πp) T i || L 2 (Ti) ≤ C(θ * )h||p|| 2,Ti . (5.22) 19
Proof. First, we split the integral on the subtriangles t i,j,k with vertices G i M j S k (see Fig. 4), where M j is the midpoint of the edge A j ⊂ ∂T i and S k one of its vertices:

Ti (Πp) T i -p 2 (x)dx = t i,j,k ⊂Ti t i,j,k (Πp) T i -p 2 (x)dx. (5.23)
Let π i,j,k p be the standard P 1 Lagrange interpolation of p on t i,j,k . We have

||p -(Πp) T i || 2 L 2 (t i,j,k ) ≤ 2||(Πp) T i -π i,j,k p|| 2 L 2 (t i,j,k ) + 2||π i,j,k p -p|| 2 L 2 (t i,j,k ) . (5.24)
It is a standard result (see Theorem 3.1.5 of [START_REF] Ciarlet | The finite element method for elliptic problems[END_REF]) that there exists a constant C, not depending on t i,j,k such that

||π i,j,k p -p|| L 2 (t i,j,k ) ≤ Cdiam(t i,j,k ) 2 ||p|| H 2 (t i,j,k ) , (5.25)
which evaluates the second term in the right-hand side of (5.24). As far as the first term in the right-hand side of (5.24) is concerned, since (Πp) T iπ i,j,k p(G i ) = 0, and since π i,j,k p is a P 1 function, we have that, for all x ∈ t i,j,k

(Πp) T i -π i,j,k p(x) = ∇π i,j,k p • (G i -x), which ensures that |(Πp) T i -π i,j,k p(x)| ≤ diam(t i,j,k )||∇π i,j,k p|| for all x ∈ t i,j,k . Thus ||π i,j,k p -(Πp) T i || L 2 (t i,j,k ) ≤ diam(t i,j,k )||∇π i,j,k p|| L 2 (t i,j,k ) . (5.26)
Using the triangular inequality

||∇π i,j,k p|| L 2 (t i,j,k ) ≤ ||∇(π i,j,k p -p)|| L 2 (t i,j,k ) + ||∇p|| L 2 (t i,j,k )
and using again Theorem 3.1.5 of [START_REF] Ciarlet | The finite element method for elliptic problems[END_REF], we have that there exists a constant C depending only the regularity parameter of the subtriangulation (and thus on θ * ) such that

||∇π i,j,k p|| L 2 (t i,j,k ) ≤ (1 + Cdiam(t i,j,k ))||p|| H 2 (t i,j,k ) . (5.27)
Gathering (5.23), (5.24), (5.25), (5.26) and (5.27) implies (5.22).

We are now able to estimate the second term in the right-hand side of (5.17)

Proposition 5.11. There exists a constant C(θ * ) that does not depend on h such that

||Πp -Πp|| T,P ≤ C(θ * ) h ||p|| 2,Ω . (5.28)
Proof. Bounds (5.21) and (5.22) imply that

|c T | ≤ C(θ * )h||p|| 2,Ω . (5.29)
Using (5.20), (5.29) and a similar bound that can be obtained in the same way for c P , we obtain (5.28).

Finally, using (5.17), (5.18) and (5.28), we obtain the following theorem:

Theorem 5.12. There exists a constant C(θ * ) that does not depend on h such that

||p -Πp|| T,P ≤ C(θ * ) h ( p 2,Ω + ω 2,Ω ).
(5.30)

Error estimate in the L 2 -norm for the velocity

Once the numerical approximation ω = (ω T i , ω P k ) i∈[1,I+J Γ ], k∈ [1,K] of the vorticity ω is known, we have to solve (4.15) in order to find the discrete velocity u = (u j ) j∈ [1,J] , with which we define a piecewise constant function u h by

u h (x) = u j , ∀x ∈ D j , ∀j ∈ [1, J].
In order to estimate the error between u h and the exact velocity û, we first introduce an intermediate discrete velocity ũ = (ũ j ) j∈ [1,J] , solution of the discrete div-curl system

                 (∇ T,P h • ũ) i,k = g T,P i,k , ∀i ∈ [1, I], ∀k ∈ [1, K], (∇ T,P h × ũ) i,k = ω T,P i,k , ∀i ∈ [1, I], ∀k / ∈ Γ, ũj • n j = σ j , ∀j ∈ Γ, (ũ • τ , 1) Γq,h = k q , ∀q ∈ [1, Q], k∈Γq |P k | (∇ P h × ũ) k = k∈Γq |P k | ω P k , ∀q ∈ [1, Q], (5.31)
where we have set

ω T i := 1 |T i | Ti ω(x) dx ∀i ∈ [1, I] and ω P k := 1 |P k | P k ω(x) dx ∀k ∈ [1, K].
Note that ũ is of course never actually computed (because the exact vorticity ω is not known), but only serves for theoretical reasons. The following triangle inequality holds

||u h -û|| 0,Ω ≤ ||u -ũ|| D + ||ũ h -û|| 0,Ω . (5.32)
A bound of the second term in the right-hand side of (5. (5.33)

The next step is the evaluation of the difference d u := uũ, which, using (4.15) and (5.31), and setting d ω := ωω , is the solution of The first term in the right-hand side of (5.40) is bounded by (5.13). Bounding the second term requires some more analysis.

               (∇ T,P h • d u ) i,k = 0, ∀i ∈ [1, I], ∀k ∈ [1, K], (∇ T,P h × d u ) i,k = (d ω ) T,P i,k , ∀i ∈ [1, I], ∀k / ∈ Γ, (d u ) j • n j = 0, ∀j ∈ Γ, (d u • τ , 1) Γq,h = 0, ∀q ∈ [1, Q], k∈Γq |P k | (∇ P h × d u ) k = k∈Γq |P k | (d ω ) P k , ∀q ∈ [1, Q], (5.34
|T i |(d ω ) T i (ψ d ) T i + 1 2 k / ∈Γ |P k |(d ω ) P k (ψ d ) P k + 1 2 k∈Γ0 |P k |(∇ T,P h × d u ) P k (ψ d ) P k + q∈[1,Q] 1 2 k∈Γq |P k |(∇
Proposition 5.15. Under hypotheses 2.1 and 5.2, there exists a constant C(θ * ) independent of h such that || ω -Πω|| T,P ≤ C(θ * ) h ω 2,Ω .

(5.41)

Proof. First, we have

|| ω -Πω|| 2 T,P = 1 2 i |T i | ω T i -(Πω) T i 2 + k |P k | ω P k -(Πω) P k 2 .
(5.42)

Let us consider the first sum in the right-hand side of (5.42); the second sum will be treated in the same way, with however an important modification that will be underlined within the developments that follow.

Since ω T i -(Πω) T i is a constant on T i , we have

|T i | ω T i -(Πω) T i 2 = Ti ω T i -(Πω) T i 2 dx ≤ 2 ||ω -ω T i || 2 L 2 (Ti) + ||ω -(Πω) T i || 2 L 2 (Ti) .
(5.43)

Using arguments similar to those that led to (5.22), we have that

||ω -(Πω) T i || L 2 (Ti) ≤ Ch||ω|| 2,Ti . (5.44)
Since T i is an open bounded set which is star-shaped with respect to G i , there exists a constant C(T i ) only depending on the shape of T i , but not on its diameter such that

ω -ω T i L 2 (Ti) ≤ C(T i )diam(T i ) ∇ω L 2 (Ti) .
(5.45)

Since T i is convex, a universal constant C(T i ) is given by 1 π , see [START_REF] Payne | An optimal Poincaré inequality for convex domain[END_REF]. Since any dual cell P k is star-shaped with respect to S k , a similar inequality holds:

ω -ω P k L 2 (P k ) ≤ C(P k )diam(P k ) ∇ω L 2 (P k ) .
(5.46) However, since P k is not necessarily convex, the estimation of C(P k ) is less obvious, but we may use explicitly computable formulas given, for example, by [START_REF] Carstensen | Constants in Clément-interpolation error and residual based a posteriori estimates in finite element methods[END_REF][START_REF] Verfürth | Error estimates for some quasi-interpolation operators[END_REF]. These formulas show that C(P k ) only depend on the angles of the subtriangulation mentioned in Hyp. 2.1. Gathering (5.43), (5.44) and (5.45), and similar inequalities on the dual cells P k , the upper bound (5.41) is obtained from (5.42). Now, we are able to estimate the first term in the right-hand side of (5.32). With (5.40), (5.13) and (5.41), we obtain from (5.35) Families of uniformly refined meshes of rectangles, or any (possibly non-conforming) union of such meshes verify such an hypothesis; see for example on Fig. 7 the second family of meshes used in section 6. In that case, each mesh of rectangles is an Ω ℓ . Another example of such families of meshes is what is called "homothetically refined triangular meshes", obtained from a coarse triangular mesh by iteratively refining each triangle into four homothetic sub-triangles by joining the midpoints of its edges. In that case, each triangle of the coarse mesh is an Ω ℓ .

Numerical results

In this section, we first test the finite volume method applied to the vorticity-velocity-pressure formulation against known analytical solutions and plot convergence curves (in log-log scale) for several quantities; then we consider the Bercovier-Engelman [START_REF] Bercovier | A finite element for the numerical solution of viscous incompressible flows[END_REF] as well as the Ruas [START_REF] Ruas | Méthodes vitesse-tourbillon pour les équations de Navier-Stokes[END_REF] test case.

Convergence curves

First, the data are chosen so that the exact solution is û = (exp(x) cos(πy), x sin(πy) + cos(πx))

T and p = xy exp(x) cos(πy).

(6.1)

The boundary conditions are those given by (4.1). Three families of increasingly fine meshes are used. The first is a family of standard triangular meshes, used on the non-simply connected domain Ω =] -1/2; 1/2[ 2 \] -1/6; 1/6[ 2 , see Fig. 6. The second and third families are nonconforming meshes used on the domain Ω =] -1/2; 1/2[ 2 . The second family has very localized non-conformities, see Fig. 7, and is obtained in the following way: the first mesh is obtained by dividing the domain into 8 × 8 identical squares, and the 4 squares at the center of the mesh are further refined into 4×4 sub-squares. Then, the subsequent meshes are obtained by dividing each cell of the previous mesh into 2 × 2 square cells. The third family has non-conformities spread over the entire domain, see Fig. 8, since every other cell is refined into 4 × 4 sub-cells. Of course, the third family of meshes is not of practical use but illustrates well the ability of the scheme to deal with heavily non-conforming meshes.

We have proved (see Sections 5.1 and 5.2) that p, ω, as well as their gradients converge to the exact solution of the Stokes problem. We are thus interested here in the numerical convergence of p, ∇ D h p, ω and ∇ D h ω, which we measure by the following errors

(e0p)(h) := ||p -Πp|| T,P ||Πp|| T,P and (e1p) 2 (h) := j |D j | |(∇ D h p) j -(Π∇p) j | 2 j |D j | |(Π∇p) j | 2
, where ∀j ∈ [1, J], (Π∇p) j = (∇p)(B j ), where B j is the center of gravity of the diamond cell D j . The same definitions hold for ω by replacing p by ω in the previous formulae. From the numerical results given in [START_REF] Domelevo | A finite volume method for the Laplace equation on almost arbitrary two-dimensional grids[END_REF], we may expect second-order accuracy for p and ω (although we were able to prove only first-order accuracy in section 5), and first-order accuracy for ∇ D h p and ∇ D h ω on general meshes. However, on meshes with diamond-cells which are almost all parallelograms, which is the case for the second family of meshes, we expect a convergence order of 1.5 for ∇ D h p and ∇ D h ω. We also proved the convergence of the velocity field u to the exact solution. The discrete relative L 2 error on the diamond cells for the velocity is measured by:

e 2 (h) := j |D j | |u j -(Πû) j | 2 j |D j | |(Πû) j | 2 ,
where (Πû) j is the value of the exact solution û at the midpoint M j of the edge A j : ∀j ∈ [1, J], (Πû) j = û(M j ).

We expect first-order convergence of the velocity field on general meshes like those of the first and third families. On the second family of meshes, since almost all diamond cells are parallelograms, we may expect from the numerical results of [START_REF] Delcourte | A discrete duality finite volume approach to Hodge decomposition and div-curl problems on almost arbitrary two-dimensional meshes[END_REF] an order of convergence of at least 1.5. In Finite Element methods, one is also usually concerned with the convergence of ∇u, since u belongs to H 1 (Ω) and since the term Ω ∇u : ∇vdx appears in the bilinear form associated with the variational formulation of the Stokes problem. In our formulation, the natural norm induced by the variational formulation is

||∇ • u|| 2 L 2 (Ω) + ||∇ × u|| 2 L 2 (Ω) 1/2
. Since ∇ T,P h

• u is always exactly imposed through the first equation of (4.15) or (4.16), we measure the errors on the derivatives of u through the error on ω = ∇ T,P h × u.

Triangular meshes

We first consider standard triangular meshes, as shown in Fig. 6. On this type of meshes, ∇ D h p and ∇ D h ω are all first-order accurate, while p, ω and u are second-order accurate, as displayed in Fig. 6. Apart from u which displays superconverge, these are the expected orders of convergence.

Locally refined meshes

On the second family of meshes (see Fig. 7), we observe a super-convergence of order 1.5 of ∇ D h p and ∇ D h ω, as expected. Moreover, as far as u is concerned, we observe in practice an order of convergence which is better than expected since it is slightly lower than 2. On the third family of meshes, see Fig. 8, we recover the same orders of convergence as those obtained on triangular meshes, except for u which has only first-order convergence.

The Bercovier-Engelman test case

The domain of computation is the unit square and the boundary conditions are u • n = 0 on the whole boundary and ω = 256y 2 (y -1) 2 for x = 0 or x = 1 and ω = 256x 2 (x -1) 2 for y = 0 or y = 1, while f = (f 1 , f 2 ) T with f 1 (x, y) = 256 x 2 (x -1) 2 (12y -6) + y(y -1)(2y -1)(12x 2 -12x + 2) + (y -1/2) f 2 (x, y) = -f 1 (y, x)

The exact solution is ω = 256 y 2 (y -1) 2 (6x 2 -6x + 1) + x 2 (x -1) 2 (6y 2 -6y + 1) , p = (x -1/2)(y -1/2) and u = (u 1 , u 2 ) T with u 1 (x, y) = -256y(y -1)(2y -1)x 2 (x -1) 2 u 2 (x, y) = -u 1 (y, x). This test case is interesting because in the source term f , the component related to ∇ × ω is much larger in magnitude as that related to ∇p, and thus one may wish to know if a numerical scheme is able to calculate the pressure field correctly. This is indeed the case as may be seen from the results displayed on Figure 9 which were computed on a triangular mesh with around 63000 triangles.

The Ruas test-case

The domain of computation is the unit circle, we use u • n = 0 on the boundary, where we also prescribe the value of the exact vorticity ω = 32 -16(x 2 + y 2 ). Setting f = ∇ × ω and g = 0 yields that p = 0. This test is interesting because a numerical method will compute a discrete pressure which will be, in general, different from 0, and thus one wishes to know the magnitude of this spurious pressure. However, in the DDFV framework, we wish to emphasize that there is a way of prescribing the right-hand side f j on each diamond-cell D j in the discrete Helmholtz-Hodge decomposition (Step 1 of our procedure) such that the discrete values p and ω computed by our scheme are exactly equal to 0 and Πω, respectively. Indeed, instead of using (4.10), we define f D j by its scalar products with n j and n ′ j in the following way and, in the same way

f D j • n j = 1 |A j | Aj f • n j dσ and f D j • n ′ j = 1 |A ′ j | A ′ j1 f • n ′ j1 dσ + A ′ j2 f • n ′ j2 dσ . ( 6
A ′ j1 f • n ′ j1 dσ + A ′ j2 f • n ′ j2 dσ = (ω(M j ) -ω(G i1 )) + (ω(G i2 ) -ω(M j )) = ω(G i2 ) -ω(G i1 ) = |A ′ j |(∇ D h × Πω) j • n ′ j ,
we have that f j = (∇ D h × Πω) j for all j; then it follows from the uniqueness of the discrete Helmholtz-Hodge decomposition of the discrete f that ω i,k = (Πω) i,k and p T,P i,k = 0. Of course, defining f D j by (6.2) requires the calculations of the 1D integrals along the edges; this may be performed by numerical quadrature, up to a user-defined precision. In the Ruas test case, since f is a P 1 function, the midpoint rule is enough to calculate the integrals exactly. We checked that the computed vorticity is exactly equal to Πω and that the computed pressure vanishes, up to the linear solver precision. The vorticity and the velocity fields are displayed on Figure 10.

Conclusion

We have proposed a discrete duality finite volume method for the two-dimensional Stokes equations with the non-standard boundary conditions (4.1) to (4.4). These non-standard boundary conditions are treated through the vorticity-velocity-pressure formulation of the Stokes equations, for which the finite volume method was successfully applied over unstructured and nonconforming meshes. Numerical results show a first-order convergence for the velocity, the pressure gradient and the vorticity gradient and a second order convergence for the pressure and the vorticity, while a superconvergence order of 1.5 for the pressure gradient and the vorticity gradient and an order two for the velocity are obtained on regular (but possibly locally non-conforming) meshes. Some of these convergence orders were proved through the theoretical analysis we led in this paper, while only suboptimal orders were obtained for pressure and vorticity on general meshes and for velocity on regular meshes.

  ) ∇ • u = g in Ω, (1.2) supplemented with one of the following non-standard sets of conditions u • n = σ over Γ , ∇ × u = ω d over Γ and Ω p(x)dx = 0 , (1.3) u • n = σ over Γ , p = p d over Γ and Ω ∇ × u(x)dx = m ω , (1.4) u • τ = σ over Γ , ∇ × u = ω d over Γ and Ω p(x)dx = 0 , (1.5) u • τ = σ over Γ , p = p d over Γ and Ω ∇ × u(x)dx = m ω . (1.6)

Figure 3 :

 3 Figure 3: Notations for the diamond cell.

3 αFigure 4 :

 34 Figure 4: Notations for Hypothesis 2.1. We can easily verify that Hypothesis 2.1 implies various weaker hypotheses: [15, Hypotheses 5.5 and 5.6], [25, Hypothesis 5.5] and [35, Definition 2.1] under which results obtained in these articles, and used in the present work, are valid.

  ) in the following way i∈[1,I]

2 k∈Γq 2 k∈Γq

 22 (4.7) with ∇ D h × ψ yields (4.13) and (4.14) again. Now, we use the fact that the first term in the left-hand side of (4.14) equals -1 |P k | (∇ P h • ∇ D h ω) k while the first term in the right-hand side of (4.14) equals 1 |P k | (∇ P h × s) k . Using the previously proved third line of (4.11), we obtain the fourth line of (4.11).

. 18 )

 18 Moreover, the right-hand sides in (4.15) and (4.16) have to verify compatibility conditions. Indeed, it is readily seen from the definition of the discrete divergence and curl operators (3.3) and (3.4) that the following two equalities hold i∈[1,I]

Lemma 4 . 4 .

 44 It was shown in[START_REF] Delcourte | A discrete duality finite volume approach to Hodge decomposition and div-curl problems on almost arbitrary two-dimensional meshes[END_REF] Proposition 5.2] that systems of the type(4.23) and(4.24) both have a unique solution.

Definition 5 . 5 .

 55 H 1 semi-norm for the pressure and the vorticity For any continuous function v, we define Πv by ∀i ∈ [1,

. 7 )

 7 Settingε p := p -Πp and ε ω := ω -Πω, (5.8)Equ. (5.7) implies the following equality:

Proposition 5 . 13 .

 513 [START_REF] Eymard | A staggered finite volume scheme on general meshes for the Navier-Stokes equations in two space dimensions[END_REF] is given by[START_REF] Delcourte | A discrete duality finite volume approach to Hodge decomposition and div-curl problems on almost arbitrary two-dimensional meshes[END_REF] Th. 5.22]: If all diamond-cells are convex and under Hypotheses 2.1 and 5.4, there exists a constant C(θ * ) independent of h such that ||ũ h -û|| 0,Ω ≤ C(θ * )h ||g|| 0,Ω + ||ω|| 0,Ω + || φ|| 2,Ω + || ψ|| 2,Ω .

Figure 6 :

 6 Figure 6: Left: Triangular mesh. Right: Errors.

Figure 7 :

 7 Figure 7: Left: Non-conforming mesh. Right: Errors.

Figure 8 :

 8 Figure 8: Left: Non-conforming mesh. Right: Errors.

Figure 9 :

 9 Figure 9: Bercovier-Engelman test. Left: pressure contours (expected extrema : -1/4 and 1/4, computed extrema: -0.24989 and 0.25003). Center: velocity field. Right: vorticity contours (expected extrema : -16 and 16, computed extrema: -16.001 and 16).

  ω) • n j = Aj (∇ω) • τ j = ω(S k1 )ω(S k2 ) = |A j |(∇ D h × Πω) j • n j

Figure 10 :

 10 Figure 10: Ruas test. Left: velocity field. Right: vorticity contours.

  ) to non-simply connected domains, with σ, p d and ω d given functions, m ω a given real number and (k q ) q∈[1,Q] a set of given real numbers:

	(1.6
	.22)
	4 Application to the Stokes equations
	In this section, we are interested in the discretization of Eqs. (1.7), (1.8) and (1.2) supplemented
	with one of the following non-standard sets of conditions, which generalize conditions (1.3) to

  to(3.17). A similar result holds for problem (4.16) using the Hodge decomposition with boundary conditions(3.19) to(3.21).

	Proposition 4.3. Problem (4.15) may be split into two independent problems:
	find

)

  Proposition 5.14. There exists a constant C(θ * ) independent of h such that ||u -ũ|| D ≤ C(θ * )||ωω || T,P . (5.35) Proof. We use the discrete Hodge decomposition (see Proposition 3.8) of d u into ∇ D h φ d + ∇ D h × ψ d . Using Proposition 4.3, we get that φ d vanishes. Next, performing like in the proof of Theorem 4.1 in [35], we obtain ||d u || 2 D = (∇ T,P h × d u , ψ d ) T,P . By definition of the scalar product (3.6), and using the second line in (5.34), it holds that

	(∇ T,P h	× d u , ψ d ) T,P =	1 2	i∈[1,I]

  × d u , ψ d ) T,P = (d ω , ψ d ) T,P . (5.39) Using (5.39), the Cauchy-Schwarz inequality (d ω , ψ d ) T,P ≤ ||d ω || T,P ||ψ d || T,P and the discrete Poincaré inequality (3.22) for ψ d , allow us to conclude. Next, it remains to estimate d ω . We first link ||ωω || T,P to a previous result by the triangle inequality ||ωω || T,P ≤ ||ω -Πω|| T,P + || ω -Πω|| T,P .

				T,P h	× d u ) P k (ψ d ) P k .	(5.36)
	Using boundary condition (3.16) for ψ d , we obtain
	|P k |(∇ T,P h	× d u ) P k (ψ d ) P k = 0 =	|P k |(d ω ) P k (ψ d ) P k .	(5.37)
	k∈Γ0		k∈Γ0
	Moreover, using boundary condition (3.17) and the last line in (5.34), we obtain
	|P k |(∇ T,P h	× d u ) P k (ψ d ) P k =	|P k |(d ω ) P k (ψ d ) P k .	(5.38)
	k∈Γq		k∈Γq
	Thus, gathering (5.36), (5.37), (5.38), and recalling that ||∇ D h × ψ d || 2 D = ||d u || 2 D = (∇ T,P h	×
	d u , ψ d ) T,P , we get
	||∇ D h × ψ d || 2 = (∇ T,P h
					(5.40)

  the following proposition: ) independent of h such that||u h -û|| 0,Ω ≤ C(θ * )h ω 2,Ω + p 2,Ω + ||g|| 0,Ω + || φ|| 2,Ω + || ψ|| 2,Ω .Remark 5.18. The convergence results obtained in this section are valid on any mesh; however, on some type of meshes, it is possible to prove some superconverge results. They will not be proved here, but it may be shown that if ω and p given by Proposition 5.1 belong to H 3 (Ω), then, on families of meshes such that there exists a finite number L of sub-domains (Ω ℓ ) ℓ∈[1,L] included in Ω and independent of the mesh step h, such that the diamond-cells which are not parallelograms are included in strips δ ℓ , having width Ch and located along the boundaries of the various Ω ℓ , where C is a constant independent of h, then ||∇ D h (p-Πp)|| D and ||∇ D h (ω-Πω)|| D both converge to 0 with the order h 3/2 , and, using the discrete Poincaré inequality, so do ||p -Πp)|| T,P and ||ω -Πω)|| T,P .

	Proposition 5.16. Under hypotheses 2.1 and 5.2, there exists a constant C(θ * ) independent
	of h such that
	||u -ũ|| D ≤ C(θ * )h ( ω 2,Ω + p 2,Ω ) .	(5.47)
	Finally, plugging (5.47) and (5.33) into (5.32) leads to the following theorem
	Theorem 5.17. If all diamond cells are convex and under hypotheses 2.1, 5.2 and 5.4, there
	exists a constant C(θ (5.48)

*

Partially supported by the European Research Council ERC Starting Grant 2009, project 239983-