# Convergence Analysis on Unstructured Meshes of a DDFV Method for Flow Problems with Full Neumann Boundary Conditions

Authors: A. Kinfack Jeutsa; A. Njifenjou; J. Nganhou
Journal: Journal of Applied Mathematics
 (2016)

---

## Abstract

A Discrete Duality Finite Volume (DDFV) method to solve on unstructured meshes the flow problems in anisotropic nonhomogeneous porous media with full Neumann boundary conditions is proposed in the present work. We start with the derivation of the discrete problem. A result of existence and uniqueness of a solution for that problem is given thanks to the properties of its associated matrix combined with adequate assumptions on data. Their theoretical properties, namely, stability and error estimates (in discrete energy norms andL
2-norm), are investigated. Numerical test is provided.


---

## Body

## 1. Introduction and the Model Problem

Efficient schemes are required for addressing flow problems in geologically complex media. The most important criteria of efficiency are mass conservation in grid blocks, accurate approximation of Darcy velocity, capability for dealing with anisotropic flow on unstructured grids and diverse heterogeneities (relevant to absolute permeability, porosity, etc.), and robust and easy implementation. Schemes well known in the literature for meeting many of the previous criteria are the following: Mixed Finite Element methods (see, e.g., [1, 2]), Control-Volume Finite Element methods (see, e.g., [3, 4]), Mimetic Finite Difference methods (see, e.g., [5, 6] and references therein), Cell-Centered Finite Volume methods (see, e.g., [7–11] and certain references therein; see also [12–14]), Multipoint Flux Approximation (see, e.g., [15–18] and some contributions to convergence analysis of MPFA O-scheme like [19]) and Discrete Duality Finite Volume methods (DDFV methods for short).The DDFV methods come in two formulations. The first formulation is based on interface flux computations for primary and dual meshes, accounting with the interface flux continuity (see, e.g., [20, 21]) and the second formulation of DDFV is based on pressure gradient reconstructions over a diamond grid (see [22–24]). Note that this second formulation attracted the attention of some mathematicians as Andreianov, Boyer, and Hubert who have greatly contributed to its mathematical development. Indeed, key ideas involved in the pressure gradient reconstruction approach have been generalized by these authors to nonlinear operators of Leray-Lions type. Motivated by the possibility of increasing the order of convergence of the pressure gradient reconstruction method for nonlinear operators, Boyer and Hubert have proposed in [25] the so-called modified DDFV.Beyond flow problems, we find the applications of DDFV methods in many areas: the numerical modeling of the surface erosion occurring at a fluid/soil interface undergoing a flow process in [26], the discretisation of partial differential equation appearing in image processing in [27], the assessment of nuclear waste repository safety in the context of simulating flow, transport in porous media in [28], and so on.For presenting our analysis of DDFV method, let us consider the 2D diffusion problem consisting of finding a functionφ which satisfies the following partial differential equation associated with nonhomogeneous Neumann boundary conditions:(1)
-
div
⁡
K
grad
⁡
φ
=
f
in
Ω
,
(2)
-
K
grad
⁡
φ
·
η
=
g
o
n
Γ
,where Ω is a given open polygonal domain, Γ is its boundary, f and g are two given functions (defined, resp., in Ω and Γ, at least almost everywhere in Lebesgue measure sense), and η is the unit normal vector to Γ outward to Ω. The permeability K
(
x
), with x
=
x
1
,
x
2
t
∈
Ω, may be a full matrix depending on solely spatial variables and obeying the following conditions.(i)
The primary mesh is such that the discontinuity ofK lies on mesh interfaces.
(ii)
Symmetry is(3)
K
i
j
x
=
K
j
i
x
a.e.  in
Ω
,
∀
1
≤
i
,
j
≤
2
.
(iii)
Uniform ellipticity and boundedness are(4)
∃
γ
min
⁡
,
γ
max
⁡
∈
R
+
∗
such that
∀
ξ
∈
R
2
,
ξ
≠
0
,
γ
min
⁡
≤
ξ
T
K
x
ξ
ξ
2
≤
γ
max
⁡
a.e.  in
Ω
,where · denotes the Euclidian norm in R
2 and where K
i
j
· are the components of K satisfying(5)
K
i
j
∈
PWC
Ω
,where PWC
(
Ω
) denotes the subspace of L
∞
(
Ω
) made of piecewise constant functions defined in Ω.One also assume thatf and g satisfy the following conditions.(i)
Compatibility condition is(6)
f
∈
L
2
Ω
,
g
∈
L
2
Γ
,
such  that
∫
Ω
f
x
d
x
-
∫
Γ
g
x
d
γ
x
=
0
.
(ii)
Null average condition is(7)
∫
Ω
φ
x
d
x
=
0
.
(iii)
Let us emphasize that the novelty of this work is that Neumann boundary conditions are imposed on the whole boundary, which causes additional difficulties in the analysis.


## 2. A Finite Volume Formulation of the Model Problem

We briefly present finite volume formulation of the model problem (1)-(2) on unstructured meshes into the same spirit as [20, 25, 29–31]. We assume thatthe diffusion matrix
K
is a piecewise constant function over
Ω.This assumption is currently used at least in industrial problems (e.g., reservoir simulation problems). Indeed, a subsurface area is made up of a collection of various geologic formations that may be characterized at intermediate scales by averaged full permeability tensors over grid blocks of the primary grid: for more details on this topic, see [32, 33].


### 2.1. Formulation of the Discrete Problem

First of all, notice that the method under consideration is analyzed in this work for general polygonal domains covered with unstructured matching primary meshesP made up of arbitrary convex polygons (see Figure 1). Let us consider some definitions needed in what follows.Figure 1
Example of a primary unstructured conforming mesh.Definition 1.
A mesh defined onΩ is compatible with the discontinuities of the permeability tensor K if these discontinuities are located along the mesh interfaces.Definition 2.
One defines an edge-point as any point located over an edge and different from the extremities of that edge.Definition 3.
Two edge-pointsI and J are named “neighboring edge-points” if they share the same vertex V in the sense that I and J belong to two different edges that intersect in V.Let us recall our main objectives in this work:(i)
Compute at the cell-points (to be defined later) and at the interior vertices from the meshP the values of the unknown function φ as a solution (expected unique) of a linear system.
(ii)
Analyze the stability and the convergence of this solution in some discrete energy norm similar to the one introduced in [31].In the context of unstructured primary meshes (including square primary meshes with cell-points chosen different from cell-centers), the definition of a discrete energy norm requires that the cell-points lie inside certain perimeters. In this connection, the main steps for defining the cell-points are as follows:(a)
Choose arbitrarily a unique point (different from a vertex) on each edge of the meshP. This process generates a finite family of edge-points.
(b)
Join every pair of neighboring edge-points by a dotted straight line.
By this way, we generate an auxiliary mesh denotedA.
(c)
Fix arbitrarily a unique point inside each intersection of a primary cell and an auxiliary cell. These points define a finite family of cell-points.Figure 2 illustrates the location of the edge-points and cell-points within both primary and auxiliary meshes.Figure 2
A primary mesh (full lines) and the associated auxiliary mesh (dotted lines), including edge-points and cell-points, respectively, in black and blue colors.Remark 4.
Note that, in the 3D framework, for a given primary mesh the corresponding auxiliary mesh is generated very easily as follows. Any primary cellC involves a finite number of faces and to each face is assigned one and only one face point lying necessarily on the boundary of C. Therefore, one could associate with these face points the smallest convex polygon containing all of them. This is the process leading to generation of the auxiliary mesh associated with the primary mesh under consideration.We should emphasize that we need to take cell-points inside auxiliary cells (see Figure2) in view of achieving the following goal which is to define a vector space for the piecewise constant solution (named “weak approximate” solution in the sequel) and to equip it with a discrete energy norm.Remark 5.
Note that there exists a DDFV approach based on pressure gradient reconstructions that can address nonlinear elliptic problems: see [25, 34] for more details. Nevertheless in these works, the choice of edge-points depends strongly on that of cell-points as each edge-point coincides (by Definition 2) with the intersection of a primary edge and the straight line joining two cell-points located in both sides of that edge.From the boundary-value problem theory (see, e.g., [35]), system (1)-(2) possesses a unique solution in H
1
Ω under assumptions (3)–(7). We have assumed that the diffusion coefficient K is a piecewise constant function over Ω. The discontinuities of K naturally divide Ω into a finite number of convex subdomains Ω
s
s
∈
S. We now make the additional assumption that the restriction over Ω
s of the exact solution to system (1)-(2), denoted by φ
Ω
s, satisfies (8)
φ
Ω
s
∈
C
2
Ω
-
s
∀
s
∈
S
.Let us now focus on a finite volume formulation of problem (1)-(2) in terms of a linear system which should be derived from the elimination of auxiliary unknowns, namely, edge-point pressures, in flux balance equations over primary cells and also dual cells (to be introduced later). This linear system should involve the real numbers u
P
P
∈
P and u
D
∗
D
∗
∈
D as discrete unknowns which are expected to be reasonable approximations of φ
P
P
∈
P (cell-point pressures) and φ
D
∗
D
∗
∈
D (interior vertex pressures), respectively, where φ
P
=
φ
(
x
1
P
,
x
2
P
) and φ
D
∗
=
φ
(
x
1
D
∗
,
x
2
D
∗
) and where D stands for dual mesh. We now give a description of the procedure leading to the linear discrete system.LetC
P be a primary cell, where P is the corresponding cell-point. We integrate the two sides of the mass balance equation (1) in C
P. Applying Ostrogradsky’s theorem to the integral in the left-hand side leads to computing the flux across the boundary of C
P. Thanks to a suitable quadrature formula this computation yields a relation involving edge-point pressures. Due to the flux continuity over the mesh interfaces, the edge-point pressures can be eliminated from the previous relation.As an illustration of this technique for computing the fluxes across primary cell boundaries, we consider the edgeA
∗
B
∗ associated with the primary cell C
P (see Figure 3).Figure 3
Two molecules for a finite volume computation of the flux across the edgeA
∗
B
∗. (a) The edge A
∗
B
∗ is lying inside the domain Ω. (b) The edge A
∗
B
∗ is part of the boundary of Ω.
(a)
(b)LetK
P be the diffusion tensor of the primary cell C
P
. Denoting ξ
A
∗
B
∗
P as the unit normal vector to A
∗
B
∗ exterior to C
P, the flux expression over the edge A
∗
B
∗ viewed as part of the boundary of C
P is given by(9)
Q
A
∗
B
∗
P
=
a
h
K
P
a
^
h
K
L
h
A
∗
B
∗
a
h
K
P
h
I
L
+
a
^
h
K
L
h
P
I
φ
P
-
φ
L
+
T
I
A
∗
,
B
∗
P
+
a
^
h
K
L
b
h
K
P
h
P
I
+
a
h
K
P
b
^
h
K
L
h
I
L
a
h
K
P
h
I
L
+
a
^
h
K
L
h
P
I
φ
B
∗
-
φ
A
∗
,where T
I
A
∗
,
B
∗
P is the truncation error and where we introduce the following necessary ingredients:(10)
σ
P
=
P
I
→
P
I
→
,
τ
h
=
A
∗
B
∗
→
A
∗
B
∗
→
,
h
=
max
⁡
size
P
,
size
D
,
size
A
,
a
h
K
P
=
1
cos
⁡
θ
h
P
,
I
ξ
A
∗
B
∗
P
t
K
P
ξ
A
∗
B
∗
P
,
b
h
K
P
=
1
cos
⁡
θ
h
P
,
I
ξ
P
I
B
∗
t
K
P
ξ
A
∗
B
∗
P
,
a
^
h
K
L
=
1
cos
⁡
θ
h
L
,
I
ξ
A
∗
B
∗
P
t
K
L
ξ
A
∗
B
∗
P
,
b
^
h
K
P
=
1
cos
⁡
θ
h
L
,
I
ξ
I
L
B
∗
t
K
L
ξ
A
∗
B
∗
P
,where θ
h
P
,
I (resp., θ
h
L
,
I) is the angle defined by the vectors σ
P and ξ
A
∗
B
∗
P (resp., -
σ
P and ξ
A
∗
B
∗
P) and where ξ
P
I
B
∗ (resp., ξ
I
L
B
∗) denotes the unit normal vector to (
P
I
) (resp., to [
I
L
]) exterior to the half-plane from R
2 containing the point B
∗ and bordered by the straight line (
P
I
) (resp., (
I
L
)). Notice that 0
≤
θ
h
P
,
I, θ
h
L
,
I
<
π
/
2. Moreover if the primary mesh (P
,
E) is regular in the sense of Definition 6, we have 0
≤
θ
h
L
,
I
<
π
/
2
-
θ and therefore 0
<
cos
⁡
(
π
/
2
-
θ
)
<
cos
⁡
(
θ
h
L
,
I
)
≤
1, where θ
∈
0
,
π
/
2 is a certain real number not depending on h.In what follows, we will need the following notations. We denote byE the set of all edge-points (from the primary mesh of course), E
i
n
t the subset of E made up of internal edge-points, that is, edge-points lying on Ω, E
e
x
t the subset of E made up of edge-points lying on edges included in the boundary of Ω, and E
P (for P
∈
P) the subset of E made up of edge-points lying on the boundary of the primary cell C
P.We now introduce the notion of regular primary mesh that should play a central role in the sequel.Definition 6.
The setP
,
E defines a regular mesh system if the following conditions are fulfilled.
There exist0
<
ϖ
≤
1 and 0
<
θ
≤
π
/
2, both of them mesh independent, such that(11)
d
P
≤
diam
⁡
P
≤
1
ϖ
d
P
∀
P
∈
P
,
(12)
0
≤
θ
h
P
,
I
<
π
2
-
θ
∀
P
∈
P
∀
I
∈
E
P
,where P stands for the set of primary cells and where d
P is the distance between a cell-point or a vertex from P and an edge-point from P.Proposition 7.
Assume that (i) the primary meshP is compatible with the discontinuities of the permeability tensor K (see Definition 1), (ii) the set P
,
E is a regular mesh system, and (iii) relation (8) is honored.
Then, there exists a strictly positive numberC mesh independent such that(13)
T
I
A
∗
,
B
∗
P
≤
C
h
2
.Using the previous notations and thanks to the consistency of the flux approximation across cell edges (see Proposition7), one reasonably can approximate the flux balance within any primary cell C
P as follows:(14)
∑
I
∈
E
P
∩
E
i
n
t
a
h
K
P
a
^
h
K
L
h
A
∗
B
∗
a
h
K
P
h
I
L
+
a
^
h
K
L
h
P
I
φ
P
-
φ
L
+
a
^
h
K
L
b
h
K
P
h
P
I
+
a
h
K
P
b
^
h
K
L
h
I
L
a
h
K
P
h
I
L
+
a
^
h
K
L
h
P
I
φ
B
∗
-
φ
A
∗
≈
∫
C
P
f
x
d
x
-
∫
Γ
P
∩
Γ
g
x
d
γ
x
=
0
∀
P
∈
P
,where L
∈
P, A
∗
∈
D, and B
∗
∈
D are such that Γ
P
∩
Γ
L
=
[
A
∗
B
∗
] and where E
P
∩
E
L
=
I. It is clear that the number of discrete unknowns φ
P
P
∈
P and φ
D
∗
D
∗
∈
D is greater than the number of equations in system (14). It is then natural to complete this system with discrete equations obtained from mass flux balance over dual cells; see Figure 4 for the definition of the dual mesh.Figure 4
A primary mesh in full lines with the associated dual mesh in red discontinuous lines and the corresponding auxiliary mesh in black dotted lines.In what follows, we will need to use the notion of pseudoedge associated with dual cells. Let us define now this notion that is illustrated in Figure5.Figure 5
A dual cell (blue discontinuous line) with its four pseudoedges “centered” on red edge-points.Definition 8.
LetP and L be two cell-points from the primary mesh (i.e., P
,
L
∈
P) such that the corresponding primary cells C
P and C
L are adjacent, and consider I
∈
E
P
∩
E
L (recall that E
P, for E
∈
P, is the set of edge-points lying in the boundary of the primary cell C
E). The line [
P
I
]
∪
[
I
L
] defines a pseudoedge denoted by [
P
I
L
] and “centered” on I, with extremities P and L.We will say that a pseudoedge is associated with a dual cellC
A
∗ if it is lying in the boundary of C
A
∗.Remark 9.
The boundary of each dual cell is a union of a finite number of pseudoedges (see Figure4).Let us now look for discrete flux balance equations over dual cells. Integrating the two sides of the balance equation (1) in a dual cell C
B
∗, applying Ostrogradsky’s theorem for the left-hand side, and exploiting Remark 9 lead to(15)
∑
I
∈
E
B
∗
-
∫
P
I
L
K
grad
⁡
φ
·
n
B
∗
d
γ
=
∫
C
B
∗
f
x
d
x
,where n
B
∗ stands for the outward unit normal vector to the boundary of C
B
∗ and where [
P
I
L
] is a pseudoedge associated with the dual cell C
B
∗. Recall that E
B
∗ is the set of edge-points lying in the boundary of the dual cell C
B
∗.Let us look for a flux approximation across the pseudoedge[
P
I
L
] viewed as part of the boundary of C
B
∗. Denoting by Q
[
P
I
L
]
B
∗ the exact flux over [
P
I
L
], it can be expressed by the relation(16)
Q
P
I
L
B
∗
=
-
∫
P
I
grad
⁡
φ
·
K
P
ξ
P
I
B
∗
d
γ
-
∫
I
L
grad
⁡
φ
·
K
L
ξ
I
L
B
∗
d
γ
.By using the same process, the computation of the flux across[
P
I
] is(17)
-
∫
P
I
grad
⁡
φ
·
K
P
ξ
P
I
B
∗
d
γ
=
c
h
K
P
a
^
h
K
L
h
P
I
a
h
K
P
h
I
L
+
a
^
h
K
L
h
P
I
φ
P
-
φ
L
+
T
P
I
B
∗
+
h
P
I
d
h
K
P
a
h
K
P
h
I
L
+
a
^
h
K
L
h
P
I
+
c
h
K
P
h
I
L
b
^
h
K
L
-
b
h
K
P
h
A
∗
B
∗
a
h
K
P
h
I
L
+
a
^
h
K
L
h
P
I
φ
B
∗
-
φ
A
∗
,where (18)
c
h
K
P
=
1
cos
⁡
θ
h
P
,
I
ξ
A
∗
B
∗
P
t
K
P
ξ
P
I
B
∗
,
d
h
K
P
=
1
cos
⁡
θ
h
P
,
I
ξ
P
I
B
∗
t
K
P
ξ
P
I
B
∗(recall that θ
h
P
,
I denotes the angle defined by the vectors σ
P and ξ
A
∗
B
∗
P).Similarly, the computation of the flux across[
I
L
] leads to(19)
-
∫
I
L
grad
⁡
φ
·
K
L
ξ
I
L
B
∗
d
γ
=
c
^
h
K
L
a
h
K
P
h
I
L
a
h
K
P
h
I
L
+
a
^
h
K
L
h
P
I
φ
P
-
φ
L
+
T
I
L
B
∗
+
h
I
L
d
^
h
K
L
a
h
K
P
h
I
L
+
a
^
h
K
L
h
P
I
+
c
^
h
K
P
h
P
I
b
h
K
P
-
b
^
h
K
L
h
A
∗
B
∗
a
h
K
P
h
I
L
+
a
^
h
K
L
h
P
I
φ
B
∗
-
φ
A
∗
,where (20)
c
^
h
K
P
=
1
cos
⁡
θ
h
L
,
I
ξ
A
∗
B
∗
P
t
K
L
ξ
I
L
B
∗
,
d
^
h
K
L
=
1
cos
⁡
θ
h
L
,
I
ξ
I
L
B
∗
t
K
L
ξ
I
L
B
∗
.We should now formulate a global estimate for the truncation errors associated with the flux approximation over the pseudoedge[
P
I
L
].Proposition 10.
Under the same assumptions as those of Proposition7, there exists a positive number C mesh independent such that(21)
T
P
I
B
∗
+
T
I
L
B
∗
≤
C
h
2
.The above inequality shows the consistency of the flux approximation across the pseudoedge[
P
I
L
]. We summarize what precedes as(22)
Q
P
I
L
B
∗
≈
c
h
K
P
a
^
h
K
L
h
P
I
+
c
^
h
K
L
a
h
K
P
h
I
L
a
h
K
P
h
I
L
+
a
^
h
K
L
h
P
I
φ
P
-
φ
L
+
ω
h
P
,
L
,
I
h
A
∗
B
∗
a
h
K
P
h
I
L
+
a
^
h
K
L
h
P
I
φ
B
∗
-
φ
A
∗
,where we have set(23)
ω
h
P
,
L
,
I
=
a
h
K
P
h
I
L
+
a
^
h
K
L
h
P
I
d
h
K
P
h
P
I
+
d
^
h
K
L
h
I
L
+
h
P
I
h
I
L
b
^
h
K
L
-
b
h
K
P
c
h
K
P
-
c
^
h
K
L
.We deduce from (22) that an approximate flux balance equation over any dual cell C
B
∗ can be formulated as follows:(24)
∑
I
∈
E
P
∩
E
i
n
t
c
h
K
P
a
^
h
K
L
h
P
I
+
c
^
h
K
L
a
h
K
P
h
I
L
a
h
K
P
h
I
L
+
a
^
h
K
L
h
P
I
φ
P
-
φ
L
+
ω
h
P
,
L
,
I
h
A
∗
B
∗
a
h
K
P
h
I
L
+
a
^
h
K
L
h
P
I
φ
B
∗
-
φ
A
∗
+
∑
I
∈
E
B
∗
∩
E
e
x
t
h
I
E
a
K
E
d
h
K
E
-
c
h
K
E
b
h
K
E
a
h
K
E
h
B
∗
D
∗
φ
B
∗
-
φ
D
∗
≈
∫
C
B
∗
f
x
d
x
-
∫
Γ
B
∗
∩
Γ
g
d
γ
+
∑
I
∈
E
B
∗
∩
E
e
x
t
c
h
K
E
h
I
E
a
K
E
h
B
∗
D
∗
∫
B
∗
D
∗
g
d
γ
∀
B
∗
∈
D
,where P
∈
P, L
∈
P, and A
∗
∈
D are such that Γ
P
∩
Γ
L
=
A
∗
B
∗ and P
L
∩
A
∗
B
∗
=
I and where E
∈
P, D
∗
∈
D a boundary-vertex are such that B
∗
D
∗
⊂
Γ
∩
Γ
E and I
∈
B
∗
D
∗.Systems (14) and (24) naturally suggest defining a finite volume formulation of the diffusion problem (1)-(2) as follows.Findφ
-
P
P
∈
P and φ
-
D
∗
D
∗
∈
D such that(25)
∑
I
∈
E
P
∩
E
i
n
t
a
h
K
P
a
^
h
K
L
h
A
∗
B
∗
a
h
K
P
h
I
L
+
a
^
h
K
L
h
P
I
φ
-
P
-
φ
-
L
+
a
^
h
K
L
b
h
K
P
h
P
I
+
a
h
K
P
b
^
h
K
L
h
I
L
a
h
K
P
h
I
L
+
a
^
h
K
L
h
P
I
φ
-
B
∗
-
φ
-
A
∗
=
∫
C
P
f
x
d
x
-
∫
Γ
P
∩
Γ
g
d
γ
=
0
∀
P
∈
P
,
(26)
∑
I
∈
E
B
∗
∩
E
i
n
t
c
h
K
P
a
^
h
K
L
h
P
I
+
c
^
h
K
L
a
h
K
P
h
I
L
a
h
K
P
h
I
L
+
a
^
h
K
L
h
P
I
φ
-
P
-
φ
-
L
+
ω
h
P
,
L
,
I
h
A
∗
B
∗
a
h
K
P
h
I
L
+
a
^
h
K
L
h
P
I
φ
-
B
∗
-
φ
-
A
∗
+
∑
I
∈
E
B
∗
∩
E
e
x
t
h
I
E
a
K
E
d
h
K
E
-
c
h
K
E
b
h
K
E
a
h
K
E
h
B
∗
D
∗
φ
-
B
∗
-
φ
-
D
∗
=
∫
C
B
∗
f
x
d
x
-
∫
C
B
∗
∩
Γ
g
d
γ
+
∑
I
∈
E
B
∗
∩
E
e
x
t
c
h
K
E
h
I
E
a
K
E
h
B
∗
D
∗
∫
B
∗
D
∗
g
d
γ
∀
B
∗
∈
D
.Remark 11.
It is useful to note that(27)
∑
B
∗
∈
D
∑
I
∈
E
B
∗
∩
E
e
x
t
c
h
K
E
h
I
E
a
K
E
h
B
∗
D
∗
∫
B
∗
D
∗
g
d
γ
=
0
.


