Afternote to “Coupling at a Distance”: Convergence Analysis and A Priori Error Estimates

2.1 Problem Setting

Consider a bounded domain Ω 0 ⊂ ℝ 2 that has a smooth parametrizable boundary that will be denoted by Γ 0 := ∂ ⁡ Ω 0 . We will denote the unbounded complement of its closure by

In this section, we will be concerned with the analysis of a discretization for the following diffusion problem:

(2.1)

The function f will be taken to be compactly supported and square integrable on Ω 0 c . The diffusion coefficient 𝜿 is a strictly positive matrix-valued function such that, denoting the identity matrix is as 𝐈 , the difference ( 𝐈 - 𝜿 ) is compactly supported in Ω 0 c . This condition implies that outside of supp ⁡ ( 𝐈 - 𝜿 ) equations (2.1a) and (2.1b) in fact coincide with Poisson’s equation. We will also require that there exist positive constants 𝜿 ¯ and 𝜿 ¯ such that, for any component function κ i ⁢ j of 𝜿 it holds that

The Dirichlet boundary data u 0 will be considered to be an element of the trace space H 1 2 ⁢ ( Γ 0 ) . The radiation condition at infinity (2.1d) is equivalent to assuming that there is a constant u ∞ such that u = u ∞ + 𝒪 ⁢ ( | 𝒙 | - 1 ) (see [18]).

Figure 1

Left: The artificial boundary Γ splits the domain of definition of problem (2.1) into an unbounded region Ω ext and a bounded annular domain Ω. Right: The computational domain Ω h is discretized by an un-fitted triangulation (blue), with boundary Γ h ∪ Γ 0 , h .

2.2 Interior and Exterior Problems

To deal with the unboundedness of the domain, later on we will make use of an integral representation that will reduce the computations to a bounded domain. To this avail, we introduce an artificial, smoothly parametrizable interface Γ enclosing Ω 0 , the support of f and the support of ( 𝐈 - 𝜿 ) . We will also require that Γ ∩ Γ 0 = ∅ . The domain interior to Γ will be denoted Ω, while the unbounded complementary region will be denoted Ω ext . The boundary of Ω will be denoted as ∂ ⁡ Ω and consists of two disjoint components: the artificial boundary Γ and the original problem boundary Γ 0 , so that ∂ ⁡ Ω = Γ ∪ Γ 0 . We will denote the unit normal vector to ∂ ⁡ Ω , pointing in the direction of Ω ext for points in Γ and in the direction of Ω 0 for points in Γ 0 , by 𝒏 . This geometric decomposition, depicted in Figure 1, splits our region of interest into two disjoint domains and allows us to rewrite the problem (2.1) in terms of an interior and an exterior problem coupled by continuity conditions at the artificial boundary Γ.

Since we aim to use an integral equation formulation, for the exterior problem we will prefer a second order formulation and will eliminate 𝒒 ext from the system. We will represent the solutions to (2.1) as the superposition

where the functions u and 𝒒 are supported in Ω, while u ext is supported in Ω ext . The pair ( u , 𝒒 ) satisfies the interior problem

(2.2)

On the other hand, the exterior function u ext satisfies

(2.3)

Above, the boundary value g ∈ H 1 2 ⁢ ( Γ ) corresponds to the trace of u tot over the artificial boundary Γ, while λ ∈ H - 1 2 ⁢ ( Γ ) is the value of the normal flux. These two functions are unknown at this point and will have to be retrieved as part the solution process. However, the knowledge of g (resp. λ) is enough to fully determine the solution to (2.2) or (2.3) considered as independent problems – as long as the equation containing λ (resp. g) is removed from the system. This observation will motivate the alternating solution scheme to be described in Section 5.

2.3 Boundary Integral Formulation for the Exterior Problem

We will now reformulate (2.3) as a boundary integral equation. To do that, we will make use of some standard results from potential theory; we refer the reader interested in further details to the classic references [13, 18] for a comprehensive account, or to [12] for a more concise treatment.

We start by introducing the single layer and double layer potentials defined respectively for η ∈ H 1 2 ⁢ ( Γ ) , μ ∈ H - 1 2 ⁢ ( Γ ) and 𝒙 ∈ ℝ 2 ∖ Γ as

where G ⁢ ( 𝒙 , 𝒚 ) is the Green function for Poisson’s equation. The functions defined by these two potentials satisfy equation (2.3a), and the following jump conditions:

