Distributed Lagrange multipliers based on fictitious domain method for second order elliptic problems

Comput. Methods Appl. Mech. Engrg. 196 (2007) 1498–1506 www.elsevier.com/locate/cma

Distributed Lagrange multipliers based on fictitious domain method for second order elliptic problems R. Glowinski *, Yu. Kuznetsov Department of Mathematics, University of Houston, 651 PGH, Houston, TX 77204-3008, USA Received 24 August 2005; received in revised form 21 April 2006; accepted 24 May 2006

Abstract In this article we further investigate the solution of linear second order elliptic boundary value problems by distributed Lagrange multipliers based fictitious domain methods. The following issues are addressed: (i) Derivation of the fictitious domain formulations. (ii) Finite element approximation. (iii) Iterative solution of the resulting finite dimensional problems (of the saddle-point type) by preconditioned conjugate gradient and Lanczos algorithms.  2006 Elsevier B.V. All rights reserved. Keywords: Second order elliptic problems; Fictitious domain methods; Lagrange multipliers; Finite element approximations; Preconditioned iterative methods; Conjugate gradient algorithms; Lanczos algorithms

1. Introduction To the best of our knowledge, fictitious domain methods (also known as domain embedding methods) have been introduced by Hyman (Ref. [18]) and Saul’ev (Refs. [27,28]) about 50 years ago (indeed the terminology is from Saul’ev) and have since then been applied to the solution of a variety of problems of increasing complexity (for a – necessary biased – historical account of fictitious domain methods see, for example, [2,3,5–8,9,Section 38,17,20,22, 24], and the references therein). Let us briefly recall the basic principle of fictitious domain methods: in order to solve a partial differential equation (or inequality) problem (P) ‘‘taking place’’ in a spatial domain X of complicated shape, we embed X in a domain G having a simpler shape and we extend (P) by ð Pe Þ such that when one restricts to X the solution ~ u of problem ð Pe Þ, it coincides (exactly or approximately) with the solution of problem (P); the main idea behind this approach is that for a domain of a simple shape one can expect (or develop) fast specialized solvers

*

Corresponding author. Tel.: +1 713 743 34 73; fax: +1 713 743 30 05. E-mail address: [email protected] (R. Glowinski).

0045-7825/$ - see front matter  2006 Elsevier B.V. All rights reserved. doi:10.1016/j.cma.2006.05.013

for a given class of problems. There are many ways to implement the above simple idea and the related literature has become sufficiently large that we have no choice but to limit the references to those works which have some relation to the present article (asking the forgiveness of those authors whose contributions are not mentioned). Some years ago (1996 to be precise) fictitious domain method with distributed Lagrange multipliers was introduced as a tool for the direct numerical simulation of particulate flow; combined with time discretizations by operator-splitting; these fictitious domain methods ‘‘opened the door’’ to the simulation of three-dimensional flow involving thousands of particles. For those situations where the boundary of the flow region varies with time the fictitious domain methodology was allowing the flow computations to be done on a fixed space region, namely, the union of the region occupied by the fluid with the region occupied by the particles (see, e.g., [12–14,9] [Chapters 8 and 9] for details). More recently, a fictitious domain method with distributed Lagrange multipliers has been (successfully) applied to the solution of the wave equation modeling propagation phenomena in heterogeneous media [4]. Despite the above successes, many fundamental issues have not been investigated in depth, concerning in particular the solution of

R. Glowinski, Yu. Kuznetsov / Comput. Methods Appl. Mech. Engrg. 196 (2007) 1498–1506

problems as basic as the linear second elliptic ones. This is a review article written for the special issue of CMAME. Our main goal is to give an extended description of the distributed Lagrange multipliers based fictitious domain methods – briefly discussed in [15,16] – with a particular attention being given to the iterative solution, by preconditioned conjugate gradient or Lanczos algorithms, of the finite dimensional saddle-point type problems resulting from the space approximation, by finite element methods, for example. As it was mentioned above the main goal of the fictitious domain method is to replace the system of mesh equations in a complex shaped domain by a system of mesh equations in a simpler shaped domain. The reason for such a replacement is to get further options for designing more efficient preconditioned algebraic solvers even if the new system has a larger dimension. Of course, the extended system should be either equivalent or ‘‘approximately’’ equivalent to the original system, i.e. we should get the solution of the original system by solving the new one. The fictitious domain method described in this paper is used to replace a finite element system for second order elliptic PDE in a complex shaped domain, with Dirichlet or mixed type boundary conditions, by a system for an elliptic equation with natural (Neumann or Robin) boundary conditions in a simpler shaped domain. The reason for the replacement of Dirichlet boundary condition by natural boundary conditions is that the variant of the fictitious domain method invented by Astrakhantsev in [1,2], was designed especially for elliptic problems with natural boundary conditions, and it works very well in numerous applications. To this end, the Astrakhantsev’s fictitious domain method can be used for the preconditioning of the matrices which appear in the Lagrange multipliers based fictitious domain method. The article is organized as follows. In Section 2, we describe the basic idea of the method on a hierarchy of elliptic problems with increasing complexity. In Section 3, we apply the method to the finite element discretizations of the second order elliptic equations and discuss the algebraic aspects of the method. In Section 4, we briefly recall the solution of saddle-point linear systems by preconditioned iterative methods. Finally, in Section 5, we describe a particular approach for designing of an efficient preconditioner for the underlying matrices and provide a condition number estimate.

1499

