A Preconditioner for FETI-DP Method of Stokes Problem with Mortar-Type CR(ℎ)−𝑃₀(ℎ) Discretization

Abstract

FETI-DP method is developed for incompressible Stokes problem discretized by mortar-type CR(ℎ)−𝑃₀(ℎ) mixed elements. A preconditioner is proposed, and it is proved that the condition number of the preconditioned FETI-DP operator is bounded by (1+log(𝐻/ℎ))², where 𝐻 and ℎ are mesh sizes. Finally, the numerical tests are presented to verify the theoretical results.

1. Introduction

FETI-DP method is a nonoverlapping domain decomposition method which has been implemented for large scale engineering applications. It has been developed for elliptic problems with heterogeneous coefficients, Stokes problems, and compressible elasticity problems (cf. [1–4]). Recently, extensions to irregular subdomains and inexact subdomain solvers have been done in [5–7]. Li [4] developed a FETI-DP method for the Stokes problem by adding primal continuity constraints to accelerate the convergence of the method. Li and Widlund [8] developed BDDC methods for incompressible Stokes equations which established a close connection between the BDDC and FETI-DP methods for the Stokes case. Kim et al. [9] developed a BDDC method discretized by mortar finite element methods for geometrically nonconforming partitions, and the BDDC preconditioner was shown to be closely related to the Neumann-Dirichlet version of the FETI-DP method.

Nonconforming discretizations are important for multiphysics simulations, contact-impact problems, generations of meshes and partitions aligned with jumps in diffusion coefficients, ℎ𝑝-adaptive methods, and special discretizations in the neighborhood of singularities. In this paper, we use the inf-sup stable CR(ℎ)−𝑃₀(ℎ) mixed elements in each subdomain. In the event of applying the mortar condition, it is necessary to know functions on the interface. However, the degrees of freedom (d.o.f) of Crouzeix-Raviart element are associated with the edge midpoints and functions on the interface depend on the nodal values corresponding to interface nodes and some subdomain interior nodes lying closest to the interface. We adopt nonstandard mortar condition introduced in [10] which is only associated with nodal values on the interface.

The inf-sup constant for mortar element space is crucial in the analysis of the approximation order. If the constant is independent of mesh size and subdomain size, then the optimal order of approximation follows independently of the number of subdomains and mesh size as in the case of elliptic problems.

The fundamental idea of this paper is the same as [4], with additional technical complications caused by nonmatching grids across subdomain interfaces and the Crouzeix-Raviart element discretization for velocity space. We propose a preconditioner and analyze the condition number bounded by (1+log(𝐻/ℎ))², where 𝐻 and ℎ are mesh sizes.

The paper is organized as follows. In Section 2, a domain decomposition procedure is described for solving incompressible Stokes problem. In Section 3, the matrix form of the mortar condition is introduced. In Section 4, we derive an FETI-DP operator for the Stokes problem and propose a preconditioner for the FETI-DP operator. Section 5 gives several technical lemmas. Section 6 is devoted to the proof of the condition number of the preconditioned FETI-DP operator. The numerical tests are presented to verify the theoretical results in Section 7.

Throughout the paper, the following notations are used: 𝑢≍𝑣,𝑢⪯𝑣, and 𝑢≽𝑣 mean that there exists positive constants 𝑐 and 𝐶 independent of the subdomain size and mesh size, such that 𝑐𝑢≤𝑣≤𝐶𝑢,𝑢≤𝐶𝑣,𝑢≥𝑐𝑣,respectively. (1)

2. A Domain Decomposition Procedure for Stokes Problems

We consider the following two-dimensional incompressible Stokes problem: −▵𝐮+∇𝑝=𝐟,inΩ,−∇⋅𝐮=0,inΩ,𝐮=𝟎,on𝜕Ω, (2) where Ω is a bounded polygonal domain in 𝑅²,𝐮 represents the velocity of fluid, 𝑝 is the pressure, and 𝐟 is the external force.

We first decompose Ω into a nonoverlapping subdomain partition {Ω_𝑖}^𝑁_𝑖=1. We consider the subdomains to be geometrically conforming. With each subdomain Ω_𝑖, we assume that the triangulation 𝑇_ℎ(Ω_𝑖) is quasiuniform. Let ℎ_𝑖=max_{𝜏∈𝑇ℎ(Ω𝑖)}diam𝜏 be the mesh size parameter. Let Γ_𝑖𝑗=𝜕Ω_𝑖∩𝜕Ω_𝑗 denote the interface between the two neighboring subdomains Ω_𝑖 and Ω_𝑗, and Γ=∪^𝑁_𝑖=1𝜕Ω_𝑖⧵𝜕Ω be the interface skeleton. We denote the sets of the Crouzeix-Raviart element nodes and the P1 conforming element nodes that are contained in Ω,𝜕Ω,Ω_𝑖,𝜕Ω_𝑖, and Γ_𝑖𝑗 by Ω^CR_ℎ, 𝜕Ω^CR_ℎ, Ω^CR_𝑖,ℎ, 𝜕Ω^CR_𝑖,ℎ, Γ^CR_𝑖𝑗,ℎ and Ω_ℎ, 𝜕Ω_ℎ, Ω_𝑖,ℎ, 𝜕Ω_𝑖,ℎ, Γ_𝑖𝑗,ℎ, respectively. The resulting triangulation can be nonmatching across the subdomain interface. We note that each interface Γ_𝑖𝑗 inherits two independent triangulations. One side of Γ_𝑖𝑗 associated with Ω_𝑖 is defined as a mortar side denoted by 𝛾_𝑚(𝑖), and the other side associated with Ω_𝑗 is defined as a nonmortar side denoted by 𝛿_𝑚(𝑗).

For each subdomain Ω_𝑖, we introduce the following spaces: 𝐻¹_𝐷Ω_𝑖=𝑣∈𝐻¹Ω_𝑖∶𝑣=0on𝜕Ω_𝑖,𝐿∩𝜕Ω₂₀Ω_𝑖=𝑝∈𝐿²Ω_𝑖∶_Ω𝑖,Π𝑝𝑑𝑥=0₀=𝑝⁰∶𝑝⁰_Ω𝑖isconstantand_Ω𝑝⁰𝑑𝑥=0. (3) Then, the variational form of the problem (2) is as follows: find (𝐮,𝑝_𝐼,𝑝⁰)∈ ∏^𝑁_𝑖=1[𝐻¹_𝐷(Ω_𝑖)]²× ∏^𝑁_𝑖=1𝐿²₀(Ω_𝑖)×Π⁰, such that ^𝑁_𝑖=1(∇𝐮,∇𝐯)_Ω𝑖−^𝑁_𝑖=1𝑝_𝐼+𝑝⁰,∇⋅𝐯_Ω𝑖=^𝑁_i=1(𝐟,𝐯)_Ω𝑖,∀𝐯∈^𝑁_𝑖=1𝐻¹_𝐷Ω_𝑖²,−^𝑁_𝑖=1∇⋅𝐮,𝑞_𝐼Ω𝑖=0,∀𝑞_𝐼∈^𝑁_𝑖=1𝐿²₀Ω_𝑖,−^𝑁_𝑖=1∇⋅𝐮,𝑞⁰_Ω𝑖=0,∀𝑞⁰∈Π⁰,𝐮|_Ω𝑖−𝐮|_Ω𝑗=0,∀Γ_𝑖𝑗⊂Γ. (4) Here, (⋅,⋅)_Ω𝑖 denotes the inner product in [𝐿²(Ω_𝑖)]^𝑛 for 𝑛=1,2.

We consider the inf-sup stable CR(ℎ_𝑖)−𝑃₀(ℎ_𝑖) elements in each subdomain Ω_𝑖. Let 𝐖^ℎ(Ω_𝑖)={𝐯_𝑖∈ [𝐿²(Ω_𝑖)]²∶𝐯_𝑖|_𝜏∈[𝑃₁(𝜏)]²,𝐯_𝑖   is  continuous  at  the  midpoints  of  the edges  of 𝜏, and 𝐯_𝑖(𝑚)=𝟎,𝑚 is  a  midpoint  on 𝜕Ω∩𝜕𝜏,∀𝜏∈𝑇_ℎ(Ω_𝑖)}, 𝑄_𝑖=𝑞_𝑖∈𝐿²₀Ω_𝑖∶𝑞_𝑖|_𝜏∈𝑃₀(𝜏),∀𝜏∈𝑇_ℎΩ_𝑖, (5) where 𝑃_𝑚(𝜏) is a set of polynomials of degree less than or equal to 𝑚 in 𝜏.

