
==== Front
BIT Numer Math
BIT Numer Math
Bit. Numerical Mathematics
0006-3835
1572-9125
Springer Netherlands Dordrecht

1034
10.1007/s10543-024-01034-9
Article
Block diagonal Calderón preconditioning for scattering at multi-screens
http://orcid.org/0000-0001-9816-6321
Cools Kristof 1
Urzúa-Torres Carolina CarolinaUrzua-Torres@tudelft.nl

2
1 Tech Lane 126, 9052 Ghent, Belgium
2 Mekelweg 6, 2628 CD Delft, The Netherlands
3 9 2024
3 9 2024
2024
64 4 346 1 2023
26 7 2024
© The Author(s) 2024
2024
https://creativecommons.org/licenses/by/4.0/ Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.
A preconditioner is proposed for Laplace exterior boundary value problems on multi-screens. To achieve this, the quotient-space boundary element method and operator preconditioning are combined. For a fairly general subclass of multi-screens, it is shown that this approach paves the way for block diagonal Calderón preconditioners which achieve a spectral condition number that grows only logarithmically with decreasing mesh size, just as in the case of simple screens. Since the resulting scheme contains many more degrees of freedom than strictly required, strategies are presented to remove almost all redundancy without significant loss of effectiveness of the preconditioner. The performance of this method is verified by providing representative numerical results. Further numerical experiments suggest that these results can be extended to a much wider class of multi-screens that cover essentially all geometries encountered in practice, leading to a significantly reduced simulation cost.

Keywords

Preconditioning
Complex screens
Galerkin boundary element method
Quotient-space boundary element method
Mathematics Subject Classification

22E46
53C35
57S20
http://dx.doi.org/10.13039/100019180 HORIZON EUROPE European Research Council 101001847 Cools Kristof NWOVI.Veni.212.253 Urzúa-Torres Carolina issue-copyright-statement© BIT Foundation 2024
==== Body
pmcIntroduction

We are interested in the behaviour of potentials near multi-screens, which are geometries composed of essentially two-dimensional piecewise smooth surfaces joined together, as shown in Fig. 1. Hence, we consider the following Dirichlet and Neumann Laplace boundary value problems (BVPs) in the exterior of the multi-screen Γ⊂R3,1 -ΔU=0inR3\Γ¯,U=gDor∂U∂n=fNonΓ,

plus the decay condition U(x)=O(‖x-1‖)as‖x‖→∞, where ‖x‖ designates the Euclidean norm of a point x in R3, O(·) the Landau symbol, and gD and fN are suitable boundary data.Fig. 1 Two examples of multi-screen geometries

Our goal is to solve these exterior BVPs efficiently by means of Galerkin boundary element methods (BEM) [29] and Calderón preconditioning [7, 31]. For this, we recast the BVPs as variational first-kind boundary integral equations (BIEs) for densities on the surface of the multi-screen.

For simple screens this approach is well established [29, Section 3.5.3]. Here, we call a simple screen an orientable, piecewise smooth two-dimensional manifold with boundary embedded in R3. For these geometries, the arising variational first-kind BIEs are known to be coercive [16, 17, 32] in Sobolev spaces of jumps of suitable field traces, in H~-12(Γ) and H~+12(Γ), respectively [27, Ch. 3]. For these trace spaces, conforming boundary element spaces are easily available, and they lead to Galerkin approximations and Calderón preconditioning whose numerical analysis is well-understood [23, 24, 28].

In contrast, the notion of jumps becomes problematic in multi-screens, since they are not globally orientable. For this reason, the tools from simple screens cannot be used straightforwardly on multi-screens. Many alternatives have been proposed to tackle this problem [5, 11–14, 35]. It is worth pointing out that at the time of writing, a rigorous analysis of these approaches in suitable trace spaces is not available. Furthermore, these approaches lead to ill-conditioned linear systems, yet are not amenable to preconditioning.

Fortunately, recent work by Claeys and Hiptmair offers the mathematical framework to overcome these difficulties [9]. The key idea is to see trace spaces from the perspective of quotient-spaces and to work with multi-valued traces. This new paradigm not only allows for a rigorous analysis, but it also paves the way for conforming Galerkin discretisation by means of quotient-space BEM, as proposed in [8]. Indeed, instead of trying to approximate jumps directly, the new approach relies on the Galerkin discretisation of multi-trace boundary element spaces. With this approach, the related BIEs give rise to Galerkin matrices with large null spaces comprised of single-trace functions. Since the right-hand-sides of the linear systems of equations are consistent, Krylov subspace iterative solvers like CG still converge to the right solution. We summarise these ideas and results in Sect. 2.

Now that the most fundamental issues have been solved, we are in the position to investigate how to improve the computational performance of quotient-space BEM for multi-screens. Indeed, one should note that the arising linear systems are ill-conditioned and that the number of conjugate gradient (CG) iteration counts increases with mesh refinement. Hence, a natural next step—and the main focus of this paper—is to devise preconditioners for multi-screen problems. In Sect. 4, we propose a simple preconditioning strategy based on opposite-order preconditioning, also known as Calderón preconditioning on closed surfaces. Moreover, we present the tools to understand the new preconditioner in the context of operator preconditioning. Numerical experiments confirm that this approach reduces considerably the number of CG iterations required to solve the system.

It is worth mentioning that an advantage of the quotient-space BEM approach is that minimal geometrical information is required. However, the disadvantage is that one pays with unnecessary computations due to the “doubling of degrees of freedom” underlying the discretisation of multi-valued traces. As an alternative, we dedicate Sect. 5 to discuss reduced quotient-space representations that require slightly more geometrical information, but lower computational effort while still rendering efficient Calderón preconditioning. Furthermore, we use the tools derived in Sect. 4 to provide some insight about the requirements that such reductions need to fulfil.

Last but not least, we should mention that another approach to preconditioning multi-screens has become available during the revision of this article [2]. We believe this confirms the problem at hand is relevant and that, as usual in mathematics, there are different ways of tackling a problem.

Quotient-space perspective

We briefly summarise the new perspective introduced in [9, Section 4-6] and the quotient-space construction of boundary element (BE) spaces from [8]. Throughout this paper we focus on three dimensions, but it is worth mentioning that the method and analysis also carries over to 2d.

Geometry

We begin by recalling the rigorous characterisation of multi-screens as given in [9, Sect. 2].

For this, the first concept we need to introduce is that of a Lipschitz (simple) screen in the sense of Buffa-Christiansen:

Definition 1

(Lipschitz Screen [9, Definition 2.1]) A Lipschitz screen is a subset Γ⊂R3 such thatits closure Γ¯ is a compact Lipschitz two-dimensional sub-manifold with boundary,

Writing ∂Γ for the the boundary of Γ¯, we have that Γ=Γ¯\∂Γ,

there exists a finite covering C of Γ¯ with cubes such that for each cube C∈C, denoting by a the length of its sides, the following holds:▪ If C contains a point in ∂Γ, there is an origin and an orthonormal basis of R3 in which the cube C can be identified with (0,a)3 and there are uniformly Lipschitz continuous functions ψ:R→R, and ϕ:R2→R, with values in (0, a) such that: 2a Γ∩C={(x,y,z)∈C:y<ψ(x),z=ϕ(x,y)},

2b ∂Γ∩C={(x,y,z)∈C:y=ψ(x),z=ϕ(x,y)}.

▪ Otherwise, Γ is the graph above (0,a)2 of an uniformly Lipschitz continuous function ϕ:R2→R, with values in (0, a).

Definition 2

(Lipschitz Partition [9, Definition 2.2]) A Lipschitz partition of R3 is a finite collection of Lipschitz open sets Ωjj=0…n such that R3=∪j=0nΩ¯j and Ωj∩Ωk=∅, if j≠k.

Definition 3

(Multi-screen [9, Definition 2.3]) A multi-screen is a subset Γ⊂R3 such that there exists a Lipschitz partition of R3 denoted Ωjj=0…n satisfying Γ⊂∪j=0n∂Ωj and such that for each j=0…n, we have that the interior of Γ¯j:=Γ¯∩∂Ωj is a Lipschitz screen.

From a numerical point of view, it will be convenient to classify multi-screens into three categories. For this, we first need to introduce the notion of irregular points on the boundary, as in [10].

Definition 4

(Irregular points [10, Sect. 6.1]) Let us consider ∂Γ:=Γ¯\int(Γ) and introduce the set of regular points of the boundary PR(∂Γ) defined asPR(∂Γ)={x∈∂Γsuch thatBx∩Γ=Bx∩Sfor some ballBxcentred atxand some simple Lipschitz screenS}.

We define the set of irregular points of the boundary asPI(∂Γ)=∂Γ\PR(∂Γ).

With this, we can classify our multi-screens as follows:Type A: Γ is a multi-screen such that there exists an underlying Lipschitz partition (Ωj)j=0...n of R3 with the property that ∂Γ¯j⊆∂Γ with Γj as in Definition 3.

Type B: Γ is a multi-screen that has irregular points and that is not of type A.

Type C: Γ is a multi-screen without irregular points and that is not of type A.

Fig. 2 Multi-screens can be classified according to the location of their irregular points

Figure 2 provides examples of multi-screens in these three different classifications. In particular, Fig. 2c depicts a Möbius strip, which will not be discussed in this paper because its analysis is more cumbersome and it does not arise in applications. Indeed typical geometries that are approximated by multi-screens in applications are antennas, tail fins in aircrafts, and heat sinks. Since all of these correspond to multi-screens of Type A and Type B, we restrict ourselves to these two types.

Trace spaces

