Revisiting generalized FEM: a Petrov–Galerkin enrichment based FEM interpolation for Helmholtz problem

A concept of Petrov–Galerkin enrichment which is appropriate for highly accurate and stable interpolation of variational solutions is introduced. In the finite element context, the setting refers to standard trial functions for the solution, while the test space will be enriched. The FEM interpolation procedure that we propose will be justified by local wavelets with vanishing moments based on Gegenbauer polynomials. For the reference Helmholtz equation, the continuous piecewise polynomial test functions are enriched using dispersion analysis on uniform meshes in 2d and 3d. From a-priori and a-posteriori numerical analysis it follows that the Petrov–Galerkin based enrichment approximates the exact interpolate solution of the Helmholtz equation with at least seventh order of accuracy.


Introduction
We investigate an interpolation problem for variational solutions obtained from finite element approximations. The specialty consists in the fact that the discrete solution should be close to the exact solution at a finite number of trial points, which are selected for topology optimization purpose. As a typical example for such topology optimization problems see [7,18,19,21].
It will suffice to approximate solutions more accurately at the selected trial points (which associate mesh nodes), rather than in the whole computational domain. To reduce the discretization error at the nodes it is necessary to address the interpolate solution. Motivated by the special need of highly accurate and stable interpolation for a class of variational solutions of scattering problems, we develop a concept of Petrov-Galerkin enrichment by the FEM.
Our motivation of the FEM interpolation come from numerous applications in the engineering sciences where topology optimization problems arise. For typical formulations of topology optimization as well as related inverse and ill-posed problems we refer to [2,10,17,20], and to [8,12,16,23,31] for common numerical approaches used in the field.
As the reference model we consider the Helmholtz equation. It is well known that generalized finite element methods (GFEM) are well suited in this case since they reduce the so-called pollution error of discretization. We refer to [13,15,26,29,33] for the theoretical background of GFEM, and to [34] for their review. Major developments for GFEM were carried out by Babuška with coauthors.
To reduce the discretization error it is possible to enrich either the trial or the test function spaces. The former approach enriching trial spaces is commonly used, but it has high computational costs and meets primarily the task of approximation rather than interpolation. The latter approach enriching test spaces is based on a Petrov-Galerkin setting, see e.g. [4,11,28]. The Petrov-Galerkin concept is also the basis for the variational multiscale method (VMS) by Hughes et al., e.g. [14].
We shall utilize necessary optimality conditions for interpolation properties which are deduced when enriching the space of test functions. Using a Petrov-Galerkin approach, we suggest low order interpolation polynomials for the trial space, and we enrich the test space with high order shape functions. The resulting Petrov-Galerkin enrichment (PGE) improves significantly the accuracy of interpolation, which is at a low price because low order finite elements are used for the discrete solution.
In the present paper a theoretical justification of PGE is provided by local wavelets with vanishing moments based on Gegenbauer polynomial approximation. We derive also practical formulas for calculation of the system matrix for the reference Helmholtz equation given over uniform meshes in 2d and 3d.
For possible extension to nonuniform meshes and isoparametric elements we refer to [25]. In the above paper, a global FE basis is enriched with polynomial bubbles based on preprocessing, where the weights multiplying these bubbles have to be found by solving local problems minimizing the dispersion. In the present context, the shape functions are to be weighted within the local optimization.

