A Trefftz Discontinuous Galerkin Method for Time Harmonic Waves with a Generalized Impedance Boundary Condition

A Trefftz Discontinuous Galerkin Method for Time Harmonic Waves with a Generalized Impedance Boundary Condition

Shelvean Kapita Corresponding author. Current address: Department of Mathematics, University of Georgia, 220 DW Brooks Dr, Athens, GA 30602. Email: shelvean.kapita@uga.edu    Peter Monk and Virginia Selgas
Institute for Mathematics and its Applications, College of Science and Engineering, 306 Lind Hall, Minneapolis, MN USA 55455; Department of Mathematical Sciences, University of Delaware, Newark DE 19716, USA; Departamento de Matemáticas, Universidad de Oviedo, EPIG, 33203 Gijón, Spain

We show how a Trefftz Discontinuous Galerkin (TDG) method for the displacement form of the Helmholtz equation can be used to approximate problems having a generalized impedance boundary condition (GIBC) involving surface derivatives of the solution. Such boundary conditions arise naturally when modeling scattering from a scatterer with a thin coating. The thin coating can then be approximated by a GIBC. A second place GIBCs arise is as higher order absorbing boundary conditions. This paper also covers both cases. Because the TDG scheme has discontinuous elements, we propose to couple it to a surface discretization of the GIBC using continuous finite elements. We prove convergence of the resulting scheme and demonstrate it with two numerical examples.

elmholtz equation, Trefftz, Discontinuous Galerkin, Generalized Impedance Boundary Condition, Error estimate, Artificial Boundary Condition


1 Introduction

The Trefftz method, in which a linear combination of simple solutions of the underlying partial differential equation on the whole solution domain are used to approximate the solution of the desired problem, dates back to the 1926 paper of Trefftz [32]. A historical discussion in relation to Ritz and Galerkin methods can be found in  [17]. From our point of view, a key paper in this area is that of Cessenat and Déspres [9] who analyzed the use of a local Trefftz space on a finite element grid to approximate the solution of the Helmholtz equation [9]. This was later shown to be a special case of the Trefftz Discontinuous Galerkin (TDG) method [8, 19] which opened the way for a more general error analysis. For more recent work in which boundary integral operators are used to contruct the Trefftz space, see for example [3, 23]. The aforementioned work all concerns the standard pressure field formulation of acoustics which results in a scalar Helmholtz equation. Indeed, TDG methods are well developed for the Helmholtz, Maxwell and Navier equations with standard boundary conditions and a recent survey can be found in [24]. For the displacement form, a TDG method has been proposed by Gabard [15] also using simple boundary conditions.

Because of the unusual boundary conditions considered in this paper, we propose to use the displacement form of the Trefftz Discontinuous Galerkin (TDG) method for approximating solutions of the Helmholtz equation governing scattering of an acoustic wave (or suitably polarized electromagnetic wave) by a bounded object. This is because the scatterer is assumed to be modeled by a Generalized Impedance Boundary condition (GIBC). These boundary conditions arise as approximate asymptotic models of thin coatings or gratings ([12, 5, 4, 33]). Importantly, they also arise as approximate absorbing boundary conditions (ABCs) and our paper shows how to handle these boundary conditions. As far as we are aware, the displacement TDG method has not been analyzed to date. We provide such an analysis and this is a one of the contributions of our paper.

In order to define the problem under consideration more precisely, let denote the region occupied by the scatterer. We assume that is an open bounded domain with connected complement having a smooth boundary . Then we can define to be the surface gradient and to be the surface divergence on (see for example [11]). In addition denotes the outward unit normal on .

Let , , denote the wave number of the field, and suppose that a given incident field impinges on the scatterer. We want to approximate the scattered field that is the solution of


The last equation, the Sommerfeld radiation condition (SRC), holds uniformly in . In addition

where is the given incident field and is assumed to be a smooth solution of the Helmholtz equation in a neighborhood of . For example, if is a plane wave then , where is the direction of propagation of the plane wave and . Alternatively could be the field due to a point source in . The coefficient functions and are used to model the thin coating on and we shall give details of the assumptions on these coefficients in the next section.

