In this paper we study the discretization of the Stokes problem (2.1)–(2.3) by the finite element pair P 1-bubble/ P 1 (the so-called mini-element), introduced by Arnold et al. (1984). This element leads to a relatively low number of degrees of freedom with a good approximate solution.

Let T h be a triangulation of Ω and T a triangle of T h . We define the space associated with the bubble by B h = { v h ∈ C 0 ( Ω ¯ ) ; ∀ T ∈ T h , v h | T = x b ( T ) } . We also defined the discrete function spaces V i h = { v h ∈ C 0 ( Ω ¯ ) ; v h | T ∈ P 1 , ∀ T ∈ T h ; v h | Γ = 0 } , i = 1 , … , d , P h = { q h ∈ C 0 ( Ω ¯ ) ; q h | T ∈ P 1 , ∀ T ∈ T h ; : ∫ Ω q h d x = 0 } , and we set X i h = V i h ⊕ B h and X h = X 1 h × ⋯ × X d h . With the above preparations, the discrete variational problem reads as follows.

Find ( u h , p h ) ∈ X h × P h such that (3.1) a ( u h , v h ) − ( p h , ∇ · v h ) Ω = ( f , v h ) Ω ∀ v h ∈ X h , (3.2) − ( q h , ∇ · u h ) Ω = 0 ∀ q h ∈ P h .

For a given triangle T, the velocity field u h and the pressure p h are approximated by linear combinations of the basis functions in the form u h ( x ) = ∑ i = 1 d + 1 ϕ i ( x ) u i + u b ϕ b ( x ) , p h ( x ) = ∑ i = 1 d + 1 ϕ i ( x ) p i , where u i and p i are nodal values of u h and p h while u b is the bubble value. The basis functions are defined by ϕ 1 ( x ) = 1 − x − y , ϕ 2 ( x ) = x , ϕ 3 ( x ) = y , ϕ b ( x ) = 27 ∏ i = 1 3 ϕ i ( x ) in 2D, and ϕ 1 ( x ) = 1 − x − y − z , ϕ 2 ( x ) = x , ϕ 3 ( x ) = y , ϕ 4 ( x ) = z , ϕ b ( x ) = 256 ∏ i = 1 4 ϕ i ( x ) in 3D.

To construct the coefficient matrices, a number of integrals involving powers of the basis functions will be computed. Integrals over a triangle T can be evaluated directly by the following formulas (3.3) ∫ T ϕ 1 α 1 ϕ 2 α 2 ϕ 3 α 3 d x = 2 | T | α 1 ! α 2 ! α 3 ! ( α 1 + α 2 + α 3 + 2 ) ! , (3.4) ∫ T ϕ 1 α 1 ϕ 2 α 2 ϕ 3 α 3 ϕ 4 α 4 d x = 6 | T | α 1 ! α 2 ! α 3 ! α 4 ( α 1 + α 2 + α 3 + α 4 + 3 ) ! , where | T | stands for the triangle area (in 2D), or the tetrahedron volume (in 3D). A useful property for the basis functions is (3.5) ∑ i = 1 d + 1 ∇ ϕ i = 0 .

