Elliptic Problems in Curved Domains Using Cubature Points Based Triangular Spectral Elements and Isoparametric Mappings

Richard Pasquetti

Lab.J.A.Dieudonn′e,Universit′e C?te d’Azur,CNRS,Inria,F-06108 Nice,France.

Abstract.Using the cubature points based triangular spectral element method and isoparametric mappings, we provide accuracy results for elliptic problems in non polygonal domains.Two regimes of convergence,associated to the bulk and to the boundary of the computational domain are clearly discerned and an efficient way to define the isoparametric mapping is proposed.

Key words:Spectral element method,triangular elements,curved domains,isoparametric mappings,cubature points,diagonal mass matrix.

1 Introduction

The present paper extends the study presented in[14],concerned with high order Finite Element Methods(FEMs)in non polygonal domains, in order(i)to provide a deeper analysis of the accuracy results and also(ii)to propose a way to improve them.More precisely,in[14]it was pointed out that when isoparametric mappings are used to approximate a curved boundary,then the convergence of the approximate solution to the exact one,with respect to the polynomial degree,shows two regimes,the first one being associated to the bulk of the domain and the other being controlled by the boundary part.Here we want to clarify this point.Moreover,even if focusing on a particular type of isoparametric mapping,a curved boundary can be approximated in many different ways.This paper provides a way to define the mapping that yields better accuracy results than the way used in[14].

As in[14],the study is restricted to isoparametric mappings,i.e.,we do not consider different treatments of curved domains like those involved with trans finite mappings,with the Isogeometric FEM or with the Nurbs Enhanced FEM,see[15]and references herein.Despite the fact that isoparametric mappings show some failures,especially in general an approximate description of the boundary,they indeed remain simple to implement and moreover preserve some important properties of Spectral Element Methods(SEMs),like the diagonal feature of the mass matrix.Moreover,from results first presented in[12],where different isoparametric mappings are compared,we restrict ourselves to isoparametric mappings based on the resolution of Partial Differential Equations(PDEs),namely the Laplace equation(harmonic extension)or the equation of linear elasticity.

With respect to some other high order FEMs developed for triangular elements,e.g.,the Fekete-Gauss triangular spectral element method(TSEM)[11],the Fekete(resp.Gauss)points of the triangle being used as interpolation(resp.quadrature)points,as in[14]here we use the “CubatureTSEM”,i.e.,only based on the so-called cubature points,see[2,4,7,10].Because with the CubatureTSEM the interpolation and quadrature points coincide,the mass matrix is diagonal,so that this interesting property of the usual SEM is recovered.Indeed,in quadrangles the usual SEM makes use of the tensorial product of Gauss-Lobatto-Legendre(GLL)points both for interpolations and quadratures.This may be useful,e.g.,to address evolution problems with explicit time schemes or to easily define high order differentiation operators[9].Indeed,in both cases the required inversion of the mass matrix can be trivially achieved.

High order FEMs are generally based on a polynomial approximation of total degreeNin each triangular element,the number of degrees of freedom being then equal to(N+1)(N+2)/2.However,it is known that in the frame of Lagrangian methods one cannot find a single set ofn=(N+1)(N+2)/2 points,with 3Npoints on the triangle boundary and(N-1)(N-2)/2 inside,such that an enough powerful quadrature rule exists[5,16].To overcome this difficulty,as first suggested in[2]one should enrich the spaceof the polynomials of total degreeNdefined on the reference triangle(master element),with bubble functions of degreeN′＞N,i.e.,to make use of the space PNUbPN′-3wherebis the unique bubble function of P3().This allows to increase the number of inner points up to(N′-1)(N′-2)/2 while keeping equal to 3Nthe number of boundary points.Up to now,satisfactory results have been obtained withN′=N+1 for 1＜N＜5,N′=N+2 forN=5 andN′=N+3 for 5＜N＜9[7].Such values ofN′are those obtained when requiring that the integrals of polynomials of degreeN+N′-2 are exactly computed,from a theoretical result that traces back to the 70’s[1],but they could be decreased from the less demanding criterion recently proposed in[3].

Note that when using the so-called “condensation technique”,one first computes the solution at the nodes located at the sides of the triangular elements and then locally its inner point values.As a result,the size of the algebraic system increases linearly withN,rather than inN2.Moreover,the system only shows aO(N)condition number and usingN′＞Nyields a negligible increase of the computational cost.Although using the CubatureTSEM,with implementation of the condensation technique,the results of the present study are not specific to this approach.

Finally,it should be mentioned that although the present study focuses on the two dimensional(2D)case,similar behaviors should be expected in the 3D one.

2 Numerical evidence of the two regimes of spectral convergence

We consider the elliptic equation-Δu+u=fin the star domain introduced in[14].In thexyplane,this domain,say Ω,has its boundary Γ parameterized by the polar angleθ,such that:

A P1mesh of this star domain is shown in Figure 1.One notices that domain Ω is especially interesting as curved domain to carry out an accuracy study,since it shows both concave and convex parts and moreover cannot be exactly approximated by polynomials.In Ω,we consider exact solutions of the formuex=cos(ax)cos(ay),whereais a parameter that allows to control the oscillatory behavior of the solution.The source termfas well of the boundary conditions,say of Dirichletu|Γ=gor of Robin type(?nu+u)|Γ=g,are defined accordingly.Note that the implementation of Robin as well as non homogeneous Neumann conditions is not straightforward when using the CubatureTSEM.Indeed,at the edges of the elements the cubature points do not coincide with Gauss points,so that a different set of points,e.g.,the GLL points,should be used to compute the boundary integrals,see[13,14]for details.

Figure 1:P1-mesh of the star shaped computational domain.

Figure 2:For a∈{5,10,20},max norm of the error versus the polynomial degree,when using the harmonic extension isoparametric mapping and for the polygonal domain,with Dirichlet(left)and Robin conditions(right).

Figure 2 shows the max norm of the error versus the polynomial approximation degree for three different values of the coefficienta,i.e.a∈{5,10,20}.In Figure 2(left)the results obtained(i)with the harmonic extension and(ii)if assuming polygonal the computational domain are compared.The polygonal domain,say Ω1of boundary Γ1,is simply defined by the P1mesh,so that no isoparametric mapping is involved.Critical values of the polynomial degree are clearly pointed out.Foragiven,below these values the convergence curves coincide,meaning that the error is associated to the bulk of the domain,whereas beyond them the influence of the isoparametric mapping is clearly pointed out.One observes that the decrease of the error remains exponential but with a smaller convergence rate.Moreover,more oscillatory is the solution greater is the critical polynomial degree,sayN=3(a=5),N=5(a=10)andN≥8(a=20).Indeed,a smooth exact solution can be easily captured,so that the criticalNis low.Beyond it,the convergence rate is controlled by the boundary fitting.

The influence of the boundary condition is shown in Figure 2(right),where a Robin condition is used.The analysis carried out for the Dirichlet case keeps true but the errors are greater,since additional sources of error are present,e.g.,the approximation of the normal to the boundary.

Figure 3:For a∈{5,10,20},max norm of the error versus the polynomial degree,when using the harmonic extension and the linear elasticity isoparametric mappings with Dirichlet(left)and Robin conditions(right).

The present results are not specific to the harmonic extension mapping.Thus,Figure 3 points out that the results obtained with the harmonic extension and with the linear elasticity mappings compare well,both in the Dirichlet and Robin cases.Recall however that,differently to the linear elasticity mapping,in some extreme cases the harmonic extension may fail to define an isoparametric mapping[8].

3 A better choice to define the isoparametric mapping

Let us now investigate a way to improve the accuracy results.For a given type of isoparametric mapping,e.g.,the harmonic extension,to fully define the mapping one has first to define the boundary nodes of the curved edges.An implicit definition of the isoparametric mapping is thus obtained.In[12,14],for a given curved triangular element the boundary nodes on the curved edge are chosen at the intersections of the lines that join the opposite vertex and the nodes obtained from the linear P1mesh,see Figure 4(left).This is relevant when transfinite mappings are involved,especially the bending procedure that we introduced in[6],but is less justified when PDEs based isoparametric mappings are chosen.Indeed,the approximation of the curved boundary at one edge of the element has not to depend on the localization of the vertex opposite to this edge.Following what is generally done for quadratic element,one may define the nodes of the curved boundary by stating that their orthogonal projections onto the sustaining straight edge coincide with the nodes that result from the linear P1-mappping,see Figure 4(right).

Figure 4:Definition of the boundary nodes by the opposite vertex(at left)and by orthogonal projection(at right).

Using the harmonic extension,let us revisit the elliptic problems previously introduced,again witha={5,10,20}in the exact solution.The improvement obtained when defining the boundary nodes by orthogonal projection is shown in Figure 5,for the Dirichlet(at left)and Robin(at right)problems.For the most oscillating exact solution,the curves coincide in the Dirichlet case and only slightly differ for the Robin one.This is however untrue for the less oscillating solutions,i.e.,when for largeNthe error is controlled by the boundary region.

Figure 5:For a∈{5,10,20}and when using the harmonic extension,max norm of the error with the boundary nodes defined by the opposite vertex(label V)and by projection(label P),for the Dirichlet problems(at left)and for the Robin problems(at right).

For a given elementTwith one edge on the boundary,let us denote byandthe images of,in the reference element,obtained with the PNand P1mappings,respectively.Then,withone has:

Introduce now the outwards unit vectororthogonal to the edge to be curved and the outwards unit vectorsparallel to the.Then,with

Defining the boundary nodes by projection means choosingfor anyj.Then,in the local coordinates system associated to the edge,the deformation amplitudedis just a polynomial of degreeN.Thus,forN=2 one obtains a symmetric arc of parabola.

