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

Nestor Sánchez , Tonatiuh Sánchez-Vizuet and Manuel E. Solano

## Abstract

In their article “Coupling at a distance HDG and BEM”, Cockburn, Sayas and Solano proposed an iterative coupling of the hybridizable discontinuous Galerkin method (HDG) and the boundary element method (BEM) to solve an exterior Dirichlet problem. The novelty of the numerical scheme consisted of using a computational domain for the HDG discretization whose boundary did not coincide with the coupling interface. In their article, the authors provided extensive numerical evidence for convergence, but the proof of convergence and the error analysis remained elusive at that time. In this article we fill the gap by proving the convergence of a relaxation of the algorithm and providing a priori error estimates for the numerical solution.

MSC 2010: 65N15; 65N30; 65R20

## 1 Introduction

The goal of this article is to conclude the work started by Cockburn, Sayas and Solano in the article Coupling at a distance [7], where an iterative solution method for a classic exterior elliptic problem was introduced. The proposed scheme amounted to a Schur complement-style algorithm that alternates between a Hybridizable Discontinuous Galerkin Method (HDG) for an interior problem and the Boundary Element Method (BEM) for an exterior problem. At the time of publication, the novelty of the method resided in the use of non-touching grids for the discretization of each of the two problems. The ready availability of two separate, uncoupled, codes for each of the discretization methods and the eagerness to show the viability of such a non-touching coupling led to the choice of an iterative alternating procedure – even though the problem in question is in fact linear.

When [7] was published, the technique for transferring information between the two grids had only been recently incorporated into the HDG literature [8] and, despite the fact that convincing numerical evidence of convergence at an optimal rate was provided, a rigorous analysis of the coupled scheme proved elusive at the time. A few years after Coupling at a distance appeared, a method for the analysis of HDG discretizations involving the transfer technique – that we now like to call the transfer path method – was developed in [5] for interior elliptic problems. Since then, both the transfer technique and the analysis method have been successfully employed for the study of linear [21, 32, 33], and non-linear [22, 25, 24, 27, 28] interior problems, as well as problems with interfaces [23, 31], however the analysis of the HDG-BEM coupling had fallen by the wayside and remained unfinished.

The current special issue honoring Francisco-Javier Sayas, one of the co-authors of the original article, seemed like the perfect venue for the missing analysis. In that sense, the present communication shall not be considered a novel contribution, but rather the conclusion, long overdue, of the original work, an after-note to the original work Coupling at a distance. With that in mind, we will stick to the iterative alternating procedure proposed in [7], even if a more efficient monolithic approach where the HDG and BEM discrete systems – along with the discrete coupling terms – are solved simultaneously is possible. The study of such a monolithic scheme applied to non-linear problems is the subject of ongoing work that will be communicated in a separate publication [26].

The method proposed in [7], rather than approaching the problem as a single coupled unit, follows the spirit of domain decomposition methods. It relies on an iterative approximation of a Dirichlet to Neumann mapping through the independent solution of an interior and an exterior problem that communicate through their Dirichlet and Neumann traces. Since these two problems are dealt with independent solvers, we will analyze their discretizations separately. After establishing the well posedness of the independent discretizations, we will then prove that, at the discrete level, the alternating solution of an interior Dirichlet and (with HDG) an exterior Neumann problem (with BEM) converges to the solution of the original unbounded problem. This latter result constitutes the main contribution of this article.

We will describe the problem setting and its reformulation as a system of coupled interior/exterior problems at the continuous level in Section 2. The discretizations of the interior problem and the boundary integral formulation for the exterior problem are described respectively in Sections 3 and 4. Finally, in Section 5, we show that it is possible to define a relaxation of the iterative process presented in [7], alternating between the solution of the interior and the boundary problems, that converges to the solution of the original problem.

## 2 Continuous Formulation

### 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

Ω 0 c := 2 Ω 0 ¯ .

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

(2.1)

(2.1a) 𝒒 tot = f in  Ω 0 c ,
(2.1b) 𝒒 tot + 𝜿 u tot = 0 in  Ω 0 c ,
(2.1c) u tot = u 0 on  Γ 0 ,
(2.1d) u tot = 𝒪 ( 1 ) as  𝒙 .

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

𝜿 ¯ κ i j ( 𝒙 ) 𝜿 ¯ for all  𝒙 Ω .

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

u tot = u + u ext and 𝒒 tot = 𝒒 + u ext ,

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

(2.2)