The Lagrange multiplier space corresponding to Γ_𝑖𝑗 is defined as follows: 𝐌𝛿_𝑚(𝑗)=𝝀∈𝐿²Γ_𝑖𝑗²∶𝝀|_𝜏isaconstantvector,∀𝜏∈𝑇_ℎ𝛿_𝑚(𝑗). (6) Then we take the Lagrange multiplier space 𝐌=_{𝛿𝑚(𝑗)⊂Γ}𝐌𝛿_𝑚(𝑗). (7) We introduce the nonstandard mortar condition (cf. [10]): 𝑄_𝑚𝐽_m𝐯_𝑖=𝑄_𝑚𝐯_𝑗,onΓ_𝑖𝑗=𝛾_𝑚(𝑖)=𝛿_𝑚(𝑗), (8) where 𝐯_𝑘∈𝐖^ℎ(Ω_𝑘)(𝑘=𝑖,𝑗),𝑄_𝑚∶[𝐿²(Γ_𝑖𝑗)]²→𝐌(𝛿_𝑚(𝑗)) is the 𝐿² orthogonal projection defined as (𝑄_𝑚𝐮,𝝀)_{[𝐿2(𝛿𝑚(𝑗))]}²=(𝐮,𝝀)_{[𝐿2(𝛿𝑚(𝑗))]}²,∀𝝀∈𝐌(𝛿_𝑚(𝑗)). The operator 𝐽_𝑚 is defined in the following. We can also write (8) as follows: _{𝛿𝑚(𝑗)}𝐽_𝑚𝐯_𝑖−𝐯_𝑗⋅𝝀_𝑖𝑗𝑑𝑠=0,∀𝝀_𝑖𝑗∈𝐌𝛿_𝑚(𝑗). (9)

We use 𝑒 to denote a triangle edge. Let 𝑍_ℎ(𝛾_𝑚(𝑖))=∏_{𝑒∈𝑇ℎ(𝛾𝑚(𝑖))}𝑃₁(𝑒) be the space of piecewise linear functions defined on the triangulation 𝑇_ℎ(𝛾_𝑚(𝑖)), and let 𝑇_ℎ/2(𝛾_𝑚(𝑖)) be the triangulation obtained as a result of dividing each edge of 𝑇_ℎ(𝛾_𝑚(𝑖)) into two equal segments. Let 𝑌_ℎ/2(𝛾_𝑚(𝑖)) be the conforming space of piecewise linear continuous functions on the triangulation 𝑇_ℎ/2(𝛾_𝑚(𝑖)). The midpoint, left and right endpoints of each edge 𝑒∈𝑇_ℎ(𝛾_𝑚(𝑖)) are denoted by 𝑥^𝑒_𝑚, 𝑥^𝑒_𝑙, and 𝑥^𝑒_𝑟, respectively. The length of 𝑒 is denoted by |𝑒|.

We now define the operator 𝐽_𝑚∶𝐖^ℎ(𝛾_𝑚(𝑖))=𝐖^ℎ(Ω_𝑖)|_{𝛾𝑚(𝑖)}→[𝑍_ℎ(𝛾_𝑚(𝑖))]². It is based on the definition of another interpolation operator 𝐼_𝑚∶𝐖^ℎ(𝛾_𝑚(𝑖))→[𝑌_ℎ/2(𝛾_𝑚(𝑖))]².

Definition 1. For 𝐮∈𝐖^ℎ(𝛾_𝑚(𝑖)),𝐼_𝑚𝐮∈[𝑌_ℎ/2(𝛾_𝑚(𝑖))]² is defined by the nodal values as 𝐼_𝑚𝐮(x)=𝐮(𝑥),𝑥∈𝛾^CR_{𝑚(𝑖),ℎ},𝑒_𝑟𝑒_𝑙+𝑒_𝑟𝐮𝑥^𝑒𝑙_𝑚+𝑒_𝑙𝑒_𝑙+𝑒_𝑟𝐮𝑥^𝑒𝑟_𝑚,𝑥∈𝛾_{𝑚(𝑖),ℎ},𝐮𝑥^𝑒𝑒_𝑚+𝑒_𝑒𝑒_𝑒+𝑒^′_𝑒𝐮𝑥^𝑒𝑒_𝑚−𝐮𝑥^𝑒′𝑒_𝑚,𝑥∈𝜕𝛾_{𝑚(𝑖),ℎ}. (10)

Here, 𝑒_𝑙 and 𝑒_𝑟 are the left- and the right-neighboring edges of 𝑥∈𝛾_{𝑚(𝑖),ℎ}, respectively. 𝑒_𝑒 represents a triangle edge of 𝑇_ℎ(𝛾_𝑚(𝑖)), touching 𝜕𝛾_𝑚(𝑖), and 𝑒^′_𝑒 is the corresponding neighboring edge.

Each component of the interpolation is done basically by first joining the neighboring edge midpoints using straight lines, and then simply extending the two end straight lines toward the end of the mortar 𝛾_𝑚(𝑖) (cf. Figure 1). The operator 𝐽_𝑚 can now be defined using 𝐼_𝑚.

Figure 1 

Here each component of 𝐮∈𝑊ℎ(𝛾𝑚(𝑖)) is a (discontinuous) piecewise linear function shown (using solid lines) in the lower figure. The dashed lines (in both upper and lower figures) correspond to each component of 𝐼𝑚𝐮. The solid lines in the upper figure correspond to each component of 𝐽𝑚𝐮.

Definition 2. For 𝐮∈𝐖^ℎ(𝛾_𝑚(𝑖)),𝐽_𝑚𝐮∈[𝑍_ℎ(𝛾_𝑚(𝑖))]² is a vector whose components are piecewise linear functions on the edges {𝑒} of 𝛾_𝑚(𝑖) defined by its values at the two endpoints 𝑥^𝑒_𝑙,𝑥^𝑒_𝑟∈𝛾_{𝑚(𝑖),ℎ} of each edge 𝑒. If 𝑒 is an interior edge of 𝛾_𝑚(𝑖), then 𝐽_𝑚𝐮(x)=𝐮𝑥^𝑒_𝑚+12𝐼_𝑚𝐮𝑥^𝑒_𝑙−𝐼_𝑚𝐮𝑥^𝑒_𝑟,𝑥=𝑥^𝑒_𝑙,𝐮𝑥^𝑒_𝑚+12𝐼_𝑚𝐮𝑥^𝑒_𝑟−𝐼_𝑚𝐮𝑥^𝑒_𝑙,𝑥=𝑥^𝑒_𝑟. (11)

It is easy to see that, if 𝑒 is a boundary edge of 𝛾_𝑚(𝑖), then 𝐽_𝑚𝐮(𝑥)=𝐼_𝑚𝐮(𝑥) for 𝑥=𝑥^𝑒_𝑙,𝑥^𝑒_𝑟.

We define the following spaces: 𝐖^ℎ(Ω)=^𝑁_𝑖=1𝐖^ℎΩ_𝑖,_∼𝐖^ℎ(Ω)=𝐯∈𝐖^ℎ(Ω)∶_Γ𝑖𝑗𝐯_𝑖−𝐯_𝑗𝑑𝑠=𝟎,∀Γ_𝑖𝑗⊂Γ,𝐕=𝐯∈𝐖^∼ℎ(Ω),𝐯satisfies(9,𝑄_𝐼)=^𝑁_𝑖=1𝑄_𝑖,𝑃=𝑞∈𝐿²₀(Ω)∶𝑞_Ω𝑖∈𝑄_𝑖,∀𝑖=1,…,𝑁, (12) where 𝐕 and 𝑃 are the spaces for velocity and pressure, respectively.

3. Matrix form of Mortar Conditions

For 𝐯=(𝐯^𝑡₁,…,𝐯^𝑡_𝑁)^𝑡∈𝐖^∼ℎ(Ω),𝐯_𝑖=(𝐯^𝑖_𝐼,𝐯^𝑖_Δ,𝐯^𝑖_𝑐)^𝑡, where the symbols 𝐼, Δ, and 𝑐 represent the d.o.f corresponding to the interior nodes, nodes related to ∫_Γ𝑖𝑗𝐯𝑑𝑠=𝟎, and nodes related to the remaining parts of edges, respectively. We define the spaces 𝐖_𝐼, 𝐖_Δ, and 𝐖_𝑐 which consists of 𝐯_𝐼, 𝐯_Δ, and 𝐯_𝑐, respectively.

Let 𝐵_𝑖^𝑖𝑗 be a matrix with entries 𝐵_𝑖^𝑖𝑗_𝑙𝑘=_{𝛿𝑚(𝑗)}𝝍_𝑙⋅𝝓_𝑘𝑑𝑠,∀𝑙=1,…,𝐿,∀𝑘=1,…,𝐾, (13) where {𝝍_𝑙}^𝐿_𝑙=1 is basis for 𝐌(𝛿_𝑚(𝑗)) and {𝝓_𝑘}^𝐾_𝑘=1 is basis for 𝐽_𝑚𝑊^ℎ(Ω_𝑖)|_{𝛿𝑚(𝑗)}. Let 𝐵_𝑗^𝑖𝑗 be a matrix with entries 𝐵_𝑗^𝑖𝑗_𝑙𝑘=−_{𝛿𝑚(𝑗)}𝝍_𝑙⋅𝝋_𝑘𝑑𝑠,∀𝑙=1,…,𝐿,∀𝑘=1,…,𝑃, (14) where {𝝍_𝑙}^𝐿_𝑙=1 is basis for 𝐌(𝛿_𝑚(𝑗)) and {𝝋_𝑘}^𝑃_𝑘=1 is basis for 𝑊^ℎ(Ω_𝑗)|_{𝛿𝑚(𝑗)}. Then we write (9) into a matrix form as 𝐵_𝑖^𝑖𝑗𝐯_{𝑖𝛿𝑚(𝑗)}+𝐵_𝑗^𝑖𝑗𝐯_{𝑗𝛿𝑚(𝑗)}=𝟎. (15)

