FEM with Trefftz trial functions on polyhedral elements

Journal of Computational and Applied Mathematics 263 (2014) 202–217

Contents lists available at ScienceDirect

Journal of Computational and Applied Mathematics journal homepage: www.elsevier.com/locate/cam

FEM with Trefftz trial functions on polyhedral elements Sergej Rjasanow, Steffen Weißer ∗ Applied Mathematics, Saarland University, Campus, 66041 Saarbrücken, Germany

article

info

Article history: Received 22 October 2012 Received in revised form 16 July 2013 MSC: 65N30 65N38 41A25 41A30 Keywords: BEM-based FEM Polyhedral mesh Convergence estimates Non-standard finite element method

abstract The goal of this paper is to generalize a variant of the BEM-based FEM for second order elliptic boundary value problems to three space dimensions. Here, the emphasis lies on polyhedral meshes with polygonal faces, where even non-convex elements are allowed. Due to an implicit definition of the trial functions in the spirit of Trefftz, the strategy yields conforming approximations and is very flexible with respect to the meshes. Thus, it gets into the line of recent developments in several areas. The arising local problems are treated by two dimensional Galerkin schemes coming from finite and boundary element formulations. With the help of a new interpolation operator and its properties, convergence estimates are proven in the H 1 - as well as in the L2 -norm. Numerical experiments confirm the theoretical results. © 2013 Elsevier B.V. All rights reserved.

1. Introduction The finite element method (FEM) is a powerful tool for the approximation of solutions of boundary value problems. Its history goes back to the mid of the twentieth century. Since that time, there have been a lot of developments for the method which is usually applied on simplicial meshes. Nowadays, the need for more general meshes is arising and thus schemes like finite volume [1], discontinuous Galerkin [2], multiscale finite element [3] or mimetic discretization methods [4] as well as virtual element methods [5] are sometimes favorable over the classic finite element approach. Other strategies like the Trefftz-discontinuous Galerkin method [6,7] or the ultra weak variational formulation [8] already make use of Trefftz trial functions. In contrast to the forthcoming presentation, however, they use non-conforming approximation spaces. All these strategies are applicable on more general meshes due to their nature. But also within the finite element method, there are some developments towards polygonal and polyhedral meshes, see [9–11]. Such meshes appear naturally in geological and biological science but also while meshing complex geometries where simplicial elements can be restrictive and deteriorate the mesh quality. In 2009, D. Copeland, U. Langer and D. Pusch [12] proposed a new strategy which goes back to boundary element domain decomposition methods [13]. This approach uses locally implicitly defined trial functions in the finite element method. These functions fulfill element-by-element the underlying differential equation and thus, we also refer to them as Trefftz trial functions. Due to this definition, the method is applicable on general polygonal meshes in two space dimensions. If the trial functions in the finite element method are treated with boundary element methods (BEM), as in the following, we call the strategy BEM-based FEM. Several variants of this promising approach have been studied concerning convergence [14–16], higher order trial functions [17] as well as in an adaptive strategy [18]. Furthermore, the studies have been extended to convection–diffusion–reaction problems [19] and H (div)-conforming approximations [20]. The main results have been

∗

Corresponding author. Tel.: +49 681 302 3811. E-mail address: [email protected] (S. Weißer).

0377-0427/$ – see front matter © 2013 Elsevier B.V. All rights reserved. http://dx.doi.org/10.1016/j.cam.2013.12.023

S. Rjasanow, S. Weißer / Journal of Computational and Applied Mathematics 263 (2014) 202–217

203

collected in two Ph.D. theses [21,22]. Already in [23], the method was formulated in three space dimensions on polyhedral meshes and this formulation has been used in the literature cited above. In this formulation, it is essential to discretize the union of all surfaces of the polyhedral elements by triangles. Afterwards, each node in the triangulation corresponds to a degree of freedom in the finite element computations. The aim of this presentation is to generalize the lowest order trial functions to arbitrary polyhedral elements with polygonal faces such that only the nodes from the polyhedral mesh produce degrees of freedom. We adapt the ideas in [17] and generalize them to three space dimensions. Compared with [17], we even relax the convexity of the elements and prove linear convergence in the H 1 -norm and quadratic convergence in the L2 -norm with respect to the mesh size. The article is organized as follows. In Section 2, we state the model problem, give a definition of regular and stable polyhedral meshes and review the boundary element method which is utilized in the numerical realization. The trial functions are introduced in Section 3 and we discuss the setup of the finite element matrix. A priori error estimates are proven in Section 4. Finally, we present numerical experiments in Section 5 which confirm the theoretical results, and we draw some conclusions in Section 6. 2. Preliminaries 2.1. Model problem For the sake of simplicity, we restrict ourselves to a model problem. Let Ω ⊂ R3 be a polyhedral and bounded domain with boundary Γ = ΓD ∪ ΓN which is split into a Dirichlet and a Neumann part. Furthermore, we assume |ΓD | > 0 and there is a given source term f ∈ L2 (Ω ), a Dirichlet datum gD ∈ H 1/2 (ΓD ) as well as a Neumann datum gN ∈ L2 (ΓN ). We consider the boundary value problem

−div(a∇ u) = f in Ω , u = gD on ΓD , a∇ u · n = gN on ΓN ,

(1)

for a ∈ L∞ (Ω ) with 0 < amin ≤ a ≤ amax on Ω . By the use of an extension uD ∈ H 1 (Ω ) of the Dirichlet datum, the Galerkin formulation for (1) with the solution u = u0 + uD ∈ H 1 (Ω ) reads Seek u0 ∈ V : aΩ (u0 , v) = ℓ(v)

∀v ∈ V

(2)

with aΩ (u, v) =

 a∇ u · ∇v Ω

and

ℓ(v) = (f , v)Ω + (gN , v)ΓN − aΩ (uD , v), where V = HD1 (Ω ) = v ∈ H 1 (Ω ) : γ0 v = 0 on ΓD .





Here, γ0 : H 1 (Ω ) → H 1/2 (Γ ) denotes the usual trace operator, see [24], and (·, ·)Ω and (·, ·)ΓN denote the L2 -scalar products over Ω and ΓN , respectively. Due to the properties of the material coefficient a ∈ L∞ (Ω ) the bilinear form aΩ (·, ·) : H 1 (Ω ) × H 1 (Ω ) → R is bounded and coercive on V , i.e. there are constants m, M > 0 such that

|aΩ (u, v)| ≤ M ∥u∥H 1 (Ω ) ∥v∥H 1 (Ω ) and aΩ (v, v) ≥ m∥v∥2H 1 (Ω ) ,

(3)

for all u, v ∈ V ⊂ H (Ω ). Since ℓ(·) : H (Ω ) → R is also bounded on V , the variational formulation (2) admits a unique solution according to the √ lemma of Lax–Milgram, see [25]. Due to the boundedness and coercivity of the bilinear form, the energy norm ∥ · ∥E = aΩ (·, ·) is equivalent to the usual Sobolev norm ∥ · ∥H 1 (Ω ) on V . 1

1

2.2. Polyhedral meshes The discretization is based on a finite partition of the domain Ω into non-overlapping polyhedral subdomains K . The collection of these subdomains is called mesh and labeled Kh . Besides the polyhedral elements K ∈ Kh , we additionally have polygonal faces F ∈ Fh , edges E ∈ Eh and nodes z ∈ Nh in the discretization. The sets of faces, edges and nodes are labeled Fh , Eh and Nh , respectively. With Nh,D , we indicate the nodes which lie on ΓD or on the transition lines between the Dirichlet and Neumann boundary. The sets of nodes, edges and faces which belong to an element K ∈ Kh are labeled N (K ), E (K ) and F (K ), respectively. In the same manner, we define N (F ) and E (F ) for F ∈ Fh as the sets of the corresponding nodes and

