A second order multipoint flux mixed finite element method on hybrid meshes

TL;DR

This paper develops a second order mixed finite element method with mass lumping for simulating single phase flow in porous media on hybrid meshes, extending previous work with new analysis and computational validation.

Contribution

It introduces a second order approximation method on hybrid meshes with a comprehensive convergence analysis and new quadrature rules, advancing numerical techniques for porous media flow.

Findings

Method achieves second order accuracy in hybrid meshes.

Convergence analysis confirms theoretical predictions.

Computational tests validate the method's effectiveness.

Abstract

We consider the numerical approximation of single phase flow in porous media by a mixed finite element method with mass lumping. Our work extends previous results of Wheeler and Yotov, who showed that mass lumping together with an appropriate choice of basis allows to eliminate the flux variables locally and to reduced the mixed problem in this way to a finite volume discretization for the pressure only. Here we construct second order approximations for hybrid meshes in two and three space dimensions which, similar to the method of Wheeler and Yotov, allows the local elimination of the flux variables. A full convergence analysis of the method is given for which new arguments and, in part, also new quadrature rules and finite elements are required. Computational tests are presented for illustration of the theoretical results.

Tables

Table 1. Table 7.1. Degrees of freedom, relative discretization errors, and convergence rates for the first order multipoint flux finite element method.

$h$	DOF $u$	DOF $p$	$\Vert u-u_h\Vert$	eoc	$\Vert {\pi }_h^{0}p-p_h\Vert$	eoc
$2^{-1}$	$108$	$28$	$0.396017$	—	$0.238983$	—
$2^{-2}$	$460$	$132$	$0.184622$	$1.10$	$0.057068$	$2.06$
$2^{-3}$	$1574$	$462$	$0.092531$	$0.99$	$0.015148$	$1.91$
$2^{-4}$	$6094$	$1822$	$0.046019$	$1.00$	$0.003804$	$1.99$
$2^{-5}$	$25050$	$7590$	$0.022949$	$1.00$	$0.000945$	$2.00$

Table 2. Table 7.2. Degrees of freedom, relative discretization errors, and convergence rates for the second order multipoint flux finite element method.

$h$	DOF $u$	DOF $p$	$\Vert u-u_h\Vert$	eoc	$\Vert p-p_h\Vert$	eoc	$\Vert {\pi }_h^0\left(p-p_h\right)\Vert$	eoc
$2^{-1}$	$164$	$84$	$0.078309$	—	$0.064026$	—	$0.033106$	—
$2^{-2}$	$724$	$396$	$0.013097$	$2.57$	$0.007741$	$3.04$	$0.002864$	$3.53$
$2^{-3}$	$2498$	$1386$	$0.002275$	$2.52$	$0.001936$	$1.99$	$0.000391$	$2.87$
$2^{-4}$	$9738$	$5466$	$0.000484$	$2.23$	$0.000469$	$2.04$	$0.000049$	$2.99$
$2^{-5}$	$40230$	$22770$	$0.000099$	$2.28$	$0.000111$	$2.06$	$0.000005$	$3.13$

Equations

$${\mathcal{K}}^{-1}u+\nabla p$$

$$\mathrm{div}u$$

$$p=g\text{on }\partial \Omega ,$$

$$-\mathrm{div}(\mathcal{K}\nabla p)$$

$$p$$

$$({\mathcal{K}}^{-1}{u}_h,v_h{)}_h-(p_h,\mathrm{div}{v}_h)$$

$$(\mathrm{div}{u}_h,q_h)$$

$${\mathsf{M}}_{\mathsf{h}}\mathsf{u}-{\mathsf{B}}^{\top }\mathsf{p}$$

$$\mathsf{Bu}$$

$$\mathsf{B}{\mathsf{M}}_{\mathsf{h}}^{-1}{\mathsf{B}}^{\top }\mathsf{p}$$

$$\Vert u-u_h{\Vert }_{L^2(\Omega )}+\Vert p-p_h{\Vert }_{L^2(\Omega )}\le Ch,$$

$$\Vert {\pi }_h^{0}p-p_h{\Vert }_{L^2(\Omega )}\le C{h}^2,$$

$$\Vert u-u_h{\Vert }_{L^2(\Omega )}+\Vert {\pi }_h^0(p-p_h){\Vert }_{L^2(\Omega )}+\Vert p-{\tilde{p}}_h{\Vert }_{L^2(\Omega )}\le C{h}^2,$$

$$\Vert {\pi }_h^0(p-p_h){\Vert }_{L^2(\Omega )}+\Vert p-{\tilde{p}}_h{\Vert }_{L^2(\Omega )}\le C{h}^3,$$

$$\underline{k}|\xi {|}^2\le {\xi }^{\top }\mathcal{K}(x)\xi \le \overline{k}|\xi {|}^2,$$

$$\Vert u{\Vert }_{H(\mathrm{div})}+\Vert p{\Vert }_{H^1(\Omega )}\le C(\Vert f{\Vert }_{L^2(\Omega )}+|g{|}_{H^{1/2}(\partial \Omega )}).$$

$$({\mathcal{K}}^{-1}u,v)-(p,\mathrm{div}v)$$

$$(\mathrm{div}u,q)$$

$$\Vert v{\Vert }_{H(\mathrm{div},\Omega )}+\Vert v{\Vert }_{L^p(\Omega )}\le C\Vert q{\Vert }_{L^2(\Omega )}.$$

$$F_T(\hat{x})=a_T+B_T\hat{x}\text{with }{a}_T\in {\mathbb{R}}^d,B_T\in {\mathbb{R}}^{d\times d}.$$

$$\Vert B_T\Vert \le c_{\gamma }{h}_T\text{and}\Vert B_T^{-1}\Vert \le c_{\gamma }{h}_T^{-1},$$

$$({\pi }_{h}p,q_h)=(p,q_h)\forall q_h\in Q_h,$$

$$\mathrm{div}{\Pi }_{h}v={\pi }_h\mathrm{div}v\text{for all }v\in H(\mathrm{div},\Omega )\cap L^p(\Omega {)}^d.$$

$$\begin{array}{rl}(u_h,v_h{)}_h& :=\sum _{T\in {\mathcal{T}}_h}(u_h,v_h{)}_{h,T}=\sum _{T\in {\mathcal{T}}_h}|T|\sum _{n=1}^{N_T}(u_h\cdot v_h)(x_{n,T})w_{n,T},\end{array}$$

$$c^{-1}(v_h,v_h)\le (v_h,v_h{)}_h\le c(v_h,v_h),$$

$${\hat{c}}^{-1}(v_h,v_h)\le ({\mathcal{K}}^{-1}{v}_h,v_h{)}_h\le \hat{c}(v_h,v_h),$$

$$({\mathcal{K}}^{-1}{w}_h,v_h{)}_h-(r_h,\mathrm{div}{v}_h)$$

$$(\mathrm{div}{w}_h,q_h)$$

$$\Vert w_h{\Vert }_{H(\mathrm{div})}+\Vert r_h{\Vert }_{L^2(\Omega )}\le C\left(\underset{v_h\in V_h}{\sup }\frac{a_h(v_h)}{\Vert v_h{\Vert }_{H(\mathrm{div})}}+\underset{q_h\in Q_h}{\sup }\frac{b_h(q_h)}{\Vert q_h{\Vert }_{L^2(\Omega )}}\right)$$