We use the following notations for the discrete velocity/pressure nodal values (3.6) u ¯ i = u i u i b , f ¯ i = f i f i b , i = 1 , … , d . Systems (3.1)–(3.2) lead to the following algebraic form, using notations (3.6) (3.7) A ¯ 0 − B ¯ 1 ⊤ 0 A ¯ − B ¯ 2 ⊤ − B ¯ 1 − B ¯ 2 0 u ¯ 1 u ¯ 2 p = f ¯ 1 f ¯ 2 0 in 2D, or (3.8) A ¯ 0 0 − B ¯ 1 ⊤ 0 A ¯ 0 − B ¯ 2 ⊤ 0 0 A ¯ − B ¯ 3 ⊤ − B ¯ 1 − B ¯ 2 − B ¯ 3 0 u ¯ 1 u ¯ 2 u ¯ 3 p = f ¯ 1 f ¯ 2 f ¯ 3 0 in 3D. In (3.7) and (3.8) A ¯ = M ¯ + R ¯, with M ¯ as the mass matrix and R ¯ as the stiffness matrix. B ¯ i is the divergence submatrix associated with the i-th partial derivative, i.e. B ¯ i ≡ ( q h , ∂ i u i h ) Ω , i = 1 , … , d . To create the algebraic system (3.7) or (3.8), the discrete system (3.1)–(3.2) is evaluated over each triangle T to obtain the element matrices and vectors M ¯ i j ( T ) = ∫ T α ϕ i ϕ j d x , R ¯ i j ( T ) = ∫ T ν ∇ ϕ i · ∇ ϕ j d x , B ¯ i j ( T ) = ∫ T ∂ 1 ϕ i ϕ j d x + ∫ T ∂ 2 ϕ i ϕ j d x , f ¯ i ( T ) = ∫ T f ϕ i d x . Then assembling operations consist of direct-summing the element matrices over the triangulation T h to obtain the global matrices M ¯ = ( M ¯ i j ), R ¯ = ( R ¯ i j ), B ¯ = ( B ¯ i j ) and f ¯ = ( f ¯ i ) M ¯ i j = ∑ T ∈ T h M ¯ i j ( T ) , R ¯ i j = ∑ T ∈ T h R ¯ i j ( T ) , B ¯ i j = ∑ T ∈ T h B ¯ i j ( T ) , f ¯ i = ∑ T ∈ T h f ¯ i ( T ) . In the next sections we detail the element matrices M ¯, R ¯, B ¯ i and the right-hand side f ¯.

For a triangle T, let { ( x i , y i ) } i = 1 , 2 , 3 be the vertices and { ϕ } i = 1 , 2 , 3 the corresponding basis functions. The gradient of ϕ i is given by (4.1) ∇ ϕ 1 t ∇ ϕ 2 t ∇ ϕ 3 t = 1 1 1 x 1 x 2 x 3 y 1 y 2 y 3 − 1 0 0 1 0 0 1 = 1 2 | T | y 2 − y 3 x 3 − x 2 y 3 − y 1 x 1 − x 3 y 1 − y 2 x 2 − x 1 , where | T | is the area of T given by 2 | T | = det x 2 − x 1 x 3 − x 1 y 2 − y 1 y 3 − y 1 = ( x 2 − x 1 ) ( y 3 − y 1 ) − ( x 3 − x 1 ) ( y 2 − y 1 ) .

A nonbubble entry of the element stiffness matrix R ¯ is given by (4.2) R ¯ i j = ∫ T ν ∇ ϕ i · ∇ ϕ j d x = | T | ν ∇ ϕ i · ∇ ϕ j , i , j = 1 , 2 , 3 . For the bubble entries R ¯ b j , for j = 1 , 2 , 3. A straightforward calculation yields to (using (3.5)) R ¯ b j = 9 4 | T | ∑ i = 1 3 ∇ ϕ i = 0 , j = 1 , 2 , 3 . For the diagonal entry corresponding to the bubble (i.e. i = j = b) we have (4.3) R ¯ b b = ν ∫ T 27 2 ∇ ( ϕ 1 ϕ 2 ϕ 3 ) · ∇ ( ϕ 1 ϕ 2 ϕ 3 ) d x = 81 10 ν | T | | ∇ ϕ 1 | 2 + | ∇ ϕ 2 | 2 + | ∇ ϕ 3 | 2 + ∇ ϕ 1 ∇ ϕ 2 + ∇ ϕ 1 ∇ ϕ 3 + ∇ ϕ 2 ∇ ϕ 3 = 81 10 ν | T | | ∇ ϕ 1 | 2 + | ∇ ϕ 2 | 2 + ∇ ϕ 1 · ∇ ϕ 2 = : ω R , using (3.5).

With the above results, the element stiffness matrix is therefore R ¯ = R 0 0 ω R .

As for the stiffness matrix, we set M = ( M ¯ i j ) i , j = 1 , … , 3 as the nonbubble part of the mass matrix. A direct calculation yields (4.4) M i j = α 6 | T | if i = j , α 12 | T | elsewhere . The bubble part of the mass matrix is given by (4.5) M ¯ b j = 3 α 20 | T | , j = 1 , 2 , 3 , M ¯ b b = 81 α 280 | T | = : ω M . The element mass matrix is therefore M ¯ = M z z ⊤ ω M , where z ⊤ = 3 20 α | T | 1 1 1 .

Finally, the element stiffness/mass matrix A ¯ is A ¯ = A z z ⊤ ω , where we have set A = R + M and ω = ω R + ω M .