### 2.2. Existence of Solutions to System (25)-(26) Conditions for Uniqueness

(i) Matrix Properties of the DDFV Problem (25)-(26). It is easily seen that the symmetry of the matrix associated with the linear system (25)-(26) essentially follows from the symmetry of the diffusion coefficient K (see assumption (3)). We assume that all the cell-points and all the vertices (with respect to the primary mesh) are numbered. Therefore one can identify P and D with two disjoint subsets of the set of positive integers. To fix ideas, let us set(28)
P
≡
1
,
…
,
n
,
D
≡
n
+
1
,
…
,
n
+
m
.Then C
a
r
d
(
P
) denotes the total number of cell-points and C
a
r
d
(
D
) stands for the total number of vertices. On the other hand, define the subvectors Φ
c
c and Φ
v
c by(29)
Φ
c
c
=
φ
-
P
P
∈
P
,
Φ
v
c
=
φ
-
D
∗
D
∗
∈
Dand set (30)
M
=
A
B
B
t
C
,where(31)
M
Φ
c
c
Φ
v
c
=
F
c
c
F
v
c
,where F
c
c is a subvector with C
a
r
d
(
P
) components defined by the right-hand side of (25) and F
v
c is a subvector with C
a
r
d
(
D
) components defined by the right-hand side of (26).Remark 12.
The matrixM satisfies the following properties:(1)
1
≤
i
,
j
≤
m
M
i
j
=
M
j
i;
(2)
1
≤
i
,
j
≤
m
∑
j
=
1
m
M
i
j
=
0;
(3)
1
≤
i
,
j
≤
m
∑
i
=
1
m
M
i
j
=
0.Let us introduce two vectors ofR
n
+
m named F
c
c and F
v
c and defined by(32)
F
c
c
i
=
1
if
1
≤
i
≤
n
0
if
n
+
1
≤
i
≤
n
+
m
,
F
v
c
i
=
0
if
1
≤
i
≤
n
1
if
n
+
1
≤
i
≤
n
+
m
.Proposition 13 (characterization of Kernel space ofM).
(i)
The matrixM is singular.
(ii)
Moreover, letker
⁡
(
M
) denote the subset of R
n
+
m defined as follows:
⁡
ker
⁡
M
=
V
∈
R
n
+
m
,
M
V
=
0 (named Kernel space of M in the sequel); then we have(33)
ker
⁡
M
=
F
c
c
,
F
v
c
.Sketch for the Proof
(i)
The singularity ofM is an immediate consequence of Remark 12.
(ii)
Define the spaceR
P
,
D by(34)
R
P
,
D
=
V
=
V
c
c
,
V
v
c
;
V
c
c
=
V
P
P
∈
P
⊂
R
,
V
v
c
=
V
D
∗
D
∗
∈
D
⊂
Rand endow it with the seminorm ·
R
P
,
D defined as(35)
V
R
P
,
D
=
∑
P
,
L
∈
P
2
,
A
∗
,
B
∗
∈
D
2
s.t.
Γ
P
∩
Γ
L
=
A
∗
B
∗
V
P
-
V
L
2
+
V
B
∗
-
V
A
∗
2
+
∑
A
∗
,
B
∗
∈
D
2
s.t.
A
∗
B
∗
∈
E
e
x
t
V
B
∗
-
V
A
∗
2
1
/
2
∀
V
∈
R
P
,
D
.Then, find the above identification of the Kernel space of M using the following Lemma.Lemma 14.
The matrixM is positive; that is,(36)
V
c
c
,
V
v
c
M
V
c
c
V
v
c
≥
0
∀
V
∈
R
P
,
D
.Proof.
LetV
=
V
c
c
,
V
v
c
∈
R
P
,
D, where V
c
c
=
V
P
P
∈
P and V
v
c
=
V
A
∗
A
∗
∈
D.
It follows from what precedes that(37)
V
c
c
,
V
v
c
M
V
c
c
V
v
c
=
∑
P
∈
P
∑
I
∈
E
P
∩
E
i
n
t
a
h
K
P
a
^
h
K
L
h
A
∗
B
∗
a
h
K
P
h
I
L
+
a
^
h
K
L
h
P
I
V
P
-
V
L
V
P
+
a
^
h
K
L
b
h
K
P
h
P
I
+
a
h
K
P
b
^
h
K
L
h
I
L
a
h
K
P
h
I
L
+
a
^
h
K
L
h
P
I
V
B
∗
-
V
A
∗
V
P
+
∑
B
∗
∈
D
∑
I
∈
E
B
∗
∩
E
i
n
t
c
h
K
P
a
^
h
K
L
h
P
I
+
c
^
h
K
L
a
h
K
P
h
I
L
a
h
K
P
h
I
L
+
a
^
h
K
L
h
P
I
V
P
-
V
L
V
B
∗
+
ω
h
P
,
L
,
I
h
A
∗
B
∗
a
h
K
P
h
I
L
+
a
^
h
K
L
h
P
I
V
B
∗
-
V
A
∗
V
B
∗
+
∑
B
∗
∈
D
∑
I
∈
E
B
∗
∩
E
e
x
t
h
I
E
a
K
E
d
h
K
E
-
c
h
K
E
b
h
K
E
a
h
K
E
h
B
∗
D
∗
V
B
∗
-
V
D
∗
V
B
∗
.Define(38)
K
11
P
L
=
a
h
K
P
a
^
h
K
L
h
A
∗
B
∗
a
h
K
P
h
I
L
+
a
^
h
K
L
h
P
I
,
K
22
P
L
=
ω
h
P
,
L
,
I
h
A
∗
B
∗
a
h
K
P
h
I
L
+
a
^
h
K
L
h
P
I
,
Π
P
I
,
A
∗
,
B
∗
=
a
K
P
d
h
K
P
-
c
h
K
P
b
h
K
P
h
I
P
a
h
K
P
h
A
∗
B
∗
,where P
,
L
∈
P are two adjacent primary cells sharing the edge A
∗
B
∗ as interface containing the edge-point I. Then, relation (37) becomes(39)
V
c
c
,
V
v
c
M
V
c
c
V
v
c
=
∑
P
,
L
∈
P
2
,
A
∗
,
B
∗
∈
D
2
s.t.
Γ
P
∩
Γ
L
=
A
∗
B
∗
K
11
P
L
V
P
-
V
L
2
+
K
22
P
L
V
B
∗
-
V
A
∗
2
+
2
K
12
P
L
V
P
-
V
L
V
B
∗
-
V
A
∗
+
∑
P
∈
P
e
x
t
,
A
∗
,
B
∗
∈
D
2
∩
Γ
2
s.t.
A
∗
B
∗
∈
E
P
Π
P
I
,
A
∗
,
B
∗
V
B
∗
-
V
A
∗
2
,where P
e
x
t denotes the subset of P made of primary cells adjacent to the domain boundary Γ. Let us prove that the homogenized symmetric matrix K
P
L is positive definite; that is,(40)
K
11
P
L
K
22
P
L
-
K
12
P
L
2
>
0
.Setting(41)
Δ
P
L
=
K
11
P
L
K
22
P
L
-
K
12
P
L
2
,it is easy to check that(42)
Δ
P
L
=
N
1
a
h
K
P
d
h
K
P
-
b
h
K
P
2
+
N
2
a
^
h
K
L
d
^
h
K
L
-
b
^
h
K
L
2
,where we have set(43)
N
1
=
a
^
h
K
L
h
I
L
a
h
K
P
h
I
L
+
a
^
h
K
L
h
P
I
2
+
a
h
K
P
a
^
h
K
L
h
P
I
h
I
L
a
h
K
P
h
I
L
+
a
^
h
K
L
h
P
I
2
,
N
2
=
a
h
K
P
h
I
L
a
h
K
P
h
I
L
+
a
^
h
K
L
h
P
I
2
+
a
h
K
P
a
^
h
K
L
h
P
I
h
I
L
a
h
K
P
h
I
L
+
a
^
h
K
L
h
P
I
2which are strictly positive numbers.
Since the diffusion matrixK is symmetric and positive definite (see assumptions (3)-(4)), Cauchy-Schwarz’s inequality for the inner product associated with K ensures that(44)
a
h
K
P
d
h
K
P
-
b
h
K
P
2
>
0
,
a
^
h
K
L
d
^
h
K
L
-
b
^
h
K
L
2
>
0as either ξ
A
∗
B
∗
P and ξ
P
I
B
∗ or ξ
A
∗
B
∗
P and ξ
I
L
B
∗ are not collinear. Therefore, Δ
P
L
>
0 and thus K
P
L is a symmetric and positive definite matrix.
It follows from what precedes that the matrixK
P
L possesses strictly positive eigenvalues. Let λ
m
i
n
P
L be its least eigenvalue. So we have(45)
V
c
c
,
V
v
c
M
V
c
c
V
v
c
≥
∑
P
,
L
∈
P
2
,
A
∗
,
B
∗
∈
D
2
s.t.
Γ
P
∩
Γ
L
=
A
∗
B
∗
λ
min
⁡
P
L
V
P
-
V
L
2
+
V
B
∗
-
V
A
∗
2
+
∑
P
∈
P
e
x
t
,
A
∗
,
B
∗
∈
D
2
∩
Γ
2
s.t.
A
∗
B
∗
∈
E
P
Π
P
I
,
A
∗
,
B
∗
V
B
∗
-
V
A
∗
2
.Remarking that there exist two real numbers δ and η strictly positive depending exclusively on the geological structure of the domain such that(46)
δ
≤
K
11
P
K
22
P
-
K
12
P
2
≤
η
∀
P
∈
P
,the following result is easily seen.
Lemma  15.(47)
Π
P
I
,
A
∗
,
B
∗
=
K
11
P
K
22
P
-
K
12
P
2
h
P
I
a
h
K
P
h
A
∗
B
∗
.
Moreover
(48)
δ
ϖ
sin
⁡
θ
γ
m
a
x
<
Π
P
I
,
A
∗
,
B
∗
<
η
ϖ
γ
m
i
n
.
Thanks to inequality (45) and to Lemma 15 the proof of Lemma 14 ends.Lemma 16.
The matrixM satisfies the relation(49)
γ
V
R
P
,
D
2
≤
V
c
c
,
V
v
c
M
V
c
c
V
v
c
∀
V
∈
R
P
,
D
,where γ is a strictly positive number mesh independent and where ·
R
P
,
D is a seminorm defined on R
P
,
D by relation (35).Proof.
It is conducted with the same arguments as those developed in the previous proof, except that one should go much farther by proving that there existsρ
>
0 mesh independent such that(50)
ρ
≤
λ
m
i
n
P
L
∀
P
,
L
∈
P
.The eigenvalues λ of the symmetric positive definite matrix K
P
L satisfy the so-called characteristic equation associated with K
P
L; that is,(51)
λ
2
-
K
11
P
L
+
K
22
P
L
λ
+
K
11
P
L
K
22
P
L
-
K
12
P
L
2
=
0
.The least eigenvalue of K
P
L denoted by(52)
λ
min
⁡
P
L
=
K
11
P
L
+
K
22
P
L
-
Δ
2
,where Δ
=
K
11
P
L
+
K
22
P
L
2
-
4
K
11
P
L
K
22
P
L
-
K
12
P
L
2, is a strictly positive number. One can easily deduce that(53)
λ
min
⁡
P
L
≥
det
⁡
K
P
L
K
11
P
L
+
K
22
P
L
+
det
⁡
K
P
L
+
1
,where det
⁡
K
P
L
=
K
11
P
L
K
22
P
L
-
K
12
P
L
2 is a strictly positive number. We should bound the quantities K
11
P
L, K
22
P
L, and det
⁡
K
P
L by mesh independent strictly positive numbers. Let us start first with det
⁡
K
P
L
. We consider a change of coordinates by moving from the initial Cartesian coordinates to a local one, namely, J
,
C
∗
D
∗
→
/
C
∗
D
∗
→
,
ξ
C
∗
D
∗
⊥, where J is the edge-point on the interface C
∗
D
∗ between the cells C
G and C
H and where ξ
C
∗
D
∗
⊥ is a vector orthogonal to C
∗
D
∗ and oriented such that the basis change matrix M is a rotation. Denoting the permeability tensor of the cell C
G by K
G
=
K
i
j
G in the initial Cartesian coordinates and by K
~
G
=
K
~
i
j
G in the local coordinates we have(54)
K
~
G
=
M
-
1
K
G
M
.Similarly, we get for the cell C
H
(55)
K
~
H
=
M
-
1
K
H
M
.Then it is easy to check that(56)
a
h
K
G
d
h
K
G
-
b
h
K
G
2
=
K
~
11
G
K
~
22
G
-
K
~
12
G
2
;that is,(57)
a
h
K
G
d
h
K
G
-
b
h
K
G
2
=
det
⁡
K
~
G
=
det
⁡
K
G
,where det
⁡
· denotes the determinant. Similarly, we have for the cell C
H
(58)
a
^
h
K
H
d
^
h
K
H
-
b
^
h
K
H
2
=
det
⁡
K
~
H
=
det
⁡
K
H
.It follows from what precedes that(59)
det
⁡
K
P
L
=
N
1
det
⁡
K
G
+
N
2
det
⁡
K
H
,where N
1 and N
2 are given, respectively, by relations (43).
On one hand, we can deduce from (4) and Definition 6 that(60)
N
i
≥
1
2
ϖ
γ
min
⁡
sin
⁡
θ
γ
max
⁡
2
∀
i
=
1,2
.On the other hand, we remark that(61)
det
⁡
K
P
≥
min
⁡
det
⁡
K
s
,
s
∈
S
∀
P
∈
P
,where the set S (introduced in (8)) depends exclusively on the lithologic structure of the medium Ω. Then we deduce from (59), (60), and (61) that(62)
det
⁡
K
P
≥
1
2
ϖ
γ
min
⁡
sin
⁡
θ
γ
max
⁡
2
min
⁡
det
⁡
K
s
,
s
∈
S
∀
P
∈
P
.Remarking that(63)
det
⁡
K
P
L
≤
K
11
P
L
K
22
P
L
+
K
12
P
L
2and exploiting again (4) and Definition 6 lead to the following inequality:(64)
det
⁡
K
P
L
≤
2
γ
max
⁡
4
max
⁡
2
,
ϖ
ϖ
3
γ
min
⁡
2
sin
⁡
θ
4
.Thanks again to (4) and Definition 6 one can easily check that(65)
K
11
P
L
+
K
22
P
L
≤
γ
max
⁡
2
ϖ
sin
⁡
θ
2
1
+
2
ϖ
+
4
γ
max
⁡
ϖ
sin
⁡
θ
.Lemma 16 follows from the combination of (53), (62), (64), and (65).Proposition 17 (discrete compatibility condition).
The right-hand side of the discrete system (25)-(26) satisfies the following discrete compatibility condition:(66)
∑
P
∈
P
∫
C
P
f
x
d
x
-
∫
Γ
P
∩
Γ
g
d
γ
=
0
,
∑
B
∗
∈
D
∫
C
B
∗
f
x
d
x
-
∫
Γ
B
∗
∩
Γ
g
d
γ
+
∑
B
∗
∈
D
∑
I
∈
E
B
∗
∩
E
ext
c
h
K
E
h
I
E
a
K
E
h
B
∗
D
∗
∫
B
∗
D
∗
g
d
γ
=
0
.Proof.
First of all, note that the double summation in the previous proposition is equal to zero thanks to Remark11. Let us consider two vectors of R
n
+
m called F
c
c and F
v
c and defined by(67)
F
c
c
i
=
1
if
1
≤
i
≤
n
0
if
n
+
1
≤
i
≤
n
+
m
,
F
v
c
i
=
0
if
1
≤
i
≤
n
1
if
n
+
1
≤
i
≤
n
+
m
.According to the matrix form of the discrete system (see relation (31)), we have(68)
F
c
c
t
M
Φ
c
c
Φ
v
c
=
F
c
c
t
F
c
c
F
v
c
.It follows from Remark 12 that(69)
F
c
c
t
M
Φ
c
c
Φ
v
c
=
Φ
c
c
Φ
v
c
t
M
F
c
c
=
0
.Therefore(70)
0
=
F
c
c
t
F
c
c
F
v
c
=
∑
P
∈
P
∫
C
P
f
x
d
x
-
∫
Γ
P
∩
Γ
g
d
γ
.Similarly we have(71)
0
=
F
v
c
t
M
Φ
c
c
Φ
v
c
=
F
v
c
t
F
c
c
F
v
c
=
∑
B
∗
∈
D
∫
C
B
∗
f
x
d
x
-
∫
C
B
∗
∩
Γ
g
d
γ
.This ends the proof of Proposition 17.An immediate consequence of the previous result is the following existence result for a solution to the discrete problem (25)-(26).Proposition 18 (existence result).
(
1) The linear system (25)-(26) possesses an infinite number of solutions.
(
2) More precisely, if s
h is a solution to the discrete system (25)-(26), then s
h
+
ker
⁡
M is the set of all solutions for this system.The question to know how to get the physical solution of the model problem has a natural answer; that is, one should follow the same way as the continuous problem analysis. Since the dimension ofker
⁡
M is 2, required are two (linearly independent) discrete versions of the null average condition (7). This naturally leads to the following result.Proposition 19 (uniqueness result).
Under the following (null average) conditions,(72)
i
∑
P
∈
P
m
e
s
C
P
φ
-
P
=
0
,
i
i
∑
D
∗
∈
D
m
e
s
C
D
∗
φ
-
D
∗
=
0
,the discrete problem (25)-(26) possesses a unique solution.Proof.
Let us consider the space(73)
R
P
,
D
=
V
=
V
c
c
,
V
v
c
;
V
c
c
=
V
P
P
∈
P
⊂
R
,
V
v
c
=
V
D
∗
D
∗
∈
D
⊂
R
.One knows from Lemma 14 that(74)
V
t
M
V
≥
0
∀
V
∈
R
P
,
D
.Then it follows that, for all V in the subspace of R
P
,
D made up of V
=
V
P
P
∈
P
,
V
D
∗
D
∗
∈
D such that conditions (72) are fulfilled, one has V
t
M
V
=
0
⇔
V
=
0.