For a multi-screen Γ⊂R3 we consider the following chains of nested Sobolev spaces13a H0,Γ1(R3)⊂H1(R3)⊂H1(R3\Γ¯),

3b H0,Γ(div,R3)⊂H(div,R3)⊂H(div,R3\Γ¯),

where the subscript X0,Γ indicates a space obtained as the closure in X of smooth functions/vectorfields compactly supported in R3\Γ¯. All inclusions in (3) should be read as “is a closed subspace of”, which describe the associated quotient-spaces Hilbert spaces. With this, we can define the multi-trace spaces [9, Sect. 5] 4a H+12(Γ):=H1(R3\Γ¯)/H0,Γ1(R3),

4b H-12(Γ):=H(div,R3\Γ¯)/H0,Γ(div,R3),

and the single-trace spaces [9, Sect. 6.1] 5a H+12([Γ]):=H1(R3)/H0,Γ1(R3),

5b H-12([Γ]):=H(div,R3)/H0,Γ(div,R3).

Remark 1

We note that H1(R3\Γ¯) and H(div,R3\Γ¯) are spaces of functions attaining different values on both sides of Γ. This implies that functions in the multi-trace spaces H+12(Γ) and H-12(Γ) are multi-valued on Γ. In other words, they can take different values on both sides of Γ.

Since the spaces H+12([Γ]) and H-12([Γ]) are closed subspaces of H+12(Γ) and H-12(Γ), respectively [9, Proposition 6.2], we can also introduce the jump spaces [9, Sect. 6.2] as 6a H~+12([Γ]):=H+12(Γ)/H+12([Γ]),

6b H~-12([Γ]):=H-12(Γ)/H-12([Γ]).

Remark 2

It is worth mentioning that single-trace spaces are a generalisation of the spaces H±12(Γ) on simple screens. Indeed, definition (5) follows from the characterisation of trace spaces as the quotient between the domain of the Dirichlet and normal traces, and their kernels.

Then, multi-trace spaces are the counterpart of single-trace spaces when starting from the multi-valued spaces H1(R3\Γ¯) and H(div,R3\Γ¯).

Finally, if one defines jump operators [·]:H±12(Γ)→H~±12([Γ]) as in [9, Def. 6.5], one gets that their kernels are the single-trace spaces. This motivates the characterisation of jump-spaces as the quotients (6).

Next, we consider the canonical surjections7 πD:H1(R3\Γ¯)→H+12(Γ)andπN:H(div,R3\Γ¯)→H-12(Γ),

and, with H1(Δ,R3\Γ)={u∈H1(R3\Γ),Δu∈L2(R3)}, we define the relevant trace operatorsDirichlet trace:γD:H1(R3\Γ)→H+12(Γ),γD:=πD,Neumann trace:γN:H1(Δ,R3\Γ)→H-12(Γ),γN:=πN∘grad.

Moreover, we remark that they map onto H+12([Γ]) and H-12([Γ]) when restricted to H1(R3) and H1(Δ,R3), respectively.2

As noted in [9, Sect. 5.1], Green’s Formula in R3 does not hold for elements of H1(R3\Γ¯) and H(div,R3\Γ¯). As these spaces underlie the definitions of H+12(Γ) and H-12(Γ), that implies we cannot use the usual L2-duality pairing.

As a remedy, we introduce a bilinear pairing on H+12(Γ)×H-12(Γ):8 ≪u,p≫:=∫[Γ]updσ:=∫Rd\Γ¯p·∇u+udiv(p)dx,

with arbitrary representatives u∈H1(R3\Γ¯) and p∈H(div,R3\Γ¯) [9, Sect. 5.1].3 Note that this pairing induces the following isometric dualities [9, Prop. 5.1 and Sect. 6.2]H-12(Γ)≅H+12(Γ)′,H~-12([Γ])≅H+12([Γ])′,H~+12([Γ])≅H-12([Γ])′.

Furthermore, the bilinear pairing (8) offers a characterisation of single-trace spaces through self-polarity:

Proposition 1

([9, Proposition 6.3]) For u∈H+12(Γ) and p∈H-12(Γ) the following equivalences hold true:u∈H+12([Γ])⟺≪u,q≫=0∀q∈H-12([Γ]),p∈H-12([Γ])⟺≪v,p≫=0∀v∈H+12([Γ]).

Remark 3

These polarity properties may seem surprising starting from the quotient space definition of the multi-trace spaces. However, thinking of the interpretation of the pairing (8) as an L2-type pairing on an inflated multi-screen, they make sense: u∈H+12([Γ]) is even when crossing Γ and p∈H-12([Γ]) is odd when crossing Γ (because the normal changes direction). The result is that contributions from opposite sides of Γ cancel.

Weakly singular and hypersingular BIEs

Let G(z):=14π‖z‖ be the fundamental solution of the Laplace equation in R3. For x∉Γ, let Gx(y):=χx(y)G(x-y) with χx:R3→R a smooth cut-off function that is 1 in a neighborhood of Γ and 0 in a neighborhood of x as in [9, Sect. 8]. This allows the definition of the single and double layer potentials by9 SLϕ(x):=≪γDGx,ϕ≫,DLv(x):=-≪γNGx,v≫.

The weakly singular and hypersingular boundary integral operator (BIO) are the continuous operators V0:=γD∘SL:H-12(Γ)→H+12(Γ),W0:=γN∘DL:H+12(Γ)→H-12(Γ). For sufficiently smooth arguments, the corresponding bilinear forms admit weakly singular representations given by [8, Sect. 3]10 ≪V0ϕ,ψ≫=∫[Γ]∫[Γ]G(y-x)ϕ(y)ψ(x)dσ(y)dσ(x),

11 ≪W0v,p≫=∫[Γ]∫[Γ]G(y-x)curlΓv(y)·curlΓp(x)dσ(y)dσ(x),

In order to solve the Dirichlet Laplace BVP, we solve the following variational BIE: Given gD∈H+12([Γ]), find ϕ∈H-12(Γ) such that12 ≪V0ϕ,ψ≫=≪gD,ψ≫∀ψ∈H-12(Γ).

To solve the Neumann Laplace BVP, we solve the variational BIE: Given fN∈H-12([Γ]), find u∈H+12(Γ) such that13 ≪W0v,p≫=≪fN,p≫∀p∈H+12(Γ).

We conclude this section by recalling some properties of these BIEs: First, as a consequence of the polarity from Proposition 1, we get that

Lemma 1

([8, Lemma 3.2]) The nullspaces of V0 and W0 agree with H-12([Γ]) and H+12([Γ]), respectively.

In analogy to the situation on simple screens, we have

Proposition 2

The operators V0:H~-12([Γ])→H+12([Γ]) and W0:H~+12([Γ])→H-12([Γ]) are elliptic, i.e.14 ≪V0q,q≫≥αV‖q‖H~-12([Γ])2∀q∈H~-12([Γ]),

15 ≪W0v,v≫≥αW‖v‖H~+12([Γ])2∀v∈H~+12([Γ]),

with αV,αW>0 depending only on Γ.

Proof

The proof is similar to that of [9, Prop. 8.7] but setting ψ to be the single layer potential for the Laplacian and using the fact that Δψ=0 and that |ψ|H1(R3\Γ¯) is equivalent to ‖ψ‖H1(Δ,R3\Γ¯) [29, Thm. 2.10.10]. □

The previous results combined with continuity of the operators give us

Proposition 3

([9, Prop. 8.9]) The operators V0:H~-12([Γ])→H+12([Γ]) and W0:H~+12([Γ])→H-12([Γ]) are isomorphisms.

Additionally, these operators remain well-defined on the multi-trace spaces H-12(Γ) and H+12(Γ), respectively. However, Lemma 1 implies that they have non-trivial nullspaces when considered on multi-trace spaces. Although this excludes uniqueness of solutions for (12) and (13), Proposition 1 still provides existence, since gD∈H+12([Γ]) and fN∈H-12([Γ]) guarantees consistency of the right-hand side linear forms: they vanish on the single-trace spaces.

Operator preconditioning on quotient-space BEM

In order to explain what changes in the quotient-space BEM setting, we recall the essential ingredients of operator preconditioning as presented in [22]: Let X, Y be Banach spaces, and consider the finite-dimensional subspaces Xh⊂X and Yh⊂Y with dimensions N:=dimXh and M:=dimYh, and bases (φi)i=1N and (ϕi)i=1M, respectively. Further, let a∈L(X×X,R), b∈L(Y×Y,R) and m∈L(X×Y,R) be continuous bilinear forms satisfying discrete inf-sup conditions:16 supvh∈Xh|a(uh,vh)|‖vh‖X≥αA‖uh‖X,∀uh∈Xh,

17 supwh∈Yh|b(qh,wh)|‖wh‖Y≥αB‖qh‖Y,∀qh∈Yh,

18 supwh∈Yh|m(vh,wh)|‖wh‖Y≥αM‖vh‖X,∀vh∈Xh.

If N=M, then [22, Theorem 2.1] implies that the associated Galerkin matrices Ah:=a(φi,φj)i,j=1N,Bh:=b(ϕi,ϕj)i,j=1N,Mh:=m(φi,ϕj)i,j=1N, satisfy19 κsp(Mh-1BhMh-TAh)≤‖a‖‖b‖‖m‖2αAαBαM2,

where κsp designates the spectral condition number and ‖·‖ denotes the corresponding operator norms.

From (19), the idea is that if we are solving (12), the Galerkin matrix of V0 would play the role of Ah and we would need to find a suitable bilinear form b such that the spectral condition number is as small as possible if we precondition the resulting linear system with Ph=Mh-1BhMh-T. Analogously, if we are solving (13), the Galerkin matrix of W0 would play the role of Ah and we would need a different choice of b.