$$({\mathcal{K}}^{-1}{w}_h,v_h{)}_h\le \overline{k}c\Vert w_h{\Vert }_{H(\mathrm{div})}\Vert v_h{\Vert }_{H(\mathrm{div})},\forall w_h,v_h\in V_h$$

$$\underset{v_h\in V_h}{\sup }\frac{(\mathrm{div}{v}_h,q_h)}{\Vert v_h{\Vert }_{H(\mathrm{div})}}\ge \beta \Vert q_h{\Vert }_{L^2(\Omega )},\text{for all }{q}_h\in Q_h,$$

$$({\mathcal{K}}^{-1}{v}_h,v_h{)}_h$$

Full text

**Keywords: porous medium equation, mixed finite elements, mass-lumping, multipoint flux method **

**AMS-classification (2000): 35F45, 65N30, 76M10 **

1. Introduction

We consider the numerical approximation of single phase flow through a saturated porous medium modeled by the Darcy law and the continuity equation

[TABLE]

Here $\Omega$ is the computational domain, $u$ denotes the flow velocity, $p$ represents the pressure in the fluid, $\mathcal{K}$ is the the hydraulic conductivity tensor, and $f$ are the source terms. For ease of notation, we assume that the pressure is known on the whole boundary, i.e.

[TABLE]

but other types of boundary conditions could be considered with minor modifications.

The simulation of problems of the form (1.1)–(1.3) is of practical relevance, i.e., in oil recovery or groundwater hydrology [10, 24]. In such applications, the conductivity $\mathcal{K}$ is usually anisotropic and spatially inhomogeneous, and mixed finite element methods, which have a local conservation property and which perform well for rough and anisotropic coefficients, seem particularly well-suited for practical computations.

A major drawback of mixed finite element approximations of (1.1)–(1.3) is that they naturally lead to algebraic saddle point problems that are larger and more difficult to solve than the symmetric positive definite systems obtained by more standard discretization schemes applied to the reduced problem

[TABLE]

which results from elimination of the velocity variable $u$ in (1.1)–(1.3). An elegant approach that allows to keep the advantages of the mixed finite element approximation and, at the same time, to eliminate the difficulties arising from the saddle-point structure, was proposed by Wheeler and Yotov [28]. They use special quadrature rules and appropriate basis functions for the approximation of the velocity, and consider a discrete variational approximation of the following form: Find $(u_{h},p_{h})\in V_{h}\times Q_{h}$ such that

[TABLE]

Here $(\cdot,\cdot)_{h}$ denotes an approximation of the $L^{2}$–scalar product $(\cdot,\cdot)$ over $\Omega$ which is obtained by appropriate numerical integration on every element. The degrees of freedom and quadrature rules considered in [28] are depicted in Figure 1.1.

A particular choice of basis functions for the spaces $V_{h}$ and $Q_{h}$ leads to a linear algebraic system of the form

[TABLE]

in which the mass matrix $\mathsf{M}_{\mathsf{h}}$ for the velocity variable is block diagonal. This allows to locally eliminate $\mathsf{u}$ from the system, very similar to the continuous level, and to obtain a reduced problem for the pressure variable of the form

[TABLE]

By the block diagonal structure, the mass matrix $\mathsf{M}_{\mathsf{h}}$ has a sparse inverse, and therefore the system matrix $\mathsf{B}\mathsf{M}_{\mathsf{h}}^{-1}\mathsf{B}^{\top}$ is sparse and has a compact stencil. Since the vector $\mathsf{p}$ here represents local cell averages of the pressure, the system (1.10) can be interpreted as a finite volume discretization of the reduced problem (1.3)–(1.4), which is closely related to certain finite difference and finite volume methods based on multipoint flux approximations [1, 21]. Due to this connection, the above approach was called a multipoint flux mixed finite element method in [28]; we refer to [4, 25, 26] for a detailed discussion of different discretization methods for subsurface flow and their mutual relations.

Besides the algebraic properties mentioned above, the paper [28] also provides a detailed numerical analysis of the approach, including error estimates of first order

[TABLE]

which are optimal in view of the pressure approximation by piecewise constants. By duality arguments, one can even obtain super-convergence for the pressure

[TABLE]

provided that the domain $\Omega$ is convex and that the conductivity $\mathcal{K}$ is sufficiently smooth. Here $\pi_{h}^{0}:L^{2}(\Omega)\to P_{0}(\mathcal{T}_{h})$ denotes the projection onto piecewise constants over the mesh. Let us note that, although linear finite elements are used for the velocity approximation, second order convergence in the velocity is in general not valid, which can be explained by the consistency error introduced by numerical integration in the mass lumping procedure.

The paper [28] covers triangular and tetrahedral grids, as well as quadrilateral meshes consisting of (slightly perturbed) parallelograms. In [19], the extension to hexahedral meshes was considered, which required the construction of an extension of the BDDF finite element [7]. The case of general quadrilateral or hexahedral grids, resulting from non-affine transformations of corresponding reference elements, was treated in [27]. A similar analysis was developed previously for the investigation of related finite difference and finite volume methods [2, 20]. In a recent preprint [3], high order approximations on quadrilateral and hexahedral grids have been considered.

In this paper, we propose and analyze a second order multipoint flux mixed finite element method on hybrid meshes which, in two dimensions, is based on local approximation spaces and quadrature rules as outlined in Figure 1.2.

Similar to before, the selection of appropriate basis functions allows the systematic reduction to a cell-centered discretization for the reduced problem (1.4)–(1.5). The algebraic problem (1.10) here amounts to an approximation of the pressure with piecewise linear but discontinuous finite elements, and our approach could therefore, in principle, also be interpreted as a particular discontinuous Galerkin approximation.

The accuracy of the quadrature formula depicted in Figure 1.2 seems to be insufficient at first glance. Indeed, the same arguments of proof used in [28] would lead to sub-optimal convergence orders. Our analysis, therefore, relies on special properties of the discrete spaces that allow to show the desired convergence rates of second order, i.e.

[TABLE]

where $\widetilde{p}_{h}^{\,1}$ is a piecewise linear or quadratic approximation for the pressure which can be obtained, e.g., by the post-processing strategy of Stenberg [23]. The resulting multipoint flux mixed finite element scheme thus yields the optimal approximation order in both variables. Using duality arguments, we will also establish super-convergence estimates for the pressure which here can be expressed as

[TABLE]

where $\widetilde{p}_{h}^{\,2}$ now is a piecewise quadratic approximation obtained by local post-processing.

We develop a general framework for the convergence analysis of mixed finite element methods with mass lumping which allows us to treat, in a unified manner, triangular and quadrilateral elements in two dimensions as well as tetrahedral, prismatic, and hexahedral elements in three dimensions. Our analysis also automatically covers hybrid meshes consisting of different element types; see Figure 1.3 for a sketch.