## 2.1. Formulation of the Discrete Problem

First of all, notice that the method under consideration is analyzed in this work for general polygonal domains covered with unstructured matching primary meshesP made up of arbitrary convex polygons (see Figure 1). Let us consider some definitions needed in what follows.Figure 1
Example of a primary unstructured conforming mesh.Definition 1.
A mesh defined onΩ is compatible with the discontinuities of the permeability tensor K if these discontinuities are located along the mesh interfaces.Definition 2.
One defines an edge-point as any point located over an edge and different from the extremities of that edge.Definition 3.
Two edge-pointsI and J are named “neighboring edge-points” if they share the same vertex V in the sense that I and J belong to two different edges that intersect in V.Let us recall our main objectives in this work:(i)
Compute at the cell-points (to be defined later) and at the interior vertices from the meshP the values of the unknown function φ as a solution (expected unique) of a linear system.
(ii)
Analyze the stability and the convergence of this solution in some discrete energy norm similar to the one introduced in [31].In the context of unstructured primary meshes (including square primary meshes with cell-points chosen different from cell-centers), the definition of a discrete energy norm requires that the cell-points lie inside certain perimeters. In this connection, the main steps for defining the cell-points are as follows:(a)
Choose arbitrarily a unique point (different from a vertex) on each edge of the meshP. This process generates a finite family of edge-points.
(b)
Join every pair of neighboring edge-points by a dotted straight line.
By this way, we generate an auxiliary mesh denotedA.
(c)
Fix arbitrarily a unique point inside each intersection of a primary cell and an auxiliary cell. These points define a finite family of cell-points.Figure 2 illustrates the location of the edge-points and cell-points within both primary and auxiliary meshes.Figure 2
A primary mesh (full lines) and the associated auxiliary mesh (dotted lines), including edge-points and cell-points, respectively, in black and blue colors.Remark 4.
Note that, in the 3D framework, for a given primary mesh the corresponding auxiliary mesh is generated very easily as follows. Any primary cellC involves a finite number of faces and to each face is assigned one and only one face point lying necessarily on the boundary of C. Therefore, one could associate with these face points the smallest convex polygon containing all of them. This is the process leading to generation of the auxiliary mesh associated with the primary mesh under consideration.We should emphasize that we need to take cell-points inside auxiliary cells (see Figure2) in view of achieving the following goal which is to define a vector space for the piecewise constant solution (named “weak approximate” solution in the sequel) and to equip it with a discrete energy norm.Remark 5.
Note that there exists a DDFV approach based on pressure gradient reconstructions that can address nonlinear elliptic problems: see [25, 34] for more details. Nevertheless in these works, the choice of edge-points depends strongly on that of cell-points as each edge-point coincides (by Definition 2) with the intersection of a primary edge and the straight line joining two cell-points located in both sides of that edge.From the boundary-value problem theory (see, e.g., [35]), system (1)-(2) possesses a unique solution in H
1
Ω under assumptions (3)–(7). We have assumed that the diffusion coefficient K is a piecewise constant function over Ω. The discontinuities of K naturally divide Ω into a finite number of convex subdomains Ω
s
s
∈
S. We now make the additional assumption that the restriction over Ω
s of the exact solution to system (1)-(2), denoted by φ
Ω
s, satisfies (8)
φ
Ω
s
∈
C
2
Ω
-
s
∀
s
∈
S
.Let us now focus on a finite volume formulation of problem (1)-(2) in terms of a linear system which should be derived from the elimination of auxiliary unknowns, namely, edge-point pressures, in flux balance equations over primary cells and also dual cells (to be introduced later). This linear system should involve the real numbers u
P
P
∈
P and u
D
∗
D
∗
∈
D as discrete unknowns which are expected to be reasonable approximations of φ
P
P
∈
P (cell-point pressures) and φ
D
∗
D
∗
∈
D (interior vertex pressures), respectively, where φ
P
=
φ
(
x
1
P
,
x
2
P
) and φ
D
∗
=
φ
(
x
1
D
∗
,
x
2
D
∗
) and where D stands for dual mesh. We now give a description of the procedure leading to the linear discrete system.LetC
P be a primary cell, where P is the corresponding cell-point. We integrate the two sides of the mass balance equation (1) in C
P. Applying Ostrogradsky’s theorem to the integral in the left-hand side leads to computing the flux across the boundary of C
P. Thanks to a suitable quadrature formula this computation yields a relation involving edge-point pressures. Due to the flux continuity over the mesh interfaces, the edge-point pressures can be eliminated from the previous relation.As an illustration of this technique for computing the fluxes across primary cell boundaries, we consider the edgeA
∗
B
∗ associated with the primary cell C
P (see Figure 3).Figure 3
Two molecules for a finite volume computation of the flux across the edgeA
∗
B
∗. (a) The edge A
∗
B
∗ is lying inside the domain Ω. (b) The edge A
∗
B
∗ is part of the boundary of Ω.
(a)
(b)LetK
P be the diffusion tensor of the primary cell C
P
. Denoting ξ
A
∗
B
∗
P as the unit normal vector to A
∗
B
∗ exterior to C
P, the flux expression over the edge A
∗
B
∗ viewed as part of the boundary of C
P is given by(9)
Q
A
∗
B
∗
P
=
a
h
K
P
a
^
h
K
L
h
A
∗
B
∗
a
h
K
P
h
I
L
+
a
^
h
K
L
h
P
I
φ
P
-
φ
L
+
T
I
A
∗
,
B
∗
P
+
a
^
h
K
L
b
h
K
P
h
P
I
+
a
h
K
P
b
^
h
K
L
h
I
L
a
h
K
P
h
I
L
+
a
^
h
K
L
h
P
I
φ
B
∗
-
φ
A
∗
,where T
I
A
∗
,
B
∗
P is the truncation error and where we introduce the following necessary ingredients:(10)
σ
P
=
P
I
→
P
I
→
,
τ
h
=
A
∗
B
∗
→
A
∗
B
∗
→
,
h
=
max
⁡
size
P
,
size
D
,
size
A
,
a
h
K
P
=
1
cos
⁡
θ
h
P
,
I
ξ
A
∗
B
∗
P
t
K
P
ξ
A
∗
B
∗
P
,
b
h
K
P
=
1
cos
⁡
θ
h
P
,
I
ξ
P
I
B
∗
t
K
P
ξ
A
∗
B
∗
P
,
a
^
h
K
L
=
1
cos
⁡
θ
h
L
,
I
ξ
A
∗
B
∗
P
t
K
L
ξ
A
∗
B
∗
P
,
b
^
h
K
P
=
1
cos
⁡
θ
h
L
,
I
ξ
I
L
B
∗
t
K
L
ξ
A
∗
B
∗
P
,where θ
h
P
,
I (resp., θ
h
L
,
I) is the angle defined by the vectors σ
P and ξ
A
∗
B
∗
P (resp., -
σ
P and ξ
A
∗
B
∗
P) and where ξ
P
I
B
∗ (resp., ξ
I
L
B
∗) denotes the unit normal vector to (
P
I
) (resp., to [
I
L
]) exterior to the half-plane from R
2 containing the point B
∗ and bordered by the straight line (
P
I
) (resp., (
I
L
)). Notice that 0
≤
θ
h
P
,
I, θ
h
L
,
I
<
π
/
2. Moreover if the primary mesh (P
,
E) is regular in the sense of Definition 6, we have 0
≤
θ
h
L
,
I
<
π
/
2
-
θ and therefore 0
<
cos
⁡
(
π
/
2
-
θ
)
<
cos
⁡
(
θ
h
L
,
I
)
≤
1, where θ
∈
0
,
π
/
2 is a certain real number not depending on h.In what follows, we will need the following notations. We denote byE the set of all edge-points (from the primary mesh of course), E
i
n
t the subset of E made up of internal edge-points, that is, edge-points lying on Ω, E
e
x
t the subset of E made up of edge-points lying on edges included in the boundary of Ω, and E
P (for P
∈
P) the subset of E made up of edge-points lying on the boundary of the primary cell C
P.We now introduce the notion of regular primary mesh that should play a central role in the sequel.Definition 6.
The setP
,
E defines a regular mesh system if the following conditions are fulfilled.
There exist0
<
ϖ
≤
1 and 0
<
θ
≤
π
/
2, both of them mesh independent, such that(11)
d
P
≤
diam
⁡
P
≤
1
ϖ
d
P
∀
P
∈
P
,
(12)
0
≤
θ
h
P
,
I
<
π
2
-
θ
∀
P
∈
P
∀
I
∈
E
P
,where P stands for the set of primary cells and where d
P is the distance between a cell-point or a vertex from P and an edge-point from P.Proposition 7.
Assume that (i) the primary meshP is compatible with the discontinuities of the permeability tensor K (see Definition 1), (ii) the set P
,
E is a regular mesh system, and (iii) relation (8) is honored.
Then, there exists a strictly positive numberC mesh independent such that(13)
T
I
A
∗
,
B
∗
P
≤
C
h
2
.Using the previous notations and thanks to the consistency of the flux approximation across cell edges (see Proposition7), one reasonably can approximate the flux balance within any primary cell C
P as follows:(14)
∑
I
∈
E
P
∩
E
i
n
t
a
h
K
P
a
^
h
K
L
h
A
∗
B
∗
a
h
K
P
h
I
L
+
a
^
h
K
L
h
P
I
φ
P
-
φ
L
+
a
^
h
K
L
b
h
K
P
h
P
I
+
a
h
K
P
b
^
h
K
L
h
I
L
a
h
K
P
h
I
L
+
a
^
h
K
L
h
P
I
φ
B
∗
-
φ
A
∗
≈
∫
C
P
f
x
d
x
-
∫
Γ
P
∩
Γ
g
x
d
γ
x
=
0
∀
P
∈
P
,where L
∈
P, A
∗
∈
D, and B
∗
∈
D are such that Γ
P
∩
Γ
L
=
[
A
∗
B
∗
] and where E
P
∩
E
L
=
I. It is clear that the number of discrete unknowns φ
P
P
∈
P and φ
D
∗
D
∗
∈
D is greater than the number of equations in system (14). It is then natural to complete this system with discrete equations obtained from mass flux balance over dual cells; see Figure 4 for the definition of the dual mesh.Figure 4
A primary mesh in full lines with the associated dual mesh in red discontinuous lines and the corresponding auxiliary mesh in black dotted lines.In what follows, we will need to use the notion of pseudoedge associated with dual cells. Let us define now this notion that is illustrated in Figure5.Figure 5
A dual cell (blue discontinuous line) with its four pseudoedges “centered” on red edge-points.Definition 8.
LetP and L be two cell-points from the primary mesh (i.e., P
,
L
∈
P) such that the corresponding primary cells C
P and C
L are adjacent, and consider I
∈
E
P
∩
E
L (recall that E
P, for E
∈
P, is the set of edge-points lying in the boundary of the primary cell C
E). The line [
P
I
]
∪
[
I
L
] defines a pseudoedge denoted by [
P
I
L
] and “centered” on I, with extremities P and L.We will say that a pseudoedge is associated with a dual cellC
A
∗ if it is lying in the boundary of C
A
∗.Remark 9.
The boundary of each dual cell is a union of a finite number of pseudoedges (see Figure4).Let us now look for discrete flux balance equations over dual cells. Integrating the two sides of the balance equation (1) in a dual cell C
B
∗, applying Ostrogradsky’s theorem for the left-hand side, and exploiting Remark 9 lead to(15)
∑
I
∈
E
B
∗
-
∫
P
I
L
K
grad
⁡
φ
·
n
B
∗
d
γ
=
∫
C
B
∗
f
x
d
x
,where n
B
∗ stands for the outward unit normal vector to the boundary of C
B
∗ and where [
P
I
L
] is a pseudoedge associated with the dual cell C
B
∗. Recall that E
B
∗ is the set of edge-points lying in the boundary of the dual cell C
B
∗.Let us look for a flux approximation across the pseudoedge[
P
I
L
] viewed as part of the boundary of C
B
∗. Denoting by Q
[
P
I
L
]
B
∗ the exact flux over [
P
I
L
], it can be expressed by the relation(16)
Q
P
I
L
B
∗
=
-
∫
P
I
grad
⁡
φ
·
K
P
ξ
P
I
B
∗
d
γ
-
∫
I
L
grad
⁡
φ
·
K
L
ξ
I
L
B
∗
d
γ
.By using the same process, the computation of the flux across[
P
I
] is(17)
-
∫
P
I
grad
⁡
φ
·
K
P
ξ
P
I
B
∗
d
γ
=
c
h
K
P
a
^
h
K
L
h
P
I
a
h
K
P
h
I
L
+
a
^
h
K
L
h
P
I
φ
P
-
φ
L
+
T
P
I
B
∗
+
h
P
I
d
h
K
P
a
h
K
P
h
I
L
+
a
^
h
K
L
h
P
I
+
c
h
K
P
h
I
L
b
^
h
K
L
-
b
h
K
P
h
A
∗
B
∗
a
h
K
P
h
I
L
+
a
^
h
K
L
h
P
I
φ
B
∗
-
φ
A
∗
,where (18)
c
h
K
P
=
1
cos
⁡
θ
h
P
,
I
ξ
A
∗
B
∗
P
t
K
P
ξ
P
I
B
∗
,
d
h
K
P
=
1
cos
⁡
θ
h
P
,
I
ξ
P
I
B
∗
t
K
P
ξ
P
I
B
∗(recall that θ
h
P
,
I denotes the angle defined by the vectors σ
P and ξ
A
∗
B
∗
P).Similarly, the computation of the flux across[
I
L
] leads to(19)
-
∫
I
L
grad
⁡
φ
·
K
L
ξ
I
L
B
∗
d
γ
=
c
^
h
K
L
a
h
K
P
h
I
L
a
h
K
P
h
I
L
+
a
^
h
K
L
h
P
I
φ
P
-
φ
L
+
T
I
L
B
∗
+
h
I
L
d
^
h
K
L
a
h
K
P
h
I
L
+
a
^
h
K
L
h
P
I
+
c
^
h
K
P
h
P
I
b
h
K
P
-
b
^
h
K
L
h
A
∗
B
∗
a
h
K
P
h
I
L
+
a
^
h
K
L
h
P
I
φ
B
∗
-
φ
A
∗
,where (20)
c
^
h
K
P
=
1
cos
⁡
θ
h
L
,
I
ξ
A
∗
B
∗
P
t
K
L
ξ
I
L
B
∗
,
d
^
h
K
L
=
1
cos
⁡
θ
h
L
,
I
ξ
I
L
B
∗
t
K
L
ξ
I
L
B
∗
.We should now formulate a global estimate for the truncation errors associated with the flux approximation over the pseudoedge[
P
I
L
].Proposition 10.
Under the same assumptions as those of Proposition7, there exists a positive number C mesh independent such that(21)
T
P
I
B
∗
+
T
I
L
B
∗
≤
C
h
2
.The above inequality shows the consistency of the flux approximation across the pseudoedge[
P
I
L
]. We summarize what precedes as(22)
Q
P
I
L
B
∗
≈
c
h
K
P
a
^
h
K
L
h
P
I
+
c
^
h
K
L
a
h
K
P
h
I
L
a
h
K
P
h
I
L
+
a
^
h
K
L
h
P
I
φ
P
-
φ
L
+
ω
h
P
,
L
,
I
h
A
∗
B
∗
a
h
K
P
h
I
L
+
a
^
h
K
L
h
P
I
φ
B
∗
-
φ
A
∗
,where we have set(23)
ω
h
P
,
L
,
I
=
a
h
K
P
h
I
L
+
a
^
h
K
L
h
P
I
d
h
K
P
h
P
I
+
d
^
h
K
L
h
I
L
+
h
P
I
h
I
L
b
^
h
K
L
-
b
h
K
P
c
h
K
P
-
c
^
h
K
L
.We deduce from (22) that an approximate flux balance equation over any dual cell C
B
∗ can be formulated as follows:(24)
∑
I
∈
E
P
∩
E
i
n
t
c
h
K
P
a
^
h
K
L
h
P
I
+
c
^
h
K
L
a
h
K
P
h
I
L
a
h
K
P
h
I
L
+
a
^
h
K
L
h
P
I
φ
P
-
φ
L
+
ω
h
P
,
L
,
I
h
A
∗
B
∗
a
h
K
P
h
I
L
+
a
^
h
K
L
h
P
I
φ
B
∗
-
φ
A
∗
+
∑
I
∈
E
B
∗
∩
E
e
x
t
h
I
E
a
K
E
d
h
K
E
-
c
h
K
E
b
h
K
E
a
h
K
E
h
B
∗
D
∗
φ
B
∗
-
φ
D
∗
≈
∫
C
B
∗
f
x
d
x
-
∫
Γ
B
∗
∩
Γ
g
d
γ
+
∑
I
∈
E
B
∗
∩
E
e
x
t
c
h
K
E
h
I
E
a
K
E
h
B
∗
D
∗
∫
B
∗
D
∗
g
d
γ
∀
B
∗
∈
D
,where P
∈
P, L
∈
P, and A
∗
∈
D are such that Γ
P
∩
Γ
L
=
A
∗
B
∗ and P
L
∩
A
∗
B
∗
=
I and where E
∈
P, D
∗
∈
D a boundary-vertex are such that B
∗
D
∗
⊂
Γ
∩
Γ
E and I
∈
B
∗
D
∗.Systems (14) and (24) naturally suggest defining a finite volume formulation of the diffusion problem (1)-(2) as follows.Findφ
-
P
P
∈
P and φ
-
D
∗
D
∗
∈
D such that(25)
∑
I
∈
E
P
∩
E
i
n
t
a
h
K
P
a
^
h
K
L
h
A
∗
B
∗
a
h
K
P
h
I
L
+
a
^
h
K
L
h
P
I
φ
-
P
-
φ
-
L
+
a
^
h
K
L
b
h
K
P
h
P
I
+
a
h
K
P
b
^
h
K
L
h
I
L
a
h
K
P
h
I
L
+
a
^
h
K
L
h
P
I
φ
-
B
∗
-
φ
-
A
∗
=
∫
C
P
f
x
d
x
-
∫
Γ
P
∩
Γ
g
d
γ
=
0
∀
P
∈
P
,
(26)
∑
I
∈
E
B
∗
∩
E
i
n
t
c
h
K
P
a
^
h
K
L
h
P
I
+
c
^
h
K
L
a
h
K
P
h
I
L
a
h
K
P
h
I
L
+
a
^
h
K
L
h
P
I
φ
-
P
-
φ
-
L
+
ω
h
P
,
L
,
I
h
A
∗
B
∗
a
h
K
P
h
I
L
+
a
^
h
K
L
h
P
I
φ
-
B
∗
-
φ
-
A
∗
+
∑
I
∈
E
B
∗
∩
E
e
x
t
h
I
E
a
K
E
d
h
K
E
-
c
h
K
E
b
h
K
E
a
h
K
E
h
B
∗
D
∗
φ
-
B
∗
-
φ
-
D
∗
=
∫
C
B
∗
f
x
d
x
-
∫
C
B
∗
∩
Γ
g
d
γ
+
∑
I
∈
E
B
∗
∩
E
e
x
t
c
h
K
E
h
I
E
a
K
E
h
B
∗
D
∗
∫
B
∗
D
∗
g
d
γ
∀
B
∗
∈
D
.Remark 11.
It is useful to note that(27)
∑
B
∗
∈
D
∑
I
∈
E
B
∗
∩
E
e
x
t
c
h
K
E
h
I
E
a
K
E
h
B
∗
D
∗
∫
B
∗
D
∗
g
d
γ
=
0
.


