Convergence of nonlinear numerical approximations for an elliptic linear problem with irregular data

. This work is devoted to the study of the approximation, using two nonlinear numerical methods, of a linear elliptic problem with measure data and heterogeneous anisotropic diffusion matrix. Both methods show convergence properties to a continuous solution of the problem in a weak sense, through the change of variable 𝑢 = 𝜓 ( 𝑣 ), where 𝜓 is a well chosen diffeomorphism between ( − 1 , 1) and R , and 𝑣 is valued in ( − 1 , 1). We first study a nonlinear finite element approximation on any simplicial grid. We prove the existence of a discrete solution, and, under standard regularity conditions, we prove its convergence to a weak solution of the problem by applying H¨older and Sobolev inequalities. Some numerical results, in 2D and 3D cases where the solution does not belong to 𝐻 1 (Ω), show that this method can provide accurate results. We then construct a numerical scheme which presents a convergence property to the entropy weak solution of the problem in the case where the right-hand side belongs to 𝐿 1 . This is achieved owing to a nonlinear control volume finite element (CVFE) method, keeping the same nonlinear reformulation, and adding an upstream weighting evaluation and a nonlinear 𝑝 -Laplace vanishing stabilisation term.


Introduction
This paper is devoted to the numerical approximation of a second order linear elliptic equation in divergence form with coefficients in ∞ (Ω) and measure data. More precisely we consider the following problem: find a function defined on Ω such that −div(Λ∇ ) = in Ω, (1.1) supplemented with the boundary condition = 0 on Ω, (1.2) under the following assumptions: − Λ ∈ ∞ (Ω) × and there exists , > 0 such that, for a.e. ∈ Ω, Λ( ) is symmetric and, for all ∈ R , | | 2 ≤ Λ( ) · ≤ | | 2 , (1.3b) − ∈ (Ω), where we denote by (Ω) the space of Radon measures, defined as the topological dual space of (︀ Ω )︀ , the space of continuous functions on Ω with its usual norm. (1.3c) This type of problem may arise for example in some models of underground oil or water resources management. Then is the source term of the quantity diffused by the system. In the case where it is possible to consider ∈ 2 (Ω), the natural framework of the problem consists in considering the standard weak formulation of (1.1) in 1 0 (Ω). But in the case where is modeling a singular source term at the scale of the domain, the right-hand side can be a measure instead of a function. For example, in underground fluid management [19], the wells, which are used for injection or production purposes, are cylindrical holes with very small diameters compared to the size of the domain. Therefore, the source terms are accurately modeled by measures supported by lines.
A proper mathematical sense to a solution of problems (1.1) and (1.2) under Assumptions 1.3 can be given as follows. The existence of a weak solution for any ≥ 2 is given in [4]. Its uniqueness is proved for = 2 in [20] for general diffusion fields: the proof relies on a regularity result [22] which holds on domains Ω with 2 boundary, extended in [21] to all domains with Lipschitz boundaries. In the case ≥ 3, this uniqueness result remains true if Λ is sufficiently regular (see [15] for such a proof when Λ = Id), but it is no longer true for general diffusion fields (as the ones introduced in the numerical examples in Sect. 4, inspired by the counter-example introduced by Serrin [24] and detailed by Prignet [23]). This paper is focused on the approximation of problem (1.5). Consistently with the space (Ω) introduced in Definition 1.1, we investigate weak/strong convergence of approximate solutions in 1, (Ω) for ∈ (︁ 1, −1 )︁ . Let us cite a few different works providing such convergence results.
In [15], a finite volume method is used. It relies on the Two-Point Flux Approximation in grids satisfying an orthogonality condition, which is restricting the kind of meshes which can be used and the kind of anisotropy that can be considered for Λ. The convergence of the approximate solution to the unique weak solution of the problem is proved in the case Λ = Id. The method of proof implies the use of nonlinear functions of the discrete unknown as test function, which is easily managed by the use of Two-Point Flux Approximation (see [18] for a discussion about this problem). Since our goal is to consider general diffusion fields Λ, this scheme no longer applies.
The standard finite element framework is considered in [8]: it consists in computing the solution to the following linear discrete problem ∈ (Ω) and ∫︁ Ω Λ ∇ · ∇ d = ∫︁ Ω d , for any ∈ (Ω), (1.6) where (Ω) is the finite dimensional space resulting from the use of the 1 finite elements on triangles or tetrahedra, and Λ is a piecewise constant approximation of Λ. This solution can always be computed, since (Ω) ⊂ (︀ Ω )︀ . But the study of its convergence cannot be done in the classical way, which consists in taking = for getting an estimate. Indeed, for > 1, this method does not yield an estimate of ‖ ‖ (Ω) , which would be necessary for passing to the limit. Nevertheless, the convergence of the approximate solution to a weak solution (which is moreover the renormalised solution or equivalently the entropy solution in the sense of Def. 3.1) can be proved in the case where the finite element scheme is similar to that used in the finite volume framework, that is when it relies on a Two-Point Flux interpretation of the finite element scheme. This requires the use of 1 finite elements on triangles or tetrahedra which satisfy strong geometrical constraints (in 2D, for two triangles sharing the same edge, the sum of the two opposite angles must be lower than in the case Λ = Id). The assumption of general diffusion fields Λ only satisfying Assumption (1.3b) is incompatible with this method of proof. We therefore propose in this paper two nonlinear numerical methods for the approximation of this problem, with different goals for the two methods. The first one concerns an accurate nonlinear finite element scheme, for which a proof of convergence to a weak solution remains available on general simplicial meshes and fields Λ. The second one concerns a nonlinear Control-Volume Finite-Element scheme (CVFE for short) for which, in addition, a proof of convergence to the entropy weak solution holds in the case where ∈ 1 (Ω).
Notes applying to the whole paper -We denote by ‖ · ‖ the norm of (Ω) for any ≥ 1.
-The Lebesgue measure of any subset of R is denoted by | |.