Hence, the first question is: which Galerkin matrix one should consider? In other words, what discretisation should we choose? The answer following quotient-space BEM is to discretise the multi-trace spaces. However, due to Lemma 1, the bilinear forms of V0 and W0 will not satisfy its inf-sup condition on their discrete multi-trace space. We therefore have to extend (19) to use operator preconditioning for (12) and (13).

As usual, bounding the spectral condition number entails bounding the largest eigenvalue λmax:=λmax(Mh-1BhMh-TAh) from above and the smallest λmin:=λmin(Mh-1BhMh-TAh) from below. However, when using quotient-space BEM, we have to consider the spectral condition number away from the kernel of Ah, i.e., κ~sp(Mh-1BhMh-TAh)=λmax/λmin~, where λmin~ is the smallest non-zero eigenvalue of Mh-1BhMh-TAh, since this quantity determines the convergence of CG on singular systems [20, 26]. In order to write these bounds, we need to introduce some notation first.

Let Ah:Xh→Xh′, Bh:Yh→Yh′ and Mh:Xh→Yh′ be the bounded linear operators associated to the bilinear forms a, b and m, respectively.

For λmax we proceed in the classical way and arrive to20 λmax≤‖Mh-1‖2‖Bh‖‖Ah‖=αM-2‖Bh‖‖Ah‖.

For λ~min we have to take a slightly different approach since we need to restrict Ah to the space where its corresponding bilinear form a satisfies a discrete inf-sup condition. Moreover, we have to establish when the discrete inf-sup condition will bound the smallest eigenvalue. We study this in the next Lemma.

Lemma 2

(Discrete inf-sup constant in the quotient space norm) Let X be a Hilbert space. Let a be a continuous bilinear form on X×X. Let X⊆X be both the left and right nullspace of a. Let Xh be a finite dimensional subspace of X and Ah:Xh→Xh′ the bounded linear operator associated to a. We assume that Xh is nullspace conforming to a in the sense that Xh:=kerAh=kerAh′ is a linear subspace of X∩Xh⊆X.

If a satisfies a discrete inf-sup condition in Xh/X×Xh/X with constant αa>0 and if the norms on X/X and Xh/Xh are equivalent, i.e. there exists ceq>0, such that for all uh∈Xh21 ‖uh‖X/X≤‖uh‖X/Xh≤ceq‖uh‖X/X,

then22 supvh∈Xh\{0}a(uh,vh)‖vh‖X≥αaceq‖uh‖X/X.

Proof

We have, with Xh⊥={uh∈Xh,∀vh∈Xh,(uh,vh)X=0},23 supvh∈Xh\{0}a(uh,vh)‖vh‖X≥supv~h∈Xh⊥\{0}a(uh,v~h)‖v~h‖X.

For v~h∈Xh⊥, we have ‖v~h‖X=‖v~h‖X/Xh (e.g. [6, Eq. (4)]), and so24 supv~h∈Xh⊥\{0}a(uh,v~h)‖v~h‖X=supv~h∈Xh/Xh\{0}a(uh,v~h)‖v~h‖X/Xh,

where we have identified Xh⊥ and Xh/Xh. Finally, by the norm equivalence (21) and the discrete inf-sup condition, we have that25 supvh∈Xh\{0}a(uh,vh)‖vh‖X≥ceq-1supv~h∈Xh/Xh\{0}a(uh,v~h)‖v~h‖X/X≥αaceq‖uh‖X/X.

□

Note that inequality (21) holds in our application by virtue of Lemma 7.

Using Lemma 2 and the properties of the related bilinear forms, we can bound λ~min as follows26 λ~min≥infuh∈Xh/Xh\{0}‖Mh-1BhMh-∗Ahuh‖X‖uh‖X≥αM-12αBαAceq≥αBαA‖M‖2ceq.

Finally, we combine (20) and (26) and get27 κ~sp(M-1BhM-TAh)≤ceq‖B‖‖A‖‖M‖2αBαAαM2.

Remark 4

For the discrete quotient norm to be bounded from below by the continuous quotient norm, it suffices that Xh⊆X∩Xh; equality is not required. This opens the door to reduction schemes for the multi-trace space variational formulation. A reduction scheme is a choice Xh∘⊂Xh⊂X with corresponding discrete left/right nullspace Xh∘ such that (i) Xh∘⊂X and (ii) Xh∘/Xh∘=Xh/Xh. Such a choice leads, on the one hand, to approximations in the jump space of equal quality, and, on the other hand, does not preclude the construction of efficient operator preconditioners. We will discuss this again in Sect. 5.

Calderón preconditioning for multi-screens

As already mentioned in the introduction, the linear systems arising from the discretisation of (12) and (13) using Quotient-space BEM are ill-conditioned, which causes that the number of CG iteration counts increases with mesh refinement. One should note that this is not a particularity of Quotient-space BEM. Indeed, we usually encounter this difficulty when using low-order BEM discretisation of first-kind integral equations on simple screens and closed surfaces. In those cases, one typically remedies the problem by using so-called Calderón preconditioning, which combines Calderón identities with operator preconditioning to build a convenient and effective preconditioner [7, 22, 31].

In this paper, we will extend this approach and devise Calderón preconditioners for the problem at hand. For the sake of presentation, we will discuss the details for the case of the Hypersingular operator W0 and note that one can proceed analogously to precondition the weakly singular operator V0.

Discretisation

When considering the Neumann problem (13), we want to find solutions v in H~+12([Γ]). The idea of Quotient-space BEM is to discretise the multi-trace space H+12(Γ) instead.

As mentioned in Remark 1, multi-trace spaces can take different values on both sides of Γ. One way to grasp this is to imagine an “infinitesimally inflated” screen, as illustrated in Fig. 3 for a 2D multi-screen. With this, one can intuitively understand the trace spaces introduced in Sect. 2.2 as follows:H+12(Γ) can be seen as a standard Dirichlet trace space on the surface of the “inflated screen”. Similarly, H-12(Γ) can be viewed as the standard Neumann trace space on the surface of the “inflated screen”.

The single-trace space H+12([Γ]) simply consists of single-valued functions on Γ. One can follow the same intuition for H-12([Γ]), however, its right interpretation as a single-valued normal component requires that one fixes a local normal n on Γ.

It is worth noting that this depiction via the “inflated screen” is only meant to help the reader to visualise the related spaces. Furthermore, every time we resort to this intuition, the thickness of the “inflated screen” is indeed zero.Fig. 3 “Inflating” a 2D multi-screen

Fig. 4 Illustration of virtual mesh for 2D multi-screen (please note that here the thickness of the virtual meshes is non-zero only to simplify the picture.)

Following this intuition, a simple way to implement the multi-trace space in an already existing BEM code, is via a triangular virtual surface mesh of Γ, as defined in [8, Sect. 4.1]. In essence, this virtual mesh is built from a triangulation T0 of Γ. Then, for every triangle K∈T0 one creates two copies K+ and K- with the same geometry but to be regarded as different entities. The reader may imagine K+ and K- as the two faces of K. In addition, one defines nK as the unit normal vector of K, and labels the copies such that nK points from K- to K+. In this way, K+ is endowed with nK and -nK is assigned to K-. Finally, the union of these oriented triangles forms the set underlying what we call the virtual surface mesh of Γ. Figure 4 illustrates this idea for a 2D multi-screen (where triangles are replaced by segments of the mesh). We refer to [8, Sect. 4.1] for further details.

Let Th be a triangular virtual surface mesh of Γ, with target element size h, and let Tˇh be its dual as realised on the barycentric refinement [4]. The BE spaces above could be chosen as [30, Sect. 2.2]Xh(Γ)=S1,0(Th): piecewise linear “continuous” functions on Th,

Yh(Γ)=S0,-1(Tˇh): piecewise constant functions on Tˇh,

yet we will need a different choice for our preconditioner to work, as we will discuss in the next section.

Implementation

Note that by construction of the virtual mesh, which can be understood as a triangulation of a closed surface (see Fig. 4c), the space S1,0(Th) has one degree of freedom at the vertices in Th∩∂Γ. Since solutions for the hypersingular equation (13) belong to H~+12([Γ]), we know they will be zero on ∂Γ [8]. Indeed, because functions in H~+12([Γ]) are only determined up to contributions in H+12([Γ]), degrees of freedom on ∂Γ (which by construction are in H+12([Γ])) can be safely deleted. Hence, instead of working with S1,0(Th), we consider S01,0(Th)⊂H+12(Γ): piecewise linear “continuous” functions on the inflated screen that are zero on ∂Γ, as shown in Fig. 4d.

When dealing with multi-screens of type A, this will have the computational advantage of allowing us to decouple the BE spaces on each side of the triangular virtual surface mesh Th, as depicted in Figs. 4d and 5.

Let us illustrate how we implemented these BE spaces on a multi-screen Γ that is the union of three simple screens Mi,i=1,2,3 meeting at a junction: We define the inflated multi-screen [Γ] as 28 [Γ]=∪l=13Il

with Il=Ml∪Ml+1. The normal on Il is chosen outward. Each simple screen Mi appears once as the front and once as the back of the multi-screen (see Fig. 5).

For i=1,2,3, we create the triangular surface mesh Mi,h of Mi with target element size h, and such that the meshes Mi,h for i=1,2,3 match up along the junction. The simple screens Il inherit this mesh. In other words, we have Il,h=Ml,h⋃Ml+1,h,l=1,2,3.