We assume that f 2 L2(X) and consider the following variational formulation of (1): Find u 2 H 10 ðXÞ, such that UX ðuÞ 6 UX ðvÞ 8v 2 H 10 ðXÞ; where UX ðvÞ ¼

1 2

Z

jrvj2 dx 

ð2Þ

Z

ð3Þ

fv dx:

X

X

The solution to problem (2) and (3) exists and is unique. 2.1. Fictitious domain method for a simply connected domain Let X be a simply connected domain, and G be a bounded domain in Rd , d = 2, 3, containing X, i.e. X  G. We shall denote G n X by x and call x a fictitious domain. We assume that x consists of m simply connected regularshaped subdomains xk, k = 1, . . . , m, which do not touch k \ x  l ¼ ;, l 5 k, k,l = 1, . . . , m. An exameach other, i.e. x ple of X and G is shown in Fig. 1. A fictitious domain formulation of problem (2) as a constrained minimization problem can be defined as follows: Find w 2 H 10 ðGÞ, such that UG ðwÞ 6 UG ðvÞ where UG ðvÞ ¼ and KF ¼



1 2

Z

8v 2 K F ;

jrvj2 dx 

ð4Þ

Z

ð5Þ

fv dx X

G

v : v 2 H 10 ðGÞ;

 2 jrvj dx ¼ 0 :

Z

ð6Þ

x

Problem (4)–(6) has unique solution w 2 H 10 ðGÞ. It is obvious that w = 0 in x and w coincides with the solution u to (2) and (3) in X. An equivalent distributed Lagrange multipliers based formulation of problem (4)–(6) is defined by Find w 2 H 10 ðGÞ, k 2 H1(x), k = 0 on Cx \ CG, such that Z Z Z rw  rv dx þ rk  rv dx ¼ fv dx; x X ZG ð7Þ rw  rl dx ¼ 0; x

2. Dirichlet problem In this section, we consider the fictitious domain method in the differential variational setting. Let X be a bounded domain in Rd , d = 2, 3, with a piecewise smooth boundary C  CX satisfying the cone condition. The model problem is defined as follows: Du ¼f

ω1 Ω ω2

in X

u ¼0 on C:

ð1Þ

Fig. 1. The domain G consists of the subdomains X, x1, and x2.

1500

R. Glowinski, Yu. Kuznetsov / Comput. Methods Appl. Mech. Engrg. 196 (2007) 1498–1506

for all v 2 H 10 ðGÞ, l 2 H1(x), l = 0 on Cx \ CG, where Cx and CG are the boundaries of x and G, respectively. The function k in (7) is said to be a distributed Lagrange multiplier. One can easily verify that the function k in (7) is a harmonic function, i.e. Dk ¼ 0 in x;

ð8Þ

and satisfies the boundary conditions k ¼ 0 on Cx \ CG ; o ðk  uÞ ¼ 0 on Cx \ C; on

ð9Þ

where u is the solution function to (1), and n is the outward unit normal to C.