Motivation and organisation
Section 2 of this paper explores a non-linear finite element formulation of the problem such that, when taking the solution of the scheme as test function, we get an estimate leading to a convergence property to some function ∈ (Ω). Hence we introduce a strictly increasing diffeomorphism : (−1, 1) → R, and we replace the function by the function ( ) where (Ω) ⊂ (−1, 1). A similar idea, used in [5] for getting such an if ≥ 0 and −1 if < 0. Unfortunately, there is no fixed value for providing the estimate ∈ 1, 0 (Ω) for all ∈ (1, /( − 1)): one has to let → 0 for covering the whole desired range of values for . So this choice cannot be kept in our framework.
Let us sketch the computations which lead to the expression of chosen in this paper. We first transform the problem (1.5) into the problem, formally given as follows: find a function : Ω → (−1, 1), such that, for any regular function vanishing at the boundary, Assuming that can be taken as a test function, and denoting for all ∈ (−1, 1) by ( ) = ∫︀ 0 √︀ ′ ( ) d , the following relation holds From this relation, using Hölder and Sobolev inequalities, an estimate on ‖∇ ( )‖ , for all ∈ (1, /( − 1)), is derived under the condition that, up to a quantity that can be controlled, ′ is bounded by . This leads to the following definition for (see Fig. 1 for a graphical representation of , ′ and ): : for a given > 0, is considered in the numerical section. Since choosing values ̸ = 1 does not change the methods of the mathematical analysis, we mainly let = 1 in this paper in order to avoid additional notations.
Note that for any ∈ (−1, 1), we have and the reciprocal function to is the function −1 : R → (−1, 1) defined for any ∈ R by Note also that this function satisfies |( −1 ) ′ ( )| ≤ 1 for any ∈ R. It is noticeable that the function ↦ → log(1+| |) also plays an important role in the analysis done in [14] or [15]. It is then possible to define a finite element formulation of the nonlinear problem (2.1), and to prove its convergence to a weak solution of the problem in the sense of Definition 1.1, hence providing a first proof for the convergence of a numerical scheme, using a general simplicial mesh, for a general diffusion field and an irregular right-hand side. Note that we are not able to prove, in the case where ∈ 1 (Ω), the convergence of this scheme to the entropy weak solution in the sense of Definition 3.1 given below, contrary to the modified scheme studied in Section 3. But, as shown in Section 4, this modification leads to a severe loss of numerical convergence order. Section 2 is organized as follows. In Section 2.2, we present the numerical scheme and Section 2.3 provides the main results that we are able to prove on this scheme. In Section 2.4, estimates of the discrete solution are established, following the ideas sketched above. These estimates allow, in Section 2.5, to prove the existence of at least one solution to the scheme, using the now classical topological degree arguments. Section 2.6 is devoted to the convergence proof of a subsequence of discrete solutions to a weak solution of the problem (recall that the uniqueness of this solution holds only if = 2). The proof of the weak convergence of a sequence of approximate solutions is based first on the compactness of the sequence of approximate solutions and then on the identification of the limit.

The numerical scheme Elements and nodes
We define a simplex in dimension ≥ 1 as the interior of the convex hull of a given set of + 1 points (called its vertices) which are not all contained in the same hyperplane (a simplex is a triangle if = 2 and a tetrahedron if = 3).
We consider a finite family of simplices (called the mesh of the domain), which satisfies the following properties.
(1) The family containing all the vertices of the elements of the mesh, called the family of the nodes of the mesh, is denoted by ( ) ∈ , and, for an element ∈ the family of the + 1 vertices of is denoted by ( ) ∈ , with ⊂ . The set is partitioned into = int ∪ ext , where for all ∈ ext , ∈ Ω (the exterior nodes) and for all ∈ int , ∈ Ω (the interior nodes).
(2) The union of the closure of all the elements of is equal to Ω.
(3) The intersection of two different elements of is empty.
Such a family is then a regular simplicial mesh of Ω in the usual sense of the finite element literature [9].
For ∈ , one denotes by the set of all ∈ such that ∈ . For any ∈ , we denote by | | the measure in R of , by ℎ the diameter of and by the diameter of the largest ball included in . Then, we define the mesh size ℎ and the mesh regularity by