A direct integration yields the element divergence matrix − B ¯ = [ − B ¯ 1 − B ¯ 2 ], where B ¯ i = [ B i − B i b ] with (4.6) B i = | T | 3 ∂ i ϕ 1 ∂ i ϕ 2 ∂ i ϕ 3 ∂ i ϕ 1 ∂ i ϕ 2 ∂ i ϕ 3 ∂ i ϕ 1 ∂ i ϕ 2 ∂ i ϕ 3 , i = 1 , 2 , and (4.7) B i b = 9 | T | 20 ∂ i ϕ 1 ∂ i ϕ 2 ∂ i ϕ 3 , i = 1 , 2 .

The contribution of the right-hand side component f i , in nonbubble terms, is given by (4.8) f i ( T ) = | T | 3 f i T , i = 1 , 2 , 3 , where f i T is a mean value of f i on T. The bubble terms of the right-hand side are (4.9) f i b ( T ) = ∫ T f i ϕ b d x = 9 20 | T | f i T , i = 1 , 2 .

With the above element matrices and vectors the 11 × 11 element system corresponding to (3.7) is A z 0 0 − B 1 ⊤ z ⊤ ω 0 0 B 1 b ⊤ 0 0 A z − B 2 ⊤ 0 0 z ⊤ ω B 2 b ⊤ − B 1 B 1 b − B 2 B 2 b 0 u 1 u 1 b u 2 u 2 b p = f 1 f 1 b f 2 f 2 b 0 . To reveal diagonal blocks, the system above can be rearranged as follows (4.10) A 0 z 0 − B 1 ⊤ 0 A 0 z − B 2 ⊤ z ⊤ 0 ω 0 B 1 b ⊤ 0 z ⊤ 0 ω B 2 b ⊤ − B 1 − B 2 B 1 b B 2 b 0 u 1 u 2 u 1 b u 2 b p = f 1 f 2 f 1 b f 2 b 0 . We can now eliminate the bubble unknowns u 1 b and u 2 b since they correspond to diagonal blocks (the ω blocks). From (4.10) 3 and (4.10) 4 , we deduce that (4.11) u i b = ( f i b − B i b ⊤ p − z ⊤ u i ) / ω , i = 1 , 2 . Substituting (4.11) into (4.10) 1 , (4.10) 2 and (4.10) 5 we obtain a linear system in ( u 1 u 2 p ) ⊤ whose matrix is (4.12) A − ω − 1 z z ⊤ 0 − B 1 ⊤ − ω − 1 z B 1 b ⊤ 0 A − ω − 1 z z ⊤ − B 2 ⊤ − ω − 1 z B 2 b ⊤ − B 1 − ω − 1 B 1 b z ⊤ − B 2 − ω − 1 B 2 b z ⊤ − ω − 1 ( B 1 b B 1 b ⊤ + B 2 b B 2 b ⊤ ) and the right-hand side (4.13) f 1 − ω − 1 z f 1 b f 2 − ω − 1 z f 2 b − ω − 1 ( B 1 b f 1 b + B 2 b f 2 b ) .

Let { x i = ( x i , y i , z i ) } i = 1 , … , 4 be the vertices of a tetrahedron T and { ϕ i } i = 1 , … , 4 the corresponding basis functions. The gradient ∇ x ϕ i over T are given by ∇ ϕ 1 ⊤ ∇ ϕ 2 ⊤ ∇ ϕ 3 ⊤ ∇ ϕ 4 ⊤ = 1 1 1 1 x 1 x 2 x 3 x 4 y 1 y 2 y 3 y 4 z 1 z 2 z 3 z 4 − 1 0 0 0 1 0 0 0 1 0 0 0 1 . An alternative formula for computing ∇ ϕ i is ∇ x ϕ i = J − 1 ∇ ξ ϕ i ( ξ ) , where J is the Jacobean matrix of the mapping ξ ↦ x ( ξ ), that is, (4.14) J = x 2 − x 1 y 2 − y 1 z 2 − z 1 x 3 − x 1 y 3 − y 1 z 3 − z 1 x 4 − x 1 y 4 − y 1 z 4 − z 1 . The volume of tetrahedron T is given, from (4.14), by | T | = det ( J ) / 6.

