MORTAR MULTISCALE FINITE ELEMENT ... - Semantic Scholar

MORTAR MULTISCALE FINITE ELEMENT METHODS FOR STOKES-DARCY FLOWS VIVETTE GIRAULT, DANAIL VASSILEV, AND IVAN YOTOV Abstract. We investigate mortar multiscale numerical methods for coupled Stokes and Darcy flows with the Beavers–Joseph–Saffman interface condition. The domain is decomposed into a series of subdomains (coarse grid) of either Stokes or Darcy type. The subdomains are discretized by appropriate Stokes or Darcy finite elements. The solution is resolved locally (in each coarse element) on a fine scale, allowing for non-matching grids across subdomain interfaces. Coarse scale mortar finite elements are introduced on the interfaces to approximate the normal stress and impose weakly continuity of the normal velocity. Stability and a priori error estimates in terms of the fine subdomain scale h and the coarse mortar scale H are established for fairly general grid configurations, assuming that the mortar space satisfies a certain inf-sup condition. Several examples of such spaces in two and three dimensions are given. Numerical experiments are presented in confirmation of the theory.

1. Introduction Mathematical and numerical modeling of coupled Stokes and Darcy flows has become a topic of significant interest in recent years. Such coupling occurs in many applications, including surface water-groundwater interaction, flows through vuggy or fractured porous media, industrial filters, fuel cells, and cardiovascular flows. The most commonly used model is based on the experimentally derived Beavers–Joseph–Saffman interface condition [10, 55], a slip with friction condition for the Stokes flow with a friction coefficient that depends on the permeability of the adjacent porous media. Existence and uniqueness of a weak solution has been studied in [46, 24]. Numerous stable and convergent numerical methods have been developed, see, e.g., [46, 24, 54, 29, 44, 28, 49, 32, 33] for methods based on different numerical discretizations suitable for each region, and [47, 4, 19, 59, 45, 9] for approaches employing unified finite elements. The full Beavers–Joseph condition was considered in [20, 2]. A coupling of Stokes-Darcy flows with transport was analyzed in [58]. The nonlinear system of coupled Navier-Stokes and Darcy flows has been studied in [35, 26, 8]. In this paper we develop multiscale mortar methods for multi-domain non-matching grid discretizations of Stokes-Darcy flows in two and three dimensions. Non-matching grids provide flexibility in meshing complex geometries with relatively simple locally constructed subdomain grids that are suitable for the choice of subdomain discretization methods. Mortar finite elements play the role of Lagrange multipliers to impose weakly interface conditions. In [46], a Lagrange multiplier approximating the normal stress was introduced to impose continuity of the normal velocity for discretizations involving mixed finite element methods for Darcy and conforming Stokes elements. With a choice of the Lagrange multiplier space as the normal trace of the Darcy velocity space, the analysis in [46] applied to non-matching grids on the Stokes–Darcy interface, although this was not explicitly noted. A similar choice was considered in subsequent mortar discretizations for Stokes–Darcy flows [54, 29, 14]. Mortar methods for mixed finite element discretizations for Darcy have been studied in [60, 5, 52, 6]. The analysis in these papers allows for the mortar grid to be Date: March 7, 2013. Key words and phrases. Stokes-Darcy flows, mortar finite element, mixed finite element, multiscale finite element. The second and third authors were partially supported by the NSF grant DMS 1115856 and the DOE grant DEFG02-04ER25618. 1

2

different from the traces of the subdomain grids with the assumption that the mortar space satisfies a suitable solvability condition that limits the number of mortar degrees of freedom. Mortar discretizations for Stokes have been developed in [11, 12]. There, the mortar grid was chosen to be the trace of one of the neighboring subdomain grids, similar to the choice in mortar methods for conforming Galerkin discretizations for second order elliptic problems [13]. In this work we allow for non-matching grid interfaces of Stokes–Darcy, Stokes–Stokes, and Darcy–Darcy types. We develop multiscale discretizations, where the subdomains are discretized on a fine scale h and the mortar space is discretized on a coarse scale H. Our method is based on saddle-point formulations in both regions and employs inf-sup stable mixed finite elements for Darcy and conforming elements for Stokes. The mortars approximate different physical variables and are used to impose different matching conditions depending on the type of interface. On Stokes–Stokes interfaces, the mortar functions represent the entire stress vector and impose weak continuity of the entire velocity vector. On Stokes–Darcy and Darcy–Darcy interfaces, the mortars approximate the normal stress, which is just the pressure in the Darcy region, and impose weak continuity of the normal velocity. The mortar spaces are assumed to satisfy suitable inf-sup conditions, allowing for very general subdomain and mortar grid configurations. We consider approximations of different polynomial degrees on the three types of interfaces and the two types of subdomains. The mortar spaces can be continuous or discontinuous, the latter providing localized mass conservation across interfaces. Our method is more general than existing Stokes–Stokes mortar methods [11, 12] and Stokes–Darcy mortar methods [46, 54, 29, 14]. On Darcy–Darcy interfaces, our condition is closely related to the solvability condition considered in [60, 5, 52, 6]. The stability and convergence analysis relies on the construction of a bounded global interpolant in the space of weakly continuous velocities that also preserves the velocity divergence in the usual discrete sense. This is done in two steps, starting from suitable local interpolants and correcting them to satisfy the interface matching conditions. The correction step requires the existence of bounded mortar interpolants. This is a very general condition that can be easily satisfied in practice. We present two examples in 2-D and one example in 3-D that satisfy this solvability condition. Our error analysis shows that the global velocity and pressure errors are bounded by the fine scale local approximation error and the coarse scale non-conforming error. Since the polynomial degrees on subdomains and interfaces may differ, one can choose higher order mortar polynomials to balance the fine scale and the coarse scale error terms and obtain fine scale asymptotic convergence. The dependence of the stability and convergence constants on the subdomain size is explicitly determined. In particular, the stability and fine scale convergence constants do not depend on the size of subdomains, while the coarse scale non-conforming error constants deteriorate when the subdomain size goes to zero. This is to be expected, as the relative effect of the non-conforming error becomes more significant in such regime. However, this dependence can be made negligible by choosing higher order mortar polynomials, as mentioned above. Our multiscale Stokes–Darcy formulation can be viewed as an extension of the mortar multiscale mixed finite element (MMMFE) method for Darcy developed in [6]. The MMMFE method provides an alternative to other multiscale methods in the literature such as the variational multiscale method [41, 3] and the multiscale finite element method [40, 22]. All three methods utilize a divide and conquer approach: solve relatively small fine scale subdomain problems that are only coupled on the coarse scale through a reduced number of degrees of freedom. The mortar multiscale approach is more flexible as it allows for employment of a posteriori error estimation to adaptively refine the mortar grids where necessary to improve the global accuracy [6]. Following the non-overlapping domain decomposition approach from [37], it can be shown that the global Stokes–Darcy problem can be reduced to a positive definite coarse scale interface problem [57]. The latter can be solved using a preconditioned Krylov space solver requiring Stokes or Darcy subdomain solves at each iteration. An alternative more efficient implementation for MMMFE discretizations for Darcy was developed in [31], where a multiscale flux basis giving the interface flux response for each coarse scale

3

mortar degree of freedom is precomputed. The multiscale flux basis is used to replace the solution of subdomain problems by a simple linear combination. The application of this methodology to the Stokes–Darcy interface problem will be discussed in a forthcoming paper. We should mention that there have been a number of papers in the literature studying domain decomposition methods for the Stokes–Darcy problem, primarily in the two-subdomain case, see, e.g., [25, 27, 39, 30, 21]. 1.1. Notation and preliminaries. Let Ω be a bounded, connected Lipschitz domain of IRn , n = 2, 3, with boundary ∂Ω and exterior unit normal vector n, and let Γ be a part of ∂Ω with positive n − 1 measure: |Γ| > 0. We do not assume that Γ is connected, but if it is not connected, we assume that it has a finite number of connected components. In the case when n = 3, we also assume that Γ is itself Lipschitz. Let 1 H0,Γ (Ω) = {v ∈ H 1 (Ω) ; v|Γ = 0}.

1 (Ω) reads: There exists a constant P depending only on Ω and Γ such Poincar´e’s inequality in H0,Γ Γ that

(1.1)

1 ∀v ∈ H0,Γ (Ω) , kvkL2 (Ω) ≤ PΓ |v|H 1 (Ω) .

The norms and spaces are made precise later on. The formula (1.1) is a particular case of a more general result (cf. [50, 15]): Proposition 1.1. Let Ω be a bounded, connected Lipschitz domain of IRn and let Φ be a seminorm on H 1 (Ω) satisfying: 1) There exists a constant P1 such that (1.2)

∀v ∈ H 1 (Ω) , Φ(v) ≤ P1 kvkH 1 (Ω) .

2) The condition Φ(c) = 0 for a constant function c holds if and only if c = 0. Then there exists a constant P2 depending only on Ω, such that
1 (1.3) ∀v ∈ H 1 (Ω) , kvkL2 (Ω) ≤ P2 |v|2H 1 (Ω) + Φ(v)2 2 .

We recall Korn’s first inequality: There exists a constant C1 depending only on Ω and Γ such that (1.4)

1 ∀v ∈ H0,Γ (Ω)n , |v|H 1 (Ω) ≤ C1 kD(v)kL2 (Ω) ,

where D(v) is the deformation rate tensor, also called the symmetric gradient tensor:
1 D(v) = ∇ v + ∇ v T . 2 This is a particular case of the following more general result (see (1.6) in [16]): Proposition 1.2. Let Ω be a bounded, connected Lipschitz domain of IRn and let Φ be a seminorm on H 1 (Ω)n satisfying: 1) There exists a constant C2 such that (1.5)

∀v ∈ H 1 (Ω)n , Φ(v) ≤ C2 kvkH 1 (Ω) .

2) The condition Φ(m) = 0 for a rigid-body motion m holds if and only if m is a constant vector. Then there exists a constant C3 depending only on Ω, such that
1 (1.6) ∀v ∈ H 1 (Ω)n , |v|H 1 (Ω) ≤ C3 kD(v)k2L2 (Ω) + Φ(v)2 2 .

In particular, Proposition 1.2 implies that there exists a constant CΩ , only depending on Ω such that (see (1.9) in [16]), 2 ! 12 Z 2 1 n , (1.7) ∀v ∈ H (Ω) , |v|H 1 (Ω) ≤ CΩ kD(v)kL2 (Ω) + curl v Ω

4

where | · | denotes the Euclidean vector norm. For any non-negative integer m, recall the classical Sobolev space (cf. [1] or [50]) n o H m (Ω) = v ∈ L2 (Ω) ; ∂ k v ∈ L2 (Ω) ∀ |k| ≤ m ,

equipped with the following seminorm and norm (for which it is a Hilbert space) 1  1  2 2 X X Z 2 k 2     |v|H k (Ω) , kvkH m (Ω) = . |∂ v| dx |v|H m (Ω) = |k|=m Ω

0≤|k|≤m

This definition is extended to any real number s = m + s′ for an integer m ≥ 0 and 0 < s′ < 1 by defining in dimension n the fractional semi-norm and norm  1 2  1 X Z Z |∂ k v(x) − ∂ k v(y)|2 2 2 2  s |v|H s (Ω) =  , kvk = kvk + |v| dx dy . H (Ω) H m (Ω) H s (Ω) n+2 s′ |x − y| Ω Ω |k|=m

1

1

2 (Γ) for a Lipschitz In particular, we shall frequently use the fractional Sobolev spaces H 2 (Γ) and H00 surface Γ when n = 3 or curve when n = 2 with the following seminorms and norms: !1  1 Z Z 2 2 |v(x) − v(y)|2 2 2 , kvk 12 (1.8) |v| 12 dxdy = = kvkL2 (Γ) + |v| 1 , n H (Γ) H (Γ) H 2 (Γ) |x − y| Γ Γ

(1.9)

|v|

1 2 (Γ) H00

=

|v|2 1 H 2 (Γ)

+

Z

Γ

|v(x)|2 dx d∂Γ (x)

!1

2

, kvk

1 2 (Γ) H00

=



kvk2L2 (Γ)

+ |v|

2

1 2 (Γ) H00

1

2

,

where d∂Γ (x) denotes the distance from x to ∂Γ. When Γ is a subset of ∂Ω with positive n − 1 1

1 2 measure, H00 (Ω). The above norms (1.8) and (Γ) is the space of traces of all functions of H0,∂Ω\Γ 1

(1.9) are not equivalent except when Γ is a closed surface or curve. The dual space of H 2 (Γ) is 1 denoted by H − 2 (Γ). In addition to these spaces, we shall use the Hilbert space  H(div; Ω) = v ∈ L2 (Ω)n ; div v ∈ L2 (Ω) , equipped with the graph norm

 1 2 kvkH(div;Ω) = kvk2L2 (Ω) + kdiv vk2L2 (Ω) .

1

The normal trace v · n of a function v of H(div; Ω) on ∂Ω belongs to H − 2 (∂Ω) (cf. [34]). The same result holds when Γ is a part of ∂Ω and is a closed surface. But if Γ is not a closed surface, then 1

2 v · n belongs to the dual of H00 (Γ). When v · n = 0 on ∂Ω, we use the space

H0 (div; Ω) = {v ∈ H(div; Ω) ; v · n = 0 on ∂Ω} . 2. Problem statement

2.1. Coupled Stokes and Darcy systems. Let Ω be partitioned into two non-overlapping regions: the region of the Darcy flow, Ωd , and the region of the Stokes flow, Ωs , each one possibly non-connected, but with a finite number of connected components, and with Lipschitz-continuous boundaries ∂Ωd and ∂Ωs : Ω = Ωd ∪ Ωs . Let Γd = ∂Ωd ∩ ∂Ω, Γs = ∂Ωs ∩ ∂Ω, Γsd = ∂Ωd ∩ ∂Ωs . The unit normal vector on Γsd exterior to Ωd , respectively Ωs , is denoted by nd , respectively ns . In dimension three, we assume that Γd , Γs ,

5

and Γsd also have Lipschitz-continuous boundaries. Let f d be the gravity force in Ωd , f s a given body force in Ωs , let νd > 0, respectively νs > 0, be the constant viscosity coefficient of the Darcy, respectively Stokes flow, let K be the rock permeability tensor in Ωd , let qd be an external source or sink term in Ωd , and let α > 0 be the slip coefficient in the Beavers-Joseph-Saffman law [10, 55], determined by experiment. As far as the data are concerned, we assume on one hand that K is bounded, symmetric and uniformly positive definite in Ωd : there exist two constants, λmin > 0 and λmax > 0 such that ∀x ∈ Ωd , ∀χ ∈ IRn , λmin |χ|2 ≤ K(x)χ · χ ≤ λmax |χ|2 ,

(2.1)

and we assume on the other hand, that the source qd satisfies the solvability condition Z (2.2) qd dx = 0. Ωd

The fluid velocity and pressure in Ωd , respectively Ωs , are denoted by ud and pd , respectively by us and ps . The stress tensor of the Stokes flow is denoted by σ(us , ps ), σ(us , ps ) = −ps I + 2νs D(us ). In the Darcy region Ωd , the pair (ud , pd ) satisfies νd K −1 ud + ∇ pd = f d in Ωd , div ud = qd in Ωd , ud · n = 0 on Γd .

(2.3) (2.4) (2.5)

In the Stokes region Ωs , the pair (us , ps ) satisfies −div σ(us , ps ) ≡ −2νs div D(us ) + ∇ ps = f s in Ωs , div us = 0 in Ωs , us = 0 on Γs .

(2.6) (2.7) (2.8)

These operators suggest to look for (ud , pd ) in H(div; Ωd )×H 1 (Ωd ) and (us , ps ) in H 1 (Ωs )n ×L2 (Ωs ). The Darcy and Stokes flows are coupled on Γsd through the following interface conditions us · ns + ud · nd = 0 on Γsd ,

(2.9) (2.10) √



− σ(us , ps )ns · ns ≡ ps − 2νs D(us )ns · ns = pd

on Γsd ,

√

Kl Kl (2.11) − σ(us , ps )ns · τ l ≡ − 2 D(us )ns · τ l = us · τ l , on Γsd , 1 ≤ l ≤ n − 1, νs α α
where τ l , 1 ≤ l ≤ n − 1 is an orthogonal system of unit tangent vectors on Γsd and Kl = Kτ l · τ l . Conditions (2.9) and (2.10) incorporate continuity of flux and normal stress, respectively. Condition (2.11) is known as the Beavers-Joseph-Saffman law [10, 55, 42] describing slip with friction, where √ Kl /α is a friction coefficient. 2.2. First variational formulation. For any functions ϕd defined in Ωd and ϕs defined in Ωs , it is convenient to define ϕ in Ω by ϕ|Ωd = ϕd and ϕ|Ωs = ϕs . With this notation, regarding the data, we assume that f ∈ L2 (Ω)n , we extend qd by zero in Ωs , i.e. we set qs = 0 and owing to (2.2), the extended function q belongs to L20 (Ω). Before setting problem (2.3)–(2.11) into an equivalent variational formulation, it is useful to interpret the interface conditions (2.10)–(2.11). First we observe from the regularity of f s that each row of σ(us , ps ) belongs to H(div; Ωs ); hence 1 σ(us , ps )ns belongs to H − 2 (∂Ωs )n , and in particular is well-defined as an element of the dual 1

2 of H00 (Γsd )n , which is a distribution space on Γsd . But without further information, it cannot

6

be multiplied directly by the normal or tangent vectors, since the boundary is only Lipschitzcontinuous. To bypass this difficulty, following [35] we define the function on Γsd g = −pd ns −

(2.12)

and replace (2.10)–(2.11) by the condition (2.13)

n−1 X l=1

ν α √s (us · τ l )τ l , Kl

σ(us , ps )ns = g

on Γsd . 1

As the traces of pd and of all components of us on Γsd belong to H 2 (Γsd ), Sobolev’s imbeddings [1] imply that g belongs to Lr (Γsd )n for any finite r when n = 2 and r = 4 when n = 3. Hence condition (2.13) makes sense. Let us check that it implies (2.10)–(2.11). First, prescribing (2.13) guarantees that σ(us , ps )ns belongs at least to L4 (Γsd )n and thus can be multiplied by the normal or tangent vectors. Then by virtue of this regularity, (2.12), (2.13) are equivalent to:





σ(us , ps )ns · ns ns +

n−1 X l=1







σ(us , ps )ns · τ l τ l = −pd ns −

n−1 X l=1

ν α √s (us · τ l )τ l , Kl

and therefore, by identifying on both sides the components of the normal and tangent vectors (that forms an orthonormal set), we derive (2.10)–(2.11). Hence (2.13) is the interpretation of (2.10)–(2.11). Now, let (ud , pd ) ∈ H(div; Ωd ) × H 1 (Ωd ) and (us , ps ) ∈ H 1 (Ωs )n × L2 (Ωs ) be a solution of (2.3)–(2.11). In order to set (2.3)–(2.11) in variational form, we take the scalar product of (2.3) and (2.6) respectively with any test function v d in H 1 (Ωd )n satisfying v d · n = 0 on Γd , and any v s in H 1 (Ωs )n satisfying v s = 0 on Γs . Then we apply Green’s formula in Ωd and Ωs . This yields Z Z Z Z −1 (2.14) νd K ud · v d − pd div v d + pd v d · nd = f d · vd, Ωd

(2.15)

2 νs

Z

Ωs

Ωd

D(us ) : D(v s ) −

Z

Ωs

Γsd

Ωd

ps div v s − hσ(us , ps )ns , v s iΓsd = 1

Z

Ωs

f s · vs,