The concept of Petrov-Galerkin enrichment
We consider a class of variational problems formulated in the following form. Let X 0 be a closed subspace of a complex Hilbert space X . For a given u 0 ∈ X find u − u 0 ∈ X 0 such that b(u, v) = 0 for all v ∈ X 0 , (2.1) where the map b : X × X → C is continuous and sesquilinear over C (respectively, bilinear over R).
Assuming the existence of a variational solution to (2.1) we look for its finite dimensional approximation in a trial space X h ⊂ X . Let X 0 h be a closed subspace of X h and u h 0 ∈ X h . For finite elements, the index h ∈ R + refers to the mesh size. Following the Petrov-Galerkin approach we set a discrete counterpart of (2.1) in the form: over a test space Y h ⊂ X , where the continuous and sesquilinear form b h : X × X → C refers to a discrete version of b. It may result from approximation of a computational domain, or may relate to the modified forms resulting from stabilization methods. In particular, b h = b is possible option. We remark that the solution u h depends on the choice of the test space Y h . The main idea resides in defining the test space Y h in (2.2). The standard Galerkin method relies on setting Y h = X 0 h which is not the best choice. In fact, we consider the following. Let B X := {φ i ∈ X |i = 1, . . . , ∞} form a FE basis of X . In this basis, we denote by X ⊥ h ⊂ X the orthogonal complement to X h with respect to (2.2), namely Noting that in general X = X h X ⊥ h we define the space We exclude the trivial case Z h = ∅. Indeed, in the subsequent considerations we have The following discussion will rely on an interpolation operator I X h : X → X h defined on a suitable subspace X ⊆ X where exact interpolation conditions at a fixed number of points are well-defined. This property is important, for instance, in topology optimization. This requires that the solution u of (2.1) belongs to X such that I X h u is well defined. In our case I X h is the classical interpolation operator and thus it suffices if X is contained in the space of continuous functions.
For any such interpolation operator, the approximation error between the solutions u of (2.1) and u h of (2.2) can be estimated as where u − I X h u X is the interpolation error, and I X h u − u h X is the error of the FE solution with respect to exact interpolation. We note that the former does not depend on the test space Y h , whereas the latter one is dependent on Y h .
Our goal is to construct the test space in such a manner that the latter error is minimized: If the minimum in (2.5) is zero, this implies that If the exact solution u and, hence, its interpolate I X h u are known, then the necessary optimality condition for (2.5) can be deduced directly from (2.6). For this task, we substitute (2.6) in (2.2) to determine the test space Y h in Z h in Algorithm 2.1 below. We assume that the union of all finite element spaces with respect to the family of meshes under consideration is dense in X . Since the FE basis is given in X we suggest the following conceptual algorithm for construction of the interpolation by means of the discrete problem (2.2).
If the test space Y h is constructed such that (2.7) in Step 2 is satisfied, then the discrete solution u h of (2.2) coincides with the exact interpolation solution I X h u. This construction needs knowledge of the exact solution u ∈ X of (2.1), or, at least, particular solutions for specific data u 0 . Otherwise, we suggest to relax (2.7) with the following approximate condition: Then the error of the FE solution of (2.2) can be estimated by We note that from (2.9) the stability estimate u h − u 0 In Sect. 3 we shall present our approach for the specific case of the Helmholtz equation. We confine ourselves to the case when X h is the space of continuous piecewise linear trial functions. For this case, it is known that the approximation error is of order O(h). We realize Algorithm 2.1 by substituting particular solutions in the form of plane waves into (2.8) and by dispersion analysis. In this way we determine the space Y h consisting of weighted continuous piecewise quadratic test functions in Z h . As a result of asymptotic analysis for κ := kh 2 → 0, where k stands for the wave number, we obtain the error function Er(h) = o(κ 7 ) in (2.8). Moreover, we prove that the system matrix associated to b h is positive definite and hence the inf-sup condition (2.9) holds. Therefore Proposition 2.2 guarantees interpolation of order o(κ 7 ).

The Helmholtz equation
In this section we first formulate the Helmholtz problem under consideration and then show how to use Algorithm 2.1.