Let 𝐖_𝑖=𝐖^ℎ(Ω_𝑖)|_𝜕Ω𝑖,𝑖=1,…,𝑁. The standard discrete Stokes harmonic extension operator 𝑆𝐻_𝑖∶ 𝐖_𝑖→𝐖^ℎ(Ω_𝑖) is defined by the following: given a velocity 𝐮_𝑖∈𝐖_𝑖, find 𝑆𝐻_𝑖𝐮_𝑖∈𝐖^ℎ(Ω_𝑖), 𝑝^𝑖_𝐼∈𝑄_𝑖 and 𝑝⁰∈Π⁰, such that on each subdomain Ω_𝑖, ^𝑁_𝑖=1∇𝐮_𝑖,∇𝐯_𝑖Ω𝑖−^𝑁_𝑖=1𝑝^𝑖_𝐼+𝑝⁰,∇⋅𝐯_𝑖Ω𝑖=𝟎,∀𝐯_𝑖∈𝐖^ℎΩ_𝑖,−^𝑁_𝑖=1∇⋅𝐮_𝑖,𝑞^𝑖_𝐼Ω𝑖=0,∀𝑞^𝑖_𝐼∈𝑄_𝑖,−^𝑁_𝑖=1∇⋅𝐮_𝑖,𝑞⁰_Ω𝑖=0,∀𝑞⁰∈Π⁰,𝑆𝐻_𝑖𝐮_𝑖=𝐮_𝑖,on𝜕Ω_𝑖. (16)

Define 𝐸_𝑖𝑗∶𝐌(𝛿_𝑚(𝑗))→𝐌 to be an extension operator by zero and 𝑅^𝑙_𝑖𝑗∶𝐖_𝑙→𝐖_𝑙|_{𝛿𝑚(𝑗)} for 𝑙=𝑖,𝑗 to be a restriction operator, and let 𝐵_𝑖=∑_{𝛿𝑚(𝑗)⊂Γ}𝐸_𝑖𝑗𝐵_𝑖^𝑖𝑗𝑅^𝑖_𝑖𝑗+∑_{𝛿𝑚(𝑖)⊂Γ}𝐸_𝑖𝑗𝐵_𝑖^𝑖𝑗𝑅^𝑖_𝑖𝑗. Then (15) becomes 𝐵𝐰=0, (17) where 𝐵=𝐵₁,…,𝐵_𝑁,𝐰=𝐰^𝑡₁,…,𝐰^𝑡_𝑁^𝑡,𝐰_𝑖=𝐯_{𝑖𝜕Ω𝑖},𝑖=1,…,𝑁. (18)

Let 𝐵_𝑖,𝑐 and 𝐵_𝑖,Δ be matrices that consist of the columns of 𝐵_𝑖 corresponding to the d.o.f related to ∫_Γ𝑖𝑗𝐰𝑑𝑠=0 and the d.o.f related to the remaining parts of the edges, respectively. Then, (17) can also be written into 𝐵_Δ𝐰_Δ+𝐵_𝑐𝐰_c=0, (19) where 𝐵_Δ=(𝐵_1,Δ,…,𝐵_𝑁,Δ) and 𝐵_𝑐=∑^𝑁_𝑖=1𝐵_𝑖,𝑐𝐿^𝑖_𝑐, where 𝐿^𝑖_𝑐∶𝐖_𝑐→𝐖^𝑖_𝑐 is the restriction map that gives 𝐿^𝑖_𝑐𝐰_𝑐=𝐰^𝑖_𝑐,∀𝑖=1,…,𝑁, for 𝐰_𝑐∈𝐖_𝑐. We will use any of these two expressions (17) and (19) for the sake of convenience.

4. FETI-DP Method

4.1. FETI-DP Operator

In this section, we formulate a FETI-DP operator with the mortar condition (19).

We define the following spaces: 𝐖=𝐯∈^𝑁_𝑖=1𝐖_𝑖∶_Γ𝑖𝑗𝐯_𝑖−𝐯_𝑗𝑑𝑠=0. (20)

To solve the Stokes problem efficiently and correctly, we will consider the following primal constraints: _{𝛿𝑚(𝑗)}𝐽_𝑚𝐯_𝑖−𝐯_𝑗𝑑𝑠=𝟎,∀𝛿_𝑚(𝑗)⊂Γ. (21)

Note that (21) holds by replacing 𝝀_𝑖𝑗=¹₀ or 01 in (9), because constant multipliers belong to the Lagrange multiplier space 𝐌(𝛿_𝑚(𝑗)). These primal constraints were introduced by Li (cf. [4]). The primal constraints enlarge the size of the coarse problem, so that it may lead to a fast convergence of the FETI-DP iteration.

We rewrite (21) as 𝑅^𝑡𝐵_Δ𝐰_Δ+𝐵_𝑐𝐰_𝑐=0, (22) where the matrix 𝑅 is a Boolean matrix, the number of columns of the matrix 𝑅 equals to twice the number of Γ_𝑖𝑗, and the number of rows of the matrix 𝑅 equals to the d.o.f of the space 𝐌. For 𝝀∈𝐌, at each interior nodal point of Γ_𝑖𝑗, 𝝀|_Γ𝑖𝑗 has two components corresponding to horizontal and vertical parts of velocity function. For 𝝀∈𝐌,𝑅^𝑡𝝀=0 means that for all Γ_𝑖𝑗, the sums of 𝝀|_Γ𝑖𝑗 corresponding to the horizontal and vertical parts of velocity function are zero.

Let 𝑁 be the Lagrange multiplier space corresponding to the constraints (22), and for 𝝁∈𝑁,𝝁|_Γ𝑖𝑗 has two components that correspond to the constraints for horizontal and vertical velocity. By introducing Lagrange multipliers 𝝀 and 𝝁 to enforce the constraints (19) and (22), the following is induced from the Galerkin approximation to (4) as follows: find (𝐮_𝐼,𝑝_𝐼,𝐮_Δ,𝐮_𝑐,𝑝⁰,𝝁,𝝀)∈ 𝐖_𝐼×𝑄_𝐼×𝐖_Δ×𝐖_𝑐× Π⁰×𝑁×𝐌, such that 𝐴_𝐼𝐼𝐺_𝐼𝐼𝐴_𝐼Δ𝐴_𝐼𝑐𝐺_𝐼0𝐺𝟎𝟎_𝑡𝐼𝐼𝟎𝐺^𝑡_Δ𝐼𝐺^𝑡_𝑐𝐼𝐴𝟎𝟎𝟎_Δ𝐼𝐺_Δ𝐼𝐴_ΔΔ𝐴_Δ𝑐𝐺_Δ0𝐵^𝑡_Δ𝑅𝐵^𝑡_Δ𝐴_𝑐𝐼𝐺_𝑐𝐼𝐴_𝑐Δ𝐴_𝑐𝑐𝐺_𝑐0𝐵^𝑡_𝑐𝑅𝐵^𝑡_𝑐𝐺^𝑡_𝐼0𝟎𝐺^𝑡_Δ0𝐺^𝑡_𝑐0𝟎𝟎𝟎𝟎𝟎𝑅^𝑡𝐵_Δ𝑅^𝑡𝐵_𝑐𝟎𝟎𝟎𝟎𝟎𝐵_Δ𝐵_𝑐𝐮𝟎𝟎𝟎^𝐼𝑝_𝐼𝐮_Δ𝐮_𝑐𝑝⁰𝝁𝝀=𝐟_𝐼𝟎𝐟_Δ𝐟_𝑐𝟎𝟎𝟎. (23)

Here, 𝐴_𝐼𝐼𝐴_𝐼Δ𝐴_𝐼𝑐𝐴_Δ𝐼𝐴_ΔΔ𝐴_Δ𝑐𝐴_𝑐𝐼𝐴_𝑐Δ𝐴_𝑐𝑐 is a stiffness matrix induced from ∑^𝑁_𝑖=1(∇𝐮,∇𝐯)_Ω𝑖, (𝐺^𝑡_𝐼𝐼𝐺^𝑡_Δ𝐼𝐺^𝑡_𝑐𝐼) is a matrix induced from ∑^𝑁_𝑖=1(−∇⋅𝐯,𝑝_𝐼)_Ω𝑖, and (𝐺^𝑡_𝐼0𝐺^𝑡_Δ0𝐺^𝑡_𝑐0) is a matrix induced from ∑^𝑁_𝑖=1(−∇⋅𝐯,𝑝⁰)_Ω𝑖. Since 𝑝⁰|_Ω𝑖 is constant and 𝐯_𝐼|_𝜕Ω𝑖=𝟎 for 𝐯_𝐼∈𝐖_𝐼, we have 𝐺_𝐼0=𝟎 by the divergence theorem.