204

S. Rjasanow, S. Weißer / Journal of Computational and Applied Mathematics 263 (2014) 202–217

edges. The nodes, edges, faces and elements are d-dimensional objects with d = 0, 1, 2, 3, respectively. For d ≥ 1, they are open sets and their boundaries decompose into (d − 1)-dimensional objects. Furthermore, we assume that the interior of a d-dimensional object does not contain any object of dimension less or equal than d. Consequently, by inserting a hanging node on a straight edge, the edge is split into two new ones. Therefore, it is possible to have, in the sense of classical FEM, hanging nodes on edges and hanging edges on faces. However, it is not allowed to have a hanging node in the interior of a face unless additional edges are introduced. The diameter of an element K ∈ Kh , a face F ∈ Fh and the length of an edge E ∈ Eh are denoted by hK , hF and hE , respectively. Definition 1. A set of faces Fh is called regular if all faces are flat and there is a constant σF such that every face F ∈ Fh is star-shaped with respect to a circle inscribed in F with radius ρF and center zF and the aspect ratio hF /ρF is uniformly bounded from above by σF for all F ∈ Fh . Definition 2. A mesh Kh is called regular if the associated set of faces Fh is regular and if there is a constant σK such that every element K ∈ Kh is star-shaped with respect to a ball inscribed in K with radius ρK and center zK and the aspect ratio hK /ρK is uniformly bounded from above by σK for all K ∈ Kh . Definition 3. A mesh Kh is called stable if there is a constant cK such that for every K ∈ Kh the estimate hK ≤ cK hE for all E ∈ E (K ) holds true. A regular mesh which only consists of tetrahedral elements is consequently shape regular in the common sense of classical finite element methods, see [25]. In the following, we always consider sequences of regular and stable polyhedral meshes Kh , where the constants in the definitions above hold uniformly. The stability ensures hE ≤ hF ≤ hK ≤ cK hE ≤ cK hF for K ∈ Kh and all F ∈ F (K ) and E ∈ E (F ). Consequently, the number of nodes, edges and faces per element is uniformly bounded for a sequence of regular and stable polyhedral meshes. In this presentation the regularity and stability conditions on the mesh are more restrictive than the assumptions in [16]. However, we have chosen them in such a way that the notion of admissible elements is rather obvious to the reader. The analysis in this paper can be generalized quite easily to the situation in [16]. 2.3. A review of the boundary element method An alternative approach to the finite element method is the so called boundary element method. An introduction to this approach can be found in [26]. In the case of a constant material coefficient a(·) and vanishing right hand side f in (1), we end up with the Laplace equation. We consider this equation on an arbitrary element K ∈ Kh and prescribe some Dirichlet data on the boundary, i.e.

−1u = 0 in K , u = g on ∂ K .

(4)

With the help of the so called fundamental solution of minus the Laplacian, which is U ∗ (x, y) =

1 4π |x − y|

for x, y ∈ R3 ,

the usual trace operator γ0K : H 1 (K ) → H 1/2 (∂ K ) as well as the conormal derivative, which takes the form

γ1K v = nK · γ0K ∇v ∈ H −1/2 (∂ K ) for sufficiently regular v and the outer unit normal vector nK to ∂ K , we have the representation formula   ∗ K u(x) = U (x, y)γ1 u(y) dsy − γ1K,y U ∗ (x, y)γ0K u(y) dsy for x ∈ K . ∂K

(5)

∂K

Applying the trace and the conormal derivative operators to this formula, it is possible to derive the relation between the Dirichlet datum γ0K u and the Neumann datum γ1K u on ∂ K . It holds 1 γ1K u = SK γ0K u with SK = V− K

1 2

I + KK ,



(6)

where the operator SK : H 1/2 (∂ K ) → H −1/2 (∂ K ) is called a Steklov–Poincaré operator. Here, we have used the standard boundary integral operators which are well studied, see e.g. [27,26]. The single-layer potential operator VK and the double-layer potential operator KK are given as

(VK ζ )(x) = for x ∈ ∂ K .

 ∂K

U ∗ (x, y)ζ (y) dsy

and (KK ξ )(x) = lim

ε→0

 y∈∂ K :|y−x|≥ε

γ1K,y U ∗ (x, y)ξ (y) dsy

S. Rjasanow, S. Weißer / Journal of Computational and Applied Mathematics 263 (2014) 202–217

205

Let Tl = Tl (∂ K ) be an admissible triangulation of ∂ K and denote its elements, the triangles, by T . The set of nodes is labeled Ml = Ml (∂ K ). The index l stands for a discretization parameter. For each T ∈ Tl , we define a function

τT0 =



1, 0,

in T else

and set Φl0 = {τT0 : T ∈ Tl },

which is a basis of the space of piecewise constant functions over the mesh Tl . It is span Φl0 ⊂ H −1/2 (∂ K ), and thus we use this discrete space for the approximation of the Neumann trace. Additionally, we define for z ∈ Ml a function 1, ϕ = linear, 0,



1 z

at z on T ∈ Tl at x ∈ Ml \ {z }

and set Φl1 = {ϕz1 : z ∈ Ml }.

(7)

These functions form a basis of the space of piecewise linear and globally continuous functions over the boundary ∂ K . Therefore, it holds that span Φl1 ⊂ H 1/2 (∂ K ) and we use this discrete space for the approximation of the Dirichlet trace. In the boundary value problem (4), the Dirichlet datum g ∈ H 1/2 (∂ K ) is given and we approximate it by gl ∈ span Φl1 . This datum is used in a discrete Galerkin formulation for (6) Seek tl ∈ span Φl0 : (VK tl , ϑ)L2 (∂ K ) =

 1 2

I + KK gl , ϑ





L2 (∂ K )

∀ϑ ∈ Φl0

to approximate the unknown Neumann datum t = γ1K u ∈ H −1/2 (∂ K ). Due to the properties of the boundary integral operators, the variational formulation has a unique solution. The representations tl =



tτ τ

and



gl =

gϕ ϕ

ϕ∈Φl1

τ ∈Φl0

yield the system of linear equations VK ,l t l =

1 2

MK ,l + KK ,l g ,



l

where the underline refers to the coefficient vector, e.g. t l = (tτ )τ ∈Φ 0 . l

The approximations gl ∈ span Φl1 and tl ∈ span Φl0 can be used in the representation formula (5) to obtain an approximation of the exact solution u of (4) in K . Because of the piecewise linear data, an analytic formula exists to evaluate the boundary integrals in (5). In the implementation, we realize the entries of the boundary element matrices by the use of a semi analytic integration scheme. The inner integral is evaluated analytically and the outer one is approximated by Gaussian quadrature. More details on this and on an efficient treatment of the matrices can be found in [28]. 3. BEM-based FEM in 3D In this section, we discuss a generalization of a variant of the BEM-based finite element method studied for the two dimensional case in [17,18]. So, we address the definition of trial functions on meshes with polyhedral elements which satisfy the regularity of Definition 2. These functions are used to construct an approximation space Vh which can be utilized in the discrete Galerkin formulation of the finite element method. The idea of the BEM-based FEM is to define the trial functions implicitly on each element as local solutions of the underlying differential equation and to treat the local problems by boundary element methods. Therefore, the material coefficients in the differential equation are approximated by piecewise constants and the right hand side is neglected, see [18]. For our model problem, there is only one material coefficient and we end up with the Laplace equation. For other boundary value problems like convection–diffusion–reaction, for example, there might be more than one material coefficient. 3.1. Generalization of the trial functions In order to get a nodal basis of Vh , we declare for each node z ∈ Nh a trial function ψz which is equal to one in z and zero in all other nodes of the mesh. In [17,18], the authors used linear data on the edges for the two dimensional case. To be more precise, they defined ψz as a unique solution of