The Helmholtz problem formulation and discretization
Let ⊂ R d , d ∈ {2, 3}, be a bounded domain with the Lipschitz boundary ∂ . We consider the following model problem: Find u ∈ H 1 ( ; C) =: X such that where the wave number k ∈ R and the Dirichlet data u 0 ∈ H 1 ( ; C) are given. Except for eigenvalues, the unique solvability of the variational equation (3.1) is argued by the Fredholm's theorem, see e.g. [6,Section 5.3]. The solution of (3.1) fulfills the homogeneous Helmholtz equation: − u−k 2 u = 0 in . Therefore, the interior C ∞ -regularity of u follows from the general theory of linear elliptic PDEs, see [22,Chapter 3]. Moreover, we assume ∂ of class C 1 and u 0 ∈ W 1/ p , p (∂ ; C) for p the conjugate of p and p > d. Then there exists U 0 ∈ W 1, p ( ) such that U 0 = u 0 on ∂ , see [35,Lemma 1.49]. The difference w := u − U 0 has zero trace and solves the inhomogeneous Helmholtz equation: − w − k 2 w = k 2 U 0 − U 0 in . Therefore, the embedding u ∈ W 1, p ( ) ⊂ C( ) holds, see [35,Theorem 3.16]. This provides well-posedness of a pointwise evaluation of u in for the interpolation I X h u (see (3.2) below).
While we consider the Dirichlet problem (3.1), our subsequent considerations can be extended to the cases of Neumann, Robin, and mixed boundary conditions as well as inhomogeneous equations.
We assume that the computational domain h ⊂ R d is polyhedral endowed with a mesh M h and coincides with the physical domain . Its boundary is denoted by ∂ h . In particular, we assume that h consists of uniform elements of size h > 0, quadrilaterals in R 2 and polyhedra in R 3 . Let N h denote the set of all mesh nodes x j ∈ M h ⊂ h , j = 1, . . . , N . We set X h ⊂ H 1 ( h ; C) as continuous piecewise d-linear (bilinear in R 2 , trilinear in R 3 ) finite element functions with respect to the mesh M h over h . The respective FE functions from X 0 h have zero traces at ∂ h . In the following by I X h we denote the classical, piecewise d-linear interpolation operator in C( ; C) =: X fulfilling This operator is stable, see e.g. [5, Lemma 2.1].
The discrete counterpart of (3.1) thus reads: Find u h ∈ X h such that 3) is to be determined according to Algorithm 2.1. For this purpose we construct in the next section a continuous piecewise polynomial basis in H 1 ( h ; C) based on Gegenbauer polynomials.

The continuous piecewise polynomial basis
To construct a continuous piecewise polynomial basis we suggest to use Gegenbauer polynomials G n ∈ P n ([−1, 1]), n = 0, 1, . . . , which are defined by the recursion More of their properties are given in Appendix A.
In the numerical framework, the Gegenbauer polynomials (3.4) are well suited as hierarchical shape functions for hp-FEM, e.g. see [32]. In our theoretical optimization context, they will be used for mother wavelets with vanishing moments in Appendix A.
Based on (3.4) the hierarchical shape functions are given in the local coordinates ξ ∈ [−1, 1] with respect to the binary parameter t ∈ {−1, 1} as In the computational domain, we apply the following affine coordinate transformation given for every fixed mesh point The 2 d adjacent polyhedra in (3.8) forming the patch centered at the point x j are as illustrated in 2d in Fig. 1. Applying (3.7) to (3.6) we get the transformed shape functions as polynomials P n (Q j (t 1 ,...,t d ) ) of degree at most n = max i=1,...,d n i on the polyhedra (3.8) in the reference patch j . We have the following lemma.