Let us note that the results presented in this paper may also be relevant for the numerical solution of wave propagation problems, where mass lumping is required to allow for an efficient realization of explicit time stepping schemes. Related mass lumping strategies have been developed for rectangular and simplicial elements in [13] and in [16, 17] for Maxwell’s equations; let us refer to [12] for an overview. In [14], we provided a full convergence analysis for a mixed finite element approximation of the acoustic wave equation with mass lumping based on the method of Wheeler and Yotov [28].

The remainder of the manuscript is organized as follows: In Section 2, we introduce our notation and some basic assumptions that guarantee the well-posedness of the problem. In Section 3, we then present a general framework for the analysis of mixed finite element approximations and show basic properties such as existence and stability. The first order estimates derived in Section 4 also cover most of the results presented in [28]. In Section 5, we then establish the improved convergence rates (1.13)–(1.14) under some general assumptions; this can be seen as the main contribution of the paper. Appropriate finite elements and quadrature rules for different element types in two and three space dimensions are discussed in Section 6 and some numerical results are presented for illustration in Section 7, including a comparison with the first order methods of [28].

2. Preliminaries

In the following, we introduce our notation and some basic assumptions that are required throughout the paper. In addition, we present a general framework for the formulation of mixed finite element approximations with mass lumping.

2.1. Problem data, function spaces, and well-posedness

Throughout the presentation, $\Omega$ denotes a bounded Lipschitz domain in two or three space dimensions. In order to ensure the well-posedness of problem (1.1)–(1.3), we assume that $\mathcal{K}\in L^{\infty}(\Omega)^{d\times d}$ and that $\mathcal{K}(x)$ is symmetric and uniformly positive definite, i.e.,

[TABLE]

for a.e. $x\in\Omega$ with some constants $0<\underline{k},\overline{k}<\infty$. For our convergence analysis, we will later impose additional conditions on $\Omega$ and $\mathcal{K}$ as required.

We denote by $L^{p}(\Omega)$ and $W^{k,p}(\Omega)$ the usual Lebesgue and Sobolev spaces, and we write $H^{k}(\Omega)=W^{k,2}(\Omega)$ for ease of notation. The scalar product and norm of $L^{2}(\Omega)$ are denoted by $(\cdot,\cdot)$ and $\|\cdot\|_{L^{2}}$. We further define $H(\mathrm{div},\Omega)=\{u\in L^{2}(\Omega)^{d}:\mathrm{div}\,u\in L^{2}(\Omega)\}$ and denote by $\|u\|_{H(\mathrm{div},\Omega)}=(\|u\|_{L^{2}(\Omega)}^{2}+\|\mathrm{div}\,u\|_{L^{2}(\Omega)}^{2})^{1/2}$ the natural norm for this space. Similar notation is used for the norm of other spaces.

Let us recall some well-known results concerning the analysis of our model problem.

Lemma 2.1**.**

For any $f\in L^{2}(\Omega)$ and $g\in H^{1/2}(\partial\Omega)$, there exists a unique solution $u\in H(\mathrm{div},\Omega)$, $p\in H^{1}(\Omega)$ of problem (1.1)–(1.3) and

[TABLE]

Proof.

Problem (1.1)–(1.3) is equivalent to the reduced problem (1.3)–(1.4), for which existence of a unique solution $p\in H^{1}(\Omega)$ follows from the Lax-Milgram theorem. Via equation (1.1), one can recover the velocity field $u=-\mathcal{K}\nabla p\in L^{2}(\Omega)^{d}$, and by inserting thisin equation (1.2), one can then deduce that $u\in H(\mathrm{div},\Omega)$. ∎

Remark 2.2**.**

Let us note that on Lipschitz domains, one can actually deduce that $u\in L^{p}(\Omega)^{d}$ for some $p>2$, which allows to define normal traces $n\cdot u$ of the velocity along sufficiently smooth $(d-1)$-dimensional manifolds; see [5] for details.

By standard arguments, one can further show the following statement.

Lemma 2.3**.**

Under the assumptions of Lemma 2.1, the solution of Problem (1.1)–(1.3) corresponds to the unique solution $u\in H(\mathrm{div},\Omega)$, $p\in L^{2}(\Omega)$ of the variational problem

[TABLE]

The key for establishing the existence of a unique solution $u\in H(\mathrm{div},\Omega)$ to the variational problem (2.2)–(2.3) via Brezzi’s splitting theorem [5, 6] is the surjectivity of the divergence operator, which we state here for later reference.

Lemma 2.4**.**

For any $q\in L^{2}(\Omega)$ there exists a function $v\in H(\mathrm{div},\Omega)$ such that $q=\mathrm{div}\,v$ and $\|v\|_{H(\mathrm{div},\Omega)}\leq C\|q\|_{L^{2}(\Omega)}$ with $C$ depending only on $\Omega$. Moreover, one can choose $v\in H(\mathrm{div},\Omega)\cap L^{p}(\Omega)^{d}$ with $p>2$ and such that

[TABLE]

The second assertion follows from regularity results for the Poisson problem on Lipschitz domains and allows to define normal traces $n\cdot v$ along smooth sub-manifolds of dimension $d-1$, which is required to prove discrete variants of Lemma 2.4; see e.g. [5] for details.

3. Mixed finite element discretization

We now introduce our notation and the general form of mixed finite element approximations that will be considered in the rest of the paper.

3.1. Mesh

Let $\mathcal{T}_{h}=\{T_{n}:n=1,\ldots,N\}$ denote a geometrically conforming partition of the domain $\Omega$ into $d$–dimensional elements in the sense of [11, 18]. In particular, neighboring elements either share a common vertex, edge, or face. Moreover, we assume that any element $T\in\mathcal{T}_{h}$ is the image $T=F_{T}(\widehat{T})$ of one of a finite number of reference elements $\widehat{T}\in\{\widehat{T}_{1},\ldots,\widehat{T}_{M}\}$ under an affine mapping

[TABLE]

A partition $\mathcal{T}_{h}$ of $\Omega$ satisfying the above conditions will simply be called a mesh in the sequel. We denote by $h_{T}$ the diameter of the element $T$ and call $h=\max_{T}h_{T}$ the global mesh size. To have some control on the shape of the elements, we further require that the mesh $\mathcal{T}_{h}$ is uniformly shape regular, i.e., there exists $c_{\gamma}>0$ such that

[TABLE]

where $h_{T}$ is the diameter of the element $T$. Without loss of generality, one may assume that $c_{\gamma}$ is chosen such that $c_{\gamma}^{-1}h_{T}^{d}\leq|\det B_{T}|\leq c_{\gamma}h_{T}^{d}$ holds as well.

Using standard convention, we denote by $W^{k,p}(\mathcal{T}_{h})=\{v\in L^{2}(\Omega):v|_{T}\in W^{k,p}(T)\}$ the broken Sobolev spaces over the mesh $\mathcal{T}_{h}$ and we use $H^{k}(\mathcal{T}_{h})=W^{k,2}(\mathcal{T}_{h})$ for abbreviation.

3.2. Approximation spaces

Let $P_{k}(\mathcal{T}_{h})=\{v\in L^{2}(\Omega):v|_{T}\in P_{k}(T)\}$ be the space of piecewise polynomials over the mesh $\mathcal{T}_{h}$ of degree smaller or equal to $k$. We assume that the approximation spaces $V_{h}$, $Q_{h}$ in the discrete variational problem (1.6)–(1.7) satisfy