basis functions and barycentric coordinates
We denote, for any ∈ , by the continuous function defined on Ω which is piecewise affine in each ∈ , continuous, and such that ( ) = 1 and ( ) = 0 for all ∈ ∖ { }. Recall that, for any ∈ and ∈ , ( ( )) ∈ is the family of the barycentric coordinates of point with respect to the vertices ( ) ∈ of .

Approximation space and scheme
We denote by (Ω) the conforming finite element space defined by We approximate Λ by some piecewise constant function Λ : Ω → R × such that where for any ∈ , Λ ∈ R × is assumed to be symmetric and to satisfy | | 2 ≤ Λ · ≤ | | 2 for any ∈ R . (2.5) The value Λ can be chosen as the mean value of Λ on , but it can also be chosen as the value of Λ at any point of if Λ is sufficiently regular (which is done in the numerical examples of Sect. 4). We then consider the following approximation of problems (1.1) and (1.2). Definition 2.2. We say that is a solution of the numerical scheme if ∈ (Ω), max ∈Ω | ( )| < 1 and We remark that, since Λ , ∇ and ∇ are constant in each element of , this numerical scheme leads to the computation, for each ∈ , of the quantity ∫︀ ′ ( ) d . We provide in Section A.1 the mathematical computation of this quantity.

Main results of Section 2
The first main result of this section is the existence of a solution to the nonlinear scheme. Once we have a discrete solution ∈ (Ω) at hand for all meshes, then we can study the convergence of the scheme when the discretisation parameter ℎ tends to zero. More precisely, consider a sequence (︀ ( ) )︀

≥1
of meshes of Ω in the sense specified in Section 2.2 such that and such that the sequence of regularity is bounded: there exists ⋆ such that ( ) ≤ ⋆ , for any ≥ 1. (2.8) Let (Λ ( ) ) ≥1 be such that (2.4) and (2.5) hold and such that Note that we do not assume that (2.9) holds for = +∞, which enables the convergence result to apply in the cases considered in the numerical section. We then consider a sequence ( ( ) ) ≥1 of solutions with respect to the sequence (︀ ( ) )︀ ≥1 and we study the convergence of this sequence to a solution of the continuous problem in the sense of Definition 1.1. Thus the second main result of this section is that this sequence converges, up to a subsequence, to a weak solution of the continuous problem in the sense of Definition 1.1, as stated by the following theorem.
be a sequence of meshes of the computational domain Ω in the sense specified in Section 2.2 such that (2.7) and (2.8) are satisfied and let (Λ ( ) ) ≥1 be such that (2.4), (2.5) and (2.9) hold. For any ≥ 1, let ( ) be an arbitrary numerical solution in the sense of Definition 2.2 in the case where = ( ) (the existence of ( ) is given by Thm. 2.3). Then, there exists ∈ (Ω), weak solution of problems (1.1) and (1.2) in the sense of Definition 1.1, such that, up to the extraction of a subsequence, ( ( ( ) )) ≥1 converges to weakly in 1, Moreover, for all > 0, ∈ 1 0 (Ω) holds. In the case = 2, the whole sequence converges to the unique solution of the problem.
Remark 2.5. Note also that = ( ) where is the limit of the subsequence ( ( ) ) ≥1 in (Ω) for any The remaining of this section is dedicated to the proof of Theorems 2.3 and 2.4.

Estimates
As done in the introduction, we denote in the whole paper by ‖ · ‖ the norm in (Ω) or in (Ω) for any ∈ [1, +∞]. Following the arguments presented in the introduction, we define the following function : (2.10) In Figure 1, we propose a graphical representation of functions , ′ and .
The following lemma is used in the course of the following computations.
Lemma 2.6. For any ∈ (0, 1), there exists strictly positive reals , and only depending on , such that Proof. We have, for all ∈ (0, 1) and ∈ (−1, 1), Consequently we obtain that Owing to (2.12), this gives which provides the conclusion of the lemma owing to a convexity argument.
We have the following estimate. (2.14) Then there exists 1 depending only on and ‖ ‖ (Ω) such that Proof. On one hand, using the fact that max ∈Ω | ( )| < 1 we have and on the other hand, accounting for (2.5), we have Hence we get The next lemma provides a lower bound on 1 − | | for any function such that the estimate (2.14) holds.

Let
∈ N which will be chosen later, and let us define )︀ . Note that we have > 0 and Since | ( )| < 1 for any ∈ , and using the fact that we get that, for a.e. ∈ Ω,
For this value of , which only depends on and , we get the conclusion of the lemma.