4 Detailed analysis of the accuracy results

Let us introduce the computational domain Ω2obtained by substituting to the edges at the boundary symmetric arcs of parabola.As explained previously,see Section 3,if we define the boundary nodes by projection,then the boundary Γ2of Ω2is perfectly described.Using Ω2as computational domain should indicate if the loss of accuracy,with respect to polygonal domains,comes from the boundary of the domain.

Figure 6:Top panel:Cubature nodes for the the domain Ω(label PN),as obtained with the harmonic extension(the boundary nodes being defined by projection).Comparison to those associated the polygonal domain Ω1(label P1).Central panel:Cubature nodes for the domains Ω (label PN)and Ω2(label P2-PN),as obtained with the PNmapping associated to the harmonic extension.For the domain Ω2,the curved boundary is just a symmetric arc of parabola.Bottom panel:Cubature nodes obtained for the domain Ω2with the P2-subparametric mapping(label P2)and with the PNmapping(label P2-PN).Inside the element the nodes differ but they coincide at the element boundary.

Moreover,for this very specific domain one can also make use of the P2-mapping usually used for quadratic finite elements.Thus,if we assume that the elementA1A2A3has its edgeA2A3on the boundary,then the P2-mapping is simply defined bywherehis the distance of the parabola summit to the straight edge and(λ1,λ2,λ3)the usual barycentric coordinates.If using the former strategy and,e.g.,the harmonic extension,the Jacobian determinant is a polynomial of degree 2N-2,whereas for the latter one it is of degree 2,independently ofN.This should indicate if the source of inaccuracy comes from a non satisfactory quadrature in the curved triangles.For a specific element,we show in Figure 6 the various distributions of nodes that can be obtained.

Using the exact solution and parametrization introduced in Section 2,the convergence curves obtained with the harmonic extension and the P2-mapping are given in Figure 7,for the Dirichlet(at left)and Robin(at right)problems.The results obtained for the domain Ω and for the polygonal domain Ω1are also presented.Clearly,for the Dirichlet problem using the P2-subparametric mapping allows to recover an optimal accuracy,at least tillN=7.In the Dirichlet case the loss of accuracy with respect to the polygonal domain Ω1is thus linked to the surface integrals inside the curved elements.For the Robin problem,that additionally involves boundary integrals,the result is different since one clearly observes an influence of the boundary region both for the P2and PNmappings.Moreover,the curves obtained with the two mappings nearly coincide.Indeed,for the piece-wise parabolic domain Ω2the numerical quadratures involved with the P2and PNmappings are exactly the same.In the Robin case,the loss of accuracy is thus governed by the evaluations of the integrals at the edges.Moreover,we have checked that this also holds(i)in the Neumann case and(ii)for the Poisson equation.

One may then think to enhance the accuracy of the computations of the boundary integrals.It can be done simply by increasing the number quadrature points,since,as explained in[13],to compute the boundary integrals we make use of Gauss points rather than the cubature points located on the edges.We have checked this approach,using 2N+1 rather thanN+1 GLL points,without finding any improvement.Indeed,the failure comes from the polynomial approximation of the kernel of the integral,but this approximation remains associated to the cubature points.

To summarize,in the Dirichlet case it is required(i)to exactly describe the geometry(like e.g.done with the isogeometric FEM)and additionally(ii)to make use of a subparametric mapping to recover an optimal accuracy,i.e.,to recover the exponential convergence rate obtained with the polygonal domain Ω1,but this remains insufficient in the Robin case.This lack of convergence is due to the point-wise evaluation of the side Jacobian which is not polynomial,except in the particular case where it is constant(polygonal domain),so that even for the piece-wise parabolic domain Ω2the convergence rate remains suboptimal.Thus,withN＞1 one can compute the surface of the domain Ω2up to the machine accuracy,but not its perimeter,see Tab.1.

Finally,for both the Dirichlet and Robin problems one observes that the curves labeledPN and P2-PN(obtained for the domains Ω and Ω2,respectively)are close to be parallel,meaning that when describing exactly the geometry the accuracy is improved but not the convergence rate.

Table 1:Relative error in the computation of the surface of Ω2and of its perimeter when using the harmonic extension.Similar results are obtained with the P2-mapping.

Figure 7:Computational domain Ω2,a=5 in the exact solution;Max norm of the error when using the harmonic extension(label P2-PN)and with the P2-subparametric mapping(label P2)for the Dirichlet problem(at left)and for the Robin problem(at right).The results obtained for the domain Ω(label PN)and for the polygonal domain Ω1(label polygon)are also presented.

Acknowledgments

The P1FEM mesh has been generated with the free software “Triangle”.We are grateful to Dr.Youshan Liu,Chinese Academy of Sciences,Beijing,for transmitting us the cubature points and weights used in[7].