The discrete primal multi-trace space Xh(Γ) is built as the product of these spaces, i.e. 29 Xh(Γ)=S01,0(Th)=⨂l=13S01,0(Il,h).

Construct the dual BE spaces on the simple screens Il following the cue from [4]. For this, let Iˇl,h denote the dual barycentric mesh to Il,h, built as in [25, Definition 2]. Then we introduce the space S0,-1(Iˇl,h)⊂H-12(Il) of piecewise constant functions supported by the dual cells of Iˇl,h that correspond to nodes not on the boundary of Il. In particular we have that dimS0,-1(Iˇl,h)=dimS01,0(Il,h).

The discrete dual multi-trace space Yh(Γ) is built as the product of these dual spaces, i.e. 30 Yh(Γ)=S0,-1(Tˇh)=⨂l=13S0,-1(Iˇl,h).

By construction, we have that N:=dim(Yh(Γ))=dim(Xh(Γ)).

Fig. 5 Back-front (decoupled) conforming mesh on the multi-screen

Remark 5

The description of discrete multi-trace spaces (29) and (30) is not valid for multi-screens of Types B and C.

Block diagonal Calderón preconditioner

Based on Calderón preconditioning for closed surfaces and its applicability to simple screens, one could think of preconditioning the hypersingular operator W0 with the weakly singular operator V0. However, it is clear from Lemma 1 that V0 will not do the job.

As an alternative, we propose to use Calderón preconditioning blockwise, under the considerations of the previous sections. Let Xh(Γ)=span{φk}k=1N and consider the Galerkin matrix for the hypersingular operator, i.e.,31 Wh[i,j]=≪W0φj,φi≫,i,j=1,⋯,N.

Then, we will build a preconditioner for Wh based on the block matrix32 BhV:=Vˇh,1000Vˇh,2000Vˇh,3,

where Vˇh,l[i,j]=⟨V0ψˇj,ψˇi⟩Il with ψˇi,ψˇj in the standard basis of S0,-1(Iˇl,h) for l=1,2,3, and ⟨·,·⟩Il denotes the usual L2-duality pairing over Il.

The motivation to consider this BhV is that the choice of discrete spaces from (29) and (30) allows us to decouple what is happening on the dual space of each simple screen Il. Furthermore, they would agree with the standard discretisation of the jump spaces on simple screens. More concretely, we have that S01,0(Il,h)⊂H~12(Il) and S0,-1(Iˇl,h)⊂H~-12(Il).

In order to analyse the impact of this preconditioner in the number of PCG iteration counts, we aim to bound the resulting condition number. For this, we first need to make some assumptions.

Let us begin by noticing that the Lipschitz partition Ωjj=0…n such that Γ⊂∪j=0m∂Ωj is not unique. We illustrate this for a two-dimensional multi-screen Γ with a triple junction in Fig. 6. Nevertheless, for the multi-screens considered in this paper, one can always find a Lipschitz partition such that the Lebesgue measure meas(Γ∩∂Ωi)>0 for all i=0,⋯,m. In other words, we can always assume we have the configuration corresponding to Fig. 6b. For simplicity of the proofs, this is the type of Lipschitz partitions that we will consider. This and the particular order of the domains is stated in the following:Fig. 6 Example of two Lipschitz partitions for Γ being a 2D multi-screen with a triple junction

Assumption 1

Let Ωjj=0…n be a Lipschitz partition in R3. We assume Γ to be a multi-screen such that Γ⊂∪j=0m∂Ωj and meas(Γ∩∂Ωi)>0 for all i=0,⋯,m. Moreover, and without loss of generality, we assume that Γ and the Lipschitz partition Ωjj=0…n are such that Ω0 is the only unbounded domain.

Next, let us recall the notation Γi=Γ∩∂Ωi for i=0,⋯,m. With this, we state one more assumption that we will need in order to show the required condition number bound:

Assumption 2

Let Γ be a multi-screen as in Assumption 1. We assume that for each j=0,..,m, the family of meshes {Γjh}h∈H,h>0 of Γj:agree at the junction(s);

are uniformly shape-regular, and locally quasi-uniform;

satisfy the following local mesh condition: For each triangle τl∈Γjh, we define the index set J(l):={k∈{1,⋯,M}:τl∩supp(φk)≠0}, where M:=dim(S01,0(Γjh)). Moreover, for each basis function φk∈S01,0(Γjh), local quasi-uniformity gives us an associated mesh size h^k that satisfies 1cQ≤h^khl≤cQ,∀lsuch thatτl∩supp(φk)≠0,k=1,…,M,

where hl is the mesh size of τl. Then, we require that 33 577-∑k1∈J(l)h^k1∑k2∈J(l)h^k2-1≥c0j>0∀τl∈Γjh,

with a global constant c0j>0.

Remark 6

The local mesh condition (33) was first introduced in [30, Chapter 2.1]. It is considered mild because it is fulfilled by a broad set of meshes used in applications, including geometrically graded meshes, algebraically 2-graded meshes and families of meshes generated by adaptive red-green algorithms [18].

Remark 7

It is worth noticing that the local mesh condition (33) also guarantees that the L2(Γ)-duality product between Xh(Γ) and Yh(Γ) as chosen in (29) and (30) is stable [30]. Hence, by choosing m to be ≪·,·≫, our implementation leads to a Galerkin matrix Mh that is bounded and invertible.

Under these considerations, we will show the following bound for the resulting condition number:

Theorem 1

Let Wh be the Galerkin matrix defined in (31), and Mh the Galerkin matrix of the duality pairing for Xh(Γ)×Yh(Γ) as chosen above.

Assume that there exists an operator Rh+:H+12(Γ)→Xh(Γ) such thatRh+ is a h-uniformly bounded projection

Rh+(H+12([Γ]))⊆Xh(Γ).

Then, under the mesh conditions from Assumption 2, we have34 κ~sp(Mh-1BhVMh-TWh)≤(1+|logh|)2‖Rh+‖αM2‖W0‖‖V0‖‖M‖2αBαa,

where κ~sp denotes the spectral condition number away from the kernel, ‖·‖ the operator norms and α(·) the corresponding inf-sup constants.

The proof will be given at the end of this Section, after we have presented some required preliminary results.

Remark 8

We refer the reader to [1, Theorem 2.1] for a construction of the operator Rh+ for the case stated in Theorem 1. It is worth mentioning that here we still state it as a requirement to make it explicit that this is a key piece in our theory. Although we believe one can follow the cue from [1] to also build the projection operator required for V0, the reader should be cautioned that story becomes considerably more complicated when considering Maxwell equations.

Remark 9

The existence of the operator Rh+ required by Theorem 1 is proven in [1] for meshes that agree on the front and back of the structure. In practice this does not pose a major limitation because in a typical usage scenario the simple screens Γj are built by fusing together two or more unique meshes for the interfaces ∂Ωi∪∂Ωj (see Sect. 4.2). The interface meshes are used both as front and back and thus necessarily agree.

Auxiliary Lemmas

In this section we prove some auxiliary results that hold for the multi-screens under consideration. We remark that for this we follow the cue from [9, Sect. 5.2] and use the properties of the associated volume-based spaces.

Lemma 3

Let Γ be a multi-screen as in Assumption 1. Then, there exist continuous embeddingsH+12(Γ)↪H12(Γ0)×⋯×H12(Γm),H-12(Γ)↪H-12(Γ0)×⋯×H-12(Γm).

Proof

We recall that H+12(Γ)=H1(R3\Γ¯)/H0,Γ1(R3) and note that H1(R3\Γ¯)⊂H1(R3\∪j=0m∂Ωj¯). This induces the injectionH+12(Γ)=H1(R3\Γ¯)/H0,Γ1(R3)↪H1(R3\∪j=0m∂Ωj¯)/H0,Γ1(R3).

Additionally, we have the natural identificationH1(R3\∪j=0m∂Ωj¯)≅H1(Ω0)×⋯H1(Ωm),

that associates u∈H1(R3\∪j=0m∂Ωj¯) with (u|Ω0,⋯,u|Ωm). From this natural identification, we get the isomorphismH1(R3\∪j=0m∂Ωj¯)/H0,Γ1(R3)≅[H1(Ω0)/H0,Γ1(Ω0)]×⋯×[H1(Ωm)/H0,Γ1(Ωm)]≅H12(Γ0)×⋯×H12(Γm).

Therefore, we have the injectionH+12(Γ)↪H12(Γ0)×⋯×H12(Γm).

H-12(Γ)↪H-12(Γ0)×⋯×H-12(Γm) follows analogously. □

Lemma 4

Let Γ be a multi-screen as in Assumption 1. Then, for w∈H+12(Γ) and φ∈H-12(Γ) such that φ|Γj∈H~-12(Γj)∀j=0,⋯,m, we have that35 ≪w,φ≫=∑l=0m⟨w|Γl,φ|Γl⟩Sl,

Proof

Let us consider w∈H+12(Γ) and φ∈H-12(Γ). By definition≪w,φ≫=∫[Γ]wφdσ=∫R3\Γ¯p·∇U+Udiv(p)dx,

for U∈H1(R3\Γ¯) and p∈H(div,R3\Γ¯) such that πD(U)=w and πN(p)=φ.

For j=0,⋯,m, we set Uj=U|Ωj and pj=p|Ωj, and let nj denote the outwards unit normal vector to Ωj. Then, by linearity of the integrals and Green’s formula, we get36 ∫R3\Γ¯p·∇U+Udiv(p)dx=∑j=0m∫Ωjpj·∇Uj+Ujdiv(pj)dx,=∑j=0m∫∂Ωj(Uj)|∂Ωjnj·(pj)|∂Ωjdσ.