2 (Γsd )n and its dual space. The validity of where h·, ·iΓsd denotes the duality pairing between H00 (2.14) and (2.15) follows from the above considerations. By summing (2.14) and (2.15), by using the fact that nd = −ns , and by applying (2.13), the term on the interface, say I, reads Z Z Z I = −hσ(us , ps )ns , v s iΓsd − pd v d · ns = − g · vs − pd v d · ns .

Γsd

Γsd

Γsd

Then the expression (2.12) for g yields Z n−1 XZ νs α √ (us · τ l )(v s · τ l ) + pd [v · n], (2.16) I= Kl Γsd l=1 Γsd where the jump [v · n] is defined by

[v · n] = v s · ns + v d · nd .

Finally, we eliminate this jump by enforcing strongly the transmission condition (2.9) on the test function v. In view of the interior terms in (2.14) and (2.15) and what remains in (2.16), we see that we can reduce the regularity of our functions and work in the space ˜ = {v ∈ H(div; Ω) ; v s ∈ H 1 (Ωs )n , v|Γs = 0, (v · n)|Γ = 0}, (2.17) X d

which is a Hilbert space equipped with the norm (2.18)

2 ˜ , kvk ˜ = kvk2 ∀v ∈ X H(div;Ω) + |v s |H 1 (Ωs ) X


1 2

.

7

Note that the restriction of v · n on Γsd belongs at least to L4 (Γsd ). Let us denote W = L20 (Ω) with the norm kwkW = kwkL2 (Ω) . Then we propose the following variational formulation : Find ˜ × W such that (u, p) ∈ X Z Z Z ˜ , νd ∀v ∈ X p div v K −1 ud · v d + 2 νs D(us ) : D(v s ) − Ωd

(2.19)

Ωs

+

n−1 XZ l=1

ν α √s (us · τ l )(v s · τ l ) = Kl

Γsd

∀w ∈ W ,

(2.20)

Ω

Z

w div u =

Z

Z

f · v,

Ω

w qd .

Ωd

Ω

˜ × W of problem Lemma 2.1. For any data f in L2 (Ω)n and qd in L20 (Ωd ), any solution (u, p) ∈ X 1 (2.19)–(2.20) satisfies pd ∈ H (Ωd ) and solves (2.3)–(2.11). Conversely, any solution of (2.3)– ˜ × W with pd ∈ H 1 (Ωd ), solves (2.19)–(2.20). (2.11), (u, p) ∈ X Proof. It stems from the above considerations that any solution (ud , pd ) ∈ H(div; Ωd ) × H 1 (Ωd ) ˜ × L2 (Ω) and solves and (us , ps ) ∈ H 1 (Ωs )n × L2 (Ωs ) of (2.3)–(2.11) is such that (u, p) belongs to X (2.19)–(2.20). Moreover as Ω is connected, (2.19) only defines p up to an additive constant and this constant can be chosen so that p belongs to W . ˜ × W solve (2.19)–(2.20), and denote its restriction to Ωd and Ωs as Conversely, let (u, p) ∈ X above. By choosing smooth test functions with compact support first in Ωd and next in Ωs , we immediately derive that (ud , pd ) is a solution of (2.3)–(2.5) and (us , ps ) is a solution of (2.6)–(2.8). Furthermore, since both f d and νd K −1 ud belong to L2 (Ωd )n , (2.3) implies that pd ∈ H 1 (Ωd ). It remains to recover the transmission conditions (2.9)–(2.11). First, (2.9) is a consequence of ˜ Next, we recover (2.14) and (2.15) by taking the scalar product of (2.3) the definition (2.17) of X. ˜ that is smooth in Ωd and in Ωs , and by applying Green’s formula and (2.6) with a function v ∈ X in both regions. By comparing with (2.19), this gives Z

Γsd

pd v d · nd − hσ(us , ps )ns , v s iΓsd =

n−1 XZ l=1

Γsd

ν α √s (us · τ l )(v s · τ l ). Kl

By taking into account the orientation of the normal, the regularity of v, and the definition (2.12) of g, this is equivalent to: Z Z n−1 XZ νs α √ (us · τ l )(v s · τ l ) = g · vs. hσ(us , ps )ns , v s iΓsd = − pd v s · ns − Kl Γsd Γsd l=1 Γsd As the trace space of v s on Γsd is large enough, this implies (2.13).



2.3. Existence and uniqueness of the solution. For any functions ud , v d in L2 (Ωd )n and us , v s in H 1 (Ωs )n , we define the bilinear form (2.21)

a ˜(u, v) = νd

Z

Ωd

K −1 ud · v d + 2 νs

Z

D(us ) : D(v s ) +

Ωs

n−1 XZ l=1

Γsd

ν α √s (us · τ l )(v s · τ l ), Kl

and for any functions v d ∈ H(div; Ωd ), v s ∈ H(div; Ωs ) and w ∈ L2 (Ω), we define the bilinear form Z Z ˜ (2.22) b(v, w) = − w div v d − w div v s . Ωd

Ωs

8

˜ ×X ˜ and ˜b(·, ·) is continuous on X ˜ × L2 (Ω): Note that a ˜(·, ·) is continuous on X (2.23)

˜ ×X ˜ , |˜ ∀(u, v) ∈ X a(u, v)| ≤

νd kud kL2 (Ωd ) kv d kL2 (Ωd ) + 2 νs k∇ us kL2 (Ωs ) k∇ v s kL2 (Ωs ) λmin n−1 X νs α √ kus · τ l kL2 (Γsd ) kv s · τ l kL2 (Γsd ) , + λmin l=1

˜ × L2 (Ω) , |˜b(v, w)| ≤kvk ˜ kwkL2 (Ω) . ∀(v, w) ∈ X X

˜ × W such that Then (2.19)–(2.20) has the familiar form : Find (u, p) ∈ X Z ˜ ˜ f · v, (2.24) ∀v ∈ X , a ˜(u, v) + b(v, p) = Ω

(2.25) Next, we set (2.26)

∀w ∈ W , ˜b(u, w) = −

Z

w qd .

Ωd

˜ 0 = {v ∈ X ˜ ; div v = 0}, X

and more generally, for a given function g ∈ W , we define the affine variety ˜ g = {v ∈ X ˜ ; div v = g}. (2.27) X ˜ q such that Then we consider the reduced problem : Find u ∈ X Z ˜0 , a f · v, (2.28) ∀v ∈ X ˜(u, v) = Ω

recall that q is qd extended by zero on Ωs . It is well known [34, 18] that showing equivalence between problems (2.28) and (2.24)–(2.25) amounts to proving the following inf-sup condition. Lemma 2.2. There exists a constant β > 0 such that ˜b(v, w) (2.29) ∀w ∈ W , sup ≥ βkwkW . ˜ ˜ kvkX v∈X Proof. Let w ∈ W . The inf-sup condition between H01 (Ω)n and L20 (Ω) implies that there exists a function v ∈ H01 (Ω)n such that div v = w in Ω and |v|H 1 (Ω) ≤

1 kwkL2 (Ω) , κ

˜ and a straightforward where κ depends only on Ω; see for example [34] or [18]. Then v belongs to X 1
2 + 2 2 , where P e constant in (1.1).  computation yields (2.29) with β ≥ κ/ P∂Ω ∂Ω is the Poincar´ Lemma 2.2 implies that (2.28) and (2.24)–(2.25) are equivalent in the following sense.

˜ × W be a solution of Proposition 2.1. Let (f , qd ) be given in L2 (Ω)n × L20 (Ωd ). Let (u, p) ∈ X ˜ (2.24)–(2.25). Then u solves (2.28). Conversely, let u ∈ Xqd be a solution of (2.28). Then there exists a unique p in W such that (u, p) satisfies (2.24)–(2.25). Now, let us prove that (2.28), and hence (2.24)–(2.25), is well-posed. This relies on the ellipticity ˜ 0 . If Ωs is connected and |Γs | > 0, a partial ellipticity result for a of a ˜(·, ·) on X ˜(·, ·) follows directly from Korn’s inequality (1.4): νd νs ˜,a kv d k2L2 (Ωd ) , (2.30) ∀v ∈ X ˜(v, v) ≥ 2 2 |v s |2H 1 (Ωs ) + λmax C1

9

with the constant C1 of (1.4). If Ωs is connected and |Γs | = 0, then Γsd = ∂Ωs up to a set of zero measure, and proving the analogue of (2.30) makes use of (1.7) and the tangential components on Γsd . Indeed, we have formally X n−1 |(v s · τ l )(x)|, a.e. on ∂Ωs , (v s × ns )(x) ≤

(2.31)

l=1

and therefore Z

Ωs

curl v s = −

Z

v s × ns

Γsd

Z ⇒

Ωs

Z curl v s ≤

Γsd

Hence (1.7) and a straightforward manipulation yield for all v s in (2.32)

2νs kD(v s )k2L2 (Ωs )

+

n−1 X l=1

 n−1 X l=1

 |v s · τ l | .

H 1 (Ωs )n :


α ν α ν √s kv s · τ l k2L2 (Γsd ) ≥ s2 min 2, √ |v s |2H 1 (Ωs ) . CΩ Kl λmax |Γsd |

As a consequence (2.30) is replaced by
νd νs α ˜,a |v s |2H 1 (Ωs ) . (2.33) ∀v ∈ X ˜(v, v) ≥ kv d k2L2 (Ωd ) + 2 min 2, √ λmax CΩ λmax |Γsd |

Finally, if Ωs is not connected, the analogue of (2.30) holds on all connected components of Ωs that are adjacent to Γs and the analogue of (2.33) holds on all connected components of Ωs that are not adjacent to Γs . It remains to establish that the mapping:
1 (2.34) v 7→ |v|X˜ 0 = |v s |2H 1 (Ωs ) + kv d k2L2 (Ωd ) 2 ˜ 0 equivalent to kvk ˜ . This is the object of the next lemma. is a norm on X X

Lemma 2.3. There exists a constant C4 such that 2 ˜ 0 , kv s kL2 (Ω ) ≤ C4 |v s |2 1 ∀v ∈ X H (Ωs ) + kv d kL2 (Ωd ) s

(2.35)


1

2

.

Proof. Let us assume that Ωs is connected; the case when Ωs is not connected is treated as above. If |Γs | > 0, (2.35) follows from Poincar´e’s inequality (1.1) applied in Ωs and does not require the norm of v d in the right-hand side. When |Γs | = 0, the proof of (2.35) is a variant of the proof of Peetre-Tartar’s Lemma [51]. Let us recall its argument. Assume that (2.35) is not true. Then ˜ 0 such that there exists a sequence (v m ) in X m m lim kv m d kL2 (Ωd ) = lim |v s |H 1 (Ωs ) = 0 and ∀m , kv s kL2 (Ωs ) = 1.

m→∞

m→∞

˜ 0 is reflexive, this implies that there exists a function v ∈ X ˜ 0 such that As X ˜ lim v m = v weakly in X. m→∞

˜ implies that Moreover v s = c, a constant vector, and v d = 0. Then the fact that v belongs to X c · ns = 0 on Γsd , and since Γsd coincides a.e. with ∂Ωs , this implies that c = 0. Thus 1 n lim v m s = 0 weakly in H (Ωs ) ,

m→∞

hence 2 n lim v m s = 0 strongly in L (Ωs ) .

m→∞

This contradicts the fact that for all m kv m s kL2 (Ωs ) = 1. Therefore (2.35) combined with either (2.30) or (2.33) yields the ellipticity of a ˜(·, ·).



10

Ωd,4

Ωd,5

Ωd,3

Ωd,6

Ωs,2 Ωs,3

Ωs,1

Ωd,2 Ωd,1

Ωs,4

Figure 2.1. Domain decomposition with multiple connected components. Lemma 2.4. There exists a constant C5 > 0 such that ˜0 , a (2.36) ∀v ∈ X ˜(v, v) ≥ C5 kvk2˜ . X

˜ 0 , and the inf-sup condition With the continuity of a ˜(·, ·) and ˜b(·, ·), the ellipticity of a ˜(·, ·) on X (2.29), the Babuˇska-Brezzi’s theory [7, 17] implies immediately that (2.24)–(2.25) is well-posed. ˜ × W and there exists a Theorem 2.1. Problem (2.24)–(2.25) has a unique solution (u, p) ∈ X constant C that depends only on Ωd , Ωs , λmin , λmax , α, νd , and νs , such that
(2.37) kukX˜ + kpkL2 (Ω) ≤ C kf kL2 (Ω) + kqd kL2 (Ωd ) . In turn, the well-posedness of problem (2.3)–(2.11) stems from Lemma 2.1.

2.4. Domain decomposition of the Darcy and Stokes regions. Let Ωs , respectively Ωd , be decomposed into Ms , respectively Md , non-overlapping, open Lipschitz subdomains: Md s Ωs = ∪M i=1 Ωs,i , Ωd = ∪i=1 Ωd,i .

Set M = Md + Ms ; according to convenience we can also number the subdomains with a single index i, 1 ≤ i ≤ M , the Darcy subdomains running from Ms + 1 to M . An example of a domain decomposition with multiple connected components is depicted in Figure 2.1. Let ni denote the outward unit normal vector on ∂Ωi . For 1 ≤ i ≤ M , let the boundary interfaces be denoted by Γi , with possibly zero measure: Γi = ∂Ωi ∩ ∂Ω, and for 1 ≤ i < j ≤ M , let the interfaces between subdomains be denoted by Γij , again with possibly zero measure: Γij = ∂Ωi ∩ ∂Ωj . In addition to Γsd , let Γdd , respectively Γss , denote the set of interfaces between subdomains of Ωd , respectively Ωs . Then, assuming that the solution (u, p) of (2.3)–(2.11) is slightly smoother, we can obtain an equivalent formulation by writing individually (2.3)–(2.11) in each subdomain Ωi , 1 ≤ i ≤ M , and complementing these systems with the following interface conditions (2.38)

(2.39)

[ud · n] = 0 ,

[us ] = 0 ,

[pd ] = 0 on Γdd ,

[σ(us , ps )n] = 0 on Γss ,

where the jumps on an interface Γij , 1 ≤ i < j ≤ M , are defined as [v · n] = v i · ni + v j · nj ,

[σn] = σ i ni + σ j nj ,

[v] = (vi − vj )|Γij ,

using the notation vi = v|Ωi . The smoothness requirement on the solution is meant to ensure that the jumps [ud · n], respectively [σ(us , ps )n], are well-defined on each interface of Γdd , respectively Γss .

11

Finally, let us prescribe weakly the interface conditions (2.38), (2.39), and (2.9) by means of Lagrange multipliers, usually called mortars. For this, it is convenient to attribute a unit normal vector nij to each interface Γij of positive measure, directed from Ωi to Ωj (recall that i < j). The basic velocity spaces are: Xd = {v ∈ L2 (Ωd )n ; v d,i := v|Ωd,i ∈ H(div; Ωd,i ), 1 ≤ i ≤ Md , 1

v · nij ∈ H − 2 (Γij ), Γij ∈ Γdd ∪ Γsd , v · n = 0 on Γd },

(2.40)

Xs = {v ∈ L2 (Ωs )n ; v s,i := v|Ωs,i ∈ H 1 (Ωs,i )n , 1 ≤ i ≤ Ms , v = 0 on Γs }, and the mortar spaces are: (2.41)


n 1 ∀Γij ∈ Γss , Λij = H − 2 (Γij )

1

∀Γij ∈ Γsd ∪ Γdd , Λij = H 2 (Γij ).

,

˜ (see (2.17)) by Then we replace X

X = {v ∈ L2 (Ω)n ; v d := v|Ωd ∈ Xd , v s := v|Ωs ∈ Xs },

(2.42)

we keep W = L20 (Ω) for the pressure, and we define the mortar spaces
n
n 1 Λs = {λ ∈ D′ (Γss ) ; λ|Γij ∈ H − 2 (Γij ) for all Γij ∈ Γss }, 1

Λsd = {λ ∈ L2 (Γsd ) ; λ|Γij ∈ H 2 (Γij ) for all Γij ∈ Γsd },

(2.43)

1

Λd = {λ ∈ L2 (Γdd ) ; λ|Γij ∈ H 2 (Γij ) for all Γij ∈ Γdd }.

We equip these spaces with broken norms: |||v|||Xd =

Md X i=1

kvk2H(div;Ωd,i )

|||λ|||Λs =

1

 X

2

Γij ∈Γss

, |||v|||Xs =

kλk2

i=1

1

2

1 H − 2 (Γij )

|||λ|||Λd =

Ms X

kvk2H 1 (Ωs,i )

, |||λ|||Λsd =

 X

Γij ∈Γdd

kλk2

H

1  2 , |||v|||X = |||v|||2Xd + |||v|||2Xs ,

1

2

 X

Γij ∈Γsd

1

2

1 2 (Γij )

kλk2

H

1

2

1 2 (Γij )

,

.

Note that in most geometrical situations, Xd (and hence X) is not complete for the above norm, but none of the subsequent proofs require its completeness. The matching condition between subdomains is weakly enforced through the following bilinear forms: X ∀v ∈ Xs , ∀µ ∈ Λs , bs (v, µ) = h[v], µiΓij , Γij ∈Γss

(2.44)

∀v ∈ Xd , ∀µ ∈ Λd , bd (v, µ) = ∀v ∈ X, ∀µ ∈ Λsd , bsd (v, µ) =

X

Γij ∈Γdd

X

Γij ∈Γsd

h[v · n], µiΓij , h[v · n], µiΓij .

12

For the velocity and pressure in the Darcy and Stokes regions, we use the following bilinear forms: Z ∀(u, v) ∈ Xs × Xs , as,i (u, v) = 2 νs D(us,i ) : D(v s,i ) Ωs,i

+ (2.45)

n−1 XZ

∂Ωs,i ∩Γsd

l=1

∀(u, v) ∈ Xd × Xd , ad,i (u, v) = νd 2

∀v ∈ X, ∀w ∈ L (Ω) , bi (v, w) = − Then we set

Z

Ωd,i

Z

Ωi

∀(u, v) ∈ X × X , a(u, v) =

ν α √s (us · τ l )(v s · τ l ) , 1 ≤ i ≤ Ms , Kl

K −1 ud,i · v d,i , 1 ≤ i ≤ Md ,

wdiv v i , 1 ≤ i ≤ M.

Ms X

as,i (u, v) +

i=1

2

Md X

ad,i (u, v),

i=1

∀(v, w) ∈ X × L (Ω) , b(v, w) =

M X

bi (v, w).

i=1

The second variational formulation reads: Find (u, p, λsd , λd , λs ) ∈ X × W × Λsd × Λd × Λs such that Z f · v, ∀v ∈ X , a(u, v) + b(v, p) + bsd (v, λsd ) + bd (v, λd ) + bs (v, λs ) = Ω Z (2.46) w qd , ∀w ∈ W , b(u, w) = − Ωd

∀µ ∈ Λsd , bsd (u, µ) = 0 , ∀µ ∈ Λd , bd (u, µ) = 0 , ∀µ ∈ Λs , bs (u, µ) = 0.

It remains to prove that (2.46) is equivalent to (2.3)–(2.11) when the solution is sufficiently smooth. Since we know from Theorem 2.1 that (2.3)–(2.11) has a unique solution, equivalence will also establish that (2.46) is uniquely solvable. ˜ ×W of (2.3)–(2.11) with pd ∈ H 1 (Ωd ) satisfies Theorem 2.2. Assume that the solution (u, p) ∈ X 1

1

∀Γij ∈ Γdd ∪ Γsd , (ud · nd )|Γij ∈ H − 2 (Γij ) , ∀Γij ∈ Γss , (σ(us , ps )ns )|Γij ∈ H − 2 (Γij )n .

