h-p Spectral element methods for three dimensional elliptic problems on non-smooth domains

https://doi.org/10.1016/j.amc.2014.02.025Get rights and content

Abstract

This paper is devoted to the study of non-conforming h-p spectral element methods for three dimensional elliptic problems on non-smooth domains using parallel computers. To overcome the singularities which arise in the neighbourhoods of vertices, edges and vertex-edges we use local systems of coordinates together with geometric meshes which become finer near corners and edges. We provide estimates for the error and the computational complexity of our method and present results of numerical simulations that have been performed to validate the theoretical estimates obtained in [22].

Introduction

It is well known that there are three different types of singularities for elliptic problems on non-smooth domains in R3, namely, the vertex, the edge and the vertex-edge. In addition, the solutions to these problems are anisotropic in the neighbourhoods of edges and vertex-edges (see [1], [2], [3], [4], [18], [19], [20], [21], [23], [26], [38], [39], [40] and references therein for further details).

Babuška and Guo [1], [2] started the study of analytic regularity of elliptic problems on non-smooth domains in R3 in the frame work of weighted Sobolev spaces with Cauchy type control of all derivatives in the so-called countably normed spaces in the neighbourhoods of vertices, edges and vertex-edges in spherical, cylindrical and Cartesian coordinates. Recently, Costabel and coworkers setteled the proof of the analytic regularity estimates [11], [13] using anisotropic weighted Sobolev spaces introduced in [1], [2] by filling the gap which was left over by Babuška and Guo.

A method for obtaining a numerical solution to exponential accuracy for elliptic problems on non-smooth domains in R2 was first proposed by Babuska and Guo [6], [7] in the frame work of the finite element method (FEM). In [16], [17], [36], [37] an exponentially accurate h-p spectral element method (SEM) was proposed for two dimensional elliptic problems on non-smooth domains using parallel computers. The key ingredients were use of auxiliary mappings and geometrical mesh near corners, proposed originally by Babuska and Guo [6], [7].

Among many approaches that have been attempted over the past three decades to provide accurate and economical solutions to elliptic boundary value problems containing singularities, the two principal approaches are adaptive finite element methods (AFEM) and the use of singular basis functions. In three dimensional problems, the AFEM approach, though very successful in many practical applications, faces the problems of increasing computational cost particularly in h-p adaptive refinements (in presence of singularities) [20] while ill-conditioning occur in the latter [20], [26]. In adaptive FEM approach desirable convergence is not guaranteed a priori and crucially depends on algorithmic details (see [10], [31] and literature therein). Recent survey by Nochetto et al. [31] summarizes a comprehensive study on adaptive FEM, their advantages and disadvantages over other standard methods.

As an extension of the second approach, the method of auxiliary mapping (MAM) was introduced by Babuška and Oh (see [5], [27], [28], [29] and references therein). The MAM has proven to be highly successful in dealing with elliptic boundary value problems containing singularities in R2 [5], [27], [28], [29] and R3 [20], [26] using p-version of the finite element method. However, the success of this method depends on optimal choices of auxiliary mappings. The solution of elliptic problems on non-smooth domains has been examined by Pathria and Karniadakis in [30] and Karniadakis and Spencer in [23], in the framework of spectral/hp element methods.

In [22] we have proposed an exponentially accurate h-p spectral element method for solving three dimensional elliptic problems on non-smooth domains by choosing as our solution the spectral element function which minimizes the sum of a weighted squared norm of the residuals in the partial differential equations and the squared norm of the residuals in the boundary conditions in fractional Sobolev spaces and enforce continuity by adding a term which measures the jump in the function and its derivatives at inter-element boundaries in fractional Sobolev norms, to the functional being minimized.

The method we had described is a least-squares collocation method and to obtain our approximate solution we solve the normal equations using Preconditioned Conjugate Gradient Method (PCGM). The residuals in the normal equations can be obtained without computing and storing mass and stiffness matrices [36], [37]. A preconditioner is defined for the method which is a block diagonal matrix, where each diagonal block corresponds to an element. It is also shown in [14], [22] that there exists a new preconditioner which can be diagonalized using separation of variables technique.