Now, let us point out that for any j=0,⋯,m we know that functions in H1(R3\Γ¯) and p∈H(div,R3\Γ¯) do not jump across ∂Ωj\Γ. This allow us to simplify (36) further as37 ∑j=0m∫∂Ωj(Uj)|∂Ωjnj·(pj)|∂Ωjdσ=∫Γ∑j=0mvjμjdσ,

where vj=(Uj)|Γ and μj=nj·(pj)|Γ.

Next, we note thatμj=nj·(pj)|∂Ωj=πN(p|Ωj)=(πN(p))|Γj=φ|Γj∈H~-12(Γj),vj=(Uj)|∂Ωj=πD(U|Ωj)=(πD(U))|Γj=w|Γj∈H12(Γj).

This implies that on the right hand side of (37), we are allowed to split the integral over Γ into the sum of the integrals over Γj for j=0,⋯,m. Furthermore, we can write it in terms of the H12(Γj)×H~-12(Γj)-duality pairings, i.e.∫Γ∑j=0mvjμjdσ=∑j=0m∫Γjvjμjdσ=∑j=0m⟨vj,μj⟩Γj.

Finally, using again that μj=φ|Γj and vj=w|Γj, we conclude that≪w,φ≫Γ=∑j=0m⟨w|Γj,φ|Γj⟩Γj.

□

The next ingredient we need is an inverse inequality on H+12(Γ). Before we can derive it, we recall an inverse inequality on standard trace spaces. Let us consider the simple screens Γj for j=0,⋯,m and let S0,-1(Γjh) be the space of piecewise constants on the mesh Γjh of Γj.

Lemma 5

([21, Lemma 2.8], [34, Sect. 5], [33, Chap. 4]) For j=0,⋯,m, the following inverse inequality holds:38 ‖φh‖H~-12(Γj)≤c2(1+|logh|)‖φh‖H-12(Γj),

for all φh∈S0,-1(Γjh)⊂H~-12(Γj), with mesh size h≤1 and c2>0 independent of h.

Remark 10

Lemma 5 also holds for all φh∈S0,-1(Γˇjh)⊂H~-12(Γj). In order to see this, one should note that: (i) all steps in the proof are valid for barycentric refinements; (ii) elements of S0,-1(Γˇjh) are linear combinations of piecewise constants on the barycentric refinement; (iii) our mesh assumptions guarantee that the number of neighbours is uniformly bounded with respect to the mesh size h. Hence, the coefficients of the linear combination can be absorbed by the constant c2, which remains independent of h.

Let S1,0(Γjh) be the space of piecewise linear functions on the primal mesh of Γj, as defined in [4]. We introduce the generalised L2-projection Q~hj:L2(Γj)→S1,0(Γjh) as39 ⟨Q~hju,ϕh⟩Γj=⟨u,ϕh⟩Γj,∀ϕh∈S0,-1(Γˇjh).

Then, from [23, Theorem 4.3], we know that under Assumption 2, we have40 ‖Q~hju‖H12(Γj)≤cQj‖u‖H12(Γj),∀u∈H12(Γj).

Moreover, [30, Thm. 2.1 and 2.2] shows that Assumption 2 implies the discrete inf-sup condition of the duality pairing of multi-trace spaces.

Next, note that Xh(Γ) and Yh(Γ) defined in Sect. 4.2 generalize to441 Xh(Γ)=⨂j=0mS01,0(Γjh),andYh(Γ)=⨂j=0mS0,-1(Γˇjh).

Lemma 6

(Inverse inequality in H-12(Γ)) Let Γ be a multi-screen as in Assumption 1 and consider finite dimensional spaces Xh(Γ)⊂H+12(Γ) and Yh(Γ)⊂H-12(Γ) defined in (41) and such that Assumption 2 is satisfied. Then, we have that for all φh∈Yh(Γ)42 ‖φh‖H-12(Γ)≤CP(1+|logh|)∑j=0m‖φh‖H-12(Γj)21/2,

with mesh size h≤1, and CP>0 independent of h.

Proof

By definition of dual norm and Lemma 4, we get43 ‖φh‖H-12(Γ)=supv∈H+12(Γ)\{0}|≪v,φh≫|‖v‖H+12(Γ)=supv∈H+12(Γ)\{0}|∑j=0m⟨vj,φhj⟩Γj|‖v‖H+12(Γ),

where we have set vj=v|Γj and φhj=(φh)|Γj. Then, let us consider the index set J of all the indices 0≤k≤m such that vk≠0. Thus, we have that44 |∑j=0m⟨vj,φhj⟩Γj|‖v‖H+12(Γ)=|∑j∈J⟨vj,φhj⟩Γj|‖v‖H+12(Γ).

Next, we derive two inequalities that will help us proceed. First, using the embedding from Lemma 3 and then Young’s inequality m-times, we obtain45 ‖v‖H+12(Γ)2≥∑j∈J‖vj‖H12(Γj)2≥1|J|∑j∈J‖vj‖H12(Γj)2

where |J| is the size of the index set J. Second, we remark that for a sum ∑i=1∗ai with all coefficients ai>0 and m∗∈N, we have that 1∑i=1∗ai≤1ak for all k=1,⋯m∗. Hence,46 m∗∑i=1∗ai≤∑k=1m∗1ak

Plugging these in (43) gives47 ‖φh‖H-12(Γ)=(44)supv∈H+12(Γ)\{0}|∑j∈J⟨vj,φhj⟩Γj|‖v‖H+12(Γ)≤(45)|J|1/2supv∈H+12(Γ)\{0}|∑j∈J⟨vj,φhj⟩Γj|∑k∈J‖vk‖H12(Sk)≤(46)1|J|1/2supv∈H+12(Γ)\{0}|∑j∈J⟨vj,φhj⟩Γj‖vj‖H12(Γj)|≤1|J|1/2∑j∈Jsupwj∈H12(Γj)\{0}|⟨wj,φhj⟩Γj|‖wj‖H12(Γj)=1|J|1/2∑j=0msupwj∈H12(Γj)\{0}|⟨wj,φhj⟩Γj|‖wj‖H12(Γj).

Then, using Q~hj from (39) and its continuity (40), we get48 ‖φh‖H-12(Γ)≤1|J|1/2∑j=0mcQjsupwj∈H12(Γj)\{0}|⟨Q~hjwj,φhj⟩Γj|‖Q~hjwj‖H12(Γj)≤1|J|1/2∑j=0mcQjsupwhj∈S0,1(Γjh)\{0}|⟨whj,φhj⟩Γj|‖whj‖H12(Γj).

Applying Cauchy-Schwarz and (38), we obtain49 ‖φh‖H-12(Γ)≤1|J|1/2∑j=0mcQj‖φhj‖H~-12(Γj)≤c3(1+|logh|)∑j=0m‖φhj‖H-12(Γj),

with c3:=c2max(cQj)|J|1/2. For convenience, we work the estimate further50 ‖φh‖H-12(Γ)2≤c32(1+|logh|)2∑j=0m‖φhj‖H-12(Γj)2≤(45)|J|c32(1+|logh|)2∑j=0m‖φhj‖H-12(Γj)2,

which gives (42) with CP=c2max(cQj).

□

Finally, recall that‖u‖H~+12([Γ])=infx∈H+12([Γ])‖u+x‖H+12(Γ),

and that S1,0(Th) is the space spanned by piecewise linear “continuous” functions on (the virtual mesh) Th. We define the norm51 ‖u‖H~+12([Γ]),D:=infxh∈H+12([Γ])∩S1,0(Th)‖u+xh‖H+12(Γ),

and study its relation with the continuous jump norm.

Lemma 7

Assume that there exists an operator Rh+:H+12(Γ)→S1,0(Th) such that (i) Rh+ is a h-uniformly bounded projection,

(ii) Rh+(H+12([Γ]))⊆H+12([Γ])∩S1,0(Th).

Then there exist constants ce,Ce>0 independent of h and such thatce‖uh‖H~+12([Γ])≤‖uh‖H~+12([Γ]),D≤Ce‖uh‖H~+12([Γ]),∀uh∈H~+12([Γ])∩S1,0(Th).

Proof

By definition, for uh∈H~+12([Γ])∩S1,0(Th), we have that52 ‖uh‖H~+12([Γ])≤‖uh‖H~+12([Γ]),D

Choose uh=Rh+u where u∈H~+12([Γ])⊂H+12(Γ). Then we get‖uh‖H~+12([Γ]),D=infxh∈H+12([Γ])∩S1,0(Th)‖uh+xh‖H+12(Γ).

Note that, by surjectivity of Rh+ and property (ii), there exists an x∈H+12([Γ]) such that xh=Rh+x. This allows us to write‖uh‖H~+12([Γ]),D=infx∈H+12([Γ])‖Rh+(u+x)‖H+12(Γ).

Since Rh+ is continuous, we further obtain53 ‖uh‖H~+12([Γ]),D≤‖Rh+‖infx∈H+12([Γ])‖u+x‖H+12(Γ)=‖Rh+‖‖uh‖H~+12([Γ])

for all uh∈H~+12([Γ])∩S1,0(Th).

From this, we conclude that the two norms are equivalent. Moreover, the fact that Rh+ is h-uniformly bounded guarantees that the related constants are h-independent.

□

Proof of Theorem 1

We are now in the position to show the ellipticity of our preconditioner and then prove Theorem 1.

Proposition 4

Let BhV:Yh(Γ)→(Yh(Γ))′ be the linear operator corresponding to BhV defined in (58). For the discrete spaces defined in this section, we have that for all h∈R it holds that54 ≪BhVuh,uh≫≥αBV(1+|logh|)-2‖uh‖H-12(Γ)2