• (A1)

$V_{h}\subset P_{k}(\mathcal{T}_{h})^{d}\cap H(\mathrm{div},\Omega)$ and $Q_{h}=\mathrm{div}\,V_{h}\supset P_{0}(\mathcal{T}_{h})^{d}$.

We denote by $\pi_{h}:L^{2}(\Omega)\to Q_{h}$ the $L^{2}$-orthogonal projection onto $Q_{h}$, defined by

[TABLE]

and use $\pi_{h}^{k}:L^{2}(\Omega)\to P_{k}(\mathcal{T}_{h})$ for the orthogonal projections onto $P_{k}(\mathcal{T}_{h})$ and $P_{k}(\mathcal{T}_{h})^{d}$. We further assume to have access to a suitable interpolation operator $\Pi_{h}$ for the velocity space with the following properties.

• (A2)

$\Pi_{h}:H(\mathrm{div},\Omega)\cap L^{p}(\Omega)^{d}\to V_{h}$, $p>2$ is a linear and bounded projection operator satisfying the commuting diagram property

[TABLE]

These assumptions will be used to establish the well-posedness of the discrete variational problem (1.6)–(1.7). In particular, they yield the following discrete variant of Lemma 2.4.

Lemma 3.1**.**

Let (A1)–(A2) hold. Then for any $q_{h}\in Q_{h}$ there exists a function $v_{h}\in V_{h}$ such that $q_{h}=\mathrm{div}\,v_{h}$ and $\|v_{h}\|_{H(\mathrm{div})}\leq C\|q_{h}\|_{L^{2}(\Omega)}$.

Proof.

By Lemma 2.4, there exists a function $v\in H(\mathrm{div},\Omega)\cap L^{p}(\Omega)^{d}$ such that $\mathrm{div}v=q_{h}$ and $\|\mathrm{div}v\|+\|v\|_{L^{p}}\leq C\|q_{h}\|$. Then the choice $v_{h}=\Pi_{h}v$ together with the properties in assumption (A2) immediately yield the result. ∎

A variety of particular approximation spaces $V_{h}$, $Q_{h}$ and corresponding interpolation and projection operators $\Pi_{h}$, $\pi_{h}$ that satisfy assumptions (A1)–(A2) can be found in the literature [7, 8, 9, 22]; we refer to [5] and to Section 6 for details and examples.

3.3. Quadrature rule

As a final ingredient for the formulation of our method, we need to specify the type of quadrature rule to be used in the definition of the discrete variational problem (1.6)–(1.7). For any $u_{h},v_{h}\in V_{h}$, let us define

[TABLE]

where $(x_{n,T},w_{n,T})$, $n=1,\ldots,N_{T}$ denotes the integration points and weights specifying the local quadrature rules on the element $T$. Note that the definition also naturally applies to scalar valued functions. In order to ensure the well-posedness of the discrete variational problem (1.6)–(1.7) defined below, we require that

• (A3)

all weights $\hat{w}_{n,T}$ are positive and the lumped scalar product $(\cdot,\cdot)_{h}$ is uniformly equivalent to the exact scalar product $(\cdot,\cdot)$ on the discrete velocity space $V_{h}$, i.e.,

[TABLE]

for all $v_{h}\in V_{h}$ with some generic constant $c>0$ independent of $V_{h}$.

If $\mathcal{K}$ is bounded from above and below, and piecewise continuous, one can even show that

[TABLE]

for all discrete functions $v_{h}\in V_{h}$, which will be utilized below. In the sequel, we always assume that $\mathcal{K}$ is piecewise continuous.

3.4. Discrete variational problem

We next show that the system (1.6)-(1.7) attains a unique solution. For later reference, we here consider a slightly more general problem, namely: Find $(w_{h},r_{h})\in V_{h}\times Q_{h}$ such that

[TABLE]

where $a_{h}:V_{h}\to\mathbb{R}$ and $b_{h}:Q_{h}\to\mathbb{R}$ are some given linear functionals. Standard arguments for the analysis of mixed variational problems yield the following result.

Lemma 3.2**.**

Let (A1)–(A3) hold and let $a_{h}:V_{h}\to\mathbb{R}$ and $b_{h}:Q_{h}\to\mathbb{R}$ be bounded linear functionals. Then problem (3.8)–(3.9) has a unique solution $w_{h}\in V_{h}$, $r_{h}\in Q_{h}$ with

[TABLE]

with a constant $C$ that only depends on the constants in the assumptions.

Proof.

It suffices to establish the conditions of Brezzi’s splitting theorem [6].

(i) One has $(\mathrm{div}\,v_{h},q_{h})\leq\|v_{h}\|_{H(\mathrm{div})}\|q_{h}\|_{L^{2}(\Omega)}$ and

[TABLE]

i.e., the bilinear forms of the discrete variational problem (3.8)-(3.9) are uniformly continuous on the spaces $V_{h}\times Q_{h}$ and $V_{h}\times V_{h}$, respectively.

(ii) Lemma 3.1 implies the discrete inf-sup condition

[TABLE]

with $\beta>0$ depending only on the constant from Lemma 3.1.

(iii) By assumption (A1), we know that $\mathrm{div}\,V_{h}=Q_{h}$. The condition $(\mathrm{div}\,v_{h},q_{h})=0$ for all functions $q_{h}$ therefore implies that $\mathrm{div}\,v_{h}=0$. From (A3), we know that

[TABLE]

Together with the above considerations, we obtain

[TABLE]

for all $v_{h}\in V_{h}$ satisfying $(\mathrm{div}\,v_{h},q_{h})=0$ for all $q_{h}\in Q_{h}$ and with $\alpha=\underline{k}c^{-1}>0$. This is the discrete kernel ellipticity condition required for Brezzi’s splitting theorem. The assertions of the Lemma are then a direct consequence of the results in [6]; see [5] for further details. ∎

From the previous Lemma one can deduce that for any $f\in L^{2}(\Omega)$ and $g\in H^{1/2}(\partial\Omega)$, the discrete variational problem (1.6)–(1.7) admits a unique solution $(u_{h},p_{h})\in V_{h}\times Q_{h}$.

4. First order estimates

Based on assumptions (A1)–(A3), we now present an abstract error analysis for the discrete variational problem (1.6)–(1.7). We start with splitting the error

[TABLE]

into approximation errors and discrete error components. Some additional assumptions will be required to bound the approximation errors and the consistency errors introduced by the numerical quadrature. In this section, we establish convergence rates of first order; the improved estimates (1.13)–(1.14) announced in the introduction will be proved in the next section.

4.1. Approximation error

In order to ensure sufficient approximation properties of the spaces $V_{h}$ and $Q_{h}$, we assume that

• (A4)

the projection operators $\Pi_{h}$ and $\pi_{h}$ are defined locally, i.e., $(\Pi_{h}v)|_{T}=\Pi_{T}v|_{T}$ and $(\pi_{h}q)|_{T}=\pi_{T}q|_{T}$ for appropriate operators $\Pi_{T}$ and $\pi_{T}$ which satisfy

[TABLE]