Existence of a solution to Scheme (2.6)
The purpose of this section is to prove Theorem 2.3, which states the existence of a solution to the numerical scheme in the sense of Definition 2.2, by applying the topological degree method [12].
Proof of Theorem 2.3. Let us define the continuous function 19) and the function where for any = ( ) ∈ int ∈ R int , ∈ [0, 1] and ∈ int , the quantity ℱ ( , ) is defined by This mapping is well defined and continuous, since, for any = ( ) ∈ int ∈ R int , we have max ∈Ω |ℒ( )( )| ≤ max ∈ int | −1 ( )| < 1 and ℒ is continuous. We also notice that the equation ℱ( , 1) = 0 gives In particular we obtain which means that ℒ( ) is a solution of the numerical scheme (2.6). Let ∈ (0, 1] and let = ( ) ∈ int ∈ R int be such that ℱ( , ) = 0. We have, for any ∈ int , Multiplying the previous identity by −1 ( ) and summing over the internal nodes leads to Dividing by and using 1− ∑︀ which means that (2.14) holds for ℒ( ). Hence, from Lemma 2.8, we then obtain We then obtain that | −1 ( )| ≤ 2 for all ∈ int , and therefore that Define the relatively compact open set For = 0, the linear equation ℱ( , 0) = 0 has the unique solution = 0. The topological degree corresponding to ℱ( , 0) and is therefore equal to 1 since = 0 belongs to . Hence for any ∈ [0, 1] any solution to ℱ( , ) = 0 necessarily belongs to . Therefore, owing to the invariance of the topological degree by homotopy, there exists at least one ∈ R int (non necessarily unique) such that ℱ( , 1) = 0, which means that ℒ( ) is a solution to Scheme (2.6) in the sense of Definition 2.2.
We can now prove the weak convergence of a discrete solution to a solution of the continuous problem.

Motivation and organisation
As already noticed, for some diffusion fields Λ satisfying Hypothesis (1.3b), there exist non-zero weak solutions to problems (1.1) and (1.2) in the sense of Definition 1.1 for = 0. One such diffusion field is precisely considered in the numerical examples presented in Section 4. But, among these solutions, there is one and only one which is in some sense the limit of more regular problems, as proved in [3]. This solution must satisfy a regularity criterion which is the basis of the notion of entropy weak solution, in the case where ∈ 1 (Ω). We give in Definition 3.1 below this entropy weak solution sense, formulated in the particular case of bounded domains and linear uniformly elliptic operators.
Definition 3.1. Entropy solution to the linear elliptic problem with right-hand side in 1 (Ω).
An entropy solution of problems (1.1) and (1.2) in the sense of [3] is a function ∈ (Ω) satisfying that for any > 0 and for any ∈ ∞ (Ω), It is proved in [3] that it is equivalent to replace (3.1) and the truncations by for any ∈ ∞ (Ω) and for any ∈ ℱ, where ℱ is defined as the set of all functions It is also proved in [3] that the entropy solution exists, is unique and is a weak solution in the sense of Definition 1.1.

