Caccioppoli-type estimates and H -Matrix approximations to inverses for FEM-BEM couplings

We consider three diﬀerent methods for the coupling of the ﬁnite element method and the boundary element method, the Bielak-MacCamy coupling, the symmetric coupling, and the Johnson-Nédélec coupling. For each coupling we provide discrete interior regularity estimates. As a consequence, we are able to prove the existence of exponentially convergent H -matrix approximants to the inverse matrices corresponding to the lowest order Galerkin discretizations of the couplings.


Introduction
Transmission problems are usually posed on unbounded domains, where a (possibly nonlinear) equation is given on some bounded domain, and another linear equation is posed on the complement of the bounded domain.While the interior problem can be treated numerically by the finite element method (FEM), the unbounded nature of the exterior problem makes this problematic.A suitable method to treat unbounded problems is provided by the boundary element method (BEM), where the differential equation in the unbounded domain is reformulated via an integral equation posed just on the boundary.In order to combine both methods for transmission problems, additional conditions on the interface have to be fulfilled, which leads to different approaches for the coupling of the FEM and the BEM.We study three different FEM-BEM couplings, the Bielak-MacCamy coupling [BM84], Costabel's symmetric coupling [Cos88,CES90], and the Johnson-Nédélec coupling [JN80].Well-posedness and unique solvability of these formulation have been studied in, e.g., [Ste11, Say13, AFF + 13], where a main observation is that the couplings are equivalent to an elliptic problem.
Elliptic problems typically feature interior regularity known as Caccioppoli estimates, where stronger norms can be estimated by weaker norms on larger domains.In this paper, we provide such Caccioppolitype estimates for the discrete problem.More precisely, we obtain simultaneous interior regularity estimates for the finite element solution as well as for the single-and double-layer potential of the boundary element solution (cf.Theorems 2.3, 2.4, 2.5).Discrete Caccioppoli-type estimates for the FEM and the BEM separately can be found in our previous works [FMP15,AFM20,FMP16,FMP17]. While the techniques for the FEM and the BEM part are similar therein, some essential modifications have to be made to treat the coupling terms on the boundary.
An important consequence of Caccioppoli-type estimates is the existence of low-rank approximants to inverses of FEM or BEM matrices, as these inverses are usually dense matrices [BH03, Bör10a, FMP15, FMP16, FMP17].In particular, FEM and BEM inverses can be approximated in the data-sparse Hmatrix format, introduced in [Hac99].In comparison with other compression methods, H-matrices have the advantage that they come with an additional approximative arithmetic that allows for addition, multiplication, inversion or LU -decompositions in the H-matrix format; for more details we refer to [Gra01,GH03,Hac09].In this work, we present an approximation result for the inverses of stiffness matrices corresponding to the lowest order FEM-BEM discretizations.On admissible blocks, determined by standard admissibility conditions, the inverses can be approximated by a low-rank factorization, where the error converges exponentially in the rank employed.
The paper is structured as follows: In Chapter 2, we present our model problem and state the main results of the article, the discrete Caccioppoli-type interior regularity estimates for each coupling, and the existence of exponentially convergent H-matrix approximants to the inverse matrices corresponding to the FEM-BEM discretizations of the couplings.Chapter 3 is concerned with the proofs of the Caccioppolitype estimates.Chapter 4 provides an abstract framework for the proof of low-rank approximability to inverse matrices, which can be applied for other model problems as well.Finally, Chapter 5 provides some numerical examples.
Remark 2.1.In the following, we consider three different variational formulations, namely, the symmetric coupling, the Bielak-MacCamy coupling, and the Johnson-Nédélec coupling for our model problem.All three are well-posed without compatibility assumptions on the data.The compatibility condition f, 1 L 2 (Ω) + ϕ 0 , 1 L 2 (Γ) = 0 for d = 2 ensures the radiation condition (2.1e); lifting the compatibility condition yields a solution that satisfies a different radiation condition, namely, u ext = b log |x|+O(|x| −1 ) as |x| → ∞ for some b ∈ R for the three coupling strategies considered.Our analysis requires only the unique solvability of the variational formulations.
The single-layer operator V is elliptic for d = 3 and for d = 2 provided diam(Ω) < 1.The hyper-singular operator W is semi-elliptic with a kernel of dimension being the number of components of connectedness of Γ.
In addition to the boundary integral operators, we need the volume potentials V and K defined by In this paper, we study discretizations of weak solutions of the model problem reformulated via three different FEM-BEM couplings: the Bielak-MacCamy coupling, Costabel's symmetric coupling, and the Johnson-Nédélec coupling.All these couplings lead to a variational formulation of finding where a : X × X → R is a bilinear form and g : X → R is continuous linear functional.
For the discretization, we assume that Ω is triangulated by a quasi-uniform mesh T h = {T 1 , . . ., T n } of mesh width h := max Tj ∈T h diam(T j ).The elements T j ∈ T h are open triangles (d = 2) or tetrahedra (d = 3).Additionally, we assume that the mesh T h is regular in the sense of Ciarlet and γ-shape regular in the sense that we have diam(T j ) ≤ γ |T j | 1/2 for all T j ∈ T h , where |T j | denotes the Lebesgue measure of T j .By K h := {K 1 , . . ., K m }, we denote the restriction of T h to the boundary Γ, which is a regular and shape-regular triangulation of the boundary.
For simplicity, we consider lowest order Galerkin discretizations in S 1,1 (T h ) × S 0,0 (K h ), where with P p (T ) denoting the space of polynomials of maximal degree p on an element T .We let B h := {ξ j : j = 1, . . ., n} be the basis of S 1,1 (T h ) consisting of the standard hat functions, and we let W h := {χ j : j = 1, . . ., m} be the basis of S 0,0 (K h ) that consists of the characteristic functions of the surface elements.These bases feature the following norm equivalences: for the isomorphisms Φ : R n → S 1,1 (T h ), x → n j=1 x j ξ j and Ψ : R m → S 0,0 (K h ), y → m j=1 y j χ j .
Finally, we need the notion of concentric boxes.
Definition 2.2.(Concentric boxes) Two (quadratic) boxes B R and B R ′ of side length R and R ′ are said to be concentric if they have the same barycenter and B R can be obtained by a stretching of B R ′ by the factor R/R ′ taking their common barycenter as the origin.
Before we can state our first main results, the interior regularity estimates, we specify the norm we are working with, an h-weighted H (2.5) We note that u will be the interior solution, v be chosen as a single-layer potential and w as a double-layer potential (which jumps across Γ), which explains the different requirements for the set ω.
The following theorem is one of the main results of our paper.It states that for the interior finite element solution and the single-layer potential of the boundary element solution, a Caccioppoli type estimate holds, i.e., the stronger H 1 -seminorm can be estimated by a weaker h-weighted H 1 -norm on a larger domain.
16 , and let B R and B (1+ε)R be two concentric boxes.Assume that the data is localized away from Then, there exists a constant C depending only on Ω, d, and the γ-shape regularity of the quasi-uniform triangulation T h , such that for the solution (u h , ϕ h ) of (2.8) we have , where the norms on the right-hand side are defined in (2.5).
With the bases B h of S 1,1 (T h ) and W h of S 0,0 (K h ), the Galerkin discretization (2.8) leads to a block matrix A bmc ∈ R (n+m)×(n+m) where A ∈ R n×n is given by