As can be seen from the second equation in (1), the GIBC involves a non-homogeneous second order partial differential equation on the boundary of the scatterer, and this complicates the implementation using a TDG method which uses discontinuous local solutions of the homogeneous equation element by element. In addition, because the problem is posed on an infinite domain we need to truncate the domain to apply the TDG method, and then apply a suitable artificial boundary condition (ABC) on the outer boundary. Because TDG methods have discontinuous basis functions, when the GIBC or ABC involve derivatives, these boundary conditions can be applied more easily if we convert them to displacement based equations, so we propose to solve (1) by converting it to a vector problem. To this end, we introduce , in which case . Using this relationship we see that v should satisfy


where the radiation condition (last equation) holds uniformly for all in directions .

The use of the displacement variable for the Helmholtz equation with standard boundary conditions in the context of plane wave methods was considered by Gabard in [15], but no error estimates were proved. In particular he used the PUFEM [2] and DEM [14] approaches, not TDG. The use of the displacement vector as the primary variable is often necessary in studies of fluid-structure interaction (see e.g. [34]). To date, no error estimates have been proved for the displacement based formulation with or without the GIBC. The vector formulation is useful in its own right. For example, using finite element methods, Brenner et al. [6] show that a vector formulation can also be advantageous for sign changing materials, although we do not consider that problem here.

Our approach to discretizing (2) is to use TDG in a bounded subdomain of , and standard finite elements or trigonometric polynomial based methods to discretize the GIBC on the boundary. The domain is truncated using the Neumann-to-Dirichlet (NtD) map on an artificial boundary that is taken to be a circle. Other truncation conditions could be used. Since it is not the focus of the paper, we assume for simplicity that the NtD map is computed exactly. The discretization of the NtD map could be analyzed using the techniques from [25, 26], and it is also possible to use an integral equation approach to approximate the NtD on a more general artificial boundary but this remains to be analyzed.

Our analysis of the discrete problem follows the pattern of the analysis of finite element methods for approximating the standard problem of scattering by an impenetrable scatterer using the Dirichlet-to-Neumann boundary condition from [29]. We first show that the GIBC can be discretized leaving the displacement equation continuous. Then we show that this semi-discrete problem can also be discretized successfully. The analysis of the error in the TDG part of the problem is motivated by the analysis of TDG for Maxwell’s equations in [22] and uses the Helmholtz decomposition of the vector field satisfying (2) as a critical tool.

The contributions of this paper are 1) a first application and analysis of TDG to the displacement Helmholtz problem; 2) a method for incorporating a discretization of the GIBC into the TDG scheme using novel numerical fluxes from [25]; 3) an error analysis of the fully discrete problem (except for the NtD map as described earlier), and the first numerical results for TDG applied to this problem.

In the remainder of the paper we use bold font to represent vector fields and we will work in . We utilize the usual gradient and divergence operators (both in the domain and on the boundary), and also a vector and scalar curl defined by

for any and .

The paper proceeds as follows. In the next section we formulate problem (2) in a variational way and show it is well posed using the theory of Buffa [7]. Then in Section 3 we describe and analyze the discretization of the GIBC using finite elements (or trigonometric basis functions). The fully discrete TDG scheme is described in Section 4 where we also prove a basic error estimate and show well-posedness of the fully discrete problem. We then prove convergence in a special mesh independent norm. In Section 5 we provide a preliminary numerical test of the algorithm, and in Section 6 we draw some conclusions.

2 Variational Formulation of the Displacement Method

In this section we give details of our assumptions on the coefficients in the GIBC, and formulate the displacement problem (2) in variational setting suitable for analysis. Then we show that the problem is well-posed. The functions in (1) are complex valued functions and we assume that there exists a constant such that


Of key importance will be the operator defined weakly as the solution operator for the boundary condition on relating the Neumann and Dirichlet boundary data there. More precisely, for each we define to be the solution of


An essential assumption is the following.

Assumption 1.

The only solution of

is .

We will show that Assumption 1 together with the conditions (3) ensure that the operator is well-defined.