(2.2a) 𝒒 = f in  Ω ,
(2.2b) 𝒒 + 𝜿 u = 0 in  Ω ,
(2.2c) u = g on  Γ ,
(2.2d) 𝒒 𝒏 = λ on  Γ ,
(2.2e) u = u 0 on  Γ 0 .

On the other hand, the exterior function u ext satisfies

(2.3)

(2.3a) - Δ u ext = 0 in  Ω ext ,
(2.3b) u ext = g on  Γ ,
(2.3c) u ext 𝒏 = - λ on  Γ ,
(2.3d) u tot = 𝒪 ( 1 ) as  | 𝒙 | .

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

𝒮 μ ( 𝒙 ) := Γ G ( 𝒙 , 𝒚 ) μ ( 𝒚 ) 𝑑 Γ 𝒚 (single layer) ,
𝒟 η ( 𝒙 ) := Γ 𝒏 ( 𝒚 ) G ( 𝒙 , 𝒚 ) η ( 𝒚 ) 𝑑 Γ 𝒚 (double layer) ,

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:

[ [ 𝒮 μ ] ] := 0 , [ [ 𝒮 μ ] ] := μ , [ [ 𝒟 η ] ] := - η , [ [ 𝒟 η ] ] := 0 ,

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

(2.4) [ [ v ] ] := lim ϵ 0 ( v ( 𝒚 - ϵ 𝒏 ) - v ( 𝒚 + ϵ 𝒏 ) ) and [ [ 𝒗 ] ] := lim ϵ 0 ( 𝒗 ( 𝒚 - ϵ 𝒏 ) - 𝒗 ( 𝒚 + ϵ 𝒏 ) ) 𝒏 ( 𝒚 ) .

In a similar fashion we can define the average operators as

(2.5) { { v } } := 1 2 lim ϵ 0 ( v ( 𝒚 - ϵ 𝒏 ) + v ( 𝒚 + ϵ 𝒏 ) ) and { { 𝒗 } } := 1 2 lim ϵ 0 ( 𝒗 ( 𝒚 - ϵ 𝒏 ) + 𝒗 ( 𝒚 + ϵ 𝒏 ) ) 𝒏 ( 𝒚 ) ,

and use them to define the following boundary integral operators:

𝒱 μ := { { 𝒮 μ } } , 𝒦 μ := { { ( 𝒮 μ ) } } , 𝒦 η := { { 𝒟 η } } and 𝒲 η := - { { ( 𝒟 η ) } } .

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

(2.6) u ext = 𝒟 g - 𝒮 λ + u

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

{ { u ext } } = 1 2 g = 𝒦 g - 𝒱 λ + 1 2 u ,

giving rise to the integral equation

(2.7)

(2.7a) ( 1 2 - 𝒦 ) g = - 𝒱 λ + 1 2 u .

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

(2.7b) λ , 1 Γ = 0 .

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

(2.8) g 1 / 2 , Γ c λ - 1 / 2 , Γ .

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

(2.9) u ext Ω C BIE λ - 1 / 2 , Γ + | u | .

### 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)

(2.10a) 𝒒 = f in  Ω ,
(2.10b) 𝜿 - 1 𝒒 + u = 0 in  Ω ,
(2.10c) u = ξ 0 on  Ω .

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