Remark 3.3.
In fact, Andrea Dall'Aglio proved in 1996 that the inequality (3.1) can be replaced by an equality [11]. This result is stated and briefly proved in Lemma A.6 in the appendix, as well as the fact that (3.2) can also be replaced by an equality.
The above entropy solution sense is reviewed in Prignet [23], as well as different mathematical senses for a solution to this problem. It is proved to be equivalent to the renormalised sense introduced in [8,10]. See [3,4,25] for more detailed definitions and properties; one can also refer to [13,14] for some extensions to the case where the problem is not coercive.
Recall that, in [8], the authors prove the convergence of the approximate solution obtained by the linear finite element approximation (1.6) to the entropy (or renormalised) solution, but only under strong restrictions on the meshes and Λ-fields that can be considered. These restrictions are not requested for the nonlinear finite element scheme studied in Section 2 of this work, which is shown to converge to a weak solution of the problem on any simplicial grid, with any diffusion field. Nevertheless, the convergence of this nonlinear finite element scheme to the entropy solution of the problem could not be proved, mainly because no inequality can be obtained by introducing nonlinear functions of the primary unknown as test functions.
We therefore propose in Section 3 of this paper a nonlinear scheme which holds for any simplicial grid and for any diffusion field Λ satisfying Hypothesis (1.3b). We construct this scheme in such a way that we can prove its convergence to the entropy weak solution of the problem, owing to three ideas.
(1) The first one is motivated in Section 2.1 above: it relies on writing = ( ), using the function defined by (2.2). Then, as in Section 2, an estimate is deduced when one takes the approximation of the bounded function as test function in the numerical scheme.
(2) The second one is the use of a control volume-finite element scheme (called for short CVFE in this work), which allows to simultaneously adopt the finite element point of view, for computing the coefficients of a rigidity matrix, and the finite volume point of view, for the computation of some nonlinear expressions with respect to the unknown. Then, following [2,6,7] and other works by the same authors, the evaluation of ′ is upstream weighted with respect to the sign of the so-called "transmissibility" between control volumes (this quantity, which only depends on the mesh and on Λ, is defined by (3.4)). Owing to this technique, it becomes possible to derive estimates with using test functions under the form of nonlinear functions of the primary unknown. For this purpose, these nonlinear functions must satisfy some specific properties (see (A.3)), which is restricting the range of the nonlinear functions which can be considered: for example, truncations do not provide such monotony inequalities, see Remark 3.16. As in Section 2, Sobolev inequalities (see Lem. A.1) play an important role for establishing these estimates; the application of these inequalities to the CVFE scheme is done through the equivalence property between continuous and discrete norms (see Lem. A.2). The discrete Sobolev inequalities as stated in [17] could also be used, up to the adaptation of the treatment of the homogeneous Dirichlet boundary conditions. (3) We notice in Remark 3.11 that this upstream weighting yields a kind of weak -Laplace stabilisation term with = 3 (in analogy with the Godunov scheme for scalar hyperbolic equation, which provides a weak BV-inequality). Nevertheless, we could not conclude the convergence proof to the entropy solution only using this weak inequality, and we had to introduce in the scheme a nonlinear vanishing stabilisation term, based on a -Laplace operator with > 3. This term is not necessary for proving, up to a subsequence, the convergence of the scheme to a weak solution; but it is strongly used in the proof of the convergence to the entropy weak solution in the case where ∈ 1 (Ω). One of the difficulties in the definition of this term is that it must be sufficiently large for bounding the expressions which must tend to zero, and sufficiently small for vanishing against regular test functions.
The organisation of Section 3 is similar to that of Section 2. In Section 3.2, we present the CVFE numerical scheme and give the main results (existence and convergence) concerning this scheme in Section 3.3. In Section 3.4, estimates on the discrete solution are established, including an inequality induced by the stabilisation term. These estimates allow, in Section 3.5, to derive the existence of a solution to the scheme, using a topological degree argument (similarly to what is done in Sect. 2). Section 3.6 is first devoted to the convergence proof of a subsequence of discrete solutions to a weak solution of the problem in the general case ∈ (Ω) (recall that the uniqueness of this solution holds only if = 2). Then the convergence of the whole sequence to the entropy weak solution is proved in the case where ∈ 1 (Ω). This last proof strongly relies on the presence of the -Laplace stabilisation term. Let us observe that the complexity of these convergence proofs is largely greater than that of the convergence to a weak solution of the nonlinear finite element scheme studied in Section 2.

Definition of the scheme
The Control Volume Finite Element (CVFE) scheme is based on 1 -finite elements on a primal simplicial mesh, and on piecewise constant functions on a dual mesh (see Fig. 2 for an illustration in the case = 2). Let us introduce the geometrical objects used in the definition of this scheme, which will be considered in the following to be all collected by the notation .

Primal mesh
We use the elements, nodes, 1 basis functions defined in Section 2.2 which define the primal mesh of Ω. the subset of ℰ containing all the edges of , and for any { , } ∈ ℰ, we denote by the set of all ∈ such that { , } ∈ ℰ . Accounting for the fact that any simplex has ( + 1)/2 edges, we define for any { , } ∈ ℰ,

Dual mesh
Once the primal triangular mesh has been built, we can define its dual barycentric mesh ℳ as follows. To each ∈ and ∈ , we define the set , of all ∈ such that ( ) > ( ) for all ∈ ∖ { } (we then have = ⋃︀ ∈ , ). Then we define = ⋃︀ ∈ , and ℳ = { , ∈ }. Note that Ω = ⋃︀ ∈ and ∩ = ∅ for ̸ = . We refer to Figure 2 for an illustration of the primary and dual barycentric meshes in the 2D case. The Lebesgue measure of is denoted by . The geometrical construction of ensures that We then denote by : Ω → R the characteristic function of the subset for any ∈ . The space of all real families ( ) ∈ is classically denoted by R . Then we set For all ∈ R and for any continuous function : R → R, we denote by ( ) the element of R such that ( ( )) = ( ) for any ∈ . Given a family = ( ) ∈ ∈ R , we denote Π ∈ (︀ Ω )︀ and Π ℳ ∈ 1 (Ω) the functions defined by

Transmissibility coefficients
In order to define the CVFE scheme, we define the transmissibility coefficients (whose sign is not specified, contrary to the case of the finite volume-type schemes) and Note that Λ = 0 unless { , } ∈ ℰ. Moreover, since ∑︀ ∈ ∇ = 0, we have that: As a consequence of (3.4) and (3.5), given and two elements of R 0 , one has
Definition 3.4. Let 0 ∈ [0, +∞) and ∈ (2, +∞) be given. We say that is a solution of the numerical scheme if there holds ∈ R 0 , max ∈ | | < 1 and we let 0 = 0 in the numerical results given in this paper, we observe suboptimal convergence orders, resulting from the upstream weighting scheme used in the definition of Ψ , ( , ).