Lemma 3.1 The continuous piecewise polynomials of degree n on the uniform polyhedral mesh M h over h ⊂ R d can be spanned with the shape functions (3.9) of degree at most n on the polyhedra.
Proof Firstly we note note that the piecewise d-linear function in (3.9), namely form the usual "hat" function supported on the patch j . Second, it holds because the Gegenbauer polynomials of degree two and more vanish at the boundary (see property (A.3) in Appendix A). Therefore, the shape functions (3.9) of degree at most n form a basis for the continuous piecewise polynomials of degree n supported on the patch j . We associate the center-point x j of j to the nodes N h = {x j } N j=1 of a uniform polyhedral mesh M h over h . Since every continuous piecewise polynomial with respect to the mesh can be partitioned continuously over the set of overlapping patches { j } N j=1 which cover h , this proves the assertion.
For the underlying Helmholtz problem (3.3), from Lemma 3.1 and Theorem A.6 which is proved in Appendix A, the main result of this section will follow.  Let u ∈ H 1 ( ; C) be the exact solution of the variational problem (3.1). As introduced above, let I X h u ∈ X h denote the interpolation over M h satisfying the interpolation condition (3.2). Following (2.7) we look for a test space Y h in Z h satisfying the following equation: (3.10) Due to Theorem 3.2, Y h should be sought within continuous piecewise quadratic polynomials on the mesh M h over h . In fact, one does not gain any benefit by expanding the basis of Y h beyond the quadratic polynomials as their contribution to the system matrix of (3.10) will be zero. Since the dimension of X h is less than Z h , substituting the span of X h and the span of Z h into (3.10) would result in an undetermined problem for the unknown coefficients by the respective basis functions. We reduce the number of quadratic test functions using appropriate simplifications and accounting for particular solutions u for specific data u 0 . As natural simplification, we suggest to rely on symmetric test functions in a finite element basis. Further we provide a dispersion analysis of (3.10) with respect to the particular solutions u specified by plane waves.