Costabel's symmetric coupling
Using the representation formula, or more precisely, both single-and double-layer potential, for the exterior solution, one obtains an expression By coupling the interior and exterior solution in a symmetric way (which uses all four boundary integral operators), this leads to Costabel's symmetric coupling, introduced in [Cos88] and [Han90].Here, the bilinear form and right-hand side are given by (2.10a) (2.10b) The Galerkin discretization leads to the problem of finding With similar arguments as for the Bielak-MacCamy coupling, [AFF + 13] prove unique solvability for the symmetric coupling for any C ell > 0.
The following theorem is similar to Theorem 2.3 and provides a simultaneous Caccioppoli-type estimate for the interior solution as well as for the single-layer potential of the boundary solution and the double-layer potential of the trace of the interior solution.Here, the double-layer potential additionally appears since all boundary integral operators, especially the hyper-singular operator appear in the coupling.
Theorem 2.4.Let ε ∈ (0, 1) and R ∈ (0, 2 diam(Ω)) be such that h R < ε 32 , and let B R and B (1+ε)R be two concentric boxes.Assume that the data is localized away from B (1+ε)R , i.e., (supp f ∪ supp v 0 ∪ supp w 0 ) ∩ B (1+ε)R = ∅.Then, there exists a constant C depending only on Ω, d, and the γ-shape regularity of the quasi-uniform triangulation T h , such that for the solution (u h , ϕ h ) of (2.11) we have where the norm on the right-hand side is defined in (2.6).
With the bases B h of S 1,1 (T h ) and W h of S 0,0 (K h ), the Galerkin discretization (2.11) leads to a block matrix A sym ∈ R (n+m)×(n+m) where A, M, K are defined in (2.9), and W ∈ R n×n is given by W ij = W ξ j , ξ i L 2 (Γ) .

The Johnson-Nédélec coupling
The Johnson-Nédélec coupling, introduced in [JN80] again uses the representation formula for the exterior solution, but differs from the symmetric coupling in the way how the interior and exterior solutions are coupled on the boundary.Instead of all four boundary integral operators, only the single-layer and the double-layer operator are needed.The bilinear form for the Johnson-Nédélec coupling is given by (2.14b) The Galerkin discretization in S 1,1 (T h ) × S 0,0 (K h ) leads to the problem of finding As in the case of the Bielak-MacCamy coupling, the Johnson-Nédélec coupling has an unique solution provided C ell > 1/4, see [AFF + 13].
The following theorem gives the analogous result to Theorem 2.3 and Theorem 2.4 for the Johnson-Nédélec coupling.Similarly to the symmetric coupling, we simultaneously control a stronger norm of the interior solution and both layer potentials by a weaker norm on a larger domain.
With the bases B h of S 1,1 (T h ) and W h of S 0,0 (K h ), the Galerkin discretization (2.15) leads to a matrix A jn ∈ R (n+m)×(n+m) where A, M, K, V are defined in (2.9).