−1ψz = 0 in K for all elements K ,  1 for x = z ψz (x) = 0 for all other nodes x, ψz is linear on each edge of K , where the elements K are convex polygons. In the case of polyhedral elements with triangulated surfaces, we can get a straightforward generalization by replacing the word ‘ edge’ by ‘triangular face’. This strategy is chosen in the present literature [23,15,14,16].

206

S. Rjasanow, S. Weißer / Journal of Computational and Applied Mathematics 263 (2014) 202–217

Fig. 1. Stepwise construction of trial functions.

To handle arbitrary polygonal faces of the polyhedral elements, we follow another strategy. If we look again into the two dimensional case, we observe that the values of the trial functions are fixed in the nodes and extended uniquely along the edges by linear functions. This linear extension is nothing other than a harmonic extension along the edge and thus the trial function is also defined on the edges according to the underlying differential equation. Therefore, we propose a stepwise construction for the trial functions in the case of polyhedral elements with polygonal faces as sketched in Fig. 1. A similar idea has been used in two dimensions for the construction of multiscale finite elements in [29]. Denoting the i-dimensional Laplace operator by ∆i , we define the trial function ψz , which is related to z ∈ Nh , as a unique solution of

−∆3 ψz = 0 in K for all K ∈ Kh , −∆2 ψz = 0 in F for all F ∈ Fh , −∆1 ψz = 0 in E for all E ∈ Eh ,  1 for x = z , ψz (x) = 0 for x ∈ Nh \ {z }, where the Laplace operators have to be understood in the corresponding linear parameter spaces. The values in the nodes are prescribed. Afterwards, we solve a Dirichlet problem for the Laplace equation on each edge. Then, we use the computed data as Dirichlet datum for the Laplace problem on each face, and finally we proceed with the Laplace problem on each element, where the solutions on the faces are used as boundary values. In the case of convex faces and elements, these problems are understood in the classical sense and we have ψz ∈ C 2 (K ) ∩ C 0 (K ). In the more general situation of non-convex elements, the weak solution is considered such that we have at least ψz ∈ H 1 (K ). The novelty of the present paper concerning the BEM-based FEM is to define the trial functions implicitly on the edges, faces, and elements and not only, as in previous publications, on the element level. In particular, considering a Helmholtz or convection–diffusion equation, the trial functions already become non-linear on the edges and the properties of the differential equation are built into the trial functions on all levels, on the edges, faces, and elements. In the case of the Helmholtz equation one has to care about the eigenvalues of −∆. However, for a sufficiently small mesh size everything works well, see [23]. Remark 1. For the model problem (1), the data of the trial functions along the edges is linear. This might limit the efficiency of the BEM-based FEM compared to other Trefftz-type approaches when it is applied to large elements, cf. [6]. However, in contrast to Trefftz-discontinuous Galerkin methods, this choice of data along the edges ensures a conforming approximation space. Furthermore, the presented approach is of first order and can be extended with additional trial functions to obtain better approximations, see [17]. 3.2. Discrete Galerkin formulation We set

Ψh = {ψz : z ∈ Nh } as well as Ψh,D = {ψz : z ∈ Nh,D }, and we introduce the trial space Vh = span Ψh,0

with Ψh,0 = Ψh \ Ψh,D .

According to the regularity of the trial functions inside the elements and their continuity across the element boundaries, the trial space is conforming in the sense of Vh ⊂ V . The discrete version of the Galerkin formulation (2) for the approximation uh of the exact solution u with uh = uh,0 + uh,D and uh,0 =

 ψ∈Ψh,0

βψ ψ ∈ Vh as well as uh,D =

 ψ∈Ψh,D

βψ ψ

S. Rjasanow, S. Weißer / Journal of Computational and Applied Mathematics 263 (2014) 202–217

207

Fig. 2. Polyhedral element with surface triangulations of level l = 0, 1, 2.

reads



βψ aΩ (ψ, φ) = (f , φ)Ω + (gN , φ)ΓN −



βψ aΩ (ψ, φ) for φ ∈ Ψh,0 ,

ψ∈Ψh,D

ψ∈Ψh,0

where uh,D is the discrete extension of the boundary data gD in the model problem (1). In the representations of uh,0 and uh,D the same symbol βψ is used for the coefficients, but there should be no confusion since Ψh,0 ∩ Ψh,D = ∅. 3.3. Construction of the trial functions If there are only triangular faces of the polyhedra, the definition of the trial space agrees with the former generalization. Otherwise, we have to solve additional Dirichlet boundary value problems on the faces of the elements. This can be done by the help of a two dimensional boundary element method. Nevertheless, we propose to use a 2D finite element method on triangular meshes of the faces. For each F ∈ Fh , we introduce a mesh Tl (F ) of level l. The coarsest mesh with l = 0 is obtained by connecting the nodes z ∈ N (F ) with the point zF , which is fixed once per face according to Definition 1. Afterwards, the meshes of level l ≥ 1 are defined recursively by splitting each triangle of the previous level into four similar triangles, see Fig. 2. The set of nodes in the triangular mesh is denoted by Ml (F ). Obviously, the discretizations of the faces can be combined to a whole conforming surface mesh of an element K ∈ Kh by setting



Tl (∂ K ) =



Tl (F ) and Ml (∂ K ) =

F ∈F (K )

Ml (F ).

F ∈F (K )

Now, we address the approximation of the trace of a trial function ψz , z ∈ Nh on a face F ∈ Fh with z ∈ N (F ). For the finite element computations on the face, we have to use a further trial space. Here, we utilize the usual basis Φl1 (F ) of piecewise linear and continuous functions. It holds

Φl1 (F ) = {ϕz1 : z ∈ Ml (F )}, cf. (7). We approximate the trace of the trial function ψz by gzF =



gϕ ϕ,

(8)

ϕ∈Φl1 (F )

where the coefficients gϕ which are related to ϕ = ϕx with x ∈ ∂ F are fixed such that gzF coincides with the piecewise linear data of ψz on the edges of the face F . Consequently, we obtain the discrete Galerkin formulation Seek gzF :





∇ϕx · ∇ϕy = 0,

gϕx

x∈Ml (F )

∀y ∈ Ml (F ) : y ̸∈ ∂ F ,

(9)

F

for the approximation of the trace of ψz on F ∈ Fh , where z ∈ N (F ). This discrete variational formulation admits a unique solution and the corresponding system of linear equations can be solved by the Cholesky decomposition or the conjugate gradient method, for example. Changing the level of the face discretization, the accuracy of the finite element approximation can be adapted. Remark 2. In the case of another differential equation like for Helmholtz or convection–diffusion problems, we should choose a mesh level l ≥ 1 for the faces. This automatically implies a discretization of the edges which can be used to handle the non-linear data of the trial functions along the edges in the 2D finite element method on the faces. To get an approximation gzK of the trace of ψz on the whole boundary of an element K ∈ Kh , the Dirichlet problems on the faces F ∈ F (K ) are solved successively. Thus, we have an approximation gzF of γ0K ψz on each face F of K . Afterwards, these approximations on the faces are combined to an approximation over the whole surface mesh by setting gzK =

 ϕ∈Φl1 (∂ K )