The dispersion analysis in 2d
In this section we describe in detail the dispersion analysis in 2d. The 3d results will be presented in Appendix B.
To choose an appropriate finite element basis for the test space Y h we assemble the linear and quadratic shape functions from (3.9) in a globally continuous and symmetric way as follows. In R 2 , at every patch consisting of four adjacent quadrilaterals which are defined for t 1 see Fig. 1, we define the node (linear-linear), edge (linear-quadratic), and bubble (quadratic-quadratic) modes with respect to two spatial dimensions, respectively, by (3.11c) These shape functions are illustrated in Fig. 2. In comparison, the trial space X h consists of bilinear shape functions corresponding to the node mode (3.11a) only. Thus, the test space Y h is enriched compared to X h . Next we apply the finite element ansatz function as Substituting (3.12) and (3.13) into (3.10) we obtain the system matrix, which we denote by A( Its entries are given by for i, j = 1, . . . , N . Henceforth the notation κ := kh 2 is used. These coefficients can be calculated explicitly and we find according to the nine-point interior stencil To determine these unknown weights α N , α E , and α B we employ particular solutions of the Helmholtz equation. In fact, utilizing interior regularity the solution u of (3.1) satisfies the following equation point-wise (3.16) The particular solutions of (3.16) are plane waves of the complex form with an arbitrary incident angle θ ∈ (−π, π]. The piecewise d-linear interpolation I X h u ∈ X h of (3.17) on every polyhedron in the patch j , j = 1, . . . , N , reads where the notation C := cos θ and S := sin θ is used. Inserting Our consideration results in the following. We note that Proposition 3.3 justifies the optimality condition (2.7) given in Step 2 of Algorithm 2.1 for the underlying Helmholtz problem.
In realistic situations, the incident angle θ is unknown a-priori, and the dispersion equation (3.19) cannot be solved for arbitrary angle. Therefore, in the following we rely on the asymptotic model of (3.19) when κ → 0. In view of the considerations of Sect. 2 this will give us an approximate optimality condition (2.8) instead of (2.7) (respectively, instead of (3.10) for the underlying Helmholtz problem).
We look for the weights in (3.13) given in asymptotic form with respect to κ 2 as with five unknown coefficients α N 1 ∈ R, and α E i , α B i ∈ R for i = 0, 1. The chosen number of five unknowns will be argued below (3.21). We substitute (3.15) and (3.20) into (3.19) and apply asymptotic relations: This substitution reduces the dispersion equation to the following asymptotic equality where the residual Er( · , θ) = o(κ 7 ) for θ ∈ (−π, π]. With the following five coefficients (3.24) The corresponding biquadratic finite element basis functions Fig. 3 in dependence of κ. These basis functions correspond to the center-point in the patch shown in Fig. 2. They are a specific combination with the coefficients (α N (κ 2 ), α E (κ 2 ), α B (κ 2 )) of the three shape modes therein. When varying κ we observe a difference between the basis functions shown in the plots (a) and (b) for κ = 2 7 (11 − √ 46) ≈ 1.0977 (see the remark after Theorem 3.4) and κ = 0.25, respectively, while the basis functions depicted in the plots (b) and (c) for κ = 0.25 and κ = 10 −10 are visually indistinguishable. This fact explains that κ need not be chosen very small for computations, in spite of the asymptotic arguments used for the analysis κ → 0. We shall discuss in the next session the choice of κ for the computational realization.
We finish this section with few remarks. It is argued in [3] that, generally, no further reduction of degree in the residual error o(κ 7 ) can be attained in the context of the dispersion equation (3.19). From (3.23) we can determine also the discrete wave number k κ such that |k κ − k| = o(κ 7 ).
In the following section we show well-posedness of (3.3) for the chosen trial and test spaces. Subsequently, we estimate the respective error of discretization.

Well posedness and a-priori error analysis
Now we consider the discrete variational problem problem (3.3) with the test space Y h = Y κ h spanned by the biquadratic finite element functions weighted according to (3.20) and (3.22). The test space is enriched in comparison with the bilinear finite elements in the trial space X h .
To guarantee well posedness of (3.3) in the test space Y κ h , positive definiteness of the system matrix A(κ 2 ) ∈ Sym(N 2 ) with coefficients A 0 (κ 2 ), A 1 (κ 2 ), A 2 (κ 2 ) from (3.24) is needed. We note that the coefficients enter the system matrix A(κ 2 ) in the following way. All nonzero elements in each row as well as in each column consist of nine elements which are 4 A 0 (κ 2 ) (once), 2 A 1 (κ 2 ) (four times), and A 2 (κ 2 ) (four times). We start our consideration with the Laplace operator which corresponds to κ = 0 in (3.24) and then extend it by continuity to κ > 0. We note that the relations A 0 (0) > 0 and A 0 (0) + 2 A 1 (0) + A 2 (0) = 0 hold which imply the usual consistency conditions for the Laplace operator. In fact, we can derive the following properties of the matrix A(0) ∈ Sym(N 2 ): is irreducible: no permutation matrix P exists such that P A(0)P can be reduced.
From (i), (iii), and (iv) it follows that A(0) is irreducibly diagonal dominant, hence nonsingular, see [36]. Since the determinant of A(κ 2 ) is a continuous function of its entries, the following existence theorem holds.

Theorem 3.4 There exists κ 2 0 > 0 (which may be small) such that the discrete variational problem (3.3) stated in the trial space X h and in the test space Y
We can evaluate an upper bound for the constant κ 2 0 in Theorem 3.4 from the necessary condition (i) by requiring: κ 0 < 2 7 (11 − √ 46) ≈ 1.0977. For larger values of κ we obtain negative diagonal entries. In diverse tests the choice κ ≤ 0.25 was successful (see Fig. 4).
Well-posedness can be assured only if k 2 is bounded away from the finitedimensional eigenvalues of the Laplace operator. Otherwise, if k 2 approaches the finite-dimensional eigenvalues it is said to enter a zone of degeneracy, see the related investigation in [9,24,27]. In 3d, in Appendix B for a specific choice of α B 1 we will get A 0 (κ 2 ) > 0 for all κ 2 , which is the improvement compared to the 2d case. Now with the help of Theorem 3.4 we investigate the error of (3.3). We consider the plane wave solution u in (3.17) and its linear interpolate I X h u given in (3.18). The optimality condition (3.10) is not satisfied exactly with the test space Y h = Y κ h . But with (3.23) we can argue that a residual term φ κ h ∈ X h exists such that

25b)
holds (compare with (2.8)). In comparison with (3.25a), the discrete counterpart u h ∈ X h solves the homogeneous equation (3.26) For κ 2 < κ 2 0 , from Theorem 3.4 the inf-sup condition (2.9) follows which in our case has the form of the LBB condition with the LBB constant γ : Therefore, analogously to Proposition 2.2 in Sect. 2, from 3.24 and (3.27) we infer the following result on the asymptotic error.
Theorem 3.5 Let κ 2 < κ 2 0 . For an arbitrary plane wave solution u, due to the asymptotic estimate of the dispersion (3.23), the error between the exact linear interpolation and the FE solution of (3.26) is given by The assertion of Theorem 3.5 holds also in the 3d case, which will be described in Appendix B.
We note that the Dirichlet problem (3.26) can be generalized to an arbitrary mixed setting in the following way. We consider the Dirichlet, Neumann, and Robin boundary conditions stated at pairwise disjoint boundaries D h ∪ N h ∪ R h = ∂ h of the computational domain h endowed with the mesh M h . In the discretized form it leads to the following variational problem: Find u h ∈ X h such that u h = I X h u 0 on D h and (3.29) where the continuous data I X h u 0 and β h are given in h and R h , respectively. Using bilinear finite elements for u h ∈ X h , and the biquadratic finite elements given by (3.11), (3.13) with the weights from (3.20), (3.22) for v h ∈ Y κ h , the stiffness and mass matrices can be calculated according to the terms in (3.29). The high order of interpolation, which was validated for the Dirichlet problem (3.26), may not hold for the variational boundary conditions appearing in (3.29). Possible ways of improvement are discussed in Sect. 4.