where the jump operator is defined for 𝒚 ∈ Γ and scalar and vector functions v and 𝒗 respectively as

In a similar fashion we can define the average operators as

and use them to define the following boundary integral operators:

We are now in a position to recast the exterior problem (2.3) in terms of boundary integral equations. To that avail, we will represent u ext in Ω ext as

and extend it by zero for 𝒙 ∈ Ω . The constant u ∞ captures the far field behavior of the function and will have to be determined. Since u ext ≡ 0 in Ω, by applying the integral operators above to the integral representation (2.6), the boundary condition (2.3b) leads to

giving rise to the integral equation

(2.7)

To ensure that u ext = u ∞ as | 𝒙 | → ∞ , we must impose the additional restriction

Equation (2.7a) will be used as part of the alternating scheme described in Section 5, where an approximation of λ will be produced by a numerical solution of the interior problem (2.2) and the density g solving (2.7a) will be then used as the Dirichlet datum for (2.2).

Therefore, if Γ has two continuous derivatives and λ ∈ H - 1 2 ⁢ ( Γ ) is problem data satisfying the constraint (2.7b), then the unique solvability of equation (2.7) and continuous dependence on problem data follow from standard results in boundary integral equations (see, for instance [14, Section 6.4]). Moreover, there exists a constant c > 0 , depending only on Γ and the norms of ( 1 2 - 𝒦 ) - 1 and 𝒱 , such that

Moreover, from this estimate and the representation formula (2.6), it follows that there exists C BIE > 0 such that

2.4 Variational Formulation for the Interior Problem

In this subsection, we will study the interior Dirichlet boundary value problem obtained from (2.2) by removing (2.2d) altogether and considering that the boundary trace g, appearing in (2.2c), is known. This yields the problem

(2.10)

Above, the source term f ∈ L 2 ⁢ ( Ω ) and the Dirichlet boundary data ξ 0 ∈ H 1 2 ⁢ ( ∂ ⁡ Ω ) is given by

To derive the weak formulation of this system, we test (2.10a) with an arbitrary w ∈ L 2 ⁢ ( Ω ) and (2.10b) with 𝒗 ∈ 𝑯 ⁢ ( 𝐝𝐢𝐯 ; Ω ) , integrate by parts and incorporate (2.10c) leading to

where ( ⋅ , ⋅ ) Ω and 〈 ⋅ , ⋅ 〉 ∂ ⁡ Ω denote the L 2 -inner products over Ω and ∂ ⁡ Ω , respectively. From the three preceding equations, we arrive at the variational problem:

The well-posedness of (2.11) follows from standard arguments of Babǔska–Brezzi theory [11, Section 2.4] and the solution satisfies

We will, however, not solve the problem as stated above and instead will consider a slightly different version posed in a subdomain. This approach, known as the transfer path method will be described in detail in Section 3.2, and will require us first to discuss the geometric setting of the discretization, which we will do next.

3.1 Geometric Setting and Notation

3.2 Transferal of Boundary Conditions

Having introduced all the necessary geometric concepts we can now return to the interior problem (2.10) which we will now pose in a polygonal computational domain Ω h ⊂ Ω satisfying the admissibility requirements discussed in the previous subsection. In addition, we will need to define a bijective^[1] mapping

assigning a point 𝒙 ¯ ∈ ∂ ⁡ Ω to every point 𝒙 ∈ ∂ ⁡ Ω h .

For any fixed computational domain Ω h , the solution pair to (2.11) satisfies the related problem

(3.2)

where the boundary condition φ 0 𝒒 can be calculated by integrating equation (2.10b) along a path connecting ∂ ⁡ Ω to ∂ ⁡ Ω h . More precisely, if we denote the distance between 𝒙 and 𝒙 ¯ by l ⁢ ( 𝒙 ) , and by 𝒕 the unit vector 𝒙 ¯ - 𝒙 | 𝒙 ¯ - 𝒙 | , the boundary conditions on Γ h can be expressed in terms of the flux 𝒒 and the trace of u on ∂ ⁡ Ω , as

Note that the required bijectivity of ϕ ⁢ ( 𝒙 ) implies that 𝒕 can not be tangent to a boundary edge. Thus, the solution of (2.11) also satisfies the abstract formulation