H-Matrix approximation of inverses
As a consequence of the Caccioppoli-type inequalities, we are able to prove the existence of H-matrix approximants to the inverses of the stiffness matrices corresponding to the discretized FEM-BEM couplings.
We briefly introduce the matrix compression format of H-matrices.For more detailed information, we refer to [Hac99,Beb08,Hac09,Bör10b].
The main idea of H-matrices is to store certain far field blocks of the matrix efficiently as a low-rank matrix.In order to choose blocks that are suitable for compression, we need to introduce the concept of admissibility.
Remark 2.7.Definition 2.6 clusters the degrees of freedom associated with triangulation T h of Ω and the triangulation K h of Γ simultaneously.
The block-partition of H-matrices is based on so-called cluster trees.
Definition 2.8 (cluster tree).A cluster tree with leaf size n leaf ∈ N is a binary tree T I with root I such that each cluster τ ∈ T I is either a leaf of the tree and satisfies |τ | ≤ n leaf , or there exist disjoint subsets τ ′ , τ ′′ ∈ T I of τ , so-called sons, with τ = τ ′ ∪ τ ′′ .Here and below, |τ | denotes the cardinality of the finite set τ .The level function level : T I → N 0 is inductively defined by level(I) = 0 and level(τ ′ ) := level(τ )+1 for τ ′ a son of τ .The depth of a cluster tree is depth(T I ) := max τ ∈TI level(τ ).Definition 2.9 (far field, near field, and sparsity constant).A partition P of I × I is said to be based on the cluster tree T I , if P ⊂ T I × T I .For such a partition P and a fixed admissibility parameter η > 0, we define the far field and the near field as P far := {(τ, σ) ∈ P : (τ, σ) is η-admissible}, P near := P \ P far . (2.18) The sparsity constant C sp of such a partition was introduced in [HK00, Gra01] as Definition 2.10 (H-matrices).Let P be a partition of I × I that is based on a cluster tree T I and η > 0. A matrix A ∈ R (n+m)×(n×m) is an H-matrix with blockwise rank r, if for every η-admissible block (τ, σ) ∈ P far , we have a low-rank factorization where Due to the low-rank structure on far-field blocks, the memory requirement to store an H matrix is given by ∼ C sp depth(T I )r(n + m).Provided C sp is bounded and the cluster tree is balanced, i.e., depth(T I ) ∼ log(n + m), which can be ensured by suitable clustering methods (e.g.geometric clustering, [Hac09]), we get a storage complexity of O(r(n + m) log(n + m)).
The following theorem shows that the inverse matrices A −1 bmc , A −1 sym , and A −1 jn corresponding to the three mentioned FEM-BEM couplings can be approximated in the H-matrix format, and the error converges exponentially in the maximal block rank employed.
Theorem 2.11.For a fixed admissibility parameter η > 0, let a partition P of I × I that is based on the cluster tree T I be given.Then, there exists an H-matrix B H with maximal blockwise rank r such that for the Bielak-MacCamy coupling.In the same way, there exists a blockwise rank-r H-matrix B H such that for the symmetric coupling and for the Johnson-Nédélec coupling.The constants C apx > 0 and b > 0 depend only on Ω, d, η, and the γ-shape regularity of the quasi-uniform triangulations T h and K h .

The Caccioppoli-type inequalities
In this section, we provide the proofs of the interior regularity estimates of Theorems 2.3-2.5.We start with some well-known facts about the volume potential operators V , K and the boundary integral operators V, K, K ′ , W .For details, we refer to [SS11, Ch. 3] and [Ste07,Ch. 6].
• With the interior trace operator γ int 0 (for Ω) and exterior trace operator γ ext 0 (for R d \Ω), we have • Similarly, with the interior γ int 1 u := γ int 0 ∇u • ν and exterior conormal derivative γ ext 1 u := γ ext 0 ∇u • ν (ν is the outward normal vector of Ω), we have and consequently the jump conditions • The potentials V ϕ and Ku are harmonic in R d \Γ and are bounded operators (see [SS11, Ch. 3.1.2]) Consequently, we have the boundedness for the boundary integral operators as In the following, the notation abbreviates ≤ up to a constant C > 0 which depends only on Ω, the dimension d, and the γ-shape regularity of T h .Moreover, we use ≃ to indicate that both estimates and hold.

