Mixed Finite Element Methods for the Poisson Equation Using Biorthogonal and Quasi-Biorthogonal Systems

Abstract

We introduce two three-field mixed formulations for the Poisson equation and propose finite element methods for their approximation. Both mixed formulations are obtained by introducing a weak equation for the gradient of the solution by means of a Lagrange multiplier space. Two efficient numerical schemes are proposed based on using a pair of bases for the gradient of the solution and the Lagrange multiplier space forming biorthogonal and quasi-biorthogonal systems, respectively. We also establish an optimal a priori error estimate for both finite element approximations.

1. Introduction

In many practical situations, it is important to compute dual variables of partial differential equations more accurately. For example, the gradient of the solution is the dual variable in case of the Poisson equation, whereas the stress or pressure variable is the dual variable in case of elasticity equation. Working with the standard finite element approach these variables should be obtained a posteriori by differentiation, which will result in a loss of accuracy. In these situations, a mixed method is often preferred as these variables can be directly computed using a mixed method.

In this paper, we introduce two mixed finite element methods for the Poisson equation using biorthogonal or quasi-biorthogonal systems. Both formulations are obtained by introducing the gradient of the solution of Poisson equation as a new unknown and writing an additional variational equation in terms of a Lagrange multiplier. This gives rise to two additional vector unknowns: the gradient of the solution and the Lagrange multiplier. In order to obtain an efficient numerical scheme, we carefully choose a pair of bases for the space of the gradient of the solution and the Lagrange multiplier space in the discrete setting. Choosing the pair of bases forming a biorthogonal or quasi-biorthogonal system for these two spaces, we can eliminate the degrees of freedom associated with the gradient of the solution and the Lagrange multiplier and arrive at a positive definite formulation. The positive definite formulation involves only the degrees of freedom associated with the solution of the Poisson equation. Hence a reduced system is obtained, which is easy to solve. The first formulation is discretized by using a quasi-biorthogonal system, whereas the second one, which is a stabilized version of the first one, is discretized using a biorthogonal system.

There are many mixed finite element methods for the Poisson equation [1–8]. However, all of them are based on the two-field formulation of the Poisson equation and hence are not amenable to the application of the biorthogonal and quasi-biorthogonal systems. We need a three-field formulation to apply the biorthogonal and quasi-biorthogonal systems which leads to a symmetric formulation (see [9] for a three-field formulation in linear elasticity). The use of biorthogonal and quasi-biorthogonal systems allows an easy static condensation of the auxiliary variables (gradient of the solution and Lagrange multiplier) leading to a reduced linear system. These variables can be recovered just by inverting a diagonal matrix. Therefore, in this paper we present three-field formulations of the Poisson equation to apply the biorthogonal and quasi-biorthogonal systems.

The structure of the rest of the paper is as follows. In Section 2, we introduce our three-field formulations and show their well-posedness. We propose finite element methods for both formulations and prove a priori error estimates in Section 3. Finally, a short conclusion is drawn in Section 4.

2. Two Three-Field Formulations of the Poisson Equation

In this section we introduce two three-field mixed formulations of the Poisson problem. Let Ω⊂ℝ^𝑑, 𝑑∈{2,3}, be a bounded convex domain with polygonal or polyhedral boundary 𝜕Ω with the outward pointing normal 𝐧 on 𝜕Ω. There are many mixed formulations of the Poisson equation −Δ𝑢=𝑓inΩ (1) with Dirichlet boundary condition 𝑢=0on𝜕Ω, (2) Neumann boundary condition 𝜕𝑢𝜕𝐧=0on𝜕Ω, (3) or a mixture of these two (see [2–6, 8, 10]).

We start with the following minimization problem for the Poisson problem: 𝐽(𝑢)=inf_{𝑣∈𝐻10(Ω)}𝐽(𝑣), (4) with 𝐽(𝑣)=12_Ω|∇𝑣|²𝑑𝑥−ℓ(𝑣), (5) where, for simplicity, we have assumed the Dirichlet boundary condition, 𝑓∈𝐿²(Ω), and ℓ(𝑣)=_Ω𝑓𝑣𝑑𝑥. (6) Let 𝑉=𝐻¹₀(Ω) and 𝑅=[𝐿²(Ω)]^𝑑, and for two vector-valued functions  𝜶∶Ω→ℝ^𝑑 and  𝜷∶Ω→ℝ^𝑑, the Sobolev inner product on the Sobolev space [𝐻^𝑘(Ω)]^𝑑 is defined as ⟨𝜶,𝜷⟩_𝑘,Ω∶=^𝑑_𝑖=1𝛼_𝑖,𝛽_𝑖𝑘,Ω, (7) where (𝜶)_𝑖=𝛼_𝑖and(𝜷)_𝑖=𝛽_𝑖 with 𝛼_𝑖,𝛽_𝑖∈𝐻^𝑘(Ω), for 𝑖=1,…,𝑑, and the norm ‖⋅‖_𝑘,Ω is induced from this inner product, where 𝑘∈ℕ∪{0}. Note that ⟨𝛼_𝑖,𝛽_𝑖⟩_𝑘,Ω is the standard inner product on the Sobolev space 𝐻^𝑘(Ω).