## 2.2. Existence of Solutions to System (25)-(26) Conditions for Uniqueness

(i) Matrix Properties of the DDFV Problem (25)-(26). It is easily seen that the symmetry of the matrix associated with the linear system (25)-(26) essentially follows from the symmetry of the diffusion coefficient K (see assumption (3)). We assume that all the cell-points and all the vertices (with respect to the primary mesh) are numbered. Therefore one can identify P and D with two disjoint subsets of the set of positive integers. To fix ideas, let us set(28)
P
≡
1
,
…
,
n
,
D
≡
n
+
1
,
…
,
n
+
m
.Then C
a
r
d
(
P
) denotes the total number of cell-points and C
a
r
d
(
D
) stands for the total number of vertices. On the other hand, define the subvectors Φ
c
c and Φ
v
c by(29)
Φ
c
c
=
φ
-
P
P
∈
P
,
Φ
v
c
=
φ
-
D
∗
D
∗
∈
Dand set (30)
M
=
A
B
B
t
C
,where(31)
M
Φ
c
c
Φ
v
c
=
F
c
c
F
v
c
,where F
c
c is a subvector with C
a
r
d
(
P
) components defined by the right-hand side of (25) and F
v
c is a subvector with C
a
r
d
(
D
) components defined by the right-hand side of (26).Remark 12.
The matrixM satisfies the following properties:(1)
1
≤
i
,
j
≤
m
M
i
j
=
M
j
i;
(2)
1
≤
i
,
j
≤
m
∑
j
=
1
m
M
i
j
=
0;
(3)
1
≤
i
,
j
≤
m
∑
i
=
1
m
M
i
j
=
0.Let us introduce two vectors ofR
n
+
m named F
c
c and F
v
c and defined by(32)
F
c
c
i
=
1
if
1
≤
i
≤
n
0
if
n
+
1
≤
i
≤
n
+
m
,
F
v
c
i
=
0
if
1
≤
i
≤
n
1
if
n
+
1
≤
i
≤
n
+
m
.Proposition 13 (characterization of Kernel space ofM).
(i)
The matrixM is singular.
(ii)
Moreover, letker
⁡
(
M
) denote the subset of R
n
+
m defined as follows:
⁡
ker
⁡
M
=
V
∈
R
n
+
m
,
M
V
=
0 (named Kernel space of M in the sequel); then we have(33)
ker
⁡
M
=
F
c
c
,
F
v
c
.Sketch for the Proof
(i)
The singularity ofM is an immediate consequence of Remark 12.
(ii)
Define the spaceR
P
,
D by(34)
R
P
,
D
=
V
=
V
c
c
,
V
v
c
;
V
c
c
=
V
P
P
∈
P
⊂
R
,
V
v
c
=
V
D
∗
D
∗
∈
D
⊂
Rand endow it with the seminorm ·
R
P
,
D defined as(35)
V
R
P
,
D
=
∑
P
,
L
∈
P
2
,
A
∗
,
B
∗
∈
D
2
s.t.
Γ
P
∩
Γ
L
=
A
∗
B
∗
V
P
-
V
L
2
+
V
B
∗
-
V
A
∗
2
+
∑
A
∗
,
B
∗
∈
D
2
s.t.
A
∗
B
∗
∈
E
e
x
t
V
B
∗
-
V
A
∗
2
1
/
2
∀
V
∈
R
P
,
D
.Then, find the above identification of the Kernel space of M using the following Lemma.Lemma 14.
The matrixM is positive; that is,(36)
V
c
c
,
V
v
c
M
V
c
c
V
v
c
≥
0
∀
V
∈
R
P
,
D
.Proof.
LetV
=
V
c
c
,
V
v
c
∈
R
P
,
D, where V
c
c
=
V
P
P
∈
P and V
v
c
=
V
A
∗
A
∗
∈
D.
It follows from what precedes that(37)
V
c
c
,
V
v
c
M
V
c
c
V
v
c
=
∑
P
∈
P
∑
I
∈
E
P
∩
E
i
n
t
a
h
K
P
a
^
h
K
L
h
A
∗
B
∗
a
h
K
P
h
I
L
+
a
^
h
K
L
h
P
I
V
P
-
V
L
V
P
+
a
^
h
K
L
b
h
K
P
h
P
I
+
a
h
K
P
b
^
h
K
L
h
I
L
a
h
K
P
h
I
L
+
a
^
h
K
L
h
P
I
V
B
∗
-
V
A
∗
V
P
+
∑
B
∗
∈
D
∑
I
∈
E
B
∗
∩
E
i
n
t
c
h
K
P
a
^
h
K
L
h
P
I
+
c
^
h
K
L
a
h
K
P
h
I
L
a
h
K
P
h
I
L
+
a
^
h
K
L
h
P
I
V
P
-
V
L
V
B
∗
+
ω
h
P
,
L
,
I
h
A
∗
B
∗
a
h
K
P
h
I
L
+
a
^
h
K
L
h
P
I
V
B
∗
-
V
A
∗
V
B
∗
+
∑
B
∗
∈
D
∑
I
∈
E
B
∗
∩
E
e
x
t
h
I
E
a
K
E
d
h
K
E
-
c
h
K
E
b
h
K
E
a
h
K
E
h
B
∗
D
∗
V
B
∗
-
V
D
∗
V
B
∗
.Define(38)
K
11
P
L
=
a
h
K
P
a
^
h
K
L
h
A
∗
B
∗
a
h
K
P
h
I
L
+
a
^
h
K
L
h
P
I
,
K
22
P
L
=
ω
h
P
,
L
,
I
h
A
∗
B
∗
a
h
K
P
h
I
L
+
a
^
h
K
L
h
P
I
,
Π
P
I
,
A
∗
,
B
∗
=
a
K
P
d
h
K
P
-
c
h
K
P
b
h
K
P
h
I
P
a
h
K
P
h
A
∗
B
∗
,where P
,
L
∈
P are two adjacent primary cells sharing the edge A
∗
B
∗ as interface containing the edge-point I. Then, relation (37) becomes(39)
V
c
c
,
V
v
c
M
V
c
c
V
v
c
=
∑
P
,
L
∈
P
2
,
A
∗
,
B
∗
∈
D
2
s.t.
Γ
P
∩
Γ
L
=
A
∗
B
∗
K
11
P
L
V
P
-
V
L
2
+
K
22
P
L
V
B
∗
-
V
A
∗
2
+
2
K
12
P
L
V
P
-
V
L
V
B
∗
-
V
A
∗
+
∑
P
∈
P
e
x
t
,
A
∗
,
B
∗
∈
D
2
∩
Γ
2
s.t.
A
∗
B
∗
∈
E
P
Π
P
I
,
A
∗
,
B
∗
V
B
∗
-
V
A
∗
2
,where P
e
x
t denotes the subset of P made of primary cells adjacent to the domain boundary Γ. Let us prove that the homogenized symmetric matrix K
P
L is positive definite; that is,(40)
K
11
P
L
K
22
P
L
-
K
12
P
L
2
>
0
.Setting(41)
Δ
P
L
=
K
11
P
L
K
22
P
L
-
K
12
P
L
2
,it is easy to check that(42)
Δ
P
L
=
N
1
a
h
K
P
d
h
K
P
-
b
h
K
P
2
+
N
2
a
^
h
K
L
d
^
h
K
L
-
b
^
h
K
L
2
,where we have set(43)
N
1
=
a
^
h
K
L
h
I
L
a
h
K
P
h
I
L
+
a
^
h
K
L
h
P
I
2
+
a
h
K
P
a
^
h
K
L
h
P
I
h
I
L
a
h
K
P
h
I
L
+
a
^
h
K
L
h
P
I
2
,
N
2
=
a
h
K
P
h
I
L
a
h
K
P
h
I
L
+
a
^
h
K
L
h
P
I
2
+
a
h
K
P
a
^
h
K
L
h
P
I
h
I
L
a
h
K
P
h
I
L
+
a
^
h
K
L
h
P
I
2which are strictly positive numbers.
Since the diffusion matrixK is symmetric and positive definite (see assumptions (3)-(4)), Cauchy-Schwarz’s inequality for the inner product associated with K ensures that(44)
a
h
K
P
d
h
K
P
-
b
h
K
P
2
>
0
,
a
^
h
K
L
d
^
h
K
L
-
b
^
h
K
L
2
>
0as either ξ
A
∗
B
∗
P and ξ
P
I
B
∗ or ξ
A
∗
B
∗
P and ξ
I
L
B
∗ are not collinear. Therefore, Δ
P
L
>
0 and thus K
P
L is a symmetric and positive definite matrix.
It follows from what precedes that the matrixK
P
L possesses strictly positive eigenvalues. Let λ
m
i
n
P
L be its least eigenvalue. So we have(45)
V
c
c
,
V
v
c
M
V
c
c
V
v
c
≥
∑
P
,
L
∈
P
2
,
A
∗
,
B
∗
∈
D
2
s.t.
Γ
P
∩
Γ
L
=
A
∗
B
∗
λ
min
⁡
P
L
V
P
-
V
L
2
+
V
B
∗
-
V
A
∗
2
+
∑
P
∈
P
e
x
t
,
A
∗
,
B
∗
∈
D
2
∩
Γ
2
s.t.
A
∗
B
∗
∈
E
P
Π
P
I
,
A
∗
,
B
∗
V
B
∗
-
V
A
∗
2
.Remarking that there exist two real numbers δ and η strictly positive depending exclusively on the geological structure of the domain such that(46)
δ
≤
K
11
P
K
22
P
-
K
12
P
2
≤
η
∀
P
∈
P
,the following result is easily seen.
Lemma  15.(47)
Π
P
I
,
A
∗
,
B
∗
=
K
11
P
K
22
P
-
K
12
P
2
h
P
I
a
h
K
P
h
A
∗
B
∗
.
Moreover
(48)
δ
ϖ
sin
⁡
θ
γ
m
a
x
<
Π
P
I
,
A
∗
,
B
∗
<
η
ϖ
γ
m
i
n
.
Thanks to inequality (45) and to Lemma 15 the proof of Lemma 14 ends.Lemma 16.
The matrixM satisfies the relation(49)
γ
V
R
P
,
D
2
≤
V
c
c
,
V
v
c
M
V
c
c
V
v
c
∀
V
∈
R
P
,
D
,where γ is a strictly positive number mesh independent and where ·
R
P
,
D is a seminorm defined on R
P
,
D by relation (35).Proof.
It is conducted with the same arguments as those developed in the previous proof, except that one should go much farther by proving that there existsρ
>
0 mesh independent such that(50)
ρ
≤
λ
m
i
n
P
L
∀
P
,
L
∈
P
.The eigenvalues λ of the symmetric positive definite matrix K
P
L satisfy the so-called characteristic equation associated with K
P
L; that is,(51)
λ
2
-
K
11
P
L
+
K
22
P
L
λ
+
K
11
P
L
K
22
P
L
-
K
12
P
L
2
=
0
.The least eigenvalue of K
P
L denoted by(52)
λ
min
⁡
P
L
=
K
11
P
L
+
K
22
P
L
-
Δ
2
,where Δ
=
K
11
P
L
+
K
22
P
L
2
-
4
K
11
P
L
K
22
P
L
-
K
12
P
L
2, is a strictly positive number. One can easily deduce that(53)
λ
min
⁡
P
L
≥
det
⁡
K
P
L
K
11
P
L
+
K
22
P
L
+
det
⁡
K
P
L
+
1
,where det
⁡
K
P
L
=
K
11
P
L
K
22
P
L
-
K
12
P
L
2 is a strictly positive number. We should bound the quantities K
11
P
L, K
22
P
L, and det
⁡
K
P
L by mesh independent strictly positive numbers. Let us start first with det
⁡
K
P
L
. We consider a change of coordinates by moving from the initial Cartesian coordinates to a local one, namely, J
,
C
∗
D
∗
→
/
C
∗
D
∗
→
,
ξ
C
∗
D
∗
⊥, where J is the edge-point on the interface C
∗
D
∗ between the cells C
G and C
H and where ξ
C
∗
D
∗
⊥ is a vector orthogonal to C
∗
D
∗ and oriented such that the basis change matrix M is a rotation. Denoting the permeability tensor of the cell C
G by K
G
=
K
i
j
G in the initial Cartesian coordinates and by K
~
G
=
K
~
i
j
G in the local coordinates we have(54)
K
~
G
=
M
-
1
K
G
M
.Similarly, we get for the cell C
H
(55)
K
~
H
=
M
-
1
K
H
M
.Then it is easy to check that(56)
a
h
K
G
d
h
K
G
-
b
h
K
G
2
=
K
~
11
G
K
~
22
G
-
K
~
12
G
2
;that is,(57)
a
h
K
G
d
h
K
G
-
b
h
K
G
2
=
det
⁡
K
~
G
=
det
⁡
K
G
,where det
⁡
· denotes the determinant. Similarly, we have for the cell C
H
(58)
a
^
h
K
H
d
^
h
K
H
-
b
^
h
K
H
2
=
det
⁡
K
~
H
=
det
⁡
K
H
.It follows from what precedes that(59)
det
⁡
K
P
L
=
N
1
det
⁡
K
G
+
N
2
det
⁡
K
H
,where N
1 and N
2 are given, respectively, by relations (43).
On one hand, we can deduce from (4) and Definition 6 that(60)
N
i
≥
1
2
ϖ
γ
min
⁡
sin
⁡
θ
γ
max
⁡
2
∀
i
=
1,2
.On the other hand, we remark that(61)
det
⁡
K
P
≥
min
⁡
det
⁡
K
s
,
s
∈
S
∀
P
∈
P
,where the set S (introduced in (8)) depends exclusively on the lithologic structure of the medium Ω. Then we deduce from (59), (60), and (61) that(62)
det
⁡
K
P
≥
1
2
ϖ
γ
min
⁡
sin
⁡
θ
γ
max
⁡
2
min
⁡
det
⁡
K
s
,
s
∈
S
∀
P
∈
P
.Remarking that(63)
det
⁡
K
P
L
≤
K
11
P
L
K
22
P
L
+
K
12
P
L
2and exploiting again (4) and Definition 6 lead to the following inequality:(64)
det
⁡
K
P
L
≤
2
γ
max
⁡
4
max
⁡
2
,
ϖ
ϖ
3
γ
min
⁡
2
sin
⁡
θ
4
.Thanks again to (4) and Definition 6 one can easily check that(65)
K
11
P
L
+
K
22
P
L
≤
γ
max
⁡
2
ϖ
sin
⁡
θ
2
1
+
2
ϖ
+
4
γ
max
⁡
ϖ
sin
⁡
θ
.Lemma 16 follows from the combination of (53), (62), (64), and (65).Proposition 17 (discrete compatibility condition).
The right-hand side of the discrete system (25)-(26) satisfies the following discrete compatibility condition:(66)
∑
P
∈
P
∫
C
P
f
x
d
x
-
∫
Γ
P
∩
Γ
g
d
γ
=
0
,
∑
B
∗
∈
D
∫
C
B
∗
f
x
d
x
-
∫
Γ
B
∗
∩
Γ
g
d
γ
+
∑
B
∗
∈
D
∑
I
∈
E
B
∗
∩
E
ext
c
h
K
E
h
I
E
a
K
E
h
B
∗
D
∗
∫
B
∗
D
∗
g
d
γ
=
0
.Proof.
First of all, note that the double summation in the previous proposition is equal to zero thanks to Remark11. Let us consider two vectors of R
n
+
m called F
c
c and F
v
c and defined by(67)
F
c
c
i
=
1
if
1
≤
i
≤
n
0
if
n
+
1
≤
i
≤
n
+
m
,
F
v
c
i
=
0
if
1
≤
i
≤
n
1
if
n
+
1
≤
i
≤
n
+
m
.According to the matrix form of the discrete system (see relation (31)), we have(68)
F
c
c
t
M
Φ
c
c
Φ
v
c
=
F
c
c
t
F
c
c
F
v
c
.It follows from Remark 12 that(69)
F
c
c
t
M
Φ
c
c
Φ
v
c
=
Φ
c
c
Φ
v
c
t
M
F
c
c
=
0
.Therefore(70)
0
=
F
c
c
t
F
c
c
F
v
c
=
∑
P
∈
P
∫
C
P
f
x
d
x
-
∫
Γ
P
∩
Γ
g
d
γ
.Similarly we have(71)
0
=
F
v
c
t
M
Φ
c
c
Φ
v
c
=
F
v
c
t
F
c
c
F
v
c
=
∑
B
∗
∈
D
∫
C
B
∗
f
x
d
x
-
∫
C
B
∗
∩
Γ
g
d
γ
.This ends the proof of Proposition 17.An immediate consequence of the previous result is the following existence result for a solution to the discrete problem (25)-(26).Proposition 18 (existence result).
(
1) The linear system (25)-(26) possesses an infinite number of solutions.
(
2) More precisely, if s
h is a solution to the discrete system (25)-(26), then s
h
+
ker
⁡
M is the set of all solutions for this system.The question to know how to get the physical solution of the model problem has a natural answer; that is, one should follow the same way as the continuous problem analysis. Since the dimension ofker
⁡
M is 2, required are two (linearly independent) discrete versions of the null average condition (7). This naturally leads to the following result.Proposition 19 (uniqueness result).
Under the following (null average) conditions,(72)
i
∑
P
∈
P
m
e
s
C
P
φ
-
P
=
0
,
i
i
∑
D
∗
∈
D
m
e
s
C
D
∗
φ
-
D
∗
=
0
,the discrete problem (25)-(26) possesses a unique solution.Proof.
Let us consider the space(73)
R
P
,
D
=
V
=
V
c
c
,
V
v
c
;
V
c
c
=
V
P
P
∈
P
⊂
R
,
V
v
c
=
V
D
∗
D
∗
∈
D
⊂
R
.One knows from Lemma 14 that(74)
V
t
M
V
≥
0
∀
V
∈
R
P
,
D
.Then it follows that, for all V in the subspace of R
P
,
D made up of V
=
V
P
P
∈
P
,
V
D
∗
D
∗
∈
D such that conditions (72) are fulfilled, one has V
t
M
V
=
0
⇔
V
=
0.


## 3. Stability and Error Estimates

We deal here with a theoretical analysis of the solution for the discrete system (25), (26), and (72). Recall that the existence and uniqueness of that solution (under explicit conditions) is proven in the previous section. We assume in what follows that the primary mesh is regular in the sense of Definition 6.