In this paper we present results of numerical experiments for a number of test problems on non-smooth domains with analytic as well as singular solutions to validate the theoretical estimates obtained in [22] for a non-conforming h-p spectral element method. In order to show the effectiveness of the method in dealing with elliptic problems on non-smooth domains we shall consider various examples of Laplace and Poisson equations on a number of polyhedral domains containing all three types of singularities namely, vertex, edge and vertex-edge singularities with different types of boundary conditions. This will include all possible cases of singularities that may arise at the corners and edges.

After obtaining our approximate solution at the Gauss–Lobatto-Legendre points consisting of non-conforming spectral element functions we can make a correction to it so that the corrected solution is conforming and is an exponentially accurate approximation to the true solution in the H1 norm over the whole domain [22]. Let uSE be the spectral element solution obtained from the minimization problem and w be the exact solution. Then the relative error is defined as||E||rel=||uSE-w||H1||w||H1.

All our computations are carried out on a distributed memory parallel computer which consists of V20Z and V40Z Sun Fire Servers. Each of the Sun Fire Servers is a 2 CPU 64 bit 2.4 GHZ Oeptron AMD machine with 4 GB RAM. The library used for inter processor communication is Message Passing Interface (MPI). In all our theoretical results the spectral element functions are defined as tensor product of monomials, however, we have used tensor product of Legendre polynomials for simulations. Throughout this paper, Ω,Ω(v),Ω(e),Ω(v-e) etc. will denote polyhedral domains in R3 and (x1,x2,x3),(ϕ,θ,ρ) and (r,θ,x3) will denote the Cartesian, the spherical and the cylindrical coordinates, respectively. We shall also denote the standard Cartesian coordinates by (x,y,z).

This paper is organized as follows. In Section 2, we recall in brief the notations, definitions, structure of various neighbourhoods in the singular regions of the polyhedron, spectral element functions and main results that have been obtained in [22]. Section 3 contains results on error estimates from [22] and description of numerical scheme in brief. Sections 4 Test Problems with smooth solutions, 5 Test problems containing singularities are dedicated to computational results for various test problems on non-smooth domains with analytic and singular solutions containing singularities of different types. Summary of the results and some outlines for future work are provided in Section 6.

Section snippets

Elliptic problems on polyhedral domains

Let Ω denote a polyhedron in R3, as shown in Fig. 1(a). We consider an elliptic boundary value problem posed on Ω with mixed Neumann and Dirichlet boundary conditions:Lw=FinΩ,w=g[0]forxΓ[0]Ω,wNA=g[1]forxΓ[1]=ΩΓ[0].It is assumed that the differential operatorLw(x)=i,j=13-xi(ai,jwxj)+i=13biwxi+cwis a strongly elliptic differential operator which satisfies the Lax–Milgram conditions. Moreover, ai,j=aj,i for all i,j and the coefficients of the differential operator are analytic.

Discretization and local systems of coordinates

In

Error estimates

It is well known that for three dimensional elliptic problems containing singularities at vertices and edges, the geometric mesh and a proper choice of element degree distribution leads to exponential convergence and efficiency of computations (see [3], [18], [21] and references therein). Our error analysis is similar to that in two dimensions [15], [24], [25], [36].

It was shown in Chapter 3 of [22] that the error obtained from the proposed method is exponentially small in N. The optimal rate

Test Problems with smooth solutions

We first demonstrate the performance of our method for various test problems on polyhedral domains on which the solution is smooth. Throughout this section N denote the number of refinements in each direction and W the degree of the polynomials used for approximation. In all our computations in this section we have employed a parallel computer and each element is mapped onto a single processor.

It is known [22] that the error in the regular (smooth) region obeys||uSE-w||H1Ce-bNdof1/3.Here, Ndof

Test problems containing singularities

In order to show the effectiveness of the proposed method in dealing with boundary value problems containing singularities we now present results of numerical simulations for model problems on several polyhedral domains having singularities of different types discussed in Section 2 which aim to verify the validity of the asymptotic error estimates and estimates of computational complexity that have been obtained in [22].

As earlier, N will always denote the number of layers in the geometric mesh

Conclusions