Main results of Section 3
The first main result of this section is similar to that of Section 2.
Theorem 3.6. There exists (at least) one solution ∈ R 0 to Scheme (3.8) (in the sense of Def. 3.4).
Therefore we can consider a sequence (︀ ( ) )︀ ≥1 of meshes of Ω in the sense specified in Section 2.2 such that (2.7) and (2.8) are satisfied. For this sequence, owing to Theorem 3.6, we can consider for all ≥ 1 a solution ( ) to Scheme (3.8) with respect to ( ) . The second main result of this section, also similar to that of Section 2, is that the functions reconstructed from this sequence converge, up to a subsequence, to a weak solution of the continuous problem in the sense of Definition 1.1, as stated by the following theorem.
Moreover, for all > 0, ∈ 1 0 (Ω) holds. In the case = 2, the whole sequence converges to the unique solution of the problem. Finally, we have the following important property in the case ∈ 1 (Ω), which is not proved for the scheme studied in Section 2. The remaining part of this section is dedicated to the proof of Theorems 3.6, 3.7 and 3.9.
We have the following estimate.
As a consequence of the previous result, we can obtain a bound away from one on a function |Π | such that the estimate (3.12) holds. Proof. Owing to Lemma A.1 with = 1 and to (A.1), we can write Consequently using Lemma 3.10 we obtain (2,1) sob which gives in particular for any ∈ , Using the fact that the function is a continuous strictly increasing one-to-one function from (−1, 1) to R, we then obtain for any ∈ int , Remark 3.13. Note that the proof of Lemma 3.12 is shorter and simpler than the proof of Lemma 2.8, owing to the finite volume point of view allowed by the CVFE scheme.

Existence of a solution to CVFE Scheme (3.8)
The purpose of this section is to prove Theorem 3.6, which states the existence of a solution to the numerical scheme in the sense of Definition 3.4, by applying the topological degree method [12]. This proof follows the same lines as that of Theorem 2.3.

Convergence to a weak solution in the general case ∈ (Ω)
The goal of this section is the proof of Theorem 3.7. In this section, it is possible to let 0 = 0, which means that the convergence results to a weak solution also hold without the stabilisation term. We could as well consider ∈ (1, 2] in the stabilisation term under a suitable definition of ( ) in the case = . The first step is the following compactness lemma. (︀ ( ) )︀ weakly converges to ∇̃︀ ( ) in 2 (Ω) and̃︀ ( ) ∈ 1 0 (Ω) which implies that ( ) ∈ 1 0 (Ω) for all > 0.
We now turn to the convergence proof to a weak solution of the continuous problem.

Convergence to the entropy solution in the case ∈ 1 (Ω)
In this section, we assume that the right hand side of (1.1) is defined by ∈ 1 (Ω), and that 0 > 0 and > 3, which are requested in the course of the convergence proof. The next lemma provides a general convergence result due to the presence of the -Laplace stabilisation term. This result is used several times in the proof of convergence to the entropy solution. Lemma 3.20. Let (︀ ( ) )︀ ≥1 be a sequence of simplicial meshes of Ω such that (2.7) and (2.8) hold. Let ( ( ) ) ≥1 be a sequence of solutions to the numerical scheme (3.8) in the sense of Definition 3.4 with 0 > 0 and > 3. Let us define, for any ∈ (0, 3], ∈ [0, +∞), and for any ≥ 1 (dropping some indices in the discrete quantities involved in the right-hand side), Proof. Applying Hölder's inequality with exponents and − where > 3 ≥ we get Applying (3.14) in Lemma 3.10, we get remains bounded for all ≥ 0. Since, from Lemma 3.17, we have that we conclude using the fact > 3 ≥ that which implies that (3.26) holds.
We now turn to the proof of convergence to the entropy solution.
Let us denote = ( ) for any ∈ . We remark that

This implies
We then have The remaining of the proof consists in studying the limit or the limitinf of each of these terms.

Term 11
Using Lemma A.4 and ′ ≥ 0, we have with, denoting by and by the function defined on Ω, equal a.e. to on ∈ , we have ( ) Definition (3.29) for is motivated, on one hand, by defining a piecewise constant function on the elements, on the other hand, by the fact that ′ (︀ )︀ ̸ = 0 implies that max ∈ | ( )| ≤ ( ), used in the proof of the convergence of the gradient. Note that, owing to the inequality | − 0 | ≤ ∑︀ ∈ | − 0 | for any , 0 ∈ and ( ) ∈ , we have the relation Remarking that | ( )− ( )| ≤ ′ ( )| − |, we can apply the Cauchy-Schwarz inequality and the inequality Owing to (A.2) in Lemma A.2 and to (3.13) in Lemma 3.10, this leads to and therefore that lim Up to the extraction of a subsequence we can assume the convergence a.e. of ′ (︁ ( ) )︁ to ′ ( ( ) − ).