Let 𝐾_ΔΔ=𝐴_𝐼𝐼𝐺_𝐼𝐼𝐴_𝐼Δ𝐺^𝑡_𝐼𝐼𝟎𝐺^𝑡_Δ𝐼𝐴_Δ𝐼𝐺_Δ𝐼𝐴_ΔΔ,𝐾_Δ𝑐=𝐴_𝐼𝑐𝐺_𝐼0𝟎𝐺^𝑡_𝑐𝐼𝐴𝟎𝟎_Δ𝑐𝐺_Δ0𝐵^𝑡_Δ𝑅,𝐾_𝑐𝑐=𝐴_𝑐𝑐𝐺_𝑐0𝐵^𝑡_𝑐𝑅𝐺^𝑡_𝑐0𝑅𝟎𝟎_𝑡𝐵_𝑐,𝟎𝟎_∼𝐵_Δ=𝟎𝟎𝐵_Δ,𝐵^∼_𝑐=𝐵_𝑐0𝟎_∼,𝐟_Δ=𝐟_𝐼𝟎𝐟_Δ,^∼𝐟_𝑐=𝐟_𝑐𝟎𝟎,𝐳_Δ=𝐮_𝐼𝑝_𝐼𝐮_Δ,𝐳_𝑐=𝐮_𝑐𝑝⁰𝝁, (24) where 𝐳_𝑐 is a primal variable, then (23) can be written as 𝐾_ΔΔ𝐾^∼_Δ𝑐𝐵^𝑡_Δ𝐾^𝑡_Δ𝑐𝐾^∼_𝑐𝑐𝐵^𝑡_𝑐∼𝐵𝐵_Δ^∼_𝑐𝟎𝐳_Δ𝐳_𝑐𝝀=^∼𝐟_Δ∼𝐟_𝑐𝟎. (25)

In Section 5, we will show that the matrix 𝐾_ΔΔ is invertible (cf. Lemma 9). After eliminating 𝐳_Δ, we obtain the following equations for 𝐳_𝑐 and 𝝀: −𝐹_𝑐𝑐𝐹_𝑐𝑙𝐹^𝑡_𝑐𝑙𝐹_𝑙𝑙𝐳_𝑐𝝀=−𝐝_𝑐𝐝_𝑙, (26) where 𝐹_𝑙𝑙=𝐵^∼_Δ𝐾^∼−1_ΔΔ𝐵^𝑡_Δ,𝐹_𝑐𝑙=𝐾^𝑡_Δ𝑐𝐾^∼−1_ΔΔ𝐵^𝑡_Δ−𝐵^∼𝑡_𝑐,𝐹_𝑐𝑐=𝐾_𝑐𝑐−𝐾^𝑡_Δ𝑐𝐾⁻¹_ΔΔ𝐾_Δ𝑐,𝐝_𝑙=𝐵^∼_Δ𝐾^∼−1_ΔΔ𝐟_Δ,𝐝_𝑐=^∼𝐟_𝑐−𝐾^𝑡_Δ𝑐𝐾^∼−1_ΔΔ𝐟_Δ. (27) The matrix 𝐹_𝑐𝑐 which is a coarse problem in the FETI-DP method is invertible (cf. Lemma 10). Hence by eliminating 𝐳_𝑐, we obtain the following equation for 𝝀∶𝐹_DP𝝀=𝐝_𝑙−𝐹^𝑡_𝑐𝑙𝐹⁻¹_𝑐𝑐𝐝_𝑐, (28) where 𝐹_DP=𝐹_𝑙𝑙+𝐹^𝑡_𝑐𝑙𝐹⁻¹_𝑐𝑐𝐹_𝑐𝑙, and we call it the FETI-DP operator. Since we add the primal constraints to the coarse problem, 𝝀 is not uniquely determined in the space 𝐌. We define a subspace of 𝐌 as follows: 𝐌^∼=𝝀∈𝐌∶𝑅^𝑡𝝀=𝟎. (29) The matrix 𝐹_DP is symmetric and positive definite (s.p.d) on 𝐌^∼, and the solution 𝝀 of (28) is uniquely determined in 𝐌^∼.

4.2. Preconditioner

For 𝐰_𝑖∈𝐖_𝑖, we define 𝑆_𝑖𝐰_𝑖 by 𝐴^𝑖_𝐼𝐼𝐺^𝑖_𝐼𝐼𝐴^𝑖_𝐼Δ𝐴^𝑖_𝐼𝑐𝐺^𝑖_𝐼𝐼^𝑡𝟎𝐺^𝑖_Δ𝐼^𝑡𝐺^𝑖_𝑐𝐼^𝑡𝐴^𝑖_Δ𝐼𝐺^𝑖_Δ𝐼𝐴^𝑖_ΔΔ𝐴^𝑖_Δ𝑐𝐴^𝑖_𝑐𝐼𝐺^𝑖_𝑐𝐼𝐴^𝑖_𝑐Δ𝐴^𝑖_𝑐𝑐𝐮^𝑖_𝐼𝑝^𝑖_𝐼𝐰^𝑖_Δ𝐰^𝑖_𝑐=𝟎𝟎𝑆_𝑖𝐰^𝑖_Δ𝐰^𝑖_𝑐, (30) where the subscript 𝑖 for a matrix denotes the submatrix corresponding to subdomain Ω_𝑖. Since the upper left 2×2 matrix 𝐴^𝑖_𝐼𝐼𝐺^𝑖_𝐼𝐼𝐺^𝑖_𝐼𝐼^𝑡𝟎 (31) represents the local Stokes problem with a Dirichlet boundary condition, it is invertible, so that the Schur complementary 𝑆_𝑖 is well defined. We then assemble local Schur complementary matrices and define 𝑆=diag𝑆₁,…,𝑆_𝑁 (32) and it can be seen easily that 𝑆 is an s.p.d operator on 𝐖. Hence, we define ‖𝐰‖_𝐖=^𝑁_𝑖=1⟨𝑆_𝑖𝐰_𝑖,𝐰_𝑖⟩^1/2 (33) as a norm for 𝐰∈𝐖. Here, ⟨⋅,⋅⟩ denotes the 𝐿²-inner product of vectors.

We define the following spaces: 𝐖⁰_𝑖𝑗=𝐰_𝑖𝑗∈𝐖^ℎΩ_{𝑖𝛿𝑚(𝑗)}∶_{𝛿𝑚(𝑗)}𝐰_𝑖𝑗,𝐖𝑑𝑠=𝟎₀=_{𝑖𝛿𝑚(𝑗)⊂Γ}𝐖⁰_𝑖𝑗. (34)

For a function 𝐰_𝑖𝑗∈𝐖⁰_𝑖𝑗, let 𝐰^∼_𝑖𝑗 be the zero extension of 𝐰_𝑖𝑗 into 𝐖_𝑖. For 𝐰∈𝐖⁰ we define an extension 𝐰^∼∈𝐖 by 𝐰^∼=𝐰^∼₁,…,𝐰^∼_𝑁with𝐰^∼_𝑖=_{𝛿𝑚(𝑗)}^∼_⊂Γ𝐰_𝑖𝑗,∀𝑖=1,…,𝑁, (35) and define a norm on 𝐖⁰ by ‖𝐰‖_𝐖⁰=𝐰^∼_𝐖. (36) We introduce the following subspaces with the norms induced from the spaces 𝐖 and 𝐖⁰: 𝐖^∼=𝐰∈𝐖∶_{𝛿𝑚(𝑗)}𝐽_𝑚𝐰_𝑖−𝐰_𝑗𝑑𝑠=𝟎,∀𝛿_𝑚(𝑗)⊂Γ_∼,𝐖_Δ=𝐰∈𝐖^∼∶_𝜕Ω𝑖𝐰_𝑖⋅𝐧_𝑖=0,𝑖=1,…,𝑁_∼,𝐖⁰=𝐰∈𝐖⁰∶𝐰^∼∈𝐖^∼, (37)(38)(39) where 𝐧_𝑖 denotes the outward unit normal vector on the subdomain boundary 𝜕Ω_𝑖.