ξ 0 = { u 0 on  Γ 0 , g on  Γ .

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

( 𝒒 , w ) Ω = ( f , w ) Ω ,
( 𝜿 - 1 𝒒 , 𝒗 ) Ω - ( u , 𝒗 ) Ω = - 𝒗 𝒏 , ξ 0 Ω ,

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

### Problem.

Find ( 𝒒 , u ) 𝑯 ( 𝐝𝐢𝐯 ; Ω ) × L 2 ( Ω ) such that

(2.11)

(2.11a) 𝒜 ~ ( 𝒒 , 𝒗 ) + ~ ( 𝒗 , u ) = ~ 1 ( 𝒗 ) for all  𝒗 𝑯 ( 𝐝𝐢𝐯 ; Ω ) ,
(2.11b) ~ ( 𝒒 , w ) = ~ 2 ( w ) for all  w L 2 ( Ω ) ,

where the bilinear forms 𝒜 ~ : 𝑯 ( 𝐝𝐢𝐯 ; Ω ) × 𝑯 ( 𝐝𝐢𝐯 ; Ω ) , ~ : 𝑯 ( 𝐝𝐢𝐯 ; Ω ) × L 2 ( Ω ) , and the functionals 1 ~ : 𝑯 ( 𝐝𝐢𝐯 ; Ω ) and 2 ~ : L 2 ( Ω ) are defined by

𝒜 ~ ( 𝒒 , 𝒗 ) := ( 𝜿 - 1 𝒒 , 𝒗 ) Ω ,
~ ( 𝒒 , w ) := - ( w , 𝒒 ) Ω ,
~ 1 ( 𝒗 ) := - ξ 0 , 𝒗 𝒏 Ω ,
~ 2 ( w ) := - ( f , w ) Ω .

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

(2.12) 𝒒 div , Ω + u 0 , Ω C stab 𝜿 ¯ 1 2 ( f 0 , Ω + ξ 0 1 / 2 , Ω ) = C stab 𝜿 ¯ 1 2 ( f 0 , Ω + g 1 / 2 , Γ + u 0 1 / 2 , Γ 0 ) .

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 HDG Discretization of the Interior Problem

### 3.1 Geometric Setting and Notation

#### The Computational Domain.

We will consider a family of polygonal subdomains Ω h Ω that approximate Ω in the sense that the Lebesgue measure μ ( Ω Ω h ) 0 , as h 0 . We will refer to any such Ω h as a computational domain and will triangulate Ω ¯ h by a shape-regular triangulation 𝒯 h as depicted in Figure 1. A generic element in 𝒯 h will be denoted by T and the mesh parameter h will be defined as diameter of a circle inscribing an element T 𝒯 h . The set 𝒯 h := { T : T 𝒯 h } , will be referred to as the skeleton of the triangulation. The set of edges, e, of 𝒯 h will be denoted by h and we will distinguish between those edges lying entirely in the computational boundary

h := { e h : e Ω h = e } ,

and those that are either interior or have at most their endpoints in the computational boundary

h := { e h : e Ω h e } .

We will refer to the former as boundary edges and to the latter as interior edges. Note that h = h h .

Just as the boundary associated to the continuous problem (2.2) has two separate connected components, the boundary of the computational domain can be split as Ω h = Γ h Γ h , 0 , where

Γ h := { e h : d ( e , Γ ) d ( e , Γ 0 ) } and Γ h , 0 := { e h : d ( e , Γ 0 ) < d ( e , Γ ) } .

We will require that the computational domain Ω h and the triangulation 𝒯 h satisfy the following local proximity condition: for any point in the computational boundary Ω h , the minimum distance between 𝒙 and the boundary Ω = Γ Γ 0 should be, at most, of the same order of magnitude as the diameter of the smallest triangle T 𝒯 h , such that 𝒙 T . In view of this condition, the process of mesh refinement should not be understood as a sequence of finer triangulations for a fixed computational domain Ω h . Instead, as the mesh diameter h 0 , the process involves the passage through a sequence of pairs domain/triangulation ( Ω h , 𝒯 h ) that satisfy the local proximity condition and exhaust the original domain Ω as the refinement progresses. We refer the reader to [24], where this condition is discoursed in more detail, and to [28] where an algorithm for building a sequence { ( Ω h , 𝒯 h ) } h is described.

#### Mesh-Dependent Subspaces and Inner Products.

For the discrete formulation we will have introduce the following mesh-dependent inner products:

( u , w ) 𝒯 h := T 𝒯 h T u w for all  u , w L 2 ( 𝒯 h ) ,
( 𝒒 , 𝒗 ) 𝒯 h := T 𝒯 h T 𝒒 𝒗 for all  𝒒 , 𝒗 𝑳 2 ( 𝒯 h ) ,
u , w 𝒯 h := T 𝒯 h T u w for all  u , w L 2 ( 𝒯 h ) ,
u , w 𝒯 h Γ h := T 𝒯 h e T Γ h e u w for all  u , w L 2 ( 𝒯 h ) ,
u , w 𝒯 h Γ h , 0 := T 𝒯 h e T Γ h , 0 e u w for all  u , w L 2 ( 𝒯 h ) .

These inner products induce mesh-dependent norms that will be denoted, respectively, by

w Ω h := ( w , w ) 𝒯 h 1 2 , w 𝒯 h := w , w 𝒯 h 1 2 and w Γ h := w , w 𝒯 h Γ h 1 2 .

The finite-dimensional discontinuous polynomial subspaces that will be used for discretization, for k 0 , are given by

𝑽 h := { 𝒗 𝑳 2 ( 𝒯 h ) : 𝒗 | T [ k ( T ) ] 2  for all  T 𝒯 h } ,
W h := { w L 2 ( 𝒯 h ) : w | T k ( T )  for all  T 𝒯 h } ,
M h := { μ L 2 ( h ) : μ | T k ( F )  for all  F h } ,

where k ( T ) denotes the space of polynomials of degree at most k defined in T 𝒯 h . Similarly, k ( e ) denotes the space of polynomials of degree at most k defined over a edge e h .

#### Extension Patches and Extrapolation.

Since the discrete spaces are defined only over the elements of the triangulation we will need to define a way to compute our approximations in the region Ω Ω h between the boundary and the computational boundary. To this purpose, we will tessellate this region as follows:

1. Let 𝒙 1 and 𝒙 2 be the endpoints of a boundary edge e Ω h .

2. Let 𝒙 ¯ 1 and 𝒙 ¯ 2 be the corresponding points in Ω – as determined by the mapping (3.1).

3. Let 𝝈 1 and 𝝈 2 be the straight segments connecting 𝒙 ¯ 1 to 𝒙 1 and 𝒙 ¯ 2 to 𝒙 2 .

We will refer to the open region of Ω Ω h delimited by e, 𝝈 1 and 𝝈 2 and the segment of Ω connecting 𝒙 ¯ 1 to 𝒙 ¯ 2 as an extension patch and will denote it by T e ext . It is clear that for every e Γ h there is one and only one such T e ext (this justifies subindex in the notation) and that

Ω Ω h ¯ = e Γ h T ¯ e ext .

It also follows from this construction that for every T e ext there is only one element T e in the triangulation such that K ¯ e ext T ¯ e = e . We will use this fact to define an extrapolation operator that will extend the value of the piecewise polynomial functions defined on T e onto the corresponding extension patch T e ext , thus extending functions the discrete spaces above into the full domain Ω. With this in mind, we will define the values of polynomial function p on T e ext by extrapolating the values of the corresponding polynomial from T e , and will denote its as E p ( 𝒙 ) for any 𝒙 T e ext .

For a given domain Ω h and corresponding triangulation 𝒯 h , the usual notion of the exterior normal vector is well defined for almost all points in the boundary, with the possible exception of the vertices of the triangulation. We will define the exterior normal vector to the computational domain, 𝒏 h in the usual manner, and extend the definition to 𝒏 h ( 𝒙 ) = 𝒕 ( 𝒙 ) , where 𝒕 ( 𝒙 ) := 𝒙 ¯ - 𝒙 | 𝒙 ¯ - 𝒙 | , for those vertices for which the standard normal vector is not well defined. On the other hand, we will define the unit normal vector exterior to each element T T h as 𝒏 h , which will coincide with the exterior normal 𝒏 h on element edges belonging to the computational boundary Γ h .

Finally, for every edge e h we will denote the ratio between its distance to the boundary and the diameter, h T e , of its parent element as

r e := d ( e , Ω ) h T e ,

and will define the boundary proximity parameter as

R h := max e h r e ,

and will assume for this work that the family of admissible domains and triangulations ( Ω h , 𝒯 h ) is such that:

1. R h 0 as h 0 ,

2. 𝒏 h - 𝒏 = o ( h 1 2 ) as h 0 , where the normal 𝒏 h should be understood as coinciding with 𝝈 for those points in which the standard normal vector is not defined.

### 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

(3.1) ϕ : Ω h Ω , 𝒙 𝒙 ¯

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)