Remark 1.

One possible condition under which Assumption 1 holds is,

Remark 2.

On the one hand, the assumptions in (3) concerning the imaginary parts of and are governed by physics, since these quantities represent absorption when our model is deduced as an approximation of the Engquist-Nédélec condition modeling the diffraction of a time-harmonic electromagnetic wave by a perfectly conducting object covered by a thin dielectric layer (see [12]).

On the other hand, the hypothesis in (3) on the real part of is technical and ensures ellipticity (see [4, Th.2.1]); however, this property is fulfilled in the example of a medium with a thin coating (see [12]). It would also be possible to allow on as might be encountered modeling meta-materials, but a sign changing coefficient would require a more elaborate study.

The role of these properties will be clarified in Lemma 2.3.

The assumptions on the coefficients in (3) together with Assumption 1 ensure that problem (1) has a unique weak solution in the space (see later and [4, Th.2.1]).

To solve (2) we first truncate the domain. We wish to analyze the error introduced in approximating a scattering problem, concentrating on the discretization of the GIBC, so we truncate the domain using a simple analytic Neumann-to-Dirichlet map. Obviously other more general truncation approaches such as integral equations could be used. Indeed, in the numerical section, we shall consider a GIBC that arises from approximating the Neumann-to-Dirichlet map (a higher order ABC).

Let denote the ball of radius centered at the origin and set be our computational domain (i.e. the bounded domain that we will mesh for the UWVF) and , where the radius is taken large enough to enclose (see Fig. 1 for a diagram illustrating the major geometric elements of the problem). The following Neumann-to-Dirichlet (NtD) map will provide the ABC on . In particular let solve the exterior problem


for some , then . Let us recall that is an isomorphism since its inverse, the Dirichlet-to-Neumann map, is also an isomorphism [11]. Obviously, the solution of (1) satisfies and, in consequence, using the fact that ,

where we denote by the outward unit normal on . In the same way, for the solution of (1) we have that

Now we can write down a weak form for the boundary value problem (2) in the usual way, multiplying the first equation in (2) by a test vector function and integrating by parts:

where the minus sign in the last term is due to the normal field pointing outward . Using the NtD map and the boundary solution map , the above equation can be rewritten as the problem of finding such that


for any . It will be convenient to associate with the left hand side of (6) the sesquilinear form defined by

Figure 1: A cartoon showing the geometric features of the problem. The bounded scatterer is covered by a thin coating giving rise to a GIBC on . An incident wave on this scatterer causes a scattered field in the exterior of . The artificial boundary is introduced to truncate the domain resulting in a bounded computational domain and is taken to be a circle for simplicity.

In order to prove the well-posedness of this variational formulation, we now summarize some of the properties of the NtD map and GIBC boundary map . For a given function the NtD map is given by


where are the Fourier coefficients of on and

According to [10, page 97] there are constants and such that

for all . Now define by

Clearly is negative definite and

Also from [10, page 97] we can obtain the asymptotic estimate

so that when . Hence,


where is well-defined and bounded, in particular is compact.

We next state some properties of the NtD map which follow from the properties of the better known DtN map.

Lemma 2.1.

For all , it holds


The first inequality follows from [29, Lemma 3.2], whereas the second is proved as follows: For any , we may write

where, as above,

are the Fourier coefficients of on . Since , taking the imaginary part

by the Wronskian formula for Bessel functions (see e.g. [1, 9.1.16]). ∎

We note that the foregoing theory provides a direct proof that is an isomorphism as a consequence of the Fredholm alternative thanks to Lemma 2.1 and the splitting (9).

Corollary 2.2.

The operator is an isomorphism.

Next we show that is well defined.

Lemma 2.3.

Under Assumption 1 and the conditions (3), the operator defined in (4) is an isomorphism. In particular, is well-defined, linear and continuous.


We start the proof by defining the bounded sesquilinear forms by

for any . Thanks to the Riesz representation theorem, we can consider the associated operators that satisfy

for any . Notice that, under assumption (3),

and, in consequence, using the Lax-Milgram theorem guarantees that is an isomorphism.