Recall the definition of 𝐌^∼ in (29), and let ⟨⋅,⋅⟩_𝑚 be a duality pairing between 𝐌^∼ and 𝐖^∼0 defined as ⟨𝝀,𝐰⟩_𝑚=_{𝑖𝛿𝑚(𝑗)⊂Γ𝛿𝑚(𝑗)}𝝀_𝑖𝑗⋅𝐰_𝑖𝑗𝑑𝑠. (40) Then we define a dual norm for 𝝀∈𝐌^∼ by ‖𝝀‖^2∼_𝐌=max_{𝐰∈∼𝐖0⧵{0}}⟨𝝀,𝐰⟩²_𝑚‖𝐰‖²_𝐖0. (41) Now, we will find an operator 𝐹_DP which gives ⟨𝐹_DP𝝀,𝝀⟩=‖𝝀‖^2∼_𝐌 (42) and propose 𝐹⁻¹_DP as a preconditioner for 𝐹_DP in (28). Define 𝑅_𝑖𝑗∶𝐖⁰→𝐖⁰_𝑖𝑗 as a restriction operator and 𝐸^𝑖_𝑖𝑗∶𝐖⁰_𝑖𝑗→𝐖_𝑖 as an extension operator by zero. Then for 𝐰∈𝐖^∼0, ‖𝐰‖²_𝐖0=𝐰^∼2_𝐖=^𝑁_𝑖=1𝑆𝐰_𝑖^∼_𝑖,𝐰^∼_𝑖=^𝑁_𝑖=1𝑆_{𝑖𝛿𝑚(𝑗)⊂Γ}𝐸^𝑖_𝑖𝑗𝑅_𝑖𝑗𝐰,_{𝛿m(j)⊂Γ}𝐸^𝑖_𝑖𝑗𝑅_𝑖𝑗𝐰. (43) Let 𝑆=∑^𝑁_𝑖=1(∑_{𝛿𝑚(𝑗)⊂Γ}𝐸^𝑖_𝑖𝑗𝑅_𝑖𝑗)^𝑡𝑆_𝑖(∑_{𝛿𝑚(𝑗)⊂Γ}𝐸^𝑖_𝑖𝑗𝑅_𝑖𝑗). Moreover, we have ⟨𝝀,𝐰⟩_𝑚=𝐵𝐰,𝝀, (44) where 𝐵=diag_{𝑖=1,…,𝑁}diag_{𝛿𝑚(𝑗)⊂Γ}𝐵_𝑖^𝑖𝑗 and 𝐵_𝑖^𝑖𝑗 is a matrix obtained from 𝐵_𝑖^𝑖𝑗 after deleting the columns corresponding to the d.o.f. of ∫_Γ𝑖𝑗𝐰𝑑𝑠=𝟎. Note that 𝐵_𝑖^𝑖𝑗 is invertible. We define the following 𝐿²-orthogonal projections: 𝑃_∼𝐖⁰∶𝐖⁰⟶𝐖^∼0,𝑃^∼_𝐌∶𝐌⟶𝐌^∼. (45)

For 𝝀∈𝐌^∼ and 𝐰∈𝐖^∼0, we have ⟨𝝀,𝐰⟩_𝑚=𝐵_𝑝𝐰,𝝀,‖𝐰‖²_𝐖0=𝑆_𝑝𝐰,𝐰, (46) where 𝑆_𝑝=𝑃_∼𝐖⁰𝑆𝑃_∼𝐖⁰,𝐵_𝑝=𝑃^∼_𝐌𝐵𝑃_∼𝐖⁰. (47)

Then it can be shown that the operators 𝑆_𝑝∶𝐖^∼0→𝐖^∼0,𝐵_𝑝∶𝐖^∼0→𝐌^∼ are invertible and 𝑆_𝑝 is s.p.d on 𝐖^∼0. Hence, using (46), the maximum in (41) occurs when 𝐰∈𝐖^∼0 satisfies 𝑆_𝑝𝐰=𝐵^𝑡_𝑝𝝀. Therefore, we have ‖𝝀‖^2∼_𝐌=𝐵_𝑝𝑆_𝑝⁻¹𝐵^𝑡_𝑝𝝀,𝝀. (48) Let 𝐹⁻¹_DP=𝐵_𝑝𝑆_𝑝⁻¹𝐵^𝑡_𝑝⁻¹=𝐵^𝑡_𝑝⁻¹𝑆_𝑝𝐵_𝑝⁻¹ (49) be the preconditioner for the 𝐹_DP.

Define the 𝐿²-orthogonal projections 𝑃^𝑖𝑗_∼𝐖0∶𝐖⁰_{𝛿𝑚(𝑗)}⟶𝐖^∼0_{𝛿𝑚(𝑗)},𝑃^𝑖𝑗_∼𝐌∶𝐌|_{𝛿𝑚(𝑗)}⟶𝐌^∼_{𝛿𝑚(𝑗)}. (50)

Then the projection operators 𝑃_∼𝐖⁰ and 𝑃^∼_𝐌 are composed of diagonal blocks of 𝑃^𝑖𝑗_∼𝐖0 and 𝑃^𝑖𝑗_∼𝐌, respectively. Moreover, it can be shown easily that 𝑃^𝑖𝑗_∼𝐌𝐵_𝑖^𝑖𝑗𝑃^𝑖𝑗_∼𝐖0∶𝐖^∼0_{𝛿𝑚(𝑗)}⟶𝐌^∼_{𝛿𝑚(𝑗)} (51) is invertible. Hence, it follows that 𝐵_𝑝⁻¹=diag_{𝑖=1,…,𝑁}diag_{𝛿𝑚(𝑗)⊂Γ}𝐵⁻¹_𝑖𝑗, (52) where 𝐵_𝑖𝑗=𝑃^𝑖𝑗_∼𝐌𝐵_𝑖^𝑖𝑗𝑃^𝑖𝑗_∼𝐖0, and we obtain 𝐹⁻¹_DP=^𝑁_{𝑖=1𝛿𝑚(𝑗)⊂Γ}𝐸^𝑖_𝑖𝑗𝐵⁻¹_𝑖𝑗𝑅_𝑖𝑗^𝑡𝑆_{𝑖𝛿𝑚(𝑗)⊂Γ}𝐸^𝑖_𝑖𝑗𝐵⁻¹_𝑖𝑗𝑅_𝑖𝑗. (53)

5. Technical Tools

In this section, we state and prove a few technical lemmas necessary for the proof of Theorem 17 of Section 6.

Let 𝑌_ℎ/2(Ω_𝑖) be the conforming space of piecewise linear continuous functions on the triangulation 𝑇_ℎ/2(Ω_𝑖) which is constructed by joining the midpoints of the edges of the elements of 𝑇_ℎ(Ω_𝑖), and let 𝑌⁰_ℎ/2(Ω_𝑖) be the subspace of 𝑌_ℎ/2(Ω_𝑖) consisting of functions with zero traces on 𝜕Ω_𝑖.

For each open edge 𝜀⊂𝜕Ω_𝑖, we introduce a local equivalence map (isomorphism) 𝑀^𝜀_𝑖∶ 𝐖^ℎ_𝜀(Ω_𝑖)→ [𝑌_ℎ/2(Ω_𝑖)]² (cf. [11]), where 𝐖^ℎ_𝜀(Ω_𝑖) is a subspace of 𝐖^ℎ(Ω_𝑖) formed by all vectors whose component are zero in 𝜕Ω^CR_𝑖,ℎ⧵𝜀^CR_ℎ. This mapping is a slightly modified version of a similar mapping introduced in [12].

Definition 3. For given 𝐮∈𝐖^ℎ_𝜀(Ω_𝑖), we introduce 𝑀^𝜀_𝑖𝐮∈[𝑌_ℎ/2(Ω_𝑖)]² by defining the values of 𝑀^𝜀_𝑖𝐮 at the nodal points of the triangulation 𝑇_ℎ/2(Ω_𝑖).(i)For 𝑝∈Ω^CR_𝑖,ℎ∪𝜕Ω^CR_𝑖,ℎ, let 𝑀^𝜀_𝑖𝐮(𝑝)=𝐮(𝑝).(ii)For 𝑝∈Ω_𝑖,ℎ, then 𝑀^𝜀_𝑖𝐮(𝑝)=(1/𝑁(𝑝))∑_𝜏ℎ𝑗𝐮|_𝜏^ℎ_𝑗(𝑝), where the sum is taken over all triangles 𝜏^ℎ_𝑗 with the common vertex 𝑝, and 𝑁(𝑝) is the number of elements with 𝑝 as an vertex.(iii)For 𝑝∈𝜕Ω_𝑖,ℎ⧵𝜀_ℎ, let 𝑀^𝜀_𝑖𝐮(𝑝)=𝟎.(iv)For 𝑝∈𝜀_ℎ, then 𝑀^𝜀_𝑖𝐮(𝑝)=(|𝑝𝑝_𝑟|/|𝑝_𝑙𝑝_𝑟|)𝐮(𝑝_𝑙)+(|𝑝_𝑙𝑝|/|𝑝_𝑙𝑝_𝑟|)𝐮(𝑝_𝑟), where 𝑝_𝑙,𝑝_𝑟 are the left- and right-neighboring CR nodal points of 𝑝 and |𝑎𝑏| is the length of the segment with 𝑎,𝑏 as its endpoints.

Note that each component of 𝑀^𝜀_𝑖𝐮 is piecewise linear between the CR nodes of 𝜀^CR_ℎ, and 𝑀^𝜀_𝑖𝐮=𝟎 on 𝜕Ω_𝑖⧵𝜀. The mapping 𝑀^𝜀_𝑖𝐮 has the following properties (cf. [11]): |𝐮|_𝐻¹_ℎ(Ω𝑖)≍𝑀^𝜀_𝑖𝐮_{𝐻1(Ω𝑖)},‖𝐮‖_𝐿²_(Ω𝑖)≍𝑀^𝜀_𝑖𝐮_{𝐿2(Ω𝑖)},𝑀^𝜀_𝑖𝐮−𝐮_{𝐿2(Ω𝑖)}⪯ℎ_𝑖|𝐮|_𝐻¹_ℎ(Ω𝑖),𝑀^𝜀_𝑖𝐮−𝐮_𝐿2(𝜀)⪯ℎ_𝑖^1/2|𝐮|_𝐻¹_ℎ(Ω𝑖). (54)