for all elements $T\in\mathcal{T}_{h}$ with a uniform constant $C$.

This condition already allows to bound the approximation errors in the above error splitting. In addition, we can show that the projection $\Pi_{h}$ is actually stable in $H^{1}(\Omega)$.

Lemma 4.1**.**

Let (A1)–(A4) hold. Then, for all $u\in H^{1}(T)$ there holds

[TABLE]

with a constant $C$ that is independent of the element $T$ and of the function $u$.

Proof.

By the triangle inequality and Assumption (A4), we obtain

[TABLE]

Now let $\pi_{T}^{0}:L^{2}(T)^{d}\to\text{P}_{0}(T)^{d}$ be the $L^{2}$-orthogonal projection. Then

[TABLE]

where we used the triangle inequality and the approximation properties of the two projection operators. Combination of the two estimates yields the assertion. ∎

4.2. Quadrature error

As a second ingredient for our convergence analysis, we require that the numerical quadrature is sufficiently accurate. In particular, we assume that

• (A5)

the quadrature rule is exact on $P_{0}(\mathcal{T}_{h})^{d}\times V_{h}$, i.e.,

[TABLE]

As a direct consequence of (A5), we then obtain the following estimates.

Lemma 4.2**.**

Let assumptions (A1)–(A5) hold and further assume that $\mathcal{K}\in W^{1,\infty}(\mathcal{T}_{h})$. Then for any $u\in H(\mathrm{div},\Omega)\cap H^{1}(\mathcal{T}_{h})^{d}$ and any $v_{h}\in V_{h}$, there holds

[TABLE]

Proof.

Using Assumption (A5) and (3.7), one can see that

[TABLE]

Using the product rule of differentiation, we can further estimate

[TABLE]

Using $h_{T}^{d/2}\|w_{h}\|_{L^{\infty}(T)}\leq c\|w_{h}\|_{L^{2}(T)}$ for $w_{h}=v_{h},\Pi_{h}u$, which follows from a scaling argument, the regularity of $\mathcal{K}^{-1}$, and summing over all elements then yields the result. ∎

4.3. First order estimates

In the sequel, let $(u,p)$ be a sufficiently smooth solution of problem (1.1)–(1.3) and $(u_{h},p_{h})$ solve the discrete variational problem (1.6)–(1.7). Then, by combination of the previous results, we already obtain the following assertion.

Theorem 4.3**.**

Let (A1)–(A5) hold and $\mathcal{K}\in W^{1,\infty}(\mathcal{T}_{h})$. Then

[TABLE]

If, additionally, $\mathrm{div}\,u\in H^{1}(\mathcal{T}_{h})$, then also $\|\mathrm{div}(u-u_{h})\|_{L^{2}(\Omega)}\leq Ch\|\mathrm{div}\,u\|_{H^{1}(\mathcal{T}_{h})}$.

Proof.

In view of the error splitting (4.1)-(4.2) and assumption (A4), it remains to estimate the discrete error components $w_{h}=\Pi_{h}u-u_{h}$ and $r_{h}=\pi_{h}p-p_{h}$. Using the discrete and continuous variational principles (1.6)–(1.7) and (2.2)–(2.3), one can see that the discrete error components $(w_{h},r_{h})$ satisfy the equations (3.8)–(3.9) with

[TABLE]

By assumption (A4) and Lemma 4.2, we immediately obtain

[TABLE]

The first assertion then follows readily from Lemma 3.2, and the estimate of the divergence error is a direct consequence of condition (A2) and the properties of the projection operators $\Pi_{h}$ and $\pi_{h}$ stated in assumption (A4). ∎

Remark 4.4**.**

Let us note that the above theorem also covers the first order estimates derived in [28] for their choice of spaces and quadrature rules; cf Figure 1.1.

5. Improved estimates

We now establish the improved estimates (1.13)–(1.14) that were announced in the introduction and which will be summarized in a theorem at the end of this section.

5.1. Estimates for the velocity

For our analysis, we require some additional assumptions that will be introduced as needed. First of all, we assume that the space $V_{h}$ has sufficient approximation properties, i.e.,

• (A6)

the local projection operators $\Pi_{T}$ introduced in (A4) satisfy

[TABLE]

with a uniform constant $C$ for all elements $T\in\mathcal{T}_{h}$.

In addition, we also require higher exactness of the quadrature rule together with an additional special assumption on the spaces. We thus assume that

• (A7)

the quadrature rule is exact on $P_{1}(\mathcal{T}_{h})^{d}\times V_{h}^{*}$, i.e.,

[TABLE]

for all functions $w_{h}\in P_{1}(\mathcal{T}_{h})^{d}$ and $v_{h}^{*}\in V_{h}^{*}$. Moreover, there exists $V_{h}^{*}\subset V_{h}$ such that $v_{h}\in V_{h}$ and $\mathrm{div}\,v_{h}\in Q_{h}^{*}:=\mathrm{div}\,V_{h}^{*}$ imply $v_{h}\in V_{h}^{*}$.

As a direct consequence, we thenobtain the following bounds for the quadrature error.

Lemma 5.1**.**

Let (A1)–(A7) hold and further assume that $\mathcal{K}\in W^{2,\infty}(\mathcal{T}_{h})$. Then for any $u\in H(\mathrm{div},\Omega)\cap H^{1}(\mathcal{T}_{h})^{d}$ and any $v_{h}^{*}\in V_{h}^{*}$, there holds

[TABLE]

Proof.

Since the quadrature rule is exact for functions in $\text{P}_{1}(\mathcal{T}_{h})^{d}\times V_{h}^{*}$, we can write

[TABLE]

The result then follows with the same arguments as used in the proof of Lemma 4.2. ∎

We now already obtain the following improved error estimates for the velocity.

Lemma 5.2**.**

Let (A1)–(A7) hold. Moreover, assume that $\mathcal{K}\in W^{2,\infty}(\mathcal{T}_{h})$. Then

[TABLE]

Moreover, if $\mathrm{div}\,u\in H^{2}(\mathcal{T}_{h})$, we have $\|\mathrm{div}(u-u_{h})\|_{L^{2}(\Omega)}\leq Ch^{2}\|\mathrm{div}\,u\|_{H^{2}(\mathcal{T}_{h})}$.

Proof.

We proceed with similar arguments as in the proof of Theorem 4.3. Denote by $w_{h}=\Pi_{h}u-u_{h}$ and $r_{h}=\pi_{h}p-p_{h}$ the discrete error components. Using assumption (A2), we deduce from equations (2.3) and (1.7) that

[TABLE]

for all $q_{h}\in Q_{h}$. Since $Q_{h}=\mathrm{div}\,V_{h}$, this implies $\mathrm{div}\,\Pi_{h}u\equiv\mathrm{div}\,u_{h}$ and therefore

[TABLE]

From the variational characterizations (2.2)–(2.3) and (1.6)–(1.7) of the continuous and the discrete solution, we can thus infer that

[TABLE]

From equation (5.6) and assumption (A3), we can further deduce that $w_{h}\in V_{h}^{*}$. Together with assumption (A6) and Lemma 5.1, we then obtain

[TABLE]

This immediately implies that

[TABLE]