gϕ ϕ

with Φl1 (∂ K ) =

 F ∈F (K )

Φl1 (F ),

208

S. Rjasanow, S. Weißer / Journal of Computational and Applied Mathematics 263 (2014) 202–217

where the coefficients gϕ are taken from (8). Here, Φl1 (∂ K ) is the nodal basis of piecewise linear and continuous functions over the surface mesh of K of level l. This basis has been used already for the boundary element method in (7). Due to the prescribed data along the edges, we obtain a piecewise linear and globally continuous approximation of the Dirichlet trace γ0K ψz over the surface triangulation Tl (∂ K ). Instead of working with the exact trial functions ψz , we introduce their approximations ψzl as unique solutions of

−1ψzl = 0 in K , ψzl = gzK on ∂ K for all K ∈ Kh , where the piecewise linear Dirichlet datum is used. Additionally, we set Vhl = span Ψhl,0

with Ψhl,0 = Ψhl \ Ψhl,D ,

where

Ψhl = {ψzl : z ∈ Nh } as well as

Ψhl,D = {ψzl : z ∈ Nh,D }.

Thus, we obtain another conforming approximation space Vhl ⊂ V which approximates Vh and depends on the auxiliary triangulations of the faces and especially on their discretization level l. 3.4. Realization and setup of the Galerkin formulation In the computational realization, we actually solve the discrete Galerkin method which is obtained by exchanging Vh , Ψh,0 and Ψh,D by Vhl , Ψhl,0 and Ψhl,D in the discrete variational formulation above. That means, we seek an approximation ulh of the exact solution u with ulh = ulh,0 + ulh,D and ulh,0 =

 ψ∈Ψh,0

βψ ψ ∈ Vhl as well as ulh,D =



βψ ψ

(10)

ψ∈Ψhl ,D

fulfilling



βψ aΩ (ψ, φ) = (f , φ)Ω + (gN , φ)ΓN −

ψ∈Ψhl ,0



βψ aΩ (ψ, φ) for φ ∈ Ψhl,0 ,

(11)

ψ∈Ψhl ,D

where ulh,D is a discrete extension of the boundary data gD . 3.4.1. L2 -scalar product over ΓN The L2 -scalar product over the Neumann boundary is rather simple. The trial functions are given as piecewise linear on a triangulation of ΓN . Therefore, a standard Gaussian quadrature can be used on each triangle. 3.4.2. L2 -scalar product over Ω To handle the L2 -scalar product over Ω , we use a quadrature over polyhedral elements. This can be done by refining the element into tetrahedra or by using directly a numerical integration scheme for polyhedral domains, see [30]. The evaluations of the trial functions are realized with the help of the representation formula (5). Thus, the Neumann traces of the trial functions have to be approximated. 3.4.3. Bilinear form and approximation of Neumann traces At this point, it becomes clear why we have chosen a finite element discretization of the faces. The surface mesh Tl (∂ K ) and the trial functions on the boundary ∂ K for the two dimensional finite element method fit into the theory of the boundary element method reviewed in Section 2. Following the ideas given there, we obtain an approximation tzK of the Neumann trace γ1K ψz in the form tzK =



tτ τ

  K 1 1 with t Kz = V− K ,l 2 MK ,l + KK ,l g . z

(12)

τ ∈Φl0 (K )

Here, we mention that the boundary element matrices have to be computed only once per element. Inverting VK ,l with an efficient LAPACK routine, these matrices are used for the setup of all Neumann traces of trial functions ψz . This can be done in a preprocessing step and, since the BEM matrices are independent from element to element, the procedure can be performed completely in parallel. For each element K ∈ Kh , the approximations of the Dirichlet and Neumann traces of the trial functions ψz with z ∈ N (K ) are gathered in the matrices



DK = g K z

 z ∈N (K )

and

 

NK = t Kz

z ∈N (K )

S. Rjasanow, S. Weißer / Journal of Computational and Applied Mathematics 263 (2014) 202–217

209

such that each column corresponds to the datum of one ψz . The relation (12) turns into

  1 1 NK = V− K ,l 2 M K ,l + K K ,l D K . Since we have the approximations of the Dirichlet and Neumann data, we can use the representation formula (5) to get an approximation of the trial functions inside the elements, if needed. Finally, we consider the approximation of the bilinear form aΩ (·, ·). We assume that the material coefficient is constant on each element such that a(x) = aK

for x ∈ K and K ∈ Kh ,

or it is approximated by a piecewise constant function. For ψ, φ ∈ Ψhl , it holds that aΩ (ψ, φ) =



 ∇ψ · ∇φ =

aK K

K ∈Kh

  aK  K ∈Kh

2

∂K

γ ψγ φ + K 1

K 0



 ∂K

γ φγ ψ K 1

K 0

according to Green’s first identity and since ψ as well as φ are harmonic on each element K ∈ Kh . Obviously, there are x, z ∈ Nh such that ψ = ψx and φ = ψz and if there is no K ∈ Kh with x, z ∈ N (K ), we have aΩ (ψ, φ) = 0. On the other hand, for each K ∈ Kh with x, z ∈ N (K ) we get the approximations gxK , gzK for the Dirichlet traces and txK , tzK for the Neumann traces of the trial functions ψx and ψz out of the matrices DK and NK , respectively. For the integrals in the representation of the bilinear form above, we choose the approximation

 ∂K

γ ψ γ ψz ≈ K 1

K x 0

 ∂K

txK gzK =



tτ

τ ∈Φl0 (K )

 ⊤



gϕ (τ , ϕ)L2 (∂ K ) = t Kx

MK ,l g K z

ϕ∈Φl1 (K )

with the mass matrix MK ,l defined as in Section 2. This yields the symmetric approximation aΩ (ψx , ψz ) ≈

 aK    ⊤ K K ∈Kh

2

tx

 ⊤

MK ,l g K + t Kz z



M K ,l g K ,

(13)

x

where the coefficient vectors are identical to zero if x ̸∈ N (K ) and z ̸∈ N (K ), respectively. Consequently, we use the local matrix NK⊤ MK ,l DK ∈ R|N (K )|×|N (K )| for each K ∈ Kh to set up the global finite element matrix. Due to the symmetric approximation (13) of the bilinear form, we obtain a symmetric system matrix in the finite element method which is sparse and positive definite. Therefore, the conjugate gradient method can be applied to get the solution of the system of linear equations. We can also go without the symmetrization step for the bilinear form aΩ (·, ·). However, in this case we end up with a non-symmetric matrix and we have to use an appropriate solver. 4. Convergence estimates Let u = u0 + uD be the solution of the model problem obtained by the variational formulation (2) and ulh = ulh,0 + ulh,D its Galerkin approximation which is gained by (10) and (11) with the approximated trial space Vhl . We assume that the extension of the Dirichlet datum gD can be chosen in such a way that uD = ulh,D ∈ span Ψhl,D . In this section, our goal is to introduce an interpolation operator Il : H 2 (Ω ) → span Ψhl such that

