A Combined Finite Element Scheme for Second Elliptic Problems Posted in Domains with Rough Boundaries

ZHOU Jianan (周迦南), HUANG Lei (黃 磊), XU Shipeng(徐世鵬)

School of Mathematics and Computer Science, Jiangxi Science and Technology Normal University, Nanchang 330038, China

Abstract: A finite element method (FEM) is proposed to solve second elliptic problems posted in domains with rough boundaries. We also present an FEM for elliptic problems as an extension of the original method. The optimal error estimates in the energy norm are also proved. Numerical results for second elliptic equations in regular and rough domains are presented to illustrate the eciency of our method, respectively, which is consistent with theoretical analysis.

Key words: finite element; multiscale problem; elliptic problem; rough domain; penalty technique

Introduction

A combined finite element method (FEM) which has been proposed in Ref.[1] is introduced to handle problems in oscillating domains[2]. It uses a finite mesh with sizehin the vicinity of oscillating boundaries and a coarse mesh with sizeH(?h) for the interior subdomains. Owing to the discontinuity of the interface of the coarse mesh and the finite mesh, the transmission conditions across the fine-coarse mesh interface are treated by the Nitsche’s method[3]. In the definition of the bilinear form, the weighted averages of the function values from the fine and the coarse girds replace the arithmetic averages respectively. Moreover, the penalty coefficient is defined asγ/(H+h) for a positive constantγ, and then the problem that the ratioH/hdebases the convergence rate can be disposed of[4-5]. An optimal convergence in terms of elements is obtained since for problems in the domain with rough boundaries their solutions are generally inHs(Ω) with 1

In this paper we also try to propose an FEM to cope with problems in rough domains, and the main idea of the FEM can be regarded as an extension of the one proposed in Ref. [1]. The Nitsche’s method and the weighted average in the definition of the bilinear form are still employed to ensure the stability of numerical scheme. The main dierence between the two consists in that our method has one more adjustable parameter, since we apply penalty technique into both the function values from the fine and the coarse girds respectively and their numerical fluxes, just as a discontinuous FEM. In addition, the penalty coefficients are defined asσ1/Handσ2Hfor positive constantsσ1andσ2, respectively. An optimal convergence in terms of elements can be also obtained.

There are many other numerical methods to deal with the complex geometrical boundaries such as homogenization theory (seen in Refs.[6-9] and the references therein), the multiscale FEM (MsFEM)[10-12], the cut FEM[13], the extended/generalized FEM[14], and the composite finite elements[15]. We can refer to Ref. [1] for more comprehensive introduction about the methods handling the complex geometrical boundaries. Readers can also refer to Refs. [16-22] for solving elliptic problems with different FEMs.

The rest of the paper is organized as follows. In section 1 we present the model problem and recall some notations which are necessary for theoretical analysis. In section 2 our FEM is formulated, and we analyze the continuity, coercivity and convergence of the proposed method. In section 3 we deal with two model problems in regular and rough domains by our method, respectively. Numerical tests illustrate the validity and the efficiency of our method. Last but not least, some conclusions are included in section 4.

1 FEM

1.1 Model problem

Consider the second order elliptic equation as follows.

-·(Au)=finΩ,

u=0 on ?Ω，

(1)

whereΩ?Rd(d=2, 3) is a bounded Lipschitz domain. Assume that

f∈L2(Ω),

and the matrix

A(x)=[aij(x)]∈(W1, ∞(Ω))d×d,

which is symmetric and satisfies the uniformly elliptic condition,i.e., there exist two positive constantsλandΛsuch that

λ|ξ|2≤aij(x)ξiξj≤Λ|ξ|2, ?ξ∈Rd.

(2)

1.2 Some notations

For the sake of completeness, we recall some notations in Ref. [1] here. First separate domainΩinto subdomainsΩ1andΩ2such that

Ω2??Ω,

and

Ω=Ω1∪Ω2∪Γ,

whereΓ=?Ω1∩?Ω2, Lipschitz continuous.

DefinehK(orHK) as diam(K), for any elementK∈h, H. DenotehandHas the maximal element diameter overhandHrespectively.nis a unit normal vector onΓpointing fromΩ1toΩ2, whileniis the outward unit normal vector ofΩi,i=1, 2.