As for 2D case, the nonbubble entries of the element stiffness matrix are given by (4.15) R ¯ b j = 9 4 | T | ∑ i = 1 3 ∇ ϕ i = 0 , j = 1 , 2 , 3 , 4 , while R ¯ b j = 0, for all j = 1 , … , 4. For the diagonal entry, using (3.4)–(3.5), we obtain (4.16) R ¯ b b = ν ∫ T 256 2 ∇ ( ϕ 1 ϕ 2 ϕ 3 ϕ 4 ) · ∇ ( ϕ 1 ϕ 2 ϕ 3 ϕ 4 ) d x = 8192 845 ν | T | ∑ ℓ = 1 3 | ∇ ϕ ℓ | 2 + ∇ ϕ 1 · ∇ ϕ 2 + ∇ ϕ 1 · ∇ ϕ 3 + ∇ ϕ 2 · ∇ ϕ 3 = : ω R .

For the 3D mass matrix, a straightforward calculation with a linear tetrahedron shows that (for nonbubble entries) (4.17) M ¯ i j = α T | T | 10 if i = j , α T | T | 20 if i ≠ j , where α T is a mean value of α on T. The bubble part of the mass matrix is given by (4.18) M ¯ b j = 8 α 105 | T | , j = 1 , 2 , 3 , M ¯ b b = 8192 α 51975 | T | = : ω M . The element mass matrix is therefore M ¯ = M z z ⊤ ω M , where z ⊤ = 8 α 105 | T | 1 1 1 1 .

If we set − B ¯ = [ − B ¯ 1 − B ¯ 2 − B ¯ 3 ] and B ¯ i = [ B i − B i b ], a straightforward calculation yields to (4.19) B i = | T | 4 ∂ i ϕ 1 ∂ i ϕ 2 ∂ i ϕ 3 ∂ i ϕ 4 ∂ i ϕ 1 ∂ i ϕ 2 ∂ i ϕ 3 ∂ i ϕ 4 ∂ i ϕ 1 ∂ i ϕ 2 ∂ i ϕ 3 ∂ i ϕ 4 ∂ i ϕ 1 ∂ i ϕ 2 ∂ i ϕ 3 ∂ i ϕ 4 , i = 1 , … , 4 , and (4.20) B i b = 32 105 | T | ∂ i ϕ 1 , ∂ i ϕ 2 ∂ i ϕ 3 ∂ i ϕ 4 .

The contribution of the right-hand side component f i , in nonbubble terms, is given by (4.21) f i ( T ) = | T | 4 f i T , i = 1 , 2 , 3 , 4 , where f i T is a mean value of f i on T. The bubble terms of the right-hand side are (4.22) f i b ( T ) = ∫ T f i ϕ b d x = 32 105 | T | f i T , i = 1 , 2 , 3 .

As in 2D, we rearrange the Stokes system and we get (4.23) u i b = ( f i b − B i b ⊤ p − z ⊤ u i ) / ω , i = 1 , … , 4 . Substituting (4.23) into the rest of the system, we obtain a linear system whose matrix is (4.24) A − ω − 1 z z ⊤ 0 0 − B 1 ⊤ − ω − 1 z B 1 b ⊤ 0 A − ω − 1 z z ⊤ 0 − B 2 ⊤ − ω − 1 z B 2 b ⊤ 0 0 A − ω − 1 z z ⊤ − B 3 ⊤ − ω − 1 z B 3 b ⊤ − B 1 − ω − 1 B 1 b z ⊤ − B 2 − ω − 1 B 2 b z ⊤ − B 3 − ω − 1 B 3 b z ⊤ − ω − 1 ∑ i = 1 3 B i b B i b ⊤ and the right-hand side (4.25) f 1 − ω − 1 z f 1 b f 2 − ω − 1 z f 2 b − ω − 1 ∑ i = 1 3 B i b B i b ⊤ .

To derive a dual (Uzawa) algorithm for (5.1), we assume that u = u ( p ) is the solution of Poisson equation (5.4) A u = f + B ⊤ p , that is, in the decomposed useful form A i u i = f i + B i ⊤ p , i = 1 , … , d . Then multiplying (5.4) by u and substituting the result in (5.2), we obtain (5.5) L ( u ( p ) , p ) = − 1 2 u ⊤ A u − 1 2 p ⊤ C p − f p ⊤ p . If we introduce the dual functional J ∗ ( p ) = − min u L ( u ( p ) , p ) with u ( p ) solution of (5.4), the saddle-point problem (5.3) becomes.