Term 111
We have

Term 112
We remark that, for all ∈ (−1, 1), since ′ ( ) ≤ ′ ( ), the function admits the Lipschitz constant ′ ( ). Besides, we have, for all , ∈ (−1, 1), This enables to write, denoting by a bound of ′′ , We obtain, using ′ ( )+ ′ ( ) ≤ ∑︀ ∈ ′ (︀ ( ) )︀ and the Young inequality We get, using ≤ ℎ , We have, by weak convergence of ∇Π in for any 1 < < /( − 1) and strong convergence of ′ ( )∇Π in ′ with 1/ + 1/ ′ = 1 (recall that ′ is bounded, hence dominated convergence applies), Turning to ( ) 212 , we have, referring to the proof of Theorem 3.7, which proves that this term tends to 0 as → +∞. Studying Following the treatment of ( ) 112 , where Young's inequality is replaced by which also shows that this term tends to 0 as → +∞. Hence the conclusion of the study of 21 is We deduce from these equations that . This allows to conclude that (3.2) holds, which implies that is the unique entropy solution of the problem. From the uniqueness property of the limit, we deduce that the convergence properties proved by Lemma 3.18 (except the almost everywhere convergence) hold for the whole sequence.
Let us now turn to the strong convergence of the gradient We let = 0, and, for a given function ∈ ℱ, we again denote by > 0 such that ′ ( ) = 0 for all | | ≥ , and define ∈ (0, 1) such that ( ) > . Letting = ( ( )) in (3.8), we consider all the terms , = 1, 11, 2, . . . as defined above, and we get (3.27), which leads, owing to (3.28) to Since we prove above that the terms This proves that and therefore that The remaining of the proof is dedicated to show that (3.36) implies the strong convergence of ∇Π ( ) ( ) to ∇ . We now denote a given representative of the functions , ∇ and for any ∈ N, of ( ) and of ∇Π ( ) ( ) , defined everywhere in Ω, by the same notation.
Step 1. Construction of a decreasing sequence ( ) ∈N of infinite subsets of N, and of a sequence (Ω ) ∈N of subsets of Ω such that -|Ω ∖ Ω | → 0 as tends to infinity.
Step 2. Diagonal process and convergence almost everywhere of the gradient.

Numerical tests
We present here some illustrations of the behavior of the numerical schemes (2.6) and (3.8) with = 2 or = 3. The implementation of Scheme (3.8) is simply done in the case 0 = 0, using a Picard iteration for handling the term Ψ ( , ).
Let us detail the implementation of Scheme (2.6). We also adopt a simple Picard iteration method for approximating the solution of the resulting systems of nonlinear equations, since we observe good convergence properties in the studied test cases. It consists in computing the sequence such that (0) = 0 and, for all = 1, . . . , , )︁ d for any ∈ . We use the formulas given in Section A.1, which implies to numerically determine which case of equality is satisfied by the values ( ) =1,..., +1 . This is done by comparing the differences between the values with 10 −3 (smaller values lead to less precise results due to the divisions by these differences).
A second problem in the implementation of the method, in the case where the measure is in fact an element of 1 (Ω) and ( ) is singular at some points of the domain, is the computation of the right-hand sides The key point for the precision of the method is to preserve the exactness of the integrals: the use of approximate quadrature formulas with Gauss points leads to very poor accuracy in this case. We complete this goal by replacing where is a dual cell associated with the vertex (the exact shape of has no influence on the precision of the computation), and by computing the exact integration of ∫︀ ( ) d . Finally, let us observe that, in the examples below, we compare numerical solutions computed with polygonal meshes to analytical solutions available on non-polygonal domains (circles, cylinder). These analytical solutions are vanishing at the boundary of these domains, which leads to an error lower than ℎ 2 , where ℎ is the size of the mesh.