(3.2a) 𝒒 = f in  Ω h ,
(3.2b) 𝜿 - 1 𝒒 + u = 0 in  Ω h ,
(3.2c) u = φ 0 𝒒 on  Ω h ,

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

(3.3) φ 0 𝒒 ( 𝒙 ) := ξ 0 ϕ ( 𝒙 ) + 0 l ( 𝒙 ) 𝜿 - 1 𝒒 ( 𝒙 + 𝒕 ( 𝒙 ) s ) 𝒕 ( 𝒙 ) 𝑑 s for all  𝒙 Ω h .

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

𝒜 ( 𝒒 , 𝒗 ) + 𝒜 T ( 𝒒 , 𝒗 ) + ( 𝒗 , u ) = 1 ( 𝒗 ) for all  𝒗 𝑯 ( 𝐝𝐢𝐯 ; Ω h ) ,
( 𝒒 , w ) = 2 ( w ) for all  w L 2 ( Ω h ) ,

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

𝒜 ( 𝒒 , 𝒗 ) := ( 𝜿 - 1 𝒒 , 𝒗 ) Ω h ,
𝒜 T ( 𝒒 , 𝒗 ) := e Ω h e ( 0 l ( 𝒙 ) 𝜿 - 1 𝒒 ( 𝒙 + 𝒕 ( 𝒙 ) s ) 𝒕 ( 𝒙 ) ) 𝒗 ( 𝒙 ) 𝒏 h 𝑑 s 𝑑 S 𝒙 ,
( 𝒒 , w ) := - ( w , 𝒒 ) Ω h ,
1 ( 𝒗 ) := - ξ 0 ϕ , 𝒗 𝒏 h Ω h ,
2 ( w ) := - ( f , w ) Ω h .

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)

