A UNIFORMLY STABLE FORTIN OPERATOR FOR THE TAYLOR ...
A UNIFORMLY STABLE FORTIN OPERATOR FOR THE TAYLOR–HOOD ELEMENT ¨ KENT–ANDRE MARDAL, JOACHIM SCHOBERL, AND RAGNAR WINTHER
Abstract. We construct a new Fortin operator for the lowest order Taylor–Hood element, which is uniformly stable both in L2 and H 1 . The construction, which is restricted to two space dimensions, is based on a tight connection between a subspace of the Taylor– Hood velocity space and the lowest order Nedelec edge element. General shape regular triangulations are allowed for the H 1 –stability, while some mesh restrictions are imposed to obtain the L2 –stability. As a consequence of this construction, a uniform inf–sup condition associated the corresponding discretizations of a parameter dependent Stokes problem is obtained, and we are able to verify uniform bounds for a family of preconditioners for such problems, without relying on any extra regularity ensured by convexity of the domain.
1. Introduction The stability of mixed finite element methods for the linear Stokes problem is intimately connected to the existence of bounded Fortin operators, i.e., bounded interpolation operators which properly commute with the divergence operator. For some discretizations the construction of locally defined Fortin operators is well–known. However, for the celebrated Taylor–Hood element the direct construction of a commuting interpolation operator is not obvious, and, as a consequence, most of the stability proofs found in the literature typically use alternative approaches. The main contribution of the present paper is to construct a new Fortin operator for the lowest order Taylor–Hood element in two space dimensions which, under proper conditions, is uniformly bounded both in L2 and H 1 . In fact, by utilizing a tight connection between a subspace of the Taylor–Hood velocity space and the lowest order Nedelec edge element, we perform a construction which resembles the corresponding well–known construction for the Mini element, cf. [1]. The discussions of this paper is motivated by the study of finite element discretizations of a singular perturbation problem related to the linear Stokes problem. More precisely, let Ω ⊂ Rn be a bounded Lipschitz domain and ∈ (0, 1] a real parameter. We will consider Date: June 7, 2012. 1991 Mathematics Subject Classification. 65N22, 65N30. Key words and phrases. Fortin operator, parameter dependent Stokes problem, uniform preconditioners. The first author was supported by the Research Council of Norway through Grant 209951 and a Centre of Excellence grant to the Centre for Biomedical Computing at Simula Research Laboratory. The third author was supported by the Research Council of Norway through a Centre of Excellence grant to the Centre of Mathematics for Applications. 1
¨ KENT–ANDRE MARDAL, JOACHIM SCHOBERL, AND RAGNAR WINTHER
2
singular perturbation problems of the form (I − 2 ∆)u − grad p = f
in Ω,
div u = g
in Ω,
(1.1)
u = 0 on ∂Ω, where the unknowns u and p are a vector field and a scalar field, respectively. For each fixed positive value of the perturbation parameter the problem behaves like the Stokes system, but formally the system approaches the mixed formulation of the scalar Laplace equation as this parameter tends to zero. In physical terms this means that we are studying fluid flow in regimes ranging from linear Stokes flow to porous medium flow. Another motivation for studying these systems is that they frequently arise as subsystems in time stepping schemes for time dependent Stokes and Navier–Stokes equations, cf. for example [2, 6, 17, 19, 20]. Therefore, the study of preconditioners for these systems is of substantial interest. Some preliminary discussions are presented in Sections 2, while the construction and analysis of the new Fortin operator will be carried out in Section 3. Finally, in Section 4 we will relate our results to the properties of a family of preconditioners for the coefficient operators of the discrete systems obtained from the Taylor–Hood method applied to (1.1). We will show that these preconditioners behave uniformly well with respect to the model parameter , as well as the discretization parameter. For the present model similar results have been derived before, but only by utilizing extra regularity ensured by convexity of the domain. 2. Preliminaries If X is a Hilbert space, then k · kX denotes its norm. We will use H m = H m (Ω) to denote the Sobolev space of functions on Ω with m derivatives in L2 = L2 (Ω). The corresponding spaces for vector fields are denoted H m (Ω; Rn ) and L2 (Ω; Rn ). Furthermore, h·, ·i is used to denote the inner–products in both L2 (Ω) and L2 (Ω; Rn ). In general, we will use H0m to denote the closure in H m of the space of smooth functions with compact support in Ω, while L20 is the space of L2 functions with mean value zero. We will use L(X, Y ) to denote the space of bounded linear operators mapping elements of X to Y , and if Y = X we simply write L(X) instead of L(X, X). If X and Y are Hilbert spaces, both continuously contained in some larger Hilbert space, then the intersection X ∩ Y and the sum X + Y are both Hilbert spaces with norms given by kxk2X∩Y = kxk2X + kxk2Y
and kzk2X+Y =
inf
x∈X,y∈Y z=x+y
(kxk2X + kyk2Y ).
The system (1.1) admits the following weak formulation: Find (u, p) ∈ H01 (Ω; Rn ) × L20 (Ω) such that (2.1)
hu, vi + 2 hDu, Dvi + hp, div vi = hf, vi, hdiv u, qi
= hg, qi,
v ∈ H01 (Ω; Rn ), q ∈ L20 (Ω)
A UNIFORMLY STABLE FORTIN OPERATOR
3
for given data f and g. Here Dv denotes the gradient of the vector field v. With the choice of H01 (Ω; Rn ) × L20 (Ω) as a solution space the bounds on the solution will blow up as tends to zero. To obtain a proper uniform bound on the solution operator we are forced to introduce dependent spaces and norms. A possible choice in the present setting is to take (u, p) ∈ (L2 ∩ H01 )(Ω; Rn ) × ((H 1 ∩ L20 ) + −1 L20 )(Ω), with more explicit norms given by kuk2L2 ∩ H 1 = kuk2L2 + 2 kDuk2L2
and kpk2H 1 +−1 L2 = hgrad(− 2 ∆ + I)−1 p, grad pi.
Here, a natural boundary condition is assumed for the inverse of the elliptic operator. Uniform bounds for the solution in these norms can be derived from the Brezzi conditions [4, 5]. In fact, the only nontrivial condition needed to obtain such bounds is the uniform inf–sup condition hdiv v, qi (2.2) sup ≥ αkqkH 1 +−1 L2 , q ∈ L20 (Ω), 2 1 1 kvk n L ∩ H v∈H0 (Ω;R ) where the positive constant α is independent of ∈ (0, 1]. This condition is closely tied to the existence of a proper right inverse of the divergence operator in L(L20 (Ω), H01 (Ω; Rn )). If there exist such a right inverse, which can be extended to an operator in L(H0−1 (Ω), L2 (Ω; Rn )), then the uniform condition (2.2) will follow. Here, H0−1 ⊃ L20 corresponds to the dual space of H 1 ∩ L20 , obtained by extending the L2 inner–product to a duality pairing. In other words, the uniform condition (2.2) will follow from the existence of an operator S with the properties (2.3)
S ∈ L(L20 , H01 ) ∩ L(H0−1 , L2 ),
and div S = I,
cf. [15] for more details. A proper operator which satisfy these conditions is the Bogovski˘i operator, see for example [10, Section III.3], or [8, 11]. On a domain Ω ⊂ Rn , which is star shaped with respect to an open ball B, this operator is explicitly given as an integral operator of the form Z Z ∞ z z g(y)K(x − y, y) dy, where K(z, y) = n Sg(x) = θ(y + r )rn−1 dr, |z| |z| |z| Ω where θ ∈ C0∞ (B) is positive, and with integral equal to one. This operator is a right inverse of the divergence operator, and it has exactly the desired mapping properties given by (2.3), cf. [8, 11]. Furthermore, the definition of S can be extended to general bounded Lipschitz domains, by using the fact that such domains can be written as a finite union of star shaped domains. The constructed operator will again satisfy the properties given by (2.3). We refer to [10, Section III.3], [11, Section 2], and [8, Section 4.3] for more details. Remark 2.1. An alternative right inverse of the divergence operator, mapping L2 (Ω) boundedly into H01 (Ω; Rn ), can be constructed from the Stokes problem, i.e., problem (1.1) with = 1. However, this operator can, in general, not be extended to an operator in L(H0−1 , L2 ), unless the solution operator for the Stokes problem admits extra regularity. Therefore, this right inverse will, in general, not lead to the uniform inf–sup condition (2.2). The construction of preconditioners for the coefficient operators of the discrete systems derived from (2.1) has been frequently studied, cf. for example [2, 6, 17, 18, 19, 20]. Several of these studies have shown that a family of preconditioners, described in Section 4 below, behave uniformly
¨ KENT–ANDRE MARDAL, JOACHIM SCHOBERL, AND RAGNAR WINTHER
4
well, both with respect to the perturbation parameter and the discretization parameter h. However, all the theoretical verifications of this result relies on additional properties of the domain which lead to extra regularity of the solution of Stokes problem. In fact, it seems that the non–optimal results of these studies partly can be traced to the use of the wrong right inverse in the continuous case. We will return to the discussion of preconditioners for discretizations of (2.1) in Section 4 below. Finally, we mention that an alternative derivation of the uniform inf–sup condition (2.2) for non–convex domains, based on domain decomposition, is presented in the recent manuscript [13]. We will consider finite element discretizations of the problem (1.1) of the form: Find (uh , ph ) ∈ Vh × Qh such that (2.4)
huh , vi + 2 hDuh , Dvi + hph , div vi = hf, vi, hdiv uh , qi
= hg, qi,
v ∈ Vh , q ∈ Qh .
Here Vh and Qh are finite element spaces such that Vh ×Qh ⊂ H01 (Ω; Rn )×L20 (Ω), and h is the discretization parameter. For the Taylor–Hood element, and the Mini element, the pressure space Qh is in fact a subspace of H 1 , and this will be assumed throughout our discussions below. The proper discrete uniform inf–sup condition, which will guarantee uniform stability of the method (2.4), is of the form hdiv v, qi (2.5) sup ≥ αkqkH 1 +−1 L2 , q ∈ Qh , v∈Vh kvkL2 ∩ H 1 where the positive constant α is independent of both and h. The technique we will use to establish the discrete inf–sup condition (2.5) is in principle rather standard. We will just rely on the corresponding continuous condition (2.2) and bounded interpolation operators into the spaces Vh and Qh which commute properly with the divergence operator. In the present setting, this means that we need to construct interpolation operators Πh : L2 (Ω; Rn ) → Vh which are uniformly bounded, with respect to and h, in the norm of L2 ∩ H01 . However, this is equivalent to the requirement that (2.6)
kΠh kL(L2 )
and kΠh kL(H 1 )
are bounded independently of h.
Furthermore, the operators Πh must satisfy the commuting relation (2.7)
hdiv Πh v, qi = hdiv v, qi v ∈ H01 (Ω; Rn ), q ∈ Qh .
Operators Πh , satisfying (2.6) and (2.7) will be referred to as a uniformly bounded Fortin operators. If such operators are constructed, the discrete uniform inf–sup condition (2.5) will follow from the corresponding continuous condition (2.2). In the case of the Mini element the construction of uniformly bounded Fortin operators, in any space dimension, is straightforward. In fact, the design of this element, originally presented in [1], was highly motivated by the desire to construct a suitable Fortin operator (cf. also [5, Chapter VI] and [17]). However, for the Taylor–Hood element the direct construction of a bounded Fortin operator is not obvious. In fact, most of the stability proofs found in the literature for this discretization typically use alternative approaches, cf. for example [5, Section VI.6], [12, Chapter II.4.2] and the discussion given in the introduction of [9]. An exception is [9], where a local interpolation operator satisfying (2.7) is constructed. However, this operator is not
A UNIFORMLY STABLE FORTIN OPERATOR
5
defined in L2 . Below we propose a new construction of an interpolation operator, satisfying (2.6) and (2.7), by utilizing a technique for the Taylor–Hood method which is similar to the standard construction for the Mini element. As for the construction carried out in [9], this analysis is restricted to two space dimensions.
3. A Fortin operator for the Taylor–Hood element This section is devoted to the construction of the new local Fortin operator Πh for the lowest order Taylor–Hood element. We will assume that the domain Ω ⊂ R2 is a polyhedral domain which is triangulated by a family of shape regular, triangular meshes {Th } indexed by decreasing values of the mesh parameter h = maxT ∈Th hT . Here hT is the diameter of the triangle T . We recall that the mesh is shape regular if there exists a positive constant γ0 such that for all values of the mesh parameter h h2T ≤ γ0 |T |,
T ∈ Th .
Here |T | denotes the area of T . For the lowest order Taylor–Hood element the velocity space, Vh , consists of continuous piecewise quadratic vector fields, while Qh is the standard space of continuous piecewise linear scalar fields with respect to the triangulation Th . For technical reasons we will assume throughout the analysis below that any T ∈ Th has at most one edge in ∂Ω. Such an assumption is frequently made for convenience when the Taylor–Hood element is analyzed, cf. for example [5, Proposition 6.1], since most approaches requires a special construction near the boundary. On the other hand, this assumption will not hold for many simple triangulations. Therefore, at the end of this section we will outline how to refine the analysis such that this assumption can be removed. We note that if we have established the discrete inf–sup condition for a pair of spaces (Vh− , Qh ), where Vh− ⊂ Vh , then this condition will also hold for the pair (Vh , Qh ). This observation will be utilized here. We let Vh− ⊂ Vh ⊂ H01 (Ω; R2 ) be the space of piecewise quadratic vector fields which has the property that on each edge of the mesh the normal components of elements in Vh− are linear. Each function v in the space Vh− can be determined from its values at each interior vertex of the mesh, and of the mean value of the tangential component along each interior edge. In fact, the space Vh− can be decomposed as Vh− = Vh1 ⊕ Vhb , where the space Vh1 is the space of continuous piecewise linear vector fields, while the space Vhb is spanned by quadratic “edge bubbles.” To define this space of bubbles we let ∆1 (Th ) = ∆i1 (Th ) ∪ ∆∂1 (Th ) be the set of the edges of the mesh Th , where ∆i1 (Th ) is the set of interior edges and ∆∂1 (Th ) is the set of edges on the boundary of Ω. Furthermore, if T ∈ Th then ∆1 (T ) is the set of edges of T , and ∆i1 (T ) = ∆1 (T ) ∩ ∆i1 (T ). For each e ∈ ∆1 (Th ) we let Ωe be the associated macroelement consisting of the union of the two triangles T ∈ Th with e ∈ ∆1 (T ). The scalar function be is the unique continuous and piecewise quadratic function on Ωe which vanish on R the boundary of Ωe , and with e be ds = |e|, where |e| denotes the length of e. The space Vhb is defined as Vhb = span{be te | e ∈ ∆i1 (Th ) },
6
¨ KENT–ANDRE MARDAL, JOACHIM SCHOBERL, AND RAGNAR WINTHER
where te is a tangent vector along e with length |e|. Alternatively, if xi and xj are the vertices corresponding to the endpoints of e then the vector field ψe = be te is determined up to a sign as ψe = 6λi λj (xj − xi ), where {λi } are the piecewise linear functions corresponding to the barycentric coordinates, i.e., λi (xk ) = δi,k for all vertices xk . In particular, Z Z ψe · (xj − xi ) ds = 6 λi λj ds|e|2 = |e|3 . e
e
The construction of the desired interpolation operator Πh utilizes a standard Clement operator Rh mapping into the space Vh1 , cf. [7]. The operator Rh : L2 (Ω; R2 ) → Vh1 is local, it preserves constants, and it is stable in L2 and H01 . More precisely, we have for any T ∈ Th that (3.1)
k(I − Rh )vkH j (T ) ≤ chk−j T kvkH k (ΩT ) ,
0 ≤ j ≤ k ≤ 1,
where the constant c is independent of h and v. Here ΩT denote the macroelement consisting of T and all elements T 0 ∈ Th such that T ∩ T 0 6= ∅. The interpolation operator Πh will be of the form Πh = Πhb (I − Rh ) + Rh , where Πhb : L2 (Ω; R2 ) → Vhb will be constructed below. A key part of the construction is that the operator Πhb , mapping into the bubble space Vhb , satisfies an analog of (2.7), i.e., (3.2)
hdiv Πhb v, qi = hdiv v, qi,
q ∈ Qh .
If this identity holds then we will also have that hdiv(I − Πh )v, qi = hdiv(I − Πhb )(I − Rh )v, qi = 0 for all q ∈ Qh , which is (2.7). Hence, it remains to construct the operator Πhb , satisfying (3.2), and such that the corresponding operator Πh satisfies (2.6). We will need to separate the triangles which have an edge on the boundary of Ω from the interior triangles. With this purpose we define Th∂ = {T ∈ Th | T ∩ ∂Ω ∈ ∆∂1 (Th ) } and Thi = Th \ Th∂ . In order to define the operator Πhb we introduce Zh as the lowest order Nedelec space with respect to the mesh Th . Hence, if z ∈ Zh then on any T ∈ Th , z is a linear vector field such that z(x) · x is also linear. Furthermore, for each e ∈ ∆1 (Th ) the tangential component of z is continuous. As a consequence, Zh ⊂ H(curl; Ω) where the operator curl denotes the two dimensional analog of the curl–operator given by curl z = curl(z1 , z2 ) = ∂x2 z1 − ∂x1 z2 . It is well known that the proper degrees of freedom for the space Zh is the mean value of the tangential components of v, v · t, with respect to each edge in ∆1 (Th ). Furthermore, we let Z 0 (3.3) Zh = {z ∈ Zh | z · t ds = 0, T ∈ Th∂ }. ∂T
Zh0
Alternatively, the elements of are those vector fields in Zh with the property that curl z|T = 0 if T ∈ Th∂ , i.e., z is a constant vector field on T for T ∈ Th∂ . It is a key observation that the mesh assumption given above, that any T ∈ Th intersects ∂Ω in at most
A UNIFORMLY STABLE FORTIN OPERATOR
7
one edge, implies that the spaces Vhb and Zh0 have the same dimension. Furthermore, we note that grad q ∈ Zh0 for any q ∈ Qh . For each e ∈ ∆1 (Th ) let φe ∈ Zh be the basis function corresponding to the Whitney form, i.e., φe satisfies Z Z (φe · te ) ds = |e|, and (φe · te0 ) ds = 0, e0 6= e, e0
e
where, as above, te is a tangent vector of length e. Hence, if e = (xi , xj ) then the vector field φe can be expressed in barycentric coordinates as φe = λi grad λj − λj grad λi . Any z ∈ Zh0 can be written uniquely on the form X X z= ae φe + e∈∆i1 (Th )
ce φe ,
e∈∆∂1 (Th )
where the coefficients ae corresponding to interior edges can be chosen arbitrarily, but where the coefficients ce for each boundary edge should be chosen such that curl z = 0 on the associated triangle in Th∂ . We note that there is a natural mapping Φh between the spaces Vhb and Zh0 given by Φh (ψe ) = φe for all interior edges, or alternatively, Z Z −2 v · te ds, e ∈ ∆i1 (Th ). Φh (v) · te ds = |e| (3.4) e
e
Vhb (T )
Zh0 (T )
Below, we will use and to denote the restriction of the spaces Vhb and Zh0 to a single element T , and ΦT will denote the corresponding restriction of the map Φh . We will define the operator Πhb : L2 (Ω; Rn ) → Vhb by, hΠhb u, zi = hu, zi
(3.5)
∀z ∈ Zh0 .
To show that this operator is well–defined, and behaves properly, the following general formula for integration of products of barycentric coordinates over a triangle T will be useful (cf. for example [14, Section 2.13]), Z 2α! (3.6) λα1 1 λα2 2 λα3 3 dx = |T |, (2 + |α|)! T P where α! = α1 !α2 !α3 ! and |α| = i αi . Lemma 3.1. Let T ∈ Thi with edges e1 , e2 , e3 . For any v ∈ Vhb (T ) we have Z 1 2 3 |a| |T | ≤ v · ΦT (v) dx ≤ |a|2 |T |, 5 5 T P P where ΦT = Φh |T , v = i ai ψei , and |a|2 = i a2i . Proof. A direct computation gives Z Z Z 3 X 3 3 X 3 X X v · ΦT (v) dx = ai aj ψei · φej dx = ai aj bei φej · tei dx = aT M a, T
i=1 j=1
T
i=1 j=1
T
8
¨ KENT–ANDRE MARDAL, JOACHIM SCHOBERL, AND RAGNAR WINTHER
where the 3 × 3 matrix M is given by Z Z 3 3 M = {Mi,j }i,j=1 = { ψi · φej dx}i,j=1 { bei φej · tei dx}3i,j=1 , T
T
The desired result will follow from the diagonal dominance of this matrix. Let xj be the vertex opposite to ej , and let λj be the corresponding barycentric coordinate on T . Then the first diagonal element M1,1 of the matrix M is given by Z Z M1,1 = 6 λ2 λ3 (λ2 grad λ3 − λ3 grad λ2 )(x3 − x2 ) dx = 6 (λ22 λ3 + λ2 λ23 ) dx = 2|T |/5, T
T
where we have used formula (3.6) in the final step. Actually, from this formula we derive that all the diagonal elements are given by Mi,i = 2|T |/5, and similar calculations for the off–diagonal elements gives |Mi,j | = |T |/10. In addition, the matrix M is symmetric. The matrix M is therefore strictly diagonally dominant, and by the Gershgorin circle theorem all eigenvalues are bounded below by |T |/5. By combining this with the fact that M is symmetric we conclude that aT M a ≥ |a|2 |T |/5, and this is the desired bound. The next lemma is a variant of the result above for T ∈ Th∂ . Lemma 3.2. Let T ∈ Th∂ with interior edges e1 , e2 , and where e3 is the edge on the boundary. For any v ∈ Vhb (T ) we have Z 1 v · ΦT (v) dx = |a|2 |T |, 2 T 2 2 2 where v = a1 ψe1 + a2 ψe2 and |a| = a1 + a2 . Proof. Let v = a1 ψe1 +a2 ψe2 = 6a1 λ2 λ3 (x3 −x2 )+6a2 λ1 λ3 (x3 −x1 ). It is a key observation that in this case ΦT (v) is simply given as Φ(v) = −a1 grad λ2 − a2 grad λ1 . Since grad λi · tei ≡ 0 we therefore obtain from (3.6) that Z Z Z 2 2 v · ΦT (v) dx = 6a1 λ2 λ3 grad λ2 · (x2 − x3 ) dx + 6a2 λ1 λ3 grad λ1 · (x1 − x3 ) dx T T ZT Z 1 = 6a21 λ2 λ3 dx + 6a22 λ1 λ3 dx = |a|2 |T |. 2 T T This completes the proof. P P P Let v = e∈∆i (Th ) ae ψe ∈ Vhb such that z = Φh (v) = e∈∆i (Th ) ae φe + e∈∆∂ (Th ) ce φe ∈ Zh0 . 1 1 1 It follows from scaling and shape regularity that the two norms X X (3.7) kzkL2 (Ω) and ( a2e )1/2 , T ∈Th e∈∆i1 (T )
are equivalent, uniformly with respect to h. Correspondingly, the two norms X X (3.8) kvkL2 (Ω) and ( |T |2 a2e )1/2 , T ∈Th
e∈∆i1 (T )
are uniformly equivalent. In fact, these equivalences hold locally for each T ∈ Th .
A UNIFORMLY STABLE FORTIN OPERATOR
9
It will be convenient to introduce the weighted L2 –norms given by X X 2 1/2 kvk1,h = ( h−2 and kzk−1,h = ( h2T kzk2L2 (T ) )1/2 . T kvkL2 (T ) ) T ∈Th
T ∈Th
In fact, for functions v ∈ Vhb the Pnorm kvk1,h is equivalent, uniformly in h, to the ordinary 1 H –norm. Furthermore, P P if v = e ae ψe then the two norms kvk1,h and kφh (v)k−1,h are both equivalent to ( T |T | e a2e )1/2 . This observation leads to the following inf–sup property. Lemma 3.3. There is a positive constant c0 , independent of h, such that for each v ∈ Vhb sup z∈Zh0
hv, zi ≥ c0 kvk1,h . kzk−1,h
P b Proof. Let v = e∈∆i1 (Th ) ae ψe ∈ Vh be given. We simply choose the corresponding z = Φh (v). As a consequence of the norm equivalences (3.7) and (3.8), combined with Lemmas 3.1 and 3.2, we obtain XZ X 1 X hv, Φh (v)i = Φh (v) · v dx ≥ |T | a2e ≥ c0 kΦh (v)k−1,h kvk1,h , 5 T i T ∈T T ∈T h
e∈∆1 (T )
h
where c0 > 0 is independent of h. This completes the proof.
The inf–sup condition just derived implies that the operator Πhb satisfies the stability estimate (3.9)
kΠhb vk1,h ≤ c−1 0 kvk1,h ,
v ∈ L2 (Ω).
If the family {Th } is quasi–uniform, i.e., there exists a positive constant γ1 , independent of h and T ∈ Th , such that hT ≥ γ1 h, then the estimate (3.9) also implies that Πhb is uniformly bounded in the standard non weighted L2 –norm. An alternative local mesh restriction, which implies L2 –stability, is the following condition. There exists a constant δ ∈ [0, 1), independent of h, such that |Ωe0 | ≤ (1 + δ)|Ωe |,
(3.10)
for all edges e, e0 ∈ ∆i1 (Th ) sharing a common vertex. The following consequence of this condition will be useful. Lemma 3.4. Assume that the local mesh condition (3.10) holds. There is a constant ∈ [0, 1/3), and for each h and T ∈ Th a positive constant αT , such that |1 −
|Ωe | | ≤ , αT
e ∈ ∆i1 (T ).
Proof. For each T ∈ Th let αT = (1 + δ/2) min |Ωe |. i e∈∆1 (T )
10
¨ KENT–ANDRE MARDAL, JOACHIM SCHOBERL, AND RAGNAR WINTHER
Due to condition (3.10) this implies that |Ωe | ≤ (1 + δ)(1 + δ/2)−1 αT for all e ∈ ∆i1 (T ). Therefore, we can conclude that 1 |Ωe | 1+δ ≤ ≤ , 1 + δ/2 αT 1 + δ/2 which is equivalent to the desired bound with = δ/(2 + δ) < 1/3.
Lemma 3.5. Assume that the local mesh condition (3.10) holds. There is a positive constant c1 , independent of h, such that for each v ∈ Vhb sup z∈Zh0
hv, zi ≥ c1 kvkL2 (Ω) . kzkL2 (Ω)
P Proof. We let and αT be as in the statement of Lemma 3.4. For a given v = e∈∆i (Th ) ae ψe ∈ 1 Vhb we choose z = Eh Φh (v), where Eh : Zh0 → Zh0 is defined by Eh φe = |Ωe |φe . As a consequence the L2 –norm of z is equivalent to the L2 –norm of v. It follows from these norm equivalences, combined with Lemmas 3.1 and 3.2, that XZ hv, zi = hv, Eh Φh (v)i = Eh Φh (v) · v dx T ∈Th
=
X T ∈Th
≥
X T ∈Th
T
Z
Z
(Eh − αT I)Φh (v) · v dx]
Φh (v) · v dx +
[αT
T
T
X 1 3 |Ωe | αT |T |[ − max (1 − )] a2e i 5 5 e∈∆1 (T ) αT i e∈∆1 (T )
X 1 3 X ≥( − ) αT |T | a2e ≥ c1 kvkL2 (Ω) kzkL2 (Ω) . 5 5 T ∈T i h
e∈∆1 (T )
This completes the proof.
It is a direct consequence of Lemma 3.5 that the operator norm of Πhb in L(L2 (Ω; R2 )) is bounded by c−1 1 . We now can conclude with the following main result of this section. Theorem 3.6. The operator Πh is uniformly bounded in H 1 and satisfies the commuting relation (2.7). Furthermore, if the triangulation is either quasi–uniform or it satisfies the local condition (3.10), then Πh is also uniformly bounded in L2 , i.e., (2.6) holds. Proof. To show that the operator Πh fulfills the condition (2.7), it is enough to show that the operator Πhb satisfies the corresponding condition (3.2). However, since grad Qh ⊂ Zh0 we have hdiv Πhb v, qi = −hΠhb v, grad qi = −hv, grad qi = hdiv v, qi.
A UNIFORMLY STABLE FORTIN OPERATOR
11
Furthermore, since Rh is bounded in H 1 , the desired H 1 boundedness of Πh will follow from the corresponding bound on Πhb (I − Rh ). However, by shape regularity, combined with (3.9) we have X X 2 kΠhb (I − Rh )vk2H 1 (Ω) ≤ c hT−2 kΠhb (I − Rh )vk2L2 (T ) ≤ c h−2 T k(I − Rh )vkL2 (T ) , T ∈Th
T ∈Th
and hence the H 1 –bound follows from (3.1). Finally, if the triangulation is quasi–uniform or satisfies (3.10) then the L2 –bound follows from the corresponding bounds for the operators Πhb and Rh . Remark 3.1. The two conditions guaranteeing L2 –stability of the Fortin operator Πh are rather different. Condition (3.10) is local, while quasi–uniformity is a global mesh condition. In some sense, this resembles two of the mesh conditions which can be used to obtain H 1 –stability of the L2 –projection onto continuous finite elements. It is a classical, and well–known, result that quasi–uniformity implies H 1 –stability of this operator, while an alternative local condition was derived in [3]. Remark 3.2. The analysis of the Taylor–Hood method given above leans heavily on the assumption that there are no triangles in Th with more than one edge on the boundary of Ω. This assumption simplifies the analysis, but it is not necessary. For more general triangulations, where some triangles may have two edges on the boundary, we will still define the interpolation operators Πh and Πhb as above. The only necessary modification is that the definition of the bubble space Vhb must be altered, since the dimension of the space Zh0 , given by (3.3), in general is larger than the number of interior edges. A possible choice is Vhb = span({ be te | e ∈ ∆i1 (Th ) } ∪ { be ne | e ∈ ∆i,2 1 (Th ) }), where te and ne are tangent and normal vectors to the edge e, and where ∆i,2 1 (Th ) is the i subset of ∆1 (Th ) consisting of the third edges of the triangles with two edges on the boundary. Hence, the space Vhb is enriched by including “the normal edge bubbles” associated the edges b 0 in ∆i,2 1 (Th ). With this definition the space Vh has the same dimension as Zh . In [15] the analysis above is extended to this more general case. However, the proof in this setting is more technical than the analysis above, and therefore it is omitted here. 4. Uniform preconditioners The interest in the construction of preconditioners for the discretizations of the singular perturbation system (1.1) is partly motivated by the role of such systems in time stepping schemes for the time dependent Stokes and Navier–Stokes equations. The system (1.1) has a coefficient operator which formally is of the form I − 2 ∆ − grad , A = div 0 with obvious weak interpretations of the differential operators derived from the weak formulation (2.1). In the notation of [16], the corresponding block diagonal operator B =
12
¨ KENT–ANDRE MARDAL, JOACHIM SCHOBERL, AND RAGNAR WINTHER
diag((I − 2 ∆)−1 , (−∆)−1 + 2 I), can be referred to as a canonical preconditioner for A , since the uniform inf–sup condition (2.2) implies that the composition B A is an isomorphism mapping the space (L2 ∩ H01 ) × ((H 1 ∩ L20 ) + −1 L20 ) to itself, with operator norms of B A and its inverse bounded independently of . We refer to [16, Section 4] and [17, Section 3] for more details. Furthermore, as explained in Section 2 above, if the Fortin operator Πh is uniformly bounded both in H 1 and L2 then the corresponding discrete inf–sup condition (2.5) holds for the Taylor–Hood method. Therefore, it follows from Theorem 3.6 that, if the triangulation {Th } is either quasi–uniform or satisfies the local condition (3.10), then (2.5) holds. In particular, the condition numbers of the corresponding discrete operators B,h A,h will be uniformly bounded with respect to and h, cf. [16, Section 5]. This result is obtained without any convexity assumption on the domain. Since all previous theoretical studies of uniform preconditioners for the model (1.1) have relied on convexity, cf. Remark 2.1 above, we will present a numerical experiment which confirm the theory above. More precisely, the experiment indicates that the condition number of B,h A,h is basically unaffected by such properties of the domain. The results below are obtained by using a canonical preconditioner, i.e., the operator B,h is composed of exact inverses of discretizations of the differential operators appearing in the definition of B . Of course, to obtain more practical preconditioners, these inverses should again be replaced by suitable elliptic preconditioners. We consider the problem (1.1) on three two dimensional domains, referred to as Ω1 , Ω2 and Ω3 . Here Ω1 is the unit square, Ω2 is the L–shaped domain obtained by cutting out the an upper right subsquare from Ω1 , while Ω3 is the slit domain where a slit of length a half is removed from from Ω1 , cf. Figure 1. Hence, only Ω1 is convex. The corresponding problems
Figure 1. The domains Ω1 , Ω2 , and Ω3 . (1.1) were discretized by the lowest order Taylor–Hood element on a uniform triangular grid. The condition numbers of the corresponding discrete operators B,h A,h , for the three domains and for different values of and h, are computed. The results are given in Table 1 below. These results clearly indicate that the condition numbers of the operators B,h A,h are not dramatically effected by lack of convexity of the domains. References [1] D.N. Arnold, F. Brezzi, and M. Fortin, A stable finite element method for Stokes equations, Calcolo 21 (1984), 337–344. [2] J.H. Bramble and J.E. Pasciak, Iterative techniques for time dependent Stokes problems, Comput. Math. Appl. 33 (1997), no. 1-2, 13–30, Approximation theory and applications. [3] J.H. Bramble, J.E. Pasciak, and O. Steinbach On the stability of the L2 projection in H 1 (Ω), Math. Comp. 71 (2001), 147–156.
A UNIFORMLY STABLE FORTIN OPERATOR
domain
Ω1
Ω2
Ω3
\h
2−2
2−3
2−4
2−5
1
14.3 15.3 18.5 22.0
0.1
10.3 11.9 12.7 13.2
0.01
6.2
6.7
8.0
9.9
1
17.1 17.2 17.1 17.1
0.1
10.3 11.8 12.7 13.2
0.01
6.1
6.7
8.0
13
9.9
1
13.2 13.4 13.5 13.6
0.1
10.3 11.9 12.8 13.2
0.01 6.1 6.7 8.1 9.9 Table 1. Condition numbers for the operators B,h A,h
[4] F. Brezzi, On the existence, uniqueness and approximation of saddle–point problems arising from Lagrangian multipliers, RAIRO. Analyse Num´erique 8 (1974), 129–151. [5] F. Brezzi and M. Fortin, Mixed and hybrid finite element methods, Springer-Verlag, 1991. [6] J. Cahouet and J.-P. Chabard, Some fast 3D finite element solvers for the generalized Stokes problem, Internat. J. Numer. Methods Fluids 8 (1988), no. 8, 869–895. [7] P. Clement, Approximation by finite element functions using local regularization, RAIRO Anal. Num´er. 9 (1975), 77–84. [8] M. Costabel and A. McIntosh, On Bogovski˘i and regularized Poincar´e integral operators for de Rham complexes on Lipschitz domains, http://arxiv.org/abs/0808.2614v1. [9] R.S Falk, A Fortin operator for two–dimensional Taylor–Hood elements, Mathematical Modelling and Numerical Analysis (M2AN) 42 (2008), 411-424. [10] Giovanni P. Galdi, An introduction to the mathematical theory of the Navier–Stokes equations, vol 1, Linearized Steady problems, vol. 38, Springer Tracts in Natural Philosophy, Springer–Verlag, New York 1994. [11] M. Geissert, H. Heck, and M. Hieber, On the equation div u = g and Bogovski˘i’s operator in Sobolev spaces of negative order, Oper. Theory Adv. Appl. 168 (2006), 113–121. [12] V. Girault and P.-A. Raviart, Finite element methods for Navier-Stokes equations, Theory and algorithms, Springer Series in Computational Mathematics, vol. 5, Springer-Verlag, Berlin, 1986. [13] T. Kolev and P. Vassilevski, Regular decomposition of H(div) spaces, Lawrence Livermore National Laboratory, 2011. [14] M.-J. Lay and L.L. Schumaker, Spline functions on triangulations, vol. 110 of Encyclopedia of Mathematics and its Applications, Cambridge University Press, 2007. [15] K.-A. Mardal, J. Sch¨ oberl, and R. Winther, A uniform inf–sup condition with applications to preconditioning, arXiv:1201.1513v1,2012. [16] K.-A. Mardal and R. Winther, Preconditioning discretizations of systems of partial differential equations, Numer. Linear Alg. Appl. 18 (2011), 1–40. [17] K.-A. Mardal and R. Winther, Uniform preconditioners for the time dependent Stokes problem, Numer. Math. 98 (2004), no. 2, 305–327. , Erratum: “Uniform preconditioners for the time dependent Stokes problem” [Numer. Math. 98 [18] (2004), no. 2, 305–327 ], Numer. Math. 103 (2006), no. 1, 171–172.
14
¨ KENT–ANDRE MARDAL, JOACHIM SCHOBERL, AND RAGNAR WINTHER
[19] M. A. Olshanskii, J. Peters, and A. Reusken, Uniform preconditioners for a parameter dependent saddle point problem with application to generalized Stokes interface equations, Numerische Mathematik 105 (2006), 159–191. [20] S. Turek, Efficient solvers for incompressible flow problems, Springer-Verlag, 1999.
Centre for Biomedical Computing at Simula Research Laboratory, Norway and Department of Informatics, University of Oslo, Norway E-mail address: [email protected] URL: http://simula.no/people/kent-and Institute for Analysis and Scientific Computing, Wiedner Hauptstrasse 8-10, 1040 Wien, Austria E-mail address: [email protected] URL: http://www.asc.tuwien.ac.at/~schoeberl Centre of Mathematics for Applications and Department of Informatics, University of Oslo, 0316 Oslo, Norway E-mail address: [email protected] URL: http://folk.uio.no/rwinther
Abstract. We construct a new Fortin operator for the lowest order Taylor–Hood element, which is uniformly stable both in L2 and H 1 . The construction, which is restricted to two space dimensions, is based on a tight connection between a subspace of the Taylor– Hood velocity space and the lowest order Nedelec edge element. General shape regular triangulations are allowed for the H 1 –stability, while some mesh restrictions are imposed to obtain the L2 –stability. As a consequence of this construction, a uniform inf–sup condition associated the corresponding discretizations of a parameter dependent Stokes problem is obtained, and we are able to verify uniform bounds for a family of preconditioners for such problems, without relying on any extra regularity ensured by convexity of the domain.
1. Introduction The stability of mixed finite element methods for the linear Stokes problem is intimately connected to the existence of bounded Fortin operators, i.e., bounded interpolation operators which properly commute with the divergence operator. For some discretizations the construction of locally defined Fortin operators is well–known. However, for the celebrated Taylor–Hood element the direct construction of a commuting interpolation operator is not obvious, and, as a consequence, most of the stability proofs found in the literature typically use alternative approaches. The main contribution of the present paper is to construct a new Fortin operator for the lowest order Taylor–Hood element in two space dimensions which, under proper conditions, is uniformly bounded both in L2 and H 1 . In fact, by utilizing a tight connection between a subspace of the Taylor–Hood velocity space and the lowest order Nedelec edge element, we perform a construction which resembles the corresponding well–known construction for the Mini element, cf. [1]. The discussions of this paper is motivated by the study of finite element discretizations of a singular perturbation problem related to the linear Stokes problem. More precisely, let Ω ⊂ Rn be a bounded Lipschitz domain and ∈ (0, 1] a real parameter. We will consider Date: June 7, 2012. 1991 Mathematics Subject Classification. 65N22, 65N30. Key words and phrases. Fortin operator, parameter dependent Stokes problem, uniform preconditioners. The first author was supported by the Research Council of Norway through Grant 209951 and a Centre of Excellence grant to the Centre for Biomedical Computing at Simula Research Laboratory. The third author was supported by the Research Council of Norway through a Centre of Excellence grant to the Centre of Mathematics for Applications. 1
¨ KENT–ANDRE MARDAL, JOACHIM SCHOBERL, AND RAGNAR WINTHER
2
singular perturbation problems of the form (I − 2 ∆)u − grad p = f
in Ω,
div u = g
in Ω,
(1.1)
u = 0 on ∂Ω, where the unknowns u and p are a vector field and a scalar field, respectively. For each fixed positive value of the perturbation parameter the problem behaves like the Stokes system, but formally the system approaches the mixed formulation of the scalar Laplace equation as this parameter tends to zero. In physical terms this means that we are studying fluid flow in regimes ranging from linear Stokes flow to porous medium flow. Another motivation for studying these systems is that they frequently arise as subsystems in time stepping schemes for time dependent Stokes and Navier–Stokes equations, cf. for example [2, 6, 17, 19, 20]. Therefore, the study of preconditioners for these systems is of substantial interest. Some preliminary discussions are presented in Sections 2, while the construction and analysis of the new Fortin operator will be carried out in Section 3. Finally, in Section 4 we will relate our results to the properties of a family of preconditioners for the coefficient operators of the discrete systems obtained from the Taylor–Hood method applied to (1.1). We will show that these preconditioners behave uniformly well with respect to the model parameter , as well as the discretization parameter. For the present model similar results have been derived before, but only by utilizing extra regularity ensured by convexity of the domain. 2. Preliminaries If X is a Hilbert space, then k · kX denotes its norm. We will use H m = H m (Ω) to denote the Sobolev space of functions on Ω with m derivatives in L2 = L2 (Ω). The corresponding spaces for vector fields are denoted H m (Ω; Rn ) and L2 (Ω; Rn ). Furthermore, h·, ·i is used to denote the inner–products in both L2 (Ω) and L2 (Ω; Rn ). In general, we will use H0m to denote the closure in H m of the space of smooth functions with compact support in Ω, while L20 is the space of L2 functions with mean value zero. We will use L(X, Y ) to denote the space of bounded linear operators mapping elements of X to Y , and if Y = X we simply write L(X) instead of L(X, X). If X and Y are Hilbert spaces, both continuously contained in some larger Hilbert space, then the intersection X ∩ Y and the sum X + Y are both Hilbert spaces with norms given by kxk2X∩Y = kxk2X + kxk2Y
and kzk2X+Y =
inf
x∈X,y∈Y z=x+y
(kxk2X + kyk2Y ).
The system (1.1) admits the following weak formulation: Find (u, p) ∈ H01 (Ω; Rn ) × L20 (Ω) such that (2.1)
hu, vi + 2 hDu, Dvi + hp, div vi = hf, vi, hdiv u, qi
= hg, qi,
v ∈ H01 (Ω; Rn ), q ∈ L20 (Ω)
A UNIFORMLY STABLE FORTIN OPERATOR
3
for given data f and g. Here Dv denotes the gradient of the vector field v. With the choice of H01 (Ω; Rn ) × L20 (Ω) as a solution space the bounds on the solution will blow up as tends to zero. To obtain a proper uniform bound on the solution operator we are forced to introduce dependent spaces and norms. A possible choice in the present setting is to take (u, p) ∈ (L2 ∩ H01 )(Ω; Rn ) × ((H 1 ∩ L20 ) + −1 L20 )(Ω), with more explicit norms given by kuk2L2 ∩ H 1 = kuk2L2 + 2 kDuk2L2
and kpk2H 1 +−1 L2 = hgrad(− 2 ∆ + I)−1 p, grad pi.
Here, a natural boundary condition is assumed for the inverse of the elliptic operator. Uniform bounds for the solution in these norms can be derived from the Brezzi conditions [4, 5]. In fact, the only nontrivial condition needed to obtain such bounds is the uniform inf–sup condition hdiv v, qi (2.2) sup ≥ αkqkH 1 +−1 L2 , q ∈ L20 (Ω), 2 1 1 kvk n L ∩ H v∈H0 (Ω;R ) where the positive constant α is independent of ∈ (0, 1]. This condition is closely tied to the existence of a proper right inverse of the divergence operator in L(L20 (Ω), H01 (Ω; Rn )). If there exist such a right inverse, which can be extended to an operator in L(H0−1 (Ω), L2 (Ω; Rn )), then the uniform condition (2.2) will follow. Here, H0−1 ⊃ L20 corresponds to the dual space of H 1 ∩ L20 , obtained by extending the L2 inner–product to a duality pairing. In other words, the uniform condition (2.2) will follow from the existence of an operator S with the properties (2.3)
S ∈ L(L20 , H01 ) ∩ L(H0−1 , L2 ),
and div S = I,
cf. [15] for more details. A proper operator which satisfy these conditions is the Bogovski˘i operator, see for example [10, Section III.3], or [8, 11]. On a domain Ω ⊂ Rn , which is star shaped with respect to an open ball B, this operator is explicitly given as an integral operator of the form Z Z ∞ z z g(y)K(x − y, y) dy, where K(z, y) = n Sg(x) = θ(y + r )rn−1 dr, |z| |z| |z| Ω where θ ∈ C0∞ (B) is positive, and with integral equal to one. This operator is a right inverse of the divergence operator, and it has exactly the desired mapping properties given by (2.3), cf. [8, 11]. Furthermore, the definition of S can be extended to general bounded Lipschitz domains, by using the fact that such domains can be written as a finite union of star shaped domains. The constructed operator will again satisfy the properties given by (2.3). We refer to [10, Section III.3], [11, Section 2], and [8, Section 4.3] for more details. Remark 2.1. An alternative right inverse of the divergence operator, mapping L2 (Ω) boundedly into H01 (Ω; Rn ), can be constructed from the Stokes problem, i.e., problem (1.1) with = 1. However, this operator can, in general, not be extended to an operator in L(H0−1 , L2 ), unless the solution operator for the Stokes problem admits extra regularity. Therefore, this right inverse will, in general, not lead to the uniform inf–sup condition (2.2). The construction of preconditioners for the coefficient operators of the discrete systems derived from (2.1) has been frequently studied, cf. for example [2, 6, 17, 18, 19, 20]. Several of these studies have shown that a family of preconditioners, described in Section 4 below, behave uniformly
¨ KENT–ANDRE MARDAL, JOACHIM SCHOBERL, AND RAGNAR WINTHER
4
well, both with respect to the perturbation parameter and the discretization parameter h. However, all the theoretical verifications of this result relies on additional properties of the domain which lead to extra regularity of the solution of Stokes problem. In fact, it seems that the non–optimal results of these studies partly can be traced to the use of the wrong right inverse in the continuous case. We will return to the discussion of preconditioners for discretizations of (2.1) in Section 4 below. Finally, we mention that an alternative derivation of the uniform inf–sup condition (2.2) for non–convex domains, based on domain decomposition, is presented in the recent manuscript [13]. We will consider finite element discretizations of the problem (1.1) of the form: Find (uh , ph ) ∈ Vh × Qh such that (2.4)
huh , vi + 2 hDuh , Dvi + hph , div vi = hf, vi, hdiv uh , qi
= hg, qi,
v ∈ Vh , q ∈ Qh .
Here Vh and Qh are finite element spaces such that Vh ×Qh ⊂ H01 (Ω; Rn )×L20 (Ω), and h is the discretization parameter. For the Taylor–Hood element, and the Mini element, the pressure space Qh is in fact a subspace of H 1 , and this will be assumed throughout our discussions below. The proper discrete uniform inf–sup condition, which will guarantee uniform stability of the method (2.4), is of the form hdiv v, qi (2.5) sup ≥ αkqkH 1 +−1 L2 , q ∈ Qh , v∈Vh kvkL2 ∩ H 1 where the positive constant α is independent of both and h. The technique we will use to establish the discrete inf–sup condition (2.5) is in principle rather standard. We will just rely on the corresponding continuous condition (2.2) and bounded interpolation operators into the spaces Vh and Qh which commute properly with the divergence operator. In the present setting, this means that we need to construct interpolation operators Πh : L2 (Ω; Rn ) → Vh which are uniformly bounded, with respect to and h, in the norm of L2 ∩ H01 . However, this is equivalent to the requirement that (2.6)
kΠh kL(L2 )
and kΠh kL(H 1 )
are bounded independently of h.
Furthermore, the operators Πh must satisfy the commuting relation (2.7)
hdiv Πh v, qi = hdiv v, qi v ∈ H01 (Ω; Rn ), q ∈ Qh .
Operators Πh , satisfying (2.6) and (2.7) will be referred to as a uniformly bounded Fortin operators. If such operators are constructed, the discrete uniform inf–sup condition (2.5) will follow from the corresponding continuous condition (2.2). In the case of the Mini element the construction of uniformly bounded Fortin operators, in any space dimension, is straightforward. In fact, the design of this element, originally presented in [1], was highly motivated by the desire to construct a suitable Fortin operator (cf. also [5, Chapter VI] and [17]). However, for the Taylor–Hood element the direct construction of a bounded Fortin operator is not obvious. In fact, most of the stability proofs found in the literature for this discretization typically use alternative approaches, cf. for example [5, Section VI.6], [12, Chapter II.4.2] and the discussion given in the introduction of [9]. An exception is [9], where a local interpolation operator satisfying (2.7) is constructed. However, this operator is not
A UNIFORMLY STABLE FORTIN OPERATOR
5
defined in L2 . Below we propose a new construction of an interpolation operator, satisfying (2.6) and (2.7), by utilizing a technique for the Taylor–Hood method which is similar to the standard construction for the Mini element. As for the construction carried out in [9], this analysis is restricted to two space dimensions.
3. A Fortin operator for the Taylor–Hood element This section is devoted to the construction of the new local Fortin operator Πh for the lowest order Taylor–Hood element. We will assume that the domain Ω ⊂ R2 is a polyhedral domain which is triangulated by a family of shape regular, triangular meshes {Th } indexed by decreasing values of the mesh parameter h = maxT ∈Th hT . Here hT is the diameter of the triangle T . We recall that the mesh is shape regular if there exists a positive constant γ0 such that for all values of the mesh parameter h h2T ≤ γ0 |T |,
T ∈ Th .
Here |T | denotes the area of T . For the lowest order Taylor–Hood element the velocity space, Vh , consists of continuous piecewise quadratic vector fields, while Qh is the standard space of continuous piecewise linear scalar fields with respect to the triangulation Th . For technical reasons we will assume throughout the analysis below that any T ∈ Th has at most one edge in ∂Ω. Such an assumption is frequently made for convenience when the Taylor–Hood element is analyzed, cf. for example [5, Proposition 6.1], since most approaches requires a special construction near the boundary. On the other hand, this assumption will not hold for many simple triangulations. Therefore, at the end of this section we will outline how to refine the analysis such that this assumption can be removed. We note that if we have established the discrete inf–sup condition for a pair of spaces (Vh− , Qh ), where Vh− ⊂ Vh , then this condition will also hold for the pair (Vh , Qh ). This observation will be utilized here. We let Vh− ⊂ Vh ⊂ H01 (Ω; R2 ) be the space of piecewise quadratic vector fields which has the property that on each edge of the mesh the normal components of elements in Vh− are linear. Each function v in the space Vh− can be determined from its values at each interior vertex of the mesh, and of the mean value of the tangential component along each interior edge. In fact, the space Vh− can be decomposed as Vh− = Vh1 ⊕ Vhb , where the space Vh1 is the space of continuous piecewise linear vector fields, while the space Vhb is spanned by quadratic “edge bubbles.” To define this space of bubbles we let ∆1 (Th ) = ∆i1 (Th ) ∪ ∆∂1 (Th ) be the set of the edges of the mesh Th , where ∆i1 (Th ) is the set of interior edges and ∆∂1 (Th ) is the set of edges on the boundary of Ω. Furthermore, if T ∈ Th then ∆1 (T ) is the set of edges of T , and ∆i1 (T ) = ∆1 (T ) ∩ ∆i1 (T ). For each e ∈ ∆1 (Th ) we let Ωe be the associated macroelement consisting of the union of the two triangles T ∈ Th with e ∈ ∆1 (T ). The scalar function be is the unique continuous and piecewise quadratic function on Ωe which vanish on R the boundary of Ωe , and with e be ds = |e|, where |e| denotes the length of e. The space Vhb is defined as Vhb = span{be te | e ∈ ∆i1 (Th ) },
6
¨ KENT–ANDRE MARDAL, JOACHIM SCHOBERL, AND RAGNAR WINTHER
where te is a tangent vector along e with length |e|. Alternatively, if xi and xj are the vertices corresponding to the endpoints of e then the vector field ψe = be te is determined up to a sign as ψe = 6λi λj (xj − xi ), where {λi } are the piecewise linear functions corresponding to the barycentric coordinates, i.e., λi (xk ) = δi,k for all vertices xk . In particular, Z Z ψe · (xj − xi ) ds = 6 λi λj ds|e|2 = |e|3 . e
e
The construction of the desired interpolation operator Πh utilizes a standard Clement operator Rh mapping into the space Vh1 , cf. [7]. The operator Rh : L2 (Ω; R2 ) → Vh1 is local, it preserves constants, and it is stable in L2 and H01 . More precisely, we have for any T ∈ Th that (3.1)
k(I − Rh )vkH j (T ) ≤ chk−j T kvkH k (ΩT ) ,
0 ≤ j ≤ k ≤ 1,
where the constant c is independent of h and v. Here ΩT denote the macroelement consisting of T and all elements T 0 ∈ Th such that T ∩ T 0 6= ∅. The interpolation operator Πh will be of the form Πh = Πhb (I − Rh ) + Rh , where Πhb : L2 (Ω; R2 ) → Vhb will be constructed below. A key part of the construction is that the operator Πhb , mapping into the bubble space Vhb , satisfies an analog of (2.7), i.e., (3.2)
hdiv Πhb v, qi = hdiv v, qi,
q ∈ Qh .
If this identity holds then we will also have that hdiv(I − Πh )v, qi = hdiv(I − Πhb )(I − Rh )v, qi = 0 for all q ∈ Qh , which is (2.7). Hence, it remains to construct the operator Πhb , satisfying (3.2), and such that the corresponding operator Πh satisfies (2.6). We will need to separate the triangles which have an edge on the boundary of Ω from the interior triangles. With this purpose we define Th∂ = {T ∈ Th | T ∩ ∂Ω ∈ ∆∂1 (Th ) } and Thi = Th \ Th∂ . In order to define the operator Πhb we introduce Zh as the lowest order Nedelec space with respect to the mesh Th . Hence, if z ∈ Zh then on any T ∈ Th , z is a linear vector field such that z(x) · x is also linear. Furthermore, for each e ∈ ∆1 (Th ) the tangential component of z is continuous. As a consequence, Zh ⊂ H(curl; Ω) where the operator curl denotes the two dimensional analog of the curl–operator given by curl z = curl(z1 , z2 ) = ∂x2 z1 − ∂x1 z2 . It is well known that the proper degrees of freedom for the space Zh is the mean value of the tangential components of v, v · t, with respect to each edge in ∆1 (Th ). Furthermore, we let Z 0 (3.3) Zh = {z ∈ Zh | z · t ds = 0, T ∈ Th∂ }. ∂T
Zh0
Alternatively, the elements of are those vector fields in Zh with the property that curl z|T = 0 if T ∈ Th∂ , i.e., z is a constant vector field on T for T ∈ Th∂ . It is a key observation that the mesh assumption given above, that any T ∈ Th intersects ∂Ω in at most
A UNIFORMLY STABLE FORTIN OPERATOR
7
one edge, implies that the spaces Vhb and Zh0 have the same dimension. Furthermore, we note that grad q ∈ Zh0 for any q ∈ Qh . For each e ∈ ∆1 (Th ) let φe ∈ Zh be the basis function corresponding to the Whitney form, i.e., φe satisfies Z Z (φe · te ) ds = |e|, and (φe · te0 ) ds = 0, e0 6= e, e0
e
where, as above, te is a tangent vector of length e. Hence, if e = (xi , xj ) then the vector field φe can be expressed in barycentric coordinates as φe = λi grad λj − λj grad λi . Any z ∈ Zh0 can be written uniquely on the form X X z= ae φe + e∈∆i1 (Th )
ce φe ,
e∈∆∂1 (Th )
where the coefficients ae corresponding to interior edges can be chosen arbitrarily, but where the coefficients ce for each boundary edge should be chosen such that curl z = 0 on the associated triangle in Th∂ . We note that there is a natural mapping Φh between the spaces Vhb and Zh0 given by Φh (ψe ) = φe for all interior edges, or alternatively, Z Z −2 v · te ds, e ∈ ∆i1 (Th ). Φh (v) · te ds = |e| (3.4) e
e
Vhb (T )
Zh0 (T )
Below, we will use and to denote the restriction of the spaces Vhb and Zh0 to a single element T , and ΦT will denote the corresponding restriction of the map Φh . We will define the operator Πhb : L2 (Ω; Rn ) → Vhb by, hΠhb u, zi = hu, zi
(3.5)
∀z ∈ Zh0 .
To show that this operator is well–defined, and behaves properly, the following general formula for integration of products of barycentric coordinates over a triangle T will be useful (cf. for example [14, Section 2.13]), Z 2α! (3.6) λα1 1 λα2 2 λα3 3 dx = |T |, (2 + |α|)! T P where α! = α1 !α2 !α3 ! and |α| = i αi . Lemma 3.1. Let T ∈ Thi with edges e1 , e2 , e3 . For any v ∈ Vhb (T ) we have Z 1 2 3 |a| |T | ≤ v · ΦT (v) dx ≤ |a|2 |T |, 5 5 T P P where ΦT = Φh |T , v = i ai ψei , and |a|2 = i a2i . Proof. A direct computation gives Z Z Z 3 X 3 3 X 3 X X v · ΦT (v) dx = ai aj ψei · φej dx = ai aj bei φej · tei dx = aT M a, T
i=1 j=1
T
i=1 j=1
T
8
¨ KENT–ANDRE MARDAL, JOACHIM SCHOBERL, AND RAGNAR WINTHER
where the 3 × 3 matrix M is given by Z Z 3 3 M = {Mi,j }i,j=1 = { ψi · φej dx}i,j=1 { bei φej · tei dx}3i,j=1 , T
T
The desired result will follow from the diagonal dominance of this matrix. Let xj be the vertex opposite to ej , and let λj be the corresponding barycentric coordinate on T . Then the first diagonal element M1,1 of the matrix M is given by Z Z M1,1 = 6 λ2 λ3 (λ2 grad λ3 − λ3 grad λ2 )(x3 − x2 ) dx = 6 (λ22 λ3 + λ2 λ23 ) dx = 2|T |/5, T
T
where we have used formula (3.6) in the final step. Actually, from this formula we derive that all the diagonal elements are given by Mi,i = 2|T |/5, and similar calculations for the off–diagonal elements gives |Mi,j | = |T |/10. In addition, the matrix M is symmetric. The matrix M is therefore strictly diagonally dominant, and by the Gershgorin circle theorem all eigenvalues are bounded below by |T |/5. By combining this with the fact that M is symmetric we conclude that aT M a ≥ |a|2 |T |/5, and this is the desired bound. The next lemma is a variant of the result above for T ∈ Th∂ . Lemma 3.2. Let T ∈ Th∂ with interior edges e1 , e2 , and where e3 is the edge on the boundary. For any v ∈ Vhb (T ) we have Z 1 v · ΦT (v) dx = |a|2 |T |, 2 T 2 2 2 where v = a1 ψe1 + a2 ψe2 and |a| = a1 + a2 . Proof. Let v = a1 ψe1 +a2 ψe2 = 6a1 λ2 λ3 (x3 −x2 )+6a2 λ1 λ3 (x3 −x1 ). It is a key observation that in this case ΦT (v) is simply given as Φ(v) = −a1 grad λ2 − a2 grad λ1 . Since grad λi · tei ≡ 0 we therefore obtain from (3.6) that Z Z Z 2 2 v · ΦT (v) dx = 6a1 λ2 λ3 grad λ2 · (x2 − x3 ) dx + 6a2 λ1 λ3 grad λ1 · (x1 − x3 ) dx T T ZT Z 1 = 6a21 λ2 λ3 dx + 6a22 λ1 λ3 dx = |a|2 |T |. 2 T T This completes the proof. P P P Let v = e∈∆i (Th ) ae ψe ∈ Vhb such that z = Φh (v) = e∈∆i (Th ) ae φe + e∈∆∂ (Th ) ce φe ∈ Zh0 . 1 1 1 It follows from scaling and shape regularity that the two norms X X (3.7) kzkL2 (Ω) and ( a2e )1/2 , T ∈Th e∈∆i1 (T )
are equivalent, uniformly with respect to h. Correspondingly, the two norms X X (3.8) kvkL2 (Ω) and ( |T |2 a2e )1/2 , T ∈Th
e∈∆i1 (T )
are uniformly equivalent. In fact, these equivalences hold locally for each T ∈ Th .
A UNIFORMLY STABLE FORTIN OPERATOR
9
It will be convenient to introduce the weighted L2 –norms given by X X 2 1/2 kvk1,h = ( h−2 and kzk−1,h = ( h2T kzk2L2 (T ) )1/2 . T kvkL2 (T ) ) T ∈Th
T ∈Th
In fact, for functions v ∈ Vhb the Pnorm kvk1,h is equivalent, uniformly in h, to the ordinary 1 H –norm. Furthermore, P P if v = e ae ψe then the two norms kvk1,h and kφh (v)k−1,h are both equivalent to ( T |T | e a2e )1/2 . This observation leads to the following inf–sup property. Lemma 3.3. There is a positive constant c0 , independent of h, such that for each v ∈ Vhb sup z∈Zh0
hv, zi ≥ c0 kvk1,h . kzk−1,h
P b Proof. Let v = e∈∆i1 (Th ) ae ψe ∈ Vh be given. We simply choose the corresponding z = Φh (v). As a consequence of the norm equivalences (3.7) and (3.8), combined with Lemmas 3.1 and 3.2, we obtain XZ X 1 X hv, Φh (v)i = Φh (v) · v dx ≥ |T | a2e ≥ c0 kΦh (v)k−1,h kvk1,h , 5 T i T ∈T T ∈T h
e∈∆1 (T )
h
where c0 > 0 is independent of h. This completes the proof.
The inf–sup condition just derived implies that the operator Πhb satisfies the stability estimate (3.9)
kΠhb vk1,h ≤ c−1 0 kvk1,h ,
v ∈ L2 (Ω).
If the family {Th } is quasi–uniform, i.e., there exists a positive constant γ1 , independent of h and T ∈ Th , such that hT ≥ γ1 h, then the estimate (3.9) also implies that Πhb is uniformly bounded in the standard non weighted L2 –norm. An alternative local mesh restriction, which implies L2 –stability, is the following condition. There exists a constant δ ∈ [0, 1), independent of h, such that |Ωe0 | ≤ (1 + δ)|Ωe |,
(3.10)
for all edges e, e0 ∈ ∆i1 (Th ) sharing a common vertex. The following consequence of this condition will be useful. Lemma 3.4. Assume that the local mesh condition (3.10) holds. There is a constant ∈ [0, 1/3), and for each h and T ∈ Th a positive constant αT , such that |1 −
|Ωe | | ≤ , αT
e ∈ ∆i1 (T ).
Proof. For each T ∈ Th let αT = (1 + δ/2) min |Ωe |. i e∈∆1 (T )
10
¨ KENT–ANDRE MARDAL, JOACHIM SCHOBERL, AND RAGNAR WINTHER
Due to condition (3.10) this implies that |Ωe | ≤ (1 + δ)(1 + δ/2)−1 αT for all e ∈ ∆i1 (T ). Therefore, we can conclude that 1 |Ωe | 1+δ ≤ ≤ , 1 + δ/2 αT 1 + δ/2 which is equivalent to the desired bound with = δ/(2 + δ) < 1/3.
Lemma 3.5. Assume that the local mesh condition (3.10) holds. There is a positive constant c1 , independent of h, such that for each v ∈ Vhb sup z∈Zh0
hv, zi ≥ c1 kvkL2 (Ω) . kzkL2 (Ω)
P Proof. We let and αT be as in the statement of Lemma 3.4. For a given v = e∈∆i (Th ) ae ψe ∈ 1 Vhb we choose z = Eh Φh (v), where Eh : Zh0 → Zh0 is defined by Eh φe = |Ωe |φe . As a consequence the L2 –norm of z is equivalent to the L2 –norm of v. It follows from these norm equivalences, combined with Lemmas 3.1 and 3.2, that XZ hv, zi = hv, Eh Φh (v)i = Eh Φh (v) · v dx T ∈Th
=
X T ∈Th
≥
X T ∈Th
T
Z
Z
(Eh − αT I)Φh (v) · v dx]
Φh (v) · v dx +
[αT
T
T
X 1 3 |Ωe | αT |T |[ − max (1 − )] a2e i 5 5 e∈∆1 (T ) αT i e∈∆1 (T )
X 1 3 X ≥( − ) αT |T | a2e ≥ c1 kvkL2 (Ω) kzkL2 (Ω) . 5 5 T ∈T i h
e∈∆1 (T )
This completes the proof.
It is a direct consequence of Lemma 3.5 that the operator norm of Πhb in L(L2 (Ω; R2 )) is bounded by c−1 1 . We now can conclude with the following main result of this section. Theorem 3.6. The operator Πh is uniformly bounded in H 1 and satisfies the commuting relation (2.7). Furthermore, if the triangulation is either quasi–uniform or it satisfies the local condition (3.10), then Πh is also uniformly bounded in L2 , i.e., (2.6) holds. Proof. To show that the operator Πh fulfills the condition (2.7), it is enough to show that the operator Πhb satisfies the corresponding condition (3.2). However, since grad Qh ⊂ Zh0 we have hdiv Πhb v, qi = −hΠhb v, grad qi = −hv, grad qi = hdiv v, qi.
A UNIFORMLY STABLE FORTIN OPERATOR
11
Furthermore, since Rh is bounded in H 1 , the desired H 1 boundedness of Πh will follow from the corresponding bound on Πhb (I − Rh ). However, by shape regularity, combined with (3.9) we have X X 2 kΠhb (I − Rh )vk2H 1 (Ω) ≤ c hT−2 kΠhb (I − Rh )vk2L2 (T ) ≤ c h−2 T k(I − Rh )vkL2 (T ) , T ∈Th
T ∈Th
and hence the H 1 –bound follows from (3.1). Finally, if the triangulation is quasi–uniform or satisfies (3.10) then the L2 –bound follows from the corresponding bounds for the operators Πhb and Rh . Remark 3.1. The two conditions guaranteeing L2 –stability of the Fortin operator Πh are rather different. Condition (3.10) is local, while quasi–uniformity is a global mesh condition. In some sense, this resembles two of the mesh conditions which can be used to obtain H 1 –stability of the L2 –projection onto continuous finite elements. It is a classical, and well–known, result that quasi–uniformity implies H 1 –stability of this operator, while an alternative local condition was derived in [3]. Remark 3.2. The analysis of the Taylor–Hood method given above leans heavily on the assumption that there are no triangles in Th with more than one edge on the boundary of Ω. This assumption simplifies the analysis, but it is not necessary. For more general triangulations, where some triangles may have two edges on the boundary, we will still define the interpolation operators Πh and Πhb as above. The only necessary modification is that the definition of the bubble space Vhb must be altered, since the dimension of the space Zh0 , given by (3.3), in general is larger than the number of interior edges. A possible choice is Vhb = span({ be te | e ∈ ∆i1 (Th ) } ∪ { be ne | e ∈ ∆i,2 1 (Th ) }), where te and ne are tangent and normal vectors to the edge e, and where ∆i,2 1 (Th ) is the i subset of ∆1 (Th ) consisting of the third edges of the triangles with two edges on the boundary. Hence, the space Vhb is enriched by including “the normal edge bubbles” associated the edges b 0 in ∆i,2 1 (Th ). With this definition the space Vh has the same dimension as Zh . In [15] the analysis above is extended to this more general case. However, the proof in this setting is more technical than the analysis above, and therefore it is omitted here. 4. Uniform preconditioners The interest in the construction of preconditioners for the discretizations of the singular perturbation system (1.1) is partly motivated by the role of such systems in time stepping schemes for the time dependent Stokes and Navier–Stokes equations. The system (1.1) has a coefficient operator which formally is of the form I − 2 ∆ − grad , A = div 0 with obvious weak interpretations of the differential operators derived from the weak formulation (2.1). In the notation of [16], the corresponding block diagonal operator B =
12
¨ KENT–ANDRE MARDAL, JOACHIM SCHOBERL, AND RAGNAR WINTHER
diag((I − 2 ∆)−1 , (−∆)−1 + 2 I), can be referred to as a canonical preconditioner for A , since the uniform inf–sup condition (2.2) implies that the composition B A is an isomorphism mapping the space (L2 ∩ H01 ) × ((H 1 ∩ L20 ) + −1 L20 ) to itself, with operator norms of B A and its inverse bounded independently of . We refer to [16, Section 4] and [17, Section 3] for more details. Furthermore, as explained in Section 2 above, if the Fortin operator Πh is uniformly bounded both in H 1 and L2 then the corresponding discrete inf–sup condition (2.5) holds for the Taylor–Hood method. Therefore, it follows from Theorem 3.6 that, if the triangulation {Th } is either quasi–uniform or satisfies the local condition (3.10), then (2.5) holds. In particular, the condition numbers of the corresponding discrete operators B,h A,h will be uniformly bounded with respect to and h, cf. [16, Section 5]. This result is obtained without any convexity assumption on the domain. Since all previous theoretical studies of uniform preconditioners for the model (1.1) have relied on convexity, cf. Remark 2.1 above, we will present a numerical experiment which confirm the theory above. More precisely, the experiment indicates that the condition number of B,h A,h is basically unaffected by such properties of the domain. The results below are obtained by using a canonical preconditioner, i.e., the operator B,h is composed of exact inverses of discretizations of the differential operators appearing in the definition of B . Of course, to obtain more practical preconditioners, these inverses should again be replaced by suitable elliptic preconditioners. We consider the problem (1.1) on three two dimensional domains, referred to as Ω1 , Ω2 and Ω3 . Here Ω1 is the unit square, Ω2 is the L–shaped domain obtained by cutting out the an upper right subsquare from Ω1 , while Ω3 is the slit domain where a slit of length a half is removed from from Ω1 , cf. Figure 1. Hence, only Ω1 is convex. The corresponding problems
Figure 1. The domains Ω1 , Ω2 , and Ω3 . (1.1) were discretized by the lowest order Taylor–Hood element on a uniform triangular grid. The condition numbers of the corresponding discrete operators B,h A,h , for the three domains and for different values of and h, are computed. The results are given in Table 1 below. These results clearly indicate that the condition numbers of the operators B,h A,h are not dramatically effected by lack of convexity of the domains. References [1] D.N. Arnold, F. Brezzi, and M. Fortin, A stable finite element method for Stokes equations, Calcolo 21 (1984), 337–344. [2] J.H. Bramble and J.E. Pasciak, Iterative techniques for time dependent Stokes problems, Comput. Math. Appl. 33 (1997), no. 1-2, 13–30, Approximation theory and applications. [3] J.H. Bramble, J.E. Pasciak, and O. Steinbach On the stability of the L2 projection in H 1 (Ω), Math. Comp. 71 (2001), 147–156.
A UNIFORMLY STABLE FORTIN OPERATOR
domain
Ω1
Ω2
Ω3
\h
2−2
2−3
2−4
2−5
1
14.3 15.3 18.5 22.0
0.1
10.3 11.9 12.7 13.2
0.01
6.2
6.7
8.0
9.9
1
17.1 17.2 17.1 17.1
0.1
10.3 11.8 12.7 13.2
0.01
6.1
6.7
8.0
13
9.9
1
13.2 13.4 13.5 13.6
0.1
10.3 11.9 12.8 13.2
0.01 6.1 6.7 8.1 9.9 Table 1. Condition numbers for the operators B,h A,h
[4] F. Brezzi, On the existence, uniqueness and approximation of saddle–point problems arising from Lagrangian multipliers, RAIRO. Analyse Num´erique 8 (1974), 129–151. [5] F. Brezzi and M. Fortin, Mixed and hybrid finite element methods, Springer-Verlag, 1991. [6] J. Cahouet and J.-P. Chabard, Some fast 3D finite element solvers for the generalized Stokes problem, Internat. J. Numer. Methods Fluids 8 (1988), no. 8, 869–895. [7] P. Clement, Approximation by finite element functions using local regularization, RAIRO Anal. Num´er. 9 (1975), 77–84. [8] M. Costabel and A. McIntosh, On Bogovski˘i and regularized Poincar´e integral operators for de Rham complexes on Lipschitz domains, http://arxiv.org/abs/0808.2614v1. [9] R.S Falk, A Fortin operator for two–dimensional Taylor–Hood elements, Mathematical Modelling and Numerical Analysis (M2AN) 42 (2008), 411-424. [10] Giovanni P. Galdi, An introduction to the mathematical theory of the Navier–Stokes equations, vol 1, Linearized Steady problems, vol. 38, Springer Tracts in Natural Philosophy, Springer–Verlag, New York 1994. [11] M. Geissert, H. Heck, and M. Hieber, On the equation div u = g and Bogovski˘i’s operator in Sobolev spaces of negative order, Oper. Theory Adv. Appl. 168 (2006), 113–121. [12] V. Girault and P.-A. Raviart, Finite element methods for Navier-Stokes equations, Theory and algorithms, Springer Series in Computational Mathematics, vol. 5, Springer-Verlag, Berlin, 1986. [13] T. Kolev and P. Vassilevski, Regular decomposition of H(div) spaces, Lawrence Livermore National Laboratory, 2011. [14] M.-J. Lay and L.L. Schumaker, Spline functions on triangulations, vol. 110 of Encyclopedia of Mathematics and its Applications, Cambridge University Press, 2007. [15] K.-A. Mardal, J. Sch¨ oberl, and R. Winther, A uniform inf–sup condition with applications to preconditioning, arXiv:1201.1513v1,2012. [16] K.-A. Mardal and R. Winther, Preconditioning discretizations of systems of partial differential equations, Numer. Linear Alg. Appl. 18 (2011), 1–40. [17] K.-A. Mardal and R. Winther, Uniform preconditioners for the time dependent Stokes problem, Numer. Math. 98 (2004), no. 2, 305–327. , Erratum: “Uniform preconditioners for the time dependent Stokes problem” [Numer. Math. 98 [18] (2004), no. 2, 305–327 ], Numer. Math. 103 (2006), no. 1, 171–172.
14
¨ KENT–ANDRE MARDAL, JOACHIM SCHOBERL, AND RAGNAR WINTHER
[19] M. A. Olshanskii, J. Peters, and A. Reusken, Uniform preconditioners for a parameter dependent saddle point problem with application to generalized Stokes interface equations, Numerische Mathematik 105 (2006), 159–191. [20] S. Turek, Efficient solvers for incompressible flow problems, Springer-Verlag, 1999.
Centre for Biomedical Computing at Simula Research Laboratory, Norway and Department of Informatics, University of Oslo, Norway E-mail address: [email protected] URL: http://simula.no/people/kent-and Institute for Analysis and Scientific Computing, Wiedner Hauptstrasse 8-10, 1040 Wien, Austria E-mail address: [email protected] URL: http://www.asc.tuwien.ac.at/~schoeberl Centre of Mathematics for Applications and Department of Informatics, University of Oslo, 0316 Oslo, Norway E-mail address: [email protected] URL: http://folk.uio.no/rwinther
