Coupling of HDG with a double-layer potential BEM
In this paper we propose and analyze a new coupling procedure for the Hybridizable Discontinuous Galerkin Method with Galerkin Boundary Element Methods based on a double layer potential representation of the exterior component of the solution of a transmission problem. We show a discrete uniform coercivity estimate for the non-symmetric bilinear form and prove optimal convergence estimates for all the variables, as well as superconvergence for some of the discrete fields. Some numerical experiments support the theoretical findings.
AMS Subject classification. 65N30, 65N38, 65N12, 65N15
In this paper we propose and analyze a new coupling procedure for the Hybridizable Discontinuous Galerkin Method (HDG)  and a Galerkin Boundary Element procedure based on a double layer potential representation. The model problem is a transmission problem in free space, coupling a linear diffusion equation with variable diffusivity () in a polygonal domain, with an exterior Laplace equation. The transmission conditions are given by imposing the value of the difference of the unknown and its flux on the interface between the interior and the exterior domains.
Much has been written on the advantages and disadvantages of the many Discontinuous Galerkin schemes. In support of this piece of work, let us briefly hint at some of the features that make the HDG –originally a derivation of the Locally Discontinuous Galerkin (LDG) method by Bernardo Cockburn and his collaborators– a family of interest. First of all, HDG uses polynomials of the same degree for discretization of all the variables –in the case of diffusion problems, a scalar unknown , the vector valued flux and a scalar unkown on the skeleton of the triangulation–, attaining optimal order in the approximation of all of them. This feature compares well with Mixed Finite Element Methods, of which HDG is a natural modification: HDG can be considered as a variant of the Raviart-Thomas and Brezzi-Douglas-Marini Mixed Elements, implemented with Lagrange multipliers on the interelement faces , eliminating the degrees of freedom used for stabilization and using a stabilization (not penalization) parameter to perform the same task. Similar to the Lagrange multiplier implementation of mixed methods by Arnold and Brezzi , HDG is implemented by hybridization, reducing the unknowns to those living on interelement faces. This reduces considerably the number of degrees of freedom and makes the method competitive with respect to other families of Finite Element Methods, while still being a mixed method, that approximates several fields simultaneously. Another interesting feature of HDG is the fact that for polynomial degrees , the two unknowns related to the scalar field (the respective approximations inside the element and on the skeleton) superconverge at rate . This allows for the application of standard postprocessing techniques  that can be traced back to Stenberg .
Among many of its good properties (for the limited set of equations where it is usable), the Boundary Element Method provides a reliable form of constructing high order absorbing boundary conditions with complete flexibility on the geometric structure of the domain (it does not require the domain to be convex or even connected). It is therefore natural to test and study the possibility of using BEM as a way of generating a coupled discretization scheme with HDG for transmission problems. The coupling of DG and BEM started with the study of LDG-BEM schemes , extended to other DG methods of the Interior Penalty (IP) family . All of these methods, plus several new ones, were presented as particular cases of a methodology for creating symmetric couplings of DG and BEM in . One of the methods in this latter paper, using HDG and BEM, was recently analyzed in . As opposed to symmetric couplings with BEM, that need two integral equations (and thus four integral operators) to reach a stable formulation, non-symmetric couplings require only one integral equation (and two integral operators). From this point of view, they provide simpler (and more natural) forms of coupling BEM with field methods (FEM, Mixed FEM or DG). The non-symmetric coupling of Mixed FEM and BEM is recent  and expands ideas used for analyzing the simple coupling of FEM-BEM . The presentation and testing of non-symmetric coupling of DG methods of the IP family with BEM , led to its analysis , providing the first rigorous proof of convergence of this kind of schemes. Note that, unlike in the case of symmetric couplings, where the stability analysis boils down to an energy argument, coercivity properties in non-symemtric couplings pose a serious analytic challenge. The present paper makes a contribution in that direction, showing how to develop a coercivity analysis for non-symmetric coupling of HDG-BEM, and pointing out some interesting facts about the size of the diffusion parameter in the interior domain.
The paper is structured as follows. In Section 2 we present the method. Sections Section 3 to Section 5 cover the analysis in progressive steps. Section 6 discusses several related extensions and possible modifications of the method. In Section 7 we include some numerical experiments. A collection of known results concerning transmission problems, potentials and integral operators, is given in Appendix A for ease of reference.
Regarding notation. Basic theory of Sobolev spaces is assumed throughout. Norms of type will be subscripted with the integration domain . All other norms will be indexed with the name of the space.
2Formulation and discretization
Let be a polygonal domain when or a Lipschitz polyhedral domain when , with boundary . For simplicity we will assume that is connected. (Connectedness plays a minor role in the amount of energy-free solutions of the transmission problem.) Let be then the domain exterior to . The outward pointing unit normal vector field on is denoted .
2.1Statement of the problem
We consider a scalar positive diffusion coefficient such that .
In principle, boundary values of and in are taken in the sense of traces for functions with local Sobolev regularity. Similarly and are defined in a weak form as elements of the space . Basic regularity requirements for data is: , , and We note that is valid data for the transmission problem, but not for the type of formulation we are going to use. Data functions are assumed to satisfy the following compatibility condition:
A more detailed discussion about this condition and solvability issues for the transmission problems is given in Section Appendix A.4. Here we just note the following:
When , condition is necessary and sufficient for the transmission problem to have a solution. In this case as .
When , problem has a unique solution even if fails to hold. If is satisfied, the solution decays as at infinity. Proposition ?(b) shows the elementary modification of data that is needed to have condition satisfied, while still providing a solution of the original problem.
2.2Coupling with a double layer potential
Condition allows us to write as a double layer potential
(See and for the definitions and Proposition ? for the representation result.) The exterior trace and normal derivative of the double layer potential can be written in terms of boundary integral operators, as explained in Section Appendix A.3. These operators are introduced in -. They allow us to write the transmission conditions in the equivalent form
At this moment, it is convenient to collect all the equations that make up the formulation we are going to discretize:
The exterior field is reconstructed as a double layer potential . The next lemma ensures that the integral boundary condition and the normalization condition for the density are adequately encoded in the system .
Integrating over and testing with , it follows that
because of Proposition ?. This proves the result.
2.3Discretization with HDG and Galerkin BEM
We start the discretization process by describing the discrete geometric elements. The domain is divided into triangles () or tetrahedra () with the usual conditions for conforming finite element meshes. A general element will be denoted , will be the set of all elements and the maximum diameter of elements of the triangulation. The triangulation is assumed to be shape-regular. Extension of the forthcoming results to non-conforming meshes is relatively simple while not straightforward, requiring the use of some finely tuned results about HDG on general grids , .
We assume the existence of a triangulation of the boundary , formed by line segments () or triangles (). We write . From the point of view of analysis, it is immaterial whether this grid is related to or not. (Note that for the experiments we will only use meshes where is the trace of on , which makes the implementation considerably simpler.) At the time of stating and proving the convergence theorems we will add some technical restrictions relating the mesh-sizes of the two grids.
The set of edges () or faces () of all elements of the triangulation will be denoted and will denote the skeleton of the triangulation. Volumetric integration will be denoted with parentheses
while integration on boundaries will be denoted with angled brackets
The unit normal outward pointing vector field on will be denoted , or simply , when there is no doubt on what the element is.
Local discrete spaces will be composed of polynomials. The set contains all -variate polynomials of degree less than or equal to and . If or , is the set of polynomials of degree less than or equal to defined on , i.e., the space of -variate polynomials on local tangential coordinates. Four global spaces will be used for discretization:
The discretization method consists of using the Hybridizable Discontinuous Galerkin method (see , ) for equations — (this part of the system can be understood as an interior Dirichlet problem) and conforming Galerkin (Boundary Element) Method for equation (this part is considered as a hypersingular integral equation for an exterior Neumann problem). We thus look for , , , and satisfying
where is a non-negative stabilization function that is constant on each edge/face of and such that it is strictly positive on at least one edge/face of each triangle. The stabilization function can be double-valued on interlement edges/faces and is in principle multiple valued, but its normal component is made to be single valued through equation . Note also that equations - are tested with the same space, but because of the integration domain, they produce together as many equations as the dimension of . Using the local solvers for the HDG method , the system can be reduced to a system with as the only unknowns. A discrete counterpart of Lemma ? holds.
Testing equation with and equation and using Proposition ?, the result follows.
The next sections deal with the analysis of the method: we first show the basic discrete coercivity arguments (Section 3), proceed to proving energy estimates based on the HDG projection (Section 4) and finally use duality arguments to analyze convergence (and superconvergence) of some of the fields (Section 5).
3Solvability of the discrete system
The aim of this section is to show that, under some conditions on , the system has a unique solution for small enough. A relevant hypothesis will be the following:
This hypothesis is related to coercivity of the underlying mixed formulation. We will discuss this hypothesis (and compare it with similar bounds for other coupled BEM-FEM schemes) in Section 7, where we will also make some experiments related to it. For some of the forthcoming arugments, we will use the piecewise constant function given by . We will pay attention to the quantities
noting that because of the shape regularity of the grid, we can bound for all . Three more polynomial spaces will appear in our arguments:
where in the last space, we have denoted .
Part (a) is straightforward. Part (b) is Lemma 4.1 of .
Using the Bramble-Hilbert lemma (or a generalized Poincaré inequality), the trace theorem and an argument about finite dimensions, it is easy to prove that
where is the reference element and is the orthogonal projector. The result follows then by a scaling argument.
Condition is equivalent to the fact that the discrete normal fluxes are single valued on internal faces/edges. Therefore
If is the orthogonal projection onto (as in Lemma ?), then
where we have used and the fact that . Therefore, by Lemma ? and Young’s inequality
for arbitrary . Therefore
where we have chosen
This finishes the proof.
We only need to prove that the homogeneous system only admits the trivial solution. Let then be a solution of with zero right hand side. By Lemma ? it follows that . We now test equations with , , , and respectively, add the equations, and simplify the result to obtain
Let now . Using then and Proposition ?, it follows that
where makes reference to the trace of on . We are now in the hypotheses of Lemma ?, which implies that , , and on for all . By Proposition ? and the fact that , this implies that .
Going back to , it follows that
By Lemma ?(b), there exists such that of . Since and on at least one , it follows by Lemma ?(a) that and therefore on . Using this information in and testing with , it follows that is piecewise constant. At the same time, on , which implies that both and are constant. Testing with implies that this constant value has to vanish. This finishes the proof of uniqueness.
4Estimates by energy arguments
Convergence analysis of the coupled HDG-BEM scheme follows the main lines of the projection-based analysis of HDG methods  (see also ). We start by recalling the local HDG projection . Given , we define as the solution of the local problems
for all . We also consider the orthogonal projections and . The error analysis is carried out by comparing the discrete solution with the projection of the exact solution, i.e., in terms of the quantities:
The error in the normal flux (recall ) is
Note also that by definition of the projections
This shows that is the best approximation of on . Using then a local trace inequality, we can easily bound for
The second inequality follows from a similar argument, comparing with instead of .
For simplicity, in some of the forthcoming arguments we will shorten . The following three discrete functionals will be relevant in the sequel as well:
By the definition of the projections and , it follows that
Subtracting from equations we obtain the error equations:
Testing these equations with , , , and respectively, adding them, and simplifying, we obtain
Let now and note that Proposition ? and imply that
Applying now the coercivity estimate of Lemma ?, it follows from that
Since the bilinear form is -coercive (Proposition ?), the result follows.
At this moment we have to introduce a technical hypothesis,
meaning that the meshsize of is controlled by the minimum local meshsize of on . This hypothesis is satisfied, for instance, if is quasiuniform in a neighbordhood of and is equal the restriction of to , or is a refinement of this inherited mesh. We will also assume that
We will derive the proof from Proposition ? and some estimates related to the functionals defined in . First of all
Using Proposition ? and a scaling argument (note the use of in the scaling argument), it follows that
An Aubin-Nitsche duality argument and properties of best approximation by piecewise polynomial functions shows that
For the third functional, we use Proposition ? and
where it is to be understood that the local projection is only applied (and needed) on . A local scaling argument (cf ) shows that
This inequality and Lemma ? can then be used in to prove that
The result is then a consequence of Proposition ? and the bounds , , , and .
Adding and subtracting , we can bound
A scaling argument (using the fact that ), plus Lemma ? and Theorem ? finish the proof.
Let be the faces (with ) or edges (when ) of . We consider the space , endowed with its product norm.
This bound is a direct consequence of Theorem ?, using well-known estimates for the best approximation by piecewise polynomials and .
5Estimates by duality arguments
Consider the problem
By Proposition ? and -, we can represent in . Therefore, using the definition of the operator –see – and Proposition ?, we can write
We will assume the following regularity estimate
Since this is a transmission problem, this hypothesis is mainly related to the regularity of the diffusion coefficient . If is smooth and equal to one in a neighborhood of , then holds (see  for a similar argument).
The gist of the proof consists of testing the error equations with , , , , and respectively, add them and manipulate the result. Let us first test the first error equation with :
Using integration by parts twice, the first condition in the definition of the HDG projection and the fact that , we can write
Adding these two equations, it follows that
where in the last equation we have applied and the fact that and are single valued on interelement faces. We next test with and apply and the fact that :
The following step consists of subtracting this equation from to obtain
Bound for . Since ,
where is the local orthogonal projection on the spaces if and . Note that the inclusion of is possible because of the definition of the HDG projection . Therefore
by regularity and .
Bound for . Let now be the Clément approximation  on a conforming finite element space on a triangulation whose restriction to the boundary is . We then use to write
where . Testing the error equation with , and using, it follows that
where we have used that
which is just after using Lemmas ? and ?. Inserting in and applying Propositions ? and ? we can bound
We now just need to bound all the terms in the right hand side of . A duality argument (Aubin-Nitsche trick), Lemma ?, and the regularity assumption show that
Well-known properties of the Clément interpolant (see in particular ) and prove also that
Finally, we just overestimate
after applying a discrete trace inquality and . Bringing the bounds , and to and using Corollary ?, we obtain the bound
The bound is now a direct consequence of , , and Theorem ?.
Note that in absence of any kind of additional regularity hypothesis, the estimate of Proposition ? can be easily repeated without the additional