The Bielak-MacCamy coupling
This section is dedicated to the proof of Theorem 2.3.The techniques employed are fairly similar to [FMP15,FMP16], where Caccioppoli-type estimates for FEM and BEM are proven.Nonetheless, in the case of the FEM-BEM couplings, the additional terms in the bilinear forms arising from the coupling on the boundary need to be treated carefully.
Then, there is C > 0 depending only on the shape-regularity of the triangulation and Γ such that for any discrete function Proof.For details, we refer to [FMP16].The main observation is that, on each element K ∈ K h , we have ∇ψ h | K ≡ 0. Therefore, the standard approximation result reduces to Since I Γ h is the L 2 -projection, we obtain an additional factor h 1/2 in the H −1/2 (Γ)-norm.A similar super-approximation result holds for the nodal interpolation operator In the proof of the Caccioppoli type inequality, we need the following inverse-type inequalities from [FMP16, Lem.3.8] and [FMP17, Lem.3.6]. . (3.9) Combining Lemma 3.1 with Lemma 3.2 (assuming supp η ⊂ B R ), we obtain estimates of the form .
(3.10) Remark 3.3.An inspection of the proof of (3.9) ([FMP17, Lem.3.6]) shows that the main observation is that Kv h is harmonic.The remaining arguments therein only use mapping properties and jump conditions for the potential K and can directly be modified such that the same result holds for the singlelayer potential as well, i.e., for every ψ h ∈ S 0,0 (T h ), we have . (3.11) Now, with the help of a local ellipticity result, the discrete variational formulation, and superapproximation, we are able to prove Theorem 2.3.
Proof of Theorem 2.3.In order to reduce unnecessary notation, we write (u, ϕ) for the Galerkin solution (u h , ϕ h ).The assumption on the support of the data implies the local orthogonality We note that this choice of δ implies that {K ∈ K h : supp η ∩ K = ∅} ⊂ B (1+δ/2)R .In the final step of the proof, we will choose two different values for δ (≤ ε) depending on ε -one of them, δ = ε 2 , explains the assumption made on ε in the theorem.
(See (3.23) for the precise form.)Since the ellipticity constant and we start with Young's inequality implies . (3.14) . (3.15) An elementary calculation shows . (3.16) Since the single-layer potential is harmonic in Ω, integration by parts (in Ω) and (3.17) Similarly, with integration by parts (in Ω and Ω ext ) and the jump condition of the single-layer potential we obtain . (3.18) Moreover, the symmetry of C implies . (3.20) Young's inequality and ∇η as well as Absorbing the gradient terms in (3.21)-(3.22) in the left-hand side of (3.20), we arrive at (3.23) Step 2: We apply the local orthogonality of (u, ϕ) to piecewise polynomials and use approximation properties.
Let I Ω h : C(Ω) → S 1,1 (T h ) be the nodal interpolation operator and I Γ h the L 2 (Γ)-orthogonal projection mapping onto S 0,0 (K h ).Then, the orthogonality (3.12) leads to (3.24) The term T 3 can be estimated in exactly the same way as in [FMP16].More explicitly, we need a second cut-off function Here, the support property of I Γ h (η 2 ϕ) follows from the assumption on δ.The trace inequality together with the super-approximation properties of I Γ h , expressed in (3.10), lead to . (3.25) With the same arguments, we obtain an estimate for T 4 (3.26) The volume term T 1 in (3.24) can be estimated as in [FMP15].Here, the super-approximation properties of I Ω h from (3.8), Young's inequality, and h δR ≤ 1 lead to (3.27) It remains to treat the coupling term T 2 involving the adjoint double-layer operator in (3.24).With the support property supp(I Ω h (η 2 u) − η 2 u) ⊂ B (1+δ/2)R , which follows from 8h ≤ δR, and (3.28) The multiplicative trace inequality for Ω, see, e.g., [BS02], the super-approximation property of I Ω h from (3.8), and h R ≤ δ 8 lead to (see also [FMP15, Eq. ( 25), ( 26)] for more details) We use estimate (3.11) and (3.29) in (3.28), which implies . (3.31) Step 3: We iterate (3.31) to obtain the claimed powers of h for the gradient terms.

The symmetric coupling
In this section, we provide the proof of Theorem 2.4.While some parts of the proof are similar to the proof of Theorem 2.3 and are therefore shortened, there are some differences as well, mainly that it does not suffices to study the single-layer potential.Indeed, one has to add a term containing the double-layer potential to the Caccioppoli inequality in order to get a localized ellipticity estimate.
Proof of Theorem 2.4.Again, we write (u, ϕ) for the Galerkin solution (u h , ϕ h ).The assumption on the support of the data implies the local orthogonality and will be chosen in the last step of the proof.
With the L 2 (Γ)-orthogonal projection I Γ h and the nodal interpolation operator I Ω h , the orthogonality (3.33) implies (3.39) The terms T 1 , T 3 , T 4 can be estimated with (3.27), (3.30) and (3.25) respectively as in the case for the Bielak-MacCamy coupling.Therefore, it remains to estimate T 2 and T 5 .We start with T 2 , which can be treated in the same way as in [FMP17].In fact, with techniques similar to (3.25), the proof of [FMP17, Lem.3.8] (taking v = Ku there and noting that [FMP17, Lemma 3.6] is employed, which does not impose orthogonality conditions on v) provides the estimate .
We finish the proof by estimating T 5 .To that end, we need another cut-off function 0 Ku, we get with a trace inequality and the approximation properties expressed in (3.10) that γ int 0 ( η Ku) . (3.40) Putting everything together in (3.39) and further in (3.38), and absorbing the terms in the left-hand side, finally yields . (3.41) Step 3: By reapplying (3.41) to the gradient terms with δ = ε 2 and suitable boxes, we get the desired result exactly as in step 3 of the proof of Theorem 2.3.