Find p such that (5.6) J ∗ ( p ) ⩽ J ∗ ( q ) , ∀ q . From (5.5), J ∗ is quadratic and coercive. From (5.4), we deduce that the mapping p ↦ u ( p ) is linear and u ( p + t d ) = u ( p ) + t w where w is the solution of the sensitivity problem (5.7) A w = B ⊤ d , or A i w i = B i ⊤ d , i = 1 , … , d . It follows that the derivative of J ∗ is (5.8) r : = ∇ J ∗ ( p ) = B u + C p + f p . With a search direction d, we compute an optimal stepsize ρ by solving ∇ J ∗ ( p + ρ d ) ⊤ d = 0 , that is, ρ = − d ⊤ ( B w + C d ) / ( r ⊤ d ) , where w is the solution of the sensitivity equation (5.7).

Since J ∗ is a quadratic function, the best search direction is a conjugate gradient direction. As a consequence, the best algorithm for (5.6) is a conjugate gradient algorithm. At each iteration k, the Fletcher–Reeves conjugate gradient direction is given by d k = r k + 1 + β k d k , β k = r k + 1 2 r k − 2 .

A dual (Uzawa) conjugate gradient algorithm for solving the saddle-point problem (5.1) is outlined in Algorithm 1. Theoretically, Algorithm 1 converges in at most n B = rank ( B ) iterations. Obviously, for large scale problems, n B is very large and it is preferable that convergence should be obtained in a number of iterations considerably less than n B .

A practical implementation of a conjugate gradient method for solving (5.1) requires a preconditioner, that is, a suitable scalar product for computing ∇ J ∗ ( p ) instead of the standard one used in (5.8). The convergence properties of the conjugate gradient method for the generalized Stokes problem are deteriorated for large values of the ratio α / ν, see e.g. (Glowinski, 2003; Glowinski and Le Tallec, 1989).

To derive a preconditioner for (5.1) following the idea of Cahouet and Chabard (1988), Glowinski and Guidoboni (2009), we need to simplify the continuous problem. We first notice that the equivalence between P 1-Bubble/ P 1 element and the stabilized formulation has been proved (Pierre, 1987, 1989; 1995; Matsumoto, 2005; Baiocchi et al., 1993). Then, if we neglect the bubble contribution in the stiffness and divergence matrices, (5.1) can be expressed in the strong form as (5.9) α u − ν Δ u + ∇ p = f , (5.10) ∇ · u − ∇ · ( ν h ∇ u ) = 0 , where ν h is an element dependent constant.

As in Glowinski and Guidoboni (2009) we define the linear operator from L 2 ( Ω ) into L 2 ( Ω ) (neglecting the constant term in (5.8)) (5.11) ϕ : = A q = ∇ · u q − ∇ · ( ν h ∇ q ) , where u q is the solution of (5.12) α u q − ν Δ u q = − ∇ q . Note that (5.12) is the (strong) sensitivity system. The idea behind preconditioning is to find a linear operator B such that B A q = q. Applying the divergence operator in (5.12), we obtain α ∇ · u q − ν ∇ · Δ u q = − Δ q or (5.13) α ∇ · u q − ν Δ ( ∇ · u q ) = − Δ q . Using (5.11) in (5.13), we obtain (5.14) − Δ q − α ∇ · ( ν h ∇ q ) + ν Δ ( ∇ · ( ν h ∇ q ) ) = α ϕ − ν Δ ϕ . Then, in practice, at each step of the preconditioned conjugate gradient algorithm, the gradient of J ∗ is computed as an approximate solution of the linear system (5.15) ( K + α C + ν K M − 1 C ) g = ( α M + ν K ) r , where K and M are the stiffness and the mass matrices, respectively. Let us introduce the mesh Reynolds number R e h = α ν h 2 , where h is the mesh size. Taking into account the CPU time and the storage requirement for computing the last term of the matrix involved in (5.15), we consider the following preconditioning system instead (5.16) ( K + α C ) g = H r , where (5.17) H = ν K + α diag ( M ) if R e h ⩽ 1 , ν K + α I if R e h > 1 . 1.\end{array}\right.\]]]> In (5.16) and (5.17), K and M are P 1 stiffness and mass matrix, respectively, and C, the bubble matrix from Functions 6.1–6.3.