Optimal preconditioning for the symmetric and non-symmetric coupling of adaptive finite elements and boundary elements

Optimal preconditioning for the symmetric and non-symmetric coupling of adaptive finite elements and boundary elements


We analyze a multilevel diagonal additive Schwarz preconditioner for the adaptive coupling of FEM and BEM for a linear 2D Laplace transmission problem. We rigorously prove that the condition number of the preconditioned system stays uniformly bounded, independently of the refinement level and the local mesh-size of the underlying adaptively refined triangulations. Although the focus is on the non-symmetric Johnson-Nédélec one-equation coupling, the principle ideas also apply to other formulations like the symmetric FEM-BEM coupling. Numerical experiments underline our theoretical findings.

Key words and phrases:
FEM-BEM coupling, preconditioner, multilevel additive Schwarz, adaptivity
2010 Mathematics Subject Classification:
65N30, 65N38, 65F08
Acknowledgment: The research of the authors MF, TF, and DP is supported through the FWF research project Adaptive Boundary Element Method, see http://www.asc.tuwien.ac.at/abem/, funded by the Austrian Science Fund (FWF) under grant P21732, as well as through the Innovative Projects Initiative of Vienna University of Technology. Moreover, MF and DP are partially funded through the FWF doctoral school Dissipation and Dispersion in Nonlinear PDEs, funded under grant W1245. This support is thankfully acknowledged.

1. Introduction

There exist plenty of works on preconditioning of FEM-BEM coupling equations, covering mainly the symmetric coupling with quasi-uniform meshes, see [CKL98, FS09, HPPS03, HMS99, HS98, KS02, MS98] and the references therein. In contrast to that, only little is known on preconditioning of the non-symmetric Johnson-Nédélec coupling, see e.g. [Med98], and also on preconditioning of adaptive FEM-BEM couplings. It is the main goal of this paper to close this gap and to extend the existing analysis to the case of the (adaptive) symmetric as well as non-symmetric Johnson-Nédélec coupling [JN80]. For the symmetric coupling [Cos88, Han90], the approach of Bramble & Pasciak [BP88] applies, which guarantees positive definiteness and symmetry of the Galerkin matrix with respect to a special inner product [CKL98, HPPS03, KS02]. Therefore, efficient iterative solvers designed for symmetric and positive definite matrices are applicable. However, due to the non-symmetry, such an approach may not work for the Johnson-Nédélec coupling in general. In [Med98], it was assumed that the coupling boundary is smooth. Hence, the double-layer integral operator is compact. The system matrix can therefore be split into a symmetric part plus a compact perturbation part (the Galerkin matrix of the double-layer integral operator). Preconditioning is done only on the symmetric part with the theory of [BP88], and convergence results for iterative solvers can then be obtained by compact perturbation theory assuming that the mesh-size of the coarsest mesh is sufficiently small. In general, however, the coupling boundary is not smooth. Therefore, the preconditioner theory must not rely on compactness of .

To state the contributions of the current work, we consider the non-symmetric stiffness matrix of the (stabilized) Johnson-Nédélec coupling, which reads in block-form


see Section 2 below. Here, denotes an appropriate stabilization vector, which ensures positive definiteness of . The matrix block is the (positive semi-definite) Galerkin matrix of the FEM part, and the matrix block is the Galerkin matrix of the simple-layer integral operator . As in [FS09, MS98], we deal with block-diagonal preconditioners of the form


Here, the appropriate operator induces a coercive, symmetric, and bounded bilinear form . The symmetric and positive definite matrices resp. are spectrally equivalent to the symmetric and positive definite Galerkin matrices resp. , i.e.


where is strongly related to the FEM block of the FEM-BEM system (1), see (23)–(25) below. Inspired by [MS98], we prove that the condition number of as well as the number of iterations to reduce the relative residual by a factor in the preconditioned GMRES algorithm with inner product depends only on .

Usually, the condition number of Galerkin matrices and on adaptively refined meshes hinges on the global mesh-ratio as well as on the number of degrees of freedom. Therefore, the construction of optimal preconditioners for iterative solvers is a necessary task. Here, optimality is understood in the sense that the condition number of the preconditioned matrix is independent of the mesh-size and the degrees of freedom.

Very recently, it was proven in [XCH10] for 2D FEM with energy space that a local multilevel additive Schwarz preconditioner is optimal, i.e. the constants in (3) are independent of the mesh-size, the number of degrees of freedom, and the number of levels. Here, “local” means, that scaling at each level is done only on newly created nodes plus neighbouring nodes, where the associated basis functions have changed. An analogous result for 2D and 3D hypersingular integral equations with energy space has been derived by the authors [FFPS13]. In [FFPS13, XCH10], the proofs rely on a stable space decomposition of the discrete subspaces in resp. . Alternatively, [XCN09] provides stable subspace decompositions in for higher order elements in any dimension on bisection grids.

In the present work, we prove the optimality of some local multilevel additive Schwarz preconditioner for 2D weakly-singular integral equations with energy space . The proof is derived by postprocessing of the corresponding result for the hypersingular integral equation [FFPS13]. Combining this with the result of [XCH10], we prove optimality of for the FEM-BEM coupling.