Let

Fors≥0 and any subsetU?h, H, denote

and

Letνbe a piecewise smooth function and

e=K1∩K2,e∈Γ,K1∈Ω1，K2∈Ω2.

Ifvi:=v|Ki,i=1, 2, then the weighted averages {ν}wand the jump [ν] ofνovereare defined as

{v}w=w1v1+w2v2, [ν]=v1-v2,

(3)

where the weight coefficients are

We also introduce the following “energy” space and corresponding finite element space:

V={ν∈L2(Ω):ν|Ωi∈H1(Ωi),i=1, 2,ν|?Ω=0,
ν|K∈H2(K), ?K∈h, H,K∩?！?},

(4)

and

Vh, H={v∈L2(Ω):ν|Ω1∈Vh,
ν|Ω2∈VH,ν|?Ω=0},

(5)

In what follows, we use the notationABandB?Ato representA≤CBandB≥CAwith a mesh-independent generic constantC>0, respectively, which can have different values in different formulas.

1.3 Formulations of the FEM

(6)

a(uh, H,vh, H)=F(vh, H), ?vh, H∈Vh, H,

(7)

where for givenσ1>0, 0<σ2≤C, ?uh, H,vh, H∈Vh, H,

Remark 1 It is also easy to see that the scheme is consistent,i.e.,

a(u,vh, H)=F(vh, H), ?vh, H∈Vh, H.

(8)

Define the corresponding energy norm for our scheme:

2 Error Estimates

In this section, we only give the proof of the stability and error estimates of the formulation with a symmetric coefficientβ=1. Forβ=1 andβ=0, the similar analysis can be obtained.

In excruciating（） agony, I often ponder this: if I could live my life again, I would never try to achieve the elimination12 of prejudices of any kind for the simple reason that there is a price to pay.

There are three inequalities( seen in Ref. [16]) used frequently in this paper, including trace inequality:

(9)

Young’s inequality:

(10)

and inverse inequality:

(11)

whereCtris a positive mesh-independent constant.

Forβ=1, we get the following results.

Lemma 1 There exists a positive mesh-independent constantσ0such that whenσ1≥σ0, it holds that

(12)

(13)

Proof Inequality (12) can be proven easily by Cauchy-Schwarz inequality. It remains to prove for inequality (13). Denotingvh=vh, H|Ω1,vH=vh, H|Ω2and using the inequalities (2) and (11), we have

(14)

Hence, applying the inequality (10) we have

a(uh, H,vh, H)=

Remark 2 As a extension of the one proposed in Ref. [1], the mesh-independent constantσ2≤Cwhich plays a role in the stability of the scheme (7) only needs to be positive.

The following lemma is analogue of theCealemma[16]for the FEM.

Lemma 2 Letuanduh, Hbe the solutions of problems (1) and (6), respectively. There exists a positive mesh-independent constantσ0independent ofHandhsuch that whenσ1≥σ0, it holds that

(15)

Proof Similar to the proof of the Lemma 5.3 in Ref. [18] and using the consistence (8) of the numerical scheme, it follows that

2[(f,uh, H-vh, H)-a(vh, H,uh, H-vh, H)]=

2[(f,uh, H-vh, H)-a(u,uh, H-vh, H)+
a(u-vh, H,uh, H-vh, H)]≤

Hence, it follows from the triangle inequality that

The proof of the lemma is completed.

In order to get a global interpolation operator, letIH:C(Ω2)VHbe the standard Lagrange one. It is well known that (seen in,e.g., Ref. [16])

|v-IHv|i, h≤Ch2-i|v|,
0≤i≤2, ?v∈H2(h).

(16)

However, since the exact solutionumay be singular near the corner points, we use the Scott-Zhang interpolation inΩ1instead of the Lagrange one. LetΠh:H1(Ω1)Vhbe the Scott-Zhang interpolation (seen in Refs. [17, 19]) so that the following estimates hold:

(17)

Define the global operatorIh, Hvia

(Ih, Hv)|Ω1=Πh(v|Ω1), (|h, Hv)|Ω2=IH(v|Ω2),