which is the desired estimate for the velocity error in $L^{2}$. From equation (5.5) and assumption (A2), we further deduce that

[TABLE]

which yields the remaining estimate for the divergence error. ∎

Remark 5.3**.**

Note that, due to assumption (A1), the velocity error only depends on the approximation properties of the velocity space $V_{h}$; see [5]. This will be the case also for the estimates of the pressure error presented in the following sections.

5.2. Estimates for the pressure

In order to prove second order estimates for the pressure, we additionally require that the spaces $V_{h}^{*}$ and $Q_{h}^{*}$ make up a stable pairing. We denote by $\pi_{h}^{*}:L^{2}(\Omega)\to Q_{h}^{*}$ the $L^{2}$-orthogonal projection onto $Q_{h}^{*}$, defined by

[TABLE]

and we now assume that there exists

• (A8)

a linear bounded projection operator $\Pi_{h}^{*}:H^{1}(\mathcal{T}_{h})^{d}\cap H(\mathrm{div},\Omega)\to V_{h}^{*}$ such that the commuting property $\pi_{h}^{*}\mathrm{div}\,v=\mathrm{div}\,\Pi_{h}^{*}v$ holds for all $v\in H(\mathrm{div},\Omega)\cap H^{1}(\mathcal{T}_{h})^{d}$.

Together with the second order estimates for the velocity error, we now obtain.

Lemma 5.4**.**

Let (A1)–(A8) hold. Moreover, assume that $\mathcal{K}\in W^{2,\infty}(\mathcal{T}_{h})$. Then

[TABLE]

Proof.

By assumption (A7), we have $\pi_{h}^{*}(p-p_{h})=\mathrm{div}\,v_{h}^{*}$ for some $v_{h}^{*}\in V_{h}^{*}$ with bound $\|v_{h}^{*}\|_{H(\mathrm{div})}\leq\beta^{-1}\|\pi_{h}^{*}(p-p_{h})\|_{L^{2}(\Omega)}$. Using the variational equations (2.2), (1.6), as well as assumption (A8), we further deduce that

[TABLE]

In the last step, we employed the equivalence of norms (A3), the bounds for the coefficients, the estimates for the approximation error in Lemma 5.2, the improved estimate for the quadrature error in (A8), and the previous estimate for the discrete velocity error. The estimate for $\|\pi_{h}^{*}(p-p_{h})\|_{L^{2}(\Omega)}$ is now obtained by the stability estimate $\|v_{h}^{*}\|_{L^{2}(\Omega)}\leq\|v_{h}^{*}\|_{H(\mathrm{div})}\leq\beta^{-1}\|\pi_{h}^{*}(p-p_{h})\|_{L^{2}(\Omega)}$ in assumption (A7). Since $P_{0}(\mathcal{T}_{h})\subset Q_{h}^{*}$, we further have $\|\pi_{h}^{0}(p-p_{h})\|_{L^{2}(\Omega)}\leq\|\pi_{h}^{*}(p-p_{h})\|_{L^{2}(\Omega)}$, which completes the proof. ∎

5.3. Super-convegence for the pressure

As a last step, we now show that even third order convergence for the pressure can be obtained, assuming the domain $\Omega$ is convex.

Lemma 5.5**.**

Let (A1)–(A8) hold, $\Omega$ be convex, and $\mathcal{K}\in W^{1,\infty}(\Omega)\cap W^{2,\infty}(\mathcal{T}_{h})$. Then

[TABLE]

Proof.

Let $\phi$ denote the solution to the auxiliary problem

[TABLE]

and note that $\|\phi\|_{H^{2}(\Omega)}\leq C\|\pi_{h}^{*}(p-p_{h})\|_{L^{2}(\Omega)}$ due to convexity of the domain. Then, by assumption (A8) for the projectors $\pi_{h}^{*}$ and $\Pi_{h}^{*}$, we obtain

[TABLE]

Using that $(u-u_{h},\nabla\phi)=-(\mathrm{div}(u-u_{h}),\phi)=-(\mathrm{div}(u-u_{h}),\phi-\pi_{h}^{1}\phi)$, which follows from (2.3) and (1.7), the first term can be further estimated by

[TABLE]

For the second term, we can use assumption (A7), similar arguments as in the proof of Lemma 5.1, and the regularity of $\mathcal{K}$, to show that locally

[TABLE]

With similar arguments as in Lemma 4.1, one can show that $\|\Pi_{T}u\|_{H^{2}(T)}\leq C\|u\|_{H^{2}(T)}$ and $\|\Pi_{h}^{*}(\mathcal{K}\nabla\phi)\|_{H^{1}(T)}\leq C\|\mathcal{K}\nabla\phi\|_{H^{1}(T)}$. From Lemma 5.4 and assumption (A6), we further obtain $\|\Pi_{T}u-u_{h}\|_{L^{2}(\Omega)}\leq Ch_{T}^{2}\|u\|_{H^{2}(\mathcal{T}_{h})}$. By summation over all elements, we then get

[TABLE]

and since $P_{0}(\mathcal{T}_{h})\subset Q_{h}^{*}$, we also have $\|\pi_{h}^{0}(p-p_{h})\|_{L^{2}(\Omega)}\leq\|\pi_{h}^{*}(p-p_{h})\|_{L^{2}(\Omega)}$. ∎

5.4. Post-processing for the pressure

In the spirit of Stenberg [23], we now define the following local post-processing procedure.

Problem 5.6** (Post-processing for the pressure).**

Find $\widetilde{p}_{h}\in\text{P}_{2}(\mathcal{T}_{h})$ such that for all $T\in\mathcal{T}_{h}$ there holds

[TABLE]

Note that $\widetilde{p}_{h}$ can be computed separately on each element $T$, rendering the method computatinally efficient. Following the analysis of [23], we obtain the following bounds.

Lemma 5.7**.**

Let (A1)–(A8) hold. Then

[TABLE]

Proof.

For convenience of the reader, we recall the basic steps of the derivation of this result. By the triangle inequality, one obtains

[TABLE]

The first term already appears in the final estimate. For the second, observe that

[TABLE]

where we used the properties of the $L^{2}$ projections and (5.10). By the optimality of the $L^{2}$–projection and the Poincaré inequality, the third term can be estimated by

[TABLE]

Using equation (2.2), we obtain for $\widetilde{q}_{h}=\pi_{h}^{2}p-\widetilde{p}_{h}$ that

[TABLE]

Dividing by $\|\nabla\widetilde{q}_{h}\|_{L^{2}(T)}$ yields the estimate for the term (iii). The assertion of the lemma now follows by combination with the previous estimates. ∎

By application of the previous estimates, we immediately obtain the following assertion.

Lemma 5.8**.**

Let (A1)–(A8) hold and $\mathcal{K}\in W^{2,\infty}(\mathcal{T}_{h})$. Then

[TABLE]

If, in addition, $\Omega$ is convex and $\mathcal{K}\in W^{1,\infty}(\Omega)$, then

[TABLE]

Proof.

By the approximation properties of assumption (A6) and Lemma 5.5, we can estimated the terms on the right hand side in Lemma 5.7, which yields the result. ∎

Remark 5.9**.**