∥v − Il v∥H 1 (Ω ) ≤ ch |v|H 2 (Ω ) for v ∈ H 2 (Ω ),

(14)

where h = max{hK : K ∈ Kh }. In particular, the constant c may depend on the mesh level l of the surface triangulation, but it is especially independent of the mesh size h. The definition of an interpolation operator with this approximation property is a non-trivial task. However, it is an important tool for theoretical considerations and it enables us to prove the following convergence estimates. Theorem 1. Let Kh be a regular and stable mesh of a bounded polyhedral domain Ω ⊂ R3 . Then, it holds

∥u − ulh ∥H 1 (Ω ) ≤ ch |u|H 2 (Ω ) for u ∈ H 2 (Ω ), where the constant c only depends on the regularity and stability parameters of the mesh, on the level l and on m and M, see (3). Proof. Since Vl ⊂ V , the statement follows directly from Céa’s lemma and (14), see [25].



For the estimate in the theorem, we have implicitly assumed that there are no errors in the evaluation of the bilinear form and the right hand side of the Galerkin formulation. The case of approximate data f , gD , gN in the discrete problem can be treated in the usual way, where the Strang lemma is used instead of Céa’s lemma, see [25]. The change of the norm in which the finite element error is measured gives an additional power of h if the adjoint variational formulation admits a sufficiently regular solution.

210

S. Rjasanow, S. Weißer / Journal of Computational and Applied Mathematics 263 (2014) 202–217

Theorem 2. Let Kh be a regular mesh of a bounded polyhedral domain Ω ⊂ R3 . Under the condition that for any g ∈ L2 (Ω ) there is a unique solution of Seek w ∈ V : aΩ (v, w) = (g , v)Ω ,

∀v ∈ V ,

with w ∈ H (Ω ) such that 2

|w|H 2 (Ω ) ≤ C ∥g ∥L2 (Ω ) , it holds

∥u − ulh ∥L2 (Ω ) ≤ ch2 |u|H 2 (Ω ) for u ∈ H 2 (Ω ), where the constant c only depends on the regularity and stability parameters of the mesh, on the level l as well as on m and M, see (3), and on C . Proof. This result follows by the use of the Aubin–Nitsche trick and Theorem 1 in the classical way, see [25].



The rest of this section is devoted to give an appropriate interpolation operator. Since H (Ω ) ⊂ C (Ω ) according to the Sobolev embedding theorem, see [31], the pointwise evaluation of a function v ∈ H 2 (Ω ) is well defined and we can utilize a nodal interpolation operator. We set 2



Il v =

0

αψ ψ with αψzl = v(z ) for z ∈ Nh .

ψ∈Ψhl

As in [17], we start studying the properties of the interpolation operator over one element K ∈ Kh . Consequently, we are just interested in the linear combination of trial functions which are related to nodes in N (K ). Lemma 1. The restriction of the interpolation operator Il to an arbitrary element K ∈ Kh fulfills Il p = p on K for each linear polynomial p ∈ P 1 (K ). Proof. Let p ∈ P 1 (K ). Obviously, the restriction of p to each face F ∈ F (K ) and all the edges E ∈ E (F ) is linear. Due to the definition of the trial functions ψzl ∈ Ψhl for z ∈ N (F ) and the nodal interpolation, p and Il p coincide on the polygonal boundary ∂ F of the face F . Since p is a linear polynomial, it fulfills the Laplace equation in the classical sense on F . Furthermore, p ∈ span Φl1 (F ) and consequently fulfills the discrete variational formulation (9) of the Laplace equation in the parameter space of the face F . Since the Galerkin formulation for the Dirichlet problem of the Laplace equation admits a unique solution, the linear function p and its interpolation Il p coincide on the face F . Repeating this argument, we see that p and Il p coincide on the whole boundary ∂ K . With the same considerations, we obtain Il p = p in K , since p and Il p agree on ∂ K and fulfill the Laplace equation in K .  The proof of the continuity of Il in Lemma 2 makes use of an auxiliary discretization of the element into tetrahedra. To this end, we connect all nodes in the triangulation Tl (∂ K ) of ∂ K with the center zK from Definition 2. The tetrahedra constructed this way are denoted by Ttet . Proposition 1. The auxiliary discretization of an element into tetrahedra, as described above, is shape regular in the usual sense of classical finite element methods. So, neighboring tetrahedra either share a common node, edge or triangular face and the aspect ratio of their diameter hTtet and the radius ρTtet of their insphere is bounded uniformly from above, i.e. hTtet /ρTtet < σtet . Here, σtet only depends on the regularity and stability parameters of Definitions 2 and 3 and on the level l in the face discretization. Proof. For an arbitrary tetrahedron, we have the relation

ρTtet =

3VTtet ATtet

,

where VTtet is the volume and ATtet is the surface area of the tetrahedron. This relation is seen as follows. The insphere is perpendicular to the faces of the tetrahedron and thus VTtet is equal to the sum of the volumes VTtet,i , i = 1, . . . , 4, of the four tetrahedrons Ttet,i obtained by connecting the vertices with the center of the insphere. Each volume is computed as VTtet,i = 31 ρTtet |Ti |, where Ti is the triangle on the surface of the initial tetrahedron Ttet . Consequently, it holds that VTtet =

4  i =1

VTtet,i =

4  1 i =1

3

ρTtet |Ti | =

1 3

ρTtet ATtet .

Let K ∈ Kh be an element of a regular and stable mesh. First, we study the case l = 0, where only one node per face is added for the triangulation of the element surface. We consider the auxiliary discretization and choose an arbitrary

S. Rjasanow, S. Weißer / Journal of Computational and Applied Mathematics 263 (2014) 202–217

211

tetrahedron Ttet with corresponding triangle T ∈ Tl (F ) in some face F ∈ F (K ) and with an edge E ∈ E (F ) such that E ⊂ ∂ T ∩ ∂ F . A rough estimate for the surface area of this tetrahedron is ATtet =

4 

|Ti | ≤

i=1

4  h2

K

i=1

2

= 2h2K .

Let dist(zK , T ) be the distance of the center zK to the triangle T and let hTE be the height of the triangle T over the edge E. For the volume of Ttet , we have VTtet =

1 3

dist(zK , T )|T | =

1 6

dist(zK , T )hTE hE .

Since the faces of the element K and the element itself are star-shaped with respect to circles and a ball according to Definitions 1 and 2, it holds that ρF ≤ hTE as well as ρK ≤ dist(zK , T ). Consequently, we obtain VTtet ≥

1 6

ρK ρF hE ≥

1 6σK σF

hK hF hE ≥

1 6σK σF

h3E .

This yields together with the stability in Definition 3 hTtet

ρTtet

=

hTtet ATtet 3VTtet

≤

4σK σF h3K h3E

3 ≤ 4σK σF cK .

In the case l ≥ 1, the volume VTtet gets smaller. The triangle T ⊂ F ∈ F (K ) is obtained by successive splitting of an initial triangle T0 of the mesh with level zero. Due to the construction, these triangles are similar and the relation |T | = |T0 |/4l holds. Taking into account this relation in the considerations above gives the general estimate hTtet

ρTtet

3 ≤ σtet with σtet = 4l+1 σK σF cK . 

Lemma 2. Let Kh be a regular and stable mesh and K ∈ Kh with hK = 1. Then, there exists a constant c which only depends on the regularity and stability parameters as well as on the level l such that

