Ferréol Soulez
email: ferreol.soulez@univ-lyon1.fr

Loïc Denis

Yves Tourneur

Éric Thiébaut

BLIND DECONVOLUTION OF 3D DATA IN WIDE FIELD FLUORESCENCE MICROSCOPY

Keywords: blind deconvolution, wide field fluorescence microscopy, image restoration, inverses problems

In this paper we propose a blind deconvolution algorithm for wide field fluorescence microscopy. The 3D PSF is modeled after a parametrized pupil function. The PSF parameters are estimated jointly with the object in a maximum a posteriori framework. We illustrate the performances of our algorithm on experimental data and show significant resolution improvement notably along the depth. Quantitative measurements on images of calibration beads demonstrate the benefits of blind deconvolution both in terms of contrast and resolution compared to non-blind deconvolution using a theoretical PSF.

INTRODUCTION

Using fluorescent dyes to identify specific cellular structures, the wide field fluorescence microscopy (WFFM) is a widely spread imaging modality in biology. It consists on imaging at its emission wavelength a cellular structure marked by fluorescent dye excited by uniform illumination. On the resulting 2D image, structures are more or less defocalized according to their distance to the focal plane. Moving this focal plane through the sample produces a 3D representation of the object. WFFM however suffers from a very coarse axial resolution compared to the lateral resolution. This blurring of the data y can be modeled by a 3D convolution between a point spread function (PSF) h and the observed 3D object x:

y = h * x + n , (1)
where n accounts for the noise. The resolution of the setup is given by the shape of the 3D PSF. To improve depth resolution, one can either change the imaging setup (e.g., using confocal or two-photon microscopy) or enhance WFFM 3D data with a deconvolution algorithm ( [START_REF] Sibarita | Deconvolution microscopy[END_REF][START_REF] Sarder | Deconvolution methods for 3-d fluorescence microscopy images[END_REF] and references therein). Most deconvolution algorithms use the theoretical diffraction-limited PSF or a more realistic experimental PSF measured using calibration beads. Blind deconvolution algorithms bypass the problem of PSF calibration by simultaneously estimating the PSF and the object. Few blind deconvolution algorithms have been proposed for WFFM, the most notable being parametric blind deconvolution (PBD) [START_REF] Markham | Parametric blind deconvolution: a robust method for the simultaneous estimation of image and blur[END_REF], AIDA [START_REF] Hom | Aida: an adaptive image deconvolution algorithm with application to multi-frame and three-dimensional data[END_REF][START_REF] Mugnier | Mistral: a myopic edgepreserving image restoration method, with application to astronomical adaptive-optics-corrected long-exposure images[END_REF] and the method proposed by Kenig[6] where the PSF is constrained to lie in a (previously learned) subspace. The proposed method is closely related to PBD [START_REF] Markham | Parametric blind deconvolution: a robust method for the simultaneous estimation of image and blur[END_REF] method as we model a PSF as a function of a parametrized pupil function. The main differences are (i) the use of a much more robust MAP framework instead of maximum likelihood and (ii) the parametrization of the PSF that enforces its normalization and can model nonsymmetric PSF.

MAP FRAMEWORK

In blind deconvolution both object x and PSF h are derived from the same measurements y. We consider maximum a posteriori estimates {x + , h + } obtained by minimizing cost function f (x, h):

{x + , h + } = arg min x∈X,h∈H J (x, h) , (2)
J (x, h) = J data (x, h) + µ J prior (x) , (3)
with J data (x, h) the neg-loglikelihood, J prior (x) a regularizing prior and µ a hyper-parameter. Prior knowledge on the PSF is implicitly enforced through PSF parametrization. Regularization J prior is thus applied only to object x.

Data-fitting term

Under Gaussian noise assumption, the neg-loglikelihood reads:

J data (x, h) = (y -h * x) ⊤ • C -1 noise • (y -h * x) , (4)
where Cnoise is the noise covariance matrix and is diagonal if noise is uncorrelated:

J data (x, h) = k w k (h * x) k -y k 2 , (5)
where w k is the inverse of the noise variance at pixel k. This model can cope with non-stationary noise and can be used to express confidence on measurements on each pixel of the data. Thus it can deal with unmeasured pixels (due to saturation, borders. . . ) by setting w k = 0 for such pixels. Furthermore, except for very low detector noise, this formulation can account for mixed Poisson + Gaussian noise by approximating it as a non-stationary uncorrelated Gaussian noise [START_REF] Mugnier | Mistral: a myopic edgepreserving image restoration method, with application to astronomical adaptive-optics-corrected long-exposure images[END_REF]:

w k def = γ max(y k , 0) + σ 2 k -1 if y k is measured, 0 otherwise, ( 6
)
where γ is the quantization factor of the detector and σ 2 k is the variance of other approximately Gaussian noises (e.g., read-out noise) at pixel k.

Object regularization term

As most observed objects are smooth with few sharp structures (e.g., edges and spikes), we use as regularizing prior a hyperbolic version of the classical 3D total variation [START_REF] Rudin | Nonlinear total variation based noise removal algorithms[END_REF]:

J prior (x) = k ∇x k 2 2 + ǫ 2 . ( 7
)
Parameter ǫ > 0 ensures differentiability of J prior at 0. When ǫ value is close to the quantization level, this function smooths out nonsignificant differences between adjacent pixels. In addition, since only intensities are measured, we constrain object x to be positive valued.

PSF PARAMETRIZATION

An abundant literature discusses PSF modeling for fluorescence microscopy. In the present work, we chose a monochromatic scalar model that defines the 3D PSF h from pupil function p. This pupil function is the in-focus point source wavefront phase and modulus at the exit pupil of the objective. To reduce the number of degrees of freedom of our model, the pupil function is parametrized by a limited number of modes. 3D PSF h(x, y, z) is defined as the squared magnitude of complex-valued amplitude PSF h(x, y, z):

h(x, y, z) = h(x, y, z) 2 . (8)
Complex-valued PSF h is, in turn, defined as the 3D Fourier transform of the complex-valued amplitude optical transfer function (OTF). This OTF is non-zero only on a spherical cap of radius n i /λ limited by the aperture angle. It can thus be fully described by a 2D function: pupil function p(κx, κy) which is non-zero on the disk defined by κ 2 1 + κ 2 2 ≤ NA/λ. 3D PSF h(x, y, z) is then related to pupil function p through:

h(x, y, z) = p(κx, κy) exp (2 i π z d(κx, κy)) . . . exp (2 i π (κx x + κy y)) dκxdκy 2 , (9)
with d(κx, κy) = (n i /λ) 2 -(κx + κy) 2 the defocus [START_REF] Stokseth | Properties of a defocused optical system[END_REF], n i the refractive index of the immersion medium and λ the emission wavelength of the fluorophore. The defocus and the 2D complex pupil function p(κx, κy) describe the properties of the optical system. As noted by Hanser [START_REF] Hanser | Phaseretrieved pupil functions in wide-field fluorescence microscopy[END_REF], the support of p(κx, κy) is a disk thus Zernike polynomials Zn provide a suitable basis to express both modulus ρ(κx, κy) and phase φ(κx, κy) of pupil function p: ρ(κx, κy) = n βnZn(κx, κy) , [START_REF] Besag | On the statistical analysis of dirty pictures[END_REF] φ(κx, κy) = n αnZn(κx, κy) .

As phase piston has no effect on the resulting PSF we set α0 = 0. Similarly, we cancel the tip-tilt parameters (α1 = α2 = 0) to keep the PSF laterally centered (this fixes the position degeneracy common with blind deconvolution problems). PSF normalization (i.e., energy conservation) can be expressed as an ℓ1 norm constraint h 1 = 1, or equivalently as an ℓ2 constraint on the complex-valued PSF h 2 = 1. This latter ℓ2 constraint can further be written as an ℓ2 constraint in Zernike orthonormal basis: β 2 = 1. By selecting only purely radial Zernike polynomials, one can easily constrain the axial symmetry of the PSF. This relatively simple PSF model offers a flexible parametrization that requires only the knowledge of the wavelength λ, the numerical aperture NA and the refractive index of the immersion medium n i .

BLIND DECONVOLUTION ALGORITHM

To solve the blind deconvolution, we have to find the optimal parameters {x + , α + , β + } that minimizes the cost function J (x, h(α, β)).