In a similar way as in Problem 5.6, one could construct a piecewise linear approximation $\widetilde{p}_{h}^{\;1}$, for which the first estimate of the previous lemma still remains valid. Let us note again that all estimates of this section, in particular also those for the pressure error, only depend on the approximation properties of the velocity space $V_{h}$.

5.5. Summary of improved error estimates

For later reference, let us summarize the results of the analysis given in this section in a single theorem.

Theorem 5.10**.**

Let (A1)–(A8) hold and $\mathcal{K}\in W^{2,\infty}(\mathcal{T}_{h})$. Then

[TABLE]

If, in addition, $\Omega$ is convex and $\mathcal{K}\in W^{1,\infty}(\Omega)$, then

[TABLE]

In the following section, we discuss appropriate finite elements and quadrature rules, for which our abstract convergence results apply.

6. Finite elements

Let us briefly comment on the general setting which is used in two and three space dimensions: We assume that $\mathcal{T}_{h}=\{T\}$ is a shape regular conforming mesh and that every element $T\in\mathcal{T}_{h}$ is the affine image of a reference triangle or square, or a reference tetrahedron, cube, or prisms, i.e., $F_{T}(\widehat{T})$ with $F_{T}(\widehat{x})=a_{T}+B_{T}\widehat{x}$. The subspaces used in our approximations are defined locally by

[TABLE]

with local spaces $V_{T}=\{v=\frac{1}{\det B_{T}}B_{T}\widehat{v}\circ F_{T}^{-1}:\widehat{v}\in\widehat{V}\}$ and $Q_{T}=\{q=\widehat{q}\circ F_{T}^{-1}:\widehat{q}\in\widehat{Q}\}$. Additionally, spaces $V_{h}^{*}$ and $Q_{h}^{*}$ are defined accordingly via mapping from the local spaces $\widehat{V}^{*}$ and $\widehat{Q}^{*}$. Only the spaces $\widehat{V}$, $\widehat{Q}$, $\widehat{V}^{*}$ and $\widehat{Q}^{*}$ on the reference elements will therefore be defined below. The local quadrature rule in all cases is of the form

[TABLE]

where $\widehat{r}_{l}$ and $\widehat{\omega}_{l}$, $l=1,\ldots,N_{T}$ are the quadrature points and weights on the reference element, and $|\widehat{T}|$ is the area of the element $\widehat{T}$; compare with equation (3.5).

In the following, we introduce ansatz spaces for different reference elements and appropriate quadrature rules for which the conditions (A1)–(A8) required for our theoretical results are valid. In addition, we discuss a choice of basis for the spaces $\widehat{V}$ such that the resulting mass matrix for the velocity variable becomes block diagonal; see Figure 1.3 and the discussion in the introduction. We start with two dimensional elements.

6.1. Triangles

As quadrature points for the triangle, we use its vertices $\widehat{r}_{1}=(0,0)$, $\widehat{r}_{2}=(1,0)$, $\widehat{r}_{3}=(0,1)$ and its barycenter $\widehat{r}_{4}=(\frac{1}{3},\frac{1}{3})$ together with weights $\omega_{1}=\omega_{2}=\omega_{3}=\frac{1}{12}$ and $\omega_{4}=\frac{3}{4}$. Note that the quadrature rule is exact for functions in $\text{P}_{2}(\widehat{T})$.

The local function spaces on the reference triangle are defined as

[TABLE]

and we denote by $\Pi_{T}$ and $\Pi_{T}^{*}$ the standard interpolation operators for these spaces; for details see [5]. For this choice, all our assumptions (A1)–(A8) are valid.

As a basis for the reference element $\widehat{V}=\text{RT}_{1}(\widehat{T})$, which complies with the quadrature points and thus leads to a block diagonal mass matrix, we choose

[TABLE]

Note that, as desired, exactly two basis functions are associated to every quadrature point while all other basis functions vanish at this point; see Figure 1.2.

6.2. Quadrilaterals

The quadrature rule on the reference rectangle is again given by its vertices $\widehat{r}_{1}=(0,0)$, $\widehat{r}_{2}=(1,0)$, $\widehat{r}_{3}=(1,0)$, and $\widehat{r}_{4}=(0,1)$ and the barycenter $\widehat{r}_{5}=(\frac{1}{2},\frac{1}{2})$ with weights $\omega_{i}=\frac{1}{12}$ for $1\leq i\leq 4$ and $\omega_{5}=\frac{2}{3}$. One can verify that this quadrature formula is exact even for functions in $\text{P}_{3}(\widehat{T})$, which seems a bit surprising.

As function spaces on the reference element, we now use

[TABLE]

and we denote by $\Pi_{T}=\Pi_{T}^{*}$ the standard interpolation operator for the $\text{BDFM}_{2}$ space; for details, see [5]. All conditions (A1)–(A8) required for our theory are valid for this choice.

As basis for $\text{BDFM}_{2}(\widehat{T})$ that complies with the above quadrature rule, we choose

[TABLE]

Note that again exactly two basis functions can be associated to every quadrature point while all other basis functions vanish at this point. As a consequence, we obtain a block diagonal mass matrix for the velocity variable; see Figure 1.3 in the introduction.

As a next step, we define appropriate finite elements in three space dimensions.

6.3. Tetrahedra

As quadrature points, we again choose the vertices $\widehat{r}_{1}=(0,0,0)$, $\widehat{r}_{2}=(1,0,0)$, $\widehat{r}_{3}=(0,1,0)$, $\widehat{r}_{4}=(0,0,1)$ of the reference tetrahedron and its barycenter $\widehat{r}_{4}=(\frac{1}{4},\frac{1}{4},\frac{1}{4})$ and we set the weights to $\omega_{i}=\frac{1}{20}$ for $1\leq i\leq 4$ and $\omega_{5}=\frac{4}{5}$. Let us note that this quadrature formula is exact for functions in $\text{P}_{2}(\widehat{T})$.

As function spaces on the reference element, we here choose

[TABLE]

and we denote by $\Pi_{T}$ and $\Pi_{T}^{*}$ the standard interpolation operators for these spaces; see again [5] for details. All conditions (A1)–(A8) are then valid for this choice.

We next construct an appropriate basis for the space $\widehat{V}=\text{RTN}_{1}(\widehat{T})$ on the reference element. We start by defining three functions

[TABLE]

associated to the interior quadrature point, and then define $12$ additional basis functions for the vertices by

[TABLE]

Let us note that exactly three basis functions are associated to every quadrature point while all other basis functions vanish at this point; see Figure 6.1. As a consequence, the resulting mass matrix for the velocity variable will be block diagonal.

6.4. Hexahedra

The quadrature rule for the reference hexahedron is again defined by its vertices $\widehat{r}_{1}=(0,0,0)$, $\widehat{r}_{2}=(1,0,0)$, $\widehat{r}_{3}=(0,1,0)$, $\widehat{r}_{4}=(1,1,0)$, $\widehat{r}_{5}=(0,0,1)$, $\widehat{r}_{6}=(1,0,1)$, $\widehat{r}_{7}=(0,1,1)$, $\widehat{r}_{8}=(1,1,1)$ and its barycenter $\widehat{r}_{4}=(\frac{1}{2},\frac{1}{2},\frac{1}{2})$ with weights $\omega_{i}=\frac{1}{24}$ for $1\leq i\leq 8$ and $\omega_{9}=\frac{2}{3}$. One can verify that this quadrature formula is exact for functions in $\text{P}_{3}(\widehat{T})\oplus\text{span}\{x^{2}yz,xy^{2}z,xyz^{3}\}$ including ceratin polynomials of fourth order.