∥Il v∥H 1 (K ) ≤ c ∥v∥H 2 (K ) for v ∈ H 2 (K ). Proof. Let v ∈ H 2 (K ). The interpolation Il v fulfills the weak formulation Find  v ∈ H 1 (K ) : γ0K v = gv

and (∇ v, ∇w)L2 (K ) = 0,

∀w ∈ H01 (K )

with a function gv = Il v|∂ K on the boundary. The data gv is piecewise linear with respect to the surface mesh of level l. To obtain homogeneous boundary data, we decompose  v = v0 +  vg , where  v0 ∈ H01 (K ) and  vg ∈ H 1 (K ) with γ0K vg = gv . According to Proposition 1, the auxiliary discretization of K into tetrahedra is regular in the classical sense of finite elements methods. Therefore, we can use the standard interpolation operator on tetrahedral meshes for linear trial functions to get some  vg . Due to this choice and since hK = 1, it holds that

∥v −  vg ∥H 1 (K ) ≤ C1 |v|H 2 (K ) , see [25], where the constant C1 only depends on the maximal aspect ratio σtet . The reverse triangular inequality yields

∥ vg ∥H 1 (K ) ≤ C1 |v|H 2 (K ) + ∥v∥H 1 (K ) ≤ max{1, C1 } ∥v∥H 2 (K ) . The function  v0 fulfills Find  v0 ∈ H01 (K ) : (∇ v0 · ∇w)L2 (K ) = −(∇ vg · ∇w)L2 (K ) ,

∀w ∈ H01 (K ).

Since K can be embedded into a cube with side length hK , the Poincaré inequality in [32] takes the form

∥w∥L2 (K ) ≤ hK |w|H 1 (K ) for w ∈ H01 (K ). According to this inequality and since hK = 1, it holds that

∥ v0 ∥2H 1 (K ) = ∥ v0 ∥2L2 (K ) + | v0 |2H 1 (K ) ≤ 2| v0 |2H 1 (K ) . Due to the variational formulation for  v0 , we find with the help of the Cauchy–Schwarz inequality that

  | v0 |2H 1 (K ) = (∇ vg , ∇ v0 )L2 (K )  ≤ | vg |H 1 (K ) | v0 |H 1 (K ) which yields

| v0 |H 1 (K ) ≤ ∥ v g ∥H 1 ( K ) .

212

S. Rjasanow, S. Weißer / Journal of Computational and Applied Mathematics 263 (2014) 202–217

The final step in the proof is to combine all estimates as follows

∥Il v∥H 1 (K ) ≤ ∥ v0 ∥H 1 (K ) + ∥ vg ∥H 1 (K ) √ v0 |H 1 (K ) + ∥ vg ∥H 1 (K ) ≤ 2 |  √  ≤ 1 + 2 ∥ vg ∥H 1 (K )  √  ≤ max{1, C1 } 1 + 2 ∥v∥H 2 (K ) .  Applying all previous considerations, we are able to prove the interpolation error estimate (14). Theorem 3. For a regular and stable mesh Kh of a bounded polyhedral domain Ω Il : H 2 (Ω ) → span Ψhl fulfills

⊂ R3 , the interpolation operator

∥v − Il v∥H 1 (Ω ) ≤ ch |v|H 2 (Ω ) for v ∈ H 2 (Ω ) where h = max{hK : K ∈ Kh } and the constant c only depends on the regularity and stability parameters of the mesh and on the level l. Proof. Without loss of generality, we assume h ≤ 1. Let us start to examine the error over one element K ∈ Kh . We scale this element such that its diameter becomes one. For this, the transformation  x → x = zK + hK x ∈ K is applied and the scaled element is denoted by  K . With the notation  v( x) = v(zK + hK x) and the variable transformation in the integrals of the norms, we obtain for v ∈ H 2 (K ) 3/2−k

|v|H k (K ) = hK

| v|H k (K ) for k = 0, 1, 2,

where | · |H 0 (K ) = ∥ · ∥L2 (K ) , and thus  v ∈ H 2 ( K ).

Let  Il be the interpolation operator with respect to  K . Due to the pointwise interpolation, it does not matter if v is first transformed into  v and then interpolated or vice versa, i.e.

 Il v = I l v. Consequently, we obtain

∥v − Il v∥2H 1 (K ) = ∥v − Il v∥2L2 (K ) + |v − Il v|2H 1 (K ) = h3K ∥ v − Il v∥2L (K ) + hK | v − Il v|2H 1 (K ) 2

≤ hK ∥ v − Il v∥2H 1 (K ) since hK ≤ 1. According to the general approximation theory in Sobolev spaces, see [33], there is a linear polynomial  p ∈ P 1 ( K ) which satisfies

| v − p|H k ( v|H 2 (K ) for k = 0, 1. K ) ≤ C (k, σK ) |

(15)

Due to the scaling, the constant C (k, σK ) is independent on the mesh size hK . Applying Lemmata 1 and 2, we obtain

 v − ∥ v − Il v∥H 1 (K ) ≤ ∥ v − p∥H 1 ( p)∥H 1 ( K ) + ∥Il ( K) ≤ (1 + c ) ∥ v − p∥H 2 ( K) ≤ (1 + c )C (σK ) | v|H 2 (K ) ,

(16)

where we also have used (15) and the fact that the second derivatives of  p vanish. Comparing the last two estimates and transforming back to the element K yields

∥v − Il v∥2H 1 (K ) ≤ chK | v|2H 2 (K ) = ch2K |v|2H 2 (K ) . In the last step of the proof, we have to sum up this inequality over all elements of the mesh and apply the square root to it. This gives

1/2

 ∥v − Il v∥H 1 (Ω ) ≤ c

 K ∈Kh

and finishes the proof.



h2K

|v|

2 H 2 (K )

≤ ch |v|H 2 (K )

S. Rjasanow, S. Weißer / Journal of Computational and Applied Mathematics 263 (2014) 202–217

213

Fig. 3. Sequence of Voronoi meshes.

Table 1 Total number of nodes when working with triangulated surfaces of different mesh levels l.

|Kh |

|Nh |

l=0

l=1

l=2

9 76 712 1 316 5 606 26 362

46 416 4 186 7 850 34 427 164 915

98 905 9 081 17 013 74 457 356 189

424 4 170 42 446 79 676 349 663 1 675 171

1 790 18 011 184 170 345 903 1 519 143 7 280 603