### 3.1. Preliminaries and a Stability Result

Let us consider the auxiliary meshA introduced in the previous section (see Figure 4). Note that each mesh cell of A contains either one cell-point or one corner point and only one which should be lying inside or on the boundary of Ω. In the sequel, anode is a cell-point or a corner point with respect to the primary mesh.Definition 20.
An auxiliary mesh cell is degenerate if the corresponding node is lying on the boundary ofΩ
.In the sequel, we will say sometimes “auxiliary cell” instead of “auxiliary mesh cell.” We denote byE
(
A) the subspace of R
P
,
D made of functions v which satisfies conditions (72). This space is obviously nonempty as there is the null function. Moreover, the solution of the discrete system (25), (26), and (72) could clearly be identified with one (and only one) function from E
(
A). In the sequel we denote by φ
-
h such a function named “cellwise constant (approximate) solution” of the diffusion problem (1)-(2). Let us endow R
P
,
D with the following discrete seminorm.For allv
≡
V
P
P
∈
P
,
V
D
∗
D
∗
∈
D
∈
R
P
,
D, define(75)
v
R
P
,
D
=
∑
P
,
L
∈
P
2
,
A
∗
,
B
∗
∈
D
2
s.t.
Γ
P
∩
Γ
L
=
A
∗
B
∗
V
P
-
V
L
2
+
V
B
∗
-
V
A
∗
2
+
∑
A
∗
,
B
∗
∈
D
2
s.t.
A
∗
B
∗
∈
E
e
x
t
V
B
∗
-
V
A
∗
2
1
/
2
∀
V
∈
R
P
,
D
,where “s.t.” means “such that.” As this seminorm is actually a discrete energy norm for the space E
(
A), we denote it by ·
A as far as functions from E
(
A) are concerned. A norm on the space R
P
,
D is defined by the mapping (76)
·
R
P
,
D
=
·
R
P
,
D
2
+
·
L
2
Ω
2
1
/
2
.We focus here on the proof of the stability of the solution to the system of (25), (26), and (72). For this purpose, we need to introduce as in [31] a preliminary result, namely, a discrete version of Poincaré inequality which is based upon the following ingredients. Consider the following spaces(77)
R
P
=
v
=
v
P
P
∈
P
;
v
P
∈
R
∀
P
∈
P
,
R
D
=
v
=
v
D
∗
D
∗
∈
D
;
v
D
∗
∈
R
∀
D
∗
∈
Dand the linear operators Π
P and Π
D defined as follows: (78)
v
∈
E
A
⟼
Π
P
v
∈
E
P
,
v
∈
E
A
⟼
Π
D
v
∈
E
Dwith Π
P
v
=
v
P
P
∈
P and Π
D
v
=
v
D
∗
D
∗
∈
D, where E
P and E
D are, respectively, subspaces of R
P and R
D made of functions satisfying conditions (72)-(i) and (72)-(ii), respectively.Let us equip the function spacesE
P and E
D with the following discrete energy norms:(79)
v
P
=
∑
P
,
L
∈
P
2
,
A
∗
,
B
∗
∈
D
2
s.t.
Γ
P
∩
Γ
L
=
A
∗
B
∗
v
P
-
v
L
2
1
/
2
,
v
D
=
∑
P
,
L
∈
P
2
,
A
∗
,
B
∗
∈
D
2
s.t.
Γ
P
∩
Γ
L
=
A
∗
B
∗
v
B
∗
-
v
A
∗
2
+
∑
A
∗
,
B
∗
∈
D
2
s.t.
A
∗
B
∗
∈
E
e
x
t
V
B
∗
-
V
A
∗
2
1
/
2
.These mappings are only seminorms, respectively, for R
P and R
D that one can transform into norms for these spaces as follows:(80)
·
R
P
=
·
R
P
2
+
·
L
2
Ω
2
1
/
2
,
·
R
D
=
·
R
D
2
+
·
L
2
Ω
2
1
/
2
,where we have set ·
R
P
=
·
P and ·
R
D
=
·
D. The following results are key ingredients for proving the stability of the discrete solution in the sense of the discrete energy norm defined previously on the space E
A.Lemma 21 (discrete versions of Poincaré inequality).
We have the following inequalities:(81)
·
L
2
Ω
2
≤
C
v
P
2
+
2
m
e
s
Ω
∫
Ω
v
x
d
x
2
∀
v
∈
R
P
,
(82)
·
L
2
Ω
2
≤
C
v
D
2
+
2
m
e
s
Ω
∫
Ω
v
x
d
x
2
∀
v
∈
R
D
,where C represents diverse strictly positive numbers mesh independent.Note that the right-hand side of inequalities (81) and (82) is, respectively, norms for R
P and R
D. Since Ω is bounded, the previous lemma permits seeing that these norms are equivalent to standard ones, namely, (80). For proving the preceding lemma one can use the same arguments as in [36].Lemma 22 (a key-result).
LetP
-
,
E
- be a regular mesh system (defined on) in the sense of Definition 6 and denote by D the corresponding dual mesh. For v
∈
E
(
P
-
), denote by v
P the value of v in the control volume P. Let γ
-
v be a discrete trace function defined a.e. (for the (d
-
1)-Lebesgue measure) by γ
-
v
=
v
P on Γ
P
∩
Γ, for all P
∈
P
-
e
x
t, where P
-
e
x
t denotes the set of mesh elements adjacent to the domain boundary Γ. Then(83)
γ
-
v
L
2
Γ
≤
C
v
P
-
+
·
L
2
Ω
∀
v
∈
E
P
-
.Similarly, we have(84)
γ
-
v
L
2
Γ
≤
C
v
D
-
+
·
L
2
Ω
∀
v
∈
E
D
-
,where C stands for diverse positive numbers mesh independent.Proof.
See, for instance, [36].Remark 23.
It is more than useful to note that(85)
v
A
2
=
v
P
P
2
+
v
D
D
2
∀
v
∈
E
A
.Let us give now one of the main results of this section.Proposition 24 (stability result).
The cellwise constant approximate solutionφ
-
h of the diffusion problem (1)-(2) satisfies the following inequality:(86)
φ
-
h
A
≤
C
f
L
2
Ω
+
g
L
2
Γ
,where C is a strictly positive real number not depending on the spatial discretization.Proof.
Multiplying (25) by φ
-
P and (26) by φ
-
D
∗ and summing over P
∈
P and D
∗
∈
D, respectively, lead to(87)
Φ
c
c
Φ
v
c
A
B
B
t
C
Φ
c
c
Φ
v
c
=
Φ
c
c
Φ
v
c
F
c
c
F
v
c
.Recall that(88)
F
c
c
=
∫
C
P
f
x
d
x
-
∫
Γ
P
∩
Γ
g
d
γ
P
∈
P
,
F
v
c
=
∫
C
B
∗
f
x
d
x
-
∫
C
B
∗
∩
Γ
g
d
γ
+
∑
A
∗
∈
N
B
∗
,
I
∈
E
B
∗
∩
E
e
x
t
with
Γ
C
A
∗
∩
Γ
C
B
∗
=
P
I
c
h
K
E
h
I
E
a
K
E
h
B
∗
D
∗
∫
B
∗
D
∗
g
d
γ
B
∗
∈
D
.Let us set (89)
LHS
=
Φ
c
c
Φ
v
c
A
B
B
t
C
Φ
c
c
Φ
v
c
,
RHS
=
Φ
c
c
Φ
v
c
F
c
c
F
v
c
.We know from Lemma 16 that there exists a strictly positive number ρ mesh independent such that(90)
ρ
φ
-
h
A
2
≤
LHS
.On the other hand, according to Remark 11 we have(91)
RHS
=
∑
P
∈
P
∫
C
P
f
φ
-
P
d
x
-
∑
P
∈
P
∫
Γ
P
∩
Γ
g
φ
-
P
d
γ
+
∑
B
∗
∈
D
∫
C
B
∗
f
x
d
x
-
∑
B
∗
∈
D
∫
Γ
B
∗
∩
Γ
g
φ
-
B
∗
d
γ
+
∑
B
∗
∈
D
∑
A
∗
∈
N
B
∗
s.t.
,
∃
I
∈
E
B
∗
∩
E
A
∗
∩
E
e
x
t
,
∃
L
∈
P
,
with
Γ
A
∗
∩
Γ
B
∗
=
L
I
φ
-
B
∗
c
h
K
L
h
L
I
a
K
L
h
A
∗
B
∗
∫
A
∗
B
∗
g
d
γ
.Cauchy-Schwarz’s inequality leads to (92)
RHS
≤
f
L
2
Ω
φ
-
h
P
L
2
Ω
+
g
L
2
Γ
γ
-
P
φ
-
h
P
L
2
Γ
+
f
L
2
Ω
φ
-
h
D
L
2
Ω
+
g
L
2
Γ
γ
-
D
φ
-
h
D
L
2
Γ
+
∑
A
∗
,
B
∗
∈
D
e
x
t
s.t.
,
∃
I
∈
E
B
∗
∩
E
A
∗
,
∃
L
∈
P
,
with
Γ
A
∗
∩
Γ
B
∗
=
L
I
φ
-
B
∗
-
φ
-
A
∗
c
h
K
L
h
L
I
a
K
L
h
A
∗
B
∗
∫
A
∗
B
∗
g
d
γ
,where γ
-
P and γ
-
D are, respectively, discrete trace operators associated with grids P and D. Thanks to relation (4), Definition 6, and Cauchy-Schwarz inequality we get(93)
RHS
≤
f
L
2
Ω
φ
-
h
P
L
2
Ω
+
φ
-
h
D
L
2
Ω
+
g
L
2
Γ
γ
-
P
φ
-
h
P
L
2
Γ
+
γ
-
D
φ
-
h
D
L
2
Γ
+
1
ϖ
γ
m
i
n
sin
⁡
θ
∑
A
∗
,
B
∗
∈
D
e
x
t
s.t.
,
E
A
∗
∩
E
B
∗
≠
∅
φ
-
B
∗
-
φ
-
A
∗
∫
A
∗
B
∗
g
d
γ
.Accounting with the fact that φ
-
h satisfies null average condition on the grids P and D, Lemma 21 ensures that(94)
φ
-
h
P
L
2
Ω
+
φ
-
h
D
L
2
Ω
≤
C
φ
-
h
P
P
+
φ
-
h
D
D
,whereas Lemma 22 guarantees that(95)
γ
-
P
φ
-
h
P
L
2
Γ
+
γ
-
D
φ
-
h
D
L
2
Γ
≤
C
φ
-
h
P
P
+
φ
-
h
D
D
.A double application of Cauchy-Schwarz’s inequality leads to(96)
∑
A
∗
,
B
∗
∈
D
e
x
t
s.t.
,
E
A
∗
∩
E
B
∗
≠
∅
φ
-
B
∗
-
φ
-
A
∗
∫
A
∗
B
∗
g
d
γ
≤
meas
⁡
Γ
1
/
2
φ
-
h
D
D
g
L
2
Γ
.Combining relations (93), (94), (95), and (96) on the one hand and using the inequality of Cauchy-Schwarz and Remark 23 on the other hand yield(97)
RHS
≤
C
f
L
2
Ω
+
g
L
2
Γ
φ
-
h
A
.Recall that C represents diverse strictly positive numbers mesh independent. The stability result follows from relations (90) and (97). The proof of Proposition 24 ends.


### 3.2. Error Estimates