The Johnson-Nédélec coupling
In this section we prove the Caccioppoli-type inequality from Theorem 2.5 for the Johnson-Nédélec coupling.Most of the appearing terms have already been treated in the previous sections.The main difference is that the double-layer potential appears naturally due to the boundary coupling terms, but the local orthogonality is not suited to provide an approximation for it, since the hypersingular operator does not appear in the bilinear form.A remedy for this problem is to localize the double-layer potential by splitting it into a local near-field and a non-local, but smooth far-field.This techniques follows [FM18], where a similar localization using commutators is employed and a more detailed description of the method can be found.
Proof of Theorem 2.5.Once again, we write (u, ϕ) for the Galerkin solution (u h , ϕ h ).The assumption on the support of the data implies the local orthogonality We note that the condition η ≡ 1 on B (1+δ/4)R is additionally imposed due to following estimate (3.43), as the localization of the double-layer operator is additionally needed in comparison with the other couplings.
Step 1: We start with a localization of the double-layer potential.More precisely, with a second cut-off function η satisfying η ≡ 1 on B R and At first, we estimate the near-field v near := η K(ηu).The mapping properties of the double-layer potential, (3.5), together with the fact that supp ∇ η ⊂ B (1+δ/4)R \ B R and the trace inequality provide Since η(1 − η) ≡ 0, the far field v far is smooth.Integration by parts using ∆ K((1 − η)u) = 0, as well as [γ 1 Ku] = 0 and η(1 − η) ≡ 0 (therefore no boundary terms appear), leads to .
With the mapping properties of K, K from (3.5), (3.6) and the multiplicative trace inequality this implies Putting the estimates for the near-field and the far-field together, we obtain . (3.43) Step 2: We provide a local ellipticity estimate, i.e., we prove a jn (u, ϕ; η 2 u, η 2 ϕ) + terms in weaker norms.
(See (3.48) for the precise form).We start with (3.43) to obtain The last two terms are already in weaker norms, and for the first two terms, we apply (3.15).Since we assumed C ell > 1/4 for unique solvability, we choose a ρ > 0 such that 1/4 < ρ/2 < C ell and set . (3.45) The first three terms can be expanded as in Theorem 2.3, where (3.16) leads to where the omitted terms (cf.(3.16)) can be estimated in weaker norms (i.e., V ϕ L 2 (B (1+δ/2)R ∩Ω) , u L 2 (B (1+δ/2)R ∩Ω) ) or lead to terms that are absorbed in the left-hand side as in the proof of Theorem 2.3 (see (3.21), (3.22)).Equations (3.37) and (3.17) give Therefore, we only have to estimate the last term in (3.45).We write in the same way as in (3.46) where, again, the omitted terms can be estimated in weaker norms (i.e., by Ku L 2 (B (1+δ/2)R \Γ) and V ϕ L 2 (B (1+δ/2)R ) or absorbed in the left-hand side.Now, integration by parts on R d \Ω and Ω together with ∆ Ku = 0 and Putting everything together into (3.45) and in turn into (3.44),we obtain . (3.48) Step 3: We apply the local orthogonality of (u, ϕ) to piecewise polynomials and use approximation properties.
Therefore, with the super-approximation properties (3.8) of I Ω h , we obtain Let {φ 1 , . . ., φ N } ⊆ X N be a basis of X N .We denote the Galerkin matrix A ∈ R N ×N by The translation of the problem of approximating matrix blocks of A −1 to the problem of approximating certain functions from low dimensional spaces essentially depends on the following crucial property (A1), the existence of a local dual basis.
We denote the coordinate mappings corresponding to the basis and the dual basis by The Hilbert space transpose of Λ is denoted by Λ T .Moreover, for τ ⊂ {1, . . ., N }, we define the sets D j (τ ) := ∪ i∈τ supp λ i,j , where λ i,j is the j-th component of λ i , and write L 2 (τ ) := k+ℓ j=1 L 2 (D j (τ )).
In the following lemma, we derive a representation formula for A −1 based on three linear operators Λ T , S N and Λ.
Lemma 4.2.([AFM20, Lem.3.10], [AFM20, Lem.3.11]) The restriction of Λ T to X N is the inverse mapping Φ −1 .More precisely, for all x, y ∈ R N and v ∈ X N , we have The mappings Λ and Λ T preserve locality, i.e, for τ ⊂ {1, . . ., N } and x ∈ R N with {i : x i = 0} ⊂ τ , we have supp(Λx) ⊂ j D j (τ ).For v ∈ L 2 , we have Moreover, there holds the representation formula Proof.For sake of completeness, we provide the derivation of the representation formula from [AFM20,Lem. 3.11].Using that Λ T = Φ −1 | XN and the definition of the discrete solution operator, we compute This lemma is the crucial step in the proof of the following lemma.
Lemma 4.3.Let A be the Galerkin matrix, Λ be the coordinate mapping for the dual basis, and S N be the discrete solution operator.Let τ × σ ⊂ {1, . . ., N } × {1, . . ., N } be an admissible block and W r ⊆ L 2 be a finite dimensional space.Then, there exist matrices Proof.We use the representation formula from Lemma 4.2 to prove the asserted estimate.With the given space W r , we define X τ σ ∈ R |τ |×r columnwise as vectors from an orthonormal basis of the space W := (Λ T W r )| τ .Then, the product X τ σ X T τ σ is the orthogonal projection onto W. Defining Y τ σ := (A −1 | τ ×σ ) T X τ σ , we can compute for all x ∈ R N with {i : Dividing both sides by x 2 , substituting f := Λx and using that the mapping Λ preserves supports, we get the desired result.
Finally, the question of approximating the whole matrix A −1 can be reduced to the question of blockwise approximation.For arbitrary matrices M ∈ R N ×N , and an arbitrary block partition P of {1, . . ., N } × {1, . . ., N } this follows from If the block partition P is based on a cluster tree T I , the more refined estimate holds, see [Gra01], [Hac09, Lemma 6.5.8], [Bör10b].
In Section 4.3, we give explicit definitions of the dual basis for the FEM-BEM coupling model problem.

Abstract setting -low dimensional approximation
We present a general framework that only uses a Caccioppoli type estimate for the construction of exponentially convergent low dimensional approximations.Let M ∈ N be fixed.For R > 0 let B R := {B i } M i=1 be a collection of boxes, i.e., B i ∈ {B R ∩ Ω, B R , B R \Γ} for all i = 1, . . ., M , where B R denotes a box of side length R. The choice, which of the three sets is taken for each index i, is determined by the application and fixed.
We write B ⊂ B ′ := {B ′ i } M i=1 meaning that B i ⊂ B ′ i for all i = 1, . . ., M .For a parameter δ > 0, we call where B R and B R+2δ are concentric boxes.Defining diam(B R ) := max{diam(B i ), i = 1, . . ., M }, we get In order to simplify notation, we drop the subscript R and write B := B R in the following abstract setting.
We use the notation H 1 (B) to abbreviate the product space Remark 4.4.For the application of the present paper, we chose boxes (or suitable subsets of those) for the sets B i .We also mention that different constructions can be employed as demonstrated in [AFM20], where a construction for non-uniform grids is presented and where the metric is not the Euclidean one but one that is based on the underlying finite element mesh.
In the following, we fix some assumptions on the collections B of interest and the norm |||•||| B on B we derive our approximation result in.In essence, we want a norm weaker than than the classical H 1 -norm that has the correct scaling (e.g., an L 2 -type norm).
(A2) Assumptions on the approximation norm |||•||| B : For each B, the Hilbertian norm |||•||| B is a norm on H 1 (B) and such that for any δ > 0 and enlarged boxes B δ and H > 0 there is a discrete with a constant C Qap > 0 that does not depend on B, B δ , δ, and N .
Finally, we require a Caccioppoli type estimate with respect to the norm from (A2).
(A3) Caccioppoli type estimate: For each B, δ > 0 and collection B δ of δ-enlarged boxes with δ ≥ C Set (N ) with a fixed constant C Set (N ) > 0 that may depend on N , there is a subspace holds.Here, the constants C Cac > 0 and α ≥ 1 do not depend on B, B δ , δ, and N .
We additionally assume the spaces H h (B δ ) to be finite dimensional and nested, i.e., By Π h,B , we denote the orthogonal projection Π h,B : H 1 (B) → H h (B) onto that space with respect to the norm |||•||| B , which is well-defined since, by assumption, H h (B) is closed.
Proof.We set ), we obtain from (A2) and (A3) that with a constant C 1 depending only on Ω since α ≥ 1 and δ ≤ 2 diam(Ω).With the choice H = δ α 2C1CQapCCac diam(B δ ) α−1 , we get the asserted error bound.Since W 1 ⊂ V H,B δ and by choice of H, we have which concludes the proof.
Iterating the single-step approximation on concentric boxes leads to exponential convergence.We iterate the approximation result of Lemma 4.5 on the sets B δℓ , ℓ = L, . . ., 1.For ℓ = L, Lemma 4.5 applied with the sets B (L−1)δ ⊂ B δL provides a subspace ), so we can use Lemma 4.5 again with the sets B (L−2)δ ⊂ B (L−1)δ , and get a subspace Continuing this process L − 2 times leads to the subspace which finishes the proof.