Then (2.46) is equivalent to (2.3)–(2.11).

Proof. The argument is similar to that used in proving Lemma 2.1. Let (u, p) be a solution of (2.3)–(2.11) satisfying the above regularity. Take the scalar product in each Ωi of (2.3) and (2.6) with a test function v in X, apply Green’s formula and add the corresponding equations. In view of (2.16) and the regularity of (u, p), this gives: Z X X (2.47) a(u, v) + b(v, p) − f · v. hσ(us , ps )nij , [v]iΓij + hpd , [v · n]iΓij = Γij ∈Γss

Γij ∈Γdd ∪Γsd

Ω

We recover the first equation in (2.46) by defining (2.48)

∀Γij ∈ Γss , λs |Γij = −σ(us , ps )|Γij nij ,

∀Γij ∈ Γdd , λd |Γij = pd |Γij

,

∀Γij ∈ Γsd , λsd |Γij = pd |Γij ,

and the remaining equations follow from the regularity of (u, p). Conversely, let (u, p, λsd , λd , λs ) in X × W × Λsd × Λd × Λs be a solution of (2.46). By choosing smooth test functions with compact support in each Ωi we recover the interior equations (2.3), (2.4), (2.6), (2.7) in each subdomain. On one hand, we easily derive from the last equation of (2.46) that u has no jump through the interfaces of Γss . Hence u ∈ H 1 (Ωs )n . On the other hand,

13

we pick an index i with 1 ≤ i ≤ Ms and an interface Γij in Γss , we take a function v in X, smooth in Ωi , zero outside Ωi , and zero on ∂Ωi \ Γij . By taking the scalar product of (2.6) in Ωi with v, applying Green’s formula, comparing with (2.46), and using the same process in Ωj , we find λs |Γij = −σ(us,i , ps,i )|Γij nij = −σ(us,j , ps,j )|Γij nij .

As λs is single-valued, this implies that σ(us , ps )|Γij nij has no jump through Γij . This is true for all interfaces in Γss . Therefore (2.6) is satisfied in Ωs . By applying a similar process to the interfaces of Γdd , we derive first that (ud , pd ) belongs to H(div; Ωd ) × H 1 (Ωd ) and next that (2.3) is satisfied in Ωd . Finally, the third equation in (2.46) implies that u belongs to H(div; Ω), therefore u is in X and the pair (u, p) solves (2.19)–(2.20); by virtue of Lemma 2.1, it also solves (2.3)–(2.11).  3. Discretization 3.1. Meshes and discrete spaces. In view of discretization, we assume from now on that Ω and all its subdomains Ωi , 1 ≤ i ≤ M , have polygonal or polyhedral boundaries. Since none of the subdomains overlap, they form a mesh, Td of Ωd and Ts of Ωs , and the union of these meshes constitutes a mesh TΩ of Ω. Furthermore, we suppose that this mesh satisfies the following assumptions: Hypothesis 3.1. (1) TΩ is conforming, i.e. it has no hanging nodes. (2) The subdomains of TΩ can take at most L different geometrical shapes, where L is a fixed integer independent of M . (3) TΩ is shape-regular in the sense that there exists a real number σ, independent of M such that diam(Ωi ) ≤ σ, (3.1) ∀i, 1 ≤ i ≤ M , diam(Bi ) where diam(Ωi ) is the diameter of Ωi and diam(Bi ) is the diameter of the largest ball contained in Ωi . Without loss of generality, we can assume that diam(Ωi ) ≤ 1.

