A discontinuous finite element formulation for Helmholtz equation

Comput. Methods Appl. Mech. Engrg. 195 (2006) 4018–4035 www.elsevier.com/locate/cma

A discontinuous finite element formulation for Helmholtz equation Gustavo Benitez Alvarez a,*, Abimael Fernando Dourado Loula a, Eduardo Gomes Dutra do Carmo b, Fernando Alves Rochinha b a

LNCC—National Laboratory of Scientific Computation, Av. Getlio Vargas 333, 25651-070, P.B. 95113, Petrpolis, RJ, Brazil b COPPE—Federal University of Rio de Janeiro, Ilha do Funda˜o, 21945-970, P.B. 68509, Rio de Janeiro, RJ, Brazil Received 17 August 2004; accepted 26 July 2005

Abstract A discontinuous finite element formulation for Helmholtz problems is presented with C0 continuity across inter-element boundaries enforced in a weak sense depending on two parameters b and k whose choice is crucial for the performance of the proposed method. These parameters are chosen through one dimension numerical experiments aiming at minimizing the pollution error. Optimal values of these parameters are then adopted in more general situations. The accuracy and stability of the proposed formulation for linear and bilinear shape functions are demonstrated in several numerical examples in one and two dimensions.  2005 Elsevier B.V. All rights reserved. Keywords: Discontinuous Galerkin; Helmholtz equation; Stabilization; Discontinuous FEM

1. Introduction The Helmholtz equation is a linear mathematical model that describes time-harmonic acoustic, elastic and electromagnetic steady state waves. One major problem in approximating this equation by classical finite element methods is linked to the loss of ellipticity with an increasing excitation frequency. Classical Galerkin finite element approximation provides good phase and amplitude accuracy as long as the mesh is fine enough with respect to the wave number in the propagation region. However, ‘‘fine enough’’ is often too expensive for adequate resolution, even for moderate wave numbers. Hence, one of the main concerns in acoustic finite element analysis is the adequacy of the finite element mesh. Acousticians often use the, so-called, ‘‘rule of the thumb’’ [29] which prescribes a relation between the minimal number of elements and the wave number. This rule takes into account only the interpolation properties of polynomial basis in a specific mesh but not the stability of the formulation. Consequently, it works well for low frequencies but for high frequencies very poor results are obtained even when the finite element mesh obeys this rule. A numerical analysis of this problem is presented for linear elements in [29] where the estimate eh 6 C1kh + C2k3h2, kh < 1, for the relative error of the pressure gradient is derived with the constants C1 and C2 independent of wave number k and element size h. The first term in the right side of the above estimate represents the interpolation error and the second one represents the pollution effect. To keep the interpolation error constant it is enough to maintain kh constant. That is the ‘‘rule of the thumb’’ which corresponds to taking a certain number of elements per wavelength. As one can clearly see, this is not sufficient to maintain the pollution error under control. Thus, even following the heuristic rule, the pollution error increases with k. *

Corresponding author. E-mail addresses: [email protected] (G.B. Alvarez), [email protected] (A.F.D. Loula), [email protected] (E.G.D. do Carmo), rochinha@adc. coppe.ufrj.br (F.A. Rochinha). 0045-7825/$ - see front matter  2005 Elsevier B.V. All rights reserved. doi:10.1016/j.cma.2005.07.013

G.B. Alvarez et al. / Comput. Methods Appl. Mech. Engrg. 195 (2006) 4018–4035

4019

The above scenario constitutes a real challenge for numerical analysts and has driven several different approaches for alleviating the pollution effect. This is exposed in a vast literature [1–27,29] and the relative merits of some of these formulations are presented in [26,27,29]. In the present article, a discontinuous finite element formulation is proposed aiming at improving the accuracy of the sound wave prediction. Refs. [7,19,23] also introduce discontinuous methods for the acoustic problem which fundamentally differ from the one proposed in the present article by the use of non polynomial basis functions. They use wave like functions in order to better capture the unique features corresponding to medium and high frequency regime. Similar goals are pursued here but the use of standard finite element basis is considered as starting point to keep the simplicity of coding the formulation. The use of discontinuous finite element methods for solving second-order differential equations has its origins in the early 1970s [30–34]. Later on, their use were successfully extended to a broader classes of problems [35–44]. In particular, the discontinuous finite element formulation developed in [44] is applied here to the Helmholtz equation. In the proposed method, the continuity across the element edges is relaxed and weakly enforced by means of interface terms collecting possible jumps of the sought field and its gradient. Those terms contain two scalar parameters that can be tuned to improve the performance of the method. Indeed, the discontinuous formulation introduced in [44] consists of a family of methods parameterized by those parameters referred to along the text as b and k. Particular choices of these parameters allow to recover other methods within the literature, like, for instance, those in [37–39] correspond to k = 1 and b = 0, which reveals to be unstable for linear polynomials. The symmetrical methods of [30,32,33] are retrieved by adopting k = 1 and b > b0 > 0. Therefore, the formulation of [44] is quite generic which constitutes a factor of attractiveness for different problems as judicious choices of k and b leads to robust methods that takes into account the special features of the target situation. This fact was particularly important in its application to the present acoustic problem, as demonstrated by the numerical simulations of Section 4. These experiments clearly show better convergence behavior of the proposed method compared to other finite element formulation. Moreover, this convergence behavior is close to the interpolation error, indicating that the accuracy of the discontinuous method is basically dependent on the pragmatic ‘‘rule of the thumb’’. The remainder of the paper is organized as follows. In Section 2, the model problem along with the standard C0 Galerkin approximation are presented. In Section 3, the discontinuous finite element formulation proposed aiming at solving the Helmholtz equation is introduced. In Section 4, the results of several numerical experiments are shown, as a way to evaluate the performance of discontinuous finite element formulation. Finally, Section 5 contains some concluding remarks. 2. The Helmholtz equation Let X  Rn (n = 1, 2) be an open bounded domain with Lipschitz continuous boundary C, what guarantees the existence of an unit normal n defined almost everywhere on C. Consider the following boundary value problem:  r  ru  k 2 u ¼ f

ð1Þ

in X;

u ¼ g on Cg ; ru  n ¼ q on Cq ;

ð2Þ ð3Þ

ru  n þ au ¼ r

ð4Þ

on Cr ;