for all uh∈H-12(Γ), and with αBV>0 independent of h.

Proof

By definition of BhV we have that55 ≪BhVuh,uh≫=∑l⟨V0uh,uh⟩Il

Let us write uhl:=uh|Il. Recall that V0 is elliptic on each Il and hence⟨V0uhl,uhl⟩Il≥c~l‖uhl‖H~-12(Il)2≥c~l‖uhl‖H-12(Il)2

is satisfied for all h, and with c~l>0 independent of h. Therefore, setting c∗=minc~l, we get56 ≪BhVuh,uh≫≥∑lc~l‖uhl‖H-12(Il)2≥c∗∑l‖uhl‖H-12(Il)2.

Finally, using the inverse inequality from Lemma 6, we obtain the desired result with αBV:=c∗C~N2. □

Proof of Theorem 1

Given the ellipticity constants from Propositions 2 and 4, our bilinear forms satisfy their inf-sup constants. From Lemma 7, norm equivalence holds with ‖uh‖X/Xh≤‖Rh+‖‖uh‖X/X for all uh in Xh. This together with Lemma 2 gives that the result follows from the bound derived in Sect. 3. □

Calderón preconditioning on reduced quotient-space BEM

The above analysis has been carried out for the case where the BE space is chosen to approximate the entire multi-trace space. Since the solution is determined in the jump space, it can be worth while to investigate whether (combinations of) basis functions can be deleted and whether the resulting method remains amenable to operator preconditioning schemes.

In this section we will introduce several ways in which the number of basis functions can be reduced, what mileage can be expected from the resulting methods, and we discuss what the ramifications are for implementations in code of these methods.

The most straightforward approach to building a well-defined BEM formulation on multi-screens is to introduce a BE space for the multi-trace space that is contained in the product space. Two key ingredients for the success of this approach are thatthe discrete left/right nullspace Xh(Γ) equals H+12([Γ])∩Xh(Γ); and that

the quotient Xh(Γ)/Xh(Γ) approximates H~+12([Γ]).

However, because we are interested in finding an approximate solution in the quotient space H~+12([Γ]), we are free to consider BE spaces Xh(Γ)∘ that do not approximate all of H+12(Γ) as long as the corresponding discrete nullspace Xh(Γ)∘=H+12([Γ])∩Xh(Γ)∘ is still a subset of H+12([Γ])∩Xh(Γ) and the quotient spaces Xh(Γ)∘/Xh(Γ)∘ and Xh(Γ)/Xh(Γ) are equal. Similar choices can be made to select a reduced dual BE space Yh(Γ)∘⊆Yh(Γ). The quality of the resulting operator preconditioning depends on the stability of the restriction of the duality form to this subspace.

How does this work in practice? In the case of nodal elements in S01,0(Th), any given basis function relates to a function in Xh(Γ) by completing it with its counterpart(s) on the opposite side(s) of the multi-screen. By removing one of the basis functions from the standard nodal basis for Xh(Γ)∘, the dimension of the discrete nullspace Xh(Γ)∘ decreases by one. The dimension of the complement of Xh(Γ)∘ remains unchanged and so necessarily Xh(Γ)∘/Xh(Γ)∘=Xh(Γ)/Xh(Γ). To put it in more physical terminology: the reduced discrete multi-trace space Xh(Γ)∘ radiates the same fields as the original one.

There are a number of reduction schemes that are fairly straightforward to implement. We will discuss three strategies that can be applied to a multi-screen Γ comprising a single junction where an odd number m simple screens meet. (i) Partial reduction: In the partition [Γ]=∪i=1mIi, basis functions based on the terms i=3,5,7,... can be discarded. This is extremely easy to implement and boils down to using Xh(Γ)∘=S01,0(I1,h)×⨂i=1⌊m/2⌋S01,0(I2i,h) instead of Xh(Γ)=⨂i=1mS01,0(Ii,h).

(ii) Single strip: The partial reduction described above in essence removes the back from part of [Γ]. This still leaves significant redundancy in Xh(Γ)∘. In our example leaving out S01,0(Ii,h) for i=3,5,7,... still leaves all the basis functions on Γ2 (excluding basis functions on the junction) that can be completed by basis functions on the other side to yield functions in Xh(Γ). As a result, we can further discard basis functions in S01,0(I1,h) that lie in the interior of Γ2. If the implementer has access to node-triangle adjacency information this approach requires minimal coding effort. The resulting BE space is Xh(Γ)∘=S01,0(I1,h)∘×⨂i=1⌊m/2⌋S01,0(I2i,h), where the ∘ superscript on the first factor denotes that this BE space is reduced by leaving out redundant basis functions linked to nodes that are in I1,h∩Γ2 but not on the junction.

(iii) Fixed overlap: Note that the efficiency of the resulting preconditioning method depends on the lower bound for the duality form ≪.,.≫ on Xh(Γ)∘×Yh(Γ)∘, which may depend on the geometry and hence may indirectly depend on h when using a single strip reduction. In those situations where this is undesirable, one can opt to retain not only those basis functions in S01,0(I1,h) that are positioned on Γ1 or on the junction, but also those on Γ2 inside a strip within a fixed mesh independent distance from the junction. Likely, this will require manipulations to the code at the level of mesh generation. The resulting method will lead to an increasing redundancy in basis functions as h tends to zero, but the user is guaranteed that the preconditioner efficiency will not be limited by degradation of the stability of the duality pairing.

We illustrate these three reduction strategies for m=3 on Fig. 7. We also point out that when m is even, all these reductions are also valid, but since one can always find a partial reduction that provides a minimal representation of the quotient space, i.e. Xh(Γ)∘=⨂i=1m/2S01,0(I2i,h), the other two proposed strategies are not computationally attractive.

Finally, it is worth mentioning that regardless the choice of reduction method and the corresponding primal BE space Xh(Γ)∘, the construction of the dual BE space Yh(Γ)∘ remains the same. The construction goes along the lines of what is described in [4], starting from the reduced surfaces and corresponding meshes57 Ii∩⋃u∈Xh(Γ)∘suppu

This means in particular that for the Neumann problem, the dual space simply consists of piecewise constants on the dual cells as detailed in [3].

Moreover, since the discrete stability of the duality pairing also implies the continuity estimate (40) used in our proofs, we have that by ensuring this stability, all results in Sect. 4 can be extended to the proposed reduced Quotient-space BEM and hence we are still within the framework of Theorem 1. For this, it is crucial to identify the reduced primal space on I1,h with the (full) space on a truncated simple screen I1,h∘, which are displayed in blue in Fig. 7c and d. So for example, one identifies S01,0(I1,h)∘ with S01,0(I1,h∘).Fig. 7 Meshes illustrating full multi-trace discretisation and three reduction strategies used here

Numerical results

Preconditioning the hypersingular operator (Neumann problem)

Consider the geometry in Fig. 7. The structure is submersed in an externally generated potential uinc(x)=x1+x3.

Linear systems are solved by the preconditioned conjugate gradient (PCG) method with the relative tolerance set to 2.0e-5. To build the preconditioners, application of the inverse Gram matrix is required. This action is computed by running a second, inner GMRES solver (the Gram matrices are between different BE spaces and thus not self-adjoint) within the outer, primal solver. Numerical experiments have shown that it is important to set the tolerance for this inner GMRES sufficiently low. In the experiments presented here the tolerance is set to 2.0e-12. Fortunately the Gram matrices are well conditioned and application of their inverses through GMRES can be computed in a small and linear number of operations, even for these very small tolerances.

Another important aspect of implementing the preconditioning strategies presented above is the use of high quality quadrature rules, especially for interactions between geometric elements that are close together. Specifically, it is important that left and right nullspaces of the discrete bilinear forms are invariant upon introduction of the quadrature error. One can either choose to adopt highly accurate quadrature rules or to use rules that are symmetric with respect to back-front mirroring across the multi-screen. Here we have opted for the highly accurate and kernel independent Sauter-Schwab rules described in [29, Chapter 5].Fig. 8 PCG iterations vs h for the Neumann problem at Γ as in Fig. 7

The obtained results are displayed in Fig. 8. For all reductions of Xh(Γ) depicted in Fig. 7, there is a clear improvement. After preconditioning there remains a slow increase in the number of iterations, commensurate with the logarithmic grow in the (34).

Preconditioning the weakly singular operator (Dirichlet problem)

We consider the same geometry, excitation and PCG tolerance used to study our preconditioner for the Neumann problem. An important difference is that basis functions for the Dirichlet problem are linked to triangles of the mesh, as opposed to vertices. The support of the primal basis functions spans only a single triangle. The reduction of the multi-trace space can be done up to the point where there is no overlap between the simple screens that support the reduced finite element spaces. For ease of comparison, we still refer to this reduction strategy as single strip.

To arrive at a non-singular preconditioner, the hypersingular operator is regularised as in [31], resulting in a block diagonal preconditioner based on58 BhW:=Wˇh,1r000Wˇh,2r000Wˇh,3r,

where Wˇh,lr[i,j]=⟨W0ψˇj,ψˇi⟩Il+⟨ψˇj,1l⟩Il⟨1l,ψˇi⟩Il with ψˇi,ψˇj in the standard basis of S1,0(Iˇl,h) for l=1,2,3 and 1l the constant function on Il taking on the value 1.Fig. 9 PCG iterations vs h for the Dirichlet problem at Γ as in Fig. 7

Essentially all conclusions drawn for the Neumann problem carry over to the study of the numerical solution of the Dirichlet problem (Fig. 9): for all reduction strategies, application of the block diagonal Calderón preconditioner leads to a much smaller number of iterations.