Following the original technique exposed in [28], we investigate in this subsection the error estimates for the finite volume approximate solution to the model problem. Recall that φ denotes the exact solution to the model problem while φ
-
h
∈
E
A is the cellwise constant approximate solution obtained from a DDFV formulation of this model problem. View also φ
-
h as a representation in terms of function of the vector Φ
=
Φ
c
c
Φ
v
c
t, where Φ
c
c
=
φ
-
P
P
∈
P and Φ
v
c
=
φ
-
D
∗
D
∗
∈
D. Recall that the vector Φ is defined above as the unique solution of the DDFV discrete system, satisfying null average conditions on grids P and D (see (25)-(26) and (72)). Let us set(98)
E
P
=
φ
x
P
-
φ
-
h
x
P
≡
φ
P
-
φ
-
P
∀
P
∈
P
,
E
D
∗
=
φ
x
D
∗
-
φ
-
h
x
D
∗
≡
φ
D
∗
-
φ
-
D
∗
∀
D
∗
∈
D
.Note that the components of the so-called error-vector E
P
P
∈
P
,
E
D
∗
D
∗
∈
D
t
∈
R
P
,
D can be viewed as values in auxiliary cells of the so-called error-function E
h defined a.e. in Ω. More precisely let us recall that any auxiliary cell A is associated with a unique cell-point or a unique vertex. Let E
A be the error corresponding to the auxiliary cell A; that is,(99)
E
A
=
φ
x
P
-
φ
-
h
x
P
if
P
is  a  cell  point  associated  with
A
φ
x
D
∗
-
φ
-
h
x
D
∗
if
D
∗
is  a  cell  point  associated  with
A
.Since the auxiliary mesh define a partition of Ω, one can set(100)
E
h
=
∑
A
∈
A
E
A
χ
A
x
,where χ
A is the characteristic function of any cell A.Our first purpose is to show that the error-vector is a solution to a square linear system of the same type as the discrete system (25)-(26) in the sense that both of them are associated with the same matrix M defined above. Note that unfortunately E
h
∈
R
P
,
D is not in E
(
A
) as it does meet the null average conditions (72) on grids P and D. So the error estimates, which are our main purpose in this section, should be investigated in the sense of the norm defined by (76) on the space R
P
,
D. In this connection we first should prove that(101)
E
h
R
P
,
D
≤
C
hand in the next step we should prove that(102)
E
h
L
2
Ω
≤
C
h
.(i) Let us start with proving that the components of the error-vector are solution to a square linear system of the same type as the discrete system (25)-(26).Accounting with truncation errors, one can see from the DDFV flux computations that the exact flux across the edgeA
∗
B
∗ viewed as part of the boundary of the primary cell C
P is given by(103)
Q
A
∗
B
∗
P
=
a
h
K
P
a
^
h
K
L
h
A
∗
B
∗
a
h
K
P
h
I
L
+
a
^
h
K
L
h
P
I
φ
P
-
φ
L
+
a
^
h
K
L
b
h
K
P
h
P
I
+
a
h
K
P
b
^
h
K
L
h
I
L
a
h
K
P
h
I
L
+
a
^
h
K
L
h
P
I
φ
B
∗
-
φ
A
∗
+
R
I
A
∗
,
B
∗
,
L
P
,where the truncation error R
I
A
∗
,
B
∗
,
L
P is given by(104)
R
I
A
∗
,
B
∗
,
L
P
=
1
+
a
h
K
P
h
I
L
a
h
K
P
h
I
L
+
a
^
h
K
L
h
P
I
T
I
A
∗
,
B
∗
P
+
a
h
K
P
h
I
L
a
h
K
P
h
I
L
+
a
^
h
K
L
h
P
I
T
I
A
∗
,
B
∗
Land meets the following inequality (see Proposition 7):(105)
R
I
A
∗
,
B
∗
,
L
P
≤
C
h
2
.On the other hand, the exact fluxes across the pseudoedge [
P
I
L
] for I
∈
Ω and across the edge [
P
I
] for I
∈
Γ, both of them acting as part of the boundary of the dual cell C
B
∗, are given by the following relations (accounting with truncation errors):(106)
Q
P
I
L
B
∗
=
c
h
K
P
a
^
h
K
L
h
P
I
+
c
^
h
K
L
a
h
K
P
h
I
L
a
h
K
P
h
I
L
+
a
^
h
K
L
h
P
I
φ
P
-
φ
L
+
ω
h
P
,
L
,
I
h
A
∗
B
∗
a
h
K
P
h
I
L
+
a
^
h
K
L
h
P
I
φ
B
∗
-
φ
A
∗
+
R
I
P
,
L
,
A
∗
B
∗
,
Q
P
I
B
∗
=
a
h
K
P
d
h
K
P
-
c
h
K
P
2
h
P
I
a
h
K
P
h
A
∗
B
∗
φ
B
∗
-
φ
A
∗
-
c
h
K
P
h
P
I
a
K
P
h
A
∗
B
∗
∫
A
∗
B
∗
g
d
γ
+
R
I
P
,
A
∗
B
∗
,where the truncation errors R
I
P
,
L
,
A
∗
B
∗ and R
I
(
P
,
A
∗
)
B
∗ meet the following inequalities (according to Propositions 7 and 10):(107)
R
I
P
,
L
,
A
∗
B
∗
≤
C
h
2
,
R
I
P
,
A
∗
B
∗
≤
C
h
2
.The system of equations satisfied by theexact nodal potentials φ
P
P
∈
P and φ
D
∗
D
∗
∈
D reads as(108)
∑
L
∈
N
P
s.t.
∃
I
∈
E
,
∃
A
∗
,
B
∗
∈
D
,
with
I
∈
E
P
∩
E
L
,
Γ
P
∩
Γ
L
=
A
∗
B
∗
a
h
K
P
a
^
h
K
L
h
A
∗
B
∗
a
h
K
P
h
I
L
+
a
^
h
K
L
h
P
I
φ
P
-
φ
L
+
a
^
h
K
L
b
h
K
P
h
P
I
+
a
h
K
P
b
^
h
K
L
h
I
L
a
h
K
P
h
I
L
+
a
^
h
K
L
h
P
I
φ
B
∗
-
φ
A
∗
=
∫
C
P
f
x
d
x
-
∫
Γ
P
∩
Γ
g
d
γ
-
∑
L
∈
N
P
s.t.
∃
I
∈
E
,
∃
A
∗
,
B
∗
∈
D
,
with
I
∈
E
P
∩
E
L
,
Γ
P
∩
Γ
L
=
A
∗
B
∗
R
I
A
∗
,
B
∗
,
L
P
∀
P
∈
P
,
∑
A
∗
∈
N
B
∗
s.t.
∃
I
∈
E
,
∃
P
,
L
∈
P
,
with
I
∈
E
A
∗
∩
E
B
∗
,
Γ
A
∗
∩
Γ
B
∗
=
P
I
L
c
h
K
P
a
^
h
K
L
h
P
I
+
c
^
h
K
L
a
h
K
P
h
I
L
a
h
K
P
h
I
L
+
a
^
h
K
L
h
P
I
φ
P
-
φ
L
+
ω
h
P
,
L
,
I
h
A
∗
B
∗
a
h
K
P
h
I
L
+
a
^
h
K
L
h
P
I
φ
B
∗
-
φ
A
∗
+
∑
A
∗
∈
N
B
∗
s.t.
∃
I
∈
E
,
∃
P
∈
P
with
I
∈
E
A
∗
∩
E
B
∗
∩
E
e
x
t
,
Γ
A
∗
∩
Γ
B
∗
=
P
I
a
K
P
d
h
K
P
-
c
h
K
P
2
h
P
I
a
h
K
P
h
A
∗
B
∗
φ
B
∗
-
φ
D
∗
=
∫
C
B
∗
f
x
d
x
-
∫
C
B
∗
∩
Γ
g
d
γ
-
∑
A
∗
∈
N
B
∗
s.t.
∃
I
∈
E
,
∃
P
,
L
∈
P
with
I
∈
E
A
∗
∩
E
B
∗
,
Γ
A
∗
∩
Γ
B
∗
=
P
I
L
R
I
P
,
L
,
A
∗
B
∗
-
∑
A
∗
∈
N
B
∗
s.t.
∃
I
∈
E
,
∃
P
∈
P
with
I
∈
E
A
∗
∩
E
B
∗
∩
E
e
x
t
,
Γ
A
∗
∩
Γ
B
∗
=
P
I
R
I
P
,
A
∗
B
∗
+
∑
A
∗
∈
N
B
∗
s.t.
∃
I
∈
E
,
∃
P
∈
P
with
I
∈
E
A
∗
∩
E
B
∗
∩
E
e
x
t
,
Γ
A
∗
∩
Γ
B
∗
=
P
I
c
h
K
P
h
P
I
a
K
P
h
A
∗
B
∗
∫
A
∗
B
∗
g
d
γ
∀
B
∗
∈
D
,where ω
h
P
,
L
,
I is defined in (23). Recall that N
· is the set of neighboring cells of a given (primary or dual) cell.Due to the conservation of our flux approximation schemes, the truncation errors naturally obey to the following relation:(109)
R
I
·
,
·
,
b
a
+
R
I
·
,
·
,
a
b
=
0
∀
I
∈
E
i
n
t
,
I
=
E
a
∩
E
b
,
(110)
R
I
·
,
b
a
+
R
I
·
,
a
b
=
0
∀
I
∈
E
e
x
t
,
I
=
E
a
∩
E
b
,where a and b are two cell-points or corner points associated with primary or dual adjacent cells. Note that, in the previous equality, the truncation errors are written in a very formal way. This is with a view to involving both primary and dual adjacent cells. Let us introduce now the set of diamond cells called in what follows the diamond mesh and denoted as M (the concept of diamond cell has been introduced for a different usage in earlier works of some authors: see, e.g., [22, 25, 34]). Each diamond cell is associated with one edge-point and vice versa (see Figure 6). A diamond cell is declared degenerate if the corresponding edge-point is lying on the boundary of Ω. An example of a degenerate diamond cell is provided in Figure 6.Figure 6
Examples of diamond cells. (a) A normal diamond cell. (b) A degenerate diamond cell.
(a)
(b)The following assumption plays a key-role in our proof of the convergence of the piecewise constant (approximate) solutionφ
-
h.(111)
∃
ν
∈
R
+
∗
such  that
ν
h
2
≤
meas
⁡
M
∀
P
∈
M
,where m
e
a
s
(
·
) is the Lebesgue measure in any spatial dimension.Recall that we have setE
P
=
φ
P
-
φ
-
P, for P
∈
P, and E
D
∗
=
φ
D
∗
-
φ
-
D
∗, for D
∗
∈
D. An adequate linear combination of equations from system (108)-(109) with those from system (25)-(26) shows that the quantities E
P
=
φ
P
-
φ
-
P
P
∈
P and E
D
∗
=
φ
D
∗
-
φ
-
D
∗
D
∗
∈
D satisfy the following equations: (112)
∑
L
∈
N
P
s.t.
∃
I
∈
E
,
∃
A
∗
,
B
∗
∈
D
,
with
I
∈
E
P
∩
E
L
,
Γ
P
∩
Γ
L
=
A
∗
B
∗
a
h
K
P
a
^
h
K
L
h
A
∗
B
∗
a
h
K
P
h
I
L
+
a
^
h
K
L
h
P
I
E
P
-
E
L
+
a
^
h
K
L
b
h
K
P
h
P
I
+
a
h
K
P
b
^
h
K
L
h
I
L
a
h
K
P
h
I
L
+
a
^
h
K
L
h
P
I
E
B
∗
-
E
A
∗
=
∑
L
∈
N
P
s.t.
∃
I
∈
E
,
∃
A
∗
,
B
∗
∈
D
,
with
I
∈
E
P
∩
E
L
,
Γ
P
∩
Γ
L
=
A
∗
B
∗
R
I
A
∗
,
B
∗
,
L
P
∀
P
∈
P
,
(113)
∑
A
∗
∈
N
B
∗
s.t.
∃
I
∈
E
,
∃
P
,
L
∈
P
,
with
I
∈
E
A
∗
∩
E
B
∗
,
Γ
A
∗
∩
Γ
B
∗
=
P
I
L
c
h
K
P
a
^
h
K
L
h
P
I
+
c
^
h
K
L
a
h
K
P
h
I
L
a
h
K
P
h
I
L
+
a
^
h
K
L
h
P
I
E
P
-
E
L
+
ω
h
P
,
L
,
I
h
A
∗
B
∗
a
h
K
P
h
I
L
+
a
^
h
K
L
h
P
I
E
B
∗
-
E
A
∗
+
∑
A
∗
∈
N
B
∗
s.t.
∃
I
∈
E
,
∃
P
∈
P
with
I
∈
E
A
∗
∩
E
B
∗
∩
E
e
x
t
,
Γ
A
∗
∩
Γ
B
∗
=
P
I
a
h
K
P
d
h
K
P
-
c
h
K
P
2
h
P
I
a
h
K
P
h
A
∗
B
∗
E
B
∗
-
E
A
∗
=
-
∑
A
∗
∈
N
B
∗
s.t.
∃
I
∈
E
,
∃
P
,
L
∈
P
with
I
∈
E
A
∗
∩
E
B
∗
,
Γ
A
∗
∩
Γ
B
∗
=
P
I
L
R
I
P
,
L
,
A
∗
B
∗
-
∑
A
∗
∈
N
B
∗
s.t.
∃
I
∈
E
,
∃
P
∈
P
with
I
∈
E
A
∗
∩
E
B
∗
∩
E
e
x
t
,
Γ
A
∗
∩
Γ
B
∗
=
P
I
R
I
P
,
A
∗
B
∗
∀
B
∗
∈
D
.(ii) Let us now prove thatE
h
R
P
,
D
≤
C
h, where C
>
0 is a mesh independent positive number.Multiplying (112) by E
P and (113) by E
B
∗ and summing over P
∈
P and B
∗
∈
D yield (thanks to Lemma 16)(114)
γ
E
h
R
P
,
D
2
≤
∑
P
∈
P
-
E
p
∑
L
∈
N
P
s.t.
∃
I
∈
E
P
∩
E
L
,
∃
A
∗
,
B
∗
∈
D
,
with
Γ
P
∩
Γ
L
=
A
∗
B
∗
R
I
A
∗
,
B
∗
,
L
P
+
∑
B
∗
∈
D
-
E
B
∗
∑
A
∗
∈
N
B
∗
s.t.
∃
I
∈
E
A
∗
∩
E
B
∗
,
∃
P
,
L
∈
P
with
Γ
A
∗
∩
Γ
B
∗
=
P
I
L
R
I
P
,
L
,
A
∗
B
∗
-
E
B
∗
∑
A
∗
∈
N
B
∗
s.t.
∃
I
∈
E
A
∗
∩
E
B
∗
∩
E
e
x
t
,
∃
P
∈
P
with
Γ
A
∗
∩
Γ
B
∗
=
P
I
R
I
P
,
A
∗
B
∗
,where γ is some strictly positive real number mesh independent.Thanks to (109)-(110) (consequence of the conservation of the flux approximation scheme) we have(115)
γ
E
h
R
P
,
D
2
≤
∑
I
∈
E
s.t.
∃
P
,
L
∈
P
with
I
∈
E
P
∩
E
L
,
∃
A
∗
,
B
∗
∈
D
with
Γ
P
∩
Γ
L
=
A
∗
B
∗
E
P
-
E
L
R
I
A
∗
,
B
∗
,
L
P
+
∑
I
∈
E
s.t.
∃
A
∗
,
B
∗
∈
D
with
I
∈
E
A
∗
∩
E
B
∗
,
∃
P
,
L
∈
P
with
Γ
A
∗
∩
Γ
B
∗
=
P
I
L
E
B
∗
-
E
A
∗
R
I
P
,
L
,
A
∗
B
∗
+
∑
I
∈
E
s.t.
∃
A
∗
,
B
∗
∈
D
with
I
∈
E
A
∗
∩
E
B
∗
∩
Γ
,
∃
P
∈
P
with
Γ
A
∗
∩
Γ
B
∗
=
P
I
E
B
∗
-
E
A
∗
R
I
P
,
A
∗
B
∗
.Therefore(116)
γ
E
h
R
P
,
D
2
≤
∑
I
P
,
L
,
A
∗
,
B
∗
∈
E
i
n
t
R
I
P
,
L
,
A
∗
,
B
∗
max
⁡
E
P
-
E
L
+
E
B
∗
-
E
A
∗
+
∑
I
P
,
A
∗
,
B
∗
∈
E
e
x
t
R
I
P
,
A
∗
,
B
∗
max
⁡
E
B
∗
-
E
A
∗
,where we have set(117)
R
I
P
,
L
,
A
∗
,
B
∗
max
⁡
=
max
⁡
R
I
A
∗
,
B
∗
,
L
P
,
R
I
P
,
L
,
A
∗
B
∗
,
R
I
P
,
A
∗
,
B
∗
=
R
I
P
,
A
∗
B
∗
.It is important to mention that, due to (105) and (107), we have(118)
0
≤
R
I
P
,
L
,
A
∗
,
B
∗
max
⁡
≤
C
h
2
,
0
≤
R
I
P
,
A
∗
,
B
∗
≤
C
h
2
.Define S
I
P
,
L
,
A
∗
,
B
∗
=
2
D Lebesgue measure of the diamond cell defined by the points P
,
L
,
A
∗
,
B
∗ and associated with I
∈
E
i
n
t and S
I
P
,
A
∗
,
B
∗
=
2
D Lebesgue measure of the diamond cell defined by the points P
,
A
∗
,
B
∗ and associated with I
∈
E
e
x
t.Therefore, thanks to Cauchy-Schwarz’s inequality, it follows from (116) that(119)
γ
E
h
R
P
,
D
2
≤
∑
I
P
,
L
,
A
∗
,
B
∗
∈
E
i
n
t
S
I
P
,
L
,
A
∗
,
B
∗
+
∑
I
P
,
A
∗
,
B
∗
∈
E
i
n
t
S
I
P
,
A
∗
,
B
∗
1
/
2
∑
I
P
,
L
,
A
∗
,
B
∗
∈
E
i
n
t
2
R
I
P
,
L
,
A
∗
,
B
∗
max
⁡
2
S
I
P
,
L
,
A
∗
,
B
∗
E
P
-
E
L
2
+
E
B
∗
-
E
A
∗
2
1
/
2
+
∑
I
P
,
A
∗
,
B
∗
∈
E
e
x
t
2
R
I
P
,
A
∗
,
B
∗
2
S
I
P
,
A
∗
,
B
∗
E
B
∗
-
E
A
∗
2
1
/
2
.One concludes with the help of assumption (111) and relation (118) that(120)
E
h
R
P
,
D
≤
C
h
,where C is a strictly positive constant that is mesh independent.(iii) Let us now prove thatE
h
L
2
Ω
≤
C
h.We start with setting(121)
A
P
=
E
∈
A
;
∃
P
∈
P
with
x
P
∈
E
,
A
D
=
E
∈
A
;
∃
D
∗
∈
D
with
x
D
∗
∈
E
.Recall A is the auxiliary grid and x
P is a cell-point while x
D
∗ is a vertex (with respect to the primary grid introduced in a preceding section).Define(122)
Ω
P
=
⋃
E
∈
A
P
E
,
Ω
D
=
⋃
E
∈
A
D
E
.Note that Ω
P and Ω
D are a partition of Ω in the sense that(123)
i
Ω
P
∩
Ω
D
=
∅
,
i
i
Ω
-
=
Ω
-
P
∪
Ω
-
D
.Two main steps are necessary in our technique to get the estimates of φ
-
φ
-
L
2
Ω; before exposing these steps, we develop some preliminaries.(124)
∫
Ω
φ
x
-
φ
-
h
x
2
d
x
=
∫
Ω
P
φ
x
-
φ
-
h
x
2
d
x
+
∫
Ω
D
φ
x
-
φ
-
h
x
2
d
x
≤
∫
Ω
φ
x
-
Π
P
φ
-
h
x
2
d
x
+
∫
Ω
φ
x
-
Π
D
φ
-
h
x
2
d
x
,where Π
P and Π
D are operators from the function space E
(
A
) into the function spaces A
P and A
D, respectively (see (78) for the definition of Π
P and Π
D).Letλ
P and λ
D be two real numbers and φ
aux
P and φ
aux
D two real functions defined almost everywhere in Ω such that(125)
i
φ
aux
P
=
φ
x
-
λ
P
,
i
i
∑
P
∈
P
meas
⁡
C
P
φ
aux
P
x
P
=
0
,
(126)
i
φ
aux
D
=
φ
x
-
λ
D
,
i
i
∑
D
∗
∈
D
meas
⁡
C
D
∗
φ
aux
D
x
D
∗
=
0
.Define(127)
E
aux
P
=
φ
aux
P
x
P
-
φ
-
P
P
∈
P
,
E
aux
D
=
φ
aux
D
x
D
∗
-
φ
-
D
∗
D
∗
∈
D
.It is then clear that the vector(128)
E
aux
=
E
aux
P
,
E
aux
D
t
∈
R
P
,
Dis identifiable with a function denoted again by E
aux which lies in the space E
(
A
) and thus in the space R
P
,
D
.Proposition 25.
The functionE
a
u
x meets the following trivial properties:(129)
E
h
R
P
,
D
=
E
a
u
x
R
P
,
D
,
(130)
E
a
u
x
R
P
,
D
2
=
E
a
u
x
P
R
P
2
+
E
a
u
x
D
R
D
2
,
(131)
E
a
u
x
P
L
2
Ω
≤
C
E
a
u
x
P
R
P
,
(132)
E
a
u
x
D
L
2
Ω
≤
C
E
a
u
x
D
R
D
,where C stands for diverse positive numbers that are mesh independent.Let us look for the estimates of
∫
Ω
φ
(
x
)
-
Π
P
φ
-
h
(
x
)
2
d
x:(133)
∫
Ω
φ
x
-
Π
P
φ
-
h
x
2
d
x
≤
3
∫
Ω
φ
x
-
φ
aux
P
x
2
d
x
+
3
∫
Ω
φ
aux
P
x
-
Π
P
φ
aux
P
x
2
d
x
+
3
∫
Ω
Π
P
φ
aux
P
x
-
Π
P
φ
-
h
x
2
d
x
.(i)  Since the primary cellsC
P (with P
∈
P) are convex and the function φ
∈
C
2
P
- the Taylor-Lagrange expansion applies and gives rise to what follows: (134)
λ
P
=
φ
x
-
φ
aux
P
x
P
+
x
-
x
P
t
grad
⁡
φ
aux
P
ξ
P
in
P
,
∀
P
∈
P
.Since φ meets the null average condition (7) and φ
aux
P aux honors condition (125)-(ii), integrating the two sides of (134) in a primary cell P and summing on P
∈
P yield(135)
∫
Ω
φ
x
-
φ
aux
P
x
2
d
x
≤
C
h
2when accounting with the fact that g
r
a
d
φ
aux
P
x
≤
C a.e. in Ω, where C represents diverse positive numbers that are mesh independent.(ii)  Thanks again to Taylor-Lagrange it is easily seen that(136)
∫
Ω
φ
aux
P
x
-
Π
P
φ
aux
P
x
2
d
x
≤
C
h
2
.(iii)  At last, we have(137)
∫
Ω
Π
P
φ
aux
P
x
-
Π
P
φ
-
h
x
2
d
x
=
∑
P
∈
P
meas
⁡
C
P
φ
aux
P
x
P
-
φ
-
h
x
P
2
=
E
aux
P
L
2
Ω
2
≤
C
E
aux
P
R
P
2
due  to
131
≤
C
E
aux
P
R
P
,
D
2
according  to
130
≤
C
E
h
R
P
,
D
2
according  to
129
≤
C
h
2
by  virtue  of
120
.Summarizing what precedes we have proven the following.Lemma 26.
One has(138)
∫
Ω
φ
x
-
Π
P
φ
-
h
x
2
d
x
≤
C
h
2
.Remarking that any dual cellC
D
∗ (with D
∗
∈
D) is actually a finite union of convex homogeneous polygons sharing D
∗ as a common vertex, the arguments that have led to the preceding lemma apply and give rise to the following.Lemma 27.
One has(139)
∫
Ω
φ
x
-
Π
D
φ
-
h
x
2
d
x
≤
C
h
2
.Let us summarize the previous developments in the following proposition in whichN
0 denotes the set of nodes (i.e., cell-points and vertices) with respect to the mesh system P
,
E
,
A
,
D defined on Ω.Proposition 28 (error estimate result).
Assume that the primary mesh is regular in the sense of Definition6 and that the discontinuities in Ω of the piecewise constant permeability tensor K generate a finite number of convex subdomains Ω
s
s
∈
S over which the exact solution φ of (1)-(2) meets the following property:(140)
φ
Ω
s
∈
C
2
Ω
-
s
∀
s
∈
S
.Under conditions (3), (4), (6), (11), and (111), the error-function E
h associated with the error-vector from R
P
,
D with components E
N
=
φ
N
-
φ
-
N
N
∈
N
0 satisfies the following estimate:(141)
E
h
R
P
,
D
2
+
E
h
L
2
Ω
2
≤
C
h
2
,where C represents mesh independent real number.


## 3.1. Preliminaries and a Stability Result