A simple way to do this is to use continuous optimization techniques. Due to non-linearity of this cost function and as these parameters have different physical meanings, the simultaneous estimation of parameters of both the object x and the PSF h(α, β) is known to very badly conditioned. This slows down convergence of optimization algorithm. For that reason, as it is classically done [START_REF] Besag | On the statistical analysis of dirty pictures[END_REF], we use an alternating minimization scheme to minimize the criterion J (x, h(α, β)). The algorithm is thus:

1. t = 0, define the initial PSF h (0) as the theoretical diffraction limited PSF of the microscope i.e., α (0) = 0, β (0) = 0.

2. t = t + 1, estimation of the optimal non-negative object x (t) according to given PSF h (t-1) :

x (t) ≈ arg min x≥0 J data (x, h (t-1) ) + µ J prior (x) . (12)
3. Improvement of the phase parameters α (t) given the other model parameters that minimize:

α (t) = arg min α J data x (t) , h(α, β (t-1) ) . (13)
4. Improvement of the modulus parameters β (t) given the other model parameters that minimize:

β (t) ≈ arg min β 2 2 =1
J data x (t) , h(α (t) , β) ,

5. Go to step 2 until convergence.

Each step involves minimization of a criterion with respect to a large number of variables. To that end, we used VMLM-B algorithm [START_REF] Thiébaut | Optimization issues in blind deconvolution algorithms[END_REF] which can account for bound constraints on the parameters to enforce object positivity. This algorithm has proved its effectiveness for image reconstruction and only requires the computation of the cost function and its gradient. The memory requirement is a few times the size of the problem. In each step, we do not exactly solve the minimization defined in Eq. ( 12), Eq. ( 13) and Eq. ( 14) but only execute few iterations (about 10) of each inner optimization.

RESULTS

We processed the dataset 1 used by Griffa et al. in [START_REF] Griffa | Comparison of deconvolution software in 3d microscopy. a user point of view part 1[END_REF][START_REF] Griffa | Comparison of deconvolution software in 3d microscopy. a user point of view part 2[END_REF] for comparative study of available deconvolution softwares. That study only considers "non-blind" deconvolution with theoretical PSF and we used it to evaluate the advantages of using a blind deconvolution technique compared to state of the art available "non blind" deconvolution methods. 1. Performance of 3 state-of-the-art deconvolution methods as reported by Griffa [START_REF] Griffa | Comparison of deconvolution software in 3d microscopy. a user point of view part 2[END_REF] compared to the proposed method (both blind and non-blind). Hyugens and AutoDeblur are commercial softwares and Deconvolution Lab is an imageJ plugin implementing [START_REF] Vonesch | A fast thresholded Landweber algorithm for wavelet-regularized multidimensional deconvolution[END_REF].

Calibration bead

The first dataset is an observation of an InSpeck green fluorescent bead with a known diameter (2.5µm) on a Olympus Cell R microscope with a ×63, 1.4 NA oil objective. The data cube is composed of 256 × 256 × 128 voxels of size 64.5 × 64.5 × 160 nm 3 . During the first step of the algorithm, the solution of Eq. ( 12) is a "nonblind" deconvolution with a theoretical PSF (no aberration). On this solution, plotted in Fig. 1 Although the in-focus axial sections of both PSF are similar, there are strong differences between axial sections. Unlike the theoretical PSF, the estimated PSF is not symmetric with respect to the plane (z = 0). This is well explained by a mismatch between the refractive indices of the immersion oil and the mounting medium.

To perform quantitative evaluation, Griffa et al. [START_REF] Griffa | Comparison of deconvolution software in 3d microscopy. a user point of view part 2[END_REF] proposed three criteria: the size of the bead given by the lateral and the axial full width at half maximum (FWHM) and the relative contrast between the center of the bead and the maximum of its radial profile. We present in Tab. 1, the values of these criteria given by [START_REF] Griffa | Comparison of deconvolution software in 3d microscopy. a user point of view part 2[END_REF] together with the values obtained with our method. We can notice that in its non-blind form (i.e. step 2 only), the proposed 3D-TV is competitive with the state of the art. As it tends to smooth the object, TV regularization may lead to an over-estimation of the diameter of the bead (its nominal width is 2500 nm). However the axial FWHM and the relative contrast shows the improvements given by our blind deconvolution method.

C Elegans embryo

In addition to this quantitative evaluation on a calibration bead, we processeded the biological sample used by Griffa et al. in their study. It is an observation of a C. Elegans embryo containing DAPI, FITC and CY3 stainings with the same microscope and objective than in previous section. The data cube is composed of 672 × 712 × 104 voxels of size 64.5 × 64.5 × 200 nm 3 in three spectral channels. The DAPI (blue 477nm) stains chromosomes in the nuclei, the FITC (green 542nm) the microtubule filaments and CY3 (red 654nm) some point wise protein aggregates. We processed each channel from this dataset individually using the same hyperparameter for each channel (µ = 2 × 10 -5 ). For each channel, the algorithm takes about 4 hours to converge using an Intel i7-975 CPU. The result of the blind deconvolution (Fig. 3) presents a significantly improved contrast without the typical haze of fluorescence wide field microscopy: the red spots are much sharper and brighter in the restoration. Furthermore, the resolution is clearly improved especially along the depth axis: the out of focus nuclei encircled in red is almost invisible in the blind deconvolution. However, as it can be seen on axial sections, the deconvolution becomes less effective as one looks deeper in the sample. This can be explained by the use of a single stationary PSF for the whole observation. Because of refractive index mismatch, it is well known that the real PSF may vary axially and even laterally.

CONCLUSION

We presented a new method for blind deconvolution of WFFM data. It uses a PSF parametrization by means of decomposition of the pupil on Zernike basis. We used a continuous optimization method to iteratively estimate both PSF parameters and the object. Except from some parameters of the setup that are generally known (wavelength, numerical aperture and refractive index of the immersion medium), it requires only tuning one hyper-parameter. Results show clear improvements of both resolution and contrast.

Possible directions for future research include extension of the PSF model to more sophisticated models (e.g. vectorial model with pixel integration), extending this method to confocal and two photon microscopy, and introducing space-variant PSF blinddeconvolution [START_REF] Denis | Fast model of space-variant blurring and its application to deconvolution in astronomy[END_REF].

Fig. 1 .

 1 Fig. 1. Sections of observed and restored images of an InSpeck green fluorescent bead of diameter 2.5 µm. Axis are in µm

Fig. 2 .

 2 Fig. 2. Sections of theoretical ((a) and (b)) and estimated PSF ((c), (d) and (e)). Axis are in µm.

Fig. 3 .

 3 Fig. 3. C. Elegans embryo observation and the blind deconvolution result. Each panels are cuts in the 3D volume along the yellow lines.

Table

	data	Hyugens AutoDeblur Deconvolution	proposed method
	parameters			Lab	non-blind blind
	transversal FWHM (in nm) 2867	2709	2709	2664	2736	2783
	axial FWHM 4760	4000	4640	4160	3054	2977
	Relative contrast 18 %	53 %	78 %	68 %	84 %	88 %

the dataset can be found at http://bigwww.epfl.ch/ deconvolution

This work has been supported by project MiTiV funded by the French National Research Agency (ANR DEFI 09-EMER-008-01).


O Angelini

C Chavant

E Chénier

R Eymard

A finite volume scheme for diffusion problems on general meshes applying monotony constraints

In order to increase the accuracy and the stability of a scheme dedicated to the approximation of diffusion operators on any type of grids, we propose a method which reduces the curvature of the discrete solution where the loss of monotony is observed. The discrete solution is shown to fulfill a variational formulation thanks to the use of Lagrange multipliers. We can then show its convergence to the solution of the continuous problem, and an error estimate is derived. A numerical method, based on Uzawa's algorithm, is shown to provide accurate and stable approximate solutions to various problems. Numerical results show the increase of precision due to the application of the method.

The initial scheme

We now describe in details the scheme, first proposed in [8] and extended in [9], used to provide the properties that are required to follow the lines given in the introduction. To this purpose, we first define the geometrical elements entering into the definition of a discretization. Definition 2.1 (Discretization) Let Ω be a polyhedral open bounded connected subset of R d , with d ∈ N ⋆ , and ∂Ω = Ω\Ω its boundary. A discretization of Ω, denoted by D, is defined as the triplet D = (M, E, P), where:

• M is a finite family of non empty connected open disjoint subsets of Ω ("'the control volumes"') such that Ω = K∈M K. For any K ∈ M , let ∂K = K\K be the boundary of K ; m K > 0 denotes the measure of K and h K denotes the diameter of K.

Introduction

A few pioneering works address the problem of monotony of numerical schemes for handling anisotropic diffusion operators on general grids. Let us quote [15; 16; 13; 14; 17; 6; 3], who attempted various methods in order to circumvent this problem. Recall that all of these methods are nonlinear, and rely on coefficients computed from the solution itself. Among them, let us recall a natural one, which consists in writing the equations associated with each unknown in a way such that each unknown can be expressed as a convex combination of neighboring ones; this is achieved in finite volume methods, writing adapted two-point expressions for the fluxes at the edges of a control volume, and ensuring the conservation of the fluxes. One difficulty which must be overcome using such methods is to achieve the convergence of the fixpoint method, necessary for solving these nonlinear schemes (see [START_REF] Lipnikov | Interpolation-free monotone finite volume method for diffusion equations on polygonal meshes[END_REF] for a discussion on this problem). We have therefore investigated another way of increasing the numerical stability of numerical schemes, using new nonlinear scheme, which modifies the curvature of the discrete solution. In order to describe the method, let us first give the classical Dirichlet problem, with a heterogeneous and anisotropic diffusion matrix. Our aim is to provide an approximation to the following problem:

-div Λ(x)∇u = f in Ω u = 0 on ∂Ω (1)
with the following hypotheses on the data:

Ω is a polyhedral open bounded connected subset of R d with d ∈ N ⋆ , Λ is a measurable function from Ω to M d (R),
There exist λ and λ such that 0 < λ ≤ λ and Sp(Λ(x)) ⊂ [λ, λ] for a.e. x ∈ Ω,

f ∈ L 2 (Ω), (2)
where we denote by M d (R) the set of symmetric positive definite d × d real matrices, and, for A ∈ M d (R), Sp(A) is the set of the eigenvalues of A. Under these hypotheses, the weak solution of (1) is the unique function u satisfying:

   u ∈ H 1 0 (Ω), Ω Λ(x)∇u(x).∇v(x)dx = Ω f (x)v(x)dx, ∀v ∈ H 1 0 (Ω) (3)
The approach that we propose in this paper is the following. We first consider the use of the hybrid finite volume scheme (see [START_REF] Eymard | A new finite volume scheme for anisotropic diffusion problems on general grids: convergence analysis[END_REF] or [START_REF] Eymard | Discretisation of heterogeneous and anisotropic diffusion problems on general non-conforming meshes. SUSHI: a scheme using stabilisation and hybrid interfaces[END_REF]) for finding an approximate solution to Problem (3). This scheme, which can be related to the mimetic finite difference method [4; 5], has the property to degenerate on the well-known two-point flux finite volume method with harmonic averaging of the diffusion [START_REF] Herbin | An error estimate for a finite volume scheme for a diffusion-convection problem on a triangular mesh[END_REF] in the particular case of triangles or rectangles and isotropic diffusion, which has the advantage to be cheap and easily implementable. In such a case, accounting from the above discussion, monotony properties are fully satisfied. The hybrid finite volume method consists in the general case in finding an element u ∈ V which satisfies the minimum value of a real α-elliptic function J on some finite dimensional Euclidean space V . The main mathematical results, concerning this scheme, are recalled in Section 2.

We then introduce the subset K ⊂ V of all v ∈ V such that G i (v) ≤ 0 for i = 1, . . . , p. These functions G i are convex regular functions, such that the interpolation in V of any smooth function belongs to K, at least for sufficiently fine discretizations. The new problem to solve is then to find u ∈ K which satisfies the minimum on K of the function J.

We show in Section 3 that this problem is well-posed, and that the discrete solution thus provided converges to the exact solution. An error estimate is also proposed. We finally show, on a numerical example proposed in Section 4, the efficiency of the method (the discrete solution being approximated thanks to Uzawa's algorithm, which has the advantage in this framework to be easily implemented in a code dedicated to the computation of the non-constrained scheme). We recall in an appendix the classical theorems concerning the existence of Lagrange multipliers and the convergence of Uzawa's algorithm.

• E is a finite family of disjoint subsets of Ω (the "'edges"' of the mesh), such that, for all σ ∈ E , σ is a non empty closed subset of a hyperplane of R d , which has a measure m σ > 0 for the (d -1) dimensional measure of σ.

We assume that, for all K ∈ M, there exists a subset E K of E such that ∂K = σ∈E K σ. We then denote by

M σ = {K ∈ M, σ ∈ E K }.
We then assume that, for all σ ∈ E, either M σ has exactly one element and then σ ⊂ ∂Ω (the set of these interfaces, called boundary interfaces, is denoted by E ext ) or M σ has exactly two elements (the set of these interfaces, called interior interfaces, is denoted by E int ). For all σ ∈ E, we denote by x σ the barycenter of σ. For all K ∈ M and σ ∈ E K , we denote by n K,σ the unit vector normal to σ outward to K.

• P is a family of points of Ω indexed by M, denoted by P = (x K ) K∈M , such that for all K ∈ M, x K ∈ K and K is assumed to be x k -star-shaped, which means that for all x ∈ K, the property [x K , x] ⊂ K holds. Denoting by d K,σ the Euclidean distance between x K and the hyperplane including σ, one assumes that d K,σ > 0.

The following notations are used. The size of the discretization is defined by: h D = sup {h K , K ∈ M}. The regularity of the mesh is measured through the parameter

θ D = max max σ∈E int,K,L∈Mσ d K,σ d L,σ , max σ∈E K,K∈Mσ h K d K,σ (4)
We have, thanks to the assumption that K is x K -star-shaped, the property:

σ∈E K m σ d K,σ = d m K (5) Let X D = {v = ((v K ) K∈M , (v σ ) σ∈E ), v K ∈ R, v σ ∈ R} = R M × R E be
the discrete space into which we look for the solution of the schemes. Let X D,0 ⊂ X D be defined by X D,0 = {u ∈ X D , u σ = 0, σ ∈ E ext }, taking into account the homogeneous Dirichlet boundary conditions.

For all φ ∈ C 0 (Ω), we denote by P D φ the element of X D defined by (φ(x K )) K∈M , (φ(x σ )) σ∈E . The space X D,0 is equipped with a Euclidean structure, defined by the following inner product:

∀(u, v) ∈ (X D,0 ) 2 , [u, v] D = K∈M σ∈E K m σ d K,σ (u σ -u K )(v σ -v K ), (6)
and the associated norm:

u 1,D = ([u, u] D ) 1 2
. Let H M (Ω) ⊂ L 2 (Ω) be the set of piecewise constant functions on the control volumes on the mesh M. For any u ∈ X D , we define

∇ K u = 1 m K σ∈E K m σ (u σ -u K )n K,σ , ∀K ∈ M, (7)
and we denote by ∇ D u ∈ H M (Ω) d the function defined by

∇ D u(x) = ∇ K u, for a.e. x ∈ K, ∀K ∈ M. (8)
Note that, thanks to the Cauchy-Schwarz inequality, we get

∇ D u L 2 (Ω) d ≤ d m Ω u 1,D . (9)
We also denote by P M u ∈ H M (Ω) the function such that P M u(x) = u K , for a.e. x ∈ K and all K ∈ M.

For a given family of strictly positive reals α = (α K ) K∈M , we consider the bilinear form defined by:

u, v D,α = K∈M m K ∇ K u • Λ K ∇ K v + α K σ∈E K m σ d K,σ R K,σ (u)R K,σ (v)n K,σ • Λ K n K,σ (10)
where we define

R K,σ (u) = u σ -u K -∇ K u • (x σ -x K ) d K,σ , ∀K ∈ M, ∀σ ∈ E K , (11)
and

Λ K = 1 m K K Λ(x)dx, ∀K ∈ M. ( 12
)
The scheme for the approximation of Problem (3) consists in finding u ∈ X D,0 such that

u, v D,α = Ω f (x)P M v(x)dx, ∀v ∈ X D,0 . (13)
We then have the following results, proven in [START_REF] Eymard | A new finite volume scheme for anisotropic diffusion problems on general grids: convergence analysis[END_REF].

Lemma 2.1 (Discrete Poincaré inequality)

Let us assume Hypothesis [START_REF] Andreianov | Discrete duality finite volume schemes for Leray-Lions-type elliptic problems on general 2D meshes[END_REF]. Let D be a discretization of Ω in the sense of Definition 2.1 and let θ ≥ θ D be given, where θ D is defined by [START_REF] Brezzi | Convergence of the mimetic finite difference method for diffusion problems on polyhedral meshes[END_REF]. Then there exists C 1 , only depending on d, Ω and on θ, such that

P M u L 2 (Ω) ≤ C 1 u 1,D . ( 14
) Lemma 2.2 (Coerciveness of •, • D,α )
Let us assume Hypothesis [START_REF] Andreianov | Discrete duality finite volume schemes for Leray-Lions-type elliptic problems on general 2D meshes[END_REF]. Let D be a discretization of Ω in the sense of Definition 2.1 and let θ ≥ θ D be given, where θ D is defined by [START_REF] Brezzi | Convergence of the mimetic finite difference method for diffusion problems on polyhedral meshes[END_REF]. Let α = (α K ) K∈M α be a family of strictly positive reals and let α ∈ (0, min{α K , K ∈ M}).

Then there exists α 0 > 0, only depending on d, α, λ and on θ, such that

α 0 [u, u] D ≤ u, u D,α , ∀u ∈ X D,0 . (15)
The two above Lemmas provide sufficient conditions to ensure the existence and uniqueness of the solution to [START_REF] Lipnikov | Monotone finite volume schemes for diffusion equations on unstructured triangular and shape-regular polygonal meshes[END_REF].

Thanks to the symmetry of •, • D,α , we get that this solution is also given by

u = argmin v∈X D,0 J D,α (v), (16)
where the notation argmin x∈E f (x) denotes an element y ∈ E, assumed to exist and to be unique, such that the minimum value of the function f on E is equal to f (y), and where J D,α is defined by

J D,α (v) = 1 2 v, v D,α - Ω f (x)P M v(x)dx, ∀v ∈ X D,0 . (17)
Lemma 2.3 (Relative compactness in X D (m) ,0 ) Let us assume Hypothesis [START_REF] Andreianov | Discrete duality finite volume schemes for Leray-Lions-type elliptic problems on general 2D meshes[END_REF]. Let (D (m) ) m∈N be a family of discretizations of Ω in the sense of Definition 2.1, such that there exists θ ≥ 0 with θ ≥ θ D (m) for all m ∈ N. Let (u m ) m∈N be a sequence such that there exists C 2 > 0 such that, for all m ∈ N, u m ∈ X D (m) ,0 and u m 1,D (m) ≤ C 2 . Then there exists u ∈ H 1 0 (Ω) and a sub-sequence of (u m ) m∈N , again denoted by (u m ) m∈N , which converges to u in L 2 (Ω). Moreover, (∇ D (m) u m ) m∈N weakly converges to ∇u in L 2 (Ω) d . Lemma 2.4 (Consistency of ∇ D u and R K,σ u) Let us assume Hypothesis [START_REF] Andreianov | Discrete duality finite volume schemes for Leray-Lions-type elliptic problems on general 2D meshes[END_REF]. Let D be a discretization of Ω in the sense of Definition 2.1 and let θ ≥ θ D be given, where θ D is defined by [START_REF] Brezzi | Convergence of the mimetic finite difference method for diffusion problems on polyhedral meshes[END_REF]. Let φ ∈ C 2 (Ω). Then there exists C 3 , only depending on d, θ and φ (through its second order partial derivatives) such that

|∇ K P D φ -∇φ(x K )| ≤ C 3 h D ,
where we denote by |x| the Euclidean norm of any x ∈ R d , and

|R K,σ (P D φ)| ≤ C 3 h D .
The following lemma is a consequence of the previous ones.

Lemma 2.5 (Consistency of •, • D,α )
Let us assume Hypothesis [START_REF] Andreianov | Discrete duality finite volume schemes for Leray-Lions-type elliptic problems on general 2D meshes[END_REF]. Let (D (m) ) m∈N be a family of discretizations of Ω in the sense of Definition 2.1, such that there exists θ ≥ 0 with θ ≥ θ D (m) for all m ∈ N. Let (u m ) m∈N be a sequence such that there exists C 2 > 0 such that, for all m ∈ N, u m ∈ X D (m) ,0 and u m 1,D (m) ≤ C 2 . We assume that there exist α ≥ α > 0, and, for all m ∈ N, a family α (m) = (α

(m) K ) K∈M (m) , such that {α (m) K , K ∈ M} ⊂ [α, α].
We also assume that there exists u ∈ H 1 0 (Ω) such that (u m ) m∈N converges to u in L 2 (Ω). Then the following holds

lim m→∞ u m , P D (m) φ D (m) ,α (m) = Ω Λ(x)∇u(x).∇φ(x)dx, ∀φ ∈ C ∞ c (Ω), (18)
where C ∞ c (Ω) is the set of elements of C ∞ (Ω) with compact support. We can finally state the convergence result, which can be easily proven using the previous lemmas (in the next section, we state and prove a convergence result, following similar steps).

Theorem 2.1 (Convergence of the non-constrained scheme) Let us assume Hypothesis [START_REF] Andreianov | Discrete duality finite volume schemes for Leray-Lions-type elliptic problems on general 2D meshes[END_REF]. Let θ > 0 and α ≥ α > 0 be given. Let D be a discretization of Ω in the sense of Definition 2.1 such that θ ≥ θ D , where θ D is defined by [START_REF] Brezzi | Convergence of the mimetic finite difference method for diffusion problems on polyhedral meshes[END_REF]. Let α = (α K ) K∈M be a family of reals such that {α K , K ∈ M} ⊂ [α, α]. Let u D be the unique solution of [START_REF] Lipnikov | Monotone finite volume schemes for diffusion equations on unstructured triangular and shape-regular polygonal meshes[END_REF]. Then u D tends in L 2 (Ω) to u, the unique solution of (3), as h D tends to 0. Moreover, ∇ D u D tends in L 2 (Ω) d to ∇u.

The scheme with constraints

Let us again assume Hypothesis [START_REF] Andreianov | Discrete duality finite volume schemes for Leray-Lions-type elliptic problems on general 2D meshes[END_REF]. Let D be a discretization of Ω in the sense of Definition 2.1 and let θ ≥ θ D be given, where θ D is defined by [START_REF] Brezzi | Convergence of the mimetic finite difference method for diffusion problems on polyhedral meshes[END_REF]. Let α = (α K ) K∈M be a family of strictly positive reals and let α ≥ α > 0 such that {α K , K ∈ M} ⊂ [α, α]. We notice that, for any u ∈ X D which is the interpolation of a linear function, the linear form R K,σ (u) vanishes, since

u σ = u K + ∇ K u • (x σ -x K )
. Hence an idea is to look for the solution of a modified scheme, given by

u = argmin v∈X (ε) D,0 J D,α (v), (19)
where, for a given ε > 0, we define

G (ε) K (v) = 1 2 σ∈E K m σ d K,σ n K,σ • Λ K n K,σ R K,σ (v) 2 -m K ε, ∀K ∈ M, (20)
and

X (ε) D,0 = {v ∈ X D,0 , G (ε)
K (v) ≤ 0, ∀K ∈ M}. (21)
Let us remark that, for any function φ ∈ C 2 c (Ω) (that are functions in C 2 (Ω) with compact support), Lemma 2.4 shows that the interpolation P D φ satisfies P D φ ∈ X (ε) D,0 provided h D be sufficiently small. We now have the following lemma. Lemma 3.1 (Characterization of the solution of the constrained problem) Let us assume Hypothesis [START_REF] Andreianov | Discrete duality finite volume schemes for Leray-Lions-type elliptic problems on general 2D meshes[END_REF]. Let D be a discretization of Ω in the sense of Definition 2.1 and let θ ≥ θ D be given, where θ D is defined by [START_REF] Brezzi | Convergence of the mimetic finite difference method for diffusion problems on polyhedral meshes[END_REF]. Let α = (α K ) K∈M α be a family of strictly positive reals and let α ∈ (0, min{α K , K ∈ M}). Let ε > 0 be given. Then there exists one and only one solution u to (19), which moreover satisfies: there exists a family of non negative reals β

= (β K ) K∈M such that (u, β) ∈ X (ε) D,0 × (R + ) M is a saddle point on X D,0 × (R + ) M of the function L defined by L(v, β) = J D,α (v) + K∈M β K G (ε) K (v), (22)
and the so-called Kuhn-Tucker relations

β K G (ε) K (u) = 0, ∀K ∈ M, (23)
hold. Finally, the following relation holds:

u, v D,α+β = Ω f (x)P M v(x)dx, ∀v ∈ X D,0 . (24)
Proof. The existence and the uniqueness of u of ( 19) is an immediate consequence of Lemma 2.2 and of the fact that

X (ε)
D,0 is a closed non empty convex set. The existence of the multipliers β ∈ (R + ) M such that (23) holds is a consequence of Theorem 5.1 given in the appendix. The saddle point property implies that

u = argmin v∈X D,0 J D,α (v) + K∈M β K G (ε) K (v) .
The variational formulation of the above relation exactly is (24). Thanks to the characterization given by Lemma 3.1, we can state the following estimate on the solution of the constrained problem.

Lemma 3.2 (Estimate on the solution of the constrained scheme)

Let us assume Hypothesis [START_REF] Andreianov | Discrete duality finite volume schemes for Leray-Lions-type elliptic problems on general 2D meshes[END_REF]. Let θ > 0, α > 0 be given. Let D be a discretization of Ω in the sense of Definition 2.1 such that θ ≥ θ D , where θ D is defined by [START_REF] Brezzi | Convergence of the mimetic finite difference method for diffusion problems on polyhedral meshes[END_REF]. Let α = (α K ) K∈M be a family of reals such that min{α K , K ∈ M} ≥ α. Let ε > 0 be given. Then, for all β ∈ (R + ) M , there exists one and only one u D solution to (24), which satisfies

u D 1,D ≤ f L 2 (Ω) C 1 α 0 , ( 25
)
where α 0 only depends on d, α, λ and on θ (see Lemma 2.2) and C 1 only depends on d, Ω and on θ (see Lemma 2.1). Moreover, in the case where (u, β) is a saddle point of the function defined by ( 22), then

K∈M β K m K ≤ f 2 L 2 (Ω) C 2 1 2α 0 ε . ( 26
) Proof. Letting v = u in (24) provides u, u D,α+β = Ω f (x)P M u(x)dx. (27)
Since u, u D,α ≤ u, u D,α+β , we get, from Lemmas 2.2 and 2.1,

α 0 [u, u] D ≤ f L 2 (Ω) C 1 ([u, u] D ) 1/2 ,
where α 0 only depends on d, α, λ and on θ (in particular, it does not depend on β). This gives (25). We now assume that (u, β) is a saddle point of the function defined by (22). We then get, from (27),

K∈M β K σ∈E K m σ d K,σ n K,σ • Λ K n K,σ R K,σ (u) 2 ≤ Ω f (x)P M u(x)dx, which provides K∈M β K (G (ε) K (u) + m K ε) ≤ f 2 L 2 (Ω) C 2 1 2α 0 .
Thanks to (23) which holds in this case, we get

K∈M β K m K ≤ f 2 L 2 (Ω) C 2 1 2α 0 ε .
We can now state the following theorem.

Theorem 3.1 (Convergence of the constrained scheme) Let us assume Hypothesis [START_REF] Andreianov | Discrete duality finite volume schemes for Leray-Lions-type elliptic problems on general 2D meshes[END_REF]. Let θ > 0, α ≥ α > 0 be given. Let D be a discretization of Ω in the sense of Definition 2.1 such that θ ≥ θ D , where θ D is defined by [START_REF] Brezzi | Convergence of the mimetic finite difference method for diffusion problems on polyhedral meshes[END_REF]. Let α = (α K ) K∈M be a family of reals such that {α

K , K ∈ M} ⊂ [α, α].
Let ε D > 0 be given, and u D be the unique solution of ( 19) for ε = ε D . Then u D tends in L 2 (Ω) to u, the unique solution of ( 3), as h D and h D / √ ε D tend to 0. Moreover,

∇ D u D tends in L 2 (Ω) d to ∇u.
Proof. Lemmas 3.2 and 2.3 allow to extract from any sequence of solutions defined by a sequence of discretizations, a sub-sequence which converges in L 2 (Ω) to some u ∈ H 1 0 (Ω). Let φ ∈ C ∞ c (Ω) be given. Let D be a discretization belonging to this extracted sub-sequence, and let us take v = P D φ in (24). We get

u D , P D φ D,α + T 1 (u D , P D φ) = Ω f (x)P M P D φ(x)dx,
where T 1 is defined by

T 1 (w, v) = K∈M β K σ∈E K m σ d K,σ n K,σ • Λ K n K,σ R K,σ (w)R K,σ (v), ∀w, v ∈ X D . ( 28
)
Let us apply the Cauchy-Schwarz inequality. We get

T 1 (u D , P D φ) 2 ≤ T 1 (u D , u D )T 1 (P D φ, P D φ).
Since, applying Lemma 2.4, we get that there exists C 4 , only depending on d, θ and φ, such that

|R K,σ (P D φ)| ≤ C 4 h D ,
we can write

T 1 (P D φ, P D φ) ≤ C 2 4 h 2 D λd K∈M β K m K .
The above relations lead to the existence of C 5 , only depending on θ, φ, Λ, d, Ω and α such that

T 1 (u D , P D φ) 2 ≤ C 5 h 2 D ε .
Hence, under the condition that h D / √ ε D tends to 0, we get that T 1 (u D , P D φ) tends to 0 as well. This completes the proof of the convergence of the scheme, since the convergence of u D , P D φ D,α is provided by Lemma 2.5, and the proof of convergence of Ω f (x)

P M P D φ(x)dx to Ω f (x)φ(x)dx is straightforward. This implies that u D , u D D,α converges
to Ω Λ(x)∇u(x).∇u(x)dx and that T 1 (u D , u D ) tends to 0. From the preceding convergence property of u D , u D D,α , we deduce as in [START_REF] Eymard | A new finite volume scheme for anisotropic diffusion problems on general grids: convergence analysis[END_REF] that

∇ D u D converges to ∇u in L 2 (Ω) d .
Let us underline that the above convergence result is proven in the general framework of diffusion matrix Λ which are only assumed to be bounded, without further assumptions. This implies that the only regularity assumed on the continuous solution of (3) is u ∈ H 1 0 (Ω). On the contrary, the following error estimate result is only stated in the particular case Λ = Id, and u ∈ C 2 (Ω). Although these assumptions are quite restrictive (for example, we could handle the cases where Λ is piecewise C 1 and u is piecewise H 2 ), they prevent from entering into too technical details but allow for getting an indication of the origin of the error (in many similar cases, a priori error estimates for finite volume methods are known to be not sharp and pessimistic compared to the results obtained numerically). Theorem 3.2 (Error estimate) Let us assume Hypothesis [START_REF] Andreianov | Discrete duality finite volume schemes for Leray-Lions-type elliptic problems on general 2D meshes[END_REF]. We also assume that Λ = Id, and that the solution u of (3) satisfies u ∈ C 2 (Ω).

Let θ > 0, α ≥ α > 0 be given. Let D be a discretization of Ω in the sense of Definition 2.1 such that θ ≥ θ D , where θ D is defined by [START_REF] Brezzi | Convergence of the mimetic finite difference method for diffusion problems on polyhedral meshes[END_REF]. Let α = (α K ) K∈M be a family of reals such that {α K , K ∈ M} ⊂ [α, α]. Let ε D > 0 be given, and u D be the unique solution of ( 19) for ε = ε D . Then there exists C 6 depending only on d, Ω, θ, α, α and u such that:

u D -P D (u) 1,D ≤ C 6 h D √ ε D + h 2 D 1/2 , ( 29
)
there exists C 7 depending only on d, Ω, θ, α, α and u such that:

P M u D -u L 2 (Ω) ≤ C 7 h D √ ε D + h 2 D 1/2 , ( 30
)
and there exists C 8 depending only on d, Ω, θ, α, α and u such that:

∇ D u D -∇u L 2 (Ω) d ≤ C 8 h D √ ε D + h 2 D 1/2 . ( 31
)
Remark 3.1 The above error estimate is in accordance with the convergence theorem 3.1, which requires that h D √ ε D tends to zero.

Proof. Let v ∈ X D,0 . We integrate (1) (which resumes to -∆u = f ) in K ∈ M, we multiply by v K and we sum the result on K ∈ M. We get

- K∈M σ∈E K (v K -v σ ) σ ∇u(x) • n K,σ ds(x) = Ω f (x)P M v(x)dx.

This expression implies that

- K∈M σ∈E K (v K -v σ )m σ ∇ K P D u • n K,σ = Ω f (x)P M v(x)dx + T 2 (v), defining T 2 (v) by T 2 (v) = - K∈M σ∈E K (v K -v σ )m σ ∇ K P D u - 1 m σ σ ∇u(x)ds(x) • n K,σ .
Hence we get

K∈M m K ∇ K P D u • ∇ K v = Ω f (x)P M v(x)dx + T 2 (v),
which leads to

P D u, v D,α = Ω f (x)P M v(x)dx + T 2 (v) + T 3 (v), defining T 3 (v) by T 3 (v) = K∈M α K σ∈E K m σ d K,σ R K,σ (P D u)R K,σ (v).
We can then write, subtracting (24),

P D u -u D , v D,α = T 2 (v) + T 3 (v) + T 1 (u D , v). ( 32
)
We get, thanks to Lemma 2.4 and to the Cauchy-Schwarz inequality,

T 2 (v) ≤ C 9 h D v 1,D ,
where C 9 only depends on u, θ, Ω and d. Similarly, thanks again to Lemma 2.4, we have

|R K,σ (P D u)| ≤ C 10 h D ,
where C 10 only depends on u, θ, Ω and d. Thanks to the Cauchy-Schwarz inequality, we then get that

T 3 (v) 2 ≤ C 2 10 h 2 D d m Ω α 2 K∈M σ∈E K m σ d K,σ R K,σ (v) 2 .
Remarking that, thanks to the Young inequality, we have

R K,σ (v) 2 ≤ 2 (v σ -v K ) 2 d 2 K,σ + |∇ K v| 2 |x σ -x K | 2 d 2 K,σ
using ( 4) and ( 9), we conclude that

T 3 (v) ≤ C 11 h D v 1,D ,
where C 11 only depends on u, θ, Ω, α and d. We then take v = P D u -u D . The last term in the right hand side of (32) can be handled as follows:

T 1 (u D , P D u -u D ) ≤ T 1 (u D , P D u) ≤ C 12 h D √ ε ,
where C 12 is computed following the same steps as in the proof of Theorem 3.1. On the other hand, we have

α 0 P D u -u D 2 1,D ≤ P D u -u D , P D u -u D D,α ,
where α 0 only depends on θ, α and d. Gathering these results, we get

α 0 P D u -u D 2 1,D ≤ (C 9 + C 11 )h D P D u -u D 1,D + C 12 h D √ ε .
Thanks to the Young inequality, we get

(C 9 + C 11 )h D P D u -u D 1,D ≤ 1 2 α 0 P D u -u D 2 1,D + C 13 h 2 D ,
where C 13 only depends on α 0 , C 9 and C 11 . We finally get

1 2 α 0 P D u -u D 2 1,D ≤ C 13 h 2 D + C 12 h D √ ε ,
which concludes the proof of (29). The proof of (30) is a straightforward consequence of Lemma 2.1, and that of (31) is an easy consequence of Lemma 2.4 and of (9).

Implementation and numerical results

In order to find an approximation of the solution u to Problem (19), let us consider the use of Uzawa's algorithm, consisting, for a given real ρ > 0 and an initial family β (0) ∈ (R + ) M (we take it null in the examples below), in the definition of the sequence (u (n) , β (n) ) n∈N such that

u (n) = argmin u∈X D,0 L(u, β (n) ), β (n+1) K = max(β (n) K + ρG (ε) K (u (n) ), 0), ∀K ∈ M, ∀n ∈ N. (33)
Indeed, this algorithm is such that, at each iteration, the minimization problem u (n) = argmin u∈X D,0 L(u, β (n) ) to be solved is given by (24), hence the structure of the code for solving the non-constrained problem is not modified (it suffices to replace α by α + β (n) ). Thanks to Estimate (25), we get that the hypotheses of Theorem 5.2 given in the appendix, stating the convergence of Uzawa's algorithm, are verified. Nevertheless, from a practical point of view, since the constraints are arbitrarily chosen by fixing the value of ε, we are not interested in finding an accurate approximation of the saddle point (which can be too much expensive). We therefore stop the iteration procedure once the family

u (n) is such that G (ε) K (u (n)
) ≤ 0 for all K ∈ M. Moreover, we do not allow the Lagrangian multipliers to decrease within the iteration procedure, which is ensured by the following slight modification of (33)

u (n) = argmin u∈X D,0 L(u, β (n) ), β (n+1) K = max(β (n) K + ρG (ε) K (u (n) ), β (n) K ), ∀K ∈ M, ∀n ∈ N.
(34) Remark 4.1 Although the above procedure does not ensure the convergence to the solution of Problem (19), we have numerically observed that the solution was only weakly modified by using (34) instead of (33).

As we show below, this procedure allows to reach significant improvement in accuracy and stability, compared to nonconstrained scheme, within a small number of iterations. Note that values for the numerical parameters must be selected for running the algorithm, and that, in particular, the selection of an optimal value for ρ appears to be a complicate problem. Although theorem 5.2 suggests an order of magnitude for ρ, namely α 2M using the notations of the statement, which could indicate an order 1/ε 2 for ρ since M should behave as the square of the gradient of the constraints, we have preferred in the following results to study the sensitivity of the results with respect to the values of the parameters.

Anisotropic and heterogeneous case with an analytical solution

This test is inspired from [2; 12], and induces numerical locking for some schemes. It has been selected as test case in the Benchmark on discretization schemes for anisotropic diffusion problems on general grids [START_REF] Herbin | Benchmark on discretization schemes for anisotropic diffusion problems on general grids[END_REF]. In the domain Ω = (0; 1) 2 , we consider a case where the diffusion operator Λ is heterogeneous and anisotropic, and given by

Λ(x 1 , x 2 ) = 1 (x 2 1 + x 2 2 ) 10 -3 x 2 1 + x 2 2 (10 -3 -1)x 1 x 2 (10 -3 -1)x 1 x 2 x 2 1 + 10 -3 x 2 2 .
The analytical solution and the right-hand-side of (3) are defined in this test case by u(x 1 , x 2 ) = sin πx 1 sin πx 2 , f = -div(Λ∇u).

We now compare the solution provided by different grids and different selections for the numerical parameters. We present in Figure 1 the two irregular grids used for this test, also performed with regular square grids (not shown). The interest Figure 1: Left: Coarse irregular grid (mesh 1). Right: Refined irregular grid (mesh 2).

of these irregular grids is that they include a few features of irregular grids generated by the geological observations in the underground engineering framework. In order to get a comparison point, we have compared the results obtained with α = 1 and β = 0 (non-constrained scheme) using mesh 1, with the results obtained with α = 10 -3 and β given by Uzawa's algorithm (constrained scheme) using mesh 1 and mesh 2. Note that the choice α = 1 and β = 0 is the one which leads to two-point fluxes on rectangular meshes and on some regular triangular meshes in the case of isotropic problems (see [START_REF] Eymard | Discretisation of heterogeneous and anisotropic diffusion problems on general non-conforming meshes. SUSHI: a scheme using stabilisation and hybrid interfaces[END_REF]). We thus provide in Figure 2 the difference between the three solutions and the analytical solution, computed along the line with equation x 2 = 0.5. We see that the respect of the constraints decreases strongly the error, which decreases again using the finest mesh. This is confirmed by the L 2 errors of the solution and of its gradient, shown in Table 1, using a variety of grids and numerical parameters. We again used "(nc)" for non-constrained scheme, with α = 1, and "(c)" for constrained scheme. We again remark that the results obtained using the constrained scheme are much more accurate than the results without constraint. Moreover we observe that, as could be expected, the values of the multipliers which are required are lower and lower as the mesh size decrease (recall that, for the size of the mesh tending to 0, the constraints are necessarily satisfied by regular solutions). We also remark that the results are not highly sensitive to the value of ε, which is not the case for the parameter ρ. In all these cases, 3 iterations where sufficient to get a point satisfying all the constraints. We observe that, for the highest values of ρ, we get more precise values of the gradient. This shows that the obtained final point verifies stronger constraints, leading to a more regular gradient. By choosing a sufficiently high value for ρ as the grid size decreases, we observe that the order of convergence is not far from 2, showing that the result proven in this paper is not sharp (note that this is the case for the great majority of the error estimates provided in the finite volume framework). We finally remark in Figure 3 that the highest values of β are obtained on the regular grid and the coarse irregular grid (mesh 1) around the point (1, 1), whereas they are obtained at the location of the highest perturbations of the refined irregular grid mesh 2.

- Figure 2: Profiles with irregular grids along the line x 2 = 0.5: Difference with the analytical solution computed using the non-constrained scheme with the coarse irregular grid, the constrained scheme with the coarse irregular grid and the constrained scheme with the refined irregular. The conclusion of this first test case is that the constrained scheme leads to accurate tuning of the Lagrange multiplier, ensuring the best precision for the solution and its gradient.

Anisotropic case without source terms

We consider a test similar to the test 3 of the Benchmark on discretization schemes for anisotropic diffusion problems on general grids [START_REF] Herbin | Benchmark on discretization schemes for anisotropic diffusion problems on general grids[END_REF]. In this test, the main directions of an anisotropic diffusion matrix are tilted with respect to the boundary conditions and the mesh. In the domain Ω = (0; 1) 2 , we consider a non-homogeneous Dirichlet problem, without right-hand-side. Let us consider Problem [START_REF] Aavatsmark | A compact multipoint flux approximation method with improved robustness[END_REF] where the diffusion Λ is homogeneous, and given by:

Λ = R θ 1 0 0 10 -1 R -1 θ ,
where R θ is the matrix of a rotation of angle θ = 40 degrees and f = 0. The non-homogeneous boundaries conditions are continuous and piecewise linear on ∂Ω, and such that

u(x) =       
1 on ((0.; 0.2) × 0. ∪ 0. × (0.; 0.2)) 0 on ((0.8; 1.) × 1. ∪ 1. × (0.8; 1.))

1 2
on ((0.3; 1.) × 0. ∪ 0. × (0.3; 1.))

1 2
on ((0.; 0.7) × 1. ∪ 1. × (0.; 0.7))

We compare the solution, again obtained with a reference 200 × 200 mesh, α = 1 and β = 0, with the results obtained using a regular 20 × 20 mesh and using the irregular grid mesh 2 (see Figure 1). Table 1: Comparison of the errors of the solution and its gradient, using different grids and values of numerical parameters.

constrained scheme, the iterations start from the value α = 10 -9 (value also used for the computations without constraint). We see in Figure 4 that the approximate solutions are in good agreement, as far as one can compare piecewise constant functions of different meshes. In any case, this qualitative comparison is not sufficient to give precise indications about the precision and the regularity of the solutions. These indications are provided by the profiles of the solutions along given lines. Indeed, we present in Figures 5 and6 the solution along the lines of equations x 2 = 0.2 and x 2 = 0.4. We first notice that the unconstrained solution with α = 10 -9 is dramatically far from the reference solution, and that it presents numerous oscillations in the case of the irregular grid. On the contrary, the constrained solution happens to be close to the reference one, and is regular even in the case of the irregular grid.

In order to get indications of the capability of the fluxes provided by the constrained scheme to be used for the coupled transport of any scalar or vector quantity, we have computed the streamlines given by the approximate solution. The method is classical in 2D: we compute a scalar potential, initialized by 0 at a corner of the domain, and increased along each edge of the mesh by the oriented value of the flux across the edge. Hence we get a value of the potential at all the vertices of the mesh, the isovalues of which are the streamlines along the velocity -Λ∇u. We observe that these streamlines, given in Figure 7, mainly follow the behavior of that provided by the reference solution (although there clearly remain effects of the distortion of the grid-blocks). Table 2 provides the number of iterations, with respect to various values of the numerical parameters, for the regular 20 × 20 coarse grid. Parameter ρ is the length of the step along the gradient in Uzawa's algorithm (see (39)). Concerning the behavior of the method with respect to the numerical parameters, we observe similar results to those obtained in Section 4.1. Again, the strategy allowing the only increase of the Lagrangian multipliers within the iterative Uzawa's procedure provides β min and β max -values strongly dependent of the parameter ρ but little sensitive to the values of ε, which determine the level for the approximate constraints. Table 3 provides the number of iterations, with respect to various values of the numerical parameters, for the irregular grid. We see that the same conclusions, applied to the case of the regular grid, hold for the irregular one. We remark that the number of iterations is higher in the case of the irregular grid, for the same choices of the numerical parameters. This seems to be compatible with the supplementary variations in the approximate solution due to the distortion of the grid blocks.

Two wells with anisotropy

We consider a test case which is very similar to test 9 described in [START_REF] Herbin | Benchmark on discretization schemes for anisotropic diffusion problems on general grids[END_REF], also provided in [START_REF] Aavatsmark | A compact multipoint flux approximation method with improved robustness[END_REF], which is focused on a problem with zero boundary flow, with two values of the unknowns imposed within two interior grid blocks of a given mesh. It is well-known that, under this form, this problem cannot be considered for a convergence study. We have therefore modified the data in order that this convergence to a reference solution could be numerically observed, providing the opportunity of exploring the behavior of the scheme on meshes with highly contrasted grid block sizes (such meshes are more and more used in oil engineering for modeling the source terms, since the classical approach by well indices fails in heterogeneous anisotropic or tilted wells cases). We consider the domain Ω = (0; 1) × (0; 1) \ (Ω ∂Ω 1 and u = 0 on ∂Ω 2 . The diffusion Λ is assumed to be homogeneous, and given by:

Λ = R θ 1 0 0 10 -3 R -1 θ ,
where R θ is the matrix of a rotation of angle θ = 67.5 degrees and f = 0. We have computed this problem with two grids. The coarse one is depicted in the left part of Figure 8. The solution, obtained with the refined one (which is the coarse grid refined by a factor 10), is given in the right part of Figure 8. Note that the solution obtained on the refined grid Figure 8: Left: The coarse grid for the two wells with anisotropy test case. Middle: the solution obtained with the coarse grid using the constrained scheme. Right: the solution obtained with the refined grid.

ensures u(x) ∈ [0; 1], which is not the case for the coarse solution, whose minimum value is -.04 and maximum value is 1.04: in this case, the maximum principle is not satisfied by the coarse solution. We again get more precise indications of the regularity and the precision of the solutions by the study of the profiles of the solution along selected lines. We have therefore selected lines which cross the domains Ω 1 and Ω 2 . We observe in Figure 9 that, again, the unconstrained solution is far from the reference one, but that the constrained solution is acceptably accurate and regular. We give in We remark in this table that the values of the Lagrange multiplier belong to the range [2; 7000] (the highest values taken around the wells) for the first three choices, and to the range [21; 70000] for the last line, showing again that the convergence is not reached. As in Section 4.1, the final value for the multiplier, provided by this algorithm, strongly depends on the value of ρ, but again does not seem to depend on ε. We observed very close numerical solutions in all the cases presented in this table. We again remark that taking high value of ρ allows to satisfy the constraints within a very small number of iterations.

One well in a distorted quadrilateral domain

We finally consider a test case which is very similar to test 8 described in [START_REF] Herbin | Benchmark on discretization schemes for anisotropic diffusion problems on general grids[END_REF], again also provided in [START_REF] Aavatsmark | A compact multipoint flux approximation method with improved robustness[END_REF]. Since, again, our aim is to examine a convergence behavior, we modify it slightly, replacing the Dirac source term by a boundary condition p = 1 on the boundary of a very small polygonal domain Ω 1 with vertices given in the trigonometric order by the coordinates (x 1 , x 2 ) equal to (.481, .0156), (.501, .0156), (.499, .0176) and (.479, .0176). Hence, the domain is given as Ω = Ω 2 \ Ω 1 , where Ω 2 is the polygonal domain with vertices given in the trigonometric order by (0., 0.), (1., 0.), (-0.0192, .0333) and (0.9808, .0333). Then the boundary condition p = 0 is prescribed on ∂Ω 2 . We show in Figure 10 the coarse mesh (the refined one is 10 times finer), and the values of the solution (comprised between -10 -3 and .947 for the coarse one using the constrained scheme, between -10 -24 and .994 for the fine one). We see that again, the maximum principle is not respected by the approximate solution. The profiles given in Figure 11 show the efficiency of Figure 10: Top: The coarse grid for the one well in a distorted quadrilateral domain test case. Middle: zoom on the solution obtained with the coarse grid using the constrained scheme. Bottom: zoom on the solution obtained with the refined grid.

the constrained scheme (note that the unconstrained one is again far from the solution, and present dramatic oscillations). We show in Table 5 the number of iterations, together with the dependence of the final value of β max on ρ.

Conclusion

The method proposed in this paper aims at tuning the parameters of a numerical scheme, in this case the hybrid finite volume scheme, with identifying them as the Lagrange multipliers in a minimization problem under constraints. A series of numerical examples shows that this identification is successful, and that the constrained scheme satisfies the required criteria of stability and accuracy, even on distorted meshes and in the case of highly anisotropic problems. Note that further studies remain to be driven for the a priori assessment of the numerical parameters.

Nevertheless, this method does not cure in all cases one of the difficult up-to-date problems, which consists in computing solutions which do not violate the maximum principle. A research direction, opened by this paper, relies in formulations of this problem by minimization problems under constraints.

Appendix: the Lagrange multipliers (Kuhn-Tucker) theorem and Uzawa's algorithm

Let us first recall, for the sake of completeness, classical results which hold for regular minimization problems under constraints. Among many possible references, we refer to [START_REF] Culioli | Introduction à l'optimisation, chapter Optimisation non-linéaire sous contraintes[END_REF] for proofs of these results.

Theorem 5.1 (Lagrange multipliers) Let V be a finite dimensional Euclidean space and let K be the convex closed subset of V , defined by

K = {v ∈ V, G i (v) ≤ 0, for 1 ≤ i ≤ p } , ( 36
)
where p ∈ N ⋆ and for all i = 1, . . . , p the real function G i : V → R is convex and continuously differentiable. We assume that the set {v ∈ V, G i (v) < 0, for all 1 ≤ i ≤ p} is non empty (sufficient condition for the qualification of the constraints). Let J : V → R be a continuously differentiable strictly convex function such that lim |u|→∞ J(u) = +∞, and let u ⋆ be the unique solution of the minimization problem u ⋆ = argmin u∈K J(u).

(37)

Then there exists β ⋆ = (β ⋆ i ) 1≤i≤p ∈ (R + ) p such that (u ⋆ , β ⋆ ) is a saddle point on V × (R + ) p of the function L : V × R p → R defined by

L(u, β) = J(u) + p i=1 β i G i (u).

This means that

L(u ⋆ , β) ≤ L(u ⋆ , β ⋆ ) ≤ L(u, β ⋆ ), ∀(u, β) ∈ V × (R + ) p .

Moreover, the so-called Kuhn and Tucker relations hold:

     ∇J(u ⋆ ) + p i=1 β ⋆ i ∇G i (u ⋆ ) = 0, β ⋆ i G i (u ⋆ ) = 0, ∀i = 1, . . . , p, (38)
are satisfied. Reciprocally, if there exists (u ⋆ , β ⋆ ) ∈ K × (R + ) p such that relations (38) are satisfied, then u ⋆ = argmin u∈K J(u) and (u ⋆ , β ⋆ ) is a saddle point on K × (R + ) p of the function L.

Let us now recall Uzawa's algorithm, to find an approximation of the solution u ⋆ of (37). Let ρ > 0 be a given real and let β (0) ∈ (R + ) p be given. We define the sequence (u (n) , β (n) ) n∈N by

u (n) = argmin u∈V L(u, β (n) ), β (n+1) i
= max(β

(n) i
+ ρG i (u (n) ), 0), ∀i = 1, . . . , p, ∀n ∈ N.

(39) Theorem 5.2 (Convergence of Uzawa's algorithm) Let V be a finite dimensional euclidean space and let K be the convex closed subset of V , defined by (36), where p ∈ N ⋆ and for all i = 1, . . . , p the real function G i : V → R is convex and continuously differentiable. We assume that the set {v ∈ V, G i (v) < 0, for all 1 ≤ i ≤ p} is non empty. Let J : V → R be a continuously differentiable function such that there exists α > 0 with

(∇J(u) -∇J(v), u -v) ≥ α u -v 2 , ∀u, v ∈ V, (40)
(then J is called "α-elliptic", which is sufficient to show that J is strictly convex and verifies lim |u|→∞ J(u) = +∞).

Let us assume that there exists B ≥ 0 such that, for all β ∈ (R + ) p , argmin u∈V L(u, β) ≤ B. Let M be defined by M = max{ p i=1 ∇G i (u) 2 , u ≤ B}. Then, for all ρ ∈ (0, α 2M ) and for all β (0) ∈ (R + ) p , the sequence defined by ( 39) is such that (u (n) ) n∈N converges to the solution u ⋆ of (37).

Figure 3 :

 3 Figure 3: Distribution of β on 80 × 80 grid (left), mesh 1 (middle), mesh 2 (right).

Remark 4 . 2

 42 In all colored figures representing the field of a scalar value, the red color stands for the highest value, the dark blue color for the lowest value.

Figure 4 :

 4 Figure 4: Left top: solution with fine grid. Middle top: constrained solution with regular coarse grid. Right top: constrained solution with irregular grid. Middle bottom: values of β with regular coarse grid. Right bottom: values of β with irregular grid.

Figure 5 :

 5 Figure 5: Profiles with coarse regular and fine regular grids along the line x 2 = 0.2 (left) and along the line x 2 = 0.4 (right).

Figure 6 :

 6 Figure 6: Profiles with irregular and fine regular grids along the line x 2 = .2 (left) and along the line x 2 = .4 (right).

Figure 7 :

 7 Figure 7: Stream lines: fine regular mesh (left), regular 20 × 20 grid (middle), irregular grid (right)

Figure 9 :

 9 Figure 9: Left: Profiles along the line x 2 = .5. Middle: Profiles along the line x 1 = .318. Right: Profiles along the line x 1 = .682.

Figure 11 :√ 3 =

 113 Figure 11: Left: Profiles along the line x 2 = .00835. Middle: Profiles along the line x 2 = .0167. Right: Profiles along the line x 1 + x2 √ 3 = .5.

Table 2 :

 2 [START_REF] Aavatsmark | A compact multipoint flux approximation method with improved robustness[END_REF] ∪ Ω 2 ), where Ω 1 = (0.308; 0.328) × (0.49; 0.51) and Ω 2 = (0.672; 0.692) × (0.49; 0.51), which represent two wells, with u = 1 on Number of iterations for the regular 20 × 20 coarse grid.

	ε	Nb. iterations	ρ	β min	β max
	10 -7	3	10 6	2	150000
	10 -7	10	10 5	0.2	15000
	10 -6	10	10 5	0.1	15000
	10 -6	22	10 4 0.01	1500
	ε	Nb. iterations	ρ	β min	β max
	10 -7	11	10 6	4.5	300000
	10 -7	18	10 5	0.5	30000
	10 -6	16	10 5 10 -9	30000
	10 -6	31	10 4 10 -9	3000

Table 3 :

 3 Number of iterations for the irregular grid.

Table 4

 4 the relation between the number of Uzawa's iterations and the range of the Lagrange multipliers, for a variety of numerical parameters.

	ε	Nb. iterations	ρ	β min β max
	10 -5	3	10 5	2	7000
	10 -7	40	10 5	2	7000
	10 -8	45	10 5	2	7000
	10 -8	3	10 6	21	70000

Table 4 :

 4 Number of iterations and range of the Lagrange multiplier.

Table 5 :

 5 Number of iterations and range of the Lagrange multiplier.