Application of the abstract framework for the FEM-BEM couplings
In this section, we specify the assumptions (A1)-(A3) for the FEM-BEM couplings.
The dual functions λ i are constructed by use of L 2 -dual bases for S 1,1 (T h ) and S 0,0 (K h ).[AFM20, Sec.3.3] gives an explicit construction of a suitable dual basis {λ Ω i : i = 1, . . ., n 1 } for S 1,1 0 (T h ).This is done elementwise in a discontinuous fashion, i.e., λ Ω i ∈ S 1,0 (T h ) ⊂ L 2 (Ω), where each λ Ω i is non-zero only on one element of T h (in the patch of the hat function ξ i ), and the function on this element is given by the push-forward of a dual shape function on the reference element.Moreover, the local stability estimate n j=1 holds for all x ∈ R n , and we have supp λ Ω i ⊂ supp ξ i .We note that the zero boundary condition is irrelevant for the construction.The same can be done for the boundary degrees of freedom, i.e., there exists a dual basis {λ Γ i : i = 1, . . ., n 2 } with the analogous stability and support properties.For the boundary degrees of freedom in S 0,0 (K h ), the dual mappings are given by µ Γ i := χ i / χ i 2 L 2 (Ω) , i.e., the dual basis coincides -up to scaling -with the given basis {χ i : i = 1, . . ., m} of S 0,0 (K h ).With (2.4b), this gives for all y ∈ R m .Now, the dual basis is defined as λ i := (λ Ω i , 0, 0) for i = 1, . . ., n 1 , λ i+n1 := (0, λ Γ i , 0) for i = 1, . . ., n 2 and λ i+n := (0, 0, µ Γ i ) for i = 1, . . ., m, and (4.9), (4.10) together with the analogous one for the λ Γ i show (A1).
)} in a piecewise fashion by (4.12) We denote the patch of an element R ∈ R H by The Scott-Zhang projection reproduces piecewise affine functions and has the following local approximation property for piecewise H s functions: with a constant C depending only on the shape-regularity of R H and d.
δ .We define the operator where I H denotes the classical Scott-Zhang operator for the mesh R H .We have Each term on the right-hand side can be estimated with the same arguments.We only work out the details for the second component.Assuming h ≤ H, and using approximation properties and stability of the Scott-Zhang projection, we get The Caccioppoli inequalities and (A3) Theorem 2.3-Theorem 2.5 provide the Caccioppoli type estimates asserted in (A3) with δ = εR/2.For the Bielak-MacCamy coupling we have α = 1 and C Set = 8h, for the symmetric coupling α = 1 and C Set = 16h.For the Johnson-Nédélec we have to take α = 2 and C Set = 16h.
where the bilinear form a(•, •) is either a sym or a jn .For the Bielak-MacCamy coupling, it suffices to require With these definitions, the closedness and nestedness of the spaces H h (B R ) clearly holds.