Let us consider the auxiliary meshA introduced in the previous section (see Figure 4). Note that each mesh cell of A contains either one cell-point or one corner point and only one which should be lying inside or on the boundary of Ω. In the sequel, anode is a cell-point or a corner point with respect to the primary mesh.Definition 20.
An auxiliary mesh cell is degenerate if the corresponding node is lying on the boundary ofΩ
.In the sequel, we will say sometimes “auxiliary cell” instead of “auxiliary mesh cell.” We denote byE
(
A) the subspace of R
P
,
D made of functions v which satisfies conditions (72). This space is obviously nonempty as there is the null function. Moreover, the solution of the discrete system (25), (26), and (72) could clearly be identified with one (and only one) function from E
(
A). In the sequel we denote by φ
-
h such a function named “cellwise constant (approximate) solution” of the diffusion problem (1)-(2). Let us endow R
P
,
D with the following discrete seminorm.For allv
≡
V
P
P
∈
P
,
V
D
∗
D
∗
∈
D
∈
R
P
,
D, define(75)
v
R
P
,
D
=
∑
P
,
L
∈
P
2
,
A
∗
,
B
∗
∈
D
2
s.t.
Γ
P
∩
Γ
L
=
A
∗
B
∗
V
P
-
V
L
2
+
V
B
∗
-
V
A
∗
2
+
∑
A
∗
,
B
∗
∈
D
2
s.t.
A
∗
B
∗
∈
E
e
x
t
V
B
∗
-
V
A
∗
2
1
/
2
∀
V
∈
R
P
,
D
,where “s.t.” means “such that.” As this seminorm is actually a discrete energy norm for the space E
(
A), we denote it by ·
A as far as functions from E
(
A) are concerned. A norm on the space R
P
,
D is defined by the mapping (76)
·
R
P
,
D
=
·
R
P
,
D
2
+
·
L
2
Ω
2
1
/
2
.We focus here on the proof of the stability of the solution to the system of (25), (26), and (72). For this purpose, we need to introduce as in [31] a preliminary result, namely, a discrete version of Poincaré inequality which is based upon the following ingredients. Consider the following spaces(77)
R
P
=
v
=
v
P
P
∈
P
;
v
P
∈
R
∀
P
∈
P
,
R
D
=
v
=
v
D
∗
D
∗
∈
D
;
v
D
∗
∈
R
∀
D
∗
∈
Dand the linear operators Π
P and Π
D defined as follows: (78)
v
∈
E
A
⟼
Π
P
v
∈
E
P
,
v
∈
E
A
⟼
Π
D
v
∈
E
Dwith Π
P
v
=
v
P
P
∈
P and Π
D
v
=
v
D
∗
D
∗
∈
D, where E
P and E
D are, respectively, subspaces of R
P and R
D made of functions satisfying conditions (72)-(i) and (72)-(ii), respectively.Let us equip the function spacesE
P and E
D with the following discrete energy norms:(79)
v
P
=
∑
P
,
L
∈
P
2
,
A
∗
,
B
∗
∈
D
2
s.t.
Γ
P
∩
Γ
L
=
A
∗
B
∗
v
P
-
v
L
2
1
/
2
,
v
D
=
∑
P
,
L
∈
P
2
,
A
∗
,
B
∗
∈
D
2
s.t.
Γ
P
∩
Γ
L
=
A
∗
B
∗
v
B
∗
-
v
A
∗
2
+
∑
A
∗
,
B
∗
∈
D
2
s.t.
A
∗
B
∗
∈
E
e
x
t
V
B
∗
-
V
A
∗
2
1
/
2
.These mappings are only seminorms, respectively, for R
P and R
D that one can transform into norms for these spaces as follows:(80)
·
R
P
=
·
R
P
2
+
·
L
2
Ω
2
1
/
2
,
·
R
D
=
·
R
D
2
+
·
L
2
Ω
2
1
/
2
,where we have set ·
R
P
=
·
P and ·
R
D
=
·
D. The following results are key ingredients for proving the stability of the discrete solution in the sense of the discrete energy norm defined previously on the space E
A.Lemma 21 (discrete versions of Poincaré inequality).
We have the following inequalities:(81)
·
L
2
Ω
2
≤
C
v
P
2
+
2
m
e
s
Ω
∫
Ω
v
x
d
x
2
∀
v
∈
R
P
,
(82)
·
L
2
Ω
2
≤
C
v
D
2
+
2
m
e
s
Ω
∫
Ω
v
x
d
x
2
∀
v
∈
R
D
,where C represents diverse strictly positive numbers mesh independent.Note that the right-hand side of inequalities (81) and (82) is, respectively, norms for R
P and R
D. Since Ω is bounded, the previous lemma permits seeing that these norms are equivalent to standard ones, namely, (80). For proving the preceding lemma one can use the same arguments as in [36].Lemma 22 (a key-result).
LetP
-
,
E
- be a regular mesh system (defined on) in the sense of Definition 6 and denote by D the corresponding dual mesh. For v
∈
E
(
P
-
), denote by v
P the value of v in the control volume P. Let γ
-
v be a discrete trace function defined a.e. (for the (d
-
1)-Lebesgue measure) by γ
-
v
=
v
P on Γ
P
∩
Γ, for all P
∈
P
-
e
x
t, where P
-
e
x
t denotes the set of mesh elements adjacent to the domain boundary Γ. Then(83)
γ
-
v
L
2
Γ
≤
C
v
P
-
+
·
L
2
Ω
∀
v
∈
E
P
-
.Similarly, we have(84)
γ
-
v
L
2
Γ
≤
C
v
D
-
+
·
L
2
Ω
∀
v
∈
E
D
-
,where C stands for diverse positive numbers mesh independent.Proof.
See, for instance, [36].Remark 23.
It is more than useful to note that(85)
v
A
2
=
v
P
P
2
+
v
D
D
2
∀
v
∈
E
A
.Let us give now one of the main results of this section.Proposition 24 (stability result).
The cellwise constant approximate solutionφ
-
h of the diffusion problem (1)-(2) satisfies the following inequality:(86)
φ
-
h
A
≤
C
f
L
2
Ω
+
g
L
2
Γ
,where C is a strictly positive real number not depending on the spatial discretization.Proof.
Multiplying (25) by φ
-
P and (26) by φ
-
D
∗ and summing over P
∈
P and D
∗
∈
D, respectively, lead to(87)
Φ
c
c
Φ
v
c
A
B
B
t
C
Φ
c
c
Φ
v
c
=
Φ
c
c
Φ
v
c
F
c
c
F
v
c
.Recall that(88)
F
c
c
=
∫
C
P
f
x
d
x
-
∫
Γ
P
∩
Γ
g
d
γ
P
∈
P
,
F
v
c
=
∫
C
B
∗
f
x
d
x
-
∫
C
B
∗
∩
Γ
g
d
γ
+
∑
A
∗
∈
N
B
∗
,
I
∈
E
B
∗
∩
E
e
x
t
with
Γ
C
A
∗
∩
Γ
C
B
∗
=
P
I
c
h
K
E
h
I
E
a
K
E
h
B
∗
D
∗
∫
B
∗
D
∗
g
d
γ
B
∗
∈
D
.Let us set (89)
LHS
=
Φ
c
c
Φ
v
c
A
B
B
t
C
Φ
c
c
Φ
v
c
,
RHS
=
Φ
c
c
Φ
v
c
F
c
c
F
v
c
.We know from Lemma 16 that there exists a strictly positive number ρ mesh independent such that(90)
ρ
φ
-
h
A
2
≤
LHS
.On the other hand, according to Remark 11 we have(91)
RHS
=
∑
P
∈
P
∫
C
P
f
φ
-
P
d
x
-
∑
P
∈
P
∫
Γ
P
∩
Γ
g
φ
-
P
d
γ
+
∑
B
∗
∈
D
∫
C
B
∗
f
x
d
x
-
∑
B
∗
∈
D
∫
Γ
B
∗
∩
Γ
g
φ
-
B
∗
d
γ
+
∑
B
∗
∈
D
∑
A
∗
∈
N
B
∗
s.t.
,
∃
I
∈
E
B
∗
∩
E
A
∗
∩
E
e
x
t
,
∃
L
∈
P
,
with
Γ
A
∗
∩
Γ
B
∗
=
L
I
φ
-
B
∗
c
h
K
L
h
L
I
a
K
L
h
A
∗
B
∗
∫
A
∗
B
∗
g
d
γ
.Cauchy-Schwarz’s inequality leads to (92)
RHS
≤
f
L
2
Ω
φ
-
h
P
L
2
Ω
+
g
L
2
Γ
γ
-
P
φ
-
h
P
L
2
Γ
+
f
L
2
Ω
φ
-
h
D
L
2
Ω
+
g
L
2
Γ
γ
-
D
φ
-
h
D
L
2
Γ
+
∑
A
∗
,
B
∗
∈
D
e
x
t
s.t.
,
∃
I
∈
E
B
∗
∩
E
A
∗
,
∃
L
∈
P
,
with
Γ
A
∗
∩
Γ
B
∗
=
L
I
φ
-
B
∗
-
φ
-
A
∗
c
h
K
L
h
L
I
a
K
L
h
A
∗
B
∗
∫
A
∗
B
∗
g
d
γ
,where γ
-
P and γ
-
D are, respectively, discrete trace operators associated with grids P and D. Thanks to relation (4), Definition 6, and Cauchy-Schwarz inequality we get(93)
RHS
≤
f
L
2
Ω
φ
-
h
P
L
2
Ω
+
φ
-
h
D
L
2
Ω
+
g
L
2
Γ
γ
-
P
φ
-
h
P
L
2
Γ
+
γ
-
D
φ
-
h
D
L
2
Γ
+
1
ϖ
γ
m
i
n
sin
⁡
θ
∑
A
∗
,
B
∗
∈
D
e
x
t
s.t.
,
E
A
∗
∩
E
B
∗
≠
∅
φ
-
B
∗
-
φ
-
A
∗
∫
A
∗
B
∗
g
d
γ
.Accounting with the fact that φ
-
h satisfies null average condition on the grids P and D, Lemma 21 ensures that(94)
φ
-
h
P
L
2
Ω
+
φ
-
h
D
L
2
Ω
≤
C
φ
-
h
P
P
+
φ
-
h
D
D
,whereas Lemma 22 guarantees that(95)
γ
-
P
φ
-
h
P
L
2
Γ
+
γ
-
D
φ
-
h
D
L
2
Γ
≤
C
φ
-
h
P
P
+
φ
-
h
D
D
.A double application of Cauchy-Schwarz’s inequality leads to(96)
∑
A
∗
,
B
∗
∈
D
e
x
t
s.t.
,
E
A
∗
∩
E
B
∗
≠
∅
φ
-
B
∗
-
φ
-
A
∗
∫
A
∗
B
∗
g
d
γ
≤
meas
⁡
Γ
1
/
2
φ
-
h
D
D
g
L
2
Γ
.Combining relations (93), (94), (95), and (96) on the one hand and using the inequality of Cauchy-Schwarz and Remark 23 on the other hand yield(97)
RHS
≤
C
f
L
2
Ω
+
g
L
2
Γ
φ
-
h
A
.Recall that C represents diverse strictly positive numbers mesh independent. The stability result follows from relations (90) and (97). The proof of Proposition 24 ends.


## 3.2. Error Estimates

Following the original technique exposed in [28], we investigate in this subsection the error estimates for the finite volume approximate solution to the model problem. Recall that φ denotes the exact solution to the model problem while φ
-
h
∈
E
A is the cellwise constant approximate solution obtained from a DDFV formulation of this model problem. View also φ
-
h as a representation in terms of function of the vector Φ
=
Φ
c
c
Φ
v
c
t, where Φ
c
c
=
φ
-
P
P
∈
P and Φ
v
c
=
φ
-
D
∗
D
∗
∈
D. Recall that the vector Φ is defined above as the unique solution of the DDFV discrete system, satisfying null average conditions on grids P and D (see (25)-(26) and (72)). Let us set(98)
E
P
=
φ
x
P
-
φ
-
h
x
P
≡
φ
P
-
φ
-
P
∀
P
∈
P
,
E
D
∗
=
φ
x
D
∗
-
φ
-
h
x
D
∗
≡
φ
D
∗
-
φ
-
D
∗
∀
D
∗
∈
D
.Note that the components of the so-called error-vector E
P
P
∈
P
,
E
D
∗
D
∗
∈
D
t
∈
R
P
,
D can be viewed as values in auxiliary cells of the so-called error-function E
h defined a.e. in Ω. More precisely let us recall that any auxiliary cell A is associated with a unique cell-point or a unique vertex. Let E
A be the error corresponding to the auxiliary cell A; that is,(99)
E
A
=
φ
x
P
-
φ
-
h
x
P
if
P
is  a  cell  point  associated  with
A
φ
x
D
∗
-
φ
-
h
x
D
∗
if
D
∗
is  a  cell  point  associated  with
A
.Since the auxiliary mesh define a partition of Ω, one can set(100)
E
h
=
∑
A
∈
A
E
A
χ
A
x
,where χ
A is the characteristic function of any cell A.Our first purpose is to show that the error-vector is a solution to a square linear system of the same type as the discrete system (25)-(26) in the sense that both of them are associated with the same matrix M defined above. Note that unfortunately E
h
∈
R
P
,
D is not in E
(
A
) as it does meet the null average conditions (72) on grids P and D. So the error estimates, which are our main purpose in this section, should be investigated in the sense of the norm defined by (76) on the space R
P
,
D. In this connection we first should prove that(101)
E
h
R
P
,
D
≤
C
hand in the next step we should prove that(102)
E
h
L
2
Ω
≤
C
h
.(i) Let us start with proving that the components of the error-vector are solution to a square linear system of the same type as the discrete system (25)-(26).Accounting with truncation errors, one can see from the DDFV flux computations that the exact flux across the edgeA
∗
B
∗ viewed as part of the boundary of the primary cell C
P is given by(103)
Q
A
∗
B
∗
P
=
a
h
K
P
a
^
h
K
L
h
A
∗
B
∗
a
h
K
P
h
I
L
+
a
^
h
K
L
h
P
I
φ
P
-
φ
L
+
a
^
h
K
L
b
h
K
P
h
P
I
+
a
h
K
P
b
^
h
K
L
h
I
L
a
h
K
P
h
I
L
+
a
^
h
K
L
h
P
I
φ
B
∗
-
φ
A
∗
+
R
I
A
∗
,
B
∗
,
L
P
,where the truncation error R
I
A
∗
,
B
∗
,
L
P is given by(104)
R
I
A
∗
,
B
∗
,
L
P
=
1
+
a
h
K
P
h
I
L
a
h
K
P
h
I
L
+
a
^
h
K
L
h
P
I
T
I
A
∗
,
B
∗
P
+
a
h
K
P
h
I
L
a
h
K
P
h
I
L
+
a
^
h
K
L
h
P
I
T
I
A
∗
,
B
∗
Land meets the following inequality (see Proposition 7):(105)
R
I
A
∗
,
B
∗
,
L
P
≤
C
h
2
.On the other hand, the exact fluxes across the pseudoedge [
P
I
L
] for I
∈
Ω and across the edge [
P
I
] for I
∈
Γ, both of them acting as part of the boundary of the dual cell C
B
∗, are given by the following relations (accounting with truncation errors):(106)
Q
P
I
L
B
∗
=
c
h
K
P
a
^
h
K
L
h
P
I
+
c
^
h
K
L
a
h
K
P
h
I
L
a
h
K
P
h
I
L
+
a
^
h
K
L
h
P
I
φ
P
-
φ
L
+
ω
h
P
,
L
,
I
h
A
∗
B
∗
a
h
K
P
h
I
L
+
a
^
h
K
L
h
P
I
φ
B
∗
-
φ
A
∗
+
R
I
P
,
L
,
A
∗
B
∗
,
Q
P
I
B
∗
=
a
h
K
P
d
h
K
P
-
c
h
K
P
2
h
P
I
a
h
K
P
h
A
∗
B
∗
φ
B
∗
-
φ
A
∗
-
c
h
K
P
h
P
I
a
K
P
h
A
∗
B
∗
∫
A
∗
B
∗
g
d
γ
+
R
I
P
,
A
∗
B
∗
,where the truncation errors R
I
P
,
L
,
A
∗
B
∗ and R
I
(
P
,
A
∗
)
B
∗ meet the following inequalities (according to Propositions 7 and 10):(107)
R
I
P
,
L
,
A
∗
B
∗
≤
C
h
2
,
R
I
P
,
A
∗
B
∗
≤
C
h
2
.The system of equations satisfied by theexact nodal potentials φ
P
P
∈
P and φ
D
∗
D
∗
∈
D reads as(108)
∑
L
∈
N
P
s.t.
∃
I
∈
E
,
∃
A
∗
,
B
∗
∈
D
,
with
I
∈
E
P
∩
E
L
,
Γ
P
∩
Γ
L
=
A
∗
B
∗
a
h
K
P
a
^
h
K
L
h
A
∗
B
∗
a
h
K
P
h
I
L
+
a
^
h
K
L
h
P
I
φ
P
-
φ
L
+
a
^
h
K
L
b
h
K
P
h
P
I
+
a
h
K
P
b
^
h
K
L
h
I
L
a
h
K
P
h
I
L
+
a
^
h
K
L
h
P
I
φ
B
∗
-
φ
A
∗
=
∫
C
P
f
x
d
x
-
∫
Γ
P
∩
Γ
g
d
γ
-
∑
L
∈
N
P
s.t.
∃
I
∈
E
,
∃
A
∗
,
B
∗
∈
D
,
with
I
∈
E
P
∩
E
L
,
Γ
P
∩
Γ
L
=
A
∗
B
∗
R
I
A
∗
,
B
∗
,
L
P
∀
P
∈
P
,
∑
A
∗
∈
N
B
∗
s.t.
∃
I
∈
E
,
∃
P
,
L
∈
P
,
with
I
∈
E
A
∗
∩
E
B
∗
,
Γ
A
∗
∩
Γ
B
∗
=
P
I
L
c
h
K
P
a
^
h
K
L
h
P
I
+
c
^
h
K
L
a
h
K
P
h
I
L
a
h
K
P
h
I
L
+
a
^
h
K
L
h
P
I
φ
P
-
φ
L
+
ω
h
P
,
L
,
I
h
A
∗
B
∗
a
h
K
P
h
I
L
+
a
^
h
K
L
h
P
I
φ
B
∗
-
φ
A
∗
+
∑
A
∗
∈
N
B
∗
s.t.
∃
I
∈
E
,
∃
P
∈
P
with
I
∈
E
A
∗
∩
E
B
∗
∩
E
e
x
t
,
Γ
A
∗
∩
Γ
B
∗
=
P
I
a
K
P
d
h
K
P
-
c
h
K
P
2
h
P
I
a
h
K
P
h
A
∗
B
∗
φ
B
∗
-
φ
D
∗
=
∫
C
B
∗
f
x
d
x
-
∫
C
B
∗
∩
Γ
g
d
γ
-
∑
A
∗
∈
N
B
∗
s.t.
∃
I
∈
E
,
∃
P
,
L
∈
P
with
I
∈
E
A
∗
∩
E
B
∗
,
Γ
A
∗
∩
Γ
B
∗
=
P
I
L
R
I
P
,
L
,
A
∗
B
∗
-
∑
A
∗
∈
N
B
∗
s.t.
∃
I
∈
E
,
∃
P
∈
P
with
I
∈
E
A
∗
∩
E
B
∗
∩
E
e
x
t
,
Γ
A
∗
∩
Γ
B
∗
=
P
I
R
I
P
,
A
∗
B
∗
+
∑
A
∗
∈
N
B
∗
s.t.
∃
I
∈
E
,
∃
P
∈
P
with
I
∈
E
A
∗
∩
E
B
∗
∩
E
e
x
t
,
Γ
A
∗
∩
Γ
B
∗
=
P
I
c
h
K
P
h
P
I
a
K
P
h
A
∗
B
∗
∫
A
∗
B
∗
g
d
γ
∀
B
∗
∈
D
,where ω
h
P
,
L
,
I is defined in (23). Recall that N
· is the set of neighboring cells of a given (primary or dual) cell.Due to the conservation of our flux approximation schemes, the truncation errors naturally obey to the following relation:(109)
R
I
·
,
·
,
b
a
+
R
I
·
,
·
,
a
b
=
0
∀
I
∈
E
i
n
t
,
I
=
E
a
∩
E
b
,
(110)
R
I
·
,
b
a
+
R
I
·
,
a
b
=
0
∀
I
∈
E
e
x
t
,
I
=
E
a
∩
E
b
,where a and b are two cell-points or corner points associated with primary or dual adjacent cells. Note that, in the previous equality, the truncation errors are written in a very formal way. This is with a view to involving both primary and dual adjacent cells. Let us introduce now the set of diamond cells called in what follows the diamond mesh and denoted as M (the concept of diamond cell has been introduced for a different usage in earlier works of some authors: see, e.g., [22, 25, 34]). Each diamond cell is associated with one edge-point and vice versa (see Figure 6). A diamond cell is declared degenerate if the corresponding edge-point is lying on the boundary of Ω. An example of a degenerate diamond cell is provided in Figure 6.Figure 6
Examples of diamond cells. (a) A normal diamond cell. (b) A degenerate diamond cell.
(a)
(b)The following assumption plays a key-role in our proof of the convergence of the piecewise constant (approximate) solutionφ
-
h.(111)
∃
ν
∈
R
+
∗
such  that
ν
h
2
≤
meas
⁡
M
∀
P
∈
M
,where m
e
a
s
(
·
) is the Lebesgue measure in any spatial dimension.Recall that we have setE
P
=
φ
P
-
φ
-
P, for P
∈
P, and E
D
∗
=
φ
D
∗
-
φ
-
D
∗, for D
∗
∈
D. An adequate linear combination of equations from system (108)-(109) with those from system (25)-(26) shows that the quantities E
P
=
φ
P
-
φ
-
P
P
∈
P and E
D
∗
=
φ
D
∗
-
φ
-
D
∗
D
∗
∈
D satisfy the following equations: (112)
∑
L
∈
N
P
s.t.
∃
I
∈
E
,
∃
A
∗
,
B
∗
∈
D
,
with
I
∈
E
P
∩
E
L
,
Γ
P
∩
Γ
L
=
A
∗
B
∗
a
h
K
P
a
^
h
K
L
h
A
∗
B
∗
a
h
K
P
h
I
L
+
a
^
h
K
L
h
P
I
E
P
-
E
L
+
a
^
h
K
L
b
h
K
P
h
P
I
+
a
h
K
P
b
^
h
K
L
h
I
L
a
h
K
P
h
I
L
+
a
^
h
K
L
h
P
I
E
B
∗
-
E
A
∗
=
∑
L
∈
N
P
s.t.
∃
I
∈
E
,
∃
A
∗
,
B
∗
∈
D
,
with
I
∈
E
P
∩
E
L
,
Γ
P
∩
Γ
L
=
A
∗
B
∗
R
I
A
∗
,
B
∗
,
L
P
∀
P
∈
P
,
(113)
∑
A
∗
∈
N
B
∗
s.t.
∃
I
∈
E
,
∃
P
,
L
∈
P
,
with
I
∈
E
A
∗
∩
E
B
∗
,
Γ
A
∗
∩
Γ
B
∗
=
P
I
L
c
h
K
P
a
^
h
K
L
h
P
I
+
c
^
h
K
L
a
h
K
P
h
I
L
a
h
K
P
h
I
L
+
a
^
h
K
L
h
P
I
E
P
-
E
L
+
ω
h
P
,
L
,
I
h
A
∗
B
∗
a
h
K
P
h
I
L
+
a
^
h
K
L
h
P
I
E
B
∗
-
E
A
∗
+
∑
A
∗
∈
N
B
∗
s.t.
∃
I
∈
E
,
∃
P
∈
P
with
I
∈
E
A
∗
∩
E
B
∗
∩
E
e
x
t
,
Γ
A
∗
∩
Γ
B
∗
=
P
I
a
h
K
P
d
h
K
P
-
c
h
K
P
2
h
P
I
a
h
K
P
h
A
∗
B
∗
E
B
∗
-
E
A
∗
=
-
∑
A
∗
∈
N
B
∗
s.t.
∃
I
∈
E
,
∃
P
,
L
∈
P
with
I
∈
E
A
∗
∩
E
B
∗
,
Γ
A
∗
∩
Γ
B
∗
=
P
I
L
R
I
P
,
L
,
A
∗
B
∗
-
∑
A
∗
∈
N
B
∗
s.t.
∃
I
∈
E
,
∃
P
∈
P
with
I
∈
E
A
∗
∩
E
B
∗
∩
E
e
x
t
,
Γ
A
∗
∩
Γ
B
∗
=
P
I
R
I
P
,
A
∗
B
∗
∀
B
∗
∈
D
.(ii) Let us now prove thatE
h
R
P
,
D
≤
C
h, where C
>
0 is a mesh independent positive number.Multiplying (112) by E
P and (113) by E
B
∗ and summing over P
∈
P and B
∗
∈
D yield (thanks to Lemma 16)(114)
γ
E
h
R
P
,
D
2
≤
∑
P
∈
P
-
E
p
∑
L
∈
N
P
s.t.
∃
I
∈
E
P
∩
E
L
,
∃
A
∗
,
B
∗
∈
D
,
with
Γ
P
∩
Γ
L
=
A
∗
B
∗
R
I
A
∗
,
B
∗
,
L
P
+
∑
B
∗
∈
D
-
E
B
∗
∑
A
∗
∈
N
B
∗
s.t.
∃
I
∈
E
A
∗
∩
E
B
∗
,
∃
P
,
L
∈
P
with
Γ
A
∗
∩
Γ
B
∗
=
P
I
L
R
I
P
,
L
,
A
∗
B
∗
-
E
B
∗
∑
A
∗
∈
N
B
∗
s.t.
∃
I
∈
E
A
∗
∩
E
B
∗
∩
E
e
x
t
,
∃
P
∈
P
with
Γ
A
∗
∩
Γ
B
∗
=
P
I
R
I
P
,
A
∗
B
∗
,where γ is some strictly positive real number mesh independent.Thanks to (109)-(110) (consequence of the conservation of the flux approximation scheme) we have(115)
γ
E
h
R
P
,
D
2
≤
∑
I
∈
E
s.t.
∃
P
,
L
∈
P
with
I
∈
E
P
∩
E
L
,
∃
A
∗
,
B
∗
∈
D
with
Γ
P
∩
Γ
L
=
A
∗
B
∗
E
P
-
E
L
R
I
A
∗
,
B
∗
,
L
P
+
∑
I
∈
E
s.t.
∃
A
∗
,
B
∗
∈
D
with
I
∈
E
A
∗
∩
E
B
∗
,
∃
P
,
L
∈
P
with
Γ
A
∗
∩
Γ
B
∗
=
P
I
L
E
B
∗
-
E
A
∗
R
I
P
,
L
,
A
∗
B
∗
+
∑
I
∈
E
s.t.
∃
A
∗
,
B
∗
∈
D
with
I
∈
E
A
∗
∩
E
B
∗
∩
Γ
,
∃
P
∈
P
with
Γ
A
∗
∩
Γ
B
∗
=
P
I
E
B
∗
-
E
A
∗
R
I
P
,
A
∗
B
∗
.Therefore(116)
γ
E
h
R
P
,
D
2
≤
∑
I
P
,
L
,
A
∗
,
B
∗
∈
E
i
n
t
R
I
P
,
L
,
A
∗
,
B
∗
max
⁡
E
P
-
E
L
+
E
B
∗
-
E
A
∗
+
∑
I
P
,
A
∗
,
B
∗
∈
E
e
x
t
R
I
P
,
A
∗
,
B
∗
max
⁡
E
B
∗
-
E
A
∗
,where we have set(117)
R
I
P
,
L
,
A
∗
,
B
∗
max
⁡
=
max
⁡
R
I
A
∗
,
B
∗
,
L
P
,
R
I
P
,
L
,
A
∗
B
∗
,
R
I
P
,
A
∗
,
B
∗
=
R
I
P
,
A
∗
B
∗
.It is important to mention that, due to (105) and (107), we have(118)
0
≤
R
I
P
,
L
,
A
∗
,
B
∗
max
⁡
≤
C
h
2
,
0
≤
R
I
P
,
A
∗
,
B
∗
≤
C
h
2
.Define S
I
P
,
L
,
A
∗
,
B
∗
=
2
D Lebesgue measure of the diamond cell defined by the points P
,
L
,
A
∗
,
B
∗ and associated with I
∈
E
i
n
t and S
I
P
,
A
∗
,
B
∗
=
2
D Lebesgue measure of the diamond cell defined by the points P
,
A
∗
,
B
∗ and associated with I
∈
E
e
x
t.Therefore, thanks to Cauchy-Schwarz’s inequality, it follows from (116) that(119)
γ
E
h
R
P
,
D
2
≤
∑
I
P
,
L
,
A
∗
,
B
∗
∈
E
i
n
t
S
I
P
,
L
,
A
∗
,
B
∗
+
∑
I
P
,
A
∗
,
B
∗
∈
E
i
n
t
S
I
P
,
A
∗
,
B
∗
1
/
2
∑
I
P
,
L
,
A
∗
,
B
∗
∈
E
i
n
t
2
R
I
P
,
L
,
A
∗
,
B
∗
max
⁡
2
S
I
P
,
L
,
A
∗
,
B
∗
E
P
-
E
L
2
+
E
B
∗
-
E
A
∗
2
1
/
2
+
∑
I
P
,
A
∗
,
B
∗
∈
E
e
x
t
2
R
I
P
,
A
∗
,
B
∗
2
S
I
P
,
A
∗
,
B
∗
E
B
∗
-
E
A
∗
2
1
/
2
.One concludes with the help of assumption (111) and relation (118) that(120)
E
h
R
P
,
D
≤
C
h
,where C is a strictly positive constant that is mesh independent.(iii) Let us now prove thatE
h
L
2
Ω
≤
C
h.We start with setting(121)
A
P
=
E
∈
A
;
∃
P
∈
P
with
x
P
∈
E
,
A
D
=
E
∈
A
;
∃
D
∗
∈
D
with
x
D
∗
∈
E
.Recall A is the auxiliary grid and x
P is a cell-point while x
D
∗ is a vertex (with respect to the primary grid introduced in a preceding section).Define(122)
Ω
P
=
⋃
E
∈
A
P
E
,
Ω
D
=
⋃
E
∈
A
D
E
.Note that Ω
P and Ω
D are a partition of Ω in the sense that(123)
i
Ω
P
∩
Ω
D
=
∅
,
i
i
Ω
-
=
Ω
-
P
∪
Ω
-
D
.Two main steps are necessary in our technique to get the estimates of φ
-
φ
-
L
2
Ω; before exposing these steps, we develop some preliminaries.(124)
∫
Ω
φ
x
-
φ
-
h
x
2
d
x
=
∫
Ω
P
φ
x
-
φ
-
h
x
2
d
x
+
∫
Ω
D
φ
x
-
φ
-
h
x
2
d
x
≤
∫
Ω
φ
x
-
Π
P
φ
-
h
x
2
d
x
+
∫
Ω
φ
x
-
Π
D
φ
-
h
x
2
d
x
,where Π
P and Π
D are operators from the function space E
(
A
) into the function spaces A
P and A
D, respectively (see (78) for the definition of Π
P and Π
D).Letλ
P and λ
D be two real numbers and φ
aux
P and φ
aux
D two real functions defined almost everywhere in Ω such that(125)
i
φ
aux
P
=
φ
x
-
λ
P
,
i
i
∑
P
∈
P
meas
⁡
C
P
φ
aux
P
x
P
=
0
,
(126)
i
φ
aux
D
=
φ
x
-
λ
D
,
i
i
∑
D
∗
∈
D
meas
⁡
C
D
∗
φ
aux
D
x
D
∗
=
0
.Define(127)
E
aux
P
=
φ
aux
P
x
P
-
φ
-
P
P
∈
P
,
E
aux
D
=
φ
aux
D
x
D
∗
-
φ
-
D
∗
D
∗
∈
D
.It is then clear that the vector(128)
E
aux
=
E
aux
P
,
E
aux
D
t
∈
R
P
,
Dis identifiable with a function denoted again by E
aux which lies in the space E
(
A
) and thus in the space R
P
,
D
.Proposition 25.
The functionE
a
u
x meets the following trivial properties:(129)
E
h
R
P
,
D
=
E
a
u
x
R
P
,
D
,
(130)
E
a
u
x
R
P
,
D
2
=
E
a
u
x
P
R
P
2
+
E
a
u
x
D
R
D
2
,
(131)
E
a
u
x
P
L
2
Ω
≤
C
E
a
u
x
P
R
P
,
(132)
E
a
u
x
D
L
2
Ω
≤
C
E
a
u
x
D
R
D
,where C stands for diverse positive numbers that are mesh independent.Let us look for the estimates of
∫
Ω
φ
(
x
)
-
Π
P
φ
-
h
(
x
)
2
d
x:(133)
∫
Ω
φ
x
-
Π
P
φ
-
h
x
2
d
x
≤
3
∫
Ω
φ
x
-
φ
aux
P
x
2
d
x
+
3
∫
Ω
φ
aux
P
x
-
Π
P
φ
aux
P
x
2
d
x
+
3
∫
Ω
Π
P
φ
aux
P
x
-
Π
P
φ
-
h
x
2
d
x
.(i)  Since the primary cellsC
P (with P
∈
P) are convex and the function φ
∈
C
2
P
- the Taylor-Lagrange expansion applies and gives rise to what follows: (134)
λ
P
=
φ
x
-
φ
aux
P
x
P
+
x
-
x
P
t
grad
⁡
φ
aux
P
ξ
P
in
P
,
∀
P
∈
P
.Since φ meets the null average condition (7) and φ
aux
P aux honors condition (125)-(ii), integrating the two sides of (134) in a primary cell P and summing on P
∈
P yield(135)
∫
Ω
φ
x
-
φ
aux
P
x
2
d
x
≤
C
h
2when accounting with the fact that g
r
a
d
φ
aux
P
x
≤
C a.e. in Ω, where C represents diverse positive numbers that are mesh independent.(ii)  Thanks again to Taylor-Lagrange it is easily seen that(136)
∫
Ω
φ
aux
P
x
-
Π
P
φ
aux
P
x
2
d
x
≤
C
h
2
.(iii)  At last, we have(137)
∫
Ω
Π
P
φ
aux
P
x
-
Π
P
φ
-
h
x
2
d
x
=
∑
P
∈
P
meas
⁡
C
P
φ
aux
P
x
P
-
φ
-
h
x
P
2
=
E
aux
P
L
2
Ω
2
≤
C
E
aux
P
R
P
2
due  to
131
≤
C
E
aux
P
R
P
,
D
2
according  to
130
≤
C
E
h
R
P
,
D
2
according  to
129
≤
C
h
2
by  virtue  of
120
.Summarizing what precedes we have proven the following.Lemma 26.
One has(138)
∫
Ω
φ
x
-
Π
P
φ
-
h
x
2
d
x
≤
C
h
2
.Remarking that any dual cellC
D
∗ (with D
∗
∈
D) is actually a finite union of convex homogeneous polygons sharing D
∗ as a common vertex, the arguments that have led to the preceding lemma apply and give rise to the following.Lemma 27.
One has(139)
∫
Ω
φ
x
-
Π
D
φ
-
h
x
2
d
x
≤
C
h
2
.Let us summarize the previous developments in the following proposition in whichN
0 denotes the set of nodes (i.e., cell-points and vertices) with respect to the mesh system P
,
E
,
A
,
D defined on Ω.Proposition 28 (error estimate result).
Assume that the primary mesh is regular in the sense of Definition6 and that the discontinuities in Ω of the piecewise constant permeability tensor K generate a finite number of convex subdomains Ω
s
s
∈
S over which the exact solution φ of (1)-(2) meets the following property:(140)
φ
Ω
s
∈
C
2
Ω
-
s
∀
s
∈
S
.Under conditions (3), (4), (6), (11), and (111), the error-function E
h associated with the error-vector from R
P
,
D with components E
N
=
φ
N
-
φ
-
N
N
∈
N
0 satisfies the following estimate:(141)
E
h
R
P
,
D
2
+
E
h
L
2
Ω
2
≤
C
h
2
,where C represents mesh independent real number.


