Abstract
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.
Similar content being viewed by others
1 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.
2 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\in X\) find \(u-u_0\in X^0\) such that
where the map \(b: X\times X\mapsto {\mathbb {C}}\) is continuous and sesquilinear over \({\mathbb {C}}\) (respectively, bilinear over \({\mathbb {R}}\)).
Assuming the existence of a variational solution to (2.1) we look for its finite dimensional approximation in a trial space \(X_h\subset X\). Let \(X^0_h\) be a closed subspace of \(X_h\) and \(u_0^h\in X_h\). For finite elements, the index \(h\in {\mathbb {R}}_+\) refers to the mesh size. Following the Petrov–Galerkin approach we set a discrete counterpart of (2.1) in the form: Find \(u^h -u_0^h\in X^0_h\) such that
over a test space \(Y_h\subset X\), where the continuous and sesquilinear form \(b_h: X\times X\mapsto {\mathbb {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 \({\mathcal {B}}_X :=\{\phi _i\in X\vert i=1,\dots ,\infty \}\) form a FE basis of X. In this basis, we denote by \(X^\perp _h\subset X\) the orthogonal complement to \(X_h\) with respect to (2.2), namely
Noting that in general \(X \not =X_h\bigoplus X^\perp _h\) we define the space
We exclude the trivial case \(Z_h =\emptyset \). Indeed, in the subsequent considerations we have the direct sum \(X =Z_h\bigoplus X^\perp _h\) in a piecewise polynomial basis \({\mathcal {B}}_X\), and \(X_h\subset Z_h\). Since \(X^\perp _h\) is contained in the kernel of the linear operator \(b_h(u^h,\,\cdot \,): X\mapsto {\mathbb {C}}\) in (2.2), this suggests to look for a test space \(Y_h\) in \(Z_h\).
The following discussion will rely on an interpolation operator \(I_{X_h}: {\widetilde{X}}\mapsto X_h\) defined on a suitable subspace \({\widetilde{X}}\subseteq 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 \({\widetilde{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 \({\widetilde{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 \(\Vert u-I_{X_h}u\Vert _X\) is the interpolation error, and \(\Vert I_{X_h}u-u^h\Vert _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).
Algorithm 2.1
Fix the form \(b_h\) in (2.2) and set the trial space \(X_h\).
Step 1. Find the orthogonal complement \(X^\perp _h\) in (2.3).
Step 2. Determine the test space \(Y_h\) by weighting the basis functions in \(Z_h\) to satisfy the necessary optimality condition
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\in {\widetilde{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:
where \(\Phi _h: X\mapsto {\mathbb {C}}\) is linear, and \(\mathrm{Er}:{\mathbb {R}}_+\mapsto {\mathbb {R}}_+\) expresses an error function.
Proposition 2.2
Let the test space \(Y_h\) be such that (2.8) holds, and let the form \(b_h\) in (2.2) satisfy the following inf-sup condition:
Then the error of the FE solution of (2.2) can be estimated by
Proof
Due to (2.9), estimate (2.10) follows by subtracting (2.2) from (2.8). \(\square \)
We note that from (2.9) the stability estimate \(\Vert u^h -u^0\Vert _X \le {\textstyle \frac{1}{\gamma }} \Vert b_h(u^0,\,\cdot \,)\Vert _{X^{-1}}\) follows, where \(X^{-1}\) denotes the dual space of X.
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 \(\mathrm{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 \(\kappa :=\frac{kh}{2}\rightarrow 0\), where k stands for the wave number, we obtain the error function \(\mathrm{Er}(h)=\mathrm{o}(\kappa ^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 \(\mathrm{o}(\kappa ^7)\).
3 The Helmholtz equation
In this section we first formulate the Helmholtz problem under consideration and then show how to use Algorithm 2.1.
3.1 The Helmholtz problem formulation and discretization
Let \(\Omega \subset {\mathbb {R}}^d\), \(d\in \{2,3\}\), be a bounded domain with the Lipschitz boundary \(\partial \Omega \). We consider the following model problem: Find \(u\in H^1(\Omega ;{\mathbb {C}}) =:X\) such that
where the wave number \(k\in {\mathbb {R}}\) and the Dirichlet data \(u_0\in H^1(\Omega ;{\mathbb {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: \(-\Delta u -k^2 u =0\) in \(\Omega \). Therefore, the interior \(C^\infty \)-regularity of u follows from the general theory of linear elliptic PDEs, see [22, Chapter 3]. Moreover, we assume \(\partial \Omega \) of class \(C^1\) and \(u_0\in W^{1/p',p}(\partial \Omega ;{\mathbb {C}})\) for \(p'\) the conjugate of p and \(p>d\). Then there exists \(U_0\in W^{1,p}(\Omega )\) such that \(U_0 =u_0\) on \(\partial \Omega \), see [35, Lemma 1.49]. The difference \(w:=u-U_0\) has zero trace and solves the inhomogeneous Helmholtz equation: \(-\Delta w -k^2 w =k^2 U_0 -\Delta U_0\) in \(\Omega \). Therefore, the embedding \(u\in W^{1,p}(\Omega )\subset C({\overline{\Omega }})\) holds, see [35, Theorem 3.16]. This provides well-posedness of a pointwise evaluation of u in \({\overline{\Omega }}\) 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 \(\Omega _h\subset {\mathbb {R}}^d\) is polyhedral endowed with a mesh \(M_h\) and coincides with the physical domain \(\Omega \). Its boundary is denoted by \(\partial \Omega _h\). In particular, we assume that \(\Omega _h\) consists of uniform elements of size \(h>0\), quadrilaterals in \({\mathbb {R}}^2\) and polyhedra in \({\mathbb {R}}^3\). Let \(N_h\) denote the set of all mesh nodes \(x^j\in M_h\subset {\overline{\Omega }}_h\), \(j=1,\dots ,N\). We set \(X_h\subset H^1(\Omega _h;{\mathbb {C}})\) as continuous piecewise d-linear (bilinear in \({\mathbb {R}}^2\), trilinear in \({\mathbb {R}}^3\)) finite element functions with respect to the mesh \(M_h\) over \(\Omega _h\). The respective FE functions from \(X^0_h\) have zero traces at \(\partial \Omega _h\).
In the following by \(I_{X_h}\) we denote the classical, piecewise d-linear interpolation operator in \(C({\overline{\Omega }};{\mathbb {C}}) =:{\widetilde{X}}\) fulfilling
This operator is stable, see e.g. [5, Lemma 2.1].
The discrete counterpart of (3.1) thus reads: Find \(u^h\in X_h\) such that
The unknown test space \(Y_h\subset H^1(\Omega _h;{\mathbb {C}})\) in (3.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(\Omega _h;{\mathbb {C}})\) based on Gegenbauer polynomials.
3.2 The continuous piecewise polynomial basis
To construct a continuous piecewise polynomial basis we suggest to use Gegenbauer polynomials \(G_n\in {\mathbb {P}}_n([-1,1])\), \(n=0,1,\dots \), 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 \(\xi \in [-1,1]\) with respect to the binary parameter \(t\in \{-1,1\}\) as
For vector coordinates \(\xi =(\xi _1,\dots ,\xi _d)\in [-1,1]^d\) and binary multi-indices \((t_1,\dots ,t_d)\in \{-1,1\}^d\), we define the d-dimensional shape functions as
which are d-dimensional polynomials \({\mathbb {P}}_n([-1,1]^d)\) of degree at most \(n={\displaystyle \max _{i=1,\dots ,d}} n_i\).
In the computational domain, we apply the following affine coordinate transformation given for every fixed mesh point \(x^j=(x^j_1,\dots ,x^j_d)\in {\mathbb {R}}^d\) and \((t_1,\dots ,t_d)\in \{-1,1\}^d\) by
which transforms the parent domain \([-1,1]^d\) to the polyhedra
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 \({\mathbb {P}}_n(Q^j_{(t_1,\dots ,t_d)})\) of degree at most \(n={\displaystyle \max _{i=1,\dots ,d}} n_i\) on the polyhedra (3.8) in the reference patch \(\Pi ^j\). We have the following lemma.
Lemma 3.1
The continuous piecewise polynomials of degree n on the uniform polyhedral mesh \(M_h\) over \(\Omega _h\subset {\mathbb {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 \(\Pi ^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 \(\Pi ^j\). We associate the center-point \(x^j\) of \(\Pi ^j\) to the nodes \(N_h=\{x^j\}_{j=1}^N\) of a uniform polyhedral mesh \(M_h\) over \(\Omega _h\). Since every continuous piecewise polynomial with respect to the mesh can be partitioned continuously over the set of overlapping patches \(\{\Pi ^j\}_{j=1}^N\) which cover \({\overline{\Omega }}_h\), this proves the assertion. \(\square \)
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.
Theorem 3.2
Let \(X_h\) be the space of continuous piecewise d-linear polynomials over \({\mathbb {C}}\) on a uniform polyhedral mesh \(M_h\) over \(\Omega _h\subset {\mathbb {R}}^d\). The orthogonal complement with respect to (3.3) (see (A.16))
can be spanned by continuous piecewise polynomials of degree \(n\ge 3\) on the mesh \(M_h\) over \(\Omega _h\). Moreover, continuous piecewise polynomials of degree at most \(n=2\) form a basis for the finite-dimensional space \(Z_h\) for fixed h.
By Theorem 3.2 we next construct the test space \(Y_h\) in \(Z_h\) according to Step 2 in Algorithm 2.1 as follows.
Let \(u\in H^1(\Omega ;{\mathbb {C}})\) be the exact solution of the variational problem (3.1). As introduced above, let \(I_{X_h}u\in 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:
Due to Theorem 3.2, \(Y_h\) should be sought within continuous piecewise quadratic polynomials on the mesh \(M_h\) over \(\Omega _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.
3.3 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 \({\mathbb {R}}^2\), at every patch \(\Pi ^j =Q^j_{(1,-1)}\cup Q^j_{(-1,-1)}\cup Q^j_{(1,1)}\cup Q^j_{(-1,1)}\) consisting of four adjacent quadrilaterals which are defined for \(t_1,t_2\in \{-1,1\}\) by
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
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
with unknown coefficients \(\{U^j\}_{j=1}^N\in {\mathbb {C}}^N\), and
with unknown weights \(\alpha ^N, \alpha ^E, \alpha ^B \in {\mathbb {R}}\).
Substituting (3.12) and (3.13) into (3.10) we obtain the system matrix, which we denote by \(A(\kappa ^2) =\{a_{ij}(\kappa ^2)\}_{i,j=1}^N \in {\mathbb {R}}^{N\times N}\). Its entries are given by
for \(i,j=1,\dots ,N\). Henceforth the notation \(\kappa :=\frac{kh}{2}\) is used. These coefficients can be calculated explicitly and we find
according to the nine-point interior stencil
with stencil coefficients
and unknown weights \(\alpha ^N, \alpha ^E, \alpha ^B\).
To determine these unknown weights \(\alpha ^N\), \(\alpha ^E\), and \(\alpha ^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
The particular solutions of (3.16) are plane waves of the complex form
with an arbitrary incident angle \(\theta \in (-\pi ,\pi ]\). The piecewise d-linear interpolation \(I_{X_h}u\in X_h\) of (3.17) on every polyhedron in the patch \(\Pi ^j\), \(j=1,\dots ,N\), reads
Substituting expressions (3.18) and (3.13) into (3.10), for the interior stencil we derive the usual dispersion equation with respect to the incident angle \(\theta \):
where the notation \(C:=\cos \theta \) and \(S:=\sin \theta \) is used. Inserting (3.15) into (3.19), one linear equation \({\mathcal {D}}(\alpha ^N, \alpha ^E, \alpha ^B; \kappa ^2, \theta )=0\) for three unknowns \(\alpha ^N\), \(\alpha ^E\), and \(\alpha ^B\) in dependence of \(\kappa ^2\) and \(\theta \) is obtained. For each \(\kappa =\frac{kh}{2}\) and for each \(\theta \), equality (3.19) implies a linear condition which can be solved for the variables
in dependence of the parameters \((\kappa ^2,\theta )\).
Our consideration results in the following.
Proposition 3.3
If the incident direction \((\cos \theta ,\sin \theta )\) in \({\mathbb {R}}^2\) is fixed a-priori, then there exist nontrivial weights \(\alpha ^N(\kappa ^2,\theta ), \alpha ^E(\kappa ^2,\theta ), \alpha ^B(\kappa ^2,\theta )\in {\mathbb {R}}\) which solve the dispersion equation (3.19). These weights determine the finite element basis (3.13) for the test space \(Y_h\). The exact interpolation \(I_{X_h}u\) from (3.18) solves problem (3.10) stated over \(Y_h\). If the solution to (3.10) is unique, then it coincides with \(I_{X_h}u\).
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 \(\theta \) 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 \(\kappa \rightarrow 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 \(\kappa ^2\) as
with five unknown coefficients \(\alpha ^N_1\in {\mathbb {R}}\), and \(\alpha ^E_i, \alpha ^B_i\in {\mathbb {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 \(\mathrm{Er}(\,\cdot \,,\theta ) =\mathrm{o}(\kappa ^7)\) for \(\theta \in (-\pi ,\pi ]\). With the following five coefficients
we eliminate the five linear combinations by the terms of order \(\mathrm{O}(\kappa ^2)\), \(\mathrm{O}(\kappa ^4)\), \(\mathrm{O}(\kappa ^6)\), \(\mathrm{O}((CS)^2\kappa ^4)\), and \(\mathrm{O}((CS)^2\kappa ^6)\) in (3.21), which provides the asymptotic equality
Inserting (3.22) into (3.15) we get finally the following stencil coefficients
The corresponding biquadratic finite element basis functions
for \(x\in Q^j_{(t_1,t_2)}\), \(t_1,t_2\in \{-1,1\}\), are depicted in Fig. 3 in dependence of \(\kappa \).
These basis functions correspond to the center-point in the patch shown in Fig. 2. They are a specific combination with the coefficients \((\alpha ^N(\kappa ^2), \alpha ^E(\kappa ^2), \alpha ^B(\kappa ^2))\) of the three shape modes therein. When varying \(\kappa \) we observe a difference between the basis functions shown in the plots (a) and (b) for \(\kappa =\sqrt{\frac{2}{7}(11-\sqrt{46})}\approx 1.0977\) (see the remark after Theorem 3.4) and \(\kappa =0.25\), respectively, while the basis functions depicted in the plots (b) and (c) for \(\kappa =0.25\) and \(\kappa =10^{-10}\) are visually indistinguishable. This fact explains that \(\kappa \) need not be chosen very small for computations, in spite of the asymptotic arguments used for the analysis \(\kappa \rightarrow 0\). We shall discuss in the next session the choice of \(\kappa \) 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 \(\mathrm{o}(\kappa ^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_\kappa \) such that \(|k_\kappa -k| =\mathrm{o}(\kappa ^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.
3.4 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^\kappa \) 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^\kappa \), positive definiteness of the system matrix \(A(\kappa ^2)\in \mathrm{Sym}(N^2)\) with coefficients \(A_0(\kappa ^2), A_1(\kappa ^2), A_2(\kappa ^2)\) from (3.24) is needed. We note that the coefficients enter the system matrix \(A(\kappa ^2)\) in the following way. All nonzero elements in each row as well as in each column consist of nine elements which are \(4A_0(\kappa ^2)\) (once), \(2A_1(\kappa ^2)\) (four times), and \(A_2(\kappa ^2)\) (four times). We start our consideration with the Laplace operator which corresponds to \(\kappa =0\) in (3.24) and then extend it by continuity to \(\kappa >0\). We note that the relations \(A_0(0)>0\) and \(A_0(0)+2A_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)\in \mathrm{Sym}(N^2)\):
-
(i)
A(0) has positive diagonal entries.
-
(ii)
A(0) is an L-matrix: \(a_{ij}(0)\le 0\) for \(j\not =i\), \(a_{ii}(0)>0\).
-
(iii)
A(0) is diagonally dominant: \(\sum _{j\not =i} |a_{ij}(0)|\le |a_{ii}(0)|\).
-
(iv)
A(0) is irreducible: no permutation matrix P exists such that \(PA(0)P^\top \) 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(\kappa ^2)\) is a continuous function of its entries, the following existence theorem holds.
Theorem 3.4
There exists \(\kappa _0^2>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_h^\kappa \) is well posed for all \(\kappa ^2\le \kappa _0^2\), \(\kappa :=\frac{kh}{2}\).
We can evaluate an upper bound for the constant \(\kappa _0^2\) in Theorem 3.4 from the necessary condition (i) by requiring: \(\kappa _0<\sqrt{\frac{2}{7}(11-\sqrt{46})}\approx 1.0977\). For larger values of \(\kappa \) we obtain negative diagonal entries. In diverse tests the choice \(\kappa \le 0.25\) was successful (see Fig. 4).
Well-posedness can be assured only if \(k^2\) is bounded away from the finite-dimensional 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 \(\alpha ^B_1\) we will get \(A_0(\kappa ^2)>0\) for all \(\kappa ^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^\kappa \). But with (3.23) we can argue that a residual term \(\phi _h^\kappa \in X_h\) exists such that
holds (compare with (2.8)). In comparison with (3.25a), the discrete counterpart \(u^h\in X_h\) solves the homogeneous equation
For \(\kappa ^2<\kappa _0^2\), 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 \(\gamma \):
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 \(\kappa ^2<\kappa _0^2\). 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 \(\Gamma _h^D\cup \Gamma _h^N\cup \Gamma _h^R =\partial \Omega _h\) of the computational domain \(\Omega _h\) endowed with the mesh \(M_h\). In the discretized form it leads to the following variational problem: Find \(u^h\in X_h\) such that \(u^h=I_{X_h}u_0\) on \(\Gamma _h^D\) and
where the continuous data \(I_{X_h}u_0\) and \(\beta _h\) are given in \({\overline{\Omega }}_h\) and \(\Gamma _h^R\), respectively. Using bilinear finite elements for \(u^h\in X_h\), and the biquadratic finite elements given by (3.10), (3.13) with the weights from (3.20), (3.22) for \(v^h\in Y_h^\kappa \), 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.
3.5 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 (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 \(\Omega _h =\Omega =[0,1]^d\subset {\mathbb {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,\,\cdot \,) =\mathrm{O}(\kappa )\) for each \(k\in {\mathbb {R}}\) recalling that \(\kappa =\frac{kh}{2}\). The errors \(c_1(8,\kappa )\) are depicted in Fig. 4a for the different FE methods as \(\kappa \) varies according to \(\kappa =4h\) with \(h\in (2^{-9},\dots ,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 \(\Omega _h\). In this case, \(c_1(k,\,\cdot \,)=\mathrm{o}(\kappa )\). 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 \(\ell ^2\)-norm
given over the mesh nodes \(\{x^j\}_{j=1}^N =N_h\). These errors are depicted in Fig. 4, firstly, in plot (b) in dependence of \(\kappa \in (2^{-7},\dots ,2^{-2})\) for fixed \(k=8\) and, second, in plot (c) with respect to varying \(k\in (2^{3},\dots ,2^{6})\) for fixed \(\kappa =0.25\). The former case in plot (b) describes the asymptotic behavior of the tested methods as \(\kappa =\frac{kh}{2}\rightarrow 0\). Plot (c) represents the pollution effect when increasing the wave number k even fixed parameter \(\kappa \).
Before detailed discussion of the curves depicted in Fig. 4b, c with respect to the finite \(\ell ^2\)-norm, it will be helpful to give another interpretation of these data. To explain it we present in Fig. 5 the error with respect to the linear interpolate solution \(I_{X_h}u\in 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,\kappa )\), with respect to varying \(\kappa \in (2^{-7},\dots ,2^{-2})\) in the plot (a) for fixed \(k=8\), as well as varying \(k\in (2^{3},\dots ,2^{6})\) in the plot (b) for fixed \(\kappa =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 \(\kappa =\frac{kh}{2}\).
-
The QSFEM method has the best performance for large \(\kappa \sim 1\). For small \(\kappa <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 \(\frac{\mathrm{O}(\kappa )}{\mathrm{O}(\kappa )}\), which results in numerical uncertainty of the division \(\frac{0}{0}\) as \(\kappa \rightarrow 0\). Similarly, we observe in Figs. 4 and 5 that the asymptotic decay of VMS fails for \(\kappa <2^{-5}\). For smaller \(\kappa<\!<1\) also PGE may become numerically unstable.
From the computation tests for moderately small \(\kappa \le 0.25\) we conclude that our PGE method has the interpolation order \(\mathrm{o}(\kappa ^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.
4 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}(\Omega ;{\mathbb {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.
References
Abramowitz, M., Stegun, I.A.: Handbook of Mathematical Functions. National Bureau of Standards, Washington, DC (1972)
Ammari, H., Kang, H., Kim, E.J., Lim, M., Louati, K.: A direct algorithm for ultrasonic imaging of internal corrosion. SIAM J. Numer. Anal. 49, 1177–1193 (2000)
Babuška, I., Ihlenburg, F., Paik, E.T., Sauter, S.A.: A generalized finite element method for solving the Helmholtz equation in two dimensions with minimal pollution. Comput. Methods Appl. Mech. Eng. 128, 325–359 (1995)
Barrenechea, G.R., Franca, L.P., Valentin, F.: A Petrov–Galerkin enriched method: a mass conservative finite element method for the Darcy equation. Comput. Methods Appl. Mech. Eng. 196, 2449–2464 (2007)
Cai, X.C.: The use of pointwise interpolation in domain decomposition methods with nonnested meshes. SIAM J. Sci. Comput. 16, 250–256 (1995)
Cakoni, F., Colton, D.: A Qualitative Approach to Inverse Scattering Theory, vol. 188. AMS. Springer, Berlin (2014)
Cakoni, F., Kovtunenko, V.A.: Topological optimality condition for the identification of the center of an inhomogeneity. Inverse Probl. 34, 035009 (2018)
Chandler-Wilde, S.N., Langdon, S.: A Galerkin boundary element method for high frequency scattering by convey polygons. SIAM J. Numer. Anal. 43, 610–640 (2007)
Demkowicz, L.: Asymptotic convergence in finite and boundary element methods: part 1: theoretical results. Comput. Math. Appl. 27, 69–84 (1994)
Fremiot, G., Horn, W., Laurain, A., Rao, M., Sokolowski, J.: On the analysis of boundary value problems in nonsmooth domains. Dissertationes Mathematicae. Inst. Math. Polish Acad. Sci, Warsaw 462 (2009)
Harari, I., Gosteev, K.: Bubble-based stabilization for the Helmholtz equation. Int. J. Numer. Methods Eng. 70, 1241–1260 (2007)
Heikkola, E., Kuznetsov, Yu.A., Neittaanmäki, P., Toivanen, J.: Fictitious domain methods for the numerical solution of two-dimensional scattering problems. J. Comput. Phys. 145, 89–109 (1998)
Hughes, T.J.R.: Multiscale phenomena: Green’s functions, the Dirichlet-to-Neumann formulation, subgrid scale models, bubbles and the origins of stabilized methods. Comput. Methods Appl. Mech. Eng. 127, 387–401 (1995)
Hughes, T.J.R., Sangalli, G.: Variational multiscale analysis: the fine scale Green’s function, projection, optimization, localization, and stabilized methods, SIAM. J. Numer. Anal. 45, 539–557 (2007)
Ihlenburg, F.: Finite Element Analysis of Acoustic Scattering. Springer, New York (1998)
Kaltenbacher, B., Neubauer, A., Scherzer, O.: Iterative Regularization Methods for Nonlinear Ill-Posed Problems. de Gruyter, Berlin (2008)
Khludnev, A.M., Kovtunenko, V.A.: Analysis of Cracks in Solids. WIT-Press, Southampton (2000)
Kovtunenko, V.A.: High-order topological expansions for Helmholtz problems in 2d. In: Bergounioux, M., et al. (eds.) Topological Optimization, vol. 17, pp. 64–122. Radon Ser. Comput. Appl Math. de Gruyter, Berlin (2017a)
Kovtunenko, V.A.: Two-parameter topological expansion of Helmholtz problems with inhomogeneity. In: Itou, H., et al. (eds.) Mathematical Analysis of Continuum Mechanics and Industrial Applications, vol. 26, pp. 51–81. Mathematics for Industry. Springer, Singapore (2017b)
Kovtunenko, V.A., Kunisch, K.: Problem of crack perturbation based on level sets and velocities. Z. Angew. Math. Mech. 87, 809–830 (2007)
Kovtunenko, V.A., Kunisch, K.: High precision identification of an object: optimality conditions based concept of imaging. SIAM J. Control Optim. 52, 773–796 (2014)
Ladyzhenskaya, O.A., Uraltseva, N.N.: Linear and Quasilinear Elliptic Equations. Academic Press, New York (1968)
Langer, U., Steinbach, O., Wendland, W.L. (eds.) Fast boundary element methods in industrial applications. Comput. Visual. Sci. 8, 3–4. Springer, Heidelberg (2005)
Li, Z.C.: The Trefftz method for the Helmholtz equation with degeneracy. Appl. Numer. Math. 58, 131–159 (2008)
Loula, A.F.D., Fernandes, D.T.: A quasi optimal Petrov–Galerkin method for Helmholtz problem. Int. J. Numer. Methods Eng. 80, 1595–1622 (2009)
Melenk, J.M., Babuška, I.: The partition of unity finite element method: basic theory and applications. Comput. Methods Appl. Mech. Eng. 139, 289–314 (1996)
Nadukandi, P., Oñate, E., Garcia, J.: A fourth-order compact scheme for the Helmholtz equation: Alpha-interpolation of FEM and FDM stencils. Int. J. Numer. Methods Eng. 86, 18–46 (2011a)
Nadukandi, P., Oñate, E., Garcia, J.: A Petrov-Galerkin formulation for the alpha interpolation of FEM and FDM stencils: applications to the Helmholtz equation. Int. J. Numer. Methods Eng. 89, 1367–1391 (2011b)
Oberai, A.A., Pinsky, P.M.: A multiscale finite element method for the Helmholtz equation. Comput. Methods Appl. Mech. Eng. 154, 281–297 (1998)
Oberai, A.A., Pinsky, P.M.: A residual-based finite element method for the Helmholtz equation. Int. J. Numer. Methods Eng. 49, 399–419 (2000)
Scherzer, O., Engl, H.W., Kunisch, K.: Optimal a posteriori parameter choice for Tikhonov regularization for solving nonlinear ill-posed problems. SIAM J. Numer. Anal. 30, 1796–1838 (1993)
Schröder, A.: Constraints coefficients in hp-FEM. In: Kunisch, K., Of, G., Steinbach, O. (eds.) Numerical Mathematics and Advanced Applications. Proceedings of ENUMATH 2007, pp. 183–190. Springer, Berlin (2008)
Strouboulis, T., Babuška, I., Hidajat, R.: The generalized finite element method for Helmholtz equation: theory, computation, and open problems. Comput. Meth. Appl. Mech. Eng. 195, 4711–4731 (2006)
Thompson, L.L.: A review of finite element methods for time-harmonic acoustics. J. Acoust. Soc. Am. 119, 1315–1330 (2006)
Troianiello, G.M.: Elliptic Differential Equations and Obstacle Problems. Plenum Press, New York (1987)
Varga, R.S.: Matrix Iterative Analysis. Springer, Berlin (2000)
Acknowledgements
Open access funding provided by Austrian Science Fund (FWF).
Author information
Authors and Affiliations
Corresponding author
Additional information
The research results were obtained with the support of the Austrian Science Fund (FWF) in the framework of the research projects P26147-N26: “Object identification problems: numerical analysis” (PION), the SFB F32 “Mathematical Optimization and Applications in Biomedical Sciences”, the ERC advanced Grant 668998: OCLOC, and the Austrian Academy of Sciences (OeAW).
Appendices
Appendix A: Gegenbauer polynomials and mother wavelets with vanishing moments
The Gegenbauer polynomials \(G_n\in {\mathbb {P}}_n([-1,1])\) of degree \(n=0,1,\dots \) are defined by the recursion
They are used for the hierarchical shape functions, see [32], given in the local coordinates \(\xi \in [-1,1]\) with respect to the binary parameter \(t\in \{-1,1\}\) as
We start with useful properties of the Gegenbauer polynomials. First we remind that \(G_n\) is the particular solution of the Gegenbauer equation
Since \(G_2(\pm 1)=G_3(\pm 1)=0\), due to the two term recursion (A.1), it follows that
Based on (A.2) and (A.3) we associate \(G_n\) to local wavelets as follows.
Lemma A.1
For \(n\ge 3\), the Gegenbauer polynomials \(\{G_n\}\) form a system of mother wavelets (up to a normalization factor) with vanishing moments such that
Proof
For fixed \(n\ge 3\), we prove (A.4) by induction over m. For \(m=0\), by (A.2), (A.3), and applying integration by parts we derive that
The above equality is satisfied with \(\int _{-1}^1 G_n(\xi )\,d\xi =0\) for \(n\ge 3\). Analogously, for \(m=1\) we deduce the equality
which is satisfied with either \(n=3\) or \(\int _{-1}^1 \xi G_n(\xi )\,d\xi =0\) for \(n\ge 4\).
If \(\int _{-1}^1 \xi ^{m-2} G_n(\xi )\,d\xi =0\), then we obtain for \(0\le m\le n-3\)
hence, either \(m=n-2\) or \(\int _{-1}^1 \xi ^m G_n(\xi )\,d\xi =0\) for \(m\le n-3\). This proves the assertion. \(\square \)
With the help of Lemma A.1 we can recover the orthogonality property of the Gegenbauer polynomials as given below (see [1, p.774]).
Lemma A.2
System \(\{G_n\}_{n=2}^\infty \) forms an orthogonal basis in \(C^\infty _0([-1,1])\) with respect to the weighted \(L^2\)-inner product
The property of vanishing moments in (A.4) is used to prove orthogonality of the Gegenbauer polynomials to arbitrary polynomials as follows.
Lemma A.3
For fixed \(m\ge 0\), the basis \(\{G_n\}_{n=m+3}^\infty \in H^1_0(-1,1)\cap C^\infty ([-1,1])\) is orthogonal to arbitrary polynomials \(P\in {\mathbb {P}}_m([-1,1])\) with respect to inner products induced by the \(L^2\)-norm and the \(H^1\)-seminorm, respectively, i.e.
Moreover, for every binary \(t\in \{-1,1\}\) a nontrivial linear combination
exists which is orthogonal to \(P\in {\mathbb {P}}_m([-1,1])\) with respect to the bilinear form
Proof
The \(L^2\)-orthogonality in (A.6) follows immediately from (A.4) for \(n\ge m+3\). The orthogonality with respect to the \(H^1\)-seminorm in (A.6) is obtained after integration by parts due to \(P^{\prime \prime }\in {\mathbb {P}}_{m-2}([-1,1])\) and (A.3), (A.4). For \(m+1\) arbitrary coefficients in the monomials forming P, equality (A.8) results in \(m+1\) linear homogeneous equations, which can be solved with \(m+2\) nontrivial coefficients \(\{c^t_i\}_{i=1}^{m+2}\) in (A.7) for every \(t\in \{-1,1\}\). This completes the proof. \(\square \)
Lemma A.3 is used for the proof of the following theorem on decomposition of the \(H^1\)-space in 1d.
Theorem A.4
For a given mesh \(M_h\) of size \(h>0\) in a segment \(\Omega _h\subset {\mathbb {R}}^1\), let \(X_h\) be the set of continuous (over \(\Omega _h\)) and piecewise (on \(M_h\)) polynomials of degree \(m\ge 1\) over \({\mathbb {C}}\). Then for given \(k\in {\mathbb {R}}\), the orthogonal complement
can be spanned by continuous (over \(\Omega _h\)) and piecewise (on \(M_h\)) polynomials of degree \(n\ge m+2\). The remaining polynomials of degree \(n\le m+1\) are dense in \(Z_h\).
Proof
For the ordered set \(N_h\) of nodes \(x^1<x^2<\dots <x^N\) of the mesh \(M_h\), we cover \({\overline{\Omega }}_h\) with the overlapping patches \(\Pi ^1 =Q^1_{1}\cup Q^1_{-1}\), \(\dots \), \(\Pi ^{N} =Q^{N}_{1}\cup Q^{N}_{-1}\), where \(Q^j_{1} :=[x^{j-1},x^{j}]\) for \(j=2,\dots ,N\), and \(Q^j_{-1} :=[x^j,x^{j+1}]\) for \(j=1,\dots ,N-1\). Depending on the binary index \(t\in \{-1,1\}\), every segment \(Q^j_t\), \(j=1,\dots ,N\), can be transformed to the parent domain \([-1,1]\) with the help of the coordinate transformation \(x=x^j+\frac{h}{2}(\xi -t)\), where \(h=x^j-x^{j-1}\) for \(t=1\) and \(h=x^{j+1}-x^j\) for \(t=-1\), yielding the Jacobian \(\frac{h}{2}\).
Given a partition of unity subordinated to the covering \(\cup _{j=1}^N \Pi ^j\cap {\overline{\Omega }}_h\), every function \(v\in H^1(\Omega _h;{\mathbb {C}})\) can be partitioned in a continuous way by \(v =\sum _{j=1}^N v^j\). Every such function \(v^j(x)\in H^1_0(\Pi ^j;{\mathbb {C}})\) is supported on the patch \(\Pi ^j\), \(j=1,\dots ,N\). It can be transformed to the function pair \(v^j(x^j+\frac{h}{2}(\xi -t)) =:v^t_j(\xi )\in H^1([-1,1];{\mathbb {C}})\) which is defined in the parent domain for \(t\in \{-1,1\}\) such that \(v^1_j(-1)=v^{-1}_j(1)=0\) and \(v^1_j(1)=v^{-1}_j(-1)\). After transformation, the equation in (A.9) implies
where \(\kappa :=\frac{kh}{2}\) and the polynomial \(P(\xi )\) is of the form \(w^h(x^j+\frac{h}{2}(\xi -t)) \in {\mathbb {P}}_{m}([-1,1];{\mathbb {C}})\).
We span \(v^t_j(\xi )\), \(t\in \{-1,1\}\), with the following shape functions:
where \(B^t_n\in {\mathbb {P}}_{n}([-1,1])\). The linear combination \(w^t_{m+2}\) as defined in (A.7) is chosen such that (A.8) is fulfilled. The choice of this basis implies that (A.10) is satisfied by (A.6) for the test functions \(v^t_j=B^t_n\) for \(n\ge m+3\), and by (A.8) for \(v^t_j=B^t_{m+2}\).
Therefore, substituting into (A.10) the ansatz \(v^t_j(\xi ) =\sum _{n=1}^\infty V^t_n B^t_n(\xi )\) with \(B^t_n\) from (A.11) and \(V^t_n\in {\mathbb {C}}\) where \(V^1_1 =V^{(-1)}_1\), the orthogonality condition on segments (A.10) takes the form
Due to Lemma A.3 this equality cannot be satisfied with nonzero set of coefficients \(\{V^t_1,\dots ,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(\xi ) =\sum _{n=m+2}^\infty V^t_n B^t_n(\xi )\) are needed to fulfill (A.10). This proves that \(X^\perp _h\) is spanned by piecewise polynomials of the local degree \(n\ge m+2\). Consequently, the piecewise polynomials of the local degree \(n\le m+1\) are dense in \(Z_h\), and the assertion of the theorem follows. \(\square \)
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 \({\mathbb {R}}^d\), \(d\in {\mathbb {N}}\). Further we define \({\mathbb {P}}_n\) as the space of d-dimensional polynomials of degree at most n with respect to every coordinate of the d-dimensional variable \(\xi =(\xi _1,\dots ,\xi _d)\), namely,
Lemma A.5
For fixed \(m\ge 0\), the d-dimensional Gegenbauer polynomials \(G_n(\xi ) :=\prod _{i=1}^d G_{n_i}(\xi _i) \in {\mathbb {P}}_n([-1,1]^d)\) of degree \(n\ge m+3\) are orthogonal to arbitrary polynomials \(P\in {\mathbb {P}}_m([-1,1]^d)\) with respect to inner products induced by the \(L^2\)-norm and the \(H^1\)-seminorm as follows
for all \(n\ge m+3\). Moreover, using the notion (for \(i=1,\dots ,d\))
for every \(t_1,\dots ,t_d\in \{-1,1\}\), there exist exactly \((m+2)^d-(m+1)^d\) nontrivial linear combinations of the shape functions with factors \(\{c^{(t_1,\dots ,t_d)}_{n_1\dots n_d} \}_{n_1,\dots ,n_d=1}^{m+2} \in {\mathbb {R}}^{(m+2)^d}\),
which satisfy the equality
Proof
Since the multiple integral over \([-1,1]^d\) can be split into the product of d single integrals over \([-1,1]\), we infer from (A.4) that
if there exists an index \(i\in \{1,\dots ,d\}\) such that \(n_i\ge m+3\) and \(m_i\le m\). Hence (A.12) holds.
We note that from (A.12) it follows that (A.15) holds with \(w^{(t_1,\dots ,t_d)}=G_n\) for all \(n\ge m+3\) and arbitrary \(\kappa \).
Now we take \(w^{(t_1,\dots ,t_d)} \in {\mathbb {P}}_{m+2}([-1,1]^d)\) in the form of (A.14) and a polynomial \(P\in {\mathbb {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,\dots ,t_d)\in \{1,1\}^d\) we get \((m+1)^d\) homogeneous equations for \((m+2)^d\) unknown coefficients \(c^{(t_1,\dots ,t_d)}_{n_1\dots n_d}\in {\mathbb {R}}\) by \((n_1,\dots ,n_d)\in \{1,\dots ,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,\dots ,t_d)}\in H^1((-1,1)^d)\) which do not fulfill (A.15) has \((m+1)^d\) degrees of freedom. The proof is complete.
\(\square \)
Based on Lemma A.5 we prove the main result generalizing Theorem A.4 to the multidimensional case.
Theorem A.6
For a uniform polyhedral mesh \(M_h\) of the size \(h>0\) over the computational domain \(\Omega _h\subset {\mathbb {R}}^d\), \(d\in {\mathbb {N}}\), let \(X_h\) be the finite dimensional set of continuous (over \(\Omega _h\)) and piecewise (on \(M_h\)) polynomials over \({\mathbb {C}}\) of degree at most m in every coordinate, with \(m\ge 1\) fixed. For given \(k\in {\mathbb {R}}\), the orthogonal complement
can be spanned by continuous (over \(\Omega _h\)) and piecewise (on \(M_h\)) polynomials of degree \(n\ge m+2\). The remaining polynomials of degree \(n\le m+1\) are dense in \(Z_h\).
Proof
We start by defining the transformation of functions depending on \(x\in \Omega _h\) to the local coordinates \(\xi \in [-1,1]^d\). Let \(N_h=\{x^j\}_{j=1}^N\) be the set of nodal points of the polyhedral mesh. We associate every point \(x^j\), \(j=1,\dots ,N\), to the patch \(\Pi ^j\) of all neighbor elements containing \(x^j\). These neighboring elements consist of \(2^d\) polyhedra \(Q^j_{(t_1,\dots ,t_d)}\) indexed by the binary numbers \((t_1,\dots ,t_d)\in \{-1,1\}^d\). The multi-index \(t=(t_1,\dots ,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 \(\Pi ^j =\{x_i\in {\mathbb {R}}: x^j_i-h\le x_i\le x^j_i+h\}_{i=1}^d\) associates the stencil of \(3^d\)-points and consists of the polyhedra
Using a partition of unity, every function \(v\in H^1(\Omega _h;{\mathbb {C}})\) can be partitioned by \(v=\sum _{j=1}^N v^j\) continuously over the family of overlapping patches \(\{\Pi ^j\}_{j=1}^N\) covering the computational domain \({\overline{\Omega }}_h\) in such a way that each part \(v^j(x)\in H^1_0(\Pi ^j;{\mathbb {C}})\) is supported on the reference patch \(\Pi ^j\). We transform each polyhedron \(Q^j_{(t_1,\dots ,t_d)}\) of the patch to the parent domain \([-1,1]^d\) with the linear coordinate transformation \(x_i=x^j_i+\frac{h}{2}(\xi _i-t_i)\), \(i=1,\dots ,d\). After the transformation we obtain \(2^d\) images \(v^j\bigl ( x^j +{\textstyle \frac{h}{2}}(\xi -t) \bigr ) =:v^{(t_1,\dots ,t_d)}_j(\xi )\in H^1((-1,1)^d;{\mathbb {C}})\) for \((t_1,\dots ,t_d)\in \{-1,1\}^d\) which are continuous across mutual boundaries (corresponding to \(\xi _i=t_i\), \(i=1,\dots ,d\)) and zero on the external boundary (corresponding to \(\xi _i=-t_i\), \(i=1,\dots ,d\)). Therefore, \(v^{(t_1,\dots ,t_d)}_j\) can be spanned by the Gegenbauer polynomial based shape functions
Here we use the notation (A.13), where \((t_1,\dots ,t_d)\in \{-1,1\}^d\) corresponds to \(Q^j_{(t_1,\dots ,t_d)}\).
Now we can represent the orthogonal complement \(X^\perp _h\) in local coordinates. After transformation of \(v^j\in H^1_0(\Pi ^j;{\mathbb {C}})\) to the parent domain, for all \(w^h\bigl ( x^j +{\textstyle \frac{h}{2}}(\xi -t) \bigr ) =:P(\xi ) \in {\mathbb {P}}_{m}([-1,1]^d;{\mathbb {C}})\) we express the respective integral equation in (A.16) as
Substituting \(v^{(t_1,\dots ,t_d)}_j =S^{(t_1,\dots ,t_d)}_{n_1\dots n_d}\) from (A.17), due to (A.12) in Lemma A.5 equation (A.18) is satisfied for \(n\ge m+3\) with \(n=\max _{i=1,\dots ,d} n_i\). Moreover, due to (A.15) in Lemma A.5, \((m+2)^d-(m+1)^d\) nontrivial linear combinations \(w^{(t_1,\dots ,t_d)} \in {\mathbb {P}}_{m+2}([-1,1]^d)\) of the shape functions from (A.14) satisfy equation (A.18). Therefore, the finite set of \((m+1)^d\) shape functions \(\{S^{(t_1,\dots ,t_d)}_{n_1\dots n_d}\}_{ \max _{i=1,\dots ,d} n_i=1}^{m+1}\) forms the span of those images which do not fulfill condition (A.18) in the parent domain.
After the inverse transformation of the images to the global coordinates we infer the assertion of the theorem. \(\square \)
In the particular case of \(m=1\), Theorem 3.2 follows from Theorem A.6.
Appendix B: The PGE approximation in 3d
Utilizing the spatial shape functions \(S^{(t_1,t_2,t_3)}_{n_1 n_2 n_3}(\xi )\), \(t_1,t_2,t_3\in \{-1,1\}\), \(n_1,n_2,n_3\in \{1,2\}\) for \(\xi =(\xi _1,\xi _2,\xi _3)\in [-1,1]^3\), which are defined in (3.6), we construct the symmetric basis functions corresponding to node, edge, facet, and bubble modes, respectively,
Here \(x=(x_1,x_2,x_3)\) live in the polyhedra \(Q^j_{(t_1,t_2,t_3)}\), \(t_1,t_2,t_3\in \{-1,1\}\), forming the patch \(\Pi ^j =Q^j_{(1,-1,-1)}\cup Q^j_{(-1,-1,-1)}\cup Q^j_{(1,1,-1)}\cup Q^j_{(-1,1,-1)}\cup Q^j_{(1,-1,1)}\cup Q^j_{(-1,-1,1)}\cup Q^j_{(1,1,1)}\cup Q^j_{(-1,1,1)} =\{x_i\in {\mathbb {R}}: x^j_i-h\le x_i\le x^j_i+h\}_{i=1}^3\) as defined in (3.8).
We use the finite element ansatz of the trial function
with unknown coefficients \(\{U^j\}_{j=1}^N\in {\mathbb {C}}^N\), and test functions for \(i=1,\dots ,N\)
The seven unknown parameters \(\alpha ^N_1, \alpha ^E_0, \alpha ^E_1, \alpha ^F_0, \alpha ^F_1, \alpha ^B_0, \alpha ^B_1 \in {\mathbb {R}}\) determine the weights by the modes in the test function \(v^h\) as
The corresponding system matrix \(A(\kappa ^2)\in \mathrm{Sym}(N^2)\) of the homogeneous Helmholtz equation consists of the entries \(a_{ij}(\kappa ^2)\), \(i,j=1,\dots ,N\), such that
according to the 27-point interior stencil \(\{(A_{\mathrm{stencil}}^{\mathrm{interior}}(\kappa ^2))_{i_1i_2i_3} \}_{i_1,i_2,i_3=1}^3\in {\mathbb {R}}^{3\times 3\times 3}\)
with the following stencil coefficients
Substituting for the trial function \(u^h=I_{X_h}u\) the plane wave solution
with arbitrary incident angle \(\theta =(\theta _1,\theta _2)\in (-\pi ,\pi ]\times [0,\pi ]\), into the discrete variational problem (3.10) leads to the dispersion equation
with factors
Further we expand these factors as \(\kappa \rightarrow 0\):
where \(C :=(C_2S_2)^2 +(C_1S_1)^2S_2^4\) and \(S :=(C_1S_1)^2C_2^2S_2^4\). Inserting these formula together with (B.2) and (B.3) into (B.4), the dispersion equation has the expansion
where the residual term \(\mathrm{Er}(\,\cdot \,,\theta ) =\mathrm{o}(\kappa ^7)\) for all \(\theta \in (-\pi ,\pi ]\times [0,\pi ]\). In the above expression the six factors in parenthesis become zero for the following choice of parameters:
depending on one free parameter \(\alpha ^B_1\in {\mathbb {R}}\). This results in the asymptotic relations
Inserting the weights (B.5) into (B.3) we obtain the stencil coefficients as
Setting \(\kappa =0\) in (B.7), which corresponds to the Laplace operator, the consistency conditions \(A_0(0)>0\) and \(A_0(0)+2A_1(0)+A_2(0)=0\) hold. Setting \(\alpha ^B_1 =-\frac{3741}{736}\approx -5.08\) we get \(A_0(\kappa ^2)>0\) for all \(\kappa ^2\). Therefore, Theorem 3.4 providing the existence of the discrete variational solution to the Helmholtz problem, and Theorem 3.5 estimating the error of the discretization, hold in the 3d case as well.
Rights and permissions
Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made.
About this article
Cite this article
Kovtunenko, V.A., Kunisch, K. Revisiting generalized FEM: a Petrov–Galerkin enrichment based FEM interpolation for Helmholtz problem. Calcolo 55, 38 (2018). https://doi.org/10.1007/s10092-018-0280-5
Received:
Accepted:
Published:
DOI: https://doi.org/10.1007/s10092-018-0280-5
Keywords
- Generalized FEM
- Petrov–Galerkin enrichment
- FEM interpolation
- Wavelets basis
- Gegenbauer polynomials
- Helmholtz equation