Lemma 4 (see [13]). Let 𝐮∈𝐖_𝑗 and 𝐮=𝟎 in 𝜕Ω^CR_𝑗,ℎ⧵𝛿^CR_{𝑚(𝑗),ℎ} for a slave edge 𝛿_𝑚(𝑗)⊂𝜕Ω_𝑗. Then we have 𝑆𝐻_𝑗𝐮_{𝐻1ℎ(Ω𝑗)}⪯𝑀^{𝛿𝑚(𝑗)}_𝑗𝑆𝐻_𝑗𝐮_{𝐻1/200(𝛿𝑚(𝑗))}, (55) where 𝑆𝐻_𝑗𝐮 is the discrete harmonic extension of 𝐮.

Lemma 5 (see [11]). Let 𝛿_𝑚(𝑗)⊂𝜕Ω_𝑗 be a slave, then for any 𝐮∈𝐖^ℎ(Ω_𝑗), we have 𝐮^{𝛿𝑚(𝑗)2}_{𝐻1ℎΩ𝑗}⪯𝐻1+log^𝑗ℎ²×𝐻_𝑗⁻²‖𝐮‖²_𝐿2Ω𝑗+|𝐮|²_{𝐻1ℎΩ𝑗}, (56) where 𝐮^{𝛿𝑚(𝑗)} is a discrete harmonic function taking the same vector values as those of 𝐮 at CR nodal points on 𝛿^CR_{𝑚(𝑗),ℎ} and is equal to zero vector at the remaining CR nodal points on 𝜕Ω_𝑗 and 𝐻_𝑗=diamΩ_𝑗.

We also need a special mortar operator defined over trace spaces: Π_{𝛿𝑚(𝑗)}∶[𝐿²(𝛿_𝑚(𝑗))]²→𝐖₀(𝛿_𝑚(𝑗)), where 𝐖₀(𝛿_𝑚(𝑗)) is a set of vectors whose components are continuous and equal to zero at the endpoints of 𝛿_𝑚(𝑗) and are piecewise linear over all segments that have their ends in 𝛿^CR_{𝑚(𝑗),ℎ}. Each Π_{𝛿𝑚(𝑗)}𝐮∈𝐖₀(𝛿_𝑚(𝑗)) is uniquely determined by its values in 𝛿^CR_{𝑚(𝑗),ℎ}.

Definition 6 (see [11]). Let Π_{𝛿𝑚(𝑗)}∶[𝐿²(𝛿_𝑚(𝑗))]²→𝐖₀(𝛿_𝑚(𝑗)) be defined by Π_{𝛿𝑚(𝑗)}𝐮(𝑥)=𝑄_𝑚𝐮(𝑥),∀𝑥∈𝛿^CR_𝑚𝑗,ℎ. (57)

The 𝐿² and 𝐻^1/2₀₀ stability of Π_{𝛿𝑚(𝑗)} is stated in the next lemma.

Lemma 7 (see [11]). For Π_{𝛿𝑚(𝑗)} defined previously, we have Π_{𝛿𝑚(𝑗)}𝐮_{𝐿2(𝛿𝑚(𝑗))}⪯‖𝐮‖_𝐿²_{(𝛿𝑚(𝑗))}∀𝐮∈𝐿²𝛿_𝑚(𝑗)²,Π_{𝛿𝑚(𝑗)}𝐮_{𝐻1/200(𝛿𝑚(𝑗))}⪯‖𝐮‖_𝐻^1/2_{00(𝛿𝑚(𝑗))}∀𝐮∈𝐻^1/2₀₀𝛿_𝑚(𝑗)². (58)

Lemma 8 (see [10]). Let 𝐮∈𝐖^ℎ(Ω_𝑗) be a discrete harmonic vector in Ω_𝑗, with 𝐮=𝟎 at 𝑥∈ 𝜕Ω^CR_𝑗,ℎ⧵ 𝛿^CR_{𝑚(𝑗),ℎ}. Then we have |𝐮|_𝐻¹_ℎ(Ω𝑗)⪯Π_{𝛿𝑚(𝑗)}𝐮_{𝐻1/200(𝛿𝑚(𝑗))}. (59)

Lemma 9 (see [14]). The matrix 𝐾_ΔΔ is invertible on 𝑊_𝐼×𝑄_𝐼×𝑊_Δ.

Lemma 10 (see [14]). Assume that the domain Ω has the triangulation to satisfy that −_𝑖𝑝⁰,∇⋅𝐯_ΔΩ𝑖+_𝑖,𝑗𝝀⋅_Γ𝑖𝑗𝐯^𝑖_Δ−𝐯^𝑗_Δ𝑑𝑠=0,∀𝐯_Δ∈𝐖_Δ,−_𝑖𝑝⁰,∇⋅𝐯_𝑐Ω𝑖+_𝑖,𝑗𝝀⋅_Γ𝑖𝑗𝐯^𝑖_𝑐−𝐯^𝑗_𝑐𝑑𝑠=0,∀𝐯_𝑐∈𝐖_𝑐 (60) for the solution 𝑝⁰𝝀. Then the coarse problem matrix 𝐹_𝑐𝑐 is invertible.

Remark 11. Most triangulations of the domain Ω satisfy the assumption of Lemma 10 because the number of velocity unknowns 𝐯_Δ and 𝐯_𝑐 is usually greater than the number of unknowns 𝑝⁰ and 𝝀.

Lemma 12 (see [14]). We have 𝐵(𝐖^∼_Δ)=𝐵(𝐖^∼)=𝐌^∼.

Lemma 13 (see [14]). For 𝝀∈𝐌^∼, we have 𝐹_DP𝝀,𝝀=max_{𝐰∈∼𝐖Δ⧵{𝟎}}⟨𝐵𝐰,𝝀⟩²‖𝐰‖²_𝐖. (61)

Lemma 14 (see [15]). For any 𝝀∈𝐌, we have 𝐹_DP𝝀,𝝀^1/2=sup_{𝐰∈∼𝐖Δ⧵{𝟎}𝛿𝑚(𝑗)⊂Γ}∫_{𝛿𝑚(𝑗)}𝐰𝝀𝑑𝑠‖𝐰‖_𝐖. (62)

Proof. The proof follows from the definitions of 𝐖 and 𝐖^∼_Δ, (40), and a standard algebraic argument; see [16].

Definition 15 (see [17]). Given 𝐮∈𝐖_𝑖, we define an operator 𝑃_𝑖∶𝐖_𝑖→𝐖^∼ as follows:(i)for 𝑝∈𝛾^CR_{𝑚(𝑖),ℎ}, let 𝑃_𝑖𝐮(𝑝)=𝐮(𝑝);(ii)for 𝑝∈𝛿^CR_{𝑚(𝑗),ℎ}, let 𝑃_𝑖𝐮(𝑝)=𝑄_𝑚𝐽_𝑚𝐮(𝑝);(iii) for 𝑝∈𝛿^CR_{𝑚(𝑖),ℎ}, for any slave side 𝛿_𝑚(𝑖)⊂𝜕Ω_𝑖, or if the edge midpoint 𝑝 is on 𝜕Ω or on any remaining master or slave side, let 𝑃_𝑖𝐮(𝑝)=𝟎.

We define a projection operator 𝑃∶𝐖^∼→𝐖^∼_Δ by (𝑃𝐮)|_{𝛿𝑚(𝑗)}=𝑄_𝑚𝐽_𝑚𝐮_𝑖−𝐮_𝑗on𝛿_𝑚(𝑗). (63) where 𝐮_𝑖 and 𝐮_𝑗 are, respectively, the restriction of 𝐮∈𝐖^∼ to the mortar side 𝛾_𝑚(𝑖) and the slave side 𝛿_𝑚(𝑗) of an interface Γ_𝑖𝑗.

Lemma 16. For all 𝐮∈𝐖^∼, we have ‖𝑃𝐮‖_𝐖⪯𝐻1+logℎ‖𝐮‖_𝐖, (64) where 𝐻=max_𝑖𝐻_𝑖 and ℎ=min_𝑖ℎ_𝑖.