The a-posteriori numerical analysis
We compare our approximation by the Petrov-Galerkin enrichment (PGE) with the standard Galerkin least squares (GLS) and other GFEM methods: the quasi-stabilized (c) pollution error c 2 (k,κ) 1/7 Fig. 4 The errors in the selected norms for various FE methods (QSFEM) method introduced in [3], and the variational multiscale (VMS) version from [30], which are well known in the literature. With u h we associate discrete solutions corresponding to these FE methods which are compared with respect to the exact solution u of the Helmholtz problem (3.1) given in the form of a plane wave. We present here the numerical result of tests computed in the unit cube h = = [0, 1] d ⊂ R d for d = 2. Our observations also hold for our numerical tests in 3d.
The methods mentioned above are linear in the H 1 -seminorm, i.e., where c 1 (k, · ) = O(κ) for each k ∈ R recalling that κ = kh 2 . The errors c 1 (8, κ) are depicted in Fig. 4a for the different FE methods as κ varies according to κ = 4h with h ∈ (2 −9 , . . . , 2 −4 ). We observe that the curves for GFEM methods are very close to each other and advantage over GLS. For comparison, we provide here also the quadratic approximation when both trial and test spaces are spanned by continuous piecewise biquadratic polynomials on the uniform quadrilateral mesh M h over h . In this case, c 1 (k, · ) = o(κ). The computational cost of the proposed PGE method is the same as the linear FE methods, while the cost for quadratic FEM is increased.
The difference between the various tested linear approximation methods become apparent when we examine the error with respect to the discrete 2 -norm given over the mesh nodes {x j } N j=1 = N h . These errors are depicted in Fig. 4, firstly, in plot (b) in dependence of κ ∈ (2 −7 , . . . , 2 −2 ) for fixed k = 8 and, second, in plot (c) with respect to varying k ∈ (2 3 , . . . , 2 6 ) for fixed κ = 0.25. The former case in plot (b) describes the asymptotic behavior of the tested methods as κ = kh 2 → 0. Plot (c) represents the pollution effect when increasing the wave number k even fixed parameter κ.
Before detailed discussion of the curves depicted in Fig. 4b, c with respect to the finite 2 -norm, it will be helpful to give another interpretation of these data. To explain  Table 1 The numerical performance of the FE methods as κ → 0 for fixed k it we present in Fig. 5 the error with respect to the linear interpolate solution I X h u ∈ X h in the H 1 -seminorm given by This quantifier expresses the accuracy of the tested FE methods with respect to exact interpolate of the solution. In Fig. 5 we depict the respective curves of c 3 (k, κ), with respect to varying κ ∈ (2 −7 , . . . , 2 −2 ) in the plot (a) for fixed k = 8, as well as varying k ∈ (2 3 , . . . , 2 6 ) in the plot (b) for fixed κ = 0.25. For convenience we present in the following table some observations which can be drawn from Figs. 4 and 5.
Below we comment on the behavior of the FE methods presented in Table 1.
• The standard GLS method has, evidently, the poorest performance.
• The high order (here, quadratic) finite elements have a good approximation property of the exact solution while lagging behind in interpolation in the mesh nodes. • The computational complexity of the PGE as compared to the other linear techniques is comparable since these methods only differ with respect to weight factors of the basis functions. • The VMS, QSFEM, and PGE methods have the smallest pollution error when increasing k even fixed κ = kh 2 . • The QSFEM method has the best performance for large κ ∼ 1. For small κ < 2 −3 , however, its error stops decaying and starts to grow in the tests. The reason is the following one. The corresponding stencil is represented by formula (5.7) in [3] with divided differences of the order O(κ) O(κ) , which results in numerical uncertainty of the division 0 0 as κ → 0. Similarly, we observe in Figs. 4 and 5 that the asymptotic decay of VMS fails for κ < 2 −5 . For smaller κ << 1 also PGE may become numerically unstable.
From the computation tests for moderately small κ ≤ 0.25 we conclude that our PGE method has the interpolation order o(κ 7 ), which justifies numerically the assertion of Theorem 3.4. It approximates the linearly interpolated solution in the most stable way among the tested FE methods.