As each subdomain Ωi is polygonal or polyhedral, it can be entirely triangulated. Let h > 0 denote a discretization parameter, and for each h, let Tih be a regular family of partitions of Ωi made of triangles or tetrahedra T in the Stokes region and triangles, tetrahedra, parallelograms, or parallelepipeds in the Darcy region, with no matching requirement at the subdomains interfaces. Thus the meshes are independent and the parameter h < 1 is allowed to vary with i, but to reduce the notation, unless necessary, we do not indicate its dependence on i. By regular, we mean that there exists a real number σ0 , independent of i and h such that hT ≤ σ0 , (3.2) ∀i, 1 ≤ i ≤ M, ∀T ∈ Tih , ρT where hT is the diameter of T and ρT is the diameter of the largest ball contained in T . In addition we assume that each element of Tih has at least one vertex in Ωi . For the interfaces, let H > 0 be another discretization parameter and for each H and each i < j, let TijH denote a regular family of partitions of Γij into segments, triangles or parallelograms of diameter bounded by H. These partitions may not match the traces of the subdomain grids. On these meshes, we define the following finite element spaces. In the Stokes region, for each h , W h ) ⊂ H 1 (Ω )n × L2 (Ω ) be a pair of finite element spaces satisfying a local Ωs,i , let (Xs,i s,i s,i s,i h ∩ H 1 (Ω )n and h = Xs,i uniform inf-sup condition for the divergence. More precisely, setting X0,s,i s,i 0 ⋆ h h 2 W0,s,i = Ws,i ∩ L0 (Ωs,i ), we assume that there exists a constant βs > 0, independent of h and the diameter of Ωs,i , such that R h h Ωs,i w div v ≥ βs⋆ . sup (3.3) ∀i, 1 ≤ i ≤ Ms , inf h| 1 hk 2 h |v kw h h wh ∈W0,s,i H (Ωs,i ) L (Ωs,i ) v ∈X 0,s,i

14 h In addition, since X0,s,i ⊂ H01 (Ωs,i )n , it satisfies a Korn inequality: There exists a constant α⋆ > 0, independent of h and the diameter of Ωs,i , such that

(3.4)

h ∀i, 1 ≤ i ≤ Ms , ∀v h ∈ X0,s,i , kD(v h )kL2 (Ωs,i ) ≥ α⋆ |v h |H 1 (Ωs,i ) .

There are well-known examples of pairs satisfying (3.3) (cf. [34]), such as the mini-element, the Bernardi-Raugel element, or the Taylor-Hood element. Similarly, in the Darcy region, for each h , W h ) ⊂ H(div; Ω ) × L2 (Ω ) be a pair of mixed finite element spaces satisfying a Ωd,i , let (Xd,i d,i d,i d,i h h ∩ H (div; Ω ) uniform inf-sup condition for the divergence. More precisely, setting X0,d,i = Xd,i 0 d,i ⋆ h h 2 and W0,d,i = Wd,i ∩ L0 (Ωd,i ), we assume that there exists a constant βd > 0 independent of h and the diameter of Ωd,i , such that R h h Ωd,i w div v ≥ βd⋆ . sup (3.5) ∀i, 1 ≤ i ≤ Md , inf hk hk 2 h kv kw h h wh ∈W0,d,i H(div;Ωd,i ) L (Ωd,i ) v ∈X 0,d,i

Furthermore, we assume that (3.6)

h h ∀i, 1 ≤ i ≤ Md , ∀v h ∈ Xd,i , div v h ∈ Wd,i .

Again, there are well-known examples of pairs satisfying (3.5) and (3.6) (cf. [18] or [34]), such as the Raviart-Thomas elements, the Brezzi-Douglas-Marini elements, the Brezzi-Douglas-FortinMarini elements, the Brezzi-Douglas-Dur´ an-Fortin elements, or the Chen-Douglas elements. Then we discretize straightforwardly Xd and Xs by h Xdh = {v ∈ L2 (Ωd )n ; v|Ωd,i ∈ Xd,i , 1 ≤ i ≤ Md , v · n = 0 on Γd },

and we set

h Xsh = {v ∈ L2 (Ωs )n ; v|Ωs,i ∈ Xs,i , 1 ≤ i ≤ Ms , v = 0 on Γs },

h h Wdh = {w ∈ L2 (Ωd ) ; w|Ωd,i ∈ Wd,i } , Wsh = {w ∈ L2 (Ωs ) ; w|Ωs,i ∈ Ws,i },

X h = {v ∈ L2 (Ω)n ; v|Ωd ∈ Xdh , v|Ωs ∈ Xsh } , W h = {w ∈ L20 (Ω) ; w|Ωd ∈ Wdh , w|Ωs ∈ Wsh }. The finite elements regularity implies that Xdh ⊂ Xd , Xsh ⊂ Xs and X h ⊂ X. Of course, W h ⊂ W . H H In the mortar region, we take finite element spaces ΛH s , Λd , and Λsd . These spaces consist of continuous or discontinuous piecewise polynomials. We will allow for varying polynomial degrees on different types of interfaces. Although the mortar meshes and the subdomain meshes so far are H H h unrelated, we need compatibility conditions between ΛH s , Λsd and Λd on one hand, and Xd and h Xs on the other hand. The following set of conditions is fairly crude and will be sharpened further on. ˜ there exists v h ∈ X h , v h = 0 on ∂Ωs,i \ Γij (1) For all Γij ∈ Γss ∪ Γsd , i < j, and for all v ∈ X, s,i satisfying Z Z (3.7) v h · nij = v · nij . Γij

Γij

˜ there exists v h ∈ X h , v h = 0 on ∂Ωs,j \ Γij (2) For all Γij ∈ Γss , i < j, and for all v ∈ X, s,j satisfying Z Z H h µH · v. µ · v = (3.8) ∀µH ∈ ΛH , s Γij

Γij

˜ there exists v h ∈ X h , v h · nj = 0 on (3) For all Γij ∈ Γdd ∪ Γsd , i < j, and for all v ∈ X, d,j ∂Ωd,j \ Γij satisfying Z Z H h H H H H µH v · nij . µ v · nij = (3.9) ∀µ ∈ Λd , ∀µ ∈ Λsd , Γij

Γij

15

Condition (3.7) is very easy to satisfy in practice and it trivially holds true for all examples of Stokes spaces considered in this paper. Conditions (3.8) and (3.9) state that the mortar space is controlled by the traces of the subdomain velocity spaces. Both conditions are easier to satisfy for coarser mortar grids. Condition (3.8) is more general than previously considered in the literature for mortar discretizations of the Stokes equations [11, 12]. Examples for which (3.8) holds are given in the Appendix: Sections 7.1 in 2-D and 7.2 in 3-D. The condition (3.9) on Γdd is closely related to the mortar condition for Darcy flow in [60, 5, 52, 6] and on Γsd , it is more general than existing mortar discretizations for Stokes-Darcy flows on [46, 54, 29, 14]. It is discussed in more detail in Section 4.4.
H H H H in ΛH Lemma 3.1. Under assumptions (3.8) and (3.9), the only solution λH sd ×Λd ×Λs sd , λd , λs to the system (3.10)

is the zero solution.

h H h H ∀v h ∈ X h , bs (v h , λH s ) + bd (v , λd ) + bsd (v , λsd ) = 0

Proof. Consider any Γij ∈ Γss with i < j; the proof for the other interfaces being the same. Take an arbitrary v in H01 (Ω)n and v h associated with v by (3.8), extended by zero outside Ωs,j . Then on one hand, Z Z Γij

and on the other hand,

λH s ·v =

Γij

h h H λH s · v = bs (v , λs ),

h H bd (v h , λH d ) = bsd (v , λsd ) = 0. Z λH ∀v ∈ H01 (Ω)n , s · v = 0,

Therefore

Γij

thus implying that

λH s

= 0.



Finally, we define h h H Vdh = {v ∈ Xdh ; ∀µ ∈ ΛH d , bd (v, µ) = 0} , Vs = {v ∈ Xs ; ∀µ ∈ Λs , bs (v, µ) = 0},

V h = {v ∈ X h ; v|Ωd ∈ Vdh , v|Ωs ∈ Vsh , ∀µ ∈ ΛH sd , bsd (v, µ) = 0},

(3.11)

Z h = {v ∈ V h ; ∀w ∈ W h , b(v, w) = 0}.

3.2. Variational formulations and uniform stability of the discrete problem. The discrete H H h h version of the second variational formulation (2.46) is: Find (uh , ph , λH sd , λd , λs ) ∈ X × W × H H ΛH sd × Λd × Λs such that Z h h h h h h h H h H h H f · vh, ∀v ∈ X , a(u , v ) + b(v , p ) + bsd (v , λsd ) + bd (v , λd ) + bs (v , λs ) = Ω Z (3.12) ∀wh ∈ W h , b(uh , wh ) = − w h qd , Ωd

H

∀µ

∈

ΛH sd ,

h

H

h H H H h H bsd (u , µ ) = 0 , ∀µH ∈ ΛH d , bd (u , µ ) = 0 , ∀µ ∈ Λs , bs (u , µ ) = 0.

The last three equations of (3.12) state that uh ∈ V h . Therefore, we can extract from (3.12) the following reduced formulation: Find uh ∈ V h , ph ∈ W h such that Z h h h h h h f · vh, ∀v ∈ V , a(u , v ) + b(v , p ) = Ω Z (3.13) ∀wh ∈ W h , b(uh , wh ) = − w h qd . Ωd

Lemma 3.2. Problems (3.12) and (3.13) are equivalent.

16

Proof. Clearly, (3.12) implies (3.13). Conversely, if the pair (uh , ph ) solves (3.13), existence of H H λH sd , λd , λs such that all these variables satisfy (3.12) is an easy consequence of Lemma 3.1 and an algebraic argument.  In view of this equivalence, it suffices to analyze problem (3.13). From the Babuˇska–Brezzi’s theory, uniform stability of the solution of (3.13) stems from an ellipticity property of the bilinear form a in Z h and an inf-sup condition of the bilinear form b. Let us prove an ellipticity property of the bilinear form a, valid when n = 2, 3. For this, we make the following assumptions on the mortar spaces: H Hypothesis 3.2. (1) On each Γij ∈ Γdd ∪ Γsd , ΛH d |Γij and Λsd |Γij contain at least IP 0 . n (2) On each Γij ∈ Γss , on each hyperplane F ⊂ Γij , ΛH s |F contains at least IP 0 . n H (3) On each Γij ∈ Γss , Λs |Γij contains at least IP 1 .

The second assumption guarantees that nij ∈ ΛH s |Γij ; it follows from the third assumption when Γij is flat. The third assumption is solely used for deriving a discrete Korn inequality; it can be relaxed, as we shall see in Section 7.2. The first two assumptions imply that all functions v h in V h satisfy M Z M Z X X XZ h h div v = v · ni = [v h · n] = 0. i=1

Ωi

i=1

∂Ωi

i<j

Γij

Therefore, the zero mean-value restriction on the functions of W h can be relaxed. Thus the condition v h ∈ Z h and (3.6) imply in particular that div v hd = 0 in Ωd,i , 1 ≤ i ≤ Md . Hence ∀v h ∈ Z h , |||v hd |||Xd = kv hd kL2 (Ωd ) .

(3.14)

First, we treat the simpler case when |Γs | > 0 and Ωs is connected. Lemma 3.3. Let |Γs | > 0 and Ωs be connected. Then under Hypothesis 3.2, we have ∀v h ∈ Z h , a(v h , v h ) ≥

(3.15)

νs νd |||v hd |||2Xd + 2 2 |||v hs |||2Xs , λmax C

where the constant C only depends on the shape regularity of Ts . h Proof. Since v hs ∈ Vsh and IP n1 ∈ ΛH ss |Γij for each Γij ∈ Γss , then P1 [v s ] = 0, where P1 is the orthogonal projection on IP n1 for the L2 norm on each Γij . Thus, as Ωs is connected and |Γs | > 0, inequality (1.12) in [16] gives

(3.16)

∀v hs

∈

Vsh ,

Ms X i=1

|v hs |2H 1 (Ωs,i )

≤C

2

Ms X i=1

kD(v hs )k2L2 (Ωs,i ) ,

where the constant C only depends on the shape regularity of Ts . Hence we have the analogue of (2.30): (3.17)

Ms νs X νd h 2 |||v ||| + 2 2 |v hs |2H 1 (Ωs,i ) . ∀v ∈ Z , a(v , v ) ≥ λmax d Xd C h

h

h

h

i=1

Finally the above argument permits to apply formula (1.3) in [15] in order to recover the full norm of Xs in the right-hand side of (3.17). In fact, it is enough that IP n0 ∈ ΛH ss |Γij for each Γij ∈ Γss .  Now we turn to the case when Ωs is connected and |Γs | = 0, consequently Γsd = ∂Ωs , up to a set of zero measure.

17

Lemma 3.4. Let |Γs | = 0 and Ωs be connected, i.e. Γsd = ∂Ωs . Then under Hypothesis 3.2, we have
h 2 α νs νd |||v s |||Xs , |||v hd |||2Xd + 2 min 2, √ (3.18) ∀v h ∈ Z h , a(v h , v h ) ≥ λmax C λmax |Γsd | where the constant C only depends on the shape regularity of Ts .

Proof. All constants in this proof only depend on the shape regularity of Ts . Since (3.14) holds, it suffices to derive a lower bound for as (v hs , v hs ). To begin with, as v hs ∈ Vsh , we apply Theorem 4.2 if n = 2 or 5.2 if n = 3 in [16]: ∀v hs

∈

Vsh ,

Ms X i=1

|v hs |2H 1 (Ωs,i )

≤C

2

Ms X i=1


2  kD(v hs )k2L2 (Ωs,i ) + Φ(v hs ) ,

where the functional Φ is a suitable seminorm. Let us choose Z (3.19) ∀v ∈ Xs , Φ(v) = v × ns . Γsd

Clearly, Φ is a seminorm on Xs . Next, considering that Z 1 n ∀v ∈ H (Ωs ) , Φ(v) =

Ωs

curl v ,

it is easy to check that if m is a rigid body motion, then Φ(m) = 0 if and only if m is a constant vector. Finally, observing that Φ behaves exactly like the functional Φ2 of Example 2.4 in [16], we see that Φ satisfies all assumptions of Theorem 4.2 or 5.3 in [16]. Thus 2 ! Z Ms Ms X X h h h 2 2 h 2 h (3.20) ∀v s ∈ Vs , |v s |H 1 (Ωs,i ) ≤ C kD(v s )kL2 (Ωs,i ) + v s × ns . i=1

Γsd

i=1

In order to recover the full norm of Xs in the left-hand side of (3.20), we apply Theorem 5.1 in [15] with the functional n−1 XZ |v s · τ l |. Φ(v) = l=1

Γsd

Again, Φ is a seminorm on Xs and since Γsd is a closed curve or surface, the condition Φ(c) = 0 for a constant vector c implies that c = 0. The remaining assumptions of this theorem easily follow by observing that Φ has the same behavior as the functional Φ1 of Example 4.2 in [15]. This yields ! Ms 2  n−1 XZ X . |v hs · τ l | (3.21) ∀v hs ∈ Vsh , kv hs k2L2 (Ωs ) ≤ C 2 |v hs |2H 1 (Ωs,i ) + i=1

l=1

Γsd

By combining (3.20) and (3.21), we obtain (3.22)

∀v hs

∈

Vsh ,

|||v hs |||2Xs

≤C

2

Ms X i=1

kD(v hs )k2L2 (Ωs,i )

Z +

Γsd

v hs

Then (3.18) follows from (3.22) by arguing as in deriving (2.31).

2  n−1 XZ × ns + l=1

Γsd

|v hs

· τ l|

2

!

. 

The case when Ωs is not connected follows from Lemmas 3.3 or 3.4 applied to each connected component of Ωs according that it is or is not adjacent to Γs . To control the bilinear form b in Ωs , we make the following assumption: There exists a linear approximation operator Θhs : H01 (Ω)n 7→ Vsh satisfying for all v ∈ H01 (Ω)n :

18

• ∀i, 1 ≤ i ≤ Ms ,

(3.23) • For any Γij in Γsd ,

Z

(3.24)

Z

Ωs,i


div Θhs (v) − v = 0.


Θhs (v) − v · nij = 0.

Γij

• There exists a constant C independent of v, h, H, and the diameter of Ωs,i , 1 ≤ i ≤ Ms , such that |||Θhs (v)|||Xs ≤ C|v|H 1 (Ω) .

(3.25)

A general construction strategy discussed in Section 4.1 gives an operator Θhs that satisfies (3.23) and (3.24). The stability bound (3.25) is shown to hold for the specific examples presented in Sections 7.1 and 7.2. Lemma 3.5. Assume that an operator Θhs satisfying (3.23)–(3.25) exists; then there exists a linear operator Πhs : H01 (Ω)n 7→ Vsh such that for all v ∈ H01 (Ω)n , Ms Z X wh div(Πhs (v) − v) = 0, (3.26) ∀wh ∈ Wsh , Ωs,i

i=1

∀Γij ∈ Γsd ,

(3.27)

Z


Πhs (v) − v · nij = 0,

Γij

and there exists a constant C independent of v, h, H, and the diameter of Ωs,i , 1 ≤ i ≤ Ms , such that |||Πhs (v)|||Xs ≤ C|v|H 1 (Ω) .

(3.28)

Proof. The operator Πhs is constructed by correcting Θhs : Πhs (v) = Θhs (v) + chs (v) h and where chs (v)|Ωs,i ∈ X0,s,i

(3.29)

h

∀w ∈

h W0,s,i ,1

≤ i ≤ Ms ,

Z

h

w div Ωs,i

chs (v)

=

Z

Ωs,i


wh div v − Θhs (v) .

Existence of chs (v) follows directly from (3.3) and with the same constant (3.30)

|chs (v)|H 1 (Ωs,i ) ≤

1 kdiv(v − Θhs (v))kL2 (Ωs,i ) . βs⋆

h is relaxed by applying (3.23) and using the fact that chs (v) belongs to The restriction wh ∈ W0,s,i h X0,s,i . Finally, (3.28) follows from the above bound and (3.25). 

The idea of constructing the operator Πhs via the interior inf-sup condition (3.3) and the simplified operator Θhs satisfying (3.23) and (3.25) is not new. It can be found for instance in [36] and [12]. To control the bilinear form b in Ωd , we make the following assumption: There exists a linear operator Πhd : H01 (Ω)n 7→ Vdh satisfying for all v ∈ H01 (Ω)n : • Md Z X
(3.31) ∀wh ∈ Wdh , wh div Πhd (v) − v = 0. i=1

Ωd,i

19

• For any Γij in Γsd , (3.32)

H

∀µ

∈

ΛH sd ,

Z

Γij


µH Πhd (v) − Πhs (v) · nij = 0.

• There exists a constant C independent of v, h, H, and the diameter of Ωd,i , 1 ≤ i ≤ Md , such that (3.33)

|||Πhd (v)|||Xd ≤ C|v|H 1 (Ω) .

The construction of the operator Πhd is presented in Section 4. In particular, the general construction strategy discussed in Section 4.1 gives an operator that satisfies (3.31) and (3.32). The stability bound (3.33) is shown to hold for various cases in Section 4.4. The next lemma follows readily from the properties of Πhs and Πhd . Lemma 3.6. Under the above assumptions, there exists a linear operator Πh ∈ L(H01 (Ω)n ; V h ) such that for all v ∈ H01 (Ω)n M Z X
h h (3.34) ∀w ∈ W , wh div Πh (v) − v = 0, i=1

(3.35)

Ωi

|||Πh (v)|||X ≤ C|v|H 1 (Ω) ,

with a constant C independent of v, h, H, and the diameter of Ωi , 1 ≤ i ≤ M .

Proof. Take Πh (v)|Ωs = Πhs (v) and Πh (v)|Ωd = Πhd (v). Then (3.34) follows from (3.26) and (3.31). The matching condition of the functions of V h at the interfaces of Γsd holds by virtue of (3.32). Finally, the stability bound (3.35) stems from (3.28) and (3.33).  The following inf-sup condition between W h and V h is an immediate consequence of a simple variant of Fortin’s Lemma [34, 18] and Lemma 3.6. Theorem 3.1. Under the above assumptions, there exists a constant β ⋆ > 0, independent of h, H, and the diameter of Γij for all i < j, such that (3.36)

∀wh ∈ W h , sup

vh ∈V h

b(v h , wh ) ≥ β ⋆ kwh kL2 (Ω) . |||v h |||X

Finally, well-posedness of the discrete scheme (3.13) follows from Lemma 3.3 or 3.4 and Theorem 3.1. Corollary 3.1. Under the above assumptions, problem (3.13) has a unique solution (uh , ph ) ∈ V h × W h and
(3.37) |||uh |||X + kph kL2 (Ω) ≤ C kf kL2 (Ω) + kqd kL2 (Ωd ) , with a constant C independent of h, H, and the diameter of Γij for all i < j.

Proof. Since (3.13) is set into finite dimension, it is sufficient to prove uniqueness, and as uniqueness follows from (3.37), it suffices to prove this stability estimate. Thus let (uh , ph ) solve (3.13), which is a typical linear problem with a non-homogeneous constraint. By virtue of the discrete inf-sup condition (3.36), there exists a function uhq ∈ V h such that Z h h h h ∀w ∈ W , b(uq , w ) = − w h qd , Ωd

(3.38)

|||uhq |||X ≤

1 kqd kL2 (Ωd ) . β⋆

20

Then uh0 = uh − uhq solves (3.13) with qd = 0 and the coercivity condition (3.15) or (3.18) and the discrete inf-sup condition (3.36) imply that |||uh0 |||X + kph kL2 (Ω) ≤ C(kf kL2 (Ω) + kqd kL2 (Ωd ) ).

With (3.38), this gives (3.37).



4. Elements of construction of Θhs and Πhd Constructing an operator Θhs with values in Vsh , satisfying (3.23)–(3.25), uniformly stable with respect to the diameter of the subdomains and interfaces, is not straightforward, particularly in 3-D. On the other hand, a general construction of Πhd in Ωd can be found in [5], and we shall adapt it so that it matches suitably Θhs on Γsd . Let us describe our strategy in each region. 4.1. General construction strategy. Let v ∈ H01 (Ω)n . In Ωs , we propose the following threestep construction. (1) Starting step. In each Ωs,i , 1 ≤ i ≤ Ms , take Θhs (v) = S h (v), where S h (v) is a Scott & Zhang [56] approximation operator, constructed so that S h (v)|∂Ωs,i only uses values of v h and S h (v) ∈ X h . restricted to ∂Ωs,i , and in particular, S h (v)|Γs = 0. Thus S h (v)|Ωs,i ∈ Xs,i s h (2) First correction step. For each Γij ∈ Γss ∪ Γsd with i < j, correct Θs (v) in Ωs,i by setting Θhs (v)|Ωs,i := Θhs (v)|Ωs,i + chi,Γij (v),

(4.1)

h , ch where chi,Γij (v) ∈ Xs,i i,Γij (v) = 0 on ∂Ωs,i \ Γij , Z Z
h v − S h (v)|Ωs,i · nij , ci,Γij (v) · nij = Γij

Γij

(4.2)

and satisfies suitable bounds. The existence of chi,Γij (v) (without bounds) is guaranteed by assumption (3.7). In particular, one can define it on Γij to satisfy (4.1), extend it by zero h . Note that the on ∂Ωs,i \ Γij , and extend it arbitrarily inside Ωs,i to a function in Xs,i h correction ci,Γij (v) only affects values at points located in the interior of Ωs,i and in the interior of Γij . It has no influence on values at points located on the other interfaces. For this reason, the discrete term that appears in the right-hand side of (4.1) is the trace on Γij of S h (v) coming from Ωs,i ; this term has not been yet corrected. Once this correction step is performed on all Γij ∈ Γss ∪ Γsd with i < j, the resulting function Θhs (v) satisfies on all these interfaces Z
Θhs (v)|Ωs,i − v · nij = 0. Γij

Note that this step does not modify the trace on Γij of S h (v) coming from Ωs,j . This trace will be modified in the next step. (3) Second correction step. For each Γij ∈ Γss with i < j, correct Θhs (v) in Ωs,j by setting Θhs (v)|Ωs,j := Θhs (v)|Ωs,j + chj,Γij (v),

(4.3)

h , ch where chj,Γij (v) ∈ Xs,j j,Γij (v) = 0 on ∂Ωs,j \ Γij , Z Z
H H H h ∀µ ∈ Λs , µ · cj,Γij (v) = µH · Θhs (v)|Ωs,i − S h (v)|Ωs,j , Γij

Γij

and satisfies suitable bounds. The existence of chj,Γij (v) (without bounds) is guaranteed by assumption (3.8). In particular, one can define it on Γij to satisfy (4.3), extend it by h . Here also, zero on ∂Ωs,j \ Γij , and extend it arbitrarily inside Ωs,j to a function in Xs,j the correction chj,Γij (v) has no influence on values at points located on the other interfaces.

21

(4.4)

Once this correction step is performed on all Γij ∈ Γss with i < j, the resulting function Θhs (v) satisfies on all these interfaces bs (Θhs (v), µH ) = 0 for all µH ∈ ΛH s , i.e. Z ∀µH ∈ ΛH [Θhs (v)] · µH = 0. s , Γij

(4.5)

Finally, if the second assumption in Hypothesis 3.2 holds, then n ∈ ΛH s and (4.2) and (4.4) imply that on all these interfaces, Z
Θhs (v)|Ωs,j − v · n = 0. Γij

Θhs (v)

Therefore satisfies (3.23) and (3.24). Specific constructions of the corrections chi,Γij (v) and chj,Γij (v) that guarantee that Θhs (v) also satisfies (3.25) are presented in the forthcoming subsections. Next, we propose the following two-step construction algorithm in Ωd . (1) Starting step. Set Pdh (v) = Rh (v) ∈ Xdh , where Rh (v) is a standard mixed approximation operator associated with Wdh . It preserves the normal component on the boundary: Z
(4.6) ∀Γij ⊂ ∂Ωd,k , 1 ≤ k ≤ Md , ∀v h ∈ Xdh , v h · nij Rh (v)|Ωd,k − v · nij = 0, Γij

and satisfies

(4.7)

∀1 ≤ i ≤ Md , ∀wh ∈ Wdh ,

Z

Ωd,i


wh div Rh (v) − v = 0.

(2) Correction step. It remains to prescribe the jump condition. For each Γij ∈ Γdd ∪ Γsd with i < j, correct Pdh (v) in Ωd,j by setting: Pdh (v)|Ωd,j := Pdh (v)|Ωd,j + chj,Γij (v),

(4.8)

h , ch h where chj,Γij (v) ∈ Xd,j j,Γij (v) · nj = 0 on ∂Ωd,j \ Γij , div cj,Γij (v) = 0 in Ωd,j , Z Z
H H H h ∀µ ∈ Λd , µ cj,Γij (v) · nij = µH Rh (v)|Ωd,i − Rh (v)|Ωd,j · nij , Γij

∀µH ∈ ΛH sd ,

Z

Γij

Γij

µH chj,Γij (v) · nij =

Z

Γij


µH Πhs (v)|Ωs,i − Rh (v)|Ωd,j · nij ,

and chj,Γij (v) satisfies adequate bounds. Existence of a non necessarily divergence-free chj,Γij (v) (without bounds) follows from (3.9); it suffices to extend suitably Rh (v)|Ωd,j and Rh (v)|Ωd,i or Πhs (v)|Ωs,i . The zero divergence will be prescribed in the examples. Note that chj,Γij (v) has no effect on interfaces other than Γij and no effect on the restriction of Pdh (v) in Ωd,i or on that of Πhs (v) in Ωs,i . Therefore these corrections can be superimposed. When step 2 is done on all Γij ∈ Γdd ∪ Γsd with i < j, the resulting function Pdh (v) has zero normal trace on Γd , it belongs to Vdh since, due to the first equation in (4.8), it satisfies for all Γij ∈ Γdd with i < j Z (4.9) ∀µH ∈ ΛH , µH [Pdh (v) · n] = 0, d Γij

and, as the corrections are assumed to be divergence-free in each subdomain, Z
h h (4.10) ∀w ∈ Wd , ∀1 ≤ i ≤ Md , wh div Pdh (v) − v = 0. Ωd,i

22

Furthermore, due to the second equation in (4.8), it satisfies for all Γij ∈ Γsd , Z
µH Πhs (v)|Ωs,i − Pdh (v)|Ωd,j · nij = 0. (4.11) ∀µH ∈ ΛH , sd Γij

Πhd (v)

Pdh (v)

Therefore, taking = in Ωd , it satisfies (3.31) and (3.32). The remainder of this section is devoted to corrections chi,Γij (v) and chj,Γij (v) where the constants in the stability bounds (3.25) and (3.33) are shown to be independent of the discretization parameters and the diameter of the subdomains. In particular, the bounds for the corrections will stem from the fact that each involves differences between v and good approximations of v. Beforehand, we need to refine the assumptions on the meshes at the interfaces and refine Hypothesis 3.1 on the mesh of subdomains. Hypothesis 4.1. For i < j, let T be any element of Tih that is adjacent to Γij , and let {Tℓ } denote the set of elements of Tjh that intersect T . The number of elements in this set is bounded by a fixed integer and there exists a constant C such that |Tℓ | |T |−1 ≤ C.

The same is true if the indices i and j of the triangulations are interchanged. These constants are independent of i, j, h, and the diameters of the interfaces and subdomains. Hypothesis 4.2. (1) Each Ωi , 1 ≤ i ≤ M , is the image of a “reference” polygonal or polyhedral domain by an homothety and a rigid body motion: ˆ i ) , x = Fi (ˆ ˆ + bi , (4.12) Ωi = Fi (Ω x ) = A i Ri x where Ai = diam(Ωi ), Ri is an orthogonal matrix with constant coefficients and bi a constant vector. (2) There exists a constant σ1 independent of M such that for any pair of adjacent subdomains Ωi and Ωj , 1 ≤ i, j ≤ M , we have (4.13)

Ai A−1 j ≤ σ1 .

ˆ i ) = 1. In addition, it follows from Hypothesis 3.1 that on one hand the By (4.12) diam(Ω ˆ i can take at most L geometrical shapes and on the other hand, reference domains Ω ˆi ) ≥ 1 , (4.14) ∀i, 1 ≤ i ≤ M , diam(B σ ˆ ˆ where Bi is the largest ball contained in Ωi and σ is the constant of (3.1). 4.2. A construction of chi,Γij (v). In this section we construct a correction chi,Γij (v) satisfying (4.1) and suitable bounds needed to establish the stability estimate (3.25). Recall the approximation properties of the Scott & Zhang operator of degree r ≥ 1. Let v ∈ H t (Ωs,i ), for some real number t ≥ 1 and let T be an element of Tih , 1 ≤ i ≤ Ms . Let ∆T denote the macro-element that is used for defining the values of S h (v) in T . Then (cf. [56]) there exists a constant C that depends only on r, t, and the shape regularity of Tih , such that the following local approximation error holds (4.15)

min(r,t)

kv − S h (v)kL2 (T ) + hT |v − S h (v)|H 1 (T ) ≤ ChT

|v|H min(r,t) (∆T ) .

Owing to the regularity of Tih , when (4.15) is summed over all T in Tih , it gives (4.16)

kv − S h (v)kL2 (Ωs,i ) + h|v − S h (v)|H 1 (Ωs,i ) ≤ Chmin(r,t) |v|H min(r,t) (Ωs,i ) .

Let v ∈ H01 (Ω)n . Consider the case when n = 3, the case n = 2 being simpler, and also consider h the case when Γij is a polyhedral surface, not necessarily a flat plane. Let Xi,Γ denote the trace ij h h h h of Xs,i on Γij , and Ti,Γij the trace of Ti on Γij , 1 ≤ i ≤ Ms , i < j; Ti,Γij is a triangular mesh of

23 h contains each flat face in Γij . In all cases considered, the restriction to each T of functions of Xs,i n h , at least IP 1 . Now, choose one of the faces, say F , of Γij . For each interior vertex ak of Ti,Γ ij h 1 ≤ k ≤ NF , on F , let Ok denote the macro-element of all triangles of Ti,Γ (i.e. faces on F ) that ij share the vertex ak . Thus F F = ∪N k=1 Ok . The set {Ok } is not a partition of F , but it can be transformed into a partition by setting k−1 ∆1 = O1 , and recursively ∆k = Ok \ ∪ℓ=1 ∆ℓ .

Note that some ∆k may be empty. By construction, we have

F F = ∪N k=1 ∆k , , ∆k ∩ ∆ℓ = ∅ , k 6= ℓ , ∆k ⊂ Ok .

This can be done for all flat faces of Γij . For each k, let bk be the piecewise IP 1 “bubble” function such that bk (aℓ ) = δk,ℓ , extended by zero to all vertices inside Ωs,i , and define chi,Γij (v) by (4.17)

chi,Γij (v) =

NF X X

F ⊂Γij k=1

ck bk , where ck = R

1

Z

Ok bk

∆k


v − S h (v)|Ωs,i , ck = 0 , if ∆k = ∅.

Lemma 4.1. The correction chi,Γij (v) defined by (4.17) satisfies (4.1), more precisely Z Z
v − S h (v)|Ωs,i , chi,Γij (v) = (4.18) ∀F ⊂ Γij , F

F

and there exists a constant C independent of h and the diameter of Ωs,i and Γij such that

(4.19)

∀v ∈ H01 (Ω)n , |chi,Γij (v)|H 1 (Ωs,i ) ≤ C|v|H 1 (Ωs,i ) ,

kchi,Γij (v)kL2 (Ωs,i ) ≤ Ch|v|H 1 (Ωs,i ) .

Proof. For each k, the support of the trace of bk on Γij is contained in a single flat plane, say F , therefore, Z Z Z X Z NF Z NF NF NF X X X
h v − S h (v)|Ωs,i ck ck bk = bk = c k bk = ci,Γij (v) = F k=1

F

=

Z

F

k=1

F

k=1

Ok

k=1

∆k


v − S h (v)|Ωs,i ,

since the set (∆k ) is a partition of Γ. To derive the estimate (4.19), we consider the faces T ′ in ∆k , pass to the reference element Tˆ, where T is the element of Tih adjacent to T ′ , apply a trace theorem in Tˆ, and revert to T . This gives Z X
h kv − S h (v)|Ωs,i kL1 (T ′ ) v − S (v)|Ωs,i ≤ ∆k

T ′ ⊂∆k

≤Cˆ

X

T ′ ⊂∆k

  1 |T ′ | |T |− 2 kv − S h (v)kL2 (T ) + hT |v − S h (v)|H 1 (T ) .

In this proof, Cˆ denotes constants that depend only on the reference element and the shape regularity of the triangulation. Then considering the regularity of Tih , the local approximation formula (4.15) with r = t = 1 implies that Z X
− 21 ′ − 12 h ˆ ≤ Cˆ |v| ≤ C|∆ | ρ h |T | |T | v − S (v)| (4.20) 1 T Ω k H (∆T ) s,i k |v|H 1 (Dk ) , ∆k

T ′ ⊂∆k

24

where Dk is the set of elements of Tih where the values of v are taken for computing S h (v), and ρk = minT ∈Dk ρT . Therefore, considering that |∆k | ≤ |Ok |, ck defined by (4.17) satisfies −1

|ck | ≤ Cˆ ρk 2 |v|H 1 (Dk ) .

(4.21) Now, |chi,Γij (v)|2H 1 (Ωs,i ) =

Z

Ωs,i

NF NF X Z X X 2 2 X X c k ∇ bk . c k ∇ bk = F ⊂Γij k=1

T

T ∈Tih

F ⊂Γij k=1

But given an element T in Tih there is at most a fixed (and small) number of indices k such that bk |T 6= 0. Therefore |chi,Γij (v)|2H 1 (T ) =

Z

T

NF X 2 X X |ck |2 |bk |2H 1 (T ) , ck ∇ bk ≤ Cˆ F ⊂Γij k=1

k

where the sum runs over all indices k such that bk |T 6= 0. By substituting (4.21) into this inequality and estimating the norm of bk , we easily derive that X 1 1 (4.22) |chi,Γij (v)|2H 1 (T ) ≤ Cˆ |T ||v|2H 1 (Dk ) . ρk ρ2T k

Since ρk ρ2T has the same order as |T |, then by summing (4.22) over all T in Tih and considering that there is at most a fixed (and small) number of repetitions in the Dk , we obtain the first inequality in (4.19). The second inequality in (4.19) follows in a similar manner from the representation (4.17) and bound (4.21).  4.3. A construction of chj,Γij (v) in Ωs . The 2-D case. Constructing a suitable correction solving (4.3) is fairly complex because it amounts to satisfying many conditions per interface. In 2-D, it can be done by following the sharp approach of Crouzeix & Thom´ee [23]. This is a general procedure that can also be used in 3-D, but as observed in Remark 4.1, it restricts the meshes. To bypass this restriction, we present a more ad hoc construction in 3-D, which is postponed to the Appendix: Section 7.2. We slightly modify the underlying Scott & Zhang operator by asking that the modified operator S˜h be continuous at the end points of all interfaces Γij in Γss ∪ Γsd , and coincide elsewhere with S h . This only concerns the set of points, say ak , 1 ≤ k ≤ P , that do not lie on Γs , since S h (v) is necessarily zero on Γs . Note that these are all interior cross points of the triangulation of subdomains TΩ . To achieve continuity at any such point a, we pick an element Ta in Ωs with vertex a and set Z 1 h ˜ v. S (v)(a) = |Ta| Ta Both corrections chi,Γij (v) and chj,Γij (v) are defined to satisfy respectively (4.1) and (4.3) with S h (v) replaced by S˜h (v). For v ∈ H 1 (Ω), this does not change the approximation properties (4.16) of S˜h and hence (4.19), except that the domain of definition of v is possibly a little larger than Ωs,i or Ωs,j near the end points of Γij .

˜ s,k the possibly slightly enlarged domain of definition Notation 4.1. From now on, we denote by Ω h ˜ of v in case S (v) is used, with the convention that it coincides with Ωs,k when S˜h (v) is not used.

25

Owing to the continuity of S˜h (v) at the end points of Γij and the fact that all previously made corrections vanish at the end points of the interfaces Γij in Γss ∪Γsd , the integrand in the right-hand side of (4.3) also vanishes at the end points of Γij . Thus the trace of Θhs (v)|Ωs,i − S˜h (v)|Ωs,j on 1

2 (Γij )2 , see Section 1.1. We shall see that this property is crucial for estimating Γij belongs to H00 uniformly the gradient of chj,Γij (v). We propose to split the construction of chj,Γij (v) into two parts: 1

1

h 2 2 ) such that for all ℓ ∈ H00 (Γij )2 , Xj,Γ (Γij )2 , (1) Construct an operator πsh ∈ L(H00 ij Z H ∀µH ∈ Λs , (πsh (ℓ) − ℓ) · µH = 0, Γij

(4.23)

πsh (ℓ) = 0, at the end points of Γij ,

|πsh (ℓ)|

1

2 (Γ ) H00 ij

≤ κs |ℓ|

1

2 (Γ ) H00 ij

,

with a constant κs that is independent of h, H, and the measure of Γij . (2) Define chj,Γij (v) to be a suitable extension of πsh (Θhs (v)|Ωs,i − S˜h (v)|Ωs,j ) inside Ωs,j so that it satisfies all uniform properties required for the stability bound (3.25). Note that chj,Γij (v) satisfies (4.3) by construction. Also, the existence of the operator πsh implies that condition (3.8) holds. We now discuss the construction of an operator satisfying (4.23). Uniformity with respect to ˆ j be the diameters of Γij and Ωs,j is obtained by reverting to the reference subdomains. Let Ω −1 −1 h ˆ = F (Γij ). Let the image by F associated with Ωs,j by Hypothesis 4.2 and let Γ of Tj,Γ be j j ij ˆ denoted by Tj . Similarly, set h ˆ j = {ˆ ˆ = {µ ˆ = µH ◦Fj−1 ; µH ∈ ΛH X v = v h ◦Fj−1 ; v h ∈ Xj,Γ , v h = 0, at the end points of Γij } , Λ s }. ij 1

2 ˆ 2 ˆ j unique solution of (Γ) , if we construct π ˆj (ˆℓ) in X Then, given ˆℓ in H00 Z Z ˆ ˆℓ · µ, ˆ ˆ= ˆ ˆ ∈ Λ, (4.24) ∀µ π ˆj (ℓ) · µ

ˆ Γ

ˆ Γ

and satisfying (4.25)

|ˆ πj (ˆℓ)|

1

2 (Γ) ˆ H00

ˆ ˆℓ| < C|

1

2 (Γ) ˆ H00

,

with a constant Cˆ independent of j and ˆℓ, then by reverting to Ωs,j and defining πjh (ℓ) by (4.26)

πjh (ℓ) ◦ Fj = π ˆj (ℓ ◦ Fj ) = π ˆj (ˆℓ),

we shall obtain an operator satisfying (4.23). Indeed, uniqueness of the solution of (4.24) will guarantee that the mapping π ˆj is linear, and the inequality in (4.23) with the same constant will 1

2 seminorm by a homothety and a rigid body motion. But follow from the invariance of the H00 since explicit constructions of π ˆj are fairly technical, we postpone them to the Appendix: Section 7.1 and discuss the extension part (2) now. First we construct a suitable function v in H 1 (Ωs,j ). 1 2 Lemma 4.2. For any Ωs,j ⊂ Ωs , any interface Γ ⊂ ∂Ωs,j and any ℓ ∈ H00 (Γ) there exists a 1 function v = E(ℓ) ∈ H (Ωs,j ) satisfying

v|Γ = ℓ , v|∂Ωs,j \Γ = 0,

26

the mapping E is linear and there exists a constant C independent of ℓ, v, and the diameter of Γ and Ωs,j such that |v|H 1 (Ωs,j ) ≤ C|ℓ|

(4.27)

1

2 (Γ) H00

, kvkL2 (Ωs,j ) ≤ CAj |ℓ|

1

2 (Γ) H00

,

where Aj is the diameter of Ωs,j . 1

2 ˆ Proof. With the notation of Hypothesis 4.2, let ℓˆ = ℓ ◦ Fj . Then ℓˆ belongs to H00 (Γ); hence it can 1 ˆ ˆ ˆ 2 be extended by zero to ∂ Ωj and the extended function, say ℓ0 belongs to H (∂ Ωj ) with

|ℓˆ0 |

(4.28)

1

ˆj) H 2 (∂ Ω

ˆ ˆ ℓ| ≤ C|

1

2 (Γ) ˆ H00

.

ˆ ∈ L(H 12 (∂ Ω ˆ j ); H 1 (Ω ˆ j )) such that There exists an extension operator E ˆ kEk

(4.29)

1

ˆ j );H 1 (Ω ˆ j )) L(H 2 (∂ Ω

ˆ ≤ C.

ˆ ℓˆ0 ) and v = vˆ ◦ F −1 . The mapping ℓ 7→ v is linear, v belongs to H 1 (Ωs,j ), its trace Take vˆ = E( j satisfies the statement of the lemma, and it remains to check (4.27). The following inequalities stem from (4.29), (4.28) and a straightforward scaling argument: ˆ ℓˆ0 | |v|H 1 (Ωs,j ) ≤ |ˆ v |H 1 (Ωˆ j ) ≤ C|

(4.30)

1

ˆj) H 2 (∂ Ω

ˆ j kˆ ˆ j |ℓ| kvkL2 (Ωs,j ) ≤ CA v kL2 (Ωˆ j ) ≤ CA

ˆ ˆ ℓ| ≤ C| 1

2 (Γ) H00

1

2 (Γ) ˆ H00

ˆ ≤ C|ℓ|

1

2 (Γ) H00

,

. 

Next, we take Γ = Γij and define chj,Γij (v) = S h (w),

(4.31)

w = E(πsh (Θhs (v)|Ωs,i − S˜h (v)|Ωs,j )),

where the Scott & Zhang interpolant S h is constructed so that S h (w) vanishes on ∂Ωs,j \ Γij . The uniform approximation properties (4.16) of S h , (4.25), and (4.27) imply, with constants C independent of h and the diameters of Γij , Ωs,i and Ωs,j : , (4.32) kchj,Γij (v)kH 1 (Ωs,j ) ≤ C|w|H 1 (Ωs,j ) ≤ C {Θhs (v)|Ωs,i − S˜h (v)|Ωs,j } 12 H00 (Γij )

and it remains to derive a uniform bound for the last norm in terms of v. This is the object of the next lemma.

Lemma 4.3. Let Γij ∈ Γss , i < j. There exists a constant C independent of h and the diameters of Γij , Ωs,i and Ωs,j , such that ≤ C|v|H 1 (Ω˜ s,i ∪Ω˜ s,j ) , (4.33) ∀v ∈ H01 (Ω)2 , {Θhs (v)|Ωs,i − S˜h (v)|Ωs,j } 21 H00 (Γij )

˜ s,k is defined in Notation 4.1. where Ω

Proof. Recall that the trace on Γij of Θhs (v)|Ωs,i is S˜h (v)|Ωs,i + chi,Γij (v). As chi,Γij (v) vanishes on ˆ i , we obtain ∂Ωs,i \ Γij , then by passing to the reference subdomain Ω |chi,Γij (v)|

1

2 (Γ ) H00 ij

= |chi,Γij (v) ◦ Fi |

1

2 (Γ) ˆ H00

ˆ hi,Γ (v) ◦ Fi | 1 ˆ ≤ C|c ˆ hi,Γ (v)|H 1 (Ω ) . ≤ C|c s,i H (Ω i ) ij ij

˜ s,i instead of Ωs,i . Thus it suffices to consider the Therefore it is bounded owing to (4.19), with Ω h ˜ jump [S (v)] through Γij .

27

ˆ k , k = i, j, and First, by writing [S˜h (v)] = [S˜h (v) − v], by passing to the reference subdomains Ω h by using the approximation properties (4.16) of S , we derive h
ˆ S˜h (v) − v) ◦ Fk k 1 ˆ {(S˜h (v) − v) ◦ Fk | ˆ } 1 S˜ (v) − v |Ω 1 = ≤ Ck( s,k H (Ω k ) Ωk 2 2 ˆ H (Γij )

(4.34)

H (Γ)



2 ˜h ≤ Cˆ |S˜h (v) − v|2H 1 (Ωs,k ) + A−2 k kS (v) − vkL2 (Ωs,k )  1 2 2 2 . ≤ Cˆ |v|2H 1 (Ω˜ ) + (h A−1 ) |v| ˜ ) k H 1 (Ω s,k

1

2

s,k

Since h < Ak , this bounds the first part of the norm and it remains to estimate Z 2 1 ˜h [S (v)] , Γij d(x)

where d(x) denotes the distance between x and the end points of Γij . This part does not have the same immediate bound as the above correction because the jump does not vanish on the other boundaries of the subdomains. Let us first consider the case when Γij is a straight line. There is no loss of generality in assuming that Γij lies on the axis y = 0 with end point at the origin: h Γij =]0, L[. Let 0 = x0 < x1 < · · · < xN < xN +1 = L denote the nodes of Ti,Γ and exceptionally ij set hn = xn+1 − xn . Then d(x) = Min(x, L − x) and since Z L Z L Z 2 2 2 1 ˜h 1 ˜h 1 ˜h [S (v)] ≤ [S (v)] + [S (v)] , 0 x 0 L−x Γij d(x) it suffices to examine the first term: Z L N Z xn+1 2 X 2 1 ˜h 1 ˜h [S (v)] = [S (v)] . x xn 0 x n=0


2 For any n ≥ 1, we have xn ≥ hn−1 , and since [S˜h (v)] belongs to H 1 (0, L) ∩ C 0 (0, L) , we can write Z xn+1 Z xn+1 Z x 2
2 1 ˜h 1 d ˜h h ˜ [S (v)] ≤ [S (v)(t)]dt dx [S (v)(xn )] + x xn xn xn xn dt (4.35)
hn d ˜h 2hn ˜h [S (v)(xn )]2 + k [S (v)]k2L2 (xn ,xn+1 ) . ≤ hn−1 2 dx An equivalence of norms yields 2hn ˜h 1 2hn ˜h [S (v)(xn )]|2 ≤ k[S (v)]k2L∞ (xn ,xn+1 ) ≤ Cˆ k[S˜h (v)]k2L2 (xn ,xn+1 ) hn−1 hn−1 hn−1 1 = Cˆ k(S˜h (v)|Ωs,i − v) − (S˜h (v)|Ωs,j − v)k2L2 (xn ,xn+1 ) . hn−1 By arguing as in the proof of Lemma 4.1, we obtain 1

(4.36)

ˆ 2 |v| 1 ˜ , kS˜h (v)|Ωs,i − vkL2 (xn ,xn+1 ) ≤ Ch T H (∆n,i )

˜ n,k denotes the set of elements in Ωs,k where T is the element of Tih adjacent to [xn , xn+1 ] and ∆ used in defining S˜h (v)|Ωs,k on [xn , xn+1 ]. The estimation of the second term is more technical because now S˜h (v) is constructed on Tjh whereas [xn , xn+1 ] is the side of an element of Tih . Let {Tℓ } denote the set of elements of Tjh that intersect [xn , xn+1 ]. The argument of Lemma 4.1 gives: X 1 hT2ℓ |v|H 1 (∆ (4.37) kS˜h (v)|Ωs,j − vkL2 (xn ,xn+1 ) ≤ Cˆ ˜ n,ℓ ) . ℓ

28

Then combining (4.36) and (4.37), using the regularity of the triangulation and Hypothesis 4.1, we obtain
2hn ˜h 2 (4.38) [S (v)(xn )]|2 ≤ Cˆ |v|2H 1 (∆ ˜ n,i ) + |v|H 1 (∆ ˜ n,j ) . hn−1

Similarly,

(4.39)

2 1 h2n 1

d ˜h

≤ Cˆ [ S (v)] k[S˜h (v)]k2L2 (xn ,xn+1 ) = Cˆ k[S˜h (v) − v]k2L2 (xn ,xn+1 )

2 hn−1 dx hn−1 hn−1 L (xn ,xn+1 )
≤ Cˆ |v|2 1 ˜ + |v|2 1 ˜ . H (∆n,i )

H (∆n,j )

When n = 0, as [S˜h (v)(0)] = 0, there only remains the second term in the first line of (4.35), and we immediately derive Z

x1 0

Z x1  Z x 2
2 d ˜h 1 1 ˜h [S (v)] dt dx [S (v)] = x x 0 dt 0 Z x1

d




d ˜h

2

2 2 h ˜ ≤ ≤ h1 ≤ Cˆ |v|2H 1 (∆ [S (v)] 2 [S (v)] 2

˜ 0,i ) + |v|H 1 (∆ ˜ 0,j ) . dx dx L (0,x) L (0,x1 ) 0

Finally, when Γij is a polygonal line, the above argument is valid on the segments that share an end point with Γij and is simpler on the other segments as d(x) is not small compared to h.  Inequalities (4.32) and (4.33) imply the following. Corollary 4.1. Let Γij ∈ Γss , i < j. There exists a constant C independent of h and the diameters of Γij , Ωs,i and Ωs,j , such that (4.40)

∀v ∈ H01 (Ω)2 , kchj,Γij (v)kH 1 (Ωs,j ) ≤ C|v|H 1 (Ω˜ s,i ∪Ω˜ s,j ) .

Corollary 4.2. The approximation operator Θhs constructed in Section 4.1 with the above corrections chi,Γij (v) and chj,Γij (v) satisfies assumption (3.25). Proof. Since Θhs (v) is constructed by correcting the Scott-Zhang interpolant S˜h with chi,Γij (v) and chj,Γij (v), assumption (3.25) follows from (4.16), (4.19), (4.40), and the Poincar´e inequality (1.1).  Remark 4.1. The above correction chj,Γij (v) has a straightforward extension in 3-D, but its use is limited because now it requires continuity of S˜h (v) on the edges of subdomains; thus the meshes must match on these edges. This is not required by the correction defined in Section 7.2. 4.4. A construction of chj,Γij (v) in Ωd . Here we construct a correction chj,Γij (v) in Ωd satisfying (4.8) and suitable continuity bounds that are needed to establish the stability estimate (3.33). Recall that the existence of chj,Γij (v) relies on (3.9). In the construction below we directly show that (3.9) holds for a wide range of mesh configurations. Let v be given in H01 (Ω)n . Recall that the mixed approximation operator Rh defined in each Ωd,i takes its values in Xdh and satisfies (4.6) on each Γij ⊂ ∂Ωd,k , 1 ≤ k ≤ Md , and (4.7) in each Ωd,i , 1 ≤ i ≤ Md : Z
v h · nij Rh (v)|Ωd,k − v · nij = 0, ∀v h ∈ Xdh , Γij

h

∀w ∈

Wdh ,

Z

Ωd,i


wh div Rh (v) − v = 0.

29

Furthermore there exists a constant C independent of h and the geometry of Ωd,i , such that (4.41)

∀v ∈ H01 (Ω)n , kRh (v)kH(div;Ωd,i ) ≤ CkvkH 1 (Ωd,i ) , 1 ≤ i ≤ Md .

This is easily established by observing that the moments defining the degrees of freedom of Rh (v) are invariant by homothety and rigid-body motion; in particular the normal vector is preserved. In addition, it satisfies (3.6): h h . , div v h ∈ Wd,i ∀i, 1 ≤ i ≤ Md , ∀v h ∈ Xd,i

The above properties also imply (3.5): for all i, 1 ≤ i ≤ Md , R h h Ωd,i w div v ≥ βd⋆ , sup inf h kv h kH(div;Ωd,i ) kwh kL2 (Ωd,i ) wh ∈W0,d,i vh ∈X h 0,d,i

βd⋆

with a constant > 0 independent of h and Ai . h Now, let Γij ∈ Γdd ∪ Γsd ; by analogy with the situation in the Stokes region, we denote by Xd,j,Γ ij h on Γ . Following [5], we define the space of normal traces the trace space of Xd,j ij n h Xj,Γ = {w · nij ; w ∈ Xd,j,Γ }, ij ij n . Then we make the following assumpand the orthogonal projection Qhj,Γij from L2 (Γij ) into Xj,Γ ij tion: There exists a constant C, independent of H, h, i, j, and the diameters of Γij and Ωd,j , such that

(4.42)

H H H h H ∀µH ∈ ΛH d , ∀µ ∈ Λsd , kµ kL2 (Γij ) ≤ CkQj,Γij (µ )kL2 (Γij ) .

It is shown in [60] that (4.42) holds for both continuous and discontinuous mortar spaces, if the h . A similar inequality for more general grid configumortar grid TijH is a coarsening by two of Tj,Γ ij rations is shown in [52]. Formula (4.42) implies that the projection Qhj,Γij is an isomorphism from H H H n the restriction of ΛH sd , respectively Λd , to Γij , say Λsd,ij respectively Λd,ij , onto its image in Xj,Γij , and the norm of its inverse is bounded by C. Then a standard algebraic argument shows that its n H dual operator, namely the orthogonal projection from Xj,Γ into ΛH sd,ij , respectively Λd,ij , is also ij n an isomorphism from the orthogonal complement in Xj,Γ of the projection’s kernel onto ΛH sd,ij , ij H respectively Λd,ij , and the norm of its inverse is also bounded by C. As a consequence, for each n such that f ∈ L2 (Γij ), there exists v h · nij ∈ Xj,Γ ij Z Z H H H h ∀µH ∈ Λd , ∀µH ∈ Λsd , µ v · nij = f µH , Γij

Γij

and there exists a constant C independent of h, and the diameter of Γij , such that kv h · nij kL2 (Γij ) ≤ Ckf kL2 (Γij ) . This implies that (3.9) holds. Furthermore, the solution v h · nij is unique in the orthogonal complement of the projection’s kernel and by virtue of this uniqueness, v h · nij depends linearly on f . This result permits to partially solve (4.8). Lemma 4.4. Let v ∈ H01 (Ω)n . Under assumption (4.42), for each Γij ∈ Γdd ∪ Γsd , there exists n wh · nij ∈ Xj,Γ such that ij Z Z   H H H h ∀µ ∈ Λd , µH Rh (v) · n , µ w · nij = Γij Γij (4.43)   h kw · nij kL2 (Γij ) ≤ Ck Rh (v) · n kL2 (Γij ) ,

30

(4.44)

∀µH ∈ ΛH sd ,

Z

Γij

µH wh · nij = h

Z

Γij


µH Θhs (v)|Ωs,i − Rh (v)|Ωd,j · nij ,


kw · nij kL2 (Γij ) ≤ Ck Θhs (v)|Ωs,i − Rh (v)|Ωd,j · nij kL2 (Γij ) ,

with the constant C of (4.42). The mapping v 7→ wh · nij is linear.

Lemma 4.4 constructs a normal trace wh · nij on Γij and we must extend it inside Ωd,j . To this end, let ℓh ∈ L2 (∂Ωd,j ) be the extension of wh · nij by zero on ∂Ωd,j . Next, we solve the problem: Find q ∈ H 1 (Ωd,j ) ∩ L20 (Ωd,j ) such that (4.45)

∆ q = 0 in Ωd,j ,

∂q = ℓh on ∂Ωd,j . ∂nj

3

Lemma 4.5. Problem (4.45) has one and only one solution q ∈ H 2 (Ωd,j ) ∩ L20 (Ωd,j ) and p |q|H 1 (Ωd,j ) ≤ C Aj k[Rh (v) · n]kL2 (Γij ) ,

(4.46)

|q|

3

H 2 (Ωd,j )

|q|H 1 (Ωd,j ) (4.47)

|q|

3

H 2 (Ωd,j )

≤ Ck[Rh (v) · n]kL2 (Γij ) , Γij ∈ Γdd , p
≤ C Aj k Θhs (v)|Ωs,i − Rh (v)|Ωd,j · nij kL2 (Γij ) ,
≤ Ck Θhs (v)|Ωs,i − Rh (v)|Ωd,j · nij kL2 (Γij ) , Γij ∈ Γsd ,

with constants C independent of h, H, q, i, j, and the diameters of Γij and Ωd,j .

Proof. Since µH contains the constant functions, ℓh satisfies on Γij ∈ Γdd , using (4.43) and (4.6), Z Z Z  h h h ℓ = w · nij = (R (v) − v) · n] = 0, ∂Ωd,j

Γij

Γij

and this implies the unique solvability of (4.45). A similar argument holds on Γij ∈ Γsd using (4.44) ˆ j associated with Ωd,j , and (3.24). In order to check (4.46), we pass to the reference subdomain Ω −1 ˆ ˆ ˆ j be its exterior unit normal vector, Γ = Fj (Γij ), qˆ = q ◦ Fj and ℓ = ℓh ◦ Fj = (wh ◦ Fj ) · n ˆj; let n 1 2 ˆ ˆ then qˆ ∈ H (Ωj ) ∩ L0 (Ωj ) is the unique solution of ˆ ˆj. ˆ qˆ = 0 in Ω ˆ j , ∂ qˆ = Aj ℓˆ on ∂ Ω ∆ ˆj ∂ˆn

ˆ j ), it follows from [43] that qˆ belongs to H 23 (Ω ˆ j ) and As ℓˆ is in L2 (∂ Ω kˆ qk

3

ˆj) H 2 (Ω

ˆ 2 ˆ , ˆ j kℓk ≤ CA L (∂ Ωj )

ˆ j . Then (4.46) is a direct consequence with a constant Cˆ that only depends on the geometry of Ω of (4.43) and the following bounds: n−1

n

ˆ 2 ˆ , |q|H 1 (Ω ) = A 2 kℓh kL2 (Γij ) = Aj 2 kℓk j d,j L (∂ Ωj )

−1

n−3

|ˆ q |H 1 (Ωˆ j ) , |q|

The argument for (4.47) is similar, using (4.44).

H

3 2 (Ωd,j )

= Aj 2 |ˆ q|

3

ˆj) H 2 (Ω

. 

1

Now define c = ∇ q in Ωd,j . Then c belongs to H(div; Ωd,j )∩H 2 (Ωd,j )n and div c = 0. Therefore Rh (c) is well defined [18] and satisfies the approximation property for divergence-free functions [48] 1 0 0.

For Γij ∈ Γsd , we modify the argument in Corollary 4.3 to use the full approximation properties (4.16) of S h or S˜h to obtain
kchj,Γij (v)kL2 (Ωd,j ) ≤ C hr kvk r+ 21 ˜ + hs kvkH s+1 (Ω˜ s,i ∪Ω˜ d,j ) , ˜ H

(Ωs,i ∪Ωd,j )

1 0 < r ≤ min(rd + 1, t − ), 0 ≤ s ≤ min(rs , t − 1). 2 A combination of the above inequalities completes the proof of (5.5).



Next, we need to approximate the functions for the pressure. For any q ∈ L2 (Ωi ), let P h q be its onto Wih := W h |Ωi ,

L2 (Ωi )-projection

(q − P h q, wh ) = 0, ∀wh ∈ Wih ,

satisfying the approximation property kq − P h qk0,Ωi ≤ C hr |q|H r (Ωi ) , 0 ≤ r ≤ ri + 1,

(5.7)

where ri is the polynomial degree in the space Wih : ri = rs − 1 in Ωs and ri = ld in Ωd . 5.1. Error estimates. From the error equation: (5.8) ∀v h ∈ Z h , ∀qh ∈ Wh , a(u − uh , v h ) = −(b(v h , p − q h ) + bsd (v h , λsd ) + bd (v h , λd ) + bs (v h , λs )), and Lemmas 3.3 and 3.4, we immediately derive  (5.9) |||u − uh |||X ≤ C inf |||u − v h |||X + vh ∈V h (q

where V (5.10)

h (q

d)

= {v ∈ V

h;

∀w ∈

Rhu = sup

vh ∈Z h

d)

W h , b(v, w)

=

R

Ωd

inf

wh ∈W h

kp − wh kW



+ CRhu ,

w qd },

1 h h h b (v , λ ) + b (v , λ ) + b (v , λ ) , s s sd sd d d h |||v |||X

is the consistency error at the interfaces. Similarly a simple variant of (5.8) with v h ∈ V h and Theorem 3.1 yield   (5.11) kp − ph kW ≤ C |||u − uh |||X + inf kp − wh kW + CRhp , wh ∈W h

34

where Rhp is given by (5.10) with Z h replaced by V h . As the bound on the approximation error of W h follows from (5.7) and Lemma 5.1 estimates the approximation error of V h (qd ), it suffices to H H H H H bound Rhp . Note that by virtue of (3.11), we have for all µH s ∈ Λs , µd ∈ Λd , µsd ∈ Λsd , Rhp = sup

(5.12)

vh ∈V h

1 h H h H h H b (v , λ − µ ) + b (v , λ − µ ) + b (v , λ − µ ) . s s sd sd d d s sd d |||v h |||X

1 The bound for bs (v h , λs − µH s ) is straightforward because v h is measured in H in each Ωs,i , and therefore its trace can be considered on each Γi,j . But deriving an optimal bound for the other terms is more intricate because v h is measured in H(div) in each Ωd,i and this only controls the 1 H trace of the normal component in H − 2 (∂Ωd,i ). Consequently, λsd − µH sd and λd − µd must belong 1 to H 2 (∂Ωd,i ). In 2-D, this is easily achieved by interpolating λsd and λd with a variant of the Scott and Zhang interpolant that is continuous at the vertices of the subdomains; a similar idea was used in Section 4.3. But in 3-D, this construction requires an additional assumption on the mortar grids in the Darcy region.

Hypothesis 5.1. For each Ωd,i in 3-D, the mortar grids TijH on all Γij ⊂ ∂Ωd,i \ ∂Ω are chosen such that their traces on the boundaries of Γij coincide. This assumption implies conformity of mortar grids along interface boundaries in each connected region in Ωd and it allows us to consider the space ΛH,c d , which is the subset of continuous functions H H H in Λd ∪ Λsd . On Γdd ∪ Γsd , let I be the Scott-Zhang interpolant in ΛH,c d . We assume that its definition depends only on function values on Γdd ∪ Γsd . On Γss , let I H be the L2 (Γij )-orthogonal projection operator in ΛH s . This operator has the following approximation properties, assuming sufficient smoothness of µ and µ (the index ⋆ stands for dd or sd): (5.13) ∀Γij ∈ Γ⋆ , ∀τ ∈ Γij ,

(5.14) ∀Γij ∈ Γss , ∀τ ∈ Γij ,

kµ − I H (µ)kH t (τ ) ≤ CH r−t |µ|H r (∆τ ) , H

r

kµ − I (µ)kL2 (τ ) ≤ CH |µ|H r (∆τ ) ,

0 ≤ t ≤ 1, t ≤ r ≤ r⋆ + 1,

0 ≤ r ≤ rss + 1,

where ∆τ is the macroelement used in defining I H (µ) restricted to τ . The union of all ∆τ for τ in ∂Ωi may slightly overlap interfaces that are not part of ∂Ωi ; we denote the overlap by O. Lemma 5.2. There exists a constant C independent of h, H, and the diameters of the subdomains such that, for all v h ∈ X h : (5.15)

|bs (v h , λs − I H (λs ))| ≤ C H s

Ms X i=1

−1

Ai 2 kv h kH 1 (Ωs,i ) |λs |H s (∂Ωs,i ∩Γss ∪O) ,

0 ≤ s ≤ rss + 1,

provided λs is sufficiently smooth. Proof. We have h

H

bs (v , λs −I (λs )) =

X Z

Γij ∈Γss

h

Γij

H

[v ]·(λs −I (λs )) ≤

Ms X i=1

kv h kL2 (∂Ωs,i ∩Γss ) kλs −I H (λs )kL2 (∂Ωs,i ∩Γss ) .

Then (5.15) follows readily from (5.14) and the trace inequality (4.51).



35

Lemma 5.3. Under Hypothesis 5.1, there exists a constant C independent of h, H, and the diameters of the subdomains such that, for all v h ∈ X h : |bsd (v h , λsd −I H (λsd )) + bd (v h , λd − I H (λd ))| ≤C (5.16)

Md X i=1

+

Md X X i=1

provided λsd to Ωd,i .

1

1

kv h kH(div;Ωd,i ) H q− 2 |λ|H q (∂Ωd,i ∩Γdd ∪O) + H r− 2 |λ|H r (∂Ωd,i ∩Γsd ∪O)

j




 −1 Aj 2 kv h kH 1 (Ωs,j ) H r |λ|H r (∂Ωd,i ∩Γsd ∪O) ,

1 1 ≤ q ≤ rdd + 1, ≤ r ≤ rsd + 1, 2 2 and λd are sufficiently smooth, and where the sum over j runs over all Ωs,j adjacent

Proof. In the following λ denotes either λd or λsd depending on the type of interface. We can write Z X h H h H bsd (v , λ − I (λ))+bd (v , λ − I (λ)) = [v h · n](λ − I H (λ)) Γij

Γij ∈Γsd ∪Γdd

=

Md Z X i=1

h

∂Ωd,i

H

v · ni (λ − I (λ)) +

X Z

Γij ∈Γsd

Γij

v h · ns (λ − I H (λ)),

where we have used that v h · n = 0 on ∂Ωd,i ∩ ∂Ω and λ − I H (λ) has been extended continuously 1 1 from H 2 (∂Ωd,i \ ∂Ω) to H 2 (∂Ωd,i ). The argument of Lemma 5.2 can be used to bound the second sum above by the last sum in (5.16). Thus it is left to bound the first sum. Consider first one subdomain Ωd,i that is not adjacent to Γsd . Let us switch to the reference ˆ i and let E ˆ denote the extension operator defined in (4.29) relative to Ω ˆ i . For any f in domain Ω 1 1 ˆ H 2 (∂Ωd,i ), define E(f ) by: E(f ) = E(f ◦ Fi ); therefore E maps H 2 (∂Ωd,i ) into H 1 (Ωd,i ). Thus, by Green’s formula Z Z h H v · ni (λ − I (λ)) = v h · ni E(λ − I H (λ)) ∂Ωd,i

∂Ωd,i

=

Z

Ωd,i

div v h E(λ − I H (λ)) +

Z

Ωd,i

v h · ∇ E(λ − I H (λ)).

By reverting to the reference domain and observing that I H is invariant under the rigid body motion Fi , this reads Z ˆ ˆ ˆ ˆ v h · ni (λ − I H (λ)) ≤ Ain−1 kˆ v kH(div; ˆ i). ˆ Ω ˆ i ) kE(λ − I(λ))kH 1 (Ω ∂Ωd,i

On one hand, the bound in (4.29) gives

ˆ − I( ˆ ˆ λ ˆ λ))k ˆ ˆ ˆ ˆ kE( ˆ i ) ≤ Ckλ − I(λ))k H 1 (Ω ≤ Cˆ



X

ˆ ˆ i \∂ Ω ˆ Γ⊂∂ Ω

1

ˆ i) H 2 (∂ Ω

ˆ − I( ˆ ˆ λ ˆ λ))k ≤ Ck