Notation. Throughout the work, we explicitly state all constants and their dependencies in all statements of results. In proofs, however, we use the notation to abbreviate with a constant which is clear from the context. Moreover, abbreviates . Furthermore, the entries of a vector or a matrix are denoted by resp. . By , we denote the duality brackets between a Hilbert space and its dual .

Outline. The remainder of this work is organized as follows: In Section 2, we recall the basic facts on the Johnson-Nédélec coupling and define the admissible mesh-refinement strategies. We also formulate the preconditioned GMRES Algorithm 3, which is required to state the main result (Theorem 4). In Section 3, we prove the spectral estimates (3). Section 4 adapts the analysis for the symmetric coupling [MS98] and contains the proof of the main result (Theorem 4). The short Section 5 deals with extensions of the developed theory to the symmetric coupling and the one-equation Bielak-MacCamy coupling. Numerical examples from Section 6 underline our theoretical predictions and conclude the work.

2. Johnson-Nédélec coupling and main result

2.1. Model problem and analytical setting

Let be a polygonal and simply connected domain with boundary . We consider the following Laplace transmission problem in free space:


Here, denotes the outer normal on , and satisfies with uniform bounds on the maximal resp. minimal eigenvalue


With resp. and its dual , we denote the usual Sobolev spaces on resp. . For given data , , and , it is well-known that the model problem (4) admits a unique solution , with finite energy , if we impose the compatibility condition


to ensure the radiation condition (4e). Here, stands for the inner product, whereas denotes the extended inner product.

2.2. Johnson-Nédélec coupling

For the formulation of the Johnson-Nédélec coupling [JN80], the exterior solution (4b) is formulated by Green’s third formula. The latter gives rise to the simple-layer and double-layer integral operator

where denotes the fundamental solution of the 2D Laplacian and is the normal derivative. Note that boundedness holds for all . In addition, our analysis requires the hypersingular integral operator

where the integral is understood as finite part integral. It is known that and are symmetric in the sense of


for all and . Moreover, the assumption ensures -ellipticity of ,


Finally, is semi-elliptic with kernel being the constant functions,


The constants depend only on .

With the definitions (7) of the layer integral operators, the model problem (4) is equivalently recast by the Johnson-Nédélec coupling: Given , , and , find such that


for all . Setting in (11), we see that the (non-stabilized) linear operator associated to the left-hand side of (11) is indefinite. However, let denote the contraction constant of the double-layer potential [SW01]. Following the analysis in [OS13, Say09, Ste11], also (stabilized) Galerkin formulations of (11) admit unique solutions, if the ellipticity constant from (5) satisfies


and if the equation is either explicitly stabilized [OS13, Ste11] or if the discrete subspace of contains the constant functions [Say09]. In [AFF13a], the result of [Say09] is reproduced with a new proof. Introducing the notion of implicit stabilization, an equivalent elliptic operator equation of (11) is derived which fits in the frame of the Lax-Milgram lemma and thus leads to non-symmetric, but positive definite Galerkin matrices. The main result of [AFF13a] reads as follows (and also holds for a strongly monotone, but nonlinear material tensor ):

Lemma 1.

Let (12) be satisfied. For , define

as well as

Let and be arbitrary closed subspaces with . Then, the pair solves the Galerkin formulation of the Johnson-Nédélec coupling (11)


for all , if and only if it solves the operator formulation


for all . The operator is non-symmetric, but linear, continuous, and elliptic, and the constants


satisfy and depend only on and from (5) as well as on . Moreover, is a continuous linear functional on . In particular, (15) (and hence also (14)) admits a unique solution , and there holds the Céa-type estimate


where denotes the unique solution of the Johnson-Nédélec coupling (11).∎

2.3. Adaptive mesh-refinement and discrete spaces


T0 \psfragT1 \psfragT2 \psfragT3 \psfragT4 \psfragT12 \psfragT34

Figure 1. For each triangle , there is one fixed reference edge, indicated by the double line (left, top). Refinement of is done by bisecting the reference edge, where its midpoint becomes a new node. The reference edges of the son triangles are opposite to this newest vertex (left, bottom). To avoid hanging nodes, one proceeds as follows: We assume that certain edges of , but at least the reference edge, are marked for refinement (top). Using iterated newest vertex bisection, the element is then split into 2, 3, or 4 son triangles (bottom).

Let be a given conforming initial triangulation of into compact and non-degenerate triangles. We suppose that a sequence of refined triangulations is obtained by newest vertex bisection, see Figure 1, where is the coarsest conforming mesh such that all marked elements have been bisected. For our analysis, the set of marked elements is arbitrary, but in practice obtained from local a posteriori refinement indicators, see e.g. [AFF13a]. We note that newest vertex bisection guarantees uniform shape regularity in the sense that


where depends only on the initial mesh , see e.g. [Ver13, KPP13] and the references therein.

Let be a given initial partition of the coupling boundary into compact line segments. We suppose that a sequence of refined partitions is obtained by bisection, where the refined elements are refined into two sons of half length, i.e. with and where at least the marked elements are refined, i.e. . In addition, we suppose that the meshes are uniformly -shape regular in the sense that