5. Numerical experiments The first numerical examples in this section are formulated on the unit cube. As discretization, we utilize Voronoi meshes which are an example of polyhedral meshes. In Fig. 3, the first three meshes of the sequence for the convergence experiments are visualized. We see that the elements are non-trivial polyhedra with arbitrary polygonal faces. The meshes have been produced by generating random points according to [34] and constructing the corresponding Voronoi diagram in accordance with [35]. It is assumed that the mesh generator provides the points zK and zF from the Definitions 1 and 2. However, for convex elements and faces we may use the center of mass instead which is computable. In the setup of the local boundary element matrices, we use a semi analytical integration scheme. The inner integral in the Galerkin matrices is evaluated analytically and the outer one is approximated by Gaussian quadrature. The symmetric systems of linear equations arising on the faces and the global FEM system are solved by the conjugate gradient method without any preconditioning. Of course, for larger problems an efficient solver is of particular interest. It is possible to use a FETI-type solver as was shown in [21]. In Table 1, we sketch the number of elements |Kh | and the number of nodes |Nh | in the different Voronoi meshes. The proposed strategy approximates the solution by a linear combination of as many trial function as nodes are in the mesh. Therefore, the number of degrees of freedom in the finite element method is |Nh | minus the number of nodes on the Dirichlet boundary ΓD . The method proposed in [12] needs to triangulate the surfaces of the elements and the number of trial functions corresponds to the total number of nodes after the triangulation. In Table 1, this total number of nodes is listed in the case that the faces are triangulated with the level l = 0, 1, 2. We recognize that in this situation many more trial functions and thus degrees of freedom are required in the global computations. Roughly speaking, the number of nodes doubles if the coarsest discretization of the faces is used. If a finer triangulation is needed, the number of nodes and thus the number of degrees of freedom increase ten times for l = 1 and even more than forty times for l = 2. Since the diameter of the elements are equal in all four situations, the approximation errors of the finite element computations are of the same order for fixed l. However, the constant in [12] might be better than the one obtained for the presented strategy. This is due to the fact that if h is fixed and only l is increased the method in [12] still converges since it is equivalent to a boundary element domain decomposition approach [13]. The method proposed in this manuscript gives for small l comparable approximations while requiring a minimal set of degrees of freedoms. In the following, we investigate the influence of the face discretization. These triangulations of the faces are required to define the approximated trial functions ψzl ∈ Ψhl on the faces with the help of local, two dimensional finite element methods. The finer the discretization is chosen the better we approximate the original trial functions ψz ∈ Ψh . Even though the face discretization does not blow up the global system matrix, the computational effort for the local problems increases if the discretization level l is raised. As one example, we pick the element K from Fig. 2 and list the number of nodes |Ml (∂ K )| and the number of triangles |Tl (∂ K )| in the surface discretization of K for different levels l in Table 2. The main tasks in the local problems are the evaluation of the boundary element matrix entries and the inversion of the single-layer potential matrix VK ,l which gives a local complexity of O (|Tl (∂ K )|3 ). The next example analyzes the rates of convergence for different values of l.

214

S. Rjasanow, S. Weißer / Journal of Computational and Applied Mathematics 263 (2014) 202–217 Table 2 Number of nodes |Ml (∂ K )| and number of triangles |Tl (∂ K )| in the surface discretization of the element in Fig. 2 for different levels.

|N (K )|

l

|Ml (∂ K )|

|Tl (∂ K )|

12

0 1 2 3 4

20 74 290 1154 4610

36 144 576 2304 9216

Fig. 4. Relative error in ∥ · ∥E (•) and ∥ · ∥L2 (Ω ) (+) with respect to h for l = 0, 1, 2 in Example 1 and triangles with slope one and two.

Example 1. Consider the Dirichlet boundary value problem

−1u = 0 in Ω = (0, 1)3 , u = gD on Γ with gD = γ0 u such that u(x) = e2

√

2π(x1 −0.3)

cos(2π (x2 − 0.3)) sin(2π (x3 − 0.3)),

x ∈ R3

is the exact solution. In Fig. 4, the relative errors ∥u − ulh ∥E /∥u∥E and ∥u − ulh ∥L2 (Ω ) /∥u∥L2 (Ω ) are given with respect to h = max{hK : K ∈ Kh } in a logarithmic plot for different discretization levels l = 0, 1, 2 of the faces. This example shows that the discretization level of the faces does not influence the rates of convergence as we expect from the theory in Section 4. Additionally, Fig. 4 indicates that the constant in the error estimate can be chosen to be independent of the level l. The coarsest face discretization with l = 0 is sufficient to analyze the convergence rates in the forthcoming numerical experiments. Due to this choice, the local complexity in the two dimensional finite element method on the faces F ∈ Fh and the local boundary element methods on the elements K ∈ Kh is rather small. Furthermore, in Fig. 4, we recognize linear convergence for the approximation error measured in the energy norm and quadratic convergence if the error is measured in the L2 -norm. This is the first numerical experiment for the BEM-based finite element method in three space dimensions which confirms these rates of convergence on Voronoi meshes with polyhedral elements and arbitrary polygonal faces. Besides the Dirichlet problem for the Laplace equation, we also give examples for the Poisson problem and the case of a non-constant material coefficient. Example 2. The function u(x) = cos(π x1 ) sin(2π x2 ) sin(3π x3 ), x ∈ R3 fulfills the boundary value problem

−1u = f in Ω = (0, 1)3 , u = gD on Γ with f = 14π 2 u and gD = γ0 u fixed. In Fig. 5, the errors ∥u − ulh ∥E /∥u∥E and ∥u − ulh ∥L2 (Ω ) /∥u∥L2 (Ω ) are shown with respect to h = max{hK : K ∈ Kh } in logarithmic scale. Example 3. We take the two functions already considered in Examples 1 and 2 and label them by u1 and u2 , respectively. They fulfill the boundary value problems

   −div 72 − x1 − x2 − x3 ∇ ui = fi in Ω = (0, 1)3 , ui = giD on Γ

S. Rjasanow, S. Weißer / Journal of Computational and Applied Mathematics 263 (2014) 202–217

215

Fig. 5. Relative error with respect to h for Example 2 with l = 0 and triangles with slope one and two.

Fig. 6. Relative error ∥ui − ulih ∥/∥ui ∥ for i = 1 (•) and i = 2 (+) with respect to h for l = 0 in Example 3 and triangles with slope one and two.

for i = 1, 2, where fi and giD have to be chosen appropriately. In Fig. 6, the relative errors ∥ui − ulih ∥E /∥ui ∥E and ∥ui − ulih ∥L2 (Ω ) /∥ui ∥L2 (Ω ) are shown with respect to h = max{hK : K ∈ Kh } in logarithmic scale for i = 1, 2. In the previous two examples, we have also obtained optimal rates of convergence for the finite element approximation which confirm the theoretical results of the previous section. The BEM-based FEM yields linear convergence in the energy norm and quadratic convergence in the L2 -norm. Finally, we consider the model problem on a L-shaped domain with a singular solution such that u ̸∈ H 2 (Ω ), but 2

u ∈ H 1+ 3 (Ω ). Thus, the preliminaries of Theorem 1 are not fulfilled. However, due to the theory of interpolation spaces we expect a convergence order of 2/3. Example 4. With the help of cylindrical coordinates (r , ϕ, x3 ), where r ≥ 0, ϕ ∈ [π /2, 2π ] and x3 ∈ R, the function u(r cos ϕ, r sin ϕ, x3 ) = r 2/3 sin

2  3

ϕ−

π



2

2

∈ H 1+ 3 (Ω )

fulfills the Laplace equation in

  Ω = (−1, 1) × (−1, 1) × (0, 1) \ [0, 1]3 with appropriate Dirichlet data. In Fig. 7, we give a polygonal mesh of the domain Ω with hanging nodes and edges. Additionally, we show the relative error ∥u − ulh ∥E /∥u∥E with respect to h = max{hK : K ∈ Kh } in logarithmic scale. As expected, we have obtained the reduced order of convergence for a sequence of uniform refined meshes. To recover the linear convergence in the energy norm for singular solutions, it is necessary to perform adaptive strategies, see [18]. 6. Conclusion The classic convergence rates of finite element methods on simplicial meshes are recovered by the BEM-based FEM. However, the new methodology benefits from the flexibility with respect to arbitrary meshes with polyhedral elements having polygonal faces. This behavior together with the conforming approximations makes the BEM-based finite element method an interesting and attractive strategy for ongoing research. Since the trial functions are defined in accordance

216