In this subsection, we consider the case when X is not a simply connected domain. Let G be a simply connected bounded domain in Rd , d = 2, 3, with boundary CG, and xk, k = 1, . . . , m, be intek \ x  l ¼ ;, l 5 k, and rior subdomains of G, such that x CG \ Ck = ;, k, l = 1, . . . , m. Here Ck denotes the boundary of xk, k = 1, . . . , m. The domain X for problem (1) is defined as follows: ! m [ k : X¼Gn ð10Þ x k¼1

An example of G with subdomains X, x1, and x2 is shown in Fig. 2. A fictitious domain formulation of (2) and (3) as a constrained minimization problem can be defined by: Find w 2 H 10 ðGÞ, such that where 1 UG ðvÞ ¼ 2

"Z

8v 2 K F ;

2

jrvj dx þ G

xk

ð13Þ Here, ak is a nonnegative constant and bk is a positive constant, k = 1, . . . , m. The choice of the coefficients ak and bk is important when the diameters of the subdomains xk (holes) are much smaller than the diameter of G. Otherwise, we can always choose ak = 0 and bk = 1, k = 1, . . . , m. We shall discuss another choice of these coefficients when considering efficient preconditioners for the problem with small size holes. An equivalent distributed Lagrange multipliers based formulation of problem (11)–(13) is defined as follows: Find w 2 H 10 ðGÞ, kk 2 H1(xk), k = 1, . . . , m, such that Z

2.2. Domains with holes

UG ðwÞ 6 UG ðvÞ

  Z 2 K F ¼ v : v 2 H 10 ðGÞ; ½jrvj þ bk v2 dx ¼ 0; k ¼ 1; . . . ; m :

ð11Þ

m X k¼1

Z ak

# 2

jrvj dx 

xk

rw  rv dx þ G

þ Z

k¼1

m Z X k¼1

m X

Z wv dx

ak xk

½rkk  rv þ bk kk v dx ¼

xk

½rw  rlk þ bk wlk  dx ¼ 0;

Z fv dx; X

k ¼ 1; . . . ; m;

ð14Þ

xk

for all v 2 H 10 ðGÞ, lk 2 H1(xk), k = 1, . . . , m. Problems (11)–(14) are completely equivalent and have a unique solution w 2 H 10 ðGÞ, kk 2 H1(xk), k = 1, . . . , m, such that w  u in X and w  0 in xk, k = 1, . . . , m. One can easily verify that the Lagrange multipliers kk are the solutions of the differential problems  Dkk þ bk kk ¼ 0 in xk ; ð15Þ o ðkk  uÞ ¼ 0 on Ck ; on where u is the solution function in X and n is the normal unit vector on Ck directed from X into xk, k = 1, . . . , m. 2.3. Mixed boundary conditions

Z fv dx

We consider the two-dimensional Poisson equation

X

ð12Þ

Du ¼ f

in X;

ð16Þ

with the boundary conditions

and

u ¼ 0 on CD ; ou ¼ 0 on CN ; on

ω1

ω2 Ω

Fig. 2. Domain G with subdomains X, x1, and x2.

ð17Þ

where CD is a closed subset of C (the Dirichlet part of C) and CN is an open subset of C (the Neumann part of C), such that CD \ CN = ; and CD [ CN ¼ C. Let G be a bounded domain in R2 , such that X  G and x ¼ G n X is a union of m nonintersecting subdomains xk, k \ x  l ¼ ;, l 5 k, k, l = 1, . . . , m. We k = 1, . . . , m, i.e. x e k ¼ Ck \ C between assume that the interface boundaries C xk and X belong to CD, k = 1, . . . , m. Thus, each of the boundaries Ck of the subdomains xk consists of the interb k , k = 1, . . . , m. face with X and the non-interface part C

R. Glowinski, Yu. Kuznetsov / Comput. Methods Appl. Mech. Engrg. 196 (2007) 1498–1506

Γ1, D

for all v 2 H1(G), v = 0 on CG,D, and lk 2 H1(xk), lk = 0 on Ck,D, k = 1, . . . , m. Again, problems (18)–(21) are completely equivalent and have a unique solution w 2 H1(G), w = 0 on CG,D, kk 2 H1(xk), kk = 0 on Ck,D, k = 1, . . . , m, such that w = u in X and w = 0 in xk, k = 1, . . . , m. One can easily verify that kk is the solution of the differential problem

ΓN

Γ1, N

ΓD

ω1 ~ Γ1

ΓD

~ Γ2 Γ2, D

ω2

Ω

b k which is assigned for the DirichLet Ck,D be a part of C let boundary condition in the further fictitious domain formulation for (16) and (17). An example of the subdomains and its boundaries is given in Fig. 3. We denote by CG,D the union of CG \ CD and Ck,D, k = 1, . . . , m, where CG is the boundary of G. Then, the fictitious domain formulation in terms of the constrained minimization problem is defined by: Find w 2 H1(G), w = 0 on CG,D, such that where 1 UG ðvÞ ¼ 2

"Z

8v 2 K F ;

2

jrvj dx þ G

m X

ð18Þ Z ak

k¼1

# 2

jrvj dx 

xk

Z fv dx X

ð19Þ and

1

v : v 2 H ðGÞ; v ¼ 0 on CG;D ;  Z 2 ½jrvj þ bk v2  dx ¼ 0; k ¼ 1; . . . ; m :

ð20Þ

xk

Here, ak and bk are nonnegative constants, k = 1, . . . , m. In the case Ck,D = ; the constant bk should be positive. In the case CG,D = ; at least one of the constants ak, k = 1, . . . , m, should be positive. The equivalent distributed Lagrange multipliers based formulation of (18)–(20) is defined as follows: Find w 2 H1(G), w = 0 on CG,D, and kk 2 H1(xk), kk = 0 on Ck,D, k = 1, . . . , m, such that Z Z m X rw  rv dx þ ak wv dx G

þ

m X k¼1

Z

We consider the differential problem 

d X o ou aij þ cu ¼ f oxi oxj i;j¼1

u¼0

in X;

on CD ;

ðaruÞ  n þ ru ¼ gR

ð23Þ

on CR ;

where a = (aij) is a symmetric matrix uniformly positive definite in X with piecewise smooth bounded from above entries aij, i,j = 1, . . . , d, c and r are bounded piecewise smooth nonnegative functions, n is the outward unit normal to C, and f 2 L2(X). The sets CD and CR are the Dirichlet and Robin/Neumann parts of the boundary, respectively, such that CD \ CR = ;, CD is a closed set, and CD [ CR ¼ C. The weak formulation of (23) is defined as follows: Find u 2 V, such that aðu; vÞ ¼ lðvÞ 8v 2 V ;

ð24Þ

where



KF ¼

ð22Þ

2.4. Generalization

Fig. 3. Partitionings of the boundaries for X, x1, and x2.

UG ðwÞ 6 UG ðvÞ

 Dkk þ bk kk ¼ 0 in xk ; o e k; ðkk  uÞ ¼ 0 on C on

where n is the outward unit normal to C, k = 1, . . . , m.

Γ2, N

ΓN

1501

Z

k¼1

xk

½rkk  rv þ bk kk v dx ¼

xk

Z fv dx; X

½rw  rlk þ bk wlk  dx ¼ 0;

xk

k ¼ 1; . . . ; m;

ð21Þ

V ¼ fv : v 2 H 1 ðXÞ; v ¼ 0 on CD g; Z Z d Z X ou ov aðu; vÞ ¼ aij dx þ cuv dx þ ruv ds; oxi oxj X CR i;j¼1 X Z Z lðvÞ ¼ fv dx þ gR v ds: X

CR

ð25Þ The bilinear form a(u, v) is symmetric. We assume that a(u, v) is V-elliptic. Let G be a bounded domain in Rd , d = 2, 3, such that X  G and x ¼ G n X is a union of m nonintersecting k \ domains xk with boundaries Ck, k = 1, . . . , m, i.e. x  l ¼ ;, l 5 k, k,l = 1, . . . , m. We shall refer to x as a fictix tious domain with respect to X. We assume that the intere k ¼ Ck \ C between xk and X belong to face boundaries C CD, k = 1, . . . , m. Thus, each of the boundaries Ck consists e k and the non-interface part C b k , k = 1, . . . , m. For each of C b k we define subsets C b k;D and denote by Ck,D the closure of C b k;D , k = 1, . . . , m. We shall impose the Dirichlet boundary C condition on Ck,D in the fictitious domain method to be proposed for problem (24) and (25). Of course, Ck,D is

1502

R. Glowinski, Yu. Kuznetsov / Comput. Methods Appl. Mech. Engrg. 196 (2007) 1498–1506

the empty set if we do plan not impose the Dirichlet boundb k in the fictitious domain method. ary condition on C To formulate the fictitious domain method we define the bilinear forms ak ðu; vÞ ¼

d Z X i;j¼1

bk ðk; vÞ ¼

xk

d Z X i;j¼1

ðkÞ aij

ðkÞ

bij

xk

ou ov dx þ oxi oxj ok ov dx þ oxi oxj

Z ck uv dx; ð26Þ

Z

k¼1

ð27Þ  ð28Þ

k¼1

where k 2 Vk and on CG;D g:

ð29Þ

Here, CG;D ¼ ðCD \ CG Þ [

m [

! ð30Þ

Ck;D :

k¼1 ðkÞ

ðkÞ

We assume that the matrices aðkÞ ¼ ðaij Þ and bðkÞ ¼ ðbij Þ are symmetric and uniformly positive definite in xk with piecewise smooth bounded from above entries, and ck and bk are nonnegative piecewise smooth bounded from above functions in xk, k = 1, . . . , m. We also assume that the bilinear forms bk in (26) are positive definite in Vk, k = 1, . . . , m, and the bilinear form aF in (28) is positive definite in VF. The fictitious domain formulation to (24) is defined as follows: Find w 2 VF, kk 2 Vk, k = 1, . . . , m, such that aF ðw; vÞ þ bðk; vÞ ¼ lðvÞ; bðw; lÞ ¼ 0

ð31Þ

for all v 2 VF, lk 2 Vk, k = 1, . . . , m. Problems (24) and (31) are completely equivalent. The solution function w of (31) coincides with the solution function u of (24) in X, and w = 0 in xk, k = 1, . . . , m. The Lagrange multiplier functions kk in (31) are the solutions of the following problem: Find kk 2 Vk, such that bk ðkk ; vÞ ¼ gk ðvÞ where gk ðvÞ ¼

Z eC k

8v 2 V k ;

in X;

ð34Þ

on C;

where C is the boundary of X. In [27,28], Saul’ev suggested to replace problem (34) by the problem

bk ðk; vÞ;

u; v 2 V F ¼ fv : v 2 H 1 ðGÞ; v ¼ 0

d X o ou aij ¼f ox ox i j i;j¼1

u¼0

bk kv dx; xk

u; v; k 2 V k ¼ fv : v 2 H 1 ðxk Þ; v ¼ 0 on Ck;D g; m X aF ðu; vÞ ¼ aðu; vÞ þ ak ðu; vÞ; m X

For the sake of simplicity, we assume that the coefficient c in (23) is equal to zero, and CD = C, i.e. we consider the differential problem 

xk

where

bðk; vÞ ¼

2.5. Connection to the method of Saul’ev

ð32Þ

d X o ou aij; ¼f oxi oxj i;j¼1

in G;

ð35Þ

u ¼ 0 on CG ; where G is a bounded domain in Rd , and X  G. We shall impose on G the same constrains as in Section 2.4. The coefficients in (35) are defined by 8 in X; < aij aij; ¼ 1 ðkÞ ð36Þ : aij in xk k ¼ 1; . . . ; m;  where   1 is a ‘‘small’’ positive parameter. In [25], Rivkind gave the estimation of the error with respect to . The proposed approach was named by Saul’ev as the fictitious domain method. The weak formulation of problem (35) and (36) is as follows: Find u 2 H 10 ðGÞ, such that d Z d Z X ou ov 1X ou ov ~aij ~aij dx þ dx ox d ox oxi oxj i j i;j¼1 G i;j¼1 x Z ¼ fv dx 8v 2 H 10 ðGÞ; ð37Þ X

where 1d ¼ 1  1 and ( aij in X; ~aij ¼ ðkÞ aij in xk k ¼ 1; . . . ; m:

ð38Þ

Then an equivalent weak formulation of (35) and (36) in terms of distributed Lagrange multipliers is given by: Find u 2 H 10 ðGÞ, k 2 H1(x), k = 0 on CG, such that Z d Z d Z X X ou ov ok ov ~aij ~aij dx þ dx ¼ fv dx; oxi oxj oxi oxj X i;j¼1 G i;j¼1 x Z d Z X ou ol ~aij dx  d k l dx ¼ 0; oxi oxj x i;j¼1 x ð39Þ

½ðaruÞ  nv ds;

ð33Þ

u is the solution function of the original problem (24), and n is the outward unit normal to C, i = 1, . . . , m.

for all v 2 H 10 ðGÞ and l 2 H1(x), l = 0 on CG. Here,   1 1  1 u : k ¼ u  ð40Þ d 

R. Glowinski, Yu. Kuznetsov / Comput. Methods Appl. Mech. Engrg. 196 (2007) 1498–1506

1503

It is obvious that in the case d = 0 problem (39) coincides with the fictitious domain method (31). Using the terminology of S. Sobolev, the distributed Lagrange multipliers based fictitious domain method (31) is the closure of the fictitious domain method by Saul’ev when  ! 0.

of P1 finite element functions defined on xh,k and vanishing on Ch,k,D, k = 1, . . . , m. The finite element counterpart of the fictitious domain method (31) for problem (41) is defined as follows: Find wh 2 Vh,F, kh,k 2 Vh,k, k = 1, . . . , m, such that

3. Fictitious domain method for finite element problems

bh ðwh ; lh Þ ¼ 0

Let Xh be a simplicial (triangular, d = 2, or tetrahedral, d = 3) mesh with the boundary Ch partitioned into subsets Ch,D and Ch,R. We assume that Ch, Ch,D, and Ch,R approximate properly the boundaries C, CD, and CR in problem (23), respectively. We also assume that the interfaces between subdomains in X with jumps of the tensor a = a(x) are approximated properly by faces and/or edges of the mesh cells in Xh. Let Vh be the space of P1 finite element functions vanishing on Ch,D. Then, the P1 finite element approximation of (24) is defined as follows: Find uh 2 Vh, such that

for all vh 2 Vh,F, lh,k 2 Vh,k, k = 1, . . . , m. Finite element problem (45) results in the system of linear algebraic equations      w g A  ¼ ; ð46Þ k 0

ah ðuh ; vÞ ¼ lh ðvÞ

8v 2 V h ;

ð41Þ

where ah ðu; vÞ ¼

d Z X i;j¼1

ah;ij

Xh

ou ov dx þ oxi oxj

Z

ch uv dx þ Xh

Z rh uv ds Ch;R

ð42Þ and lh ðvÞ ¼

Z Xh

fh v dx þ

Z gh;R v ds:

ð43Þ

Ch;R

Here, ah,ij, ch, rh, fh, and gh,R are properly defined prolongations of aij, c, r, f, and gR, respectively. The latter assumption is typical in finite element methods for problems with curved boundaries and interfaces. In the case of polygonal/polyhedral domains and boundaries, we usually do not need the prolongations. The finite element problem (41) leads to the following system of the linear algebraic equations K u ¼ f ;

ð44Þ

with the symmetric positive definite matrix K and the right hand side vector f . To introduce the distributed Lagrange multipliers based fictitious domain method for problem (41) we have to perform similar to the above substitutions for the domains G, e k, C b k;R , and bilinear forms aF, b, ak, xk, boundaries Ck, C and bk, k = 1, . . . , m. We assume that all the substitutions are done properly. We shall indicate these substitutions by additional index h in the notations. For instance, Xh and xh,k, k = 1, . . . , m, are proper mesh subdomains of a simplicial mesh Gh, and Ch,G,D is a proper approximation to CG,D. Let Vh,F be the space of P1 finite element functions defined on Gh and vanishing on Ch,G,D, and Vh,k be the trace of Vh,F onto xh,k, i.e. Vh,k be the space

ah;F ðwh ; vh Þ þ bh ðkh ; vh Þ ¼ lh ðvh Þ;

with a symmetric n · n matrix ! A BT ; A¼ B 0

ð45Þ

ð47Þ

where A is a symmetric positive definite nF · nF matrix and B is a rectangular nk · nF matrix. To describe the block structure of matrices A and B, and the values of nF and nk we introduce the following set of matrices. We define the nk · nk matrices Ak and Bk by  k ; vk Þ ¼ ah;k ðwh;k ; vh;k Þ 8wh;k ; vh;k 2 V h;k ; ðAk w k Þ ¼ bh;k ðwh;k ; lh;k Þ 8wh;k ; lh;k 2 V h;k ; k; l ðBk w

ð48Þ

 k , vk , l k 2 Rnk , and nk is equal to the dimension of where w Vh,k, k = 1, . . . , m. It is obvious that the matrices Ak and Bk are symmetric and positive definite "k = 1, . . . , m. We define the n0 · n0 matrix A0 by  0 ; v0 Þ ¼ aðwh;0 ; vh;0 Þ ðA0 w

8wh;0 ; vh;0 2 V h;0 ;

ð49Þ

 0 , v0 2 Rn0 , Vh,0 is the restriction of Vh onto Xh, and where w n0 is the dimension of Vh,0. The matrix A is defined by the assembling of the matrices A0, A1, . . . , Am with an appropriate choice of assembling matrices Pk , k = 0,1, . . . , m, i.e. A¼

m X

Pk Ak PTk :

ð50Þ

k¼0

We assume that the matrices Pk , k = 1, . . . , m, are chosen such that 0 1 0 0 0 B0 A m 0 C X 1 B C C; ð51Þ Pk Ak PTk ¼ B .. B C @ A . k¼1 0 0 Am is an (m + 1) · (m + 1) block diagonal matrix with blocks A1, . . . , Am on the diagonal. The first diagonal block in (51) is the square null matrix. Under the above assumptions the matrix B is an (m + 1) · m block matrix. Namely

1504

R. Glowinski, Yu. Kuznetsov / Comput. Methods Appl. Mech. Engrg. 196 (2007) 1498–1506

0

0

B1

B0 B B¼B B .. @. 0

0 .. . 0

0 B2 .. . 0

...

0

... 0 . .. . .. . . . Bm

1 C C C: C A

v0 2 Rn ; ð52Þ

ð54Þ

where AFF and BF are symmetric positive definite nk · nk matrices. It is obvious that the submatrix A00 is equal to the matrix K in (44) and the subvector g0 is equal to the vector f in (44).  F ¼ 0. From the last matrix equation in (53) we get w  0 satisfies Then, system (53) implies that the subvector w the system  0 ¼ g0 ; A00 w

ð55Þ

which coincides with system (44):  0 ¼ f ; Kw

ð56Þ

and the system for the subvector  k:  0: k ¼ gF  AF 0 w BF 

ð59Þ

where e1 = 0,

The matrices A and B in (47) as well as system (46) can be also represented in another convenient block form. We  and g in (46) into partition the components of the vectors w two groups. In one group we include the components which correspond to the nodes of the mesh XhnCh,D, and into the other group the rest of the components. Then system (46) can be written in the following equivalent block form 0 10 1 0 1 0 A00 A0F 0 g0 w B CB C B C  F A ¼ @ gF A: ð53Þ @ AF 0 AFF BF A@ w  k 0 BF 0 0 Here,     A00 A0F 0 A¼ ; BT ¼ ; BF AF 0 AFF     0 g0 w ¼ ; g ¼ ; w F w gF

1 ½HAHnk  ek1 ðvk  vk1 Þ; qk

vkþ1 ¼ vk 

ð57Þ

0 ¼  u, and system (44) is equivalent to the fictitious Thus, w domain system (46) with respect to the solution subvector  0. w 4. Preconditioned iterative solvers 4.1. Preconditioned conjugate gradient method

qk ¼

kAHnk k2H  ek1 ; 2 knk k

e k ¼ qk

H

knkþ1 k2H ; 2 knk k

ð60Þ

H

k = 0, 1, . . . Here, nk ¼ Avk  g is the residual vector and  2 ¼ ðHw;  wÞ,  w  2 Rn . kwk H 1 Let v ¼ A g be the solution vector to (58) and zk ¼ vk  v , k P 0. Then, the convergence estimate kzk kH1 6

Ck



1 c2 þc1 c2 c1

 kz0 kH1 ;

ð61Þ

holds "k P 1. Here, C k ðxÞ ¼ coshðkcosh1 xÞ;

ð62Þ

are the Chebyshev polynomials, c1 is the minimal eigen2 value of the matrix ðHAÞ , and c2 is its maximal eigenvalue. In the case, when c1 and c2 are positive and 2 estimate the eigenvalues of ðHAÞ from below and from above, respectively, estimate (61) holds as well. 4.2. Preconditioned Lanczos method Let H 2 Rnn be again a symmetric positive definite matrix. Then, the preconditioned Lanczos method [19,21], is another good candidate for the iterative solution of system (58). The preconditioned Lanczos method reads as follows: v0 2 Rn ; 8 k ¼ 1; Hn0 ; < k p ¼ k ¼ 2; Hp1  a2 p1 ; : k1 k1 k2 Hp  ak p  bk p ; k P 3; vk ¼ vk1  ck pk ;

ð63Þ ð64Þ

where ak ¼

ðAHpk1 ; Apk1 Þ kApk1 k2H

;

2

bk ¼

kApk2 kH

ð65Þ

; kApk1 k2H ðnk1 ; Apk Þ ck ¼ ; 2 kApk kH k P 1, and nk ¼ Avk  g, k P 0.

nn

Let H 2 R be a symmetric positive definite matrix. Then, the minimal errors based preconditioned conjugate gradient method [21] is a good candidate for the efficient iterative solution of a system Av ¼ g;

ð58Þ

with a symmetric nonsingular matrix A 2 Rnn and a given vector g 2 Rn . This method reads as follows:

For the preconditioned Lanczos method the convergence estimate kn2k kH 6

Ck



1

c2 þc1 c2 c1

 kn0 kH ;

ð66Þ

holds "k P 1, with the same values of c2 and c1 as used in estimate (61).

R. Glowinski, Yu. Kuznetsov / Comput. Methods Appl. Mech. Engrg. 196 (2007) 1498–1506

1505

5. Preconditioning

S 1 F AF is bounded from below by a positive constant 1/(1 + c3). Here

5.1. Block diagonal preconditioner

S F ¼ AFF  AF 0 A1 00 A0F ;

ð70Þ

In paper [15], the authors proposed the preconditioner H to be a 2 · 2 block diagonal matrix, i.e.   H 0 H¼ ; ð67Þ 0 Hk

is the Schur complement matrix, and m X AF ¼ Pk Ak PTk :

ð71Þ

with symmetric positive definite diagonal submatrices H 2 RnF nF and H k ¼ Bk 2 Rnk nk . If we assume that the eigenvalues of the matrices HA and Hk(BA1BT) belong to a positive segment [a; b] then the eigenvalues of the matrix ðHAÞ2 belong to a positive segment [c1; c2] with the endpoints pffiffiffi pffiffiffi 2 2 2 ð 5 þ 1Þ 2 ð 5  1Þ and c2 ¼ b : ð68Þ c1 ¼ a 4 4 2

Better estimates for the eigenvalues of ðHAÞ can be found in [26]. In this paper, we do not discuss preconditioners H for the matrix A in (47). In fact, the major goal of the fictitious domain methods is to provide a better framework for the construction of an efficient preconditioner H. In the case of rectangular domains G and rectangular meshes Gh a preconditioner H can be a separable matrix [3,7,10,11]. In another case, when the mesh Gh is fully or logically hierarchical one, the multigrid and multilevel methods provide good tools for designing efficient preconditioners H for the matrix A with the spectral boundaries a and b independent of the mesh. Such preconditioners H1 are said to be spectrally equivalent to A. It is also an option to use Astrakhantsev’s fictitious domain method [1,2] for elliptic problems with natural (Neumann or Robin) boundary conditions to design matrix H in (67) if Gh is still a complex shaped domain. To this end, we concentrate ourselves on the analysis of the eigenvalues of the matrix Hk(BA1BT) with Hk = Bk. 5.2. Eigenvalue estimates The eigenvalue estimates for the matrix Hk(BA1BT) are based on the finite element norm preserving extension theorem [1,23,29]. In our particular case the latter result means that for any finite element functions vh,k 2 Vh,k, k = 1, . . . , m, there exists a finite element function vh,0 2 Vh,0, vh,0 = vh,k on oXh \ oxh,k, i.e. vh,0 = vh,k on the interface boundaries between Xh and xh,k, k = 1, . . . , m, such that the inequality ah;0 ðvh;0 ; vh;0 Þ 6 c3

m X

ah;k ðvh;k ; vh;k Þ;

ð69Þ

k¼1

We observe that the maximal eigenvalue of the matrix S 1 F AF is equal to one. Another statement to be used in the analysis of eigenvalues is that the matrices Ak and Bk are spectrally equivalent due to the equivalence of the bilinear forms ak and bk on Vk, k = 1, . . . , m. It means that positive constants c4 and c5 exist such that c4 ðAk v; vÞ 6 ðBk v; vÞ 6 c5 ðAk v; vÞ 8v 2 Rnk ;

ð72Þ

k = 1, . . . , n. The eigenvalue problem H k ðBA1 BT Þz ¼ mz;

ð73Þ

is equivalent (due to the block partitionings of A and B in (54)) to the eigenvalue problem z ¼ mz: S 1 F BF 

ð74Þ

On the basis of the above statement we conclude that the eigenvalues in (74) are bounded from below by the positive constant c4 aF ¼ ð75Þ 1 þ c3 and the maximal eigenvalue in (74) is bounded from above by the positive constant bF = c5. Both constants aF and bF are independent of the meshes Xh and xh,k, k = 1, . . . , m. 6. Conclusion This article is a rather simple introduction to distributed Lagrange multipliers based fictitious domain methods and preconditioning of the underlying algebraic systems. The point is that the constants in the condition number estimates for the preconditioned matrix HA may strongly depend on the shape of the original domain X and the extended domain G. For instance, the condition number strongly depends on the shapes of the holes as well as the distribution of the holes in X. The same problem arises for elliptic equations with strong discontinuities in the coefficients. The problems may be overcome by either designing more efficient preconditioners for the fictitious domain systems or using more sophisticated versions of the distributed Lagrange multipliers technique. The authors are in favor of the latter approach. References

k¼1

holds with a positive constant c3 independent of vh,k, k = 1, . . . , m. The straightforward consequence of the latter statement is that the minimal eigenvalue of the matrix

[1] G. Astrakhantsev, Iterative methods for solving variational difference schemes for two dimensional second order elliptic equations, Doctoral Thesis. LOMI Akad. Nauk SSSR, Leningrad, 1972 (in Russian).

1506

R. Glowinski, Yu. Kuznetsov / Comput. Methods Appl. Mech. Engrg. 196 (2007) 1498–1506

[2] G. Astrakhantsev, Method of fictitious domains for a second-order elliptic equation with natural boundary conditions, USSR Comput. Math. Math. Phys. 18 (1978) 114–121. [3] A. Bespalov, Yu. Kuznetsov, O. Pironneau, M.-G. Vallet, Fictitious domains with separable preconditioners versus unstructured adapted meshes, IMPACT Comput. Sci. Engrg. 4 (1992) 217–249. [4] V.A. Bokil, R. Glowinski, An operator splitting scheme with a distributed Lagrange multiplier based fictitious domain method for wave propagation problems, J. Comput. Phys. 205 (2005) 242–268. [5] B.L. Buzbee, G.H. Golub, C.W. Nielson, On direct methods for solving Poisson’s equations, SIAM J. Numer. Anal. 7 (1970) 627–656. [6] Q.V. Dinh, R. Glowinski, J. He, T.W. Pan, J. Pe´riaux, Lagrange multiplier approach to fictitious domain methods: application to fluid dynamic and electromagnetics, in: D.E. Keyes, T.F. Chan, G.A. Meurant, J.S. Scroggs, R.G. Voigt, (Eds.), Proceedings of the Fifth International Symposium on Domain Decomposition Methods for PDEs, (Norfolk, VA, May 1991), SIAM, Philadelphia (1992) pp. 151–194. [7] S. Finogenov, Yu. Kuznetsov, Two-stage fictitious components method for solving the Dirichlet boundary value problem, Soviet J. Numer. Anal. Math. Model. 3 (1988) 301–323. [8] V. Girault, R. Glowinski, Error analysis of a fictitious domain method applied to a Dirichlet problem, Jpn. J. Ind. Appl. Math. 12 (1995) 487–514. [9] R. Glowinski, Numerical methods for fluids, in: P.G. Ciarlet, J.-L. Lions (Eds.), Handbook of Numerical Analysis, vol. IX, Elsevier, Amsterdam, 2003. [10] R. Glowinski, T.I. Hesla, D.D. Joseph, T.W. Pan, J. Pe´riaux, Distributed Lagrange multiplier methods for particulate flow, in: M.O. Bristeau, G. Etgen, W. Fitzgibbon, J.-L. Lions, J. Pe´riaux, M. Wheeler (Eds.), Computational Sciences for the 21st Century, J. Wiley, Chichester, 1997, pp. 270–279. [11] R. Glowinski, T.I. Hesla, D.D. Joseph, T.W. Pan, J. Pe´riaux, A fictitious domain method with distributed Lagrange multipliers for the numerical simulation of particulate flow, in: J. Mandel, C. Farhat, X.-C. Cai (Eds.), Domain Decomposition Methods, 10, Contemporary Mathematics, 218, AMS, Providence, 1998, pp. 121–137. [12] R. Glowinski, T.I. Hesla, D.D. Joseph, T.W. Pan, J. Pe´riaux, A distributed Lagrange multiplier/fictitious domain method for flows around moving rigid bodies: application to particulate flow, Int. J. Numer. Methods Fluids 30 (1999) 1043–1066. [13] R. Glowinski, T.I. Hesla, D.D. Joseph, T.W. Pan, J. Pe´riaux, A distributed Lagrange multiplier/fictitious domain method for the simulation of flow around moving rigid bodies: Application to particulate flow, Comput. Methods Appl. Mech. Engrg. 184 (2000) 241–267. [14] R. Glowinski, T.I. Hesla, D.D. Joseph, T.W. Pan, J. Pe´riaux, A fictitious domain approach to the direct numerical simulation of incompressible viscous fluid flow past moving rigid bodies: application to particulate flow, J. Comput. Phys. 169 (2001) 363–426. [15] R. Glowinski, Yu. Kuznetsov, On the solution of the Dirichlet problem for linear elliptic operators by a distributed Lagrange multiplier method, C. R. Acad. Sci. Paris, t. 327, Se´rie I, (1998) 693– 698. [16] R. Glowinski, Yu. Kuznetsov, T. Rossi, J. Toivanen, A fictitious domain method with distributed Lagrange multipliers, in: P. Nei-

[17]

[18] [19]

[20]

[21]

[22]

[23]

[24]

[25]

[26] [27]

[28]

[29]

ttaanma¨ki, T. Tiihonen, P. Tarvainen (Eds.), ENUMATH 99, Proceedings Third European Conf. on Numer. Math. and Advanced Appl., (University of Jyva¨skyla¨, Jyva¨skyla, Finland, July 1999), World Scientific, 2000, pp. 733–742. E. Heikkola, Yu. Kuznetsov, T. Rossi, P. Tarvainen, Efficient preconditioners based on fictitious domains for elliptic FE-problems with Lagrange multipliers, Report 11, University of Jyva¨skyla¨, Dept. of Math., (1996), in: H.G. Bock, F. Brezzi, R. Glowinski, G. Kanschat, Yu. Kuznetsov, J. Pe´riaux, R. Rannacher, (Eds.), ENUMATH 97, Proceedings Second European Conf. on Numer. Math. and Advanced Appl., (Heidelberg, Germany, September–October 1997), World Scientific (1998), pp. 646–661. M.A. Hyman, Non-iterative numerical solution of boundary value problems, Appl. Sci. Res. Sec. B2 (1952) 325–351. Yu. Kuznetsov, Efficient iterative solvers for elliptic finite element problems on non-matching grids, Russian J. Numer. Anal. Math. Model. 10 (1995) 187–211. Yu. Kuznetsov, A. Matsokin, On the optimization of fictitious component method, in: G.I. Marchuk, (Ed.), Vychislitel’nye metody algebry, (Vychisl. Tsentr Sib. Otdel. Akad. Nauk SSSR, Novosibirsk, 1977), pp. 79–86 (in Russian). G.I. Marchuk, Yu. Kuznetsov, Iterative methods and quadratic functionals, Nauka, Novosibirsk, 1972, (in Russian). Translation: Me´thodes ite´ratives et fonctionnelles quadratiques. In: Sur les me´thodes nume´riques en sciences physique et e´conomique. (Paris, Dunod, 1974), pp. 3–132. G.I. Marchuk, Yu. Kuznetsov, A. Matsokin, Fictitious domain and domain decomposition methods, Soviet J. Numer. Anal. Math. Model. 1 (1986) 3–35. S. Nepomnyaschikh, Decomposition and fictitious domain methods for elliptic boundary value problems, in: D.E. Keyes, T.F. Chan, G.A. Meurant, J.S. Scroggs, R.G. Voigt (Eds.), Proceedings Fifth International Symposium on Domain Decomposition Methods for PDEs, (Norfolk, VA, May 1991), SIAM, Philadelphia, PA, 1992, pp. 62–72. W. Proskurowski, O. Widlund, On the numerical solution of Helmholtz’s equation by the capacitance matrix method, Math. Comput. 30 (1976) 434–468. V. Rivkind, Estimates of the rate of convergence of solutions of difference equations to solutions of elliptic equations with discontinuous coefficients and a numerical method of solving the Dirichlet problem, Dokl. Akad. Nauk SSSR 149 (1963), 1264–1268 (in Russian). English translation in J. Soviet Math., 4 (1964) pp. 546–550. T. Rusten, R. Winther, A preconditioned iterative method for saddle point problems, SIAM J. Matrix Anal. Appl. 13 (1992) 1352–1367. V. Saul’ev, A method for automatization of the solution of boundary value problems on high performance computers, Dokl. Akad. Nauk SSSR 144 (1962), 497–500 (in Russian). English translation in Soviet Math. Dokl. 3 (1963), pp. 763–766. V. Saul’ev, On solution of some boundary value problems on high performance computers by fictitious domain method, Siberian Math. J. 4 (1963) 912–925, in Russian. O. Widlund, An extension theorem for finite element spaces with three applications, in: W. Hackbusch, K. Witsch, (Eds.), Notes on Numerical Fluid Mechanics 16, (Friedr. Vieweg und Sohn, 1987), pp. 110–122.