Proof of Theorem 2.11
As a consequence of the above discussions, the abstract framework of the previous sections can be applied and it remains to put everything together.
The following proposition constructs the finite dimensional space required from Lemma 4.3, from which the Galerkin solution can be approximated exponentially well.Proposition 4.7 (low dimensional approximation for the symmetric coupling).Let (τ, σ) be a cluster pair with bounding boxes B Rτ and B Rσ that satisfy for given η > 0 Then, for each L ∈ N, there exists a space The constants C low , C box depend only on Ω, d, η, and the γ-shape regularity of the quasi-uniform triangulation T h and K h .
Proof.For given L ∈ N, we choose δ := Rτ 2ηL .Then, we have With B Rτ = {B Rτ ∩ Ω, B Rτ , B Rτ \Γ} and B δL Rτ = {B Rτ +2δL ∩ Ω, B Rτ +2δL , B Rτ +2δL \Γ} from (4.11), the assumption on the support of the data therefore implies the local orthogonality imposed in the space H h (B δL Rτ ).In order to define the space W L , we distinguish two cases.Case δ > 2C Set : Then, Lemma 4.6 applied with the sets B δ Rτ and B δL Rτ provides a space W L of dimension with the approximation properties for Therefore, it remains to estimate the norm |||•||| B from above and below.
With h 1, the mapping properties of V and K from (3.5), and the trace inequality we can estimate The stabilized form a sym (u, ϕ; ψ, ζ) Moreover, [AFF + 13, Thm.18] prove that the Galerkin solution also solves . Therefore, we have The stabilization term can be estimated with the mapping properties of V and K from (3.6) and the trace inequality by Inserting this in (4.17), using the trace inequality and an inverse estimate we further estimate With Young's inequality and inserting this in (4.16), we obtain the upper bound The jump conditions of the single-layer potential and Lemma 3.2 provide for arbitrary ϕ ∈ S 0,0 (K h ) Then, for each L ∈ N, there exists a space W L ⊂ S 1,1 (T h ) × S 0,0 (K h ) with dimension dim W L ≤ C low L 2d+1 such that for arbitrary right-hand sides f ∈ L 2 (Ω), ϕ 0 ∈ L 2 (Γ), and u 0 ∈ L 2 (Γ) with supp f ∪ supp ϕ 0 ∪ supp u 0 ⊂ B Rσ , the corresponding Galerkin solution (u h , ϕ h ) of (2.8) satisfies min The constants C low , C box depend only on Ω, d, η, and the γ-shape regularity of the quasi-uniform triangulation T h and K h .
Proof.The proof is essentially identical to the proof of Proposition 4.7.We stress that the bound of the dimension dim W L ≤ C low L 2d+1 is better, since no approximation for the double-layer potential is needed, i.e., we can choose M = 2 in the abstract setting.
Then, for each L ∈ N, there exists a space W L ⊂ S 1,1 (T h ) × S 0,0 (K h ) with dimension dim W L ≤ C low L 6d+1 , such that for arbitrary right-hand sides f ∈ L 2 (Ω), ϕ 0 ∈ L 2 (Γ), and w 0 ∈ L 2 (Γ) with supp f ∪ supp ϕ 0 ∪ supp w 0 ⊂ B Rσ , the corresponding Galerkin solution (u h , ϕ h ) of (2.15) satisfies min The constants C low , C box depend only on Ω, d, η, and the γ-shape regularity of the quasi-uniform triangulation T h and K h .
Proof.The proof is essentially identical to the proof of Proposition 4.7.We stress that the bound of the dimension dim W L ≤ C low L 6d+1 is worse than for the other couplings, since in the abstract setting, we have to choose M = 3 and α = 2, and the bound follows from Lemma 4.6.
Finally, we can prove the existence of H-Matrix approximants to the inverse FEM-BEM stiffness matrix.
Proof of Theorem 2.11.We start with the symmetric coupling.As H matrices are low rank only on admissible blocks, we set B H | τ ×σ = A −1 sym | τ ×σ for non-admissible cluster pairs and consider an arbitrary admissible cluster pair (τ, σ) in the following.
With a given rank bound r, we take L := ⌊(r/C low ) 1/(3d+1) ⌋.With this choice, we apply Proposition 4.7, which provides a space W L ⊂ S 1,1 (T h ) × S 0,0 (K h ) and use this space in Lemma 4.3, which produces matrices X τ σ , Y τ σ of maximal rank dim W L , which is by choice of L bounded by dim W L = C low L 3d+1 ≤ r.
Proposition 4.7 can be rewritten in terms of the discrete solution operator of the framework of Section 4.1.Let f = (f, v 0 , w 0 ) ∈ L 2 be arbitrary with supp(f ) ⊂ j D j (σ).Then, the locality of the dual functions implies supp f ∪ supp v 0 ∪ supp w 0 ⊂ B Rσ , and we obtain Defining This finishes the proof for the symmetric coupling.
The approximations to A −1 bmc and A −1 jn are constructed in exactly the same fashion.The different exponentials appear due to the different dimensions of the low-dimensional space W L in Proposition 4.8 and Proposition 4.9.