(3.4a) ( 𝜿 - 1 𝒒 h , 𝒗 ) 𝒯 h - ( u h , 𝒗 ) 𝒯 h + u ^ h , 𝒗 𝒏 h 𝒯 h = 0 ,
(3.4b) ( 𝒒 h , w ) 𝒯 h + τ u h , w 𝒯 h - τ u ^ h , w 𝒯 h = ( f , w ) 𝒯 h ,
(3.4c) μ , 𝒒 ^ h 𝒏 h 𝒯 h Ω h = 0 ,
(3.4d) u ^ h , μ Ω h = φ 0 𝒒 h , μ Ω h ,

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

(3.4e) φ 0 𝒒 h ( 𝒙 ) := ξ 0 ϕ ( 𝒙 ) + 0 l ( 𝒙 ) 𝜿 - 1 E 𝒒 h ( 𝒙 + 𝒕 ( 𝒙 ) s ) 𝒕 ( 𝒙 ) 𝑑 s for  𝒙 Ω h ,

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

(3.5) 𝒒 ^ h 𝒏 h = 𝒒 h 𝒏 h + τ ( u h - u ^ h ) on  𝒯 h ,

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

u ^ h , 𝒗 𝒏 h 𝒯 h = u ^ h , 𝒗 𝒏 h 𝒯 h Ω h + φ 0 𝒒 h , 𝒗 𝒏 h Ω h ,
τ u ^ h , w 𝒯 h = τ u ^ h , w 𝒯 h Ω h + τ φ 0 𝒒 h , w Ω h .

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

μ , 𝒒 h 𝒏 h 𝒯 h Ω h + μ , τ ( u h - u ^ h ) 𝒯 h Ω h = 0 .

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

0 = μ , 𝒒 h 𝒏 h 𝒯 h Ω h + μ , τ u h 𝒯 h Ω h - μ , τ u ^ h 𝒯 h Ω h
= T 𝒯 h e T Ω h e ( μ 𝒒 h 𝒏 h + τ μ u h - τ μ u ^ h )
= e h e ( [ [ 𝒒 h ] ] μ + 2 τ { { u h } } μ - 2 τ u ^ h μ )
= h ( [ [ 𝒒 h ] ] + 2 τ { { u h } } - 2 τ u ^ h ) μ for all  μ M h .

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

μ = [ [ 𝒒 h ] ] + 2 τ { { u h } } - 2 τ u ^ h M h

in the expression above, we deduce that

u ^ h = 1 2 τ - 1 [ [ 𝒒 h ] ] + { { u h } } on  h .

We make use of this identity to obtain

u ^ h , 𝒗 𝒏 h 𝒯 h Ω h = u ^ h , [ [ 𝒗 ] ] h = 1 2 τ - 1 [ [ 𝒒 h ] ] , [ [ 𝒗 ] ] h + { { u h } } , [ [ 𝒗 ] ] h

and

τ u ^ h , w 𝒯 h Ω h = 2 τ { { w } } , u ^ h h = [ [ 𝒒 h ] ] , { { w } } h + 2 τ { { w } } , { { u h } } h .

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)

(3.6a) 𝒜 h ( 𝒒 h , 𝒗 ) + 𝒜 T ( 𝒒 h , 𝒗 ) + h ( 𝒗 , u h ) = 1 , h ( 𝒗 ) for all  𝒗 𝑽 h ,
(3.6b) T ( 𝒒 h , w ) + h ( 𝒒 h , w ) - 𝒞 h ( u h , w ) = 2 , h ( w ) for all  w W h ,

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)