Concluding remarks
Here we discuss possible further developments of our approach.
In our application to inverse scattering problem, see [21], the Dirichlet data u 0 in (3.1) represents boundary measurement for the reconstruction of u by means of a solution to the Helmholtz equation.
We comment on extension of the formulas in Sect. 3 to Neumann and Robin boundary conditions. The weights (3.20) and (3.22) are found from the dispersion equation (3.19) at interior mesh points. In polyhedra adjacent to mesh points on the boundary, the weights need to be determined according to the dispersion equation for the boundary stencil of the mesh points in the respective patch. The corresponding dispersion equation depends on the boundary condition.
The inhomogeneous Helmholtz equation with a right-hand side f can be transformed to a homogeneous one for u − u 0 with suitable u 0 , even in the case when f resides in the dual space H −1 ( ; C).
The algorithm given in Sect. 2 has the potential of application to various PDEs expressed in variational form. We note that the result of Sect. 3 includes the Laplace equation by choosing k = 0. For vector valued problems, the algorithm can be extended to Raviart-Thomas like finite elements. Higher order finite elements can be employed within the Hermite interpolation as well.
We start with useful properties of the Gegenbauer polynomials. First we remind that G n is the particular solution of the Gegenbauer equation The above equality is satisfied with 1 −1 G n (ξ ) dξ = 0 for n ≥ 3. Analogously, for m = 1 we deduce the equality which is satisfied with either n = 3 or With the help of Lemma A.1 we can recover the orthogonality property of the Gegenbauer polynomials as given below (see [1, p.774]). The property of vanishing moments in (A.4) is used to prove orthogonality of the Gegenbauer polynomials to arbitrary polynomials as follows. 1]) is orthogonal to arbitrary polynomials P ∈ P m ([−1, 1]) with respect to inner products induced by the L 2 -norm and the H 1 -seminorm, respectively, i.e.