Our new formulation is obtained by introducing an auxiliary variable  𝝈=∇𝑢 such that the minimization problem (4) is rewritten as the following constrained minimization problem: argmin_{(𝑢,𝝈)∈[𝑉×𝑅]𝝈=∇𝑢}12‖𝝈‖²_0,Ω−ℓ(𝑢), (8) where the natural norm for the product space 𝑉×𝑅 is ‖𝑢‖²_1,Ω+‖𝝈‖²_0,Ω for (𝑢,𝝈)∈𝑉×𝑅. We write a weak variational equation for equation  𝝈=∇𝑢 in terms of the Lagrange multiplier space 𝑅 to obtain the saddle-point problem of the minimization problem (8). The saddle-point formulation is to find (𝑢,𝝈,𝝓)∈𝑉×𝑅×𝑅 such that 𝑎),())(𝑣),()𝑏))((𝑢,𝝈𝑣,𝝉))+𝑏((𝑣,𝝉,𝝓=ℓ𝑣,𝝉∈𝑉×𝑅,((𝑢,𝝈,𝝍=0,𝝍∈𝑅, (9) where 𝑎((𝑢,𝝈),(𝑣,𝝉))=_Ω𝑏))=𝝈⋅𝝉𝑑𝑥,((𝑢,𝝈,𝝍_Ω(𝝈−∇𝑢)⋅𝝍𝑑𝑥. (10) A similar three-field formulation—called the Hu-Washizu formulation—is very popular in designing finite element methods to alleviate locking effect in linear elasticity [11–13].

In order to show that the saddle-point problem (9) has a unique solution, we want to apply a standard saddle-point theory [3, 4, 10]. To this end, we need to show the following  three conditions of well-posedness. (1)The linear form ℓ(⋅), the bilinear forms 𝑎(⋅,⋅) and 𝑏(⋅,⋅) are continuous on the spaces on which they are defined. (2)The bilinear form 𝑎(⋅,⋅) is coercive on the kernel space 𝐾 defined as 𝐾={(𝑣,𝝉)∈𝑉×𝑅∶𝑏((𝑣,𝝉),𝝍)=0,𝝍∈𝑅}. (11)(3)The bilinear form 𝑏(⋅,⋅) satisfies the inf-sup condition: sup_{(𝑣,𝝉)∈𝑉×𝑅}𝑏((𝑣,𝝉),𝝍)‖𝑣‖²_1,Ω+‖𝝉‖²_0,Ω≥𝛽‖𝝍‖_0,Ω,𝝍∈𝑅 (12)for a constant 𝛽>0.

The linear form ℓ(⋅)  and the two bilinear forms 𝑎(⋅,⋅) and 𝑏(⋅,⋅) are continuous by the Cauchy-Schwarz inequality. The coercivity of the bilinear form 𝑎(⋅,⋅) on the kernel space 𝐾 follows from Poincaré inequality: 𝑎((𝑢,𝝈),(𝑢,𝝈))=‖𝝈‖²_0,Ω≥12‖𝝈‖²_0,Ω+‖𝑢‖²_1,Ω. (13) It remains to show that the bilinear form 𝑏(⋅,⋅) satisfies the inf-sup condition.

Lemma 1. The bilinear form 𝑏(⋅,⋅) satisfies the inf-sup condition; that is, sup_{(𝑣,𝝉)∈𝑉×𝑅}𝑏((𝑣,𝝉),𝝓)|𝑣‖²_1,Ω+‖𝝉‖²_0,Ω≥‖𝝓‖_0,Ω,𝝓∈𝑅. (14)

Proof. Using 𝑣=0 in the expression on the left above, we have sup_{(𝑣,𝝉)∈𝑉×𝑅}𝑏((𝑣,𝝉),𝝓)|𝑣‖²_1,Ω+‖𝝉‖²_0,Ω≥sup_𝝉∈𝑅∫_Ω𝝉⋅𝝓𝑑𝑥‖𝝉‖²_0,Ω=‖𝝓‖_0,Ω. (15)

To summarize we have proved the following theorem.

Theorem 2. The saddle-point problem (9) has a unique solution (𝑢,𝝈,𝝓)∈𝑉×𝑅×𝑅 and ‖𝑢‖_1,Ω+‖𝝈‖_0,Ω+‖𝝓‖_0,Ω≤‖ℓ‖_0,Ω. (16)

The main difficulty in this formulation is that the bilinear form 𝑎(⋅,⋅) is not coercive on the whole space 𝑉×𝑅. It is only coercive on the subspace 𝐾 of 𝑉×𝑅. One has to consider this difficulty in developing a finite element method for this problem. One approach to make the bilinear form 𝑎(⋅,⋅) coercive on the whole space 𝑉×𝑅 is to stabilize it. There are many approaches to stabilize this formulation (see [2, 5, 6]). Here we modify the bilinear form 𝑎(⋅,⋅) to get a consistent stabilization as proposed in [14] for the Mindlin-Reissner plate so that it is now coercive on the whole space 𝑉×𝑅. This is obtained by adding the term _Ω(𝝈−∇𝑢)⋅(𝝉−∇𝑣)𝑑𝑥 (17) to the bilinear form 𝑎(⋅,⋅). Thus our problem is to find (𝑢,𝝈,𝝓)∈𝑉×𝑅×𝑅 such that 𝐴),())(𝑣),()𝐵))((𝑢,𝝈𝑣,𝝉))+𝐵((𝑣,𝝉,𝝓=ℓ𝑣,𝝉∈𝑉×𝑅,((𝑢,𝝈,𝝍=0,𝝍∈𝑅, (18) where 𝐴((𝑢,𝝈),(𝑣,𝝉))=_Ω+𝝈.𝝉𝑑𝑥_Ω()⋅()𝐵))𝝈−∇𝑢𝝉−∇𝑣𝑑𝑥,((𝑢,𝝈,𝝍=_Ω(𝝈−∇𝑢)⋅𝝍𝑑𝑥. (19)

Theorem 3. The saddle-point problems (9) and (18) yield the same solution.

Proof. By construction the solution to (18) is also the solution to (9). Moreover, as the second equation of (18) yields 𝝈=∇𝑢, we can substitute this into the first equation and set 𝑣=0 to get 𝝓=−∇𝑢. This proves that the solution to (18) is also the solution to (9).

Here the bilinear form 𝐴(⋅,⋅) is coercive on the whole space 𝑉×𝑅 from the triangle and Poincaré inequalities.

Lemma 4. The bilinear form 𝐴(⋅,⋅) is coercive on 𝑉×𝑅. That is, 𝐴((𝑢,𝝈),(𝑢,𝝈))≥𝐶‖𝑢‖²_1,Ω+‖𝝈‖²_0,Ω,(𝑢,𝝈)∈𝑉×𝑅. (20)

Proof. We start with the triangle inequality ‖∇𝑢‖²_0,Ω≤‖∇𝑢−𝝈‖²_0,Ω+‖𝝈‖²_0,Ω. (21) From Poincaré inequality, there exists a constant 𝐶>0 such that 𝐶‖𝑢‖²_1,Ω≤‖∇𝑢‖²_0,Ω (22) and hence we have 𝐶‖𝑢‖²_1,Ω≤‖∇𝑢−𝝈‖²_0,Ω+‖𝝈‖²_0,Ω=𝐴((𝑢,𝝈),(𝑢,𝝈)). (23)

Moreover, the other conditions of well-posedness for this new saddle-point formulation (18) can also be shown exactly as for (9).

In the following, we present finite element approximations for both proposed three-field formulations. In the first step, we propose a finite element method for the nonstabilized formulation (9), which is based on quasi-biorthogonal systems. In the second step, we present a finite element method for the stabilized formulation (18), which is based on biorthogonal systems. The first method works only for simplicial elements, whereas the second method works for meshes of parallelotopes and simplices.

3. Finite Element Approximation and a Priori Error Estimate

Let 𝒯_ℎ be a quasi-uniform partition of the domain Ω in simplices or parallelotopes having the mesh size ℎ. Let 𝑇 be a reference simplex or 𝑑-cube, where the reference simplex is defined as 𝑇∶=𝐱∈ℝ^𝑑∶𝑥_𝑖>0,𝑖=1,…,𝑑,^𝑑_𝑖=1𝑥_𝑖<1, (24) and a 𝑑-cube 𝑇∶=(−1,1)^𝑑. Note that for 𝑑=2 a simplex is a triangle, and for 𝑑=2 it is a tetrahedron. Similarly, a 2-cube is a square, whereas, 3-cube is a cube.

The finite element space is defined by affine maps 𝐹_𝑇 from a reference element 𝑇 to a physical element 𝑇∈𝒯_ℎ. Let 𝒬(𝑇) be the space of bilinear/trilinear polynomials in 𝑇 if 𝑇 is a reference square/cube or the space of linear polynomials in 𝑇 if 𝑇 is a reference simplex. Then the finite element space based on the mesh 𝒯_ℎ is defined as the space of continuous functions whose restrictions to an element 𝑇 are obtained by maps of bilinear, trilinear, or linear functions from the reference element; that is, 𝑆_ℎ∶=𝑣_ℎ∈𝐻¹(Ω)∶𝑣_ℎ|_𝑇=𝑣_ℎ∘𝐹_𝑇⁻¹,𝑣_ℎ∈𝒬𝑇,𝑇∈𝒯_ℎ,𝑉_ℎ∶=𝑆_ℎ∩𝐻¹₀(Ω). (25) (See [3, 10, 15]). We start with some abstract assumptions for two discrete spaces [𝐿_ℎ]^𝑑⊂𝑅 and [𝑀_ℎ]^𝑑⊂𝑅 to discretize the gradient of the solution and the Lagrange multiplier spaces, where 𝐿_ℎ and 𝑀_ℎ both are piecewise polynomial spaces with respect to the mesh 𝒯_ℎ and are both subsets of 𝐿²(Ω). Explicit construction of local basis functions for these spaces for the mixed finite element method (9) is given in Section 3.1 and for the mixed finite element method (18) is given in Section 3.2.

We assume that 𝐿_ℎ and 𝑀_ℎ satisfy the following assumptions (see also [16, 17]).

Assumption 1. (i)  dim𝑀_ℎ=dim𝐿_ℎ.
(ii) There is a constant 𝛽>0 independent of the mesh 𝒯_ℎ such that 𝜙_ℎ𝐿2(Ω)≤𝛽sup_{𝜇ℎ∈𝑀ℎ⧵{0}}∫_Ω𝜇_ℎ𝜙_ℎ𝑑𝑥𝜇_ℎ𝐿2(Ω),𝜙_ℎ∈𝐿_ℎ. (26)
(iii) The space 𝑀_ℎ has the approximation property: inf_{𝜆ℎ∈𝑀ℎ}𝜙−𝜆_ℎ𝐿2(Ω)≤𝐶ℎ𝜙_𝐻1(Ω),𝜙∈𝐻¹(Ω). (27)
(iv) The space 𝐿_ℎ has the approximation property: inf_{𝜙ℎ∈𝐿ℎ}𝜙−𝜙_ℎ𝐿2(Ω)≤𝐶ℎ𝜙_𝐻1(Ω),𝜙∈𝐻¹(Ω). (28)
Now we define biorthogonality and quasi-biorthogonality (see [16, 18, 19]).

Definition 5. The pair of bases {𝜇_𝑖}_{1≤𝑖≤𝑛} of 𝑀_ℎ and {𝜑_𝑖}_{1≤𝑖≤𝑛} of 𝐿_ℎ is called  biorthogonal if the resulting Gram matrix 𝙶 is diagonal. The pair of bases {𝜇_𝑖}_{1≤𝑖≤𝑛} of 𝑀_ℎ and {𝜑_𝑖}_{1≤𝑖≤𝑛} of 𝐿_ℎ is called  quasi-biorthogonal and the resulting Gram matrix  𝙶  is called  quasidiagonal if  𝙶  is of the form 𝙳𝙶=¹0𝚁𝙳₂𝙳or𝙶=¹𝚁0𝙳₂, (29) where 𝙳₁ and 𝙳₂ are diagonal matrices and 𝚁 is a sparse rectangular matrix.

We recall that a Gram matrix 𝙶 of two sets of basis functions {𝜇_𝑖}_{1≤𝑖≤𝑛} and {𝜑_𝑖}_{1≤𝑖≤𝑛} is the matrix 𝙶 whose 𝑖𝑗th entry is _Ω𝜇_𝑖𝜑_𝑗𝑑𝑥. (30) It is worth noting that the quasi-diagonal matrix is inverted by inverting two diagonal matrices and scaling the rectangular matrix.

3.1. Finite Element Approximation for the Nonstabilized Formulation (9)

We now turn our attention to a finite element approximation for the nonstabilized formulation (9). This approximation works only for simplicial meshes and is based on a quasi-biorthogonal system as in [19]. Let 𝐵_ℎ=𝑏_ℎ∣𝑏_ℎ|𝑇=𝐶^𝑑+1_𝑖=1𝜆^𝑇_𝑖,𝑇∈𝒯_ℎ (31) be the space of bubble functions, where 𝐶 is a constant and {𝜆^𝑇_𝑖}^𝑑+1_𝑖=1 is the set of barycentric coordinates on 𝑇. Let 𝐿_ℎ=𝑆_ℎ⊕𝐵_ℎ. Explicit construction of basis functions of 𝐿_ℎ and 𝑀_ℎ on a reference element is provided below.

The finite element problem is to find (𝑢_ℎ,𝝈_ℎ,𝝓_ℎ)∈𝑉_ℎ×[𝐿_ℎ]^𝑑×[𝑀_ℎ]^𝑑 such that 𝑎𝑢_ℎ,𝝈_ℎ,𝑣_ℎ,𝝉_ℎ+𝑏𝑣_ℎ,𝝉_ℎ,𝝓_ℎ=ℓ𝑣_ℎ,𝑣_ℎ,𝝉_ℎ∈𝑉_ℎ×𝐿_ℎ^𝑑,𝑏𝑢_ℎ,𝝈_ℎ,𝝍_ℎ=0,𝝍_ℎ∈𝑀_ℎ^𝑑. (32) We note that the standard finite element space 𝑆_ℎ is enriched with elementwise defined bubble functions to obtain the space 𝐿_ℎ, which is done to obtain the coercivity of the bilinear form 𝑎(⋅,⋅) on the kernel space 𝐾_ℎ defined as (see Lemma 7) 𝐾_ℎ=𝑣_ℎ,𝝉_ℎ∈𝑉_ℎ×𝐿_ℎ^𝑑∶𝑏𝑣_ℎ,𝝉_ℎ,𝝍_ℎ=0,𝝍_ℎ∈𝑀_ℎ^𝑑. (33) To simplify the numerical analysis of the finite element problem, we introduce a quasiprojection operator: 𝑄_ℎ∶𝐿²(Ω)→𝐿_ℎ, which is defined as  𝝍_ℎΩ𝑄_ℎ𝑣𝜇_ℎ𝑑𝑥=_Ω𝑣𝜇_ℎ𝑑𝑥,𝑣∈𝐿²(Ω),𝜇_ℎ∈𝑀_ℎ. (34) This type of operator was introduced in [9] to obtain the finite element interpolation of nonsmooth functions satisfying boundary conditions and is used in [20] in the context of mortar finite elements. The definition of 𝑄_ℎ allows us to write the weak gradient as 𝝈_ℎ=𝑄_ℎ∇𝑢_ℎ, (35) where operator 𝑄_ℎ is applied to the vector ∇𝑢_ℎ component-wise. We see that 𝑄_ℎ is well defined due to Assumption 1(i) and (ii). Furthermore, the restriction of 𝑄_ℎ to 𝐿_ℎ is the identity. Hence 𝑄_ℎ is a projection onto the space 𝐿_ℎ. We note that 𝑄_ℎ is not the orthogonal projection onto 𝐿_ℎ but an oblique projection onto 𝐿_ℎ see [21, 22]. The stability and approximation properties of 𝑄_ℎ in 𝐿² and 𝐻¹-norm can be shown as in [16, 23]. In the following, we will use a generic constant 𝐶, which will take different values at different places but will be always independent of the mesh size ℎ.

Lemma 6. Under Assumption 1(i)–(iii) 𝑄_ℎ𝑣_𝐿2(Ω)≤𝐶‖𝑣‖_𝐿²_(Ω),𝑓𝑜𝑟𝑣∈𝐿²(Ω),𝑄_ℎ𝑤_𝐻1(Ω)≤𝐶|𝑤|_𝐻¹_(Ω),𝑓𝑜𝑟𝑤∈𝐻¹(Ω), (36) and for 0<𝑠≤1 and 𝑣∈𝐻^1+𝑠(Ω)𝑣−𝑄_ℎ𝑣_𝐿2(Ω)≤𝐶ℎ^1+𝑠|𝑣|_𝐻^𝑠+1_(Ω),𝑣−𝑄_ℎ𝑣_𝐻1(Ω)≤𝐶ℎ^𝑠|𝑣|_𝐻^𝑠+1_(Ω). (37)

We need to construct another finite element space 𝑀_ℎ⊂𝐿²(Ω) that satisfies Assumption 1(i)–(iii) and also leads to an efficient numerical scheme. The main goal is to choose the basis functions for 𝑀_ℎ and 𝐿_ℎ so that the matrix 𝙳 associated with the bilinear form ∫_Ω𝝈_ℎ∶𝝍_ℎ𝑑𝑥 is easy to invert. To this end, we consider the algebraic form of the weak equation for the gradient of the solution, which is given by −𝙱𝑢+𝙳𝝈=0, (38) where 𝑢 and  𝝈are the solution vectors for 𝑢_ℎ and  𝝈_ℎ and 𝙳 is the Gram matrix between the bases of [𝑀_ℎ]^𝑑 and [𝐿_ℎ]^𝑑. Hence if the matrix 𝙳 is diagonal or quasi-diagonal, we can easily eliminate the degrees of freedom corresponding to  𝝈_ℎ and  𝝓_ℎ. This then leads to a formulation involving only one unknown 𝑢_ℎ. Statically condensing out variables 𝝈_ℎ and 𝝓_ℎ we arrive at the variational formulation of the reduced system given by (68) (see Section 3.3). Note that the matrix 𝙳 is the Gram matrix between the bases of [𝑀_ℎ]^𝑑 and [𝐿_ℎ]^𝑑.

We cannot use a biorthogonal system in this case as the stability will be lost (see the Appendix). This motivates our interest in a quasi-biorthogonal system. We now show the construction of the local basis functions for 𝐿_ℎ and 𝑀_ℎ so that they form a quasi-biorthogonal system. The construction is shown on the reference triangle 𝑇={(𝑥,𝑦)∶0≤𝑥,0≤𝑦,𝑥+𝑦≤1} in the two-dimensional case and on the reference tetrahedron 𝑇={(𝑥,𝑦,𝑧)∶0≤𝑥,0≤𝑦,0≤𝑧,𝑥+𝑦+𝑧≤1} in the three-dimensional case.

We now start with the reference triangle in the two-dimensional case. Let {𝜑₁,…,𝜑₄} defined as 𝜑₁=1−𝑥−𝑦,𝜑₂=𝑥,𝜑₃𝜑=𝑦,₄=27𝑥𝑦1−𝑥−𝑦 (39) be the local basis functions of 𝐿_ℎ on the reference triangle associated with three vertices (0,0), and (1,0), (0,1) and one barycenter (1/3,1/3). Note that 𝜑₄ is the standard bubble function used, for example, in the minielement for the Stokes equations [24].

Then if we define the local basis functions {𝜇₁,…,𝜇₄} of 𝑀_ℎ as 𝜇₁=349−4𝑥−4𝑦−1403,𝜇𝑥𝑦1−𝑥−𝑦₂2=−9+4𝑥−1403,𝜇𝑥𝑦1−𝑥−𝑦₃2=−9+4𝑦−1403,𝜇𝑥𝑦1−𝑥−𝑦₄=27𝑥𝑦1−𝑥−𝑦, (40) the two sets of basis functions {𝜑₁,…,𝜑₄} and {𝜇₁,…,𝜇₄} form a quasi-biorthogonal system on the reference triangle. These four basis functions {𝜇₁,…,𝜇₄} of 𝑀_ℎ are associated with three vertices (0,0), and (1,0), (0,1) and the barycenter (1/3,1/3) of the reference triangle.

Similarly, on the reference tetrahedron, the local basis functions {𝜑₁,…,𝜑₅} of 𝐿_ℎ defined as 𝜑₁=1−𝑥−𝑦−𝑧,𝜑₂=𝑥,𝜑₃𝜑=𝑦,₄=𝑧,𝜑₅=256𝑥𝑦𝑧1−𝑥−𝑦−𝑧 (41) and the local basis functions {𝜇₁,…,𝜇₅} of 𝑀_ℎ defined as 𝜇₁=40192−5𝑥−5𝑦−5𝑧−6930,𝜇23𝑥𝑦𝑧1−𝑥−𝑦−𝑧₂=−5992+5𝑥−6930,𝜇23𝑥𝑦𝑧1−𝑥−𝑦−𝑧₃=−5992+5𝑦−6930,𝜇23𝑥𝑦𝑧1−𝑥−𝑦−𝑧₄=−5992+5𝑧−6930,𝜇23𝑥𝑦𝑧1−𝑥−𝑦−𝑧₅=256𝑥𝑦𝑧1−𝑥−𝑦−𝑧 (42) are quasi-biorthogonal. These five basis functions {𝜑₁,…,𝜑₅} of 𝐿_ℎ or {𝜇₁,…,𝜇₅} of 𝑀_ℎ are associated with four vertices (0,0,0), (1,0,0), (0,1,0), and (0,0,1) and the barycenter (1/4,1/4,and1/4) of the reference tetrahedron. Using the finite element basis functions for 𝐿_ℎ, the finite element basis functions for 𝑀_ℎ are constructed in such a way that the Gram matrix on the reference element 𝑇 is given by 𝙳_𝑇=16010006100006033403404081560 (43) for the two-dimensional case and 𝙳_𝑇=101240000124000001240000004244315431543153154096155925 (44) for the three-dimensional case. After ordering the degrees of freedom in 𝝈_ℎ and 𝝓_ℎ so that the degrees of freedom associated with the barycenters of elements come after the degrees of freedom associated with the vertices of elements, the global Gram matrix 𝙳 is quasi-diagonal.

Thus in the two-dimensional case, the local basis functions 𝜇₁,𝜇₂, and 𝜇₃ are associated with the vertices of the reference triangle, and the function 𝜇₄ is associated with the barycenter of the triangle and defined elementwise. That means 𝜇₄ is supported only on the reference element. The situation is similar in the three-dimensional case. Hence using a proper ordering, the support of a global basis function 𝜑_𝑖 of 𝐿_ℎ is the same as that of the global basis function 𝜇_𝑖 of 𝑀_ℎ for 𝑖=1,…,𝑛. However, the global basis functions 𝜇₁,…,𝜇_𝑛 of 𝑀_ℎ are not continuous on interelement boundaries.

Now we show that 𝑀_ℎ satisfies Assumption 1(i)–(iii). As the first assumption is satisfied by construction, we consider the second assumption. Let 𝜙_ℎ=∑^𝑛_𝑘=1𝑎_𝑘𝜙_𝑘∈𝑆_ℎ and set 𝜇_ℎ=∑^𝑛_𝑘=1𝑎_𝑘𝜇_𝑘∈𝑀_ℎ. By using the quasibiorthogonality relation (1) and the quasiuniformity assumption, we get _Ω𝜙_ℎ𝜇_ℎ𝑑𝐱=^𝑛_𝑖,𝑗=1𝑎_𝑖𝑎_𝑗Ω𝜙_𝑖𝜇_𝑗=𝑑𝐱_𝑛𝑖=1𝑎²_𝑖𝑐_𝑖≥𝐶^𝑛_𝑖=1𝑎²_𝑖ℎ^𝑑_𝑖≥𝐶𝜙_ℎ²_0,Ω, (45) where ℎ_𝑖 denotes the meshsize at 𝑖th vertex. Taking into account the fact that ‖𝜙_ℎ‖²_0,Ω≡‖𝜇_ℎ‖²_0,Ω≡∑^𝑛_𝑖=1𝑎²_𝑖ℎ^𝑑_𝑖, we find that Assumption 1(ii) is satisfied. Since the sum of the local basis functions of 𝑀_ℎ is one, Assumption 1(iii) can be proved as in [23, 25]. Assumption 1(iv) follows by standard arguments.

There are (𝑑+2) local basis functions of 𝐿_ℎ, where (𝑑+1) of them are associated with vertices and one is associated with the element barycenter. Then the local basis functions of 𝑀_ℎ are also divided into two parts: (𝑑+2) basis functions associated with the vertices and 1 basis function associated with the element. Let 𝜇_𝑇 be the local basis function associated with the element barycenter of 𝑇∈𝒯_ℎ for 𝑀_ℎ. Let 𝑌_ℎ be the space of piecewise constant functions with respect to 𝒯_ℎ, and 𝑃_ℎ∶𝐿²(Ω)→𝑌_ℎ is defined as 𝑃_ℎ𝑤∣_𝑇=∫_𝑇𝑤𝜇_𝑇𝑑𝑥∫_𝑇𝜇_𝑇𝑑𝑥. (46) Then 𝑃_ℎ is a projection operator onto 𝑌_ℎ. Moreover, ∇𝑣_ℎ=𝑃_ℎ𝑄_ℎ∇𝑣_ℎ, (47) and this implies that ∇𝑣_ℎ0,Ω=𝑃_ℎ𝑄_ℎ∇𝑣_ℎ0,Ω≤𝑄_ℎ∇𝑣_ℎ0,Ω. (48) This yields the following coercivity result.

Lemma 7. The bilinear form 𝑎(⋅,⋅) is coercive on the kernel space 𝐾_ℎ. In fact, 𝑎𝑢_ℎ,𝝈_ℎ,𝑢_ℎ,𝝈_ℎ≥12𝝈_ℎ²_0,Ω+∇𝑢_ℎ²_0,Ω,𝑢_ℎ,𝝈_ℎ∈𝐾_ℎ. (49)

Proof. Since 𝝈_ℎ=𝑄_ℎ∇𝑢_ℎ on the kernel space 𝐾_ℎ, we have 𝑎𝑢_ℎ,𝝈_ℎ,𝑢_ℎ,𝝈_ℎ=12𝝈_ℎ²_0,Ω+𝑄_ℎ∇𝑢_ℎ²_0,Ω. (50) The proof follows from the result ∇𝑢_ℎ0,Ω=𝑃_ℎ𝑄_ℎ∇𝑢_ℎ0,Ω≤𝑄_ℎ∇𝑢_ℎ0,Ω. (51)

Remark 8. Using a quadrilateral or hexahedral meshes, we cannot have ∇𝑣_ℎ=𝑃_ℎ𝑄_ℎ∇𝑣_ℎ, and hence the coercivity fails. One needs to modify the finite element method for these meshes.

The inf-sup condition follows exactly as in the continuous setting using Assumption 1(ii).

Lemma 9. There exists a constant 𝛽^∼>0 such that sup_{(𝑣ℎ,𝝉ℎ)∈𝑉ℎ×𝑅ℎ}𝑏𝑣_ℎ,𝝉_ℎ,𝝓_ℎ𝑣_ℎ²_1,Ω+𝝉_ℎ²_0,Ω≥𝛽^∼𝝓_ℎ0,Ω,𝝓_ℎ∈𝑅_ℎ. (52)

Proof. The proof follows exactly as in the continuous setting. We use 𝑣_ℎ=0 in the expression sup_{(𝑣ℎ,𝝉ℎ)∈𝑉ℎ×𝑅ℎ}𝑏𝑣_ℎ,𝝉_ℎ,𝝓_ℎ𝑣_ℎ²_1,Ω+𝝉_ℎ²_0,Ω (53) to get sup_{(𝑣ℎ,𝝉ℎ)∈𝑉ℎ×𝑅ℎ}𝑏𝑣_ℎ,𝝉_ℎ,𝝓_ℎ𝑣_ℎ²_1,Ω+𝝉_ℎ²_0,Ω≥sup_{𝝉ℎ∈𝑅ℎ}∫_Ω𝝉_ℎ⋅𝝓_ℎ𝑑𝑥𝝉_ℎ²_0,Ω. (54) Assumption 1(ii) then yields the final result.

Hence we have the following stability and approximation result.

Theorem 10. The discrete saddle-point problem (32) has a unique solution (𝑢_ℎ,𝝈_ℎ,𝝓_ℎ)∈𝑉_ℎ×[𝐿_ℎ]^𝑑×[𝑀_ℎ]^𝑑, and 𝑢_ℎ1,Ω+𝝈_ℎ0,Ω+𝝓_ℎ0,Ω≤‖ℓ‖_0,Ω. (55) Moreover, the following estimate holds for the solution (𝑢,𝝈,𝝓)∈𝑉×𝑅×𝑅 of (32) 𝑢−𝑢_ℎ1,Ω+𝝈−𝝈_ℎ0,Ω+𝝓−𝝓_ℎ0,Ω≤𝐶inf_{𝑣ℎ∈𝑉ℎ}𝑢−𝑢_ℎ1,Ω+inf_{𝝉ℎ∈[𝐿ℎ]𝑑}𝝈−𝝉_ℎ0,Ω+inf_{𝝍ℎ∈[𝑀ℎ]𝑑}𝝓−𝝍_ℎ0,Ω. (56)

Owing to Assumption 1(iii)-(iv) the discrete solution converges optimally to the continuous solution.

3.2. Finite Element Approximation for the Stabilized Formulation (18)

After stabilizing formulation (9), the bilinear form 𝐴(⋅,⋅) (18) is coercive on the whole product space 𝑉×𝑅. Therefore, we can use the standard finite element space 𝑆_ℎ to discretize the gradient of the solution for the saddle-point problem (18) so that 𝐿_ℎ=𝑆_ℎ. Thus our discrete saddle-point formulation of (18) is to find (𝑢_ℎ,𝝈_ℎ,𝝓_ℎ)∈𝑉_ℎ×[𝐿_ℎ]^𝑑×[𝑀_ℎ]^𝑑 so that 𝐴𝑢_ℎ,𝝈_ℎ,𝑣_ℎ,𝝉_ℎ+𝐵𝐯_ℎ,𝝉_ℎ,𝝓_ℎ=ℓ𝑣_ℎ,𝑣_ℎ,𝝉_ℎ∈𝑉_ℎ×𝐿_ℎ^𝑑,𝐵𝑢_ℎ,𝝈_ℎ,𝝍_ℎ=0,𝝍_ℎ∈𝑀_ℎ^𝑑. (57) We can now verify the three conditions of well-posedness for the discrete saddle point problem. The continuity of 𝐴(⋅,⋅), 𝐵(⋅,⋅), and ℓ(⋅) follows as in the continuous setting, and coercivity is immediate from Lemma 4 as 𝑉_ℎ×[𝐿_ℎ]^𝑑⊂𝑉×𝑅. Explicit construction of basis functions of 𝑀_ℎ on a reference element will be given below.

Assumption 1(ii) guarantees that the bilinear form 𝑏(⋅,⋅) satisfies the inf-sup condition as in the case of the previous finite element method.

We have thus proved the following theorem from the standard theory of saddle-point problem  (see [4]).

Theorem 11. The discrete saddle-point problem (57) has a unique solution (𝑢_ℎ,𝝈_ℎ,𝝓_ℎ)∈𝑉_ℎ×[𝐿_ℎ]^𝑑×[𝑀_ℎ]^𝑑 and 𝑢_ℎ1,Ω+𝝈_ℎ0,Ω+𝝓_ℎ0,Ω≤‖ℓ‖_0,Ω. (58) Moreover, the following estimate holds for the solution (𝑢,𝝈,𝝓)∈𝑉×𝑅×𝑅 of (9) 𝑢−𝑢_ℎ1,Ω+𝝈−𝝈_ℎ0,Ω+𝝓−𝝓_ℎ0,Ω≤𝐶inf_{𝑣ℎ∈𝑉ℎ}𝑢−𝑢_ℎ1,Ω+inf_{𝝉ℎ∈[𝐿ℎ]𝑑}𝝈−𝝉_ℎ0,Ω+inf_{𝝍ℎ∈[𝑀ℎ]𝑑}𝝓−𝝍_ℎ0,Ω. (59)

The optimal convergence of the discrete solution is then guaranteed by the standard approximation property of 𝑉_ℎ and 𝐿_ℎ and Assumption 1(iii).

In the following, we give these basis functions for linear simplicial finite elements in two and three dimensions [18]. The parallelotope case is handled by using the tensor product construction of the one-dimensional basis functions presented in [20, 26]. These basis functions are very useful in the mortar finite element method and extensively studied in one- and two-dimensional cases [16, 23, 26].

Note that we have imposed the condition dim𝑀_ℎ=dim𝐿_ℎ to get that our Gram matrix is a square matrix. Therefore, the local basis functions for linear/bilinear/trilinear finite spaces are associated with the vertices of triangles, tetrahedra, quadrilaterals, and hexahedra. For example, for the reference triangle 𝑇∶={(𝑥,𝑦)∶0<𝑥,0<𝑦,𝑥+𝑦<1}, we have 𝜇₁∶=3−4𝑥−4𝑦,𝜇₂𝜇∶=4𝑥−1,₃∶=4𝑦−1, (60) where the basis functions 𝜇₁, 𝜇₂ and 𝜇₃ are associated with three vertices (0,0), (1,0), and (0,1) of the reference triangle. For the reference tetrahedron 𝑇∶={(𝑥,𝑦,𝑧)∶0<𝑥,0<𝑦,0<𝑧,𝑥+𝑦+𝑧<1}, we have 𝜇₁∶=4−5𝑥−5𝑦−5𝑧,𝜇₂𝜇∶=5𝑥−1,₃∶=5𝑦−1,𝜇₄∶=5𝑧−1, (61) where the basis functions 𝜇₁, 𝜇₂, 𝜇₃, and 𝜇₄, associated with four vertices (0,0,0), (1,0,0), (0,1,0), and (0,0,1) of the reference tetrahedron (see also [17, 18]). Note that the basis functions are constructed in such a way that the biorthogonality condition is satisfied on the reference element.

The global basis functions for the test space are constructed by gluing the local basis functions together following exactly the same procedure of constructing global finite element basis functions from the local ones. These global basis functions then satisfy the condition of biorthogonality with global finite element basis functions of 𝑆_ℎ. These basis functions also satisfy Assumption 1(i)–(iv) (see [17, 18]).

3.3. Static Condensation

Here we want to eliminate  𝝈_ℎ and Lagrange multiplier  𝝓_ℎ from the saddle-point problem (57). We make use of the operator 𝑄_ℎ defined previously. Using the biorthogonality relation between the basis functions of 𝐿_ℎ and 𝑀_ℎ, the action of operator 𝑄_ℎ on a function 𝑣∈𝐿²(Ω) can be written as 𝑄_ℎ𝑣=^𝑛_𝑖=1∫_Ω𝜇_𝑖𝑣𝑑𝑥𝑐_𝑖𝜑_𝑖, (62) which tells that the operator 𝑄_ℎ is local in the sense to be given below (see also [27]).

Using the property of operator 𝑄_ℎ, we can eliminate the degrees of freedom corresponding to  𝝈_ℎ and  𝝓_ℎ so that our problem is to find 𝑢_ℎ∈𝑆_ℎ such that 𝐽𝑢_ℎ=min_{𝑣ℎ∈𝑉ℎ}𝐽𝑣_ℎ, (63) where 𝐽𝑣_ℎ=∇𝑄_ℎ∇𝑣_ℎ²_0,Ω−2ℓ𝑣_ℎ (64) for the nonstabilized formulation (9) and 𝐽𝑣_ℎ=∇𝑄_ℎ∇𝑣_ℎ²_0,Ω+𝑄_ℎ∇𝑣_ℎ−∇𝑣_ℎ²_0,Ω−2ℓ𝑣_ℎ (65) for the stabilized formulation (18).

Let the bilinear form 𝐴^∼(⋅,⋅) be defined as 𝐴^∼𝑢_ℎ,𝑣_ℎ=_Ω𝝈_ℎ∶𝝉_ℎ𝑑𝑥 (66) for the non-stabilized formulation (9) and 𝐴^∼𝑢_ℎ,𝑣_ℎ=_Ω𝝈_ℎ∶𝝉_ℎ+𝑑𝑥_Ω𝝈_ℎ−∇𝑢_ℎ⋅𝝉_ℎ−∇𝑣_ℎ𝑑𝑥 (67) for the stabilized formulation (18) with  𝝈_ℎ=𝑄_ℎ(∇𝑢_ℎ) and  𝝉_ℎ=𝑄_ℎ(∇𝑣_ℎ). Due to the biorthogonality or quasi-biorthogonality condition, the action of 𝑄_ℎ is computed by inverting a diagonal matrix.

Since the bilinear form 𝐴^∼(⋅,⋅) is symmetric, the minimization problem (63) is equivalent to the variational problem of finding 𝑢_ℎ∈𝑆_ℎ such that [3, 15] 𝐴^∼𝑢_ℎ,𝑣_ℎ=ℓ𝑣_ℎ,𝑣_ℎ∈𝑉_ℎ. (68)

The coercivity of 𝐴(⋅,⋅) and the stability result of Lemma 9 immediately yield the coercivity of the bilinear form 𝐴^∼(⋅,⋅) on the space 𝑉_ℎ with the norm ‖𝑢_ℎ‖²_0,Ω+‖∇𝑄_ℎ(∇𝑢_ℎ)‖²_0,Ω for 𝑢_ℎ∈𝑉_ℎ. That is, 𝐴^∼𝑢_ℎ,𝑢_ℎ≥𝐶𝑢_ℎ²_0,Ω+∇𝑄_ℎ∇𝑢_ℎ²_0,Ω,𝑢_ℎ∈𝑉_ℎ. (69)

4. Conclusion

We have considered two three-field formulations of Poisson problem and shown the well-posedness of them. We have proposed efficient finite element schemes for both of these formulations: one for the stabilized three-field formulation based on biorthogonal systems and the other for the non-stabilized three-field formulation based on quasi-biorthogonal systems. Optimal a priori error estimates are proved for both approaches, and a reduced discrete problem is presented.

Appendix

If we want to enrich the finite element space with the bubble function to attain stability, it is not possible to use a biorthogonal system. We consider the local basis functions of 𝐿_ℎ for the reference triangle 𝑇={(𝑥,𝑦)∶0≤𝑥,0≤𝑦,𝑥+𝑦≤1}. Let {𝜑₁,…,𝜑₄} be defined as in (39). Then using the biorthogonal relation for the four local basis functions {𝜇₁,…,𝜇₄} of 𝑀_ℎ with the four basis functions {𝜑₁,…,𝜑₄} of 𝐿_ℎ for the triangle and ∑³_𝑖=1𝜑_𝑖=1 on 𝑇, 𝜇₄ should satisfy _𝑇𝜇₄³_𝑖=1𝜑_𝑖𝑑𝑥=_𝑇𝜇₄𝑑𝑥=0 (A.1) since it is associated with the barycenter of the triangle. This leads to instability in the formulation as we suggest that the function 𝜇₄ should have nonzero integral over 𝑇 (see (46)). This explains why we cannot use a biorthogonal system. We have a similar issue in the three-dimensional case as ∑⁴_𝑖=1𝜑_𝑖=1, resulting in _𝑇𝜇₅𝑑𝑥=0 (A.2) if we impose biorthogonality.

Acknowledgments

Support from the new staff Grant of the University of Newcastle is gratefully acknowledged. The author is grateful to the anonymous referees for their valuable suggestions to improve the quality of the earlier version of this work.