(3.7a) 𝒜 h ( 𝒒 h , 𝒗 ) := ( 𝜿 - 1 𝒒 h , 𝒗 ) 𝒯 h + 1 2 τ - 1 [ [ 𝒒 h ] ] , [ [ 𝒗 ] ] h ,
(3.7b) 𝒜 T ( 𝒒 , 𝒗 ) := e Ω h e ( 0 l ( 𝒙 ) 𝜿 - 1 𝒒 ( 𝒙 + 𝒕 ( 𝒙 ) s ) 𝒕 ( 𝒙 ) ) 𝒗 ( 𝒙 ) 𝒏 h 𝑑 s 𝑑 S 𝒙 ,
(3.7c) h ( 𝒒 h , w ) := - ( w , 𝒒 h ) 𝒯 h + [ [ 𝒒 h ] ] , { { w } } h
(3.7d) T ( 𝒒 h , w ) := e Ω h e τ ( 0 l ( 𝒙 ) 𝜿 - 1 𝒒 h ( 𝒙 + 𝒕 ( 𝒙 ) s ) 𝒕 ( 𝒙 ) ) w ( 𝒙 ) 𝑑 s 𝑑 S 𝒙 ,
(3.7e) 𝒞 h ( u h , w ) := τ u h , w 𝒯 h - 2 τ { { u h } } , { { w } } h ,
(3.7f) 1 , h ( 𝒗 ) := - ξ 0 ϕ , 𝒗 𝒏 h Ω h ,
(3.7g) 2 , h ( w ) := - ( f , w ) 𝒯 h - τ ξ 0 ϕ , w Ω h .

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 :

|| | 𝒗 | || e := ( e 0 l ( 𝒙 ) | 𝒗 ( 𝒙 + s 𝒕 ( 𝒙 ) ) | 2 𝑑 s 𝑑 S 𝒙 ) 1 2 .

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

(3.8) C 1 e || | 𝒗 | || e 𝒗 T e ext C 2 e || | 𝒗 | || e .

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

(3.9) C ext e := 1 r e sup 𝝌 𝒱 k || | 𝝌 | || e 𝝌 T e and C inv e := h e sup 𝝌 𝒱 k || | 𝝌 | || T e 𝝌 T e ,

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].

We now proceed to derive an energy inequality that will lead to the well-posedness of (3.6).

### Lemma 1.

Let α h = R h 𝛋 ¯ - 1 ( 𝛋 ¯ - 𝛋 ¯ 1 2 h 1 2 τ ¯ 1 2 ) and β h = 𝛋 ¯ - 1 𝛋 ¯ 1 2 R h h 1 2 τ ¯ 1 2 . It holds

( 1 - α h ) 𝜿 - 1 2 𝒒 h 0 , Ω h 2 + ( 1 - β h ) τ 1 2 u h Ω h 2 + τ 1 2 ( u h - { { u h } } ) 𝒯 h Ω h 2 + τ - 1 2 [ [ 𝒒 h ] ] h 2
(3.10) 𝜿 1 2 h - 1 2 ξ 0 ϕ Ω h 2 + f 0 , Ω u h 0 , Ω h .

### Proof.

By taking 𝒗 = 𝒒 h and w = u h in (3.6), and subtracting the resulting expressions we obtain

𝜿 - 1 2 𝒒 h 0 , Ω h 2 + 1 2 τ - 1 2 [ [ 𝒒 h ] ] h 2 + 𝒜 T ( 𝒒 h , 𝒒 h ) + T ( 𝒒 h , u h ) + 𝒞 h ( u h , u h )
(3.11) = - ξ 0 ϕ , 𝒒 h 𝒏 h Ω h - ( f , u h ) 𝒯 h - τ ξ 0 ϕ , u h Ω h .

First of all, after performing algebraic calculations, we observe that 𝒞 is a semi-definite operator from W h × W h to . In fact,

(3.12) 𝒞 h ( u h , u h ) = τ u h , u h 𝒯 h - 2 τ { { u h } } , { { u h } } h = τ 1 2 ( u h - { { u h } } ) 𝒯 h Ω h 2 + τ 1 2 u h Ω h 2 .

We will now obtain a lower bound for the non-positive terms of left-hand side of (3.11). In this direction, the operator 𝒜 T can be bounded as follows. Let e Ω h and 𝒙 e . By the Cauchy–Schwarz inequality and the definition in (3.9),

0 l ( 𝒙 ) 𝜿 - 1 𝒒 h ( 𝒙 + 𝒕 ( 𝒙 ) s ) 𝒕 ( 𝒙 ) 𝑑 s l ( 𝒙 ) 1 2 || | 𝜿 - 1 𝒒 h | || e h T e 1 2 r e 𝜿 ¯ - 1 𝜿 ¯ 1 2 C ext e 𝜿 - 1 2 𝒒 h T e ,

where we have used the bound l ( 𝒙 ) h T e r e . Then, by the discrete trace inequality, we have

- 𝒜 T ( 𝒒 h , 𝒒 h ) | 𝒜 T ( 𝒒 h , 𝒒 h ) | 𝜿 ¯ - 1 𝜿 ¯ 1 2 e Ω h h T e 1 2 r h 𝜿 - 1 2 𝒒 h e 𝒒 h 𝒏 h e
(3.13) R h 𝜿 ¯ - 1 𝜿 ¯ 𝜿 - 1 2 𝒒 h Ω h 2 .