Computational results for a number of test problems on non-smooth domains with analytic and singular solutions confirm the estimates obtained for the error and computational complexity. In the edge neighbourhoods the behavior of the solution is very different in the direction of the edge and perpendicular to the edge. Thus, to obtain the exponential accuracy we have to refine the mesh geometrically perpendicular to the edge and optimal choice of polynomial degree distribution is needed.

References (39)

  • I. Babuška et al.

    The p-version of the finite element method for domains with corners and for infinite domains

    Numer. Methods PDEs

    (1990)
  • I. Babuška et al.

    Regularity of the solutions of elliptic problems with piecewise analytic data, part I: boundary value problems for linear elliptic equation of second order

    SIAM J. Math. Anal.

    (1988)
  • I. Babuška et al.

    The h-p version of the finite element method on domains with curved boundaries

    SIAM J. Numer. Anal.

    (1988)
  • I. Babuška et al.

    The h-p Version of the finite element method, part II: general results and applications

    Comput. Mech.

    (1986)
  • C. Bernardi et al.

    The mortar element method applied to spectral discretizations

  • M. Costabel, M. Dauge, S. Nicaise, Corner singularities and analytic regularity for linear elliptic systems, part I:...
  • P.K. Dutt et al.

    Non-conforming h-p spectral element methods for elliptic problems

    Proc. Indian Acad. Sci. (Math. Sci.)

    (2007)
  • P. Dutt et al.

    Stability estimates for h-p spectral element methods for general elliptic problems on curvilinear domains

    Proc. Indian Acad. Sci. (Math. Sci.)

    (2003)
  • P. Dutt et al.

    Stability estimates for h-p spectral element methods for elliptic problems

    Proc. Indian Acad. Sci. (Math Sci.)

    (2002)
  • Cited by (9)

    • An efficient numerical technique for estimating eigenvalues of second-order non-self-adjoint Sturm–Liouville problems

      2022, Mathematics and Computers in Simulation
      Citation Excerpt :

      Unluckily, there are not many numerical techniques available to compute eigenvalues of this type of problems efficiently. It is well known that the spectral techniques, collocation, Galerkin and Tau, are one of the most popular among existing numerical methods for solving the self-adjoint problems (see, for example [1,2,11,12,19,21,22,24,28,30,33,37,41]). The principal advantage of these techniques lies in their higher accuracy for a given number of unknowns [8,9,29].

    • A hk mortar spectral element method for the p-Laplacian equation

      2018, Computers and Mathematics with Applications
    • Spectral element method for three dimensional elliptic problems with smooth interfaces

      2017, Computer Methods in Applied Mechanics and Engineering
    • H-p spectral element methods for three dimensional elliptic problems on non-smooth domains, Part-III: Error estimates, preconditioners, computational techniques and numerical results

      2016, Computers and Mathematics with Applications
      Citation Excerpt :

      The second paper [16] addresses proof of the stability theorem. Results of numerical experiments that have been performed to validate the theoretical estimates are presented in [17]. In [15,16] we had further divided each of the elements in the regular region, vertex neighbourhoods, edge neighbourhoods and vertex-edge neighbourhoods into still smaller elements as curvilinear hexahedrons, tetrahedrons and prisms using a geometric mesh (Fig. 2) and by virtue of the fact that a tetrahedron can be split into four hexahedrons [21,22] and a prism can be split into three hexahedrons we can assume that all our elements are hexahedrons to keep the presentation simple.

    • Least-squares spectral element method for three dimensional Stokes equations

      2016, Applied Numerical Mathematics
      Citation Excerpt :

      In [10], Cockburn and Gopalkrishnan introduced a method which gives exactly incompressible velocity approximations to Stokes flow in three dimensions. In this paper we propose an exponentially accurate spectral element method for solving three dimensional Stokes equations by choosing (as our solution) the spectral element functions which minimize the sum of a squared norm of the residuals in momentum equations, continuity equation and the squared norm of the residuals in the boundary conditions in fractional Sobolev spaces and enforce continuity by adding terms which measure the jumps in the velocity and pressure components along inter-element boundaries in appropriate fractional Sobolev norms, to the functional being minimized as in [14]. This paper is organized as follows.

    View all citing articles on Scopus

    This work was carried out at the Department of Mathematics, Indian Institute of Technology Kanpur, Kanpur 208016, India.

    View full text