with Cg [ Cq [ Cr = C and Cg \ Cq = Cg \ Cr = Cq \ Cr = ;. The above problem is the well known Helmholtz equation that describes the time-harmonic sound response in a non viscous medium. The field u stands for the sound pressure, k is the wave number, f is the volume source term and g, q and r are the boundary conditions prescribed for the problem. 2.1. Continuous Galerkin finite element method Introducing   U ¼ u 2 H 1 ðXÞ : u ¼ g on Cg ;   V ¼ v 2 H 1 ðXÞ : v ¼ 0 on Cg in which H1(X) is the standard Hilbert space, the boundary value problem (1)–(4) is rephrased in a weak form as follows: find u 2 U such that aðu; vÞ ¼ f ðvÞ with aðu; vÞ ¼

Z X

8v 2 V

ð5Þ

2

½ru  rv  k uv dX þ

Z auv dC; Cr

ð6Þ

4020

G.B. Alvarez et al. / Comput. Methods Appl. Mech. Engrg. 195 (2006) 4018–4035

f ðvÞ ¼

Z

qv dC þ

Z

Cq

rv dC þ Cr

Z

ð7Þ

fv dX. X

Let us consider Mh = {X1, . . . , Xne} a finite element partition of X consisting of standard C0 elements (triangular or quadrilateral) with boundary Ce, such that: X¼

ne [

Xe ;

X ¼ X [ C;

Xe ¼ Xe [ Ce

and Xe \ Xe0 ¼ ;

if e 6¼ e0

e¼1

and the finite element subset Uh,l  U and the subspace Vh,l  V: U h;l ¼ fuh 2 H 1 ðXÞ : uhe 2 P l ðXe Þ

and uh ¼ g

V h;l ¼ fvh 2 H 1 ðXÞ : vhe 2 P l ðXe Þ and vh ¼ 0

on Cg g; on Cg g;

where h is the mesh size parameter, Pl(Xe) is the space of polynomials of degree l, gh is the interpolation of g and uhe denotes the restriction of uh to each element Xe. The Galerkin finite element approximation to (5) consists of finding uh 2 Uh,l that satisfies aðuh ; vh Þ ¼ f ðvh Þ 8vh 2 V h;l .

ð8Þ

Although convergent, this Galerkin finite element method presents a well known pollution effect on its rates of convergence. This pollution effect is very well documented in the literature and is mainly characterized by a significant phase lag spread all over the domain (dispersion error). Numerous techniques have been proposed to overcome this misbehavior ranging from simple rules to very sophisticate formulations. Given the oscillatory behavior of the solution a ‘‘rule of the thumb’’, frequently adopted in engineering practice, states that a minimal number of elements is required per wave length such that kh < 1, which is equivalent to approximately 6 elements per wave length. For large wave numbers very poor results are obtained even when this rule is obeyed, as it is clearly demonstrated by the a priori error estimates presented in [29]. To reduced the pollution effects we introduce the following finite element formulation based on discontinuous interpolation across the element edges. 3. A discontinuous finite element method The present finite element method is grounded on a weak version of the boundary value (1)–(4) in which the continuity across the element edges is imposed in a weak sense. To this end, we define new function spaces: H ðX; M h Þ ¼ fw 2 L2 ðXÞ : we 2 H 1 ðXe Þ and   U DG ¼ w 2 H ðX; M h Þ : w ¼ g on Cg ;   V DG ¼ v 2 H ðX; M h Þ : v ¼ 0 on Cg ;

r  rwe 2 L2 ðXe Þg;

and redefine our boundary value problem as: find u 2 H(X, Mh) that satisfies: r  rue  k 2e ue ¼ fe

ð9Þ

in Xe ;

with boundary conditions ue ¼ ge on Cg \ Ce ; rue  ne ¼ qe on Cq \ Ce ;

ð10Þ ð11Þ

rue  ne þ ae ue ¼ re

ð12Þ

on Cr \ Ce ;

and interface conditions ue  ue0 ¼ 0;

ðrue  rue0 Þ  ne ¼ 0 a.e. on Ce \ Ce0 ;

ð13Þ

with ue the restriction of the function u to the subdomain Xe. The family of discontinuous methods introduced here relies on a variational formulation to problem (9)–(13) including jump terms across the element edges, as follows: find u 2 UDG that satisfies: ak ðu; vÞ ¼ aG ðu; vÞ þ aDG ðu; vÞ ¼ fG ðvÞ

8v 2 V DG ;

ð14Þ

G.B. Alvarez et al. / Comput. Methods Appl. Mech. Engrg. 195 (2006) 4018–4035

4021

with aG ðu; vÞ ¼

ne Z X

fG ðvÞ ¼

"Z

aDG ðu; vÞ ¼

Cr \Ce

fe ve dX þ

Z

qe ve dC þ

"

ne X ne Z X e¼1 e0 >e

Cee0

#

Z

Cq \Ce

Xe

e¼1

 ae ue ve dC ;

Z

Xe

e¼1 ne X

½rue  rve  k 2 ue ve  dX þ

re ve dC ; Cr \Ce

1  ðrue þ rue0 Þ  ne ðve  ve0 Þ: 2

# 0 b kee ðue  ue0 Þðrve þ rve0 Þ  ne dC; ðue  ue0 Þðve  ve0 Þ þ þ hee0 2 ee0

where hee0 ¼ minfhe ; he0 g; ee0

ð15Þ

ee0

b and k are functions to be determined aiming at reducing the pollution effects compared to standard Galerkin formulations. This class of weak formulation is introduced in [44] where a detailed mathematical analysis including equivalence between the corresponding differential and variational formulations, existence and uniqueness of solution is presented for strongly elliptic problems. The equivalence between the strong form (9)–(13) and the weak form (14) is established following the same steps used in [44]. Indeed, the well posedeness of the proposed problem represented by the existence and uniqueness of solution to (14) is obtained considering the equivalence between (9)–(14) and avoiding resonance by considereing values of k for which the only solution of the homogeneous problem: find u 2 VDG such that ak ðu; vÞ ¼ 0 8v 2 V DG ; is the trivial one u = 0. Avoiding resonance we prove the following Ga¨rding condition: there exist c1 > 0 and c2 > 0 such that 2

ak ðv; vÞ þ c1 ðv; vÞ P c2 kvkV DG

8v 2 V DG \ H 1 ðXÞ;

which expresses the stability of the proposed variational formulation. This Ga¨rding condition is a consequence of the choice c1 > k2 + c2 and the VDG \ H1(X)-ellipticity of the bilinear form ak(Æ, Æ) for k = 0, i.e., there exists c2 > 0, independent of he, such that 2

ak ðv; vÞ P c2 kvkV DG

8v 2 V DG \ H 1 ðXÞ and k ¼ 0;

with 2 kvkV DG

¼

ne Z X

2

ðjrve j þ

v2e Þ dX

Xe

e¼1

þ

ne X ne Z X e¼1 e0 >e

Cee0

1 2 ðve  ve0 Þ dC. hee0

Introducing the discontinuous finite dimension spaces of degree l P 1, 2 l U h;l DG ¼ fu 2 L ðXÞ : ue 2 P ðXe Þ

V

h;l DG

2

and u ¼ gh

on Cg g;

l

¼ fv 2 L ðXÞ : ve 2 P ðXe Þ and v ¼ 0 on Cg g;

the corresponding finite element approximation of Helmholtz problem yields: find uh 2 U h;l DG that satisfies aG ðuh ; vh Þ þ aDG ðuh ; vh Þ ¼ fG ðvh Þ 8vh 2 V h;l DG ;

ð16Þ

with aG(u, v), aDG(u, v) and fG(v) as defined before. 0 0 The above formulation consists of a family of methods parameterized by the pair bee ; kee . For strongly elliptic problems, optimal values of those parameters can be determined by a priori error estimates. As this is not the case for Helm0 0 holtz equation in the high-frequency regime, optimal values for the pair bee ; kee will be determined through numerical experiments in a one dimensional problems as presented in the next section. All choices of k but one (k = 1) lead to non-symmetric formulations what, obviously, has an impact on the computational costs which is compensated by the stability and accuracy of the discontinuous formulation as presented in the next section.

4022

G.B. Alvarez et al. / Comput. Methods Appl. Mech. Engrg. 195 (2006) 4018–4035

4. Numerical results A number of examples to illustrate the main features and potential of discontinuous finite element methods to solve 0 0 Helmholtz problem is presented. For simplicity, from now on the parameters bee ; kee will have fixed values b and k for all interelement edges since we will consider only uniform meshes. We remind that depending on the choice of parameters b and k, the global matrix of our problem is not symmetric. Moreover, although sparse, the global matrices of the discontinuous formulation have larger band width than those corresponding to C0 Galerkin methods with the same mesh and polynomial interpolation. These two aspects directly impact the computational costs of the proposed methodology. However, as we will show later, the impact on computational cost is largely compensated by the gain in accuracy of the proposed method (see Table 3). At that point, before presenting the numerical results, it is convenient to recall that Helmholtz problem is not well posed for a number of exciting frequencies corresponding to the resonant modes. The corresponding values of the wave number were avoided in the numerical experiments presented below. Firstly, optimal values for b and k are determined through numerical experiments carried out on a one-dimensional model problem. Then, using b and k functions previously computed, more general situations are addressed to test the performance of the proposed discontinuous finite element method. In all examples linear and bilinear shape functions and exact 2 or 2 · 2 Gaussian integration are adopted. 4.1. Numerical determination of stabilization parameters Each choice of b and k conveys to a finite element solution uhDG ðb; kÞ. This flexibility of the proposed formulation will be explored aiming at improving the finite element approximations. The choice of b and k is such that the relative errors in the L2-norm and H 1J -seminorm are minimized in a set of numerical experiments. The optimal choice, which is not restricted to a single pair (b, k), will be adopted in more general situations. Consider the numerical solution of Eq. (1) in the interval (0, 1), with k2 = constant, f(x) = 0 and Dirichlet boundary conditions given by: u(0) = 1 and u(1) = cos(k). In this case, the exact solution is uex(x) = cos(kx). The relative errors in the L2norm and H 1J -seminorm, 2

2 kREDG ðb; kÞkL2 ðXÞ

¼

kuex  uhDG ðb; kÞkL2 ðXÞ 2

kuex kL2 ðXÞ

;

2

2 jREDG ðb; kÞjH 1 ðXÞ J

¼

juex  uhDG ðb; kÞjH 1 ðXÞ

with 2

kuex  uhDG ðb; kÞkL2 ðXÞ ¼

J

2

juex jH 1 ðXÞ ne Z X e¼1

2

juex  uhDG ðb; kÞjH 1 ðXÞ ¼ J

ne Z X e¼1

;

2

ðuex  uh;e DG ðb; kÞÞ dXe ;

Xe

Xe

2

jruex  ruh;e DG ðb; kÞj dXe þ

ne Z X e0 >e

Ce0 e

2

ðue  ue0 Þ dCee0 h2ee0

are computed for different choices of b and k. Fig. 1 shows the 3-D plot of relative errors versus b and k for k2 = 400, obtained using a uniform mesh of 40 elements and same values of both parameters for all interfaces. To highlight the performance of the proposed method, we observe that the relative errors in the L2-norm of the interpolant and continuous Galerkin finite element solution are 0.0207 and 0.125, respectively. For some regions of the bk plane, the accuracy of the discontinuous solution is very close to that achieved by the interpolant, which can be considered as target to a finite element method free of pollution. It is important to observe that those regions encompass positive and/or negative values of the parameters. Similar conclusion can be drawn from Fig. 1(b) that depicts H 1J -seminorm of the error. The only point to be mention is that in this last case the regions of good accuracy are larger. Fig. 2 presents the same study for a higher wave number, k2 = 4000, for which 160 elements were used. Observing Figs. 1 and 2, it is possible to note that as the wave number k increases the region of optimal choice, in the bk plane, in which the error of the DG approximation is close to the error of the interpolant, becomes narrower. Restricting our search for optimal values to the region where b and k are positive, the numerical experiments indicate that bopt and kopt, the optimal values for those functions, depend on two dimensionless parameters: l1 ¼ k e hee0 and l2 = 0 0 keL, where hee0 ¼ minfhei ; hei g; hei ¼ maxx;y2Xe jxi  y i j; hei ¼ maxx;y2Xe0 jxi  y i j and L = max{Li}, Li = maxx,y2X jxi  yij (i = 1, . . . , domain dimension). The dependence of the optimal values on l1 and l2 is initially investigated by fixing l1 ¼ k e hee0 ¼ p5  0:628 and plotting both bopt and kopt as a functions of l2 in Fig. 3. The curve presented in that figure is divided in three regions leading to the following interpolation using Lagrange polynomials:

G.B. Alvarez et al. / Comput. Methods Appl. Mech. Engrg. 195 (2006) 4018–4035

4023

Fig. 1. Relative error REDG for k2 = 400, 300 6 b, k 6 300, Db = 10, Dk = 10. (a) L2-norm, (b) H 1J -seminorm.

Fig. 2. Relative error REDG for k2 = 4000, 300 6 b, k 6 300, Db = 10, Dk = 10. (a) L2-norm, (b) H 1J -seminorm.

kopt

8 2 P 1 > > > C k ðiÞgi ðnÞ; > > i¼1 > > >

i¼1 > > > > > 2 > > P C 3 ðiÞg ðnÞ; : k

i

l2 ¼ 20:5;

Dl2 ¼ 39

if 1 6 l2 6 40;

l2 ¼ 160;

Dl2 ¼ 240

if 40 < l2 6 280;

l2 ¼ 650;

Dl2 ¼ 740

if 280 < l2 6 1020;

l2 ¼ 20:5;

Dl2 ¼ 39

if 1 6 l2 6 40;

l2 ¼ 160;

Dl2 ¼ 240

if 40 < l2 6 280;

l2 ¼ 650;

Dl2 ¼ 740

if 280 < l2 6 1020;

ð17Þ

i¼1

bopt

8 2 P 1 > > C b ðiÞgi ðnÞ; > > > i¼1 > > >

i¼1 > > > > > 2 > > P C 3 ðiÞg ðnÞ; : b

i

ð18Þ

i¼1

where n ¼ Dl2 2 ðl2  l2 Þ and gi(n) are the usual Lagrange polynomials. It should be emphasized that although l2 = 1020 was used as superior limit to calculate the coefficients, Fig. 3 indicates that bigger values of l2 should verify the same linear

4024

G.B. Alvarez et al. / Comput. Methods Appl. Mech. Engrg. 195 (2006) 4018–4035

30 25 20 15 10

Beta Lambda

5 0

0

200

400

µ2

600

800

1000

Fig. 3. Dependence of bopt and kopt with l2.

Table 1 Coefficients determined by least squares fitting i

C 1k

C 2k

C 3k

C 1b

C 2b

C 3b

CP

Cfb

1 2 3 4 5 6

0.8856 10.8663 – – – –

10.7138 20.8501 21.7504 22.9263 23.8373 24.1142

24.4312 25.7294 – – – –

10.6091 12.0328 – – – –

11.6889 21.3005 22.2751 23.4213 24.3291 24.6018

24.9312 26.2294 – – – –

1.0435 0.9774 1.0002 – – –

9.2814 2.2220 2.0108 7.2884 0.5120 –

dependence. The coefficients C jk ðiÞ and C jb ðiÞ determined by least squares fitting are presented in Table 1. To find the dependency of b and k functions with l1, define: Polðl1 Þ ¼

k=kopt ; b=bopt

bðl1 ; l2 Þ ¼ bopt ðl2 Þfb ðl1 Þ;

ð19Þ

kðl1 ; l2 Þ ¼ kopt ðl2 Þfk ðl1 Þ ¼ kopt ðl2 Þfb ðl1 Þ Polðl1 Þ;

ð20Þ

where Pol(l1) and fb(l1) are fitted by Lagrange polynomials through numerical experiments. The numerical fitting showed that degree of polynomial lager than 2 for Pol(l1) or lager than 4 for fb(l1), had very little influence on the relative error in L2-norm. Then the following approximation is introduced Pol(l1) and fb(l1) as Polðl1 Þ ¼

3 X

C P ðiÞgi ðnÞ;

l1 ¼ 0:33148806;

Dl1 ¼ 0:60018177;

ð21Þ

C f b ðiÞgi ðnÞ;

l1 ¼ 0:33148806;

Dl1 ¼ 0:60018177.

ð22Þ

i¼1

fb ðl1 Þ ¼

5 X i¼1

4.2. Convergence study Figs. 4 and 5 present the plot of the relative errors of continuous (CG) and discontinuous (DG) approximations versus the mesh parameter h for three values of k. The curves (CI) corresponding to the error of the nodal interpolant are also included for comparison. In all cases b and k are determined by (17)–(22). The pollution effect, for large values of k, is detected by the presence of a preasymptotic range in the convergence curves, much more significant for the C0 Galerkin approximations. The critical number of degrees of freedom for the Galerkin finite element error can be estimated by qffiffiffiffi k3 [3], whereas for the interpolant this value is N I ¼ pk . It is important to note that in the preasymptotic N CG ¼ 24 range the error of continuous finite element is not ruled by the magnitude of kh. Fig. 5 indicates that the error behavior of the nodal interpolant (CI) and of the discontinuous Galerkin solution (CG) are very close for this three values of k.

G.B. Alvarez et al. / Comput. Methods Appl. Mech. Engrg. 195 (2006) 4018–4035

4025

1.E+02

1.E+02

1.E+01 1.E+01

1.E-01

1.E-02

1.E-03

1.E-04

1.E-05 1.E-01

1.E+00

relative error

relative error

1.E+00

1.E-01 CG-k1 CG-k2 CG-k3 CI-k1 CI-k2 CI-k3

1.E-02

1.E-02

(a)

1.E-03

CG-k1 CG-k2 CG-k3 CI-k1 CI-k2 CI-k3

1.E-03 1.E-01

1.E-04

1.E-02

1.E-03

(b)

h

1.E-04

h

Fig. 4. Convergence behavior of the relative errors of continuous Galerkin method for k2 equal to 400 (CG-k1), 4000 (CG-k2) and 40,000 (CG-k3) compared to the corresponding error of the continuous interpolant (CI): (a) L2-norm, (b) H 1J -seminorm.

1.E+01

10

1.E+00 1

relative error

relative error

1.E-01

1.E-02

0.1

1.E-03

1.E-04

1.E-05 0.1 (a)

CI-k1 CI-k2 CI-k3 DG-k1 DG-k2 DG-k3

0.01

0.01 h

0.001

CI-k1 CI-k2 CI-k3 DG-k1 DG-k2 DG-k3

0.001 0.1

0.01

(b)

0.001

h 2

Fig. 5. Convergence behavior of the relative errors of discontinuous Galerkin method for k equal to 400 (DG-k1), 4000 (DG-k2) and 40,000 (DG-k3) compared to the corresponding error of the continuous interpolant (CI): (a) L2-norm, (b) H 1J -seminorm.

In Figs. 6 and 7 the same convergence study with respect to the mesh refinement is presented, but now parameterized by B = k2h2. This allows a better comparison between the discontinuous solution and the nodal interpolation, as the last satisfies the following error estimates if u 2 H2(X) [28,29]:

4026

G.B. Alvarez et al. / Comput. Methods Appl. Mech. Engrg. 195 (2006) 4018–4035

1.E+02

1.E+02 1.E+01

1.E+01

relative error

1.E+00

relative error

1.E-01 1.E-02

1.E+00

1.E-01

1.E-03

CI-k1 CI-k2 CI-k3 CG-k1 CG-k2 CG-k3

1.E-04 1.E-05

1.E-06 1.E+03 1.E+01 1.E-01

1.E-02

1.E-03

1.E-03 1.E+03 1.E+01

1.E-05

(a)

CI-k1 CI-k2 CI-k3 CG-k1 CG-k2 CG-k3

1.E-01

1.E-03

1.E-05

(b)

Fig. 6. Convergence behavior of the relative errors of continuous Galerkin method for k2 equal to 400 (CG-k1), 4000 (CG-k2) and 40,000 (CG-k3) compared to the corresponding error of the Continuous Interpolant (CI): (a) L2-norm, (b) H 1J -seminorm.

1.E+01

1.E+01

1.E+00 1.E+00

relative error

relative error

1.E-01

1.E-02

DG-k1 DG-k2 DG-k3 CI-k1 CI-k2 CI-k3

1.E-03

1.E-04

1.E-01

1.E-02

DG-k1 DG-k2 DG-k3 CI-k1 CI-k2 CI-k3

1.E-03

1.E-05 1.00E+02 1.00E+00 1.00E-02 1.00E-04 (a)

1.E+02

1.E+00

1.E-02

1.E-04

(b)

Fig. 7. Convergence behavior of the relative errors of discontinuous Galerkin method for k2 equal to 400 (DG-k1), 4000 (DG-k2) and 40,000 (DG-k3) compared to the corresponding error of the continuous interpolant (CI): (a) L2-norm, (b) H 1J -seminorm.

kuex  uhI kL2 ðXÞ kuex kL2 ðXÞ juex  uhI jH 1 ðXÞ J

juex jH 1 ðXÞ

6 C 1 k 2 h2 ; 6 C 2 kh;

J

where C1 and C2 are constants not depending on k or h. From the curves in Figs. 6 and 7, the following expressions, which are verified for both nodal interpolation and the discontinuous finite element errors, are determined:

G.B. Alvarez et al. / Comput. Methods Appl. Mech. Engrg. 195 (2006) 4018–4035

kuex  uhI kL2 ðXÞ kuex kL2 ðXÞ

kuex  uhDG kL2 ðXÞ kuex kL2 ðXÞ

juex  uhI jH 1 ðXÞ

A

¼ C 1I ðk 2 h2 Þ 1I ;

J

juex jH 1 ðXÞ

4027

A

¼ C 2I ðk 2 h2 Þ 2I ;

J

¼ C 1DG ðk 2 h2 Þ

A1DG

juex  uhDG jH 1 ðXÞ J

;

juex jH 1 ðXÞ

¼ C 2DG ðk 2 h2 Þ

A2DG

.

J

The coefficients present in the above expressions are shown in Table 2. It is worth to call attention to the outstanding performance of the discontinuous method in the present example. No significant pollution is observed in the convergence curves and optimal rates of convergence are achieved for reasonable meshes. Table 2 Convergence rates determined by numerical experiments k2

C1I

A1I

C2I

A2I

C1DG

A1DG

C2DG

A2DG

400 4000 40,000

0.0831 0.0828 0.0827

0.9996 0.9988 0.9977

0.2937 0.2893 0.2867

0.4998 0.4993 0.4986

0.0831 0.0875 0.0830

0.9997 1.0103 0.9987

0.2937 0.2896 0.2867

0.4998 0.4994 0.4986

1.05

1.01

0.525

0.505

0

0 0

0.2

0.4

0.6

1

0.8

0

0.2

0.4

0.6

1

0.8

-0.505

-0.525

-1.01

-1.05 EXACT

CG

EXACT

DG

Fig. 8. Solution of homogeneous problem in one dimension k2 = 400. Continuous Galerkin (CG) kh = 0.5 (left) and discontinuous Galerkin (DG) kh = 0.625 (right).

10

1.01

8 6 0.505 4 2 0

0

0.2

0.4

0.6

0.8

1

0

0

0.2

0.4

0.6

0.8

1

-2 -4 -0.505 -6 -8 -10

EXACT

CG

-1.01

EXACT

DG

Fig. 9. Solution of homogeneous problem in one dimension k2 = 4000. Continuous Galerkin (CG) kh = 0.5 (left) and discontinuous Galerkin (DG) kh = 0.625 (right).

4028

G.B. Alvarez et al. / Comput. Methods Appl. Mech. Engrg. 195 (2006) 4018–4035

Figs. 8–10 show the pressure distribution along the domain computed by both continuous and discontinuous finite elements, for k2 = 400, 4000 and 40,000. The corresponding exact solutions are also plotted for comparison. In all approximations the rule of the thumb was followed, with the worst mesh resolution, only used for the discontinuous elements, given by kh  0.625, what is equivalent to ten elements per wavelength. As expected, spurious dispersion of the continuous method is observed, this degradation of the solution quality is a manifestation of the pollution effect. It is important to reinforce that the proposed method is able to capture the correct phase and amplitude even for very coarser meshes. All results presented up to now were obtained using b and k given by (17)–(22). However, in Figs. 1 and 2 is observed that when b and k increase keeping a certain relationship between them, the region where the functional of the error attain its minimum is enlarged. Fig. 11 shows the relationship between b and k for which jREDG  REIj 6 103, using a mesh that verifies the rule of thumb. Observing Fig. 11, it is possible to conclude that given a b P bopt fixed and k1 6 k 6 k2 then the error of the discontinuous Galerkin method will be very close to the corresponding error of the nodal interpolant. Fig. 12 shows the dependence of Dk = k2  k1 with b for two values of k: k2 = 4000, 40,000. 1.01

30

20 0.505 10

0

0 0

0.2

0.4

0.6

0.8

0

1

0.2

0.4

0.6

0.8

1

-10 -0.505 -20

-1.01

-30 EXACT

CG

EXACT

DG

2

Fig. 10. Solution of homogeneous problem in one dimension k = 40,000. Galerkin kh = 0.395 (left) and DGFEM kh = 0.625 (right).

1.E+06 lambda

lambda

1.E+06

1.E+06

1.E+06 lambda1

lambda1

lambda2

lambda2

8.E+05

8.E+05

6.E+05

6.E+05

4.E+05

4.E+05

2.E+05

2.E+05 beta

beta

(a)

0.E+00 0.E+00 2.E+05 4.E+05 6.E+05 8.E+05 1.E+06

(b)

0.E+00 0.E+00 2.E+05 4.E+05 6.E+05 8.E+05 1.E+06

Fig. 11. Relation between b and k where the minimum error is obtained in L2-norm, k2 = 4000 (a) and k2 = 40,000 (b).

G.B. Alvarez et al. / Comput. Methods Appl. Mech. Engrg. 195 (2006) 4018–4035

4029

1.E+04

delta lambda

1.E+05

k2=40000 k2=4000

1.E+03

1.E+02

1.E+01 beta 1.E+00 1.E+00 1.E+01 1.E+02 1.E+03 1.E+04 1.E+05 1.E+06

Fig. 12. Dependence of Dk = k2  k1 with b: k2 = 4000, 40,000.

Despite leading to a non-symmetric matrix and to larger number of degrees of freedom and bandwidth compared to the continuous formulation, the discontinuous formulation is competitive in terms of computational costs, especially for large values of the wave number k. This can be attributed to the reduction of the pollution effect. This fact was observed in a

Table 3 Number of degrees of freedom to obtain a 1% relative error in L2-norm and 10% in H 1J -seminorm k2

NI (L2 and H 1J )

NDG (L2 and H 1J )

NCG (L2)

N CG ðH 1J Þ

400 4000 40,000

57 184 576

114 368 1152

139 792 >3000

66 316 1500

1.E+02

1.E+02

1.E+01 1.E+01

relative error

relative error

1.E+00

1.E-01

1.E-02

1.E+00

1.E-01

1.E-03

1.E-04

1.E-05 1.0E-01 (a)

CG-k1 CG-k2 CG-k3 CI-k1 CI-k2 CI-k3

1.E-02

1.0E-02 h

1.E-03 1.E-01

1.0E-03 (b)

CG-k1 CG-k2 CG-k3 CI-k1 CI-k2 CI-k3 1.E-02 h

1.E-03

Fig. 13. Inhomogeneous problem. Relative errors of continuous Galerkin method for k2 equal to 400 (CG-k1), 4000 (CG-k2) and 40,000 (CG-k3) compared to the corresponding error of the continuous interpolant (CI): (a) L2-norm, (b) H 1J -seminorm.

4030

G.B. Alvarez et al. / Comput. Methods Appl. Mech. Engrg. 195 (2006) 4018–4035

1.E+01

1.E+01

1.E+00

1.E+00

relative error

relative error

number of numerical experiments. Representative results of this study are summarized in Table 3, where, for some values of the wave number, the minimum number of degrees of freedom needed to achieve a fixed a-priori error is presented. Each column corresponds to a different approximation, namely: interpolant (NI), discontinuous Galerkin (NDG) and continuous Galerkin (NCG). Now, we consider an inhomogeneous problem still in one dimension with k2 constant, the inhomogeneity generated by the distributed source f(x) = 2  k2x2 and the following Dirichlet boundary conditions: u(0) = 1 and u(1) = 1 + cos(k). The analytical solution to this problem is uex(x) = x2 + cos(kx). Fig. 13 presents plots of the relative errors of the continuous Galerkin method in L2-norm and H 1J -seminorm versus h for k2 = 400 (CG-k1), k2 = 4000 (CG-k2) and k2 = 40,000 (CG-k3) compared to the corresponding errors of the continuous interpolant (CI). Again, we observe the strong pollution efect on the continuous Galerkin solution. Same plots are presented in Fig. 14 for the discontinuous Galerkin solution with

1.E-01

1.E-02

(a)

1.E-02

CI-k1 CI-k2 CI-k3 DG-k1 DG-k2 DG-k3

1.E-03 1.E-01

1.E-01

1.E-02

1.E-03 1.E-01

1.E-03

h

(b)

CI-k1 CI-k2 CI-k3 DG-k1 DG-k2 DG-k3

1.E-02

1.E-03

h

Fig. 14. Inhomogeneous problem. Relative errors of discontinuous Galerkin method for k2 equal to 400 (DG-k1), 4000 (DG-k2) and 40,000 (DG-k3) compared to the corresponding error of the continuous interpolant (CI): (a) L2-norm, (b) H 1J -seminorm.

2

2

1.5

1.5

1

1

0.5

0.5

0

0

-0.5

-0.5

-1

-1 0

0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 EXACT CG

0

0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 EXACT DG

Fig. 15. Solution of inhomogeneous problem in one dimension k2 = 4000. Continuous Galerkin (left) and DG FEM (right).

G.B. Alvarez et al. / Comput. Methods Appl. Mech. Engrg. 195 (2006) 4018–4035

4031

b and k functions determined by (17)–(22). Observing Fig. 14, one can verify once more that the error of the nodal interpolant and discontinuous Galerkin solution are very close for the three values of k. Fig. 15 presents the solutions for k2 = 4000, obtained with the continuous (left) and discontinuous (right) Galerkin methods with a mesh of 80 elements (kh = 0.79). 4.3. Plane wave propagation—minimizing the pollution effect The discontinuous finite element formulation was tested in a more challenging situation concerning the propagation of a plane wave in a arbitrary direction not necessarily aligned to the mesh. This test is decisive to examine the dispersion properties of the proposed formulation. The problem given by Eq. (1) is considered now in a square domain of unity sides, k2 constant, f(x, y) = 0 and the following Dirichlet boundary conditions: uð0; yÞ ¼ cosðkðy sin hÞÞ; uðx; 0Þ ¼ cosðkðx cos hÞÞ;

uð1; yÞ ¼ cosðkðcos h þ y sin hÞÞ; uðx; 1Þ ¼ cosðkðx cos h þ sin hÞÞ;

whose exact solution is the real part of a plane wave propagating in the h-direction: u(x, y) = cos(k(x cos h + y sin h)). Fig. 16 presents the relative errors in L2-norm and H 1J -seminorm for k2 = 400, as a function of h, corresponding to the 2.50

1.60 CI DG1 DG2 CG GLS

1.20 1.00

1.50

relative error

relative error

2.00

CI DG1 DG2 CG GLS

1.40

1.00

0.80 0.60 0.40

0.50 0.20 0.00

0.00 (a)

0

22.5

45

67.5

90

(b)

2.5

22.5

45

67.5

90

1.6 CI DG1 DG2 CG GLS

CI DG1 DG2 CG GLS

1.4 1.2 relative error

2.0

1.5 relative error

0

1.0

1.0 0.8 0.6 0.4

0.5 0.2 0.0

0.0 (c)

0

22.5

45

67.5

90

(d)

0

22.5

45

67.5

90

Fig. 16. Relative error of the discontinuous Galerkin solution (DG) compared to the continuous interpolant (CI), continuous Galerkin (CG) and (GLS) in L2-norm (a, b) and H 1J -norm (c, d) as a function of h-direction: k2 = 400, (a, c) kh = 1, coarse mesh, (b, d) kh = 0.62, resolvable mesh.

4032

G.B. Alvarez et al. / Comput. Methods Appl. Mech. Engrg. 195 (2006) 4018–4035

continuous interpolant (CI), continuous Galerkin (CG), Galerkin Least Squares (GLS) and the discontinuous Galerkin method with bopt and kopt determined by (17)–(22) (DG1), and b and k chosen such that b1 = 105  bopt and k1 = 108,685  kopt verifying the stability relation established in Fig. 11 (DG2). In case (a) the mesh is coarse kh > 0.62, in case (b) the mesh verifies the rule of thumb kh = p/5  0.62. Same results are presented in Fig. 17 for k2 = 4000. Whenever the mesh verifies the rule of thumb (case (b)) starting from a b = bopt and k = kopt, one can choose b  bopt and k  kopt keeping the stability relation between them. That is, b and k belong to the region where the functional of the errors in L2-norm and H 1J -seminorm attain small values. For the coarse mesh (case (a)) the previous statement is not valid as show the Figs. 16(a) and 17(a). In this case, it is necessary what b < bopt and k < koptto get the relative error of the discontinuous Galerkin solution close to the relative error of the nodal interpolant. Fig. 18 shows sound pressure distributions in sections x = 0.505 (top) and y = 0.505 (bottom) corresponding to the continuous (left) and discontinuous (right) Galerkin finite element methods compared to the exact solution. The results were obtained with (101 · 101) mesh for h ¼ 3p , that is the h-direction, which corresponds to the largest ‘‘phase’’ error for the 8 discontinuous Galerkin method. The b and k functions were, once again, chosen by using (17)–(22). Those results confirm the good performance of the proposed method.

0.80

3.50

0.70

3.00

0.60

relative error

relative error

2.50

2.00 CI DG1 DG2 CG GLS

1.50

1.00

0.20 0.10

0.00 0

22.5

45

67.5

90

0.00 0 0.80

3.0

0.70

22.5

45

67.5

90

67.5

90

0.60 relative error in H1-norm

relative error

(b)

3.5

2.5 CI DG1 DG2 CG GLS

2.0

1.5

1.0

0.50 CI DG1 DG2 CG GLS

0.40 0.30 0.20

0.5

0.10 0.00

0.0 (c)

CI DG1 DG2 CG GLS

0.40 0.30

0.50

(a)

0.50

0

22.5

45

67.5

90

(d)

0

22.5

45

Fig. 17. Relative error of the discontinuous Galerkin solution (DG) compared to the continuous interpolant (CI), continuous Galerkin (CG) and (GLS) in L2-norm (a, b) and H 1J -norm (c, d) as a function of h-direction: k2 = 4000: (a, c) kh = 0.79, coarse mesh, (b, d) kh = 0.62, resolvable mesh.

G.B. Alvarez et al. / Comput. Methods Appl. Mech. Engrg. 195 (2006) 4018–4035

1.3

4033

1.01

0.78 0.505

0.26 0 0

0.2

0.4

0.6

0.8

1

0

0.2

0.4

0.6

0.8

1

-0.26 -0.505 -0.78

-1.01

-1.3 EXACT

CG

1.3

EXACT

DG

1.01

0.78 0.505

0.26 0 0

0.2

0.4

0.6

0.8

0

1

0.2

0.4

0.6

0.8

1

-0.26 -0.505 -0.78

-1.01

-1.3 EXACT

CG

EXACT

DG 2

Fig. 18. Solution of homogeneous problem in two dimension at x = 0.505 (top) and y = 0.505 (bottom) for k = 4000 and h ¼ 3p . Continuous Galerkin 8 (left) and discontinuous Galerkin (right).

5. Conclusions A discontinuous finite element method for Helmholtz equation is presented with C0 continuity on the interelement boundaries enforced in a weak sense depending on two parameters b and k which are crucial for the performance of the proposed method. These parameters are first determined numerically by solving a one dimension homogeneous Helmholtz equation with constant coefficient and Dirichlet boundary conditions. Then, the optimal choice of these parameters are adopted in more general situations in a two dimensional model problem. The one dimensional tests with appropriated values of the parameter b and k show a great accuracy of the proposed finite element formulation without phase error. The error is controlled by the magnitude of kh. That is, the practical ‘‘rule of the thumb’’ guarantees good accuracy of the approximate solution, independently of the wave number. A numerical study of the dispersion properties demonstrates the good performance of the discontinuous finite element method in two dimension problems. Using the optimal values of the parameters obtained in the one dimensional experiments, the pollution error is clearly reduced compared to other stabilized finite element methods. Numerical experiments show that the convergence behavior of the proposed method is very similar to the interpolant, in both L2-norm and H 1J -seminorm, even for large values of the wave number.

4034

G.B. Alvarez et al. / Comput. Methods Appl. Mech. Engrg. 195 (2006) 4018–4035

The numerical results presented indicate a good potential of the proposed formulation to solve Helmholtz problem in the mid and high frequency regime. It should be highlight that the good results are obtained within the framework of standard finite element interpolation. Moreover, stability and accuracy stem only from the added interelement edges terms. No additional stabilization term was required. Acknowledgement The authors wish to thank the Brazilian research funding agencies FAPERJ and CNPq for their support to this work. References [1] I. Harari, T.J.R. Hughes, Finite element method for the Helmholtz equation in an exterior domain: model problems, Comput. Methods Appl. Mech. Engrg. 87 (1991) 59–96. [2] Bayliss, C.I. Goldstein, E. Turkel, On accuracy conditions for the numerical computation of waves, J. Comp. Phys. 59 (1985) 396–404. [3] F. Ihlenburg, I. Babusˇka, Finite element solution of the Helmholtz equation with high wave number Part I: The h-version of the FEM, Comput. Math. Appl. 30 (9) (1995) 9–37. [4] A.K. Aziz, R.B. Kellogg, A.B. Stephens, A two point boundary value problem with a rapidly oscillating solution, Numer. Math. 53 (1988) 107–121. [5] J. Douglas Jr., J.E. Santos, D. Sheen, L. Schreiyer, Frequency domain treatment of one-dimensional scalar waves, Math. Models Methods Appl. Sci. 3 (2) (1993) 171–194. [6] I. Harari, T.J.R. Hughes, Galerkin/least squares finite element methods for the reduced wave equation with non-reflecting boundary conditions in unbounded domains, Comput. Methods Appl. Mech. Engrg. 98 (1992) 411–454. [7] C.I. Goldstein, The weak element method applied to Helmholtz type equations, Appl. Numer. Math. 2 (1986) 409–426. [8] Y.K. Cheung, W.G. Jin, O.C. Zienkiewicz, Solution of Helmholtz equation by Trefftz method, Int. J. Numer. Methods Engrg. 32 (1) (1991) 63–78. [9] I. Harari, T.J.R. Hughes, A cost comparison of boundary element and finite element methods for problems of time-harmonic acoustics, Comput. Methods Appl. Mech. Engrg. 97 (1992) 77–102. [10] L.L. Thompson, P.M. Pinsky, A Galerkin least squares finite element method for the two-dimensional Helmholtz equation, Int. J. Numer. Methods Engrg. 38 (3) (1995) 371–397. [11] F. Ihlenburg, I. Babusˇka, Dispersion analysis and error estimation of Galerkin finite element methods for the numerical computation of wave, Int. J. Numer. Methods Engrg. 38 (1995) 3745–3774. [12] I. Babusˇka, F. Ihlenburg, E.T. Paik, S.A. Sauter, A generalized finite element method for solving the Helmholtz equation in two dimensions with minimal pollution, Comput. Methods Appl. Mech. Engrg. 128 (1995) 325–359. [13] I. Babusˇka, F. Ihlenburg, T. Strouboulis, S.K. Gangaraj, A posteriori error estimation for finite element solutions of Helmholtz equation—Part II: Estimation of the pollution error, Int. J. Numer. Methods Engrg. 40 (21) (1997) 3883–3900. [14] I. Babusˇka, S.A. Sauter, Is the pollution effect of the FEM avoidable for the Helmholtz equation considering high wave numbers, SIAM J. Numer. Anal. 34 (1997) 2392–2423. [15] L.P. Franca, C. Farhat, A.P. Macedo, M. Lesoinne, Residual-free bubbles for the Helmholtz equation, Int. J. Numer. Methods Engrg. 40 (1997) 4003–4009. [16] I. Harari, Reducing spurious dispersion, anisotropy and reflection in finite element analysis of time-harmonic acoustics, Comput. Methods Appl. Mech. Engrg. 140 (1/2) (1997) 39–58. [17] L.P. Franca, A.P. Macedo, A two-level finite element method and its application to the Helmholtz equation, Int. J. Numer. Methods Engrg. 43 (1) (1998) 23–32. [18] M. Stojek, Least-squares Trefftz-type elements for the Helmholtz equation, Int. J. Numer. Methods Engrg. 41 (5) (1998) 831–849. [19] 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 (1998) 255–299. [20] P. Bouillard, Influence of the pollution on the admissible/field error estimation for FE solutions of the Helmholtz equation, Int. J. Numer. Methods Engrg. 45 (7) (1999) 783–800. [21] K. Gerdes, F. Ihlenburg, On the pollution effect in FE solutions of the 3D-Helmholtz equation, Comput. Methods Appl. Mech. Engrg. 170 (1/2) (1999) 155–172. [22] J.L. Cipolla, Subgrid modeling in a Galerkin method for the Helmholtz equation, Comput. Methods Appl. Mech. Engrg. 177 (1999) 35–49. [23] P. Monk, Da-Qing Wang, A least-squares method for the Helmholtz equation, Comput. Methods Appl. Mech. Engrg. 175 (1999) 121–136. [24] A.A. Oberai, P.M. Pinsky, A residual-based finite element method for the Helmholtz equation, Int. J. Numer. Methods Engrg. 49 (2000) 399–419. [25] C. Farhat, A.P. Macedo, M. Lesoinne, A two-level domain decomposition method for the iterative solution of high frequency exterior Helmholtz problems, Numer. Math. 85 (2000) 283–308. [26] C. Farhat, I. Harari, L. Franca, The discontinuous enrichment method, Comput. Methods Appl. Mech. Engrg. 190 (2001) 6455–6479. [27] C. Farhat, I. Harari, U. Hetmaniuk, A discontinuous Galerkin method with Lagrange multipliers for the solution of Helmholtz problems in the midfrequency regime, Comput. Methods Appl. Mech. Engrg. 192 (2003) 1389–1419. [28] G. Strang, G.J. Fix, An Analysis of the Finite Element Method, Prentice Hall, Englewood Cliffs, NJ, 1973. [29] F. Ihlenburg, Finite Element Analysis of Acoustic Scattering, Applied Mathematical Sciences, vol. 132, Springer-Verlag, New York, 1998. [30] J. Douglas Jr., T. Dupont, Interior Penalty Procedures for Elliptic and Parabolic Galerkin Methods, Lecture Notes in Physics, vol. 58, SpringerVerlag, Berlin, 1976. [31] G.A. Baker, Finite element methods for elliptic equations using nonconforming elements, Math. Comput. 31 (1977) 45–59. [32] M.F. Wheeler, An elliptic collocation-finite element method with interior penalties, SIAM J. Numer. Anal. 15 (1) (1978) 152–161. [33] P. Percell, M.F. Wheeler, A local residual finite element procedure for elliptic equations, SIAM J. Numer. Anal. 15 (4) (1978) 705–714. [34] L.M. Delves, C.A. Hall, An implicit matching principle for global element calculations, J. Inst. Math. Appl. 23 (1979) 223–234. [35] D.N. Arnold, An interior penalty finite element method with discontinuous elements, SIAM J. Numer. Anal. 19 (4) (1982) 742–760. [36] J.T. Oden, I. Babusˇka, C.E. Baumann, A discontinuous hp finite element method for diffusion problems, J. Comp. Phys. 146 (1998) 491–519.

G.B. Alvarez et al. / Comput. Methods Appl. Mech. Engrg. 195 (2006) 4018–4035

4035

[37] I. Babusˇka, C.E. Baumann, J.T. Oden, A discontinuous hp finite element method for diffusion problems: 1-D analysis, Comput. Math. Appl. 37 (10) (1999) 103–122. [38] C.E. Baumann, J.T. Oden, A discontinuous hp finite element method for convection-diffusion problems, Comput. Methods Appl. Mech. Engrg. 175 (1999) 311–341. [39] C.E. Baumann, J.T. Oden, A discontinuous hp finite element method for the Euler and Navier–Stokes equations, Int. J. Numer. Methods Fluids 31 (1999) 79–95. [40] D.N. Arnold, F. Brezzi, B. Cockburn, D. Marini, Discontinuous Galerkin methods for elliptic problems, Technical Report of Instituto di Analisi Numerica del C.N.R., Pavia, Italy (1999). [41] E.G. Dutra do Carmo, A.V.C. Duarte, A discontinuous finite element-based domain decomposition method, Comput. Methods Appl. Mech. Engrg. 190 (2000) 825–843. [42] A.V.C. Duarte, E.G. Dutra do Carmo, F.A. Rochinha, Consistent discontinuous finite elements in elastodynamics, Comput. Methods Appl. Mech. Engrg. 190 (2000) 193–223. [43] A.V.C. Duarte, F.A. Rochinha, E.G. Dutra do Carmo, Discontinuous finite element formulations applied to cracked elastic domains, Comput. Methods Appl. Mech. Engrg. 185 (2000) 21–36. [44] E.G. Dutra do Carmo, A.V.C. Duarte, New formulations and numerical analysis of discontinuous Galerkin methods, Comput. Appl. Math. 21 (3) (2002) 661–715.