The same arguments yield to

- T ( 𝒒 h , u h ) | T ( 𝒒 h , u h ) | 𝜿 ¯ - 1 𝜿 ¯ 1 2 R h h 1 2 τ ¯ 1 2 𝜿 - 1 2 𝒒 h 0 , Ω h τ 1 2 u h 0 , Ω h
(3.14) 𝜿 ¯ - 1 𝜿 ¯ 1 2 R h h 1 2 τ ¯ 1 2 ( 1 2 𝜿 - 1 2 𝒒 h 0 , Ω h 2 + 1 2 τ 1 2 u h 0 , Ω h 2 ) .

Therefore, combining the above estimates and (3.11), we deduce that

( 1 - R h 𝜿 ¯ - 1 𝜿 ¯ - 𝜿 ¯ - 1 𝜿 ¯ 1 2 R h h 1 2 τ ¯ 1 2 ) 𝜿 - 1 2 𝒒 h 0 , Ω h 2 + ( 1 - 𝜿 ¯ - 1 𝜿 ¯ 1 2 R h h 1 2 τ ¯ 1 2 ) τ 1 2 u h Ω h 2
+ τ - 1 2 [ [ 𝒒 h ] ] h 2 + τ 1 2 ( u h - { { u h } } ) 𝒯 h Ω h 2
| ξ 0 ϕ , 𝒒 h 𝒏 h Ω h | + | ( f , u h ) 𝒯 h | + | τ ξ 0 ϕ , u h Ω h | .

Finally, the result follows by the discrete trace inequality applied to the boundary terms on the right-hand side, Young’s inequality and the definition of α h and β h . ∎

### Corollary 1.

The HDG scheme (3.6) is well-posed for h sufficiently small.

### Proof.

Let f 0 and ξ 0 = 0 . By (1) we obtain that 𝒒 h = 𝟎 . Moreover, since τ > 0 , we have that u h = 0 on the boundary Ω h and u h = { { u h } } on 𝒯 h ; therefore u h is continuous. These facts, together with (3.6b) lead to

0 = - ( u h , 𝒗 ) 𝒯 h + [ [ 𝒗 ] ] , u h h = ( u h , 𝒗 ) for all  𝒗 𝑽 h .

Thus, taking 𝒗 = u h we conclude that u h = 0 since it vanishes at the boundary. ∎

The energy estimate in Lemma 1 provides the stability bound for the vector-valued unknown 𝒒 h . On the other hand, the stability for the scalar approximation u h can be obtained by a duality argument that we omit since it is not need it for the analysis of the coupled problem. We refer the reader to the proof of [5, Lemma 3.5] or the proof of [33, Theorem 3.1] for details regarding the duality argument employed in this type of unfitted HDG methods. Therefore, it is possible to conclude that there is a constant C HDG > 0 , independent of h, such that

(3.15) J ( 𝐪 h , u h ) + u h Ω h C HDG ( f Ω h + 𝜿 1 2 h - 1 2 ξ 0 ϕ Ω h ) ,

where, for convenience of notation of the forthcoming analysis, we have denoted

(3.16) J ( 𝐪 h , u h ) := ( 𝜿 - 1 2 𝒒 h Ω h 2 + τ 1 2 u h Ω h 2 + τ 1 2 ( u h - { { u h } } ) 𝒯 h Ω h 2 + τ - 1 2 [ [ 𝒒 h ] ] h 2 ) 1 2 .

Having established the well posedness of the discrete formulation, in the following section we will study the behavior of the discretization error.

### 3.4 A Priori Error Analysis

To establish a priori error bounds for the HDG discretization, we will make use of a tool introduced by Francisco-Javier Sayas, Jay Gopalakrishnan and Bernardo Cockburn in [4]. The idea is to use a projection, known as the HDG projection, to decompose the discretization into a component involving the approximation properties of the discrete spaces 𝑽 h and W h , and another component involving the error introduced by projecting into these spaces. The HDG projection over 𝐕 h × W h , denoted by 𝚷 ( 𝐪 , u ) := ( 𝚷 v 𝐪 , Π w u ) , is the unique element-wise solution pair of

(3.17)

(3.17a) ( 𝚷 v 𝐪 , 𝐯 ) T = ( 𝐪 , 𝐯 ) T for all  𝐯 [ k - 1 ( T ) ] e ,
(3.17b) ( Π w u , w ) T = ( u , w ) T for all  w k - 1 ( T ) ,
(3.17c) 𝚷 v 𝐪 𝐧 + τ Π w u , μ e = 𝐪 𝐧 + τ u , μ e for all  μ k ( e ) ,

