Augmented high order finite volume element method for elliptic PDEs in non-smooth domains: Convergence study

1. Background The Finite Volume ElementMethod (FVEM) is a numerical method for approximating the solution of a partial differential equation (PDE) in a trial function space spanned by piecewise polynomial basis functions, similar to the Finite Element Method (FEM). The coefficients of the linear combination of the basis functions are obtained by imposing the PDE through integrations over control volumes, similar to the finite volume method. Since there are many names and variations for this method, the origin of this method is not entirely clear; however, the paper in [1] is usually cited as one of the first papers on this topic. About a decade after Bank and Rose’s publication, this method was extended to a quadratic method in [2]. Later Li et al. [3] and Xu and Zou [4] have carried out a convergence analysis for the quadraticmethod; see also [5–7]. Error estimates in the L2 norm were given in [8]. Recently, Plexousakis and Zouraris [9] have proven a priori error estimates for the high order FVEM for one-dimensional problems (the ordinary differential equation case) with order higher than two. However, to the best of our knowledge, theoretical convergence results on the FVEM, with arbitrary order piecewise polynomial basis functions for linear elliptic problems in two-dimensional domains, are currently not yet available in the literature. It is well known that a solution of the Poisson equation with an analytic right hand side and analytic Dirichlet boundary data is not necessarily analytic up to the boundary if the boundary is not smooth. Again, it is not easy to pin down the first discovery of such singular behaviour; we refer to the books in [10,11] for detailed discussion. The idea of augmentation of the trial function space for the finite element method (FEM) using results from regularity analysis can be found in [10] (Chapter 8.4.2) aswell as in [12] (Chapter 8). Although it has been shown that the augmentation of the trial function space can recover the optimal convergence rate of the FEM, additional singular basis functions in the trial function space introduce ∗ Corresponding author. E-mail address: [email protected] (Y. Aoki). 0377-0427/$ – see front matter© 2011 Elsevier B.V. All rights reserved. doi:10.1016/j.cam.2011.05.004 Author's personal copy 5178 Y. Aoki, H. De Sterck / Journal of Computational and Applied Mathematics 235 (2011) 5177–5187 complications. Integration of the singular basis function must be done analytically or using special quadrature rules, and special care is needed to preserve the sparsity of the stiffness matrix. Error estimates for the first order FVEM for elliptic PDEs with a derivative blow-up singularity in a non-convex domain are presented in [13]. They show that the rate of error convergence decreases when a singularity is present. Djadel et al. [14] have employed a grid refinement technique, similar to what is presented in Chapter 8.4.1 of [10], for the first order FVEM to improve the rate of convergence. On the other hand, to the best of our knowledge, augmentation of the trial function space has not been studied for the FVEM. In this paper, we demonstrate the augmentation of the trial function space for the arbitrary high order FVEM by presenting numerical convergence studies for the Poisson equationwith derivative blow-up singularity at a reentrant corner. We first perform a numerical convergence study of the FVEM with arbitrary order piecewise polynomial basis functions on triangular grids for non-singular solutions using the systematic way to construct control volumes that was proposed in [15], and that is a generalization of second order approaches in [2–4]; see also [16–18].Well-posedness of the FVEM on triangular grids for this particular way of constructing the control volumes was proved for second order methods in [4] under certain conditions on the angles of the triangular elements. For orders higher than two, it can be observed that H1 convergence orders are the same as the optimal orders exhibited by the Galerkin FEM, but that L2 convergence order is sub-optimal for polynomials of even order, similar to what has been observed before for the Discontinuous Galerkin method [19–21] and for the one-dimensional FVEM [9] (see also [18]). We then show numerically that the presence of a derivative blow-up singularity at a reentrant corner canpollute the numerical solution and the rate of convergence. Finally,we shownumerically how the rate of convergence can be recovered by augmentation of the trial function space, leading to an elegant and efficient augmented high order FVEM. The rest of this paper is organized as follows. Section 2 discusses model problems, and Section 3 presents the numerical method indicating the choice of control volumes and the augmentation approach. Section 4 gives numerical results and Section 5 concludes the paper. 2. Model problems In this paper, we consider the Poisson equation with Dirichlet boundary condition ∆u = f (x, y) in Ω, (1) u = g(x, y) on ∂Ω, (2) where Ω ⊂ R2 is an open polygonal domain with a finite number of vertices, ∂Ω is the boundary of the polygonal domain, f (x, y) is a function in L2(Ω)∩C0(Ω) and g(x, y) is a function inH2(∂Ω). Note that, for simplicity, we only consider classical solutions of the PDE in this paper (f ∈ C0(Ω)). The FVEM can also be used to approximate non-classical solutions (see [1]). 2.1. Regularity of solutions of the Poisson problem Since a polygonal domain is a Lipschitz domain, showing the existence and the interior smoothness of the solution of boundary value problem (BVP), (1)–(2) is straightforward. It can be proven by the Lax–Milgram theorem that there exists a unique solution u ∈ H1(Ω) of BVP (1)–(2) (See Lemma 4.4.3.1 of [10]). Also, it is known that boundary derivative blow-up singularities can occur at the vertices of the domain. We now give a known regularity result for the solution at the vertices. Let (xsi , ysi) ∈ ∂Ω be the vertices with singularities of the boundary of domainΩ , and let Nvert be the number of vertices with singularities. For each vertex (xsi , ysi), local polar coordinates ri and θi are defined as depicted in Fig. 1. Define the interior angle αi for each vertex (xsi , ysi) so that (ri, θi) ∈ Ω if 0 < θi < αi for sufficiently small ri > 0. Using the above notation, the modified shift theorem can be stated as follows. Theorem 2.1 (Modified Shift Theorem). Let u(x, y) be the solution of boundary value problem (1)–(2). Then for all m ∈ N, there exist constants ki,j ∈ R such that u(x, y) − − 0<λi,j<m+1 ki,jψi,j ∈ Hm+1(Ω), (3)