where depends only on the initial partition . Possible choices include the 1D bisection algorithms from [AFF13b]. A further choice is to consider the partition of which is induced by the triangulation of . Formally, such a coupling of and is not required for the analysis, but simplifies the implementation and is therefore used in the numerical experiments of Section 6.

In this work, we consider lowest-order Galerkin elements. We approximate functions by functions and functions by functions , where


Let denote the set of nodes of the triangulation , and let denote the set of nodes of the triangulation . For resp. , we define the patch resp. . For the construction of optimal multilevel preconditioners on adaptively refined triangulations, we need the following subsets of resp. :


The sets resp. thus consist of the new nodes plus the old nodes, where the corresponding patches have changed. A visualization of is given in Figure 2. For each node , denotes the associated hat-function with for all , where is Kronecker’s delta.

Figure 2. The left figure shows a FEM mesh , where the two elements (green) are marked for refinement. Bisection of these two elements provides the mesh (right), where two new nodes are created. The set consists of these new nodes plus their immediate neighbours (red), where the corresponding patches have changed. The union of the support of the basis functions associated to nodes in is given by the light- and dark-green areas in the right figure.

2.4. Galerkin system and block-diagonal preconditioning

Let with denote the nodal basis of , and let with denote a basis of , where denotes the characteristic function on . The Galerkin matrix of the operator has the form

where the block matrices , , , and as well as the stabilization (column) vector are given by


We stress that as well as are sparse, whereas is dense. Note that the number of non-zeros in the matrix is bounded by . Moreover, the application of the rank- stabilization matrix can be implemented efficiently with complexity for use with an iterative solver.

The discrete variational formulation (15) is equivalent to solving the following linear system of equations: Find such that


where the right-hand side vector reads

To formulate our block-diagonal preconditioner, we require an appropriate operator which is related to the FEM-domain part of . The next lemma follows from a Rellich compactness argument, since for all constants . Details are analogous to, e.g., [AFF13a, Lemma 10] and therefore left to the reader.

Lemma 2.

For , define

Then, the operator is linear, symmetric, continuous, and elliptic, and the constants

satisfy and depend only on and from (5) as well as on .∎

In this work, we investigate block-diagonal preconditioners of the form


where is a “good” approximation of the Galerkin matrix corresponding to the operator from Lemma 2 with respect to the nodal basis of , and is a “good” approximation of the Galerkin matrix . Our construction below ensures that and hence are symmetric and positive definite.

Instead of (24), we solve the preconditioned system


For this non-symmetric system of linear equations, we use a preconditioned GMRES algorithm [SS86], which will be discussed in the following subsection. The preconditioned Galerkin matrix reads in block form


2.5. Preconditioned GMRES algorithm

Let denote a symmetric and positive definite matrix and let denote a (possibly) non-symmetric, but positive definite matrix. Let, denote the standard unit vector with entries . The preconditioned GMRES algorithm reads as follows.

Algorithm 3 (Gmres).

Input: Matrices , right-hand side vector , initial guess , relative tolerance , and maximum number of iterations with .

  • Allocate memory for the matrix , the vectors , , , , and .

  • Compute initial residual and .

  • Set counter , and initialize for all .

Iterate the following steps (i)–(vii):

  • Compute .

  • For all compute

  • Compute .

  • Define the sub-matrix with entries for and compute

  • Compute and .

  • Stop iteration if or .

  • Otherwise, compute , update counter , and goto (iv).

Output: , , and .

Note that for being the identity matrix, Algorithm 3 is the usual GMRES algorithm with inner product , see e.g. [SS86]. The main memory consumption is given by the vectors , while the matrix in step (iv) is a sub-block of the matrix and thus does not need to be stored explicitly.

As is often the case for multilevel preconditioners, the application of on a vector is known, whereas the application of is unknown. We therefore note that the GMRES Algorithm 3 can be implemented without using to compute the inner products and the norms . To this end, one replaces the computation of by , where . Then, . In step (i), we compute instead of . Step (ii) is replaced by and . Instead of step (iii), we then compute . Note that in step (vi). Lastly, in step (vii) we replace the computation of by and .

2.6. Local multilevel preconditioner and main result

For both the FEM part and BEM part in (26), we will use local multilevel preconditioners which are optimal in the sense that the condition numbers of the preconditioned systems are independent of the number of levels and the mesh-size .

For the preconditioner , we use an additive Schwarz multilevel diagonal preconditioner similar to the one in [XCH10]. Recall from (22) Define


Let denote the canonical embedding with matrix representation and . Furthermore, let denote the diagonal of the Galerkin matrix with respect to the local set of nodes , i.e. for and . Then, our local multilevel diagonal preconditioner is defined via


From the definition, we see that this preconditioner corresponds to a diagonal scaling on each level, where scaling is done on the local subset only.

For all boundary nodes , let


denote the boundary hat-function with for all . To construct an efficient preconditioner for the weakly-singular integral operator in 2D, we use the Haar-basis functions for all . Recall from (22). Let . Define the local subspaces