Numerical results
In this section, we provide a numerical example that supports the theoretical results from Theorem 2.11, i.e, we compute an exponentially convergent H-matrix approximant to an inverse FEM-BEM coupling matrix.
If one is only interested in solving a linear system with one (or few) different right-hand sides, rather than computing the inverse -and maybe even its low-rank approximation -it is more beneficial to use an iterative solver.The H-matrix approximability of the inverse naturally allows for black-box preconditioning of the linear system.[Beb07] constructed LU -decompositions in the H-matrix format for FEM matrices by approximating certain Schur-complements under the assumption that the inverse can be approximated with arbitrary accuracy.Theorem 2.11 provides such an approximation result and the techniques of [Beb07,FMP15,FMP16,FMP17] can also be employed to prove the existence of H-LU-decompositions for the whole FEM-BEM matrices for each couplings.
Here, we additionally present a different, computationally more efficient approach by introducing a black-box block diagonal preconditioner for the FEM-BEM coupling matrices.
We choose the 3d-unit cube Ω = (0, 1) 3 as our geometry, and we set C = I.In the following, we only consider the Johnson-Nédélec coupling, the other couplings can be treated in exactly the same way.

4. 3 . 2
Low dimensional approximationThe sets B, B δ and the norm |||•||| BWe take M = 3 and choose collections B = B R := {B R ∩ Ω, B R , B R \Γ}, where B R is a box of side length R. For ℓ ∈ N the enlarged sets B δℓ then have the formB δℓ = B δℓ R := {B R+2δℓ ∩ Ω, B R+2δℓ , B R+2δℓ \Γ} (4.11)with the concentric boxes B R+2δℓ of side length R + 2δℓ.For v = (u, v, w), we use the norm from (2.6)|||v||| B := |||(u, v, w)||| h,R in (A2).For the Bielak-MacCamy coupling, taking M = 2 and choosing collections B R := {B R ∩ Ω, B R } would suffice, however, in order to keep the notation short, we can use M = 3 for this coupling as well by setting the third component to zero, i.e., v = (u, v, 0).The operator Q H and (A2) For the operator Q H , we use a combination of localization and Scott-Zhang interpolation, introduced in [SZ90], on a coarse grid.Since the double-layer potential is discontinuous across Γ, we need to employ a piecewise Scott-Zhang operator.Let R H be a quasi-uniform (infinite) triangulation of R d (into open simplices R ∈ R H ) with mesh width H that conforms to Ω, i.e., every R ∈ R H satisfies either R ⊂ Ω or R ⊂ Ω ext and the restrictions R H | Ω and R H | Ω ext are γ-shape regular, regular triangulations of Ω and Ω ext of mesh size H, respectively.With the Scott-Zhang projections I int H , I ext H for the grids R H | Ω and R H | Ω c , we define the operator B H | τ ×σ := X τ σ Y T τ σ , the estimates (4.3) and Λ h −d/2 together with Lemma 4.3 then give the error boundA −1 sym − B H 2 ≤ C sp depth(T I ) max{ A −1 − B H | τ ×σ 2 : (τ, σ) ∈ P } ≤ C sp depth(T I ) Λ 2 max (τ,σ)∈P far sup f ∈L 2 : supp(f )⊂ j Dj (σ) inf w∈ WL S N f − w L 2 (τ ) f L 2 C sp depth(T I )h −(d+2) 2 −L ≤ C apx C sp depth(T I )h −(d+2) exp(−br 1/(3d+1) ).