S. Rjasanow, S. Weißer / Journal of Computational and Applied Mathematics 263 (2014) 202–217

Dirichlet 1.3 1.2 0.80 0.40 0.0

Fig. 7. L-shaped domain with polyhedral mesh and Dirichlet data (left), relative error ∥u − ulh ∥E /∥u∥E with respect to h for l = 0 in Example 4 and triangle with slope 2/3 (right).

with the underlying differential equation, the system of linear equations as well as the approximations contain already some information of the solution. This property is advantageous when considering other differential equations like convection–diffusion, for example, and it has to be investigated further. References [1] I.D. Mishev, Finite volume methods on Voronoi meshes, Numer. Methods Partial Differential Equations 14 (2) (1998) 193–212. [2] D.N. Arnold, F. Brezzi, B. Cockburn, L.D. Marini, Unified analysis of discontinuous Galerkin methods for elliptic problems, SIAM J. Numer. Anal. 39 (5) (2002) 1749–1779. [3] Y. Efendiev, T.Y. Hou, Multiscale Finite Element Methods—Theory and Applications, in: Surveys and Tutorials in the Applied Mathematical Sciences, vol. 4, Springer, New York, 2009. [4] F. Brezzi, K. Lipnikov, M. Shashkov, Convergence of the mimetic finite difference method for diffusion problems on polyhedral meshes, SIAM J. Numer. Anal. 43 (5) (2005) 1872–1896. [5] L. Beirão da Veiga, F. Brezzi, A. Cangiani, G. Manzini, L. Marini, A. Russo, Basic principles of virtual element methods, Math. Models Methods Appl. Sci. 23 (1) (2013) 199–214. [6] R. Hiptmair, A. Moiola, I. Perugia, Plane wave discontinuous Galerkin methods for the 2D Helmholtz equation: analysis of the p-version, SIAM J. Numer. Anal. 49 (1) (2011) 264–284. [7] R. Hiptmair, A. Moiola, I. Perugia, Error analysis of Trefftz-discontinuous Galerkin methods for the time-harmonic Maxwell equations, Math. Comp. 82 (281) (2013) 247–268. [8] O. Cessenat, B. Despres, Application of an ultra weak variational formulation of elliptic PDEs to the two-dimensional Helmholtz problem, SIAM J. Numer. Anal. 35 (1) (1998) 255–299. [9] M. Floater, K. Hormann, G. Kós, A general construction of barycentric coordinates over convex polygons, Adv. Comput. Math. 24 (2006) 311–331. [10] A. Gillette, A. Rand, C. Bajaj, Error estimates for generalized barycentric interpolation, Adv. Comput. Math. 37 (2012) 417–439. [11] N. Sukumar, A. Tabarraei, Conforming polygonal finite elements, Internat. J. Numer. Methods Engrg. 61 (12) (2004) 2045–2066. [12] D. Copeland, U. Langer, D. Pusch, From the boundary element domain decomposition methods to local Trefftz finite element methods on polyhedral meshes, in: Domain Decomposition Methods in Science and Engineering XVIII, in: Lect. Notes Comput. Sci. Eng., vol. 70, Springer, Berlin, Heidelberg, 2009, pp. 315–322. [13] G.C. Hsiao, W.L. Wendland, Domain decomposition in boundary element methods, in: Fourth International Symposium on Domain Decomposition Methods for Partial Differential Equations (Moscow, 1990), SIAM, Philadelphia, PA, 1991, pp. 41–49. [14] C. Hofreither, U. Langer, C. Pechstein, Analysis of a non-standard finite element method based on boundary integral operators, Electron. Trans. Numer. Anal. 37 (2010) 413–436. [15] C. Hofreither, L2 error estimates for a nonstandard finite element method on polyhedral meshes, J. Numer. Math. 19 (1) (2011) 27–39. [16] C. Pechstein, C. Hofreither, A rigorous error analysis of coupled FEM-BEM problems with arbitrary many subdomains, Lect. Notes Appl. Comput. Mech. 66 (2013) 109–132. [17] S. Rjasanow, S. Weißer, Higher order BEM-based FEM on polygonal meshes, SIAM J. Numer. Anal. 50 (5) (2012) 2357–2378. [18] S. Weißer, Residual error estimate for BEM-based FEM on polygonal meshes, Numer. Math. 118 (4) (2011) 765–788. [19] C. Hofreither, U. Langer, C. Pechstein, A non-standard finite element method for convection–diffusion–reaction problems on polyhedral meshes, AIP Conf. Proc. 1404 (1) (2011) 397–404. [20] Y. Efendiev, J. Galvis, R. Lazarov, S. Weißer, Mixed FEM for second order elliptic problems on polygonal meshes with BEM-based spaces, in: Large-Scale Scientific Computing, in: Lecture Notes in Computer Science, Springer, Berlin, Heidelberg, 2014, in press. [21] C. Hofreither, A non-standard finite element method using boundary integral operators, Ph.D. Thesis, Johannes Kepler University, Linz, Austria, December 2012. [22] S. Weißer, Finite element methods with local trefftz trial functions, Ph.D. Thesis, Universität des Saarlandes, Saarbrücken, Germany, September 2012. [23] D.M. Copeland, Boundary-element-based finite element methods for Helmholtz and Maxwell equations on general polyhedral meshes, Int. J. Appl. Math. Comput. Sci. 5 (1) (2009) 60–73. [24] R.A. Adams, Sobolev Spaces, Academic Press, 1975. [25] P.G. Ciarlet, The Finite Element Method for Elliptic Problems, North-Holland, Amsterdam, 1978. [26] O. Steinbach, Numerical Approximation Methods for Elliptic Boundary Value Problems: Finite and Boundary Elements, Springer, New York, 2007. [27] W.C.H. McLean, Strongly Elliptic Systems and Boundary Integral Equations, Cambridge University Press, Cambridge, 2000. [28] S. Rjasanow, O. Steinbach, The Fast Solution of Boundary Integral Equations, in: Mathematical and Analytical Techniques with Applications to Engineering, Springer, 2007. [29] T.Y. Hou, X.-H. Wu, A multiscale finite element method for elliptic problems in composite materials and porous media, J. Comput. Phys. 134 (1) (1997) 169–189. [30] S.E. Mousavi, N. Sukumar, Numerical integration of polynomials and discontinuous functions on irregular convex polygons and polyhedrons, Comput. Mech. 47 (2011) 535–554. [31] V.I. Burenkov, Sobolev Spaces on Domains, Vol. 137, BG Teubner, 1998. [32] D. Braess, Finite Elements: Theory, Fast Solvers, and Applications in Elasticity Theory, third ed., Cambridge University Press, Cambridge, 2007. Translated from the German by Larry L. Schumaker. [33] S. Brenner, L. Scott, The Mathematical Theory of Finite Element Methods, second ed., in: Texts in Applied Mathematics, vol. 15, Springer, New York, 2002.

S. Rjasanow, S. Weißer / Journal of Computational and Applied Mathematics 263 (2014) 202–217

217

[34] M.S. Ebeida, S.A. Mitchell, A. Patney, A. Davidson, J.D. Owens, A simple algorithm for maximal Poisson-disk sampling in high dimensions, Comput. Graph. Forum 31 (2012) 785–794. [35] M.S. Ebeida, S.A. Mitchell, Uniform random Voronoi meshes, in: Proceedings of the 20th International Meshing Roundtable, Springer, Berlin, Heidelberg, 2012, pp. 273–290.