ˆ − I( ˆ 2 ˆ λ))k kλ

H

1

2

1 ˆ 2 (Γ)

1

ˆ i \∂ Ω) ˆ H 2 (∂ Ω

,

with constants here and below independent of h, H, and the diameters of Ωd,i and Γij . By space ˆ and H 1 (Γ), ˆ the estimates (5.13) summed on Γ ˆ with t = 0 and t = 1 interpolation between L2 (Γ)

36

readily imply that ˆ − I( ˆ ˆ λ))k kλ

(5.17)

H

1 ˆ 2 (Γ)

r− 12

ˆ ≤ Cˆ H i

ˆ H r (∆ ) , |λ| ˆ Γ

1 ≤ r ≤ rdd + 1, 2

ˆ i on ∂ Ω ˆ i is related to the mesh size Hi on ∂Ωi by where the mesh size H ˆ i. Hi = Ai H ˆ we choose the In the case when r is not an integer, since (5.17) is stated in the reference region Γ, fractional seminorm in the right-hand side to be defined by (1.8) and incorporate the equivalence ˆ The motivation for this choice is that the intrinsic norm is easily transformed by constant into C. Fi . Therefore ˆ − I( ˆ ˆ λ))k kλ

1

H

1 ˆ i) 2 (∂ Ω

1

ˆ r ˆ ˆ ˆ ≤ Cˆ H ˆ r− 2 |λ| ˆ r− 2 Ar− ≤ Cˆ H i i i H (∂ Ωi \∂ Ω ∪O)

n−1 2

|λ|H r (∂Ωd,i \∂Ω ∪O) .

On the other hand, −n 2

kˆ v kH(div; ˆ Ω ˆ i ) = Ai



A2i kdiv v h k2L2 (Ωd,i ) + kv h k2L2 (Ωd,i )

1 2

.

By collecting these inequalities, we derive Z
1 ˆ i r− 2 kv h kH(div;Ω ) |λ|H r (∂Ω \∂Ω ∪O) , 1 ≤ r ≤ rdd + 1. v h · ni (λ − I H (λ)) ≤ Cˆ Ai H d,i d,i ∂Ωd,i 2

By applying (5.17) to a portion of Γsd , using (5.13) the same argument handles the case when Ωd,i is adjacent to Γsd . Then (5.16) follows easily from these estimates. 

Remark 5.1. We can skip Hypothesis 5.1, work locally on each Γij , and use the L2 norm for both
1 factors v h · ni and λ − I H (λ). But this gives rise to a factor of the form Hi /h 2 in (5.16), that is clearly not optimal when h is very small compared to Hi . Let Rh denote either Rhu or Rhp . The next lemma estimates Rh .

Lemma 5.4. Under Hypothesis 5.1, the consistency error Rh satisfies  1 1 1 1 + A− 2 H q− 2 kpd k q+ 12 Rh ≤ C A− 2 H r− 2 kpd k r+ 21 H

H

(Ωd )

−1

+A

(5.18)

s

(Ωd )

H (kus k

3 H s+ 2 (Ωs )

+ kps k

1 H s+ 2 (Ωs )

 ) ,

1 1 ≤ q ≤ rdd + 1, ≤ r ≤ rsd + 1, 0 < s ≤ rss + 1, 2 2 where A = min1≤i≤M Ai . Proof. By combining (5.15) and (5.16), we easily derive h



R ≤C H

q− 21

Md X i=1

|λ|2H q (∂Ωd,i ∩Γdd ∪O)

+ Hs

Ms X i=1


1 2

+H

r− 21

Md X i=1

2 A−1 i |λs |H s (∂Ωs,i ∩Γss ∪O)


1 2

|λ|2H r (∂Ωd,i ∩Γsd ∪O) Ms X

+ Hr

i=1

−1

1

2

2 A−1 i |λ|H r (∂Ωs,i ∩Γsd ∪O)

1 1 ≤ q ≤ rdd + 1, ≤ r ≤ rsd + 1, 0 ≤ s ≤ rss + 1. 2 2 In view of the representation (2.47), local trace theorems yield: |λsd |H r (∂Ωd,i ∩Γsd ) ≤ CAi 2 kpd k


1

H r+ 2 (Ωd,i )

,


1  2 ,

37 −1

|λd |H q (∂Ωd,i ∩Γdd ) ≤ CAi 2 kpd k

1

H q+ 2 (Ωd,i )

− 12

|λs |H s (∂Ωs,i ∩Γss ) ≤ CAi

kus k

H

s+ 3 2

(Ωs,i )

+ kps k

,

H

1 s+ 2

(Ωs,i )




.

Then (5.18) follows from these bounds, that are easily extended to include O, and the fact that −1 . A−1  i 0 1 ˆ and α. Recall that | · | denotes the ˆj (ℓ) and note that D 2 is well-defined. The next lemma relates π Euclidean vector norm. Lemma 7.1. For each component ℓˆ of ˆℓ we have ˆ 2 ˆ ≤ 2|D 12 α|. kˆ πj (ℓ)k L (Γ)

(7.6)

Proof. Considering the support of the basis functions ϕˆn , we have ˆ 22 = kˆ πj (ℓ)k ˆ L (Γ)

Z

x ˆ2

2

(α1 ϕˆ1 + α2 ϕˆ2 ) + x ˆ0

+

N −2 Z x ˆn+1 X n=2

Z

x ˆN +1

(αn ϕˆn + αn+1 ϕˆn+1 )2

x ˆn

(αN −1 ϕˆN −1 + αN +1 ϕˆN +1 )2 .

x ˆN −1

By substituting the expression of these integrals and rearranging terms, we derive N −2  X ˆ n−1 + h ˆ n) ˆ0 + h ˆ 1 ) + α2 (h ˆ0 + h ˆ1 + h ˆ 2) + ˆ 2 2 ≤ 2 α2 (h αn2 (h kˆ πj (ℓ)k 2 ˆ L (Γ) 3 1 n=3  2 2 ˆ ˆ ˆ ˆ ˆN ) . + αN −1 (hN −2 + hN −1 + hN ) + αN (hN −1 + h

In view of the diagonal terms dn of D, this can be written N −1  ˆ X ˆ ˆ ˆ  2 2 hN −1 + hN ˆ 2 2 ≤ 2 α2 h0 + h1 + kˆ πj (ℓ)k α d + α . n 1 n N ˆ L (Γ) 3 3 n=2

But

ˆ N −1 + h ˆN ˆ0 + h ˆ1 ˆ0 + h ˆ1 h h h = 2d1 < 2dN . < 2d1 , and similarly ˆ 0 + 2h ˆ1 3 3 h

Hence (7.7)

N 1 X 2 2 ˆ αn dn , kˆ πj (ℓ)kL2 (Γ) ˆ ≤2 n=1

whence (7.6).

 1 2

1 2

ˆ relies on a bound for |D α|. But D α has also This lemma shows that an L2 bound for π ˆj (ℓ) the following expression 1 1 1
−1 1
D 2 α = I + D 2 KD − 2 D− 2 b . Therefore Lemma 7.1 implies that
ˆ 2 ˆ ≤ 2k I + D 12 KD − 21 −1 k2 |D − 21 b|, (7.8) kˆ πj (ℓ)k L (Γ)

43

where k · k2 denotes the matrix norm subordinated by the Euclidean norm, and more generally k · kp is the matrix norm subordinated by the lp vector norm. The next lemma gives a bound for 1 |D − 2 b|. Lemma 7.2. We have 1

|D − 2 b| ≤

(7.9)

√

ˆ 2ˆ. 3kℓk L (Γ)

Proof. By integrating the basis functions ϕˆn , we derive (7.10) 1

|D − 2 b|2 = ≤ +

N X

2 d−1 n bn

n=1 N −2 ˆ X ˆ0 + h ˆ1 + h ˆ2 ˆ0 + h ˆ1 ˆn 1h h hn−1 + h ˆ 22 ˆ 22 ˆ 22 kℓk k ℓk kℓk + + x0 ,ˆ x2 ) x0 ,ˆ x3 ) xn−1 ,ˆ xn+1 ) L (ˆ L (ˆ L (ˆ 3 d1 d2 dn n=3

 ˆ N −1 + h ˆN ˆ N −2 + h ˆ N −1 + h ˆN h h ˆ 22 ˆ 22 kℓk + k ℓk L (ˆ xN −2 ,ˆ xN +1 ) L (ˆ xN −1 ,ˆ xN +1 ) . dN −1 dN

Then proceeding as in Lemma 7.1, we obtain |D

− 12

ˆ 22 ˆ2 b| ≤ kℓk L (ˆ x0 ,ˆ x1 ) + kℓkL2 (ˆ x0 ,ˆ x2 ) + 2

whence (7.9). 1

N X

n=1

ˆ 22 ˆ2 ˆ2 kℓk L (ˆ xn−1 ,ˆ xn+1 ) + kℓkL2 (ˆ xN ,ˆ xN +1 ) + kℓkL2 (ˆ xN −1 ,ˆ xN +1 ) , 


− 1 −1

It remains to evaluate k I + D 2 KD 2 ˆ n: following constraint on the mesh length h

k2 . Following Crouzeix & Thom´ee, we introduce the

Hypothesis 7.1. There exist two constants c0 > 0 and γ > 0 independent of N such that (7.11)

∀n, m, 0 ≤ n, m ≤ N ,

ˆn h ≤ c0 γ |n−m| . ˆm h

This condition holds for quasi-uniform meshes, but it is also satisfied by much more general meshes. Proposition 7.1. If Hypothesis 7.1 holds with suitable constants c0 and γ, then there exists a constant C independent of N and j such that 1
−1 1 k2 ≤ C. (7.12) k I + D 2 KD − 2 Hence with the same constant C,

(7.13)

√ ˆ ˆ 2 , kˆ ∀ˆℓ ∈ L2 (Γ) πj (ˆℓ)kL2 (Γ) ˆ ≤ 2 3CkℓkL2 (Γ) ˆ .

ˆ j vanish Proof. The proof follows the ideas of [23] but is more complex because the functions of X ˆ whereas those of Λ ˆ do not. Let us prove convergence of the series at the end points of Γ ∞ X r=1

this will yield (7.14)

1

1

k I + D 2 KD − 2

1

1

kD 2 K r D − 2 k2 ;


−1

k2 ≤ 1 +

∞ X r=1

1

1

kD 2 K r D − 2 k2 .

44

As K has three diagonals, then K r has 2r + 1 diagonals. In addition, since the coefficients of these diagonals are non negative, this implies that 1

1

1

1

kD 2 K r D − 2 k1 ≤ (2r + 1)kD 2 K r D − 2 k∞ , and by interpolation, 1

Therefore (7.15)

1

1

1

k2 ≤

 d 1

1

kD 2 K r D − 2 k2 ≤ (2r + 1) 2 kD 2 K r D − 2 k∞ . 1 2

r

kD K D

− 21

n

sup |n−m|≤2r

dm

2

1

(2r + 1) 2 kKkr∞ .

Thus a sufficient condition for the series convergence is:  δ r  d 1 n 2 ≤c , (7.16) sup kKk∞ |n−m|≤2r dm

where c > 0 and 0 < δ < 1 are two constants independent of r. First, assuming Hypothesis 7.1, a simple but tedious computation gives for 0 < γ < 2 dn ≤ 2c0 (1 + γ + γ 2 )γ 2r . |n−m|≤2r dm

(7.17)

sup

Hence the series converges if γ < Min 2,

1  . kKk∞

It is easy to check that h ˆ0 + h ˆ1 h ˆ0 ˆN ˆN + h ˆ N −1 1
1
 h h kKk∞ = Max , 1+ , , 1+ . ˆ 0 + 2h ˆ1 h ˆ N + 2h ˆ N −1 2 ˆ0 + h ˆ1 + h ˆ2 2 ˆ N −2 + h ˆ N −1 + h ˆN h h h

To simplify the notation, set θ = c0 γ. Owing to Hypothesis 7.1, for θ ≥ 1, we have 
 1  θ θ2 θ2 1 1 + Max = 1 + . , (7.18) kKk∞ ≤ 2 2 + θ 1 + θ + θ2 2 1 + θ + θ2

Take for instance θ = 32 ; then (7.18) yields kKk∞ < 43 . Therefore the series converges for γ = 43  and in turn c0 = 98 . Then (7.13) is an immediate consequence of Lemmas 7.1 and 7.2. ˆ 2 for ˆℓ ∈ H 1 (Γ) ˆ 2 . First, we have the analogue of Lemma 7.1. Now we estimate π ˆj (ˆℓ) in H 1 (Γ) 0

Lemma 7.3. We keep the notation of Lemma 7.1. Under Hypothesis 7.1 with γ > 1, we have r 2 c0 γ 2
12 − 1 ˆ c γ(1 + γ) + 2 + (7.19) |ˆ πj (ℓ)|H 1 (Γ) ≤ |D 2 α|. 0 ˆ 3 1+γ ˆ with a prime. Then Proof. Let us denote the derivation on Γ ˆ′ = π ˆj (ℓ)

N X

αn ϕˆ′n ,

n=1

and a straightforward computation gives

N −2  1
1 1
X 2 1 1
′ 2 2 1 2 ˆ αn kˆ πj (ℓ) kL2 (Γ) + α2 + + + ˆ ≤ 2 α1 ˆ + ˆ ˆ ˆ ˆ ˆ ˆ h0 h1 h0 + h1 h2 hn−1 hn n=3
1 1 1
 1 2 2 + α + + . + αN N ˆ −1 ˆ ˆ N −1 + h ˆN ˆN hN −2 h hN −1 h

45

Hypothesis 7.1 yields c0 γ 2
1 1 1 1 1 1 1 2 1 c0 γ(1 + γ) + 2 + , (1 + c0 γ) , (1 + c0 γ), + ≤ + ≤ + ≤ ˆ ˆ ˆ ˆ ˆ ˆ ˆ 2d 3d 1 + γ 3d 1 2 n h0 h1 h0 + h1 h2 hn−1 hn 1 ˆ N −2 h

+

1 ˆ N −1 + h ˆN h

≤

1 3dN −1

c0 γ(1 + γ) + 2 +

1 1 c0 γ 2
1 + ≤ , (1 + c0 γ). ˆ ˆ 1+γ 2dN hN −1 hN

Then (7.19) follows readily from the fact that for γ > 1, c0 γ(1 + γ) + 2 +

c0 γ 2 > 2(1 + c0 γ). 1+γ 

Comparing with (7.8) and the proof of Lemma 7.2, we see that a direct estimate of π ˆj (ˆℓ)′ relies 3 on a bound for |D − 2 b|, which we can hardly expect to be sharp. This difficulty can be bypassed ˆ 2 and arguing as in [23]. Uniqueness of the solution of by using the fact that ˆℓ belongs to H01 (Γ) ˆ j implies that π ˆ ˆℓ)) = I( ˆ ˆℓ) where Iˆ denotes the standard Lagrange interpolant in X ˆj , (4.24) in X ˆj (I( ˆ is a line segment. This permits to write that is well-defined because Γ ′

ˆ ˆℓ))′ + (I( ˆ ˆℓ) − ˆℓ)′ + ˆℓ . π ˆj (ˆℓ)′ = π ˆj (ˆℓ − I(

(7.20) First, we easily derive that

ˆ ˆℓ) − ˆℓ| 1 ˆ ≤ 2|ˆℓ| 1 ˆ . |I( H (Γ) H (Γ) Therefore (7.21)

ˆ ˆ ˆℓ))| 1 ˆ , |ˆ πj (ˆℓ)|H 1 (Γ) πj (ˆℓ − I( ˆ ≤ 3|ℓ|H 1 (Γ) ˆ + |ˆ H (Γ)

and it suffices to derive a bound for this last term. Let M α = b be the system (4.24) for a generic ˆ ˆℓ). Applying (7.12) and the analogue of (7.8), we obtain component of ˆℓ − I( 1
−1 3 3 −1 ˆ ˆ ℓ))| (7.22) |ˆ πj (ℓˆ − I( k2 |D − 2 b| ≤ C|D − 2 b|, ˆ ≤ k I + D 2 KD 2 H 1 (Γ)

with the constant C of (7.12), provided c0 and γ are chosen as in the proof of Proposition 7.1. Here 1 3 1 we use the fact that (7.15) is also valid for kD − 2 K r D 2 k2 . It remains to estimate |D − 2 b|. ˆ 2 . If Hypothesis 7.1 holds with the constants c0 and γ Proposition 7.2. Let ˆℓ belong to H01 (Γ) chosen as in the proof of Proposition 7.1, then each component ℓˆ of ˆℓ satisfies √ 3 ˆ 1ˆ, (7.23) |D − 2 b| < 30c|ℓ| H (Γ) ˆ Therefore, where c depends only on the interpolant I. (7.24)

ˆ 2 , |ˆ ∀ˆℓ ∈ H01 (Γ) πj (ˆℓ)|H 1 (Γ) ˆ < 3+

with the constants c and C of (7.23) and (7.12).

√


30cC |ˆℓ|H 1 (Γ) ˆ ,

ˆ Arguing as in Lemma 7.2, we recover ˆ ℓ). Proof. We sketch the proof. To simplify, set w = ℓˆ − I( 3 a bound for |D − 2 b| by replacing dn by d3n and ℓˆ by w in (7.10). The key point here is that the properties of the interpolant Iˆ imply
ˆ 2 |ℓ| ˆ2 ˆ 2 |ℓ| ˆ2 1 kwk2L2 (ˆx0 ,ˆx1 ) ≤ c h + h 0 H 1 (ˆ 1 x0 ,ˆ x1 ) H (ˆ x1 ,ˆ x2 ) ,

46

for a constant c independent of N , with analogous expressions for the norms in the other subinterˆ n in |D − 23 b| can be bounded by applying (7.11) with c0 γ = 3 , as in vals. The factors involving h 2 the proof of Proposition 7.1. This gives 3

ˆ2 ˆ2 1 |D − 2 b|2 < 10c |ℓ| x0 ,ˆ x2 ) + H (ˆ x0 ,ˆ x1 ) + |ℓ|H 1 (ˆ and (7.23) follows.

N X

n=1


ˆ2 ˆ2 ˆ2 1 |ℓ| xN −1 ,ˆ xN +1 ) , xN ,ˆ xN +1 ) + |ℓ|H 1 (ˆ H (ˆ xn−1 ,ˆ xn+1 ) + |ℓ|H 1 (ˆ



Finally, space interpolation between (7.13) and (7.24) gives 1

(7.25)

2 ˆ 2 (Γ) , |ˆ πj (ˆℓ)| ∀ˆℓ ∈ H00

1

2 (Γ) ˆ H00

ˆ ˆℓ| < C|

1

2 (Γ) ˆ H00

,

with a constant Cˆ independent of N . ˆ j have the same dimension, but of course this ˆ and X 7.1.2. Second example. In the first example, Λ is not necessary. For the unique solvability of (3.8), it is sufficient that the set of degrees of freedom ˆ j , sufficiently rich to guarantee a good approximation ˆ be an adequate subset of those of X of Λ ˆ ˆ with roughly property of the Lagrange interpolant I. Let us briefly give another simple choice of Λ half as many degrees of freedom as in the first example. Let N be an odd integer, keep the same ϕˆ0 and ϕˆN +1 as in (7.1), and choose 1 1 ϕˆ1 (ˆ x)|[ˆx0 ,ˆx1 ] = (ˆ x−x ˆ0 ) , ϕˆ1 (ˆ x)|[ˆx1 ,ˆx2 ] = (ˆ x −x ˆ), ˆ0 ˆ1 2 h h 1 1 (ˆ x−x ˆN −1 ) , ϕˆN (ˆ x)|[ˆxN ,ˆxN +1 ] = (ˆ xN +1 − x ˆ), ϕˆN (ˆ x)|[ˆxN −1 ,ˆxN ] = ˆ ˆ hN −1 hN (7.26) 1 ϕˆ2n (ˆ x)|[ˆx2n−2 ,ˆx2n ] = (ˆ x−x ˆ2n−2 ), ˆ 2n−2 + h ˆ 2n−1 h N 1 (ˆ x2n+2 − x ˆ), 1 ≤ n ≤ , ˆ ˆ 2 h2n + h2n+1 extended by zero elsewhere. We take ϕˆ2n (ˆ x)|[ˆx2n ,ˆx2n+2 ] =

(7.27)

ˆ j = {ˆ ˆ (ˆ ˆ 1 ϕˆ1 (ˆ X v ∈ H01 (0, 1)2 ; v x) = v x) +

N/2 X

ˆ (ˆ ˆ N ϕˆN (ˆ v x2n )ϕˆ2n (ˆ x) + v x)},

n=1

(N +1)/2

(7.28)

ˆ = {µ ˆ ∈ H 1 (0, 1)2 ; µ(ˆ ˆ x) = Λ

X

ˆ x2n )ϕˆ2n (ˆ µ(ˆ x)}.

n=0

With this choice, (4.24) is uniquely solvable and it can be checked that all the results established for the first example carry over to this second example, with different constants. 7.2. The 3-D case. In addition to restricting the mesh in 3-D (see Remark 4.1), the above proofs are complex because they apply to a situation where the matrix of the system is irreducible. The proofs are much easier, when (3.8) represents local and independent conditions, but in general, this can only be achieved by suitably modifying the discrete velocity spaces. In this section we present an example for which (3.8) holds by explicitly constructing the solution of (4.3). The correction satisfies suitable continuity bounds that are needed to establish the stability estimate (3.25). We follow the approach of BenBelgacem [12] who presents a local construction in 3-D by h in Ω ′ h adding to the space Xs,j s,j a stabilizing bubble function in each face T of the trace Tj,Γij of the triangulation Tjh on Γij and choosing for ΛH constant vectors on each face T ′ . More precisely, for each j, 1 ≤ j ≤ Ms , and for each T in Tjh , let P(T ) denote a polynomial space such as the

47

mini-element or Bernardi-Raugel element, used in approximating the velocity of the Stokes problem h , let T be the element of T h with face T ′ , and define with order one. For each face T ′ of Tj,Γ j ij bT ′ | T = λ 1 λ 2 λ 3 ,

where λk , k = 1, 2, 3, denote the barycentric coordinates of the three vertices of T ′ ; it is extended by zero outside T . Then set (7.29)

h Xs,j = {v ∈ C 0 (Ωs,j )3 ;∀T ∈ Tjh , T not adjacent to Γij , v|T ∈ P(T ),

∀T ∈ Tjh , T adjacent to Γij , v|T ∈ P(T ) + bT ′ IR3 },

and choose (7.30)

h ΛH = {µH ∈ L2 (Γij )3 ; ∀T ′ ∈ Tj,Γ , µH |T ′ ∈ IP 30 }. ij

As bT ′ vanishes at all vertices of T , it does not change the approximating properties of P(T ). Note that ΛH does not satisfy Hypothesis 3.2 (3). Nevertheless a discrete Korn inequality holds in Vsh , see Propositions 7.3 and 7.4 below. With this choice, we can show that (3.8) holds. Indeed, for ˜ it is easy to see that v ∈ X, Z X 1 vh = v T ′ bT ′ , where v T ′ = R v ′ T′ T ′ bT h ′ T ∈Tj,Γ

satisfies (3.8). We now define X (7.31) chj,Γij (v) =

h T ′ ∈Tj,Γ

ij

cT ′ bT ′ , where cT ′ ij

1 =R ′ T ′ bT

Z

T′


Θhs (v)|Ωs,i − S h (v)|Ωs,j ,

where Θhs (v)|Ωs,i is computed in the preceding subsection by correcting S h (v)|Ωs,i with chi,Γij (v) defined in (4.17). Lemma 7.4. The correction chj,Γij (v) defined by (7.31) satisfies (4.3) and there exists a constant C independent of h and the diameter of Ωs,i , Ωs,j , and Γij such that for all v in H01 (Ω)n (7.32)

|chj,Γij (v)|H 1 (Ωs,j ) ≤ C|v|H 1 (Ωs,i ∪Ωs,j ) ,

kchj,Γij (v)kL2 (Ωs,j ) ≤ Ch|v|H 1 (Ωs,i ∪Ωs,j ) .

h , we can write Proof. For any face T ′ of Tj,Γ ij Z
1 Θhs (v)|Ωs,i − v − (S h (v)|Ωs,j − v) . cT ′ = R ′ T′ T ′ bT

But on Γij ,

Θhs (v)|Ωs,i = S h (v)|Ωs,i + chi,Γij (v), owing to the support of each correction. Therefore Z
1 cT ′ = R chi,Γij (v) + (S h (v)|Ωs,i − v) − (S h (v)|Ωs,j − v) = A1 + A2 + A3 . ′ T′ T ′ bT

Let T be the element of Tjh adjacent to T ′ . By arguing as in the proof of Lemma 4.1, we easily derive ˆ H 1 (∆ ) , (7.33) |A3 bT ′ |H 1 (T ) ≤ C|v| T

S h (v)

Tjh

required to define in T . The estimation of A2 is similar, where ∆T is the macro-element of but more technical because T ′ belongs to Tjh whereas S h (v)|Ωs,i is constructed on Tih . Following the argument used in the proof of Lemma 4.3, we derive X |v|H 1 (∆ ) , (7.34) |A2 bT ′ |H 1 (T ) ≤ Cˆ Tℓ

ℓ

48

where {Tℓ } denote the set of elements of Tih that intersect T . The bound for A1 follows the same lines. With the notation of Lemma 4.1, Z X 1 c k bk , A1 = R ′ T′ T ′ bT k

where the sum runs over the indices k for which Ok ∩ T ′ 6= ∅. According to Hypothesis 4.1, the number of terms in this sum is bounded by a fixed integer. Thus X −1 Cˆ ρk 2 |v|H 1 (Dk ) , |A1 | ≤ k

and Hypothesis 4.1 implies that

ˆ |A1 bT ′ |H 1 (T ) ≤ C|v| ˜ T ), H 1 (∆

(7.35)

ℓ

˜ T is a macro-element in Ωs,i and again the number of elements in ∆ ˜ T is bounded by a where ∆ ℓ ℓ fixed integer. All constants are independent of h and the diameter of Ωs,i , Ωs,j , and Γij . Finally, the first inequality in (7.32) follows readily by summing (7.33)–(7.35) and applying Hypothesis 4.1. The argument for the second inequality in (7.32) is similar.  7.2.1. A discrete Korn inequality in Vsh . Since all assumptions except Hypothesis 3.2 (3) hold, it suffices to examine the jump of functions of Vsh through the interfaces Γij ∈ Γss . The next two lemmas show that the projection of this jump on polynomials of IP 1 is small. Lemma 7.5. Let Γij ∈ Γss , i < j. There exists a constant C, independent of h and the diameters of Γij , Ωs,i , and Ωs,j such that Z
1 3 [v h ] · p ≤ Ch 2 |Γij | 2 |v h |H 1 (Ds,i ) + |v h |H 1 (Ds,j ) , (7.36) ∀v h ∈ Vsh , ∀p ∈ IP 31 , Γij

where Ds,k denotes the layer of elements in Ωs,k adjacent to Γij .

Proof. All constants in this proof are independent of h and the diameters of Γij , Ωs,i , and Ωs,j . h , which is the mesh of ΛH . There is no loss of generality Let T ′ be an element (i.e. a face) of Tj,Γ ij in assuming that T ′ lies on the x1 − x2 plane. Consider a generic component v h of v h . Since the jump already satisfies Z [v h ] = 0,

T′

it suffices to bound the product of the jump by xk , k = 1, 2, say x1 . Then for any constant c, Z Z Z

1 h h x1 . [v ]x1 = [v ] − c x1 − ′ |T | T ′ T′ T′ A scaling argument gives

Z

ˆ ′ | 21 hT ′ .

x1 − 1 x1 ≤ C|T

′ |T | T ′ L2 (T ′ ) Similarly,
ˆ T ′ |[v h ]|H 1 (T ′ ) ≤ Ch ˆ T ′ |v h |Ω |H 1 (T ′ ) + |v h |Ω |H 1 (T ′ ) . inf k[v h ] − ckL2 (T ′ ) ≤ Ch s,i

c∈IR

s,j

On one hand, equivalence of norms gives

1

− |v h |Ωs,j |H 1 (T ′ ) ≤ Cˆ ρT 2 |v h |Ωs,j |H 1 (T ) ,

where T is the element of Tjh adjacent to T ′ . On the other hand, arguing as in the proof of Lemma 7.4, we prove that −1 |v h |Ω |H 1 (T ′ ) ≤ Cˆ ρ 2 |v h |Ω |H 1 (∆ ) , s,i

T

s,i

T

49

where ∆T is the union of all elements of Tih that intersect T ′ . Then taking the product and summing h , we obtain over all faces T ′ of Tj,Γ ij Z X

1 ˆ 32 |T ′ | 2 |v h |H 1 (Ds,i ) + |v h |H 1 (Ds,j ) , [v h ]x1 ≤ Ch Γij h ′ T ∈Tj,Γ

ij

whence (7.36).



Lemma 7.6. Let Γij ∈ Γss , i < j, and let Π([v h ]) ∈ IP 31 denote the projection of [v h ] onto IP 31 on Γij . There exists a constant C, independent of h and the diameters of Γij , Ωs,i , and Ωs,j such that
3 (7.37) ∀v h ∈ Vsh , kΠ([v h ])kL2 (Γij ) ≤ Ch 2 |v h |H 1 (Ds,i ) + |v h |H 1 (Ds,j ) . Proof. By definition Π([v h ]) ∈ IP 31 solves Z 3 ∀p ∈ IP 1 ,

h

Γij

Π([v ]) · p =

Z

Γij

[v h ] · p.

ˆ j , this reads By passing to the reference subdomain Ω Z Z 3 h −1 ˆ ˆ = |Γ| |Γij | ∀ˆ p ∈ IP 1 , Π([v ]) ◦ Fj · p ˆ Γ

Γij

[v h ] · p.

The coefficients of the polynomial of degree one Π([v h ]) ◦ Fj solve a linear system whose matrix ˆ Therefore, in view of (7.36), only depends on Γ.
1 h ˆ 3 h kΠ([v h ])kL2 (Γij ) ≤ |Γij | 2 kΠ([v h ]) ◦ Fj kL∞ (Γ) ˆ ≤ Ch 2 |v |H 1 (Ds,i ) + |v |H 1 (Ds,j ) .



This last result enables us to prove Korn’s inequality in Vsh . Proposition 7.3. Let |Γs | > 0 and Ωs be connected. Then there exists h0 > 0 such that for all h ≤ h0 , ∀v h ∈ Vsh ,

(7.38)

Ms X i=1

|v h |2H 1 (Ωs,i ) ≤ C

Ms X i=1

kD(v h )k2L2 (Ωs,i ) ;

the constants C and h0 are independent of h and the diameters of the interfaces Γij and the subdomains Ωs,k . Proof. Let v h ∈ Vsh . Formula (1.12) in [16] gives, with a constant C1 that depends only on the shape regularity of the mesh of subdomains TΩ , Ms X i=1

|v h |2H 1 (Ωs,i ) ≤ C1 ≤ C1

Ms X i=1

Ms X

kD(v

i=1

h

kD(v h )k2L2 (Ωs,i ) +

)k2L2 (Ωs,i )

+ Cˆ

X i<j

X i<j

 1 kΠ([v h ])k2L2 (Γij ) diam(Γij )


 1 3 h 2 h 2 h |v |H 1 (Ωs,i ) + |v |H 1 (Ωs,j ) , diam(Γij )

owing to (7.37). But h < diam(Γij ) and there exists an h0 > 0 such that ˆ −1 ≤ 1 . ∀h ≤ h0 , h2 (C1 C) 2

50

Since C1 and Cˆ are independent of h and the diameters of the interfaces Γij and the subdomains, then so is h0 . Hence for all h ≤ h0 , Ms X i=1

|v h |2H 1 (Ωs,i )

≤ 2C

Ms X i=1

kD(v h )k2L2 (Ωs,i ) . 

By arguing as in the proof of Lemma 3.4, we easily derive Korn’s inequality when |Γs | = 0 and Ωs is connected. Proposition 7.4. Let |Γs | = 0 and Ωs be connected, i.e. Γsd = ∂Ωs . Then there exists h0 > 0 such that for all h ≤ h0 , ! Ms Ms 2 Z 2 X X X ; |v hs · τ l | kD(v h )k2L2 (Ωs,i ) + (7.39) ∀v h ∈ Vsh , |v h |2H 1 (Ωs,i ) ≤ C i=1

i=1

l=1

Γsd

the constants C and h0 are independent of h and the diameters of the interfaces Γij and the subdomains Ωs,k . The case when Ωs is not connected follows from these two propositions applied to each connected component of Ωs according that it is or not adjacent to Γs . References

[1] R. A. Adams, Sobolev Spaces, Academic Press, New York, 1975. [2] P. Angot, On the well-posed coupling between free fluid and porous viscous flows, Appl. Math. Lett., 24 (2011), pp. 803–810. [3] T. Arbogast, Analysis of a two-scale, locally conservative subgrid upscaling for elliptic problems, SIAM J. Numer. Anal., 42 (2004), pp. 576–598. [4] T. Arbogast and D. S. Brunson, A computational method for approximating a Darcy-Stokes system governing a vuggy porous medium, Comput. Geosci., 11 (2007), pp. 207–218. [5] T. Arbogast, L. C. Cowsar, M. F. Wheeler, and I. Yotov, Mixed finite element methods on non-matching multiblock grids, SIAM J. Numer. Anal., 37 (2000), pp. 1295–1315. [6] T. Arbogast, G. Pencheva, M. F. Wheeler, and I. Yotov, A multiscale mortar mixed finite element method, Multiscale Model. Simul., 6 (2007), pp. 319–346. [7] I. Babuˇ ska, The finite element method with Lagrangian multipliers, Numer. Math., 20 (1973), pp. 179–192. [8] L. Badea, M. Discacciati, and A. Quarteroni, Numerical analysis of the Navier-Stokes/Darcy coupling, Numer. Math., 115 (2010), pp. 195–227. [9] S. Badia and R. Codina, Unified stabilized finite element formulations for the Stokes and the Darcy problems, SIAM J. Numer. Anal., 47 (2009), pp. 1971–2000. [10] G. S. Beavers and D. D. Joseph, Boundary conditions at a naturally impermeable wall, J. Fluid. Mech, 30 (1967), pp. 197–207. [11] F. Ben Belgacem, The mixed mortar finite element method for the incompressible Stokes problem: convergence analysis, SIAM J. Numer. Anal., 37 (2000), pp. 1085–1100 (electronic). [12] , A stabilized domain decomposition method with nonmatching grids for the Stokes problem in three dimensions, SIAM J. Numer. Anal., 42 (2004), pp. 667–685 (electronic). [13] C. Bernardi, Y. Maday, and A. T. Patera, A new nonconforming approach to domain decomposition: the mortar element method, in Nonlinear partial differential equations and their applications. Coll`ege de France Seminar, Vol. XI (Paris, 1989–1991), vol. 299 of Pitman Res. Notes Math. Ser., Longman Sci. Tech., Harlow, 1994, pp. 13–51. [14] C. Bernardi, T. C. Rebollo, F. Hecht, and Z. Mghazli, Mortar finite element discretization of a model coupling Darcy and Stokes equations, M2AN Math. Model. Numer. Anal., 42 (2008), pp. 375–410. [15] S. C. Brenner, Poincar´e–Friedrichs inequalities for piecewise H 1 functions, SIAM J. Numer. Anal., 41 (2003), pp. 306–324. [16] , Korn’s inequalities for piecewise H 1 vector fields, Math. Comp., 73 (2004), pp. 1067–1087. [17] F. Brezzi, On the existence, uniqueness and approximation of saddle-point problems arising from Lagrange multipliers, RAIRO, Anal. Num., 2 (1974), pp. 129–151. [18] F. Brezzi and M. Fortin, Mixed and hybrid finite element methods, Springer-Verlag, New York, 1991.

51

[19] E. Burman and P. Hansbo, A unified stabilized method for Stokes’ and Darcy’s equations, J. Comput. Appl. Math., 198 (2007), pp. 35–51. [20] Y. Cao, M. Gunzburger, X. Hu, F. Hua, X. Wang, and W. Zhao, Finite element approximations for Stokes-Darcy flow with Beavers-Joseph interface conditions, SIAM J. Numer. Anal., 47 (2010), pp. 4239–4256. [21] W. Chen, M. Gunzburger, F. Hua, and X. Wang, A parallel Robin-Robin domain decomposition method for the Stokes-Darcy system, SIAM J. Numer. Anal., 49 (2011), pp. 1064–1084. [22] Z. Chen and T. Y. Hou, A mixed multiscale finite element method for elliptic problems with oscillating coefficients, Math. Comp., 72 (2003), pp. 541–576. [23] M. Crouzeix and V. Thomee, The stability in Lp and Wp1 of the l2 projection onto finite element function spaces, Mathematics of Computation, 48 (1987), pp. 521–532. [24] M. Discacciati, E. Miglio, and A. Quarteroni, Mathematical and numerical models for coupling surface and groundwater flows, Appl. Numer. Math., 43 (2002), pp. 57–74. 19th Dundee Biennial Conference on Numerical Analysis (2001). [25] M. Discacciati and A. Quarteroni, Convergence analysis of a subdomain iterative method for the finite element approximation of the coupling of Stokes and Darcy equations, Comput. Vis. Sci., 6 (2004), pp. 93–103. [26] , Navier-Stokes/Darcy coupling: modeling, analysis, and numerical approximation, Rev. Mat. Complut., 22 (2009), pp. 315–426. [27] M. Discacciati, A. Quarteroni, and A. Valli, Robin-Robin domain decomposition methods for the StokesDarcy coupling, SIAM J. Numer. Anal., 45 (2007), pp. 1246–1268 (electronic). [28] V. J. Ervin, E. W. Jenkins, and S. Sun, Coupled generalized nonlinear Stokes flow with flow through a porous medium, SIAM J. Numer. Anal., 47 (2009), pp. 929–952. [29] J. Galvis and M. Sarkis, Non-matching mortar discretization analysis for the coupling Stokes-Darcy equations, Electron. Trans. Numer. Anal., 26 (2007), pp. 350–384. [30] , FETI and BDD preconditioners for Stokes-Mortar-Darcy systems, Commun. Appl. Math. Comput. Sci., 5 (2010), pp. 1–30. [31] B. Ganis and I. Yotov, Implementation of a mortar mixed finite element method using a multiscale flux basis, Comput. Methods Appl. Mech. Eng., 198 (2009), pp. 3989–3998. ´a, A conforming mixed finite-element method for the coupling of [32] G. N. Gatica, S. Meddahi, and R. Oyarzu fluid flow with porous media flow, IMA J. Numer. Anal., 29 (2009), pp. 86–108. ´a, and F.-J. Sayas, Analysis of fully-mixed finite element methods for the Stokes[33] G. N. Gatica, R. Oyarzu Darcy coupled problem, Math. Comp., 80 (2011), pp. 1911–1948. [34] V. Girault and P.-A. Raviart, Finite Element Methods for Navier-Stokes Equations, no. 5 in Springer series in Computational Mathematics, Springer-Verlag, Berlin, Heidelberg, New York, Tokyo, 1986. `re, DG approximation of coupled Navier–Stokes and Darcy equations by Beaver– [35] V. Girault and B. Rivie Joseph–Saffman interface condition, SIAM J. Numer. Anal., 47 (2009), pp. 2052–2089. [36] V. Girault and L. R. Scott, A quasi-local interpolation operator preserving the discrete divergence, Calcolo, (2003), pp. 1–19. [37] R. Glowinski and M. F. Wheeler, Domain decomposition and mixed finite element methods for elliptic problems, in First International Symposium on Domain Decomposition Methods for Partial Differential Equations, R. Glowinski, G. H. Golub, G. A. Meurant, and J. Periaux, eds., SIAM, Philadelphia, 1988, pp. 144–172. [38] P. Grisvard, Elliptic problems in nonsmooth domains, Pitman, Boston, 1985. [39] R. H. W. Hoppe, P. Porta, and Y. Vassilevski, Computational issues related to iterative coupling of subsurface and channel flows, Calcolo, 44 (2007), pp. 1–20. [40] T. Y. Hou and X. H. Wu, A multiscale finite element method for elliptic problems in composite materials and porous media, J. Comput. Phys., 134 (1997), pp. 169–189. [41] T. J. R. Hughes, Multiscale phenomena: Green’s functions, the Dirichlet-to-Neumann formulation, subgrid scale models, bubbles and the origins of stabilized methods, Comput. Methods Appl. Mech. Engrg., 127 (1995), pp. 387–401. ¨ger and A. Mikelic ´, On the interface boundary condition of Beavers, Joseph, and Saffman, SIAM J. [42] W. Ja Appl. Math., 60 (2000), pp. 1111–1127. [43] D. S. Jerrison and C. Kenig, The Neumann problem on Lipschitz domains, Bull. AMS, 4 (1981), pp. 203–207. `re, A strongly conservative finite element method for the coupling of Stokes and [44] G. Kanschat and B. Rivie Darcy flow, J. Comput. Phys., 229 (2010), pp. 5933–5943. [45] T. Karper, K.-A. Mardal, and R. Winther, Unified finite element discretizations of coupled Darcy-Stokes flow, Numer. Methods Partial Differential Equations, 25 (2009), pp. 311–326. [46] W. J. Layton, F. Schieweck, and I. Yotov, Coupling fluid flow with porous media flow, SIAM J. Numer. Anal., 40 (2003), pp. 2195–2218. [47] K. A. Mardal, X.-C. Tai, and R. Winther, A robust finite element method for Darcy-Stokes flow, SIAM J. Numer. Anal., 40 (2002), pp. 1605–1631 (electronic).

52

[48] T. P. Mathew, Domain decomposition and iterative refinement methods for mixed finite element discretizations of elliptic problems, PhD thesis, Courant Institute of Mathematical Sciences, New York University, 1989. Tech. Rep. 463. [49] M. Mu and J. Xu, A two-grid method of a mixed Stokes-Darcy model for coupling fluid flow with porous media flow, SIAM J. Numer. Anal., 45 (2007), pp. 1801–1813. ˇas, Les M´ethodes directes en th´eorie des ´equations elliptiques, Masson, Paris, 1967. [50] J. Nec [51] J. Peetre, Espaces d’interpolation et th´eor`eme de soboleff, Ann. Inst. Fourier, 16 (1966), pp. 279–317. [52] G. Pencheva and I. Yotov, Balancing domain decomposition for mortar mixed finite element methods, Numer. Linear Algebra Appl, 10 (2003), pp. 159–180. [53] G. Pencheva and I. Yotov, Interior superconvergence in mortar and non-mortar mixed finite element methods on non-matching grids, Comput. Methods Appl. Mech. Engrg., 197 (2008), pp. 4307–4318. `re and I. Yotov, Locally conservative coupling of Stokes and Darcy flows, SIAM J. Numer. Anal., 42 [54] B. Rivie (2005), pp. 1959–1977. [55] P. G. Saffman, On the boundary condition at the surface of a porous media, Stud. Appl. Math., L, (1971), pp. 93–101. [56] L. R. Scott and S. Zhang, Finite element interpolation of nonsmooth functions satisfying boundary conditions, Math. Comp., 54 (1990), pp. 483–493. [57] D. Vassilev and I. Yotov, Domain decomposition for coupled Stokes and Darcy flows. Preprint. [58] , Coupling Stokes-Darcy flow with transport, SIAM J. Sci. Comput., 31 (2009), pp. 3661–3684. [59] X. Xie, J. Xu, and G. Xue, Uniformly-stable finite element methods for Darcy-Stokes-Brinkman models, J. Comput. Math., 26 (2008), pp. 437–455. [60] I. Yotov, Mixed finite element methods for flow in porous media, PhD thesis, Rice University, Houston, Texas, 1996. TR96-09, Dept. Comp. Appl. Math., Rice University and TICAM report 96-23, University of Texas at Austin. (Vivette Girault) UPMC–Paris 6 and CNRS, UMR 7598, F-75230 Paris Cedex 05, France and Department of Mathematics, TAMU, College Station, TX 77843, USA; email: [email protected] (Danail Vassilev) Mathematics Research Institute, College of Engineering, Mathematics and Physical Sciences, University of Exeter, Exeter, EX4 4QF, UK; email: [email protected] (Ivan Yotov) Department of Mathematics, University of Pittsburgh, 301 Thackeray Hall, Pittsburgh, Pennsylvania 15260, USA; email: [email protected]