Case = 2, measure
We consider the case where Ω = (0, 1) and (4.1) This heterogeneous and anisotropic diffusion field corresponds in 2D to the case described by Prignet [23] and Serrin [24]. We replace the average values of Λ in the elements by the values at the center of gravity of the elements. We let = (0,0) , which corresponds to the analytical solution given by ( ) = − 1 2 log | | (see the right part of Fig. 4 for a representation of the numerical solution).
In this test case, the solution is no longer in 1 0 (Ω) nor in ∞ (Ω). We use triangular meshes which are refined around the point (0, 0) (see the left part of Fig. 4), and we compute ( ) solution to Scheme (2.6), Π solution to the control-volume finite-element scheme (3.8) with  not to the entropy solution; in practice, increasing 0 leads to an increase of the observed numerical error) and ∈ solution to the linear scheme (1.6). These solutions are respectively denoted by I, II and III in the tables below.
We  Table 1, using 5 meshes whose the total number of vertices is denoted by . In this table, we define the meshsize as −1/2 and we compute the order of convergence with respect to the preceding line. The errors on the gradients are computed at the mid-edges of the triangles.
We observe that the non-linear method I provides slightly more accurate results than the linear method III, and that a numerical order, approximately equal to 2, is observed for and an order 1 is observed for ∇ . Method II seems to provide an order 1, both for and ∇ . In this test case, the solution is no longer in 1 0 (Ω) nor in ∞ (Ω) (recall that it belongs to 1 0 (Ω) only for ∈ (︀ 0, 1 2 )︀ ). The right-hand side is then in 1 (Ω), but this is not the case for the product which is not locally integrable around the point 0 (this prevents from using the solution as test function at the continuous level).
We use triangular meshes which are refined around the point (0, 0) (see Fig. 5 for an example of solution computed with such a mesh), and we again compute the solutions to Schemes I given by (2.6), II given by (3.8)  with 0 = 0 and III given by (1.6). Letting = 10 and using the same notations as in the preceding section, we obtain the results provided in Table 2.
The orders of convergence are similar to the ones observed in Table 1.
We then denote, for all = ( 1 , 2 , 3 ), by ( )︀)︀ , with = 3/4, and by for all ∈ Ω. Then the function given by ( ) = (︀ 1 − 2 3 )︀ ( ) is a weak solution in the sense of Definition 1.1, as well as + 1, 2 − 1, 2 for any ( 1 , 2 ) ∈ R 2 ∖ {(0, 0)}. In this test case, none of these solutions belongs to 1 0 (Ω) or to ∞ (Ω). We know from Prignet [23] that, for ( 1 , 2 ) ̸ = (0, 0), the truncated functions 1, 2 , for any > 0, are not in 1 (Ω). This is not the case for . We extend the 2D refined triangular meshes of (︀ 0, 1 2 )︀ (see the preceding section) on Ω by generating prisms with constant height, which are then shared into 3 tetrahedra. Letting = 1, the following 1 (Ω) errors compared to that given by the linear method are provided in Table 3, with respect to the weak solution (the references to the 3 schemes compared in this numerical section being the same as in the 2D case). In this table, we define the meshsize as −1/3 for computing the order of convergence.
One notices that all schemes seem to converge to (recall that we proved that the truncations of the limit belong to 1 0 (Ω), which excludes any solution including a non-zero term 1, 2 − 1 , 2 ). The non-linear method I seems to provide slightly more accurate results than the linear method III, with a numerical order of convergence close to 2 for the values and 1 for the gradients considering the finest meshes. The linear scheme was not expected in this case to numerically converge to any weak solution different from , since in the case where = 0, the only solution of the linear scheme is 0. The numerical order of convergence for method II seems to be closer to 1/2 than to 1, both for and ∇ .

Some open and closed problems
The two nonlinear methods presented in this paper are proved to converge to a weak solution of the linear elliptic problem with general measure data and heterogeneous anisotropic diffusion fields, which does not seem to have been proved for earlier schemes. The convergence of the nonlinear CVFE scheme to the entropy weak solution seems also to be the first one in this framework. But it was necessary to assess the numerical performance of these schemes, which is done in Section 4, where we present numerical results for 2D and 3D cases obtained with these two schemes and with the simpler linear scheme (1.6) (for which no convergence properties are proved for general meshes and diffusion field). We observe that the accuracy of the nonlinear finite element scheme is comparable to that of the simpler linear scheme (1.6), but that of the CVFE scheme is disappointing.
Two questions (at least) remain open problems.
The first one is to prove the convergence for a stronger topology of the gradient of the numerical solution obtained by the nonlinear finite element scheme, since it seems to be observed in the numerical results.
The second one is to prove that the nonlinear finite element scheme or the nonlinear CVFE scheme without the stabilisation term converge to the entropy solution, which could be expected since the entropy solution is in fact designed to be the limit of the solutions to regularised problems.

A.2. Equivalence of norms
Partial proofs of the following lemma are done in the literature, for example in [6], and we only provide the sketch of the proof for the sake of completeness.  > 0, which is increasing with respect to and also depending on , such that, for any ∈ R , 1 where ‖ · ‖ 1, ,ℳ is defined by (3.16).
Proof. We denote by 0 the reference simplex with vertices 0 and all the extremities of the canonical unit vectors. Let be the affine mapping which transforms 0 into ∈ with vertices ( 0 , . . . , ). Then Π ∘ is an affine function on 0 with values ( 0 , . . . , ) at the vertices of 0 . By equivalence of norms in finite dimension, there exists > 0 only depending on and such that 2) follows from the equivalence of norms in finite dimension spaces, from a bound of the coefficients of by ℎ , from ℎ ≤ and from a bound of the coefficients of ( ) −1 by ℎ −1 /| | involving .

A.3. Comparison lemmas
The next lemma plays an important role in the convergence properties of Scheme (3.8). where the function is defined by (2.2) and the functions ,̃︀ are defined by (3.9).