Now let $\mathcal{Q}_{1}(\widehat{T})$ denote the space of trilinear functions on the reference hexahedron $\widehat{T}$. As local ansatz spaces for the hexahedron, we then choose

[TABLE]

and we set $\widehat{Q}=\widehat{Q}^{*}=\mathrm{div}\,\widehat{V}$; note that the function $xyz$ is not contained in $\widehat{Q}$. In the usual manner [5], an interpolation operator $\Pi_{T}=\Pi_{T}^{*}$ for $\widehat{V}$ can be defined by

[TABLE]

This operator has the required approximation properties required for our theory and also satisfies the commuting diagram property with the $L^{2}$–projection on the space $\widehat{Q}=\widehat{Q}^{*}$. In summary, the conditions (A1)–(A8) required for our error estimates are valid.

In order to simplify the defintion of a basis for $\widehat{V}$, we introduce two auxiliary functions

[TABLE]

Then the basis for $\widehat{V}$ for the reference element $\widehat{T}$ is defined by the $24$ functions

[TABLE]

which are associated to the $8$ vertices, as well as the $3$ basis functions

[TABLE]

for the barycenter. Note that again exactly three basis functions are associated to every quadrature point while all other basis functions vanish at this point; see Figure 6.1. As a consequence, the mass matrix for the velocity variable will be block diagonal.

6.5. Prisms

The quadrature points for the reference prism are defined by its vertices $\widehat{r}_{1}=(0,0,0)$, $\widehat{r}_{2}=(1,0,0)$, $\widehat{r}_{3}=(0,1,0)$, $\widehat{r}_{4}=(0,0,1)$, $\widehat{r}_{5}=(1,0,1)$, $\widehat{r}_{6}=(0,1,1)$ and two additional interior points $\widehat{r}_{8}=(\frac{1}{3},\frac{1}{3},\frac{1}{3})$ and $\widehat{r}_{7}=(\frac{1}{3},\frac{1}{3},\frac{2}{3})$. Let us note that, together with the weights $\omega_{i}=\frac{1}{24}$ for $1\leq i\leq 6$ and $\omega_{7}=\omega_{8}=\frac{3}{8}$, this quadrature formula is exact for functions in $\text{P}_{1}(\widehat{T})\oplus\text{span}\{1,x,y,z,xz,yz\}$.

The auxiliary space $\widehat{V}^{*}$ here is spanned by the following set of basis functions

[TABLE]

and we set $\widehat{Q}^{*}=\mathrm{div}\,\widehat{V}^{*}=\text{P}_{1}(\widehat{T})$ as before. In order to define the space $\widehat{V}$, we enrich the space $\widehat{V}^{*}$ by additional $6$ basis functions

[TABLE]

associated to the two inner quadrature points. In the usual manner [3, 5], one can define appropriate interpolation operators $\Pi_{T}$, $\Pi_{T}^{*}$ for the spaces $\widehat{V}$, $\widehat{V}^{*}$ such that all conditions (A1)–(A8) required for our theory are valid.

As a basis for $\widehat{V}$ that complies with the quadrature rule given above, we choose

[TABLE]

for the vertices and $\Phi_{19},\ldots,\Phi_{24}$ for the two inner quadrature points. Note that again exactly three basis functions are associated to each quadrature point while all other basis functions vanish at this point; see Figure 6.2. As a consequence, we again obtain a block diagonal mass matrix for the velocity variable.

Remark 6.1**.**

For the above choices of reference elements, all assumptions (A1)–(A8) required for our theory are valid. A quick inspection shows that the degrees of freedom of different elements are compatible, such that we can consider hybrid meshes without further considerations. By construction of the basis functions, the mass matrices for the velocity variable will be block diagonal. In summary, we thus obtain a second order multipoint flux approximation in the spirit of Wheeler and Yotov [28].

7. Numerical illustration

We now illustrate our theoretical results by a simple numerical example.

7.1. Test problem

As computational domain for our tests, we choose $\Omega=\Omega_{1}\cup\Omega_{2}$, where $\Omega_{1}=(-1,1)\times(-1,0)$ and $\Omega_{2}=\{(x,y):x^{2}+y^{2}<1\text{ and }y\geq 0\}$. The exact pressure and conductivity matrix are defined as in [28] by

[TABLE]

and we set $u=-\mathcal{K}\nabla p$. These functions fulfill all regularity assumptions made throughout this manuscript, so we expect to observe the full convergence rates. We then partition $\Omega_{1}$ into triangles and $\Omega_{2}$ into rectangles in a conforming way, i.e. such that no hanging nodes are present at the interface.

7.2. First order elements

For approximation of the solution, we first consider the spaces $\text{BDM}_{1}$–$\text{P}_{0}$ used in [28]; compare with Figure 1.1 in the introduction and the correpsonding discussion. For this choice, the assumptions (A1)–(A5) are valid and the theoretical results of Section 4 apply. In Table 7.1, we display the corresponding results obtained in our computational tests on a sequence of non-nested meshes.

As predicted by our theory, we observe first order convergence for the velocity and second order convergence for the projected pressure.

7.3. Second order approximation

For the second test, we consider the approximation in $\text{RT}_{1}$–$\text{P}_{1}$ for triangles and $\text{BDFM}_{2}$–$\text{P}_{1}$ for quadrilaterals. Following the results of Section 5, we here expect second order convergence for the velocity and pressure, and third order for the projected pressure. The results of our computations are depicted in Table 7.2.

Again, the predicted convergence rates are also observed in our numerical tests. A graphical presentation of the solution is given in Figure 7.1.

8. Discussion

In this paper, we extended the multipoint flux mixed finite element method of Wheeler and Yotov [28] to second order elements and hybrid meshes in two and three space dimensions. The conditions for our abstract convergence theory were verified by construction of suitable finite elements and quadrature rules for various element types. Appropriate basis functions were defined that lead to block-diagonal mass matrices for the $H(\mathrm{div})$ variable.

Before closing, let us mention some directions for further extension and application of our results: The mass-lumping strategy proposed in this paper can also be used for the efficient numerical approximation of acoustic wave propagation. We refer to [12, 14] for results in this direction. In a similar manner as presented here, also the construction of mass-lumping strategies for $H(\mathrm{curl})$ finite elements is possible; see  [15]. Moreover, the main arguments in our analysis can in principle also be generalized to higher order. The construction of appropriate quadrature formulas and corresponding basis functions, of course, becomes more tedious when further increasing the polynomial degree.

Acknowledgements

The authors are grateful for financial support by the “Excellence Initiative” of the German Federal and State Governments via the Graduate School of Computational Engineering GSC 233 at Technische Universität Darmstadt and by the German Research Foundation (DFG) via grants TRR 146, TRR 154, and Eg-331/1-1.