Proof. Let 𝑃𝐮=(𝑃₁𝐮,…,𝑃_𝑁𝐮)∈𝐖^∼, and we have ‖𝑃𝐮‖²_𝐖=∑^𝑁_𝑖=1‖𝑃_𝑖𝐮‖²_𝐖𝑖. Next we will estimate the term ‖𝑃_𝑖𝐮‖²_𝐖𝑖.
Note that 𝑃_𝑖𝐮 can be nonzero only over Ω_𝑖 and the neighboring subdomain Ω_𝑗 that share a common edge Γ_𝑖𝑗 with Ω_𝑖, such that the slave side 𝛿_𝑚(𝑗) is associated with Ω_𝑗. Let 𝑁_𝑖 denote the set of indices of such subdomains. Thus 𝑃_𝑖𝐮²_𝐖𝑖=𝑆𝐻_𝑖𝑃_𝑖𝐮²_{𝐻1ℎΩ𝑖}+_{𝑗∈𝑁𝑖}𝑆𝐻_𝑗𝑃_𝑖𝐮²_{𝐻1ℎΩ𝑗}. (65)
Then considering Ω_𝑖, we have 𝑃_𝑖𝐮=𝐮−∑_{𝛿𝑚(𝑖)⊂𝜕Ω𝑖}𝐮^{𝛿𝑚(𝑖)} with 𝐮^{𝛿𝑚(𝑖)} being a discrete harmonic function which equals to 𝐮 at the CR nodes of 𝛿^CR_{𝑚(𝑖),ℎ} and to zero at all remaining nodes of 𝜕Ω^CR_𝑖,ℎ. Thus, by Lemma 5, 𝑆𝐻_𝑖𝑃_𝑖𝐮²_{𝐻1ℎΩ𝑖}⪯𝑆𝐻_𝑖𝐮²_{𝐻1ℎΩ𝑖}+_{𝛿𝑚(𝑖)⊂𝜕Ω𝑖}𝐮^{𝛿𝑚(𝑖)2}_{𝐻1ℎΩ𝑖}⪯𝐻1+log^𝑖ℎ_𝑖²𝐻_𝑖⁻²𝑆𝐻_𝑖𝐮²_𝐿2Ω𝑖+𝑆𝐻_𝑖𝐮²_{𝐻1ℎΩ𝑖}. (66)Next we consider those Ω_𝑗 for which 𝑗∈𝑁_𝑖, consequently the slave side 𝛿_𝑚(𝑗) of the edge Γ_𝑖𝑗 is an edge of Ω_𝑗, and the mortar side 𝛾_𝑚(𝑖) is an edge of Ω_𝑖. Then by Lemma 4, we get 𝑆𝐻_𝑗𝑃_𝑖𝐮_{𝐻1ℎ(Ω𝑗)}⪯𝑀^{𝛿𝑚(𝑗)}_𝑗𝑆𝐻_𝑗𝑃_𝑖𝐮_{𝐻1/200(𝛿𝑚(𝑗))}. (67) We also need to consider the function Π_{𝛿𝑚(𝑗)}𝑀^{𝛾𝑚(𝑖)}_𝑖𝐮^{𝛾𝑚(𝑖)}. From Definitions 3 and 15, we have 𝑀^{𝛿𝑚(𝑗)}_𝑗𝑆𝐻_𝑗𝑃_𝑖𝐮𝑝=𝑆𝐻_𝑖𝑃_𝑖𝐮𝑝=𝑄_𝑚𝑆𝐻_𝑖𝑃_𝑖𝐮𝑝=𝑄_𝑚𝐽_𝑚𝑆𝐻_𝑖𝐮𝑝,∀𝑝∈𝛿^CR_𝑚𝑗,ℎ. (68) Thus we get 𝑀^{𝛿𝑚(𝑗)}_𝑗𝑆𝐻_𝑗𝑃_𝑖𝐮_{𝐻1/200(𝛿𝑚(𝑗))}≤𝑀^{𝛿𝑚(𝑗)}_𝑗𝑆𝐻_𝑗𝑃_𝑖𝐮−Π_{𝛿𝑚(𝑗)}𝑀^{𝛾𝑚(𝑖)}_𝑖𝐮^{𝛾𝑚(𝑖)}_{𝐻1/200(𝛿𝑚(𝑗))}+Π_{𝛿𝑚(𝑗)}𝑀^{𝛾𝑚(𝑖)}_𝑖𝐮^{𝛾𝑚(𝑖)}_{𝐻1/200(𝛿𝑚(𝑗))}. (69) The pervious second term can be estimated as Π_{𝛿𝑚(𝑗)}𝑀^{𝛾𝑚(𝑖)}_𝑖𝐮^{𝛾𝑚(𝑖)2}_{𝐻1/200𝛿𝑚(𝑗)}⪯𝑀^{𝛾𝑚(𝑖)}_𝑖𝐮^{𝛾𝑚(𝑖)2}_{𝐻1/200𝛿𝑚(𝑗)}⪯𝑀^{𝛾𝑚(𝑖)}_𝑖𝐮^{𝛾𝑚(𝑖)2}_𝐻1Ω𝑖⪯𝐮^{𝛾𝑚(𝑖)2}_{𝐻1ℎΩ𝑖}⪯𝐻1+log^𝑖ℎ_𝑖²𝐻_𝑖⁻²𝑆𝐻_𝑖𝐮²_𝐿2Ω𝑖+𝑆𝐻_𝑖𝐮²_{𝐻1ℎΩ𝑖}, (70) by Lemmas 5 and 7 and the trace theorem. The first term is bounded using an inverse inequality, (68) and Definition 6 as follows: 𝑀^{𝛿𝑚(𝑗)}_𝑗𝑆𝐻_𝑗𝑃_𝑖𝐮−Π_{𝛿𝑚(𝑗)}𝑀^{𝛾𝑚(𝑖)}_𝑖𝐮^{𝛾𝑚(𝑖)2}_{𝐻1/200𝛿𝑚(𝑗)}⪯ℎ_𝑗⁻¹𝑀^{𝛿𝑚(𝑗)}_𝑗𝑆𝐻_𝑗𝑃_𝑖𝐮−Π_{𝛿𝑚(𝑗)}𝑀^{𝛾𝑚(𝑖)}_𝑖𝐮^{𝛾𝑚(𝑖)2}_{𝐿2𝛿𝑚(𝑗)}⪯_{𝑝∈𝛿CR𝑚(𝑗),ℎ}𝑀^{𝛿𝑚(𝑗)}_𝑗𝑆𝐻_𝑗𝑃_𝑖𝐮𝑝−Π_{𝛿𝑚(𝑗)}𝑀^{𝛾𝑚(𝑖)}_𝑖𝐮^{𝛾𝑚(𝑖)}𝑝²=_{𝑝∈𝛿CR𝑚(𝑗),ℎ}𝑄_𝑚𝐽_𝑚𝑆𝐻_𝑖𝐮−𝑀^{𝛾𝑚(𝑖)}_𝑖𝐮^{𝛾𝑚(𝑖)}𝑝²⪯ℎ_𝑗⁻¹𝑄_𝑚𝐽_𝑚𝑆𝐻_𝑖𝐮−𝑀^{𝛾𝑚(𝑖)}_𝑖𝐮^{𝛾𝑚(𝑖)2}_{𝐿2𝛿𝑚(𝑗)}. (71)Next, by using the fact that 𝑄_𝑚 is an orthogonal 𝐿² projection, we get 𝑄_𝑚𝑀^{𝛾𝑚(𝑖)}_𝑖𝐮^{𝛾𝑚(𝑖)}−𝐽_𝑚𝑆𝐻_𝑖𝐮_{𝐿2(𝛿𝑚(𝑗))}⪯𝑀^{𝛾𝑚(𝑖)}_𝑖𝐮^{𝛾𝑚(𝑖)}−𝐽_𝑚𝑆𝐻_𝑖𝐮_{𝐿2(𝛿𝑚(𝑗))}. (72) Finally, the fact that 𝐽_𝑚𝑆𝐻_𝑖𝐮=𝐽_𝑚𝐮^{𝛾𝑚(𝑖)} yields 𝑀^{𝛾𝑚(𝑖)}_𝑖𝐮^{𝛾𝑚(𝑖)}−𝐽_𝑚𝑆𝐻_𝑖𝐮_{𝐿2(𝛿𝑚(𝑗))}≤𝑀^{𝛾𝑚(𝑖)}_𝑖𝐮^{𝛾𝑚(𝑖)}−𝐮^{𝛾𝑚(𝑖)}_{𝐿2(𝛿𝑚(𝑗))}+𝐮^{𝛾𝑚(𝑖)}−𝐽_𝑚𝐮^{𝛾𝑚(𝑖)}_{𝐿2(𝛿𝑚(𝑗))}. (73) Applying the property of 𝑀^{𝛾𝑚(𝑖)}_𝑖𝐮 to the first term and Lemma 3.2 from [10] to the second term, Lemma 5 and the quasiuniform triangulation, we get ℎ_𝑗⁻¹𝑀^{𝛾𝑚(𝑖)}_𝑖𝐮^{𝛾𝑚(𝑖)}−𝐽_𝑚𝑆𝐻_𝑖𝐮²_{𝐿2𝛿𝑚(𝑗)}≤ℎ_𝑖ℎ_𝑗𝐮^{𝛾𝑚(𝑖)}_{𝐻1ℎ(Ω𝑖)}⪯𝐻1+log^𝑖ℎ_𝑖²𝐻_𝑖⁻²𝑆𝐻_𝑖𝐮²_𝐿2Ω𝑖+𝑆𝐻_𝑖𝐮²_{𝐻1ℎΩ𝑖}⪯𝐻1+log^𝑖ℎ_𝑖²𝑆𝐻_𝑖𝐮²_{𝐻1ℎΩ𝑖}. (74) Finally, summing over all edges in 𝜕Ω_𝑖 and then over all subdomain ends the proof.

