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.
The goal of this article is to conclude the work started by Cockburn, Sayas and Solano in the article Coupling at a distance , 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  was published, the technique for transferring information between the two grids had only been recently incorporated into the HDG literature  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  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 , 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 .
The method proposed in , 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 , 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 that has a smooth parametrizable boundary that will be denoted by . 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:
The function f will be taken to be compactly supported and square integrable on . The diffusion coefficient is a strictly positive matrix-valued function such that, denoting the identity matrix is as , the difference is compactly supported in . This condition implies that outside of 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 of it holds that
The Dirichlet boundary data will be considered to be an element of the trace space . The radiation condition at infinity (2.1d) is equivalent to assuming that there is a constant such that (see ).
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 , the support of f and the support of . We will also require that . The domain interior to Γ will be denoted Ω, while the unbounded complementary region will be denoted . The boundary of Ω will be denoted as and consists of two disjoint components: the artificial boundary Γ and the original problem boundary , so that . We will denote the unit normal vector to , pointing in the direction of for points in Γ and in the direction of for points in 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 from the system. We will represent the solutions to (2.1) as the superposition
where the functions u and are supported in Ω, while is supported in . The pair satisfies the interior problem
On the other hand, the exterior function satisfies
Above, the boundary value corresponds to the trace of over the artificial boundary Γ, while 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  for a more concise treatment.
We start by introducing the single layer and double layer potentials defined respectively for , and as
where 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 in as
and extend it by zero for . The constant captures the far field behavior of the function and will have to be determined. Since 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
To ensure that 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 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 , depending only on Γ and the norms of and , such that
Moreover, from this estimate and the representation formula (2.6), it follows that there exists 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
Above, the source term and the Dirichlet boundary data is given by
where and denote the -inner products over Ω and , respectively. From the three preceding equations, we arrive at the variational problem:
Find such that
where the bilinear forms , , and the functionals and are defined by
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 that approximate Ω in the sense that the Lebesgue measure , as . We will refer to any such as a computational domain and will triangulate by a shape-regular triangulation as depicted in Figure 1. A generic element in will be denoted by T and the mesh parameter h will be defined as diameter of a circle inscribing an element . The set , will be referred to as the skeleton of the triangulation. The set of edges, e, of will be denoted by and we will distinguish between those edges lying entirely in the computational boundary
and those that are either interior or have at most their endpoints in the computational boundary
We will refer to the former as boundary edges and to the latter as interior edges. Note that .
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 , where
We will require that the computational domain and the triangulation satisfy the following local proximity condition: for any point in the computational boundary , the minimum distance between and the boundary should be, at most, of the same order of magnitude as the diameter of the smallest triangle , such that . 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 . Instead, as the mesh diameter , the process involves the passage through a sequence of pairs domain/triangulation that satisfy the local proximity condition and exhaust the original domain Ω as the refinement progresses. We refer the reader to , where this condition is discoursed in more detail, and to  where an algorithm for building a sequence is described.
Mesh-Dependent Subspaces and Inner Products.
For the discrete formulation we will have introduce the following mesh-dependent inner products:
These inner products induce mesh-dependent norms that will be denoted, respectively, by
The finite-dimensional discontinuous polynomial subspaces that will be used for discretization, for , are given by
where denotes the space of polynomials of degree at most k defined in . Similarly, denotes the space of polynomials of degree at most k defined over a edge .
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 between the boundary and the computational boundary. To this purpose, we will tessellate this region as follows:
Let and be the endpoints of a boundary edge .
Let and be the corresponding points in – as determined by the mapping (3.1).
Let and be the straight segments connecting to and to .
We will refer to the open region of delimited by e, and and the segment of connecting to as an extension patch and will denote it by . It is clear that for every there is one and only one such (this justifies subindex in the notation) and that
It also follows from this construction that for every there is only one element in the triangulation such that . We will use this fact to define an extrapolation operator that will extend the value of the piecewise polynomial functions defined on onto the corresponding extension patch , 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 by extrapolating the values of the corresponding polynomial from , and will denote its as for any .
For a given domain and corresponding triangulation , 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, in the usual manner, and extend the definition to , 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 as , which will coincide with the exterior normal on element edges belonging to the computational boundary .
Finally, for every edge we will denote the ratio between its distance to the boundary and the diameter, , of its parent element as
and will define the boundary proximity parameter as
and will assume for this work that the family of admissible domains and triangulations is such that:
as , where the normal 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 satisfying the admissibility requirements discussed in the previous subsection. In addition, we will need to define a bijective mapping
assigning a point to every point .
For any fixed computational domain , the solution pair to (2.11) satisfies the related problem
where the boundary condition can be calculated by integrating equation (2.10b) along a path connecting to . More precisely, if we denote the distance between and by , and by the unit vector , the boundary conditions on 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 , , and the functionals and 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 , introduced by the transfer of boundary condition. The well posedness of problems of this form was established in . 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 , seeks an approximation satisfying
where E denotes the extrapolation operator. The numerical flux in the normal direction is defined as
where τ stabilization function. Throughout this analysis we will only require , where denotes the maximum value of τ.
In order to apply known results from functional analysis, we rewrite the numerical trace in terms of averages and jumps. For this, we use the equation (3.4c) and separate the term featuring as
Above, we have used the fact that the hybrid variable 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
where the bilinear forms , , , and the functionals and are defined by
The unique solvability of the scheme (3.6) will be proved by an energy argument. To that end, for and , it is convenient to define the following norm on the extension patch :
This norm is equivalent to the standard -norm as shown first in  for the two-dimensional and later extended to three dimensions in . That is, there exist positive constants and , independent of h, such that
We also introduce the element-wise constants
where . These constants are independent of h, but depend on the polynomial degree k and the mesh regularity parameter as shown in .
We now proceed to derive an energy inequality that will lead to the well-posedness of (3.6).
Let and . It holds
By taking and in (3.6), and subtracting the resulting expressions we obtain
First of all, after performing algebraic calculations, we observe that is a semi-definite operator from to . In fact,
We will now obtain a lower bound for the non-positive terms of left-hand side of (3.11). In this direction, the operator can be bounded as follows. Let and . By the Cauchy–Schwarz inequality and the definition in (3.9),
where we have used the bound . Then, by the discrete trace inequality, we have
The same arguments yield to
Therefore, combining the above estimates and (3.11), we deduce that
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 and . ∎
The HDG scheme (3.6) is well-posed for h sufficiently small.
Thus, taking we conclude that since it vanishes at the boundary. ∎
The energy estimate in Lemma 1 provides the stability bound for the vector-valued unknown . On the other hand, the stability for the scalar approximation 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 , independent of h, such that
where, for convenience of notation of the forthcoming analysis, we have denoted
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 . 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 and , and another component involving the error introduced by projecting into these spaces. The HDG projection over , denoted by , is the unique element-wise solution pair of
for every element , and . The approximation properties of are stated in Section 6. Using this projection we can then define
where is the HDG projector onto , and is the HDG projector onto . The terms and are known as the projections of the errors and the terms and are the errors of the projections. The full discretization error can then be split as
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 . We start by noting that from (3.6a) and the decompositions above, it follows that
However, since and u satisfy (2.10) in a distributional sense, we have that and therefore in . This also implies that since . Hence,
Analogously, from (3.6b) we have