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

Authors: Bishnu P. Lamichhane
Journal: Advances in Numerical Analysis
 (2013)

---

## 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.


---

## Body

## 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 Section2, 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Ω⊂ℝd, d∈{2,3}, be a bounded convex domain with polygonal or polyhedral boundary ∂Ω with the outward pointing normal n on ∂Ω. There are many mixed formulations of the Poisson equation
(1)-Δu=finΩ
with Dirichlet boundary condition
(2)u=0on∂Ω,
Neumann boundary condition
(3)∂u∂n=0on∂Ω,
or a mixture of these two (see [2–6, 8, 10]).We start with the following minimization problem for the Poisson problem:(4)J(u)=infv∈H01(Ω)J(v),
with
(5)J(v)=12∫Ω‍|∇v|2dx-ℓ(v),
where, for simplicity, we have assumed the Dirichlet boundary condition, f∈L2(Ω), and
(6)ℓ(v)=∫Ω‍fvdx.
Let V=H01(Ω) and R=[L2(Ω)]d, and for two vector-valued functions  α:Ω→ℝd and  β:Ω→ℝd, the Sobolev inner product on the Sobolev space [Hk(Ω)]d is defined as
(7)〈α,β〉k,Ω:=∑i=1d‍〈αi,βi〉k,Ω,
where (α)i=αiand(β)i=βi with αi,βi∈Hk(Ω), for i=1,…,d, and the norm ∥·∥k,Ω is induced from this inner product, where k∈ℕ∪{0}. Note that 〈αi,βi〉k,Ω is the standard inner product on the Sobolev space Hk(Ω).Our new formulation is obtained by introducing an auxiliary variableσ=∇u such that the minimization problem (4) is rewritten as the following constrained minimization problem:
(8)arg min(u,σ)∈[V×R]σ=∇u(12∥σ∥0,Ω2-ℓ(u)),
where the natural norm for the product space V×R is ∥u∥1,Ω2+∥σ∥0,Ω2 for (u,σ)∈V×R. We write a weak variational equation for equation  σ=∇u in terms of the Lagrange multiplier space R to obtain the saddle-point problem of the minimization problem (8). The saddle-point formulation is to find (u,σ,ϕ)∈V×R×R such that
(9)a((u,σ),(v,τ))+b((v,τ),ϕ)=ℓ(v),(v,τ)∈V×R,b((u,σ),ψ)=0,ψ∈R,
where
(10)a((u,σ),(v,τ))=∫Ω‍σ·τdx,b((u,σ),ψ)=∫Ω‍(σ-∇u)·ψdx.
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 a(·,·) and b(·,·) are continuous on the spaces on which they are defined.(2)
The bilinear forma(·,·) is coercive on the kernel space K defined as
(11)K={(v,τ)∈V×R:b((v,τ),ψ)=0,ψ∈R}.(3)
The bilinear formb(·,·) satisfies the inf-sup condition:
(12)sup(v,τ)∈V×Rb((v,τ),ψ)∥v∥1,Ω2+∥τ∥0,Ω2≥β∥ψ∥0,Ω,ψ∈Rfor a constant β>0.The linear formℓ(·)  and the two bilinear forms a(·,·) and b(·,·) are continuous by the Cauchy-Schwarz inequality. The coercivity of the bilinear form a(·,·) on the kernel space K follows from Poincaré inequality:
(13)a((u,σ),(u,σ))=∥σ∥0,Ω2≥12(∥σ∥0,Ω2+∥u∥1,Ω2).
It remains to show that the bilinear form b(·,·) satisfies the inf-sup condition.Lemma 1.
The bilinear formb(·,·) satisfies the inf-sup condition; that is,
(14)sup(v,τ)∈V×Rb((v,τ),ϕ)|v∥1,Ω2+∥τ∥0,Ω2≥∥ϕ∥0,Ω,ϕ∈R.Proof.
Usingv=0 in the expression on the left above, we have
(15)sup(v,τ)∈V×Rb((v,τ),ϕ)|v∥1,Ω2+∥τ∥0,Ω2≥supτ∈R∫Ω‍τ·ϕdx∥τ∥0,Ω2=∥ϕ∥0,Ω.To summarize we have proved the following theorem.Theorem 2.
The saddle-point problem (9) has a unique solution (u,σ,ϕ)∈V×R×R and
(16)∥u∥1,Ω+∥σ∥0,Ω+∥ϕ∥0,Ω≤∥ℓ∥0,Ω.The main difficulty in this formulation is that the bilinear forma(·,·) is not coercive on the whole space V×R. It is only coercive on the subspace K of V×R. One has to consider this difficulty in developing a finite element method for this problem. One approach to make the bilinear form a(·,·) coercive on the whole space V×R is to stabilize it. There are many approaches to stabilize this formulation (see [2, 5, 6]). Here we modify the bilinear form a(·,·) to get a consistent stabilization as proposed in [14] for the Mindlin-Reissner plate so that it is now coercive on the whole space V×R. This is obtained by adding the term
(17)∫Ω‍(σ-∇u)·(τ-∇v)dx
to the bilinear form a(·,·). Thus our problem is to find (u,σ,ϕ)∈V×R×R such that
(18)A((u,σ),(v,τ))+B((v,τ),ϕ)=ℓ(v),(v,τ)∈V×R,B((u,σ),ψ)=0,ψ∈R,
where
(19)A((u,σ),(v,τ))=∫Ω‍σ.τdx+∫Ω‍(σ-∇u)·(τ-∇v)dx,B((u,σ),ψ)=∫Ω‍(σ-∇u)·ψdx.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 σ=∇u, we can substitute this into the first equation and set v=0 to get ϕ=-∇u. This proves that the solution to (18) is also the solution to (9).Here the bilinear formA(·,·) is coercive on the whole space V×R from the triangle and Poincaré inequalities.Lemma 4.
The bilinear formA(·,·) is coercive on V×R. That is,
(20)A((u,σ),(u,σ))≥C(∥u∥1,Ω2+∥σ∥0,Ω2),(u,σ)∈V×R.Proof.
We start with the triangle inequality(21)∥∇u∥0,Ω2≤∥∇u-σ∥0,Ω2+∥σ∥0,Ω2.
From Poincaré inequality, there exists a constant C>0 such that
(22)C∥u∥1,Ω2≤∥∇u∥0,Ω2
and hence we have
(23)C∥u∥1,Ω2≤(∥∇u-σ∥0,Ω2+∥σ∥0,Ω2)=A((u,σ),(u,σ)).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𝒯h be a quasi-uniform partition of the domain Ω in simplices or parallelotopes having the mesh size h. Let T^ be a reference simplex or d-cube, where the reference simplex is defined as
(24)T^:={x∈ℝd:xi>0,i=1,…,d,∑i=1d‍xi<1},
and a d-cube T^:=(-1,1)d. Note that for d=2 a simplex is a triangle, and for d=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 mapsFT from a reference element T^ to a physical element T∈𝒯h. Let 𝒬(T^) be the space of bilinear/trilinear polynomials in T^ if T^ is a reference square/cube or the space of linear polynomials in T^ if T^ is a reference simplex. Then the finite element space based on the mesh 𝒯h is defined as the space of continuous functions whose restrictions to an element T are obtained by maps of bilinear, trilinear, or linear functions from the reference element; that is,
(25)Sh:={vh∈H1(Ω):vh|TSh:v==v^h∘FT-1,v^h∈𝒬(T^),T∈𝒯h},Vh:=Sh∩H01(Ω).
(See [3, 10, 15]). We start with some abstract assumptions for two discrete spaces [Lh]d⊂R and [Mh]d⊂R to discretize the gradient of the solution and the Lagrange multiplier spaces, where Lh and Mh both are piecewise polynomial spaces with respect to the mesh 𝒯h and are both subsets of L2(Ω). 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 thatLh and Mh satisfy the following assumptions (see also [16, 17]).Assumption 1.
(i)dimMh=dimLh.
(ii) There is a constantβ>0 independent of the mesh 𝒯h such that
(26)∥ϕh∥L2(Ω)≤βsupμh∈Mh∖{0}∫Ω‍μhϕhdx∥μh∥L2(Ω),ϕh∈Lh.
(iii) The spaceMh has the approximation property:
(27)infλh∈Mh∥ϕ-λh∥L2(Ω)≤Ch|ϕ|H1(Ω),ϕ∈H1(Ω).
(iv) The spaceLh has the approximation property:
(28)infϕh∈Lh∥ϕ-ϕh∥L2(Ω)≤Ch|ϕ|H1(Ω),ϕ∈H1(Ω).
Now we define biorthogonality and quasi-biorthogonality (see [16, 18, 19]).Definition 5.
The pair of bases{μi}1≤i≤n of Mh and {φi}1≤i≤n of Lh is called  biorthogonal if the resulting Gram matrix G is diagonal. The pair of bases {μi}1≤i≤n of Mh and {φi}1≤i≤n of Lh is called  quasi-biorthogonal and the resulting Gram matrix  G  is called  quasidiagonal if  G  is of the form
(29)G=[D10RD2]orG=[D1R0D2],
where D1 and D2 are diagonal matrices and R is a sparse rectangular matrix.We recall that a Gram matrixG of two sets of basis functions {μi}1≤i≤n and {φi}1≤i≤n is the matrix G whose ijth entry is
(30)∫Ω‍μiφjdx.
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
(31)Bh={bh∣bh|T=C∏i=1d+1‍λiT,T∈𝒯h}
be the space of bubble functions, where C is a constant and {λiT}i=1d+1 is the set of barycentric coordinates on T. Let Lh=Sh⊕Bh. Explicit construction of basis functions of Lh and Mh on a reference element is provided below.The finite element problem is to find(uh,σh,ϕh)∈Vh×[Lh]d×[Mh]d such that
(32)a((uh,σh),(vh,τh))+b((vh,τh),ϕh)=ℓ(vh),annnnnnm((vhτh)nϕh)(vh,τh)∈Vh×[Lh]d,b((uh,σh),ψh)=0,ψh∈[Mh]d.
We note that the standard finite element space Sh is enriched with elementwise defined bubble functions to obtain the space Lh, which is done to obtain the coercivity of the bilinear form a(·,·) on the kernel space Kh defined as (see Lemma 7)
(33)Kh={(vh,τh)∈Vh×[Lh]d:b((vh,τh),ψh)=0,ψh∈[Mh]d}.
To simplify the numerical analysis of the finite element problem, we introduce a quasiprojection operator: Qh:L2(Ω)→Lh, which is defined as  ψh(34)∫Ω‍Qhvμhdx=∫Ω‍vμhdx,v∈L2(Ω),μh∈Mh.
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 Qh allows us to write the weak gradient as
(35)σh=Qh(∇uh),
where operator Qh is applied to the vector ∇uh component-wise. We see that Qh is well defined due to Assumption 1(i) and (ii). Furthermore, the restriction of Qh to Lh is the identity. Hence Qh is a projection onto the space Lh. We note that Qh is not the orthogonal projection onto Lh but an oblique projection onto Lh see [21, 22]. The stability and approximation properties of Qh in L2 and H1-norm can be shown as in [16, 23]. In the following, we will use a generic constant C, which will take different values at different places but will be always independent of the mesh size h.Lemma 6.
Under Assumption1(i)–(iii)
(36)∥Qhv∥L2(Ω)≤C∥v∥L2(Ω),forv∈L2(Ω),|Qhw|H1(Ω)≤C|w|H1(Ω),forw∈H1(Ω),
and for 0<s≤1 and v∈H1+s(Ω)(37)∥v-Qhv∥L2(Ω)≤Ch1+s|v|Hs+1(Ω),∥v-Qhv∥H1(Ω)≤Chs|v|Hs+1(Ω).We need to construct another finite element spaceMh⊂L2(Ω) that satisfies Assumption 1(i)–(iii) and also leads to an efficient numerical scheme. The main goal is to choose the basis functions for Mh and Lh so that the matrix D associated with the bilinear form ∫Ω‍σh:ψhdx 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
(38)-Bu→+Dσ→=0,
where u→ and  σ→are the solution vectors for uh and  σh and D is the Gram matrix between the bases of [Mh]d and [Lh]d. Hence if the matrix D is diagonal or quasi-diagonal, we can easily eliminate the degrees of freedom corresponding to  σh and  ϕh. This then leads to a formulation involving only one unknown uh. Statically condensing out variables σh and ϕh we arrive at the variational formulation of the reduced system given by (68) (see Section 3.3). Note that the matrix D is the Gram matrix between the bases of [Mh]d and [Lh]d.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 forLh and Mh so that they form a quasi-biorthogonal system. The construction is shown on the reference triangle T^={(x,y):0≤x,0≤y,x+y≤1} in the two-dimensional case and on the reference tetrahedron T^={(x,y,z):0≤x,0≤y,0≤z,x+y+z≤1} in the three-dimensional case.We now start with the reference triangle in the two-dimensional case. Let{φ^1,…,φ^4} defined as
(39)φ^1=1-x-y,φ^2=x,φ^3=y,φ^4=27xy(1-x-y)
be the local basis functions of Lh on the reference triangle associated with three vertices (0,0), and (1,0), (0,1) and one barycenter (1/3,1/3). Note that φ^4 is the standard bubble function used, for example, in the minielement for the Stokes equations [24].Then if we define the local basis functions{μ^1,…,μ^4} of Mh as
(40)μ^1=349-4x-4y-1403xy(1-x-y),μ^2=-29+4x-1403xy(1-x-y),μ^3=-29+4y-1403xy(1-x-y),μ^4=27xy(1-x-y),
the two sets of basis functions {φ^1,…,φ^4} and {μ^1,…,μ^4} form a quasi-biorthogonal system on the reference triangle. These four basis functions {μ^1,…,μ^4} of Mh 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{φ^1,…,φ^5} of Lh defined as
(41)φ^1=1-x-y-z,φ^2=x,φ^3=y,φ^4=z,φ^5=256xyz(1-x-y-z)
and the local basis functions {μ^1,…,μ^5} of Mh defined as
(42)μ^1=40192-5x-5y-5z-693023xyz(1-x-y-z),μ^2=-5992+5x-693023xyz(1-x-y-z),μ^3=-5992+5y-693023xyz(1-x-y-z),μ^4=-5992+5z-693023xyz(1-x-y-z),μ^5=256xyz(1-x-y-z)
are quasi-biorthogonal. These five basis functions {φ^1,…,φ^5} of Lh or {μ^1,…,μ^5} of Mh 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 Lh, the finite element basis functions for Mh are constructed in such a way that the Gram matrix on the reference element T^ is given by
(43)DT^=[16000016000016034034034081560]
for the two-dimensional case and
(44)DT^=[124000001240000012400000124043154315431543154096155925]
for the three-dimensional case. After ordering the degrees of freedom in σh and ϕh 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 D is quasi-diagonal.Thus in the two-dimensional case, the local basis functionsμ^1,μ^2, and μ^3 are associated with the vertices of the reference triangle, and the function μ^4 is associated with the barycenter of the triangle and defined elementwise. That means μ^4 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 φi of Lh is the same as that of the global basis function μi of Mh for i=1,…,n. However, the global basis functions μ1,…,μn of Mh are not continuous on interelement boundaries.Now we show thatMh satisfies Assumption 1(i)–(iii). As the first assumption is satisfied by construction, we consider the second assumption. Let ϕh=∑k=1n‍akϕk∈Sh and set μh=∑k=1n‍akμk∈Mh. By using the quasibiorthogonality relation (1) and the quasiuniformity assumption, we get
(45)∫Ω‍ϕhμhdx=∑i,j=1n‍aiaj∫Ω‍ϕiμjdx=∑i=1n‍ai2ci≥C∑i=1n‍ai2hid≥C∥ϕh∥0,Ω2,
where hi denotes the meshsize at ith vertex. Taking into account the fact that ∥ϕh∥0,Ω2≡∥μh∥0,Ω2≡∑i=1n‍ai2hid, we find that Assumption 1(ii) is satisfied. Since the sum of the local basis functions of Mh is one, Assumption 1(iii) can be proved as in [23, 25]. Assumption 1(iv) follows by standard arguments.There are(d+2) local basis functions of Lh, where (d+1) of them are associated with vertices and one is associated with the element barycenter. Then the local basis functions of Mh are also divided into two parts: (d+2) basis functions associated with the vertices and 1 basis function associated with the element. Let μT be the local basis function associated with the element barycenter of T∈𝒯h for Mh. Let Yh be the space of piecewise constant functions with respect to 𝒯h, and Ph:L2(Ω)→Yh is defined as
(46)Phw∣T=∫T‍wμTdx∫T‍μTdx.
Then Ph is a projection operator onto Yh. Moreover,
(47)∇vh=PhQh∇vh,
and this implies that
(48)∥∇vh∥0,Ω=∥PhQh∇vh∥0,Ω≤∥Qh∇vh∥0,Ω.
This yields the following coercivity result.Lemma 7.
The bilinear forma(·,·) is coercive on the kernel space Kh. In fact,
(49)a((uh,σh),(uh,σh))≥12(∥σh∥0,Ω2+∥∇uh∥0,Ω2),((uh,σh),(uh,σh))≥12(∥σh∥0,Ω2)(uh,σh)∈Kh.Proof.
Sinceσh=Qh∇uh on the kernel space Kh, we have
(50)a((uh,σh),(uh,σh))=12(∥σh∥0,Ω2+∥Qh∇uh∥0,Ω2).
The proof follows from the result
(51)∥∇uh∥0,Ω=∥PhQh∇uh∥0,Ω≤∥Qh∇uh∥0,Ω.Remark 8.
Using a quadrilateral or hexahedral meshes, we cannot have∇vh=PhQh∇vh, 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 Assumption1(ii).Lemma 9.
There exists a constantβ~>0 such that
(52)sup(vh,τh)∈Vh×Rhb((vh,τh),ϕh)∥vh∥1,Ω2+∥τh∥0,Ω2≥β~∥ϕh∥0,Ω,ϕh∈Rh.Proof.
The proof follows exactly as in the continuous setting. We usevh=0 in the expression
(53)sup(vh,τh)∈Vh×Rhb((vh,τh),ϕh)∥vh∥1,Ω2+∥τh∥0,Ω2
to get
(54)sup(vh,τh)∈Vh×Rhb((vh,τh),ϕh)∥vh∥1,Ω2+∥τh∥0,Ω2≥supτh∈Rh∫Ω‍τh·ϕhdx∥τh∥0,Ω2.
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 (uh,σh,ϕh)∈Vh×[Lh]d×[Mh]d, and
(55)∥uh∥1,Ω+∥σh∥0,Ω+∥ϕh∥0,Ω≤∥ℓ∥0,Ω.
Moreover, the following estimate holds for the solution (u,σ,ϕ)∈V×R×R of (32)
(56)∥u-uh∥1,Ω+∥σ-σh∥0,Ω+∥ϕ-ϕh∥0,Ω≤C(infvh∈Vh∥u-uh∥1,Ω+infτh∈[Lh]d∥σ-τh∥0,Ω+infψh∈[Mh]d∥ϕ-ψh∥0,Ω).Owing to Assumption1(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 A(·,·) (18) is coercive on the whole product space V×R. Therefore, we can use the standard finite element space Sh to discretize the gradient of the solution for the saddle-point problem (18) so that Lh=Sh. Thus our discrete saddle-point formulation of (18) is to find (uh,σh,ϕh)∈Vh×[Lh]d×[Mh]d so that
(57)A((uh,σh),(vh,τh))+B((vh,τh),ϕh)=ℓ(vh),((uh,σh),(vh,τh))+cccB(vh,τh)∈Vh×[Lh]d,B((uh,σh),ψh)=0,ψh∈[Mh]d.
We can now verify the three conditions of well-posedness for the discrete saddle point problem. The continuity of A(·,·), B(·,·), and ℓ(·) follows as in the continuous setting, and coercivity is immediate from Lemma 4 as Vh×[Lh]d⊂V×R. Explicit construction of basis functions of Mh on a reference element will be given below.Assumption1(ii) guarantees that the bilinear form b(·,·) 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 (uh,σh,ϕh)∈Vh×[Lh]d×[Mh]d and
(58)∥uh∥1,Ω+∥σh∥0,Ω+∥ϕh∥0,Ω≤∥ℓ∥0,Ω.
Moreover, the following estimate holds for the solution (u,σ,ϕ)∈V×R×R of (9)
(59)∥u-uh∥1,Ω+∥σ-σh∥0,Ω+∥ϕ-ϕh∥0,Ω≤C(infvh∈Vh∥u-uh∥1,Ω+infτh∈[Lh]d∥σ-τh∥0,Ω+infψh∈[Mh]d∥ϕ-ψh∥0,Ω).The optimal convergence of the discrete solution is then guaranteed by the standard approximation property ofVh and Lh 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 conditiondimMh=dimLh 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 T^:={(x,y):0<x,0<y,x+y<1}, we have
(60)μ^1:=3-4x-4y,μ^2:=4x-1,μ^3:=4y-1,
where the basis functions μ^1, μ^2 and μ^3 are associated with three vertices (0,0), (1,0), and (0,1) of the reference triangle. For the reference tetrahedron T^:={(x,y,z):0<x,0<y,0<z,x+y+z<1}, we have
(61)μ^1:=4-5x-5y-5z,μ^2:=5x-1,μ^3:=5y-1,μ^4:=5z-1,
where the basis functions μ^1, μ^2, μ^3, and μ^4, 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 ofSh. These basis functions also satisfy Assumption 1(i)–(iv) (see [17, 18]).


### 3.3. Static Condensation

Here we want to eliminateσh and Lagrange multiplier  ϕh from the saddle-point problem (57). We make use of the operator Qh defined previously. Using the biorthogonality relation between the basis functions of Lh and Mh, the action of operator Qh on a function v∈L2(Ω) can be written as
(62)Qhv=∑i=1n‍∫Ω‍μivdxciφi,
which tells that the operator Qh is local in the sense to be given below (see also [27]).Using the property of operatorQh, we can eliminate the degrees of freedom corresponding to  σh and  ϕh so that our problem is to find uh∈Sh such that
(63)J(uh)=minvh∈VhJ(vh),
where
(64)J(vh)=∥∇(Qh(∇vh))∥0,Ω2-2ℓ(vh)
for the nonstabilized formulation (9) and
(65)J(vh)=∥∇(Qh(∇vh))∥0,Ω2+∥Qh(∇vh)-∇vh∥0,Ω2-2ℓ(vh)
for the stabilized formulation (18).Let the bilinear formA~(·,·) be defined as
(66)A~(uh,vh)=∫Ω‍σh:τhdx
for the non-stabilized formulation (9) and
(67)A~(uh,vh)=∫Ω‍σh:τhdx+∫Ω‍(σh-∇uh)·(τh-∇vh)dx
for the stabilized formulation (18) with  σh=Qh(∇uh) and  τh=Qh(∇vh). Due to the biorthogonality or quasi-biorthogonality condition, the action of Qh is computed by inverting a diagonal matrix.Since the bilinear formA~(·,·) is symmetric, the minimization problem (63) is equivalent to the variational problem of finding uh∈Sh such that [3, 15]
(68)A~(uh,vh)=ℓ(vh),vh∈Vh.The coercivity ofA(·,·) and the stability result of Lemma 9 immediately yield the coercivity of the bilinear form A~(·,·) on the space Vh with the norm ∥uh∥0,Ω2+∥∇Qh(∇uh)∥0,Ω2 for uh∈Vh. That is,
(69)A~(uh,uh)≥C(∥uh∥0,Ω2+∥∇Qh(∇uh)∥0,Ω2),uh∈Vh.


## 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
(31)Bh={bh∣bh|T=C∏i=1d+1‍λiT,T∈𝒯h}
be the space of bubble functions, where C is a constant and {λiT}i=1d+1 is the set of barycentric coordinates on T. Let Lh=Sh⊕Bh. Explicit construction of basis functions of Lh and Mh on a reference element is provided below.The finite element problem is to find(uh,σh,ϕh)∈Vh×[Lh]d×[Mh]d such that
(32)a((uh,σh),(vh,τh))+b((vh,τh),ϕh)=ℓ(vh),annnnnnm((vhτh)nϕh)(vh,τh)∈Vh×[Lh]d,b((uh,σh),ψh)=0,ψh∈[Mh]d.
We note that the standard finite element space Sh is enriched with elementwise defined bubble functions to obtain the space Lh, which is done to obtain the coercivity of the bilinear form a(·,·) on the kernel space Kh defined as (see Lemma 7)
(33)Kh={(vh,τh)∈Vh×[Lh]d:b((vh,τh),ψh)=0,ψh∈[Mh]d}.
To simplify the numerical analysis of the finite element problem, we introduce a quasiprojection operator: Qh:L2(Ω)→Lh, which is defined as  ψh(34)∫Ω‍Qhvμhdx=∫Ω‍vμhdx,v∈L2(Ω),μh∈Mh.
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 Qh allows us to write the weak gradient as
(35)σh=Qh(∇uh),
where operator Qh is applied to the vector ∇uh component-wise. We see that Qh is well defined due to Assumption 1(i) and (ii). Furthermore, the restriction of Qh to Lh is the identity. Hence Qh is a projection onto the space Lh. We note that Qh is not the orthogonal projection onto Lh but an oblique projection onto Lh see [21, 22]. The stability and approximation properties of Qh in L2 and H1-norm can be shown as in [16, 23]. In the following, we will use a generic constant C, which will take different values at different places but will be always independent of the mesh size h.Lemma 6.
Under Assumption1(i)–(iii)
(36)∥Qhv∥L2(Ω)≤C∥v∥L2(Ω),forv∈L2(Ω),|Qhw|H1(Ω)≤C|w|H1(Ω),forw∈H1(Ω),
and for 0<s≤1 and v∈H1+s(Ω)(37)∥v-Qhv∥L2(Ω)≤Ch1+s|v|Hs+1(Ω),∥v-Qhv∥H1(Ω)≤Chs|v|Hs+1(Ω).We need to construct another finite element spaceMh⊂L2(Ω) that satisfies Assumption 1(i)–(iii) and also leads to an efficient numerical scheme. The main goal is to choose the basis functions for Mh and Lh so that the matrix D associated with the bilinear form ∫Ω‍σh:ψhdx 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
(38)-Bu→+Dσ→=0,
where u→ and  σ→are the solution vectors for uh and  σh and D is the Gram matrix between the bases of [Mh]d and [Lh]d. Hence if the matrix D is diagonal or quasi-diagonal, we can easily eliminate the degrees of freedom corresponding to  σh and  ϕh. This then leads to a formulation involving only one unknown uh. Statically condensing out variables σh and ϕh we arrive at the variational formulation of the reduced system given by (68) (see Section 3.3). Note that the matrix D is the Gram matrix between the bases of [Mh]d and [Lh]d.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 forLh and Mh so that they form a quasi-biorthogonal system. The construction is shown on the reference triangle T^={(x,y):0≤x,0≤y,x+y≤1} in the two-dimensional case and on the reference tetrahedron T^={(x,y,z):0≤x,0≤y,0≤z,x+y+z≤1} in the three-dimensional case.We now start with the reference triangle in the two-dimensional case. Let{φ^1,…,φ^4} defined as
(39)φ^1=1-x-y,φ^2=x,φ^3=y,φ^4=27xy(1-x-y)
be the local basis functions of Lh on the reference triangle associated with three vertices (0,0), and (1,0), (0,1) and one barycenter (1/3,1/3). Note that φ^4 is the standard bubble function used, for example, in the minielement for the Stokes equations [24].Then if we define the local basis functions{μ^1,…,μ^4} of Mh as
(40)μ^1=349-4x-4y-1403xy(1-x-y),μ^2=-29+4x-1403xy(1-x-y),μ^3=-29+4y-1403xy(1-x-y),μ^4=27xy(1-x-y),
the two sets of basis functions {φ^1,…,φ^4} and {μ^1,…,μ^4} form a quasi-biorthogonal system on the reference triangle. These four basis functions {μ^1,…,μ^4} of Mh 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{φ^1,…,φ^5} of Lh defined as
(41)φ^1=1-x-y-z,φ^2=x,φ^3=y,φ^4=z,φ^5=256xyz(1-x-y-z)
and the local basis functions {μ^1,…,μ^5} of Mh defined as
(42)μ^1=40192-5x-5y-5z-693023xyz(1-x-y-z),μ^2=-5992+5x-693023xyz(1-x-y-z),μ^3=-5992+5y-693023xyz(1-x-y-z),μ^4=-5992+5z-693023xyz(1-x-y-z),μ^5=256xyz(1-x-y-z)
are quasi-biorthogonal. These five basis functions {φ^1,…,φ^5} of Lh or {μ^1,…,μ^5} of Mh 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 Lh, the finite element basis functions for Mh are constructed in such a way that the Gram matrix on the reference element T^ is given by
(43)DT^=[16000016000016034034034081560]
for the two-dimensional case and
(44)DT^=[124000001240000012400000124043154315431543154096155925]
for the three-dimensional case. After ordering the degrees of freedom in σh and ϕh 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 D is quasi-diagonal.Thus in the two-dimensional case, the local basis functionsμ^1,μ^2, and μ^3 are associated with the vertices of the reference triangle, and the function μ^4 is associated with the barycenter of the triangle and defined elementwise. That means μ^4 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 φi of Lh is the same as that of the global basis function μi of Mh for i=1,…,n. However, the global basis functions μ1,…,μn of Mh are not continuous on interelement boundaries.Now we show thatMh satisfies Assumption 1(i)–(iii). As the first assumption is satisfied by construction, we consider the second assumption. Let ϕh=∑k=1n‍akϕk∈Sh and set μh=∑k=1n‍akμk∈Mh. By using the quasibiorthogonality relation (1) and the quasiuniformity assumption, we get
(45)∫Ω‍ϕhμhdx=∑i,j=1n‍aiaj∫Ω‍ϕiμjdx=∑i=1n‍ai2ci≥C∑i=1n‍ai2hid≥C∥ϕh∥0,Ω2,
where hi denotes the meshsize at ith vertex. Taking into account the fact that ∥ϕh∥0,Ω2≡∥μh∥0,Ω2≡∑i=1n‍ai2hid, we find that Assumption 1(ii) is satisfied. Since the sum of the local basis functions of Mh is one, Assumption 1(iii) can be proved as in [23, 25]. Assumption 1(iv) follows by standard arguments.There are(d+2) local basis functions of Lh, where (d+1) of them are associated with vertices and one is associated with the element barycenter. Then the local basis functions of Mh are also divided into two parts: (d+2) basis functions associated with the vertices and 1 basis function associated with the element. Let μT be the local basis function associated with the element barycenter of T∈𝒯h for Mh. Let Yh be the space of piecewise constant functions with respect to 𝒯h, and Ph:L2(Ω)→Yh is defined as
(46)Phw∣T=∫T‍wμTdx∫T‍μTdx.
Then Ph is a projection operator onto Yh. Moreover,
(47)∇vh=PhQh∇vh,
and this implies that
(48)∥∇vh∥0,Ω=∥PhQh∇vh∥0,Ω≤∥Qh∇vh∥0,Ω.
This yields the following coercivity result.Lemma 7.
The bilinear forma(·,·) is coercive on the kernel space Kh. In fact,
(49)a((uh,σh),(uh,σh))≥12(∥σh∥0,Ω2+∥∇uh∥0,Ω2),((uh,σh),(uh,σh))≥12(∥σh∥0,Ω2)(uh,σh)∈Kh.Proof.
Sinceσh=Qh∇uh on the kernel space Kh, we have
(50)a((uh,σh),(uh,σh))=12(∥σh∥0,Ω2+∥Qh∇uh∥0,Ω2).
The proof follows from the result
(51)∥∇uh∥0,Ω=∥PhQh∇uh∥0,Ω≤∥Qh∇uh∥0,Ω.Remark 8.
Using a quadrilateral or hexahedral meshes, we cannot have∇vh=PhQh∇vh, 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 Assumption1(ii).Lemma 9.
There exists a constantβ~>0 such that
(52)sup(vh,τh)∈Vh×Rhb((vh,τh),ϕh)∥vh∥1,Ω2+∥τh∥0,Ω2≥β~∥ϕh∥0,Ω,ϕh∈Rh.Proof.
The proof follows exactly as in the continuous setting. We usevh=0 in the expression
(53)sup(vh,τh)∈Vh×Rhb((vh,τh),ϕh)∥vh∥1,Ω2+∥τh∥0,Ω2
to get
(54)sup(vh,τh)∈Vh×Rhb((vh,τh),ϕh)∥vh∥1,Ω2+∥τh∥0,Ω2≥supτh∈Rh∫Ω‍τh·ϕhdx∥τh∥0,Ω2.
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 (uh,σh,ϕh)∈Vh×[Lh]d×[Mh]d, and
(55)∥uh∥1,Ω+∥σh∥0,Ω+∥ϕh∥0,Ω≤∥ℓ∥0,Ω.
Moreover, the following estimate holds for the solution (u,σ,ϕ)∈V×R×R of (32)
(56)∥u-uh∥1,Ω+∥σ-σh∥0,Ω+∥ϕ-ϕh∥0,Ω≤C(infvh∈Vh∥u-uh∥1,Ω+infτh∈[Lh]d∥σ-τh∥0,Ω+infψh∈[Mh]d∥ϕ-ψh∥0,Ω).Owing to Assumption1(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 A(·,·) (18) is coercive on the whole product space V×R. Therefore, we can use the standard finite element space Sh to discretize the gradient of the solution for the saddle-point problem (18) so that Lh=Sh. Thus our discrete saddle-point formulation of (18) is to find (uh,σh,ϕh)∈Vh×[Lh]d×[Mh]d so that
(57)A((uh,σh),(vh,τh))+B((vh,τh),ϕh)=ℓ(vh),((uh,σh),(vh,τh))+cccB(vh,τh)∈Vh×[Lh]d,B((uh,σh),ψh)=0,ψh∈[Mh]d.
We can now verify the three conditions of well-posedness for the discrete saddle point problem. The continuity of A(·,·), B(·,·), and ℓ(·) follows as in the continuous setting, and coercivity is immediate from Lemma 4 as Vh×[Lh]d⊂V×R. Explicit construction of basis functions of Mh on a reference element will be given below.Assumption1(ii) guarantees that the bilinear form b(·,·) 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 (uh,σh,ϕh)∈Vh×[Lh]d×[Mh]d and
(58)∥uh∥1,Ω+∥σh∥0,Ω+∥ϕh∥0,Ω≤∥ℓ∥0,Ω.
Moreover, the following estimate holds for the solution (u,σ,ϕ)∈V×R×R of (9)
(59)∥u-uh∥1,Ω+∥σ-σh∥0,Ω+∥ϕ-ϕh∥0,Ω≤C(infvh∈Vh∥u-uh∥1,Ω+infτh∈[Lh]d∥σ-τh∥0,Ω+infψh∈[Mh]d∥ϕ-ψh∥0,Ω).The optimal convergence of the discrete solution is then guaranteed by the standard approximation property ofVh and Lh 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 conditiondimMh=dimLh 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 T^:={(x,y):0<x,0<y,x+y<1}, we have
(60)μ^1:=3-4x-4y,μ^2:=4x-1,μ^3:=4y-1,
where the basis functions μ^1, μ^2 and μ^3 are associated with three vertices (0,0), (1,0), and (0,1) of the reference triangle. For the reference tetrahedron T^:={(x,y,z):0<x,0<y,0<z,x+y+z<1}, we have
(61)μ^1:=4-5x-5y-5z,μ^2:=5x-1,μ^3:=5y-1,μ^4:=5z-1,
where the basis functions μ^1, μ^2, μ^3, and μ^4, 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 ofSh. These basis functions also satisfy Assumption 1(i)–(iv) (see [17, 18]).


## 3.3. Static Condensation

Here we want to eliminateσh and Lagrange multiplier  ϕh from the saddle-point problem (57). We make use of the operator Qh defined previously. Using the biorthogonality relation between the basis functions of Lh and Mh, the action of operator Qh on a function v∈L2(Ω) can be written as
(62)Qhv=∑i=1n‍∫Ω‍μivdxciφi,
which tells that the operator Qh is local in the sense to be given below (see also [27]).Using the property of operatorQh, we can eliminate the degrees of freedom corresponding to  σh and  ϕh so that our problem is to find uh∈Sh such that
(63)J(uh)=minvh∈VhJ(vh),
where
(64)J(vh)=∥∇(Qh(∇vh))∥0,Ω2-2ℓ(vh)
for the nonstabilized formulation (9) and
(65)J(vh)=∥∇(Qh(∇vh))∥0,Ω2+∥Qh(∇vh)-∇vh∥0,Ω2-2ℓ(vh)
for the stabilized formulation (18).Let the bilinear formA~(·,·) be defined as
(66)A~(uh,vh)=∫Ω‍σh:τhdx
for the non-stabilized formulation (9) and
(67)A~(uh,vh)=∫Ω‍σh:τhdx+∫Ω‍(σh-∇uh)·(τh-∇vh)dx
for the stabilized formulation (18) with  σh=Qh(∇uh) and  τh=Qh(∇vh). Due to the biorthogonality or quasi-biorthogonality condition, the action of Qh is computed by inverting a diagonal matrix.Since the bilinear formA~(·,·) is symmetric, the minimization problem (63) is equivalent to the variational problem of finding uh∈Sh such that [3, 15]
(68)A~(uh,vh)=ℓ(vh),vh∈Vh.The coercivity ofA(·,·) and the stability result of Lemma 9 immediately yield the coercivity of the bilinear form A~(·,·) on the space Vh with the norm ∥uh∥0,Ω2+∥∇Qh(∇uh)∥0,Ω2 for uh∈Vh. That is,
(69)A~(uh,uh)≥C(∥uh∥0,Ω2+∥∇Qh(∇uh)∥0,Ω2),uh∈Vh.


## 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.



---

*Source: 189045-2013-04-11.xml*