6. Condition Number Estimate

In this section, we give the condition number estimate of the preconditioned operator in the following main theorem of this paper.

Theorem 17. For any 𝝀∈𝐌, it holds that 𝐹_DP𝝀,𝝀⪯𝐹_DP⪯𝐻𝝀,𝝀1+logℎ²𝐹_DP𝝀,𝝀, (75) where 𝐻=max_𝑖𝐻_𝑖 and ℎ=min_𝑖ℎ_𝑖.

Proof. In the proof of this theorem we use results from this section and the algebraic arguments from [16, Section 6.4].
Lower Bound. From Lemma 13, the zero extension 𝐸(𝐰₀)∈𝐖^∼_Δ for 𝐰₀∈𝐖^∼0, and (41), we obtain 𝐹_DP𝝀,𝝀=sup_{𝐰∈∼𝐖Δ⧵{𝟎}}⟨𝐵𝐰,𝝀⟩²‖𝐰‖²_𝐖≥sup_{𝐰0∈∼𝐖0⧵{𝟎}}𝐵𝐸(𝐰₀),𝝀²𝐸(𝐰₀)²_𝐖=𝐹_DP𝝀,𝝀. (76)Upper Bound. For any 𝐰∈𝐖^∼, ⟨𝐵𝐰,𝝀⟩=_{𝛿𝑚(𝑗)⊂Γ𝛿𝑚(𝑗)}𝐽_𝑚𝐮_𝑖−𝐮_𝑗⋅𝝀_𝑖𝑗𝑑𝑠=_{𝛿𝑚(𝑗)⊂Γ𝛿𝑚(𝑗)}𝑄_𝑚𝐽_𝑚𝐮_𝑖−𝐮_𝑗⋅𝝀_𝑖𝑗𝑑𝑠=_{𝛿𝑚(𝑗)⊂Γ𝛿𝑚(𝑗)}(𝑃𝐰)|_{𝛿𝑚(𝑗)}⋅𝝀_𝑖𝑗𝑑𝑠. (77)Hence by Lemmas 13, 14, and 16, we conclude that 𝐹_DP𝝀,𝝀=sup_{𝐰∈∼𝐖Δ⧵{𝟎}}⟨𝐵𝐰,𝝀⟩²‖𝐰‖²_𝐖=sup_{𝐰∈∼𝐖Δ⧵{𝟎}}∑_{𝛿𝑚(𝑗)⊂Γ}∫_{𝛿𝑚(𝑗)}(𝑃𝐰)|_{𝛿𝑚(𝑗)}⋅𝝀_𝑖𝑗𝑑𝑠²‖𝐰‖²_𝐖⪯𝐹_DP𝝀,𝝀sup_{𝐰∈∼𝐖Δ⧵{𝟎}}‖𝑃𝐰‖²_𝐖‖𝐰‖²_𝐖⪯𝐻1+logℎ²𝐹_DP𝝀,𝝀. (78)

7. Numerical Tests

Consider the following Stokes problem: −Δ𝐮+∇𝑝=𝐟inΩ,−∇⋅𝐮=0inΩ,𝐮=𝟎on𝜕Ω, (79) where Ω=[0,1]×[0,1]⊂𝑅² and 𝐟 is given by the exact solution of the problem (79) as follows: 𝐮=sin³(𝜋𝑥)sin²𝜋𝑦cos𝜋𝑦−sin²(𝜋𝑥)sin³𝜋𝑦cos(𝜋𝑥),𝑝=𝑥²−𝑦². (80)

Let 𝑁 denote the number of subdomains. We only consider the uniform partition of Ω.𝑁=4×4 means that Ω is partitioned into 4×4 square subdomains. For all subdomains, we take the same number of nodes 𝑛, including endpoints, in horizontal and vertical edges with 𝑛=4𝑘+1 for some positive integer 𝑘. We solve (79) on both matching and nonmatching grids. For matching grids, we make uniform triangulations in each subdomain with (𝑛−1)/2+1 nodes on horizontal and vertical edges of subdomain. For nonmatching grids, we take (𝑛−1)/2+1 random quasi-uniform nodes on each horizontal and vertical edges of subdomain and generate nonuniform structured triangulations.

Now, we solve the FETI-DP operator with and without preconditioner varying 𝑁 and 𝑛 on both matching and nonmatching triangulations. The Conjugate Gradient (CG) iteration is stopped when the relative residual is less than 10⁻⁶.

In Tables 1–3, the number of CG iterations and the corresponding condition numbers are shown with varying 𝑁 and 𝑛. In Table 1, 𝑁=4×4 and 𝑛−1 increase by double. On both matching and nonmatching grids, the preconditioner performs well and the condition numbers seem to behave log²-growth as 𝑛 increases. Especially on nonmatching grids, the preconditioner is much more efficient. In Tables 2 and 3, the CG iteration becomes stable as 𝑁 increases with 𝑛=5 and 𝑛=9.

Table 1 

CG iterations (condition numbers) when 𝑁=4×4.

Table 2 

CG iterations (condition numbers) when 𝑛=5.

Table 3 

CG iterations (condition numbers) when 𝑛=9.

Moreover, we have observed the convergence behaviors of the approximated solution. The 𝐻¹ and 𝐿²-errors for velocity and pressure are examined. 𝐮^ℎ and 𝑝^ℎ denote the approximated solutions for the velocity and pressure, and ‖𝐮−𝐮^ℎ‖_1,∗ means the square root of ∑^𝑁_𝑖=1‖𝐮−𝐮^ℎ‖²_1,Ω𝑖. The errors are shown in Table 4 for various 𝑁 and 𝑛 with matching grids. Three cases are considered: when 𝑛−1 increases by double with 𝑁=4×4, when 𝑁 increases by double in both edges of Ω with 𝑛=5 and 𝑛=9. For all cases, we can see that the 𝐻¹-error for velocity ‖𝐮−𝐮^ℎ‖_1,∗,𝐿²-error for pressure ‖𝑝−𝑝^ℎ‖₀ reduce by half and 𝐿²-error for velocity ‖𝐮−𝐮^ℎ‖₀ reduces by quarter. For the finite elements CR(ℎ)−𝑃₀(ℎ), these convergence behaviors are optimal.

Table 4 

𝐻¹- and 𝐿²-errors on matching grids.

For the case of nonmatching grids, the errors are shown in Tables 5–7 with various 𝑁 and 𝑛. In Table 5, we observe that the error ‖𝐮−𝐮^ℎ‖_1,∗ and ‖𝑝−𝑝^ℎ‖₀ reduce by half and the error ‖𝐮−𝐮^ℎ‖₀ reduces by quarter as 𝑛−1 increases by double with 𝑁=4×4. When 𝑛=5 and 𝑛=9, as 𝑁 increases, the errors also show the optimal convergence behaviors in Tables 6 and 7.

Table 5 

𝐻¹- and 𝐿²-errors on nonmatching grids: 𝑁=4×4.

Table 6 

𝐻¹- and 𝐿²-errors on nonmatching grids: 𝑛=5.

Table 7 

𝐻¹- and 𝐿²-errors on nonmatching grids: 𝑛=9.

If the inf-sup constant for the space 𝑉×𝑃 is independent of 𝑁 and 𝑛, then the optimality of approximation can be shown. Let 𝛽^∗ and 𝛽 be the inf-sup constants for the space 𝑉×𝑃 and the CR(ℎ)−𝑃₀(ℎ) finite elements, respectively, and 𝛽₀ be the inf-sup constant for the space 𝑉×Π⁰. Then the constant 𝛽^∗ depends on 𝛽 and 𝛽₀ from the trick conceived by Boland and Nicolaides [18]. Hence, if the constant 𝛽₀ is independent of 𝑛 and 𝑁, then the same holds for 𝛽^∗. In [19], for 𝑉×Π⁰ which is obtained from the Hood-Taylor finite elements, it was shown that the constant 𝛽₀ is independent of 𝑛, but not shown for 𝑁. Following the proofs in [19], we can obtain the same results for the space 𝑉×Π⁰ of the CR(ℎ)−𝑃₀(ℎ) finite elements. We have no proof that 𝛽₀ is independent of 𝑁. Instead, we compute the constant 𝛽₀ numerically as 𝑁 increases. The results are given in Table 8 both for matching and nonmatching grids when 𝑛=5 and 𝑛=9. We observe that the constant 𝛽₀ becomes stable as 𝑁 increases. Table 9 gives the constant 𝛽₀ as 𝑛 increases with 𝑁=4×4. This confirms that the constant 𝛽₀ is independent of 𝑛.

Table 8 

Inf-sup constant 𝛽₀ when 𝑛=5 and 𝑛=9.

Table 9 

Inf-sup constant 𝛽₀ when 𝑁=4×4.

Acknowledgments

The author of this paper would like to express her gratitude to Professor Jinru Chen, for his guidance and encouragement, and to the reviewers, for valuable comments and suggestions. This work was supported by the National Science Foundation (NSF) of China (Grants No. 11071124 and 11226334) and the Program of Natural Science Research of Jiangsu Higher Education Institutions of China (Grant No. 12KJB110013).