for every element T 𝒯 h , and e T . The approximation properties of 𝚷 are stated in Section 6. Using this projection we can then define

𝜺 𝒒 := 𝚷 𝑽 𝒒 - 𝒒 h , ε u := Π W u - u h and 𝑰 𝒒 := 𝒒 - 𝚷 𝑽 𝒒 , I u := u - Π W u ,

where 𝚷 𝑽 is the HDG projector onto 𝐕 h , and Π W is the HDG projector onto W h . The terms 𝜺 𝒒 and ε u are known as the projections of the errors and the terms 𝑰 𝒒 and I u are the errors of the projections. The full discretization error can then be split as

𝒒 - 𝒒 h = 𝜺 𝒒 + 𝑰 𝒒 and u - u h = ε u + I u .

We will now show that the scheme (3.6) is consistent and the discretization error is driven solely by the approximation properties of the discrete spaces, as encoded by 𝑰 𝒒 , and I u . We start by noting that from (3.6a) and the decompositions above, it follows that

(3.18) 𝒜 h ( 𝒒 - 𝜺 𝒒 - 𝑰 𝒒 , 𝒗 ) + 𝒜 T ( 𝒒 - 𝜺 𝒒 - 𝑰 𝒒 , 𝒗 ) + h ( 𝒗 , u - ε u - I u ) = 1 , h ( 𝒗 ) for all  𝒗 𝑽 h .

However, since 𝒒 and u satisfy (2.10) in a distributional sense, we have that 𝒒 𝑯 ( 𝐝𝐢𝐯 ; Ω h ) and therefore [ [ 𝒒 ] ] = 0 in h . This also implies that u H 1 ( Ω h ) since u = - 𝜿 - 1 𝒒 L 2 ( Ω h ) . Hence,

𝒜 h ( 𝒒 , 𝒗 ) + 𝒜 T ( 𝒒 , 𝒗 ) + h ( 𝒗 , u ) - 1 , h ( 𝒗 ) = ( 𝜿 - 1 𝒒 , 𝒗 ) 𝒯 h + e Ω h e ( 0 l ( 𝒙 ) 𝜿 - 1 𝒒 ( 𝒙 + 𝒕 ( 𝒙 ) s ) 𝒕 ( 𝒙 ) ) 𝒗 ( 𝒙 ) 𝒏 h 𝑑 s 𝑑 S 𝒙
+ [ [ 𝒗 ] ] , { { u h } } h - ( u , 𝒗 ) 𝒯 h + ξ 0 ϕ , 𝒗 𝒏 h Ω h
= ( 𝜿 - 1 𝒒 , 𝒗 ) 𝒯 h + [ [ 𝒗 ] ] , { { u } } h - ( u , 𝒗 ) 𝒯 h + u , 𝒗 𝒏 h Ω h ,

where in the last equality we have used the fact that 𝒒 satisfies the transfer equation (3.3) and u satisfies (3.2c). Then, by integrating by parts and considering equation (3.2b), we obtain that

𝒜 h ( 𝒒 , 𝒗 ) + 𝒜 T ( 𝒒 , 𝒗 ) + h ( 𝒗 , u ) - 1 , h ( 𝒗 ) = [ [ 𝒗 ] ] , { { u } } h - u , 𝒗 𝒏 h 𝒯 h + u , 𝒗 𝒏 h Ω h = 0 .

Analogously, from (3.6b) we have

(3.19) T ( 𝒒 - 𝜺 𝒒 - 𝑰 𝒒 , w ) + h ( 𝒒 - 𝜺 𝒒 - 𝑰 𝒒 , w ) - 𝒞 h ( u - ε u - I u , w ) = 2 , h ( w ) .

Analyzing the terms above that involve 𝒒 𝑯 ( 𝐝𝐢𝐯 ; Ω h ) and u H 1 ( Ω h ) , and using again the facts that 𝒒 satisfies the transfer equation (3.3) and u satisfies (3.2c), it is easy to verify that

T ( 𝒒 , w ) + h ( 𝒒 , w ) - 𝒞 h ( u , w ) - 2 , h ( w ) = - τ u , w 𝒯 h + 2 τ u , { { w } } h + τ u , w Ω h = 0 .

Putting these arguments together it follows from (3.18) and (3.19) that the scheme is consistent and the following error equations for ( 𝜺 𝒒 , ε u )