## 4. Numerical Test

The triangular coarse and fine meshes (see Figure7) are also from FVCA5 Benchmark (see, e.g., [37]). Let us consider a diffusion problem formulated as (1)-(2), where the permeability tensor K is defined by the following relation:(142)
K
=
1
0
0
10
5
.The exact solution is taken to be u
x
,
y
=
sin
⁡
(
2
π
x
)
e
-
2
π
y
1
/
10
5. Furthermore, we shall ensure the uniqueness of the solution to (1)-(2) by enforcing the condition(143)
∫
Ω
φ
x
d
x
=
0
.
Notations
(i)⁡
n
u
n
k
w: number of unknowns.
(ii)⁡
n
n
m
a
t: number of nonzero terms in the matrix.
(iii)⁡
s
u
m
f
l
u
x: the discrete flux balance; that is, s
u
m
f
l
u
x
=
f
l
u
x
0
+
f
l
u
x
1
+
f
l
u
y
0
+
f
l
u
y
1, where f
l
u
x
0, f
l
u
x
1, f
l
u
y
0, and f
l
u
y
1 are, respectively, the outward numerical fluxes through the boundaries x
=
0, x
=
1, y
=
0, and y
=
1 (e.g., flux0 is an approximation of ∫
x
=
0
K
∇
u
·
η
d
s), and sum
f
=
∑
K
∈
T
K
f
(
x
K
), where x
K denotes some point of the control volume K; note that the residual s
u
m
f
l
u
x is a measure of the global conservation of the scheme.
(iv)⁡
u
m
i
n: value of the minimum of the approximate solution.
(v)⁡
u
m
a
x: value of the maximum of the approximate solution.
(vi)⁡
e
n
e
r
1, e
n
e
r
2: approximations of the energy following the two expressions:E
1
=
∫
Ω
K
∇
u
·
∇
u
d
x
,
E
1
=
∫
Ω
K
∇
u
·
η
u
d
x
.
Let us denote byu the exact solution, by T the mesh, and by u
T
=
u
K
K
∈
T the piecewise constant approximate solution.
(vii)⁡
e
r
l
2: relative discrete L
2-norm of the error; that is, for instance,(144)
e
r
l
2
=
∑
P
∈
P
meas
⁡
C
P
φ
x
P
-
φ
-
h
x
P
2
∑
P
∈
P
meas
⁡
C
P
φ
x
P
2
1
/
2
.
(viii)⁡
e
r
g
r
a
d: discrete L
2-norm of the error on the gradient; that is, for instance,(145)
ergrad
=
∑
σ
∈
ξ
P
i
n
t
m
σ
d
K
L
φ
-
h
x
P
-
φ
-
h
x
L
2
+
∑
σ
∗
∈
ξ
D
i
n
t
m
σ
∗
d
A
∗
B
∗
φ
-
h
x
A
∗
-
φ
-
h
x
B
∗
2
1
/
2
.
(ix)⁡
r
a
t
i
o
l
2: for i
≥
2
,
(146)
r
a
t
i
o
l
2
i
=
-
2
ln
⁡
e
r
l
2
i
⁡
-
ln
⁡
e
r
l
2
i
-
1
⁡
ln
⁡
n
u
n
k
w
i
⁡
-
ln
⁡
n
u
n
k
w
i
-
1
.
(x)⁡
ratiograd: for i
≥
2, same formula as above with ergrad instead of erl2.
(xi)⁡
e
r
f
l
x
0, e
r
f
l
x
1, e
r
f
l
y
0, e
r
f
l
y
1: relative error between f
l
u
x
0, f
l
u
x
1, f
l
u
y
0, f
l
u
y
1 and the corresponding flux of the exact solution:(147)
e
r
f
l
x
0
=
f
l
u
x
0
+
∫
x
=
0
K
∇
u
·
η
d
s
∫
x
=
0
K
∇
u
·
η
d
s
.
(xii)⁡
o
c
v
l
2: order of convergence of the method for L
2-norm of the solution as defined by erl2 with respect to the mesh size:(148)
o
c
v
l
2
=
l
n
⁡
e
r
l
2
i
max
⁡
-
l
n
⁡
e
r
l
2
i
max
⁡
-
1
l
n
⁡
h
i
max
⁡
-
ln
⁡
h
i
max
⁡
-
1
,
whereh is the maximum of the diameter of the control volume.
(xiii)⁡
o
c
v
e
n
e
r
d
i
s
c: order of convergence of the method for the norm defined by Remark 23.
(xiv)⁡
ocvgradl2: order of convergence of the method inL
2-norm of the gradient as defined by ergrad with respect to the mesh size, same formula as above with ergrad instead of erl2.
Comments about the Numerical Test Results. The results of numerical computations of pressure and its gradient show a convergence of order two in L
2-norm and a convergence of order 1.4 in discrete energy norm (see, e.g., Table 1). Moreover the similarity of curves erl2 and ergrad (Figure 8) confirm the closeness of their order of convergence. There is no discordance with the theoretical result where a linear convergence is announced (see Proposition 28). Indeed, the presence of discontinuities in the permeability tensor coefficients prevents the exact solution from being regular enough in the whole domain. More precisely in presence of such discontinuities, the exact pressure is in H
1
Ω and never in H
2
Ω no matter how regular may be the data f
,
g and the domain boundary. So a much slower convergence should take place.Table 1
Numerical results.
(a)
i
nunkw
nnmat
sumflux
erl2
ergrad
Ratiol2
ratiograd
1
21
68
−1.16E − 10
4.7E − 01
1.62E − 01
0E00
0E00
2
109
944
−6.03E − 11
1.09E − 01
2.01E02
1.78E00
−8.66E00
3
385
3576
3.50E − 10
2.12E − 02
3.80E01
2.60E00
2.64E00
4
1441
13880
−1.02E − 10
4.82E − 03
9.70E00
2.24E00
2.07E00
5
5569
54648
−6.55E − 11
1.18E − 03
2.45E00
2.09E00
2.03E00
(b)
i
erflx0
erflx1
erfly0
erfly1
erflm
umin
umax
1
1.43E − 16
1.43E − 16
5.13E02
5.11E01
4.14E − 01
−1.79E00
1.79E00
2
2.62E − 13
6.65E − 14
4.1E02
9.22E03
4.59E02
−9.42E − 01
9.23E − 01
3
1.91E − 13
1.37E − 13
1.23E04
4.61E03
8.56E01
−9.88E − 01
9.88E − 01
4
−0E00
0E00
6.14E03
2.78E01
−9.97E − 01
−9.98E − 01
9.97E − 01
5
−0E00
0E00
7.23E02
7.37E00
−9.99E − 01
−9.99E − 01
9.99E − 01
(c)
i
flux0
flux1
fluy0
fluy1
ener1
ener2
eren
enerdisc
1
6.22E00
−6.22E00
−5.82E − 11
−5.82E − 11
4.15E01
1.12E00
7.3E − 01
2.9905
2
6.22E00
−6.22E00
4.66E − 10
−5.24E − 10
3.03E01
2.46E01
1.89E − 01
1.2247
3
6.22E00
−6.22E00
−1.75E − 10
5.24E − 10
3.62E01
3.30E01
8.86E − 02
0.3368
4
6.22E00
−6.22E00
−1.31E − 10
2.91E − 11
3.80E01
3.69E01
2.83E − 02
0.1046
5
6.22E00
−6.22E00
1.24E − 10
−1.89E − 10
3.85E01
3.82E01
7.86E − 03
0.0389
(d)
ocvl2
ocvgradl2
ocvenerdisc
2.04E00
1.98E00
1.42E00Figure 7
Triangular mesh with acute angles. (a) Coarse mesh and (b) fine mesh.
(a)
(b)Figure 8
Relative discreteL
2-norm of the error (erl2) and of the error on the gradient (ergrad).



---

*Source: 5891064-2016-09-18.xml*