Also notice that, by Rellich’s theorem we know that is compactly embedded into , so that is compact.

Moreover, under Assumption 1, is injective. Therefore, by the Fredholm alternative, is an isomorphism. ∎

The next lemma shows that the impedance boundary condition does not cause a loss of uniqueness for the scattering problem.

Lemma 2.4.

For any , it holds that


Using the variational definition of for , and choosing the test function , gives

This implies

The last inequality follows from the assumptions and in (3). ∎

Starting our analysis of (6) we show that any solution is unique.

Lemma 2.5.

Problem (6) has at most one solution.


Let us consider any solution of its homogeneous counterpart, that is, such that


for all . Since is connected, by the Helmholtz decomposition theorem (see [18, Th.2.7-Ch.I]) we can rewrite v as

for some and , where

Then, the homogeneous problem (10) may be rewritten as


for all . In particular, taking we deduce that

Noticing that leads to

Hence, by uniqueness of the solution (up to a constant) of the interior Neumann problem for Laplace operator in , we have that in for some constant ; in particular, in . Furthermore, satisfies

so that in and we deduce from (11) that

In consequence, by the invertibility of and and the uniqueness of solution of the forward problem with GIBC (see [4, Th.2.1]), we have that in ; that is to say, in . Summing up, we conclude that

Using this uniqueness result and a suitable stable splitting of , we will be able to apply [7, Theorem 1.2] to prove the well-posedness of the continuous problem. In particular, we write



and is endowed with the inner product

Notice that the orthogonality of the above splitting implies that if, and only if, and for all . We also need to define the duality pairing between and its dual space , with respect to the pivot space , so that (note: this is defined without conjugation):

According to the above splitting, any has the form for some and . By the orthogonality of the splitting, and the fact that , we have that

and, in particular, the splitting is stable. Moreover, it allows us to define the linear continuous operator by .

Next we define such that if then is given via the Riesz representation theorem by

recalling the definition of in (7).

We can now state and prove the following result.

Theorem 2.6.

Problem (6) is well-posed and the Babuška-Brezzi condition is satisfied.


Let be split into for some and , and similarly . Then


We expand the troublesome term

So we can define the sesquilinear form

and use the remaining terms in (12) to define the sesquilinear form

On the one hand, since is negative definite (see Lemma 2.1) and the splitting of the space is stable, we have that there is a constant independent of such that

Now define the operator by

Notice that is compact because each sesquilinear form in its definition is compact. For example, the sesquilinear form

is compact by [27, Theorem 1.3], because the trace of functions in defined into is a compact operator; indeed, is a subset of due to our assumption of a smooth boundary , and the normal derivative operator is compact from into . The remaining sesquilinear forms are also compact by the same reasoning. Hence is compact. Then we conclude that

Hence all the conditions of [7, Assumption 1] are satisfied and the existence of a unique solution to (6) is shown by [7, Theorem 2.1]. In addition this theorem shows that there is an isomorphism such that

This in turn implies that the Babuška-Brezzi condition is satisfied. ∎

3 A Semidiscrete Problem

In this section we consider a semidiscrete problem in which the GIBC boundary operator is discretized but the space where we search for the solution in is not. As discussed in the introduction, we shall not consider the truncation of the NtD map here.

We shall need an additional assumption on the boundary operator . In particular we need to know that it smooths the solution on the boundary , so we make the following second assumption.

Assumption 2.

For each , it holds that and there exists such that for any .

Remark 3.

Note that if Assumption 2 holds, since , it also holds for .

Notice that Assumption 2 further constrains the choice of the coefficients and in the generalized impedance boundary condition on . Its role will be clarified in Lemma 3.1, where we apply Schatz’s analysis [31] in order to show that the finite element approximation of , defined shortly, converges.

On the inner boundary we consider a finite dimensional subspace of continuous piecewise polynomials of degree at least (with ) on a mesh . We assume that the mesh consists of segments of the boundary of maximum length , and that it is regular and quasi-uniform: the latter means that there exists a constant such that

where denotes the arc length of the edge