Lemma A.2 System {G
Moreover, for every binary t ∈ {−1, 1} a nontrivial linear combination exists which is orthogonal to P ∈ P m ([−1, 1]) with respect to the bilinear form can be spanned by continuous (over h ) and piecewise (on M h ) polynomials of degree n ≥ m + 2. The remaining polynomials of degree n ≤ m + 1 are dense in Z h .
Proof For the ordered set N h of nodes x 1 < x 2 < · · · < x N of the mesh M h , we cover h with the overlapping patches 1 where κ := kh 2 and the polynomial P(ξ ) is of the form We span v t j (ξ ), t ∈ {−1, 1}, with the following shape functions: Due to Lemma A.3 this equality cannot be satisfied with nonzero set of coefficients {V t 1 , . . . , V t m+1 }. Otherwise, the shape functions in (A.11) would be linearly dependent, which would lead to a contradiction.
We conclude that test functions v t j of the form v t j (ξ ) = ∞ n=m+2 V t n B t n (ξ ) are needed to fulfill (A.10). This proves that X ⊥ h is spanned by piecewise polynomials of the local degree n ≥ m + 2. Consequently, the piecewise polynomials of the local degree n ≤ m + 1 are dense in Z h , and the assertion of the theorem follows.
We note that the mesh in 1d does not need be uniform in Theorem A.4. In the following we extend Theorem A.4 to the multidimensional case R d , d ∈ N. Further we define P n as the space of d-dimensional polynomials of degree at most n with respect to every coordinate of the d-dimensional variable ξ = (ξ 1 , . . . , ξ d ), namely, for all n ≥ m + 3. Moreover, using the notion (for i = 1, . . . , d) which satisfy the equality if there exists an index i ∈ {1, . . . , d} such that n i ≥ m + 3 and m i ≤ m. Hence (A.12) holds. We note that from (A.12) it follows that (A.15) holds with w (t 1 ,...,t d ) = G n for all n ≥ m + 3 and arbitrary κ. Now we take w (t 1 ,...,t d ) ∈ P m+2 ([−1, 1] d ) in the form of (A.14) and a polynomial P ∈ P m ([−1, 1] d ) with (m +1) d degrees of freedom. Inserting the linear combination (A.14) into the left hand side of (A.15), for every fixed (t 1 , . . . , t d ) ∈ {1, 1} d we get (m + 1) d homogeneous equations for (m + 2) d unknown coefficients c (t 1 ,...,t d ) n 1 ...n d ∈ R by (n 1 , . . . , n d ) ∈ {1, . . . , m + 2} d . The respective matrix has full rank because the shape functions in (A.11) are linearly independent. Therefore, we find exactly (m + 2) d − (m + 1) d nontrivial linear combinations satisfying (A.15).
We note that the residual set of polynomial functions w (t 1 ,...,t d ) ∈ H 1 ((−1, 1)  Proof We start by defining the transformation of functions depending on x ∈ h to the local coordinates ξ ∈ [−1, 1] d . Let N h = {x j } N j=1 be the set of nodal points of the polyhedral mesh. We associate every point x j , j = 1, . . . , N , to the patch j of all neighbor elements containing x j . These neighboring elements consist of 2 d polyhedra Q j (t 1 ,...,t d ) indexed by the binary numbers (t 1 , . . . , t d ) ∈ {−1, 1} d . The multi-index t = (t 1 , . . . , t d ) is chosen in such a way that it corresponds to the "hat" function centered at the node x j . Namely, the reference patch j = {x i ∈ R : associates the stencil of 3 d -points and consists of the polyhedra Using a partition of unity, every function v ∈ H 1 ( h ; C) can be partitioned by v = N j=1 v j continuously over the family of overlapping patches { j } N j=1 covering the computational domain h in such a way that each part v j (x) ∈ H 1 0 ( j ; C) is supported on the reference patch j . We transform each polyhedron Q with the following stencil coefficients Substituting for the trial function u h = I X h u the plane wave solution u(x) = e ık(x 1 C 1 S 2 +x 2 S 1 S 2 +x 3 C 2 ) , C i := cos θ i , S i := sin θ i , i = 1, 2, with arbitrary incident angle θ = (θ 1 , θ 2 ) ∈ (−π, π] × [0, π], into the discrete variational problem (3.11) leads to the dispersion equation with factors K 1 := cos(2κC 1 S 2 ) + cos(2κ S 1 S 2 ) + cos(2κC 2 ), K 2 := cos(2κC 1 S 2 ) cos(2κ S 1 S 2 ) + cos(2κC 1 S 2 ) cos(2κC 2 ) + cos(2κ S 1 S 2 ) cos(2κC 2 ), K 3 := cos(2κC 1 S 2 ) cos(2κ S 1 S 2 ) cos(2κC 2 ).