where the bilinear forms 𝒜 : 𝑯 ⁢ ( 𝐝𝐢𝐯 ; Ω h ) × 𝑯 ⁢ ( 𝐝𝐢𝐯 ; Ω h ) → ℝ , ℬ : 𝑯 ⁢ ( 𝐝𝐢𝐯 ; Ω h ) × L 2 ⁢ ( Ω h ) → ℝ , and the functionals ℱ 1 : 𝑯 ⁢ ( 𝐝𝐢𝐯 ; Ω h ) → ℝ and ℱ 2 : L 2 ⁢ ( Ω h ) → ℝ are defined by

Beyond the difference in the domain of definition, the system above differs from the original problem (2.11) in the presence of the term 𝒜 T , introduced by the transfer of boundary condition. The well posedness of problems of this form was established in [20]. On the interest of brevity, we shall not repeat the argument here and instead will now discuss the discretization of this problem along with that of the integral equation (2.7).

3.3 Discrete Variational Formulation

Having defined all the required notation, we can now state the HDG discretization of (2.10) which, for Dirichlet data ξ 0 ∈ H 1 2 ⁢ ( ∂ ⁡ Ω ) , seeks an approximation ( 𝒒 h , u h , u ^ h ) ∈ 𝑽 h × W h × M h satisfying

(3.4)

for any test ( 𝒗 , w , μ ) ∈ 𝑽 h × W h × M h . Following [8], the approximate boundary data on ∂ ⁡ Ω h appearing on the right-hand side of (3.4d) is given by

where E denotes the extrapolation operator. The numerical flux in the normal direction 𝒒 ^ h ⋅ 𝒏 h is defined as

where τ stabilization function. Throughout this analysis we will only require 0 < τ ≤ τ ¯ < ∞ , where τ ¯ denotes the maximum value of τ.

Note that the terms 〈 u ^ h , 𝒗 ⋅ 𝒏 h 〉 ∂ ⁡ 𝒯 h and 〈 τ ⁢ u ^ h , w 〉 ∂ ⁡ 𝒯 h , given in (3.4a) and (3.4b), respectively, can be split into the contributions of the interior edges and of the boundary edges as

Replacing now the numerical flux (3.5) in (3.4c) results in

In order to apply known results from functional analysis, we rewrite the numerical trace u ^ h in terms of averages and jumps. For this, we use the equation (3.4c) and separate the term featuring u ^ h as

Above, we have used the fact that the hybrid variable u ^ h is single valued, and the average { { ⋅ } } and jump [ [ ⋅ ] ] operators are defined for every edge e in a fashion analogous to (2.4) and (2.5). Then, taking as test function

in the expression above, we deduce that

We make use of this identity to obtain

and

In this way, replacing the definition of φ 0 𝒒 h – see (3.4e) – in (3.4a) and (3.4b), together with the foregoing identities, we obtain that (3.4) is equivalent to finding ( 𝒒 h , u h ) ∈ 𝑽 h × W h such that

(3.6)

where the bilinear forms 𝒜 h : 𝑽 h × 𝑽 h → ℝ , ℬ h , ℬ T : 𝑽 h × W h → ℝ , 𝒞 h : W h × W h → ℝ , and the functionals ℱ 1 , h : 𝑽 h → ℝ and ℱ 2 , h : W h → ℝ are defined by

(3.7)

The unique solvability of the scheme (3.6) will be proved by an energy argument. To that end, for e ∈ ∂ ⁡ Ω h and 𝒗 ∈ 𝑳 2 ⁢ ( T e ext ) , it is convenient to define the following norm on the extension patch T e ext :

This norm is equivalent to the standard 𝑳 2 ⁢ ( T e ext ) -norm as shown first in [20] for the two-dimensional and later extended to three dimensions in [21]. That is, there exist positive constants C 1 e and C 2 e , independent of h, such that

This equivalence holds true under certain conditions on the transferring vectors 𝒕 ⁢ ( 𝒙 ) (cf. [21, 20])) ensuring, roughly speaking, that they cannot deviate too much from the vector normal to e.

We also introduce the element-wise constants

where 𝒱 k := { 𝒑 ∈ [ ℙ k ⁢ ( T e e ⁢ x ⁢ t ∪ T e ) ] 2 : 𝒑 ≠ 𝟎 } . These constants are independent of h, but depend on the polynomial degree k and the mesh regularity parameter as shown in [5].