(18)

In view of Lemma 2, we have the following convergence result by the interpolation error of the operatorIh, H.

Theorem 1 Letuanduh, Hbe the solutions of problems (1) and (7), respectively. Assume that

and

u∈H1+s(K), ?K∈h, 0

Then there exists a positive mesh-independent constantσ0such that whenσ1≥σ0, it holds that

(19)

σ2H‖[A

By inequalities (16) and (17), we have

Next, by use of trace inequality (9), we have

H‖

Similarly, we obtain

CH‖

and

Combining the termsT1,T2,T3andT4together, we get the estimate

Finally, by Lemma 2, the proof of the theorem is completed.

3 Numerical Tests

In this section, we will report some numerical results for our proposed FEM to verify the efficiency of the theoretical result established in section 3. Since our method is very similar to that in Ref. [1], we only apply the FEM to solve the model problem in the domain with regular boundaries and the oscillating boundaries on all sides, respectively. In our numerical tests, we also consider the following Poisson equation[1]:

Δu=finΩ,

u=gon ?Ω,

(20)

wheref=1 andu=g:= (x2+y2)/4.

For our FEM, the triangulation may be done by the same way as that in Ref. [17] (seen in Fig. 1 and Fig. 2).

Fig. 1 Regular domain

Fig. 2 General domain

Denotinguh, Has the numerical solutions computed by the methods considered in this section, we measure the relative error inL2,L∞and energy norms as

3.1 Regular domain

To compare with the method in Ref. [1], we verify the eciency of our FEM. In order to do this, we consider the model problem (20) in the squared domainΩ=[0, 1]×[0, 1], the subdomains of which are

and

Ω1=ΩΩ2.

In our first test, we fixH/h=23and then compute the model problem with a series ofH(correspondingh). The results are listed in Table 1. We observe that the optimal convergence rate of the relative energy error may be attained. It is also noted that our results are similar to those in Ref. [1] to some extent.

The second test is to further verify the efficiency of our FEM. To do this, we fixh=2-10and carry on the numerical test with differentH. The results are shown in Table 2. It is easy to see that the optimal convergence rate of the relative energy error may be attained. Also note that our results are similar to those in Ref. [1] to some extent.

Table 1 Relative energy errors of our FEM for H/h=23, β=1, σ1=20, and σ2=0.1

Table 2 Relative errors of our FEM for h=2-10, β=1, σ1=20, and σ2=0.1

3.2 General domain

In this example, we consider two arbitrary oscillating domains with general boundaries on all sides (seen in Fig. 2 and Fig. 1 in Ref. [2]). For the first one, the oscillating extent of the boundary is 1/64. In both simulations, we always fix the fine mesh sizeh=2-9, and use our FEM to test the convergence rate aboutH. The results are illustrated in Table 3 and Table 4, which show the efficiency of our method. It turns out that the convergence rate of the relative energy error aboutHis optimal, which is not subject to the oscillation of boundaries. Notice that the results of Table 4 are less optimal than those in Table 3 due to the lower regularity. We also remark that our results are similar to those in Ref. [1] to some extent.

Table 3 Relative errors of our FEM for h=2-9, β=1, σ1=20, and σ2=0.1

Table 4 Relative energy errors of our FEM for h=2-9, β=1, σ1=20, and σ2= 0.1

4 Conclusions

In this paper, we have developed an FEM similar to that in Ref. [1]. A rigorous and careful analysis has been given for the elliptic equation with general coefficient and oscillating boundaries to show the consistence, stability, and convergence of the scheme. Since the global regularities of the solutions are low, the results of the convergence are expanded as the elements, from which we actually get the quasi-optimal convergence.

Numerical experiments are carried out for Poisson equation with regular boundaries and oscillating boundaries on all sides to verify the efficiency of the proposed method. Results show the efficiency of our proposed method.

We remark that the condition number of our FEM also depends onH,h,d, and the thickness of the fine-mesh regionΩ1, the similar proof of which we can refer to Appendix A in Ref. [1]. Moreover, developing the best choices of penalty parameters for the two FEMs to optimize the computational efficiency is our future work.