Application to multi-screens of Type B

Fig. 10 Two possible coverings for [Γ]. The most economic covering on the left precludes the definition of BE spaces of direct product type (left). Allowing part of [Γ] to be multiply covered resolves this problem (right)

For multi-screens of Type B, slight modifications to the choice of BE spaces are required in order to arrive at a linear system requiring only few iterations for its solution. It may seem most natural to write [Γ] as the union of the simple screens depicted on the left in Fig. 10. Unfortunately, this partitioning does not allow the construction of a BE space Xh(Γ) that can be written as the product of BE spaces supported by the Ii. The issue is that the solution for the Neumann problem in general will not be in ⨂iH~12(Ii) and, as a result, basis functions along the segment from (0, 0, 0) to (0,-0.5,0) cannot be discarded.

Allowing overlapping coverings of [Γ] as depicted on the right of Fig. 10 solves this problem. For l=1,..,m, let I¯l denote either Il with overlap when needed, or without overlap. We can use the BE space ⨂i=1mS1,0(I¯i,h), which in a sense is larger than what we need but still leads to the correct quotient space.

We use this approach to solve the Neumann problem for the geometry in Fig. 10 and for the excitation uinc(x)=x1+x3. Figure 11 shows that upon preconditioning the number of iterations is much lower than what is required to solve the original linear system when solving the Neumann problem.

Figure 12 shows the results for the Dirichlet problem. They are in line with those from Fig. 9: the higher offset in the iteration count results in a cross-over point at smaller values of h, but asymptotically the preconditioner leads to a more efficient algorithm.

The numerical results presented in this section have been produced with the boundary element package BEAST.jl.5 The scripts to reproduce them can be found in a public Github repository.6

Note that this strategy, where we allow part of [Γ] to be multiply covered, can also be applied to enable the modelling of potential problems near Type C multi-screens such as the Möbius strip.Fig. 11 PCG iterations vs h for the Neumann problem for scattering by a geometry of type B

Fig. 12 PCG iterations vs h for the Dirichlet problem for scattering by a geometry of type B

Supplement: results for the Helmholtz equation

The above analysis does not hold verbatim for BIEs for the solution of Helmholtz BVPs: the possibility of resonances downgrades ellipticity estimates for the BIOs to mere Gårding inequalities and the number of required GMRES iterations is only loosely connected to the spectral condition number. Notwithstanding the method described in this paper performs quite well for low and moderate frequencies. In this section a set of representative numerical results are presented.

Consider the geometry in Fig. 7. The structure is illuminated by an externally generated wave uinc(x)=expiκx3, where κ is the wave number.

Results are displayed in Tables 1, 2, 3 and 4. For each reduction scheme, the columns ’unprec’ refer to the solution without preconditioner and the columns ’prec’ to the solution with preconditioner. In particular, Tables 1 and  3 summarise the GMRES iteration count for the Neumann and Dirichlet problem, respectively, at κ=1.0. Results are consistent with those for the Laplace problem (κ=0).

Tables 2 and 4 presents the GMRES iteration count for the Neumann and Dirichlet problem, respectively, at κ=10.0. Even though asymptotically preconditioning leads to a more efficient method, benefits in practice show up at much smaller values for h. For the Dirichlet problem in conjunction with the single strip reduction strategy, cross-over was not recorded within the range for h explored here.Table 1 GMRES iterations at κ=1.0, Neumann problem

h	Full	Partial	Overlap	Strip	
	Unprec	Prec	Unprec	Prec	Unprec	Prec	Unprec	Prec	
0.1	8	5	10	7	14	11	9	9	
0.05	10	6	13	8	21	11	12	10	
0.025	15	6	18	9	29	11	16	11	
0.0125	21	7	26	9	41	12	21	12	

Table 2 GMRES iterations at κ=10.0, Neumann problem

h	Full	Partial	Overlap	Strip	
	Unprec	Prec	Unprec	Prec	Unprec	Prec	Unprec	Prec	
0.1	14	8	19	10	26	17	15	12	
0.05	17	8	23	10	36	17	18	13	
0.025	22	8	28	10	48	17	23	14	
0.0125	28	9	37	12	66	18	29	16	

Table 3 GMRES iterations at κ=1.0, Dirichlet problem

h	Full	Partial	Overlap	Strip	
	Unprec	Prec	Unprec	Prec	Unprec	Prec	Unprec	Prec	
0.1	17	120	19	18	17	17	17	21	
0.05	21	15	24	19	23	19	21	23	
0.025	27	16	30	20	29	20	27	26	
0.0125	34	17	37	22	35	21	34	30	

Table 4 GMRES iterations at κ=10.0, Dirichlet problem

h	Full	Partial	Overlap	Strip	
	Unprec	Prec	Unprec	Prec	Unprec	Prec	Unprec	Prec	
0.1	27	35	31	41	29	42	27	39	
0.05	33	37	36	44	35	44	33	43	
0.025	40	39	44	45	42	46	40	46	
0.0125	48	39	53	47	50	46	48	49	

Conclusions

We have presented an effective Calderón-type preconditioner for potential problems in the exterior of multi-screens that builds on quotient-space BEM and operator preconditioning. Moreover, we have proved and confirmed numerically that it performs as standard Calderón preconditioning does on simple screens.

From a computational point of view, quotient-space BEM considering the full discretisation of multi-valued traces has the advantage of requiring minimal geometrical information but the disadvantage of doubling the number of basis functions. As an alternative, we proposed different strategies to work with reduced multi-trace discretisations that use less basis functions but require more adaptations when using a standard BEM code. We gave details regarding the additional data requirements in the implementation of all these strategies, and used the developed framework to identify the requirements that reduced spaces need to meet in order to still have efficient Calderón-type preconditioning.

Finally, we briefly presented a heuristic strategy to precondition multi-screens that also appear in applications but that are not covered by our theory. Although in essence our approach follows the same principles of our Calderón-type preconditioner for type A multi-screens, rigorous analysis has been elusive and therefore has not been treated in this article. Indeed, the key missing piece is an extension of Lemma 6 for this case. Nevertheless, we offered numerical experiments to investigate its effectivity.

Current and future work also involves extending the analysis of this preconditioning approach to Maxwell equations, where numerical results are promising [15], yet the construction of the jump aware projection Rh is considerably more challenging.

Acknowledgements

The authors would like to thank the referees whose comments helped to significantly improve this manuscript.

Funding

This project has received funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (Grant agreement No. 101001847) and from the Dutch Research Council (NWO) under the NWO-Talent Programme Veni with the project number VI.Veni.212.253.

Declarations

Conflict of interest

The research leading to these results received funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (Grant agreement No. 101001847) and from the Dutch Research Council (NWO) under the NWO-Talent Programme Veni with the project number VI.Veni.212.253.

1 We remind the reader that H(div,Ω):={v∈L2(Ω)3:divv∈L2(Ω)3} for Ω=R3 or Ω=R3\Γ¯, and refer to [19, Sect. 1.1] for definitions of the relevant Sobolev spaces.

2 Indeed, this definition of Neumann trace generalises the fact that for a Lipschitz domain Ω⊂R3, the normal derivative of v∈H1(Δ,Ω) can be computed in H-12(Ω) by projecting ∇v∈H(div,Ω) to the normal vector field n of ∂Ω.

3 Expressions ∫[Γ]⋯dσ should not be read as integrals with respect to the Lebesgue measure on [Γ]. This is simply notation introduced in [9] and that we kept for later convenience.

4 Here we follow the numbering given by the underlying Lipschitz partition, while in Sect. 4.2 we let the user number Il differently. This is harmless because one can always find a Lipschitz partition such that there is a one-to-one correspondence between Il and Γj.

5 https://github.com/krcools/BEAST.jl.

6 https://github.com/krcools/Junctions_KC_CUT.jl.

Publisher's Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
==== Refs
References

1. Averseng, M.: A stable and jump-aware quasi-interpolant onto a discrete multi-trace space (arXiv:2211.08223 [math.NA]) (2022). ArXiv:2211.08223 (2012)
2. Averseng, M., Claeys, X., Hiptmair, R.: Boundary Element Methods for the Laplace Hypersingular Integral Equation on Multiscreens: a two-level Substructuring Preconditioner. arXiv e-prints arXiv:2310.09204 (2023). 10.48550/arXiv.2310.09204
3. Buffa A Christiansen SH The electric field integral equation on Lipschitz screens: definitions and numerical approximation Numer. Math. 2003 94 2 229 267 10.1007/s00211-002-0422-0
Buffa, A., Christiansen, S.H.: The electric field integral equation on Lipschitz screens: definitions and numerical approximation. Numer. Math. 94(2), 229–267 (2003). 10.1007/s00211-002-0422-010.1007/s00211-002-0422-0
4. Buffa A Christiansen SH A dual finite element complex on the barycentric refinement Math. Comput. 2007 76 260 1743 1769 10.1090/S0025-5718-07-01965-5
Buffa, A., Christiansen, S.H.: A dual finite element complex on the barycentric refinement. Math. Comput. 76(260), 1743–1769 (2007)10.1090/S0025-5718-07-01965-5
5. Carr M Topsakal E Volakis J A procedure for modeling material junctions in 3-d surface integral equation approaches IEEE Trans. Antennas Propagat. 2004 52 5 1374 1378 10.1109/TAP.2004.827247
Carr, M., Topsakal, E., Volakis, J.: A procedure for modeling material junctions in 3-d surface integral equation approaches. IEEE Trans. Antennas Propagat. 52(5), 1374–1378 (2004)10.1109/TAP.2004.827247
6. Chandler-Wilde SN Hewett DP Moiola A Sobolev spaces on non-Lipschitz subsets of Rn with application to boundary integral equations on fractal screens Integral Equ. Operat. Theory 2017 87 2 179 224 10.1007/s00020-017-2342-5
Chandler-Wilde, S.N., Hewett, D.P., Moiola, A.: Sobolev spaces on non-Lipschitz subsets of with application to boundary integral equations on fractal screens. Integral Equ. Operat. Theory 87(2), 179–224 (2017). 10.1007/s00020-017-2342-510.1007/s00020-017-2342-5
7. Christiansen SH Nédélec JC Des préconditionneurs pour la résolution numérique des équations intégrales de frontière de l’acoustique Comptes Rendus de l’Académie des Sciences-Series I-Mathematics 2000 330 7 617 622
Christiansen, S.H., Nédélec, J.C.: Des préconditionneurs pour la résolution numérique des équations intégrales de frontière de l’acoustique. Comptes Rendus de l’Académie des Sciences-Series I-Mathematics 330(7), 617–622 (2000)
8. Claeys, X., Giacomel, L., Hiptmair, R., Urzúa-Torres, C.: Quotient-space boundary element methods for scattering at complex screens. BIT Numerical Mathematics pp. 1–29 (2021)
9. Claeys X Hiptmair R Integral equations on multi-screens Integral Equ. Operat. Theory 2013 77 2 167 197 10.1007/s00020-013-2085-x
Claeys, X., Hiptmair, R.: Integral equations on multi-screens. Integral Equ. Operat. Theory 77(2), 167–197 (2013). 10.1007/s00020-013-2085-x10.1007/s00020-013-2085-x
10. Claeys X Hiptmair R Integral equations for electromagnetic scattering at multi-screens Integral Equ. Operat. Theory 2016 84 1 33 68 10.1007/s00020-015-2242-5
Claeys, X., Hiptmair, R.: Integral equations for electromagnetic scattering at multi-screens. Integral Equ. Operat. Theory 84(1), 33–68 (2016). 10.1007/s00020-015-2242-510.1007/s00020-015-2242-5
11. Cools, K.: Mortar boundary elements for the efie applied to the analysis of scattering by pec junctions. In: 2012 Asia-Pacific Symposium on Electromagnetic Compatibility, pp. 165–168 (2012)
12. Cools, K., Andriulli, F.P.: Accuracy of the calderon preconditioned efie for the scattering by pec junctions. In: 2015 USNC-URSI Radio Science Meeting (Joint with AP-S Symposium), pp. 144–144 (2015)
13. Cools, K., Andriulli, F.P.: A regularised electric field integral equation for scattering by perfectly conducting junctions. In: 2015 9th European Conference on Antennas and Propagation (EuCAP), pp. 1–4 (2015)
14. Cools, K., Andriulli, F.P.: Well-conditioned saddle point description for scattering by a metallic junction. In: 2015 International Conference on Electromagnetics in Advanced Applications (ICEAA), pp. 1349–1352 (2015)
15. Cools, K., Urzúa-Torres, C.: Preconditioners for multi-screen scattering. In: 2022 International Conference on Electromagnetics in Advanced Applications (ICEAA), pp. 172–173 (2022)
16. Ervin VJ Stephan EP A boundary element Galerkin method for a hypersingular integral equation on open surfaces Math. Methods Appl. Sci. 1990 13 4 281 289 10.1002/mma.1670130402
Ervin, V.J., Stephan, E.P.: A boundary element Galerkin method for a hypersingular integral equation on open surfaces. Math. Methods Appl. Sci. 13(4), 281–289 (1990). 10.1002/mma.167013040210.1002/mma.1670130402
17. Ervin VJ Stephan EP El-Seoud SA An improved boundary element method for the charge density of a thin electrified plate in R3 Math. Methods Appl. Sci. 1990 13 4 291 303 10.1002/mma.1670130403
Ervin, V.J., Stephan, E.P., El-Seoud, S.A.: An improved boundary element method for the charge density of a thin electrified plate in . Math. Methods Appl. Sci. 13(4), 291–303 (1990). 10.1002/mma.167013040310.1002/mma.1670130403
18. Gimperlein H Stocek J Urzúa-Torres C Optimal operator preconditioning for pseudodifferential boundary problems Numer. Math. 2021 148 1 1 41 10.1007/s00211-021-01193-9
Gimperlein, H., Stocek, J., Urzúa-Torres, C.: Optimal operator preconditioning for pseudodifferential boundary problems. Numer. Math. 148(1), 1–41 (2021). 10.1007/s00211-021-01193-910.1007/s00211-021-01193-9
19. Girault V Raviart P Finite element methods for Navier-Stokes equations 1986 Berlin Springer
Girault, V., Raviart, P.: Finite element methods for Navier-Stokes equations. Springer, Berlin (1986)
20. Hayami, K.: Convergence of the conjugate gradient method on singular systems (2020)
21. Heuer, N.: Preconditioners for the p-version of the boundary element galerkin method in . Habilitation thesis, University of Hannover (1998). http://webdoc.sub.gwdg.de/ebook/dissts/Hannover/Heuer2000.pdf
22. Hiptmair R Operator preconditioning Comput. Math. Appl. 2006 52 5 699 706 10.1016/j.camwa.2006.10.008
Hiptmair, R.: Operator preconditioning. Comput. Math. Appl. 52(5), 699–706 (2006)10.1016/j.camwa.2006.10.008
23. Hiptmair, R., Jerez-Hanckes, C., Urzúa-Torres, C.: Optimal operator preconditioning for boundary elements on open curves. Tech. Rep. 2013-48, Seminar for Applied Mathematics, ETH Zürich, Switzerland (2013)
24. Hiptmair R Jerez-Hanckes C Urzúa-Torres C Optimal operator preconditioning for Galerkin boundary element methods on 3-dimensional screens SIAM J. Numer. Anal. 2020 58 1 834 857 10.1137/18M1196029
Hiptmair, R., Jerez-Hanckes, C., Urzúa-Torres, C.: Optimal operator preconditioning for Galerkin boundary element methods on 3-dimensional screens. SIAM J. Numer. Anal. 58(1), 834–857 (2020)10.1137/18M1196029
25. Hiptmair, R., Urzúa-Torres, C.: Dual Mesh Operator Preconditioning On 3D Screens: Low-Order Boundary Element Discretization. Tech. Rep. 2016-14, Seminar for Applied Mathematics, ETH Zürich, Switzerland (2016)
26. Kaasschieter E Preconditioned conjugate gradients for solving singular systems J. Comput. Appl. Math. 1988 24 1 265 275 10.1016/0377-0427(88)90358-5
Kaasschieter, E.: Preconditioned conjugate gradients for solving singular systems. J. Comput. Appl. Math. 24(1), 265–275 (1988). 10.1016/0377-0427(88)90358-510.1016/0377-0427(88)90358-5
27. McLean W Strongly Elliptic Systems and Boundary Integral Equations 2000 Cambridge, UK Cambridge University Press
McLean, W.: Strongly Elliptic Systems and Boundary Integral Equations. Cambridge University Press, Cambridge, UK (2000)
28. McLean W Steinbach O Boundary element preconditioners for a hypersingular integral equations on an interval Adv. Comp. Math. 1999 11 4 271 286 10.1023/A:1018944530343
McLean, W., Steinbach, O.: Boundary element preconditioners for a hypersingular integral equations on an interval. Adv. Comp. Math. 11(4), 271–286 (1999)10.1023/A:1018944530343
29. Sauter, S.A., Schwab, C.: Boundary element methods, Springer Series in Computational Mathematics, vol. 39. Springer-Verlag, Berlin, Berlin (2011). 10.1007/978-3-540-68093-2. Translated and expanded from the 2004 German original
30. Steinbach, O.: Stability estimates for hybrid coupled domain decomposition methods. Lecture Notes in Mathematics, vol. 1809. Springer-Verlag, Berlin (2003)
31. Steinbach O Wendland W The construction of some efficient preconditioners in the boundary element method Adv. Comput. Math 1998 9 191 216 10.1023/A:1018937506719
Steinbach, O., Wendland, W.: The construction of some efficient preconditioners in the boundary element method. Adv. Comput. Math 9, 191–216 (1998)10.1023/A:1018937506719
32. Stephan E Boundary integral equations for screen problems in R3 Integral Equ. Operat. Theory 1987 10 2 236 257 10.1007/BF01199079
Stephan, E.: Boundary integral equations for screen problems in . Integral Equ. Operat. Theory 10(2), 236–257 (1987)10.1007/BF01199079
33. Toselli, A., Widlund, O.: Domain decomposition methods—algorithms and theory, Springer Series in Computational Mathematics, vol. 34. Springer-Verlag, Berlin (2005). 10.1007/b137868
34. Xu J Zou J Some nonoverlapping domain decomposition methods SIAM Rev. 1998 40 4 857 914 10.1137/S0036144596306800
Xu, J., Zou, J.: Some nonoverlapping domain decomposition methods. SIAM Rev. 40(4), 857–914 (1998). 10.1137/S003614459630680010.1137/S0036144596306800
35. Yla-Oijala P Taskinen M Sarvas J Surface integral equation method for general composite metallic and dielectric structures with junctions Progr. Electromagn. Res. 2005 52 81 108 10.2528/PIER04071301
Yla-Oijala, P., Taskinen, M., Sarvas, J.: Surface integral equation method for general composite metallic and dielectric structures with junctions. Progr. Electromagn. Res. 52, 81–108 (2005)10.2528/PIER04071301
