ANALYSIS OF FINITE ELEMENT DOMAIN ... - Semantic Scholar

ANALYSIS OF FINITE ELEMENT DOMAIN EMBEDDING METHODS FOR CURVED DOMAINS USING UNIFORM GRIDS SHENG ZHANG Abstract. We analyze the error of a finite element domain embedding method for elliptic equations on a domain ω with curved boundary. The domain is embedded in a rectangular domain R on which uniform mesh and linear continuous elements are employed. The numerical scheme is based on an extension of the differential equation from ω to R by regularization with a small parameter  (for Neumann and Robin problems), or penalty with a large parameter −1 (for Dirichlet problem), or a mixture of these (for mixed boundary value problem). For Neumann and Robin problems, we prove that when  ≤ h (that is the mesh size), the error in the H 1 (ω) norm is of the optimal order O(h). For Dirichlet problem, when  ≤ h1/2 , the error is O(h1/2 ) that is not optimal. If the mesh is adjusted around ∂ω by moving the near-by nodes onto it and reconnecting some nodes such that a polygonal interpolation of ∂ω is formed in the mesh, then the optimal convergence rate O(h) holds for Dirichlet problem if ω is convex and  ≤ h. If ω is not convex, then the convergence rate can only be improved to O(h2/3 ) by such mesh adjustment with the parameter being  = h2/3 . In this latter case, a parameter smaller than h2/3 thwarts the convergence rate. Key words. Domain embedding, fictitious domain, curved boundary, uniform mesh, fast Poisson solver. Subject classification. 65N30, 65N45, 65T50.

1. Introduction In recent years, domain embedding, or fictitious domain, methods have developed into a general methodology of numerical computations of partial differential equations that has the advantage of efficiently dealing with complex domain boundaries by using uniform, or slightly adjusted, grids on a larger rectangular domain [10]. For a second order elliptic boundary value problem on a complicated domain, the ultimate goal of domain embedding methods is to quickly obtain accurate numerical solution by mainly solving the Poisson equation on uniform meshes on a rectangle by fast solvers, through, say, preconditioned iterations. A first step in these methods is to extend the boundary value problem to a larger rectangle. There are many ways to do this, one of which is using regularization for Neumann or Robin problem, or penalty for Dirichlet problem, or a combination of these for mixed boundary value problem [16, 14, 22]. This results in an extended boundary value problem on the rectangle with large jumps in the coefficients of the differential equation. A finite element domain embedding method then is obtained as a straightforward discretization of the extended equation on a uniform mesh on the rectangle. The error of the finite domain embedding method has two sources. One is due to the replacement of the original problem by the one on the Department of Mathematics, Wayne State University, Detroit, Michigan 48202, USA. Email: [email protected]. This work was partially supported by NSF grant DMS-0513559. 1

2

SHENG ZHANG

larger rectangular domain. This has been analyzed by many authors and in [22, 23] the sharp estimates on the regularization/penalty error can be found. The other is due to the discretization of the extended equation, which seems not have been completely understood. In the literature, there are numerous works devoted to this subject. In [16], a thorough analysis is presented under the assumption that the original domain boundary aligns with the uniform mesh. This assumption seems a stringent restriction on the arbitrariness of the domain. For example, the domain boundary can not be smoothly curved. And the regularity and the alignment required in this analysis are often not compatible. Furthermore, the results are actually not valid for general curved domains for Dirichlet problems. In [13, 14], an analysis for Neumann problem can be found, which did not impose such restriction on the domain. But the error estimates established there are quite inaccurate. A difficulty in this kind of error analysis for curved and complex domains is that the solution of the extended equation is not globally smooth. It has jumps in its normal derivative across the original domain boundary that inevitably cuts through some element if the mesh is uniform. In this paper, we present an analysis on such discretization errors for various boundary value problems. We first establish the estimates under the mere assumption that the original domain is Lipschitz. We then assume the domain boundary is smoothly curved and prove sharp estimates, and use the aforementioned regularization/penalty error estimate to determine the balance between the mesh-size and regularization/penalty parameter such that the overall accuracy of the finite element domain embedding method is optimized [19, 18]. We shall see that the balance is rather delicate in some cases. Preconditioning by fast Poisson solvers is a crucial ingredient in the success of finite element domain embedding methods. In this view, a slight distortion of the uniform mesh is acceptable, as long as the discrete Laplacian on the distorted mesh is spectrally equivalent to that on the uniform mesh. It turns out that no such adjustment is needed for Neumann or Robin problem since the full accuracy is already achieved by uniform meshes. But for Dirichlet problem the accuracy offered by uniform mesh is rather poor, and a slight adjustment of the uniform mesh around the original domain boundary by moving the near-by nodes onto the boundary efficiently reduces a locking effect and renders the full accuracy to the domain embedding finite element method when the domain is convex, or significantly improves the accuracy when the domain is not convex. For the latter, some local mesh refinement around the concave portion of domain boundary could be employed to achieve the full accuracy. We shall give some guides for such refinement as well. The results of this paper is valid for general self-adjoint second order elliptic equations. Similar results can be derived for finite elements of higher orders. The theory is also possible to extend to some higher order equations like those of clamped and free (but not simply supported) plates. The results could be useful to numerical analysis of other differential equations with highly discontinuous coefficients like those arising in study of composite materials. For simplicity, we only present the results for Poisson equation defined on a two-dimensional domain, and for linear continuous finite element. In the remainder of this introduction, we briefly summarize our results and compare them with existing results we found in the literature.

FINITE ELEMENT DOMAIN EMBEDDING METHODS

3

Let ω ⊂ R2 be a bounded connected domain such that its boundary Γ is of C 2 class [12]. For a given function f ∈ L2 (ω), we consider the boundary value problem −∆u = f

in ω

subject to homogeneous Dirichlet boundary condition u = 0 on Γ, or homogeneous Neumann condition ∂u/∂n = 0 with n being the unit outward normal of Γ, or Robin condition ku + ∂u/∂n = 0 with k being a smooth, bounded, and strictly positive function on Γ, or mixed boundary condition. The homogeneity is not a real restriction for Neumann and Robin conditions. A non-homogeneous Dirichlet condition could be reduced to a homogeneous one by defining a function satisfying the boundary condition, and subtracting it from u. The finite element domain embedding method is based on an approximation to u obtained by solving an -dependent boundary value problem on a larger rectangular domain R ⊂ R2 such that ω ⊂ R. The formulation of the problem on R depends on the boundary condition in the original problem. The finite element domain embedding method is then a straightforward discretization of the extended one on R with mesh size h. We shall assume that ω ⊂⊂ R. Then the boundary condition on ∂R is at one’s disposal. We impose the homogeneous Dirichlet condition. However, Neumann, Robin, mixed, or periodical boundary condition on ∂R works equally well. (If ∂R ∩ Γ 6= ∅, the condition on ∂R must be subject to some restrictions [23].) We let Ω = R \ ω be the fictitious domain, and extend the function f to a function f¯ on R by defining f¯ = 0 on Ω. Henceforth, we shall use the notation P . Q which means that there exists a constant C independent of P, Q, , and h such that P ≤ CQ. The notation P ' Q means P . Q and Q . P. We use H s (D) to denote the L2 based Sobolev space of order s on a domain D, in which the norm and semi-norm are denoted by k · ks,D and | · |s,D , respectively. In particular, the L2 norm is denoted by k · k0,D or | · |0,D . We start with the Neumann problem of which the weak formulation seeks u ∈ H 1 (ω)/R such that Z (1.1) (∇u, ∇v)[L2(ω)]2 = f vdx ∀ v ∈ H 1 (ω). ω

Here ∇ is the gradient operator, and the parentheses stand for the inner products in Hilbert spaces indicated by the subscripts. This problem is well-posed under the assumption that R f dx = 0. For  > 0, we determine u ∈ H01 (R) by a regularization such that ω Z   (1.2) (∇u , ∇v)[L2(ω)]2 + (∇u , ∇v)[L2(Ω)]2 = f¯vdx ∀ v ∈ H01 (R). R

This is an elliptic problem that yields a unique u ∈ H01 (R). As  → 0, u converges to a limit u0 ∈ H01 (R) at the sharp rate of  as long as ω is a Lipschitz domain (which is assured by our assumption on Γ). The limit u0 is harmonic on Ω and u0 |ω solves the Neumann problem (1.1). The equivalent estimate [23] ku − u0 k1,R ' |u − u0 |1,ω ' 

holds with the rare exception that u ≡ u0 , which occurs if and only if u0 = 0 on Γ, i.e., a solution of the homogeneous Neumann problem also satisfies the homogeneous Dirichlet condition on Γ. Incidentally, this is also equivalent to that u0 ∈ H 2 (R). In this case, we can take  = 1 in (1.2). And the Poisson equation on R produces the exact solution of the

4

SHENG ZHANG

original equation on ω. Except for this special case, we do not have u0 ∈ H 2 (R), since there is a jump in its normal derivative across Γ.

Figure 1. A domain ω embedded in a rectangle R with a uniform triangulation and a triangulation adjusted around Γ. We introduce a uniform triangulation of mesh size h on R, see the left figure in Figure 1, and let Hh ⊂ H01 (R) be the subspace of piecewise linear continuous functions subordinate to this triangulation. The finite element domain embedding method is a straightforward discretization of (1.2). It determines uh ∈ Hh such that Z   ¯ h dx ∀ v ∈ Hh . (1.3) (∇uh , ∇vh )[L2 (ω)]2 + (∇uh , ∇vh )[L2 (Ω)]2 = fv R

uh

on ω is the numerical solution offered by the finite domain emThen the restriction of bedding method. We prove that
√ |uh − u0 |1,ω .  + inf |u0 − vh |1,ω + |u0 − vh |1,Ω . vh ∈Hh

This estimate is valid as long as ω is a Lipschitz domain. When Γ ∈ C 2 , we have the piecewise smoothness property of u0 that u0 |ω ∈ H 2 (ω) and u0 |Ω ∈ H 2 (Ω). Using this regularity, we √ 0 0 0 construct an interpolation vh of u such that |u − vh |1,ω . h and |u − vh |1,Ω . h, thus to prove that √ (1.4) |uh − u0 |1,ω .  +h +  h.

From this we see that when  . h, we have a domain embedding finite element method that achieves the full accuracy of linear continuous finite elements with an error √ of O(h). √ This result is better than that of [13, 14], in which an estimate of the form  + h + , instead of that in (1.4), was established. They, therefore, require  = Ch2 to achieve the full accuracy. The discrete system arising from (1.3) can be effectively solved by iterative methods using the discrete Laplacian on the uniform mesh on R as a preconditioner (which obviously removes the effect of h on the condition number) and starting with special initials such as those that are zero on nodes in Ω (which removes the effect of , see [4]). Theoretically, the convergence of such preconditioned iteration remains unchanged for various  values. But from a computational point of view, one would like to avoid very small , as the aforementioned  = Ch2 . In this sense,  = Ch seems the best choice.

FINITE ELEMENT DOMAIN EMBEDDING METHODS

5

The piecewise smoothness of u0 is a consequence of the assumption that Γ ∈ C 2 and f ∈ L2 (ω) [12]. It is not valid when ω is a polygon, since either u0 |ω ∈ H 2 (ω) or u0 |Ω ∈ H 2 (Ω) will be broken, depending on the convexity of ω. However, when ω is convex and Γ aligns with the mesh-line, as assumed in [16], the result (1.4) remains valid, which needs a different argument. If we fix h and let  → 0, we have uh → u0h ∈ Hh . The restriction of the limit u0h on ω is identical to the solution given by a method proposed by Friedrichs and Keller in 1966 [11], of which an analogue for Robin problem (see below) was used in the 1980’s by the Boeing program tranair to circumvent the mesh generation problem around a complete aircraft and was also included in the software freefem3d [10]. The positive value of  renders to the equation (1.3) the advantages of better data structure and faster resolution by using fast Poisson solvers on the rectangle R as preconditioners. This method also offers an alternative of dealing with the singularity of discretized Neumann problem. See [5], where a survey of various such techniques can be found. The weak formulation of Robin problem seeks u ∈ H 1 (ω) such that Z f vdx ∀ v ∈ H 1 (ω). (1.5) (∇u, ∇v)[L2 (ω)]2 + (ku, v)L2 (Γ) = ω

This equation can be extended to R in the following way. We determine u ∈ H01 (R) such that Z    (1.6) (∇u , ∇v)[L2(ω)]2 + (ku , v)L2 (Γ) + (∇u , ∇v)[L2 (Ω)]2 = f¯vdx ∀ v ∈ H01 (R). R



This, once again, is a well defined problem on R, and lim→0 u = u0 ∈ H01 (R). The limit u0 is harmonic on Ω, and u0 |ω solves the Robin problem (1.5). The equivalent estimate (1.7)

ku − u0 k1,R ' ku − u0 k1,ω + ku − u0 k0,Γ ' 

holds as long as ω is Lipschitz [23]. There is an exception to the equivalent estimate (1.7) in which u ≡ u0 . This happens if and only if u0 = 0 on Γ. This is the only case in which u0 ∈ H 2 (R). Otherwise, u0 only has the piecewise smoothness as for the Neumann problem. The finite element domain embedding method that determines uh ∈ Hh is a straightforward discretization of (1.6) on a uniform mesh. Using the very same argument as for Neumann problem, we obtain √ kuh − u0 k1,ω .  +h +  h. The weak formulation of the Dirichlet problem seeks u ∈ H01 (ω) such that Z (∇u, ∇v)[L2(ω)]2 = f vdx ∀ v ∈ H01 (ω). ω

If ω is simply connected, then on the continuous level, the domain embedding method determines u ∈ H01 (R) such that Z  −1  (1.8) (∇u , ∇v)[L2(ω)]2 +  (∇u , ∇v)[L2(Ω)]2 = f¯vdx ∀ v ∈ H01 (R). R



We have that lim→0 u = u ∈ that the limit is characterized by u0 |ω = u and u0 |Ω = 0, and that the equivalent estimate (1.9)

0

H01 (R),

ku − u0 k1,R ' ku k1,Ω ' 

6

SHENG ZHANG

holds, under the assumption that Ω is Lipschitz [23]. (This is true since we have assumed Γ ∈ C 2 and ω ⊂⊂ R; It is also true when, for example, a L-shaped polygon is tightly embedded in a rectangle; It, however, is not true when a circular domain is tightly embedded in a square.) Furthermore, if we let u1 ∈ H01 (R) be the unique function determined by ∆u1 = 0 in ω and Z 1 (1.10) (∇u , ∇v)[L2 (Ω)]2 = f vdx − (∇u0 , ∇v)[L2 (ω)]2 ∀ v ∈ H01(R). ω

Then we have ku − u0 −  u1 k1,R ' 2 .

(1.11)

An exception to (1.9) and (1.11) is that u ≡ u0 and u1 = 0. This occurs if and only if u also satisfies the homogeneous Neumann condition on Γ. This is equivalent to say u0 ∈ H 2 (R). Except for this special case, we do not have u0 ∈ H 2 (R), since the normal derivative of u0 is not continuous across Γ. The finite element domain embedding method determines uh ∈ Hh such that Z  −1  f¯vh dx ∀ v ∈ Hh . (1.12) (∇uh , ∇vh )[L2 (ω)]2 +  (∇uh , ∇vh )[L2 (Ω)]2 = R

We prove that (1.13) kuh − uk1,ω .  + 1inf

vh ∈Hh

 ku1 − vh1 k1,ω +


√ 1 ku − vh1 k1,Ω   1 0 0 0 + 0inf ku − vh k1,ω + √ kvh k1,Ω .  vh ∈Hh

0 problem, We construct a vh1 to interpolate u1 in the √ same√way as interpolating u in Neumann 2 and bound the above second term by  h + h . Our assumption of Γ ∈ C and f ∈ L2 (ω) assures that u = u0 |ω ∈ H 2 (ω) ∩ H01(ω). We construct an interpolation vh0 in such a way that it is zero on any open triangular √ element that is not entirely contained in ω. √This makes kvh0 k1,Ω = 0 and ku0 − vh0 k1,ω . h. The above third term is then bounded by h, which is a sharp estimate. We thus have √ √ √ (1.14) kuh − uk1,ω .  +  h + h  + h. √ So the accuracy of the√finite element domain embedding method is limited by O( h), which is assured when  . h. The rather low order of accuracy is a kind of locking phenomena that is not unusual in numerical computation of equations involving large parameters, like the Reissner–Mindlin plate and nearly incompressible elasticity [2]. The accuracy of the finite element domain embedding method for Dirichlet problem can be significantly improved by a slight adjustment of the nodes around Γ. One can move the near-by nodes onto Γ, and reconnects some of the affected nodes such that no mesh-line segment has one end in ω and the other in Ω. This way, a polygonal interpolation of Γ will be formed in the mesh, see the right figure in Figure 1. This can be done by using the B¨orger’s algorithm [6] that results in a triangulation that is quasi-uniform and shape regular, on which the discrete Laplacian is spectrally equivalent to that on the uniform mesh. On the adjusted triangulation, we then construct vh0 by simply interpolating u0 at all the nodes.

FINITE ELEMENT DOMAIN EMBEDDING METHODS

7

When ω is convex, we have ku0 − vh0 k1,ω . h, see [9], and kvh0 k1,Ω = 0 in (1.13). Thus to obtain the estimate √ √ (1.15) kuh − uk1,ω .  +  h + h  + h. So when  . h the error of the finite element domain embedding method in the H 1 (ω) norm is of the optimal order O(h). If ω is not convex, we still have ku0 − vh0 k1,ω . h. However, kvh0 k1,Ω is no longer zero. But it is bounded by O(h). The error estimate becomes (1.16)

√ √ h kuh − uk1,ω .  +  h + h  + h + √ . 

This estimate is also sharp. From this estimate we see that if we take  = Ch2/3 , the domain embedding finite element method has the error estimate O(h2/3 ) in H 1 (ω) norm. It is interesting to note that in this case a smaller  thwarts the convergence rate. This should be compared with Neumann problem and Dirichlet problem on convex domains, in which smaller  does not improve but not hurts the convergence rate of finite element domain embedding method. In conclusion, a slight adjustment of the finite element mesh to accommodate Γ appears which enhances the √ to be worthwhile for the Dirichlet √ problem, 2/3 convergence rate from h to h if ω is convex, and from h to h if ω is not convex. It is important to note that the discrete system (1.12) based on such adjusted mesh can be efficiently preconditioned by using the discrete Laplacian on the uniform mesh with the preconditioned iteration initiated with a vector of zero entries corresponding to nodes in R \ ω [6, 4]. If one wishes to achieve full accuracy for non-convex ω, then the mesh could be refined around the concave portion of Γ. A local mesh size O(h1.5 ) would lead to the full accuracy of order O(h) of the domain embedding finite element method in which  = Ch. This refinement is only needed around the concave portion of the boundary of ω. See [17] for such an example. If ω is polygonal and Γ aligns with the finite element mesh, then if u0 |ω ∈ H 2 (ω) (which is true when ω is a convex polygon, but is only an assumption otherwise), we construct vh0 and vh1 by nodal interpolation to u0 and u1 in (1.13) to obtain the estimate √ √ (1.17) kuh − uk1,ω .  +  h +  h + h. Here we have used the fact that u1 |ω ∈ H 1.5 (ω) and u1 |Ω ∈ H 1.5 (Ω) [15]. Thus when  . h one has a domain embedding finite element method that has the optimal accuracy of O(h). This was mentioned in [16] where a rigorous analysis was presented for Neumann problem under this alignment assumption. If ω is not simply connected, a term like −1 (u , v)L2 (Ω) needs to be added to the left hand side of (1.8). Then the estimate (1.9) holds and u0 |ω = u regardless of the connectivity of ω. Actually, in case that ω is multiply connected, adding a term of the form −1 (u , v)L2 (Ω0 ) in the left hand side of (1.8) is sufficient to guarantee the validity of the domain embedding method on the continuous level. Here Ω0 is the union of isolated components of Ω. The finite element domain embedding methods then should be based on such a modification of (1.8). Finally, we make some remarks on mixed boundary value problem. To make the presentation sufficiently general, we assume that Γ is split into three parts such that Γ = ΓD ∪ΓN ∪ΓR ,

8

SHENG ZHANG

and on ΓD , ΓN and ΓR , homogeneous Dirichlet, Neumann, and Robin conditions are imposed, respectively. Corresponding to this splitting, we divide Ω into three parts such that Ω = ΩD ∪ ΩN ∪ ΩR , and ∂ΩD ∩ Γ = ΓD , ∂ΩN ∩ Γ = ΓN , and ∂ΩR ∩ Γ = ΓR . The weak for1 mulation of the mixed problem determines u ∈ HD (ω) (that is the space of H 1 (ω) functions that vanish on ΓD ), such that Z 1 (1.18) (∇u, ∇v)[L2(ω)]2 + (ku, v)L2(ΓR ) = f vdx ∀ v ∈ HD (ω). ω

One the continuous level, the domain embedding method determines u ∈ H01 (R) such that (1.19) (∇u , ∇v)[L2 (ω)]2 + (ku , v)L2 (ΓR ) −1





+  (∇u , ∇v)[L2(ΩD )]2 + (∇u , ∇v)[L2(Ω\ΩD )]2 =

Z

f¯vdx ∀ v ∈ H01 (R).

R

When  → 0, u converges to a limit u0 ∈ H01 (R). The limit satisfies 1) u0 |ω solves (1.18) (assuming ω is simply connected), 2) u0 = 0 on ΩD , and 3) u0 is harmonic on ΩN and ΩR . Under the assumption that both ΩD and ω ∪ ΓD ∪ ΩD are Lipschitz domains, the equivalent estimate ku − u0 k1,R ' ku − u0 k1,ω + ku k1,ΩD ' 

(1.20)

holds with the exception that u ≡ u0 , which occurs when and only when u0 ∈ H 2 (R) [22]. The finite element domain embedding method then is again a straightforward discretization of (1.18) that determines uh ∈ Hh such that (1.21) (∇uh , ∇vh )[L2 (ω)]2 + (kuh , vh )L2 (ΓR ) −1

+

(∇uh , ∇vh )[L2 (ΩD )]2 0

+

(∇uh , ∇vh )[L2 (Ω\ΩD )]2 uh k1,ω

=

Z

R

¯ h dx ∀ vh ∈ Hh . fv

can be obtained under some regularity The convergence rate estimates on ku − assumption on the solution of the original problem (1.18), which is not guaranteed by our assumptions on Γ and f . For uniform meshes, ΓD generally cuts through some elements √ and the ends of ΓD lies in the interior of some elements, the rate is only bounded by h. By moving near-by nodes onto ΓD and its ends, and reconnecting some mesh lines, as for Dirichlet problem, we can improve the convergence rate to O(h) or O(h2/3 ), depending on the convexity of ΓD with respect to ω. As for Neumann and Robin problems, no treatment is needed to accommodate ΓN or ΓR . The paper is organized as follows. In Section 2, as a preparation, we derive some estimates on the finite element interpolation error on a strip around Γ that cuts through some elements. We also include a result about -dependent variational problem, which is a corner stone in estimating the regularization/penalty error. This result is also needed in the discretization error estimate for Dirichlet problem. In Section 3, we estimate the error of the finite element domain embedding methods for Neumann and Dirichlet problems, and determine the balance between the value of the regularization/penalty parameter and the finite element mesh-size. Since no new technique is needed for Robin problem and mixed problem, we shall not go into details to prove the results mentioned above.

FINITE ELEMENT DOMAIN EMBEDDING METHODS

9

2. Technical preliminaries We need to estimate the error of piecewise linear interpolation to functions whose restriction on both ω and Ω are in H 2 , but the functions themselves do not belong to H 2(R). Their normal derivatives are not continuous across the curve Γ, and Γ cuts through some of the triangular elements. For this purpose we need to estimate the H 1 norms on thin strips of piecewise H 2 functions and their finite element interpolations. The following lemmas will be used in the next section. Let y = γ(x) be a C 2 function on [0, L]. We consider a strip domain γδ that is bounded by x = 0, x = L, y = γ(x), and y = γ(x) + δ, see Figure 2. We assume that D is a given domain such that γδ ⊂⊂ D. We have the following estimates on functions restricted on γδ . δ

Figure 2. A strip γδ with a cluster of shape regular triangles. Lemma 2.1. Let w ∈ H 1 (D), then we have

kwk2L2(γδ ) ≤ Cδkwk21,D .

Here, C is a constant independent of δ and w. Proof. In terms of the curvilinear coordinates X = x, Y = y − γ(x), we have kwk20,γδ ' Z L Z δZ L 2 w 2 dX ≤ Ckwk21,D . w dXdY . For each Y ∈ (0, δ), by a trace theorem [1], we have 0

0

0

The constant C can be chosen to be independent of Y . Integrating this estimate with respect to Y from 0 to δ gives the desired estimate.  From this result, the following lemma easily follows. Lemma 2.2. Let w ∈ H 2 (D), then we have

kwk21,γδ ≤ Cδkwk22,D .

Here, C is independent of δ and w. Lemma 2.3. Let Tδ be the set of a cluster of shape regular, but not necessarily quasi-uniform, open triangles contained in γδ , and suppose that any point in γδ can only be covered by at most a certain number of triangles in Tδ . Let w ∈ H 2 (D). For each τ ∈ Tδ , we let Iτ w be the linear interpolation of w on the vertices of τ . Then we have X |Iτ w|21,τ ≤ Cδkwk22,D . τ ∈Tδ

Here, C is a constant independent of δ and w.

10

SHENG ZHANG

Proof. Let τ ∈ Tδ be one such triangle. We affine map it onto a standard reference triangle T of unit size through the mapping Fτ . The restriction on τ of a function w ∈ H 2 (D) is mapped to a function wˆ on T such that w◦F ˆ τ = w. Let I be the linear interpolation operator on the vertices of T . For any constant p0 , we have |Iτ w|21,τ = |Iτ (w − p0 )|21,τ ≤ C|I(wˆ − p0 )|21,T ≤ Ckw ˆ − p0 k2L∞ (T ) ≤ Ckw ˆ − p0 k22,T .

Thus, by the Bramble–Hilbert Lemma

|Iτ w|21,τ ≤ C(|w| ˆ 21,T + |w| ˆ 22,T ).

Scaling the right hand side back from T to τ , we get

|Iτ w|21,τ ≤ C(|w|21,τ + h2τ |w|22,τ ).

Here hτ is the diameter of τ , which does not exceed δ. Summing up the above inequality for all triangles in the cluster Tδ , we get X X X |Iτ w|21,τ ≤ C |w|21,τ + Cδ 2 |w|22,τ ≤ C|w|21,γδ + Cδ 2 |w|22,γδ . τ ∈Tδ

τ ∈Tδ

τ ∈Tδ

In the last step, we used the assumption that any point in γδ can only be covered by at most a certain number of triangles in Tδ . The desired result follows from Lemma 2.2.  We shall apply these estimates to strips that can be covered by a finite number of (translated and rotated) strips as defined here. Strips of width δ referred to in the later sections are understood in this sense. As we mentioned in the introduction, there are two sources of errors in the finite element domain embedding methods. The first one is due to the extension of the original boundary value problem from ω to R. The sharp estimates on this kind of errors can be obtained by using the following equivalent estimates on an abstract -dependent equation. This estimate will also be directly used in estimating the finite element interpolation error for Dirichlet problem in the next section. Let H, U, and V be Hilbert spaces, A : H → U a bounded linear operator, and B : H → V a bounded linear operator with closed range. We assume that (2.1) Thus the bilinear form

kAvkU + kBvkV ' kvkH ∀ v ∈ H.

(u, v)H := (Au, Av)U + (Bu, Bv)V defines an equivalent inner product on H. Furnished with this new inner product, the space H will be denoted by H. Let ker B ⊂ H be the kernel of the operator B. Let f be a linear continuous functional on H such that f |ker B = 0. There is a unique u ∈ H such that (2.2)

(Au , Av)U + (Bu , Bv)V = hf , vi ∀v ∈ H.

Without loss of generality, we assume that the operator B maps H onto V . (If necessary, we replace V by the range of B in it.) Then B is an isomorphism between (ker B)⊥ H (the orthogonal complement of ker B with respect to the H-norm) and V . Since we have assumed f |ker B = 0, according the the closed range theorem and Riesz representation theorem, there exists a unique u0 ∈ (ker B)⊥ H such that (2.3)

(Bu0 , Bv)V = hf , vi ∀ v ∈ H.

FINITE ELEMENT DOMAIN EMBEDDING METHODS

11

Lemma 2.4. Under the assumptions that f |ker B = 0 and B has closed range in V , as  → 0, u , the solution of (2.2), converges to the limit u0 ∈ (ker B)⊥ H defined by (2.3), and we have the equivalence estimate ku − u0 kH ' kB(u − u0 )kV '  kAu0 kU .

(2.4)

Therefore, if Au0 = 0 then u ≡ u0 . Otherwise u converges to u0 at the sharp rate of . Proof. From (2.2) and (2.3), we see that (2.5)

(A(u − u0 ), Av)U + (B(u − u0 ), Bv)V = − (Au0 , Av)U ∀v ∈ H.

Taking v = u − u0 , we get (2.6)

 kA(u − u0 )k2U + kB(u − u0 )k2V = − (Au0 , A(u − u0 ))U .

 0 It is easy to see that u lies in the subspace (ker B)⊥ H . So does u − u . Because B is an ⊥ isomorphism between (ker B)H and V , we have

kA(u − u0 )kU . ku − u0 kH ' kB(u − u0 )kV .

(2.7)

By the Cauchy–Schwarz inequality, there exists a constant C such that 1 | (Au0 , A(u − u0 ))U | ≤ C 2 kAu0 k2U + kB(u − u0 )k2V . 2 This estimate and the equation (2.6) show that kB(u − u0 )kV .  kAu0 k.

(2.8)

The equivalence (2.1) and (2.7) then lead to ku − u0 kH .  kAu0 kU . To see the lower bound, in (2.5) we take v = u0 to get (A(u − u0 ), Au0 )U + (B(u − u0 ), Bu0 )V = −  kAu0 k2U . If Au0 6= 0, then  kAu0 kU .  kA(u − u0 )kU +

kBu0 kV kB(u − u0 )kV . kB(u − u0 )kV . kAu0 kU

The equivalence (2.4) then follows.



3. Analysis of domain embedding methods In this section, we derive the error estimates for the finite element domain embedding methods for various boundary conditions, as described in the Introduction. The main issues appear in Neumann and Dirichlet boundary value problems. In Subsection 3.1, we analyze Neumann problem and assume that the grid on R is uniform. However, all the results hold for shape regular triangulations. We analyze the Dirichlet problem with uniform mesh in Subsection 3.2. The results of this subsection motivate mesh adjustment around Γ for Dirichlet problem, to which we devote Subsection 3.3.

12

SHENG ZHANG

3.1. Neumann boundary condition. We first consider the Neumann problem. Let ω ⊂ R2 be a bounded domain whose boundary Γ is in the C 2 class. Let f ∈ L2 (ω), we seek function u on ω solving the homogeneous Neumann problem whose weak formulation seeks u ∈ H 1 (ω)/R such that Z (3.1) (∇u, ∇v)[L2(ω)]2 = f vdx ∀ v ∈ H 1 (ω). ω

We assume that ω f (x)dx = 0, so that this is well-posed in the quotient space H 1 (ω)/R. Let R ⊂ R2 be a larger rectangular domain such that ω ⊂⊂ R. Let Ω = R \ ω be the fictitious domain. We extend the function f to a function f¯ on R by defining f¯ = 0 on Ω. The first step in the finite element domain embedding method is replacing (3.1) by the following problem defined on the rectangle R. For  > 0, one determines u ∈ H01 (R) such that Z   (3.2) (∇u , ∇v)[L2(ω)]2 + (∇u , ∇v)[L2(Ω)]2 = f¯vdx ∀ v ∈ H01 (R). R

R

This equation fits in the form of (2.2) in an obvious manner. The condition that B has closed range is assured as long as ω is a Lipschitz domain. All the other conditions of Lemma 2.4 are also satisfied, see [23] for details. We have lim→0 u = u0 . The limit u0 ∈ H01 (R) is characterized by the following two equations. Z 0 (3.3) (∇u , ∇v)[L2(ω)]2 = f vdx ∀ v ∈ H 1 (ω). ω

I.e., the restriction of u0 on ω solves the Neumann problem (3.1). And (3.4)

0

∆u = 0 in Ω,

∂u0 h − , 1iH −1/2 (Γ)×H 1/2 (Γ) = 0. ∂n

Here n− is the outward normal to Γ viewed as a boundary of Ω. By Lemma 2.4 we have (3.5)

ku − u0 k1,R ' |u − u0 |1,ω '  |u0 |1,Ω .

It is easy to see that |u0|1,Ω = 0 if and only if u0 = 0 on Γ, in which case u ≡ u0. This is the exceptional case in which u0 ∈ H 2 (R). Generally, u0 is only piecewise smooth with a normal derivative jump across Γ. We introduce a uniform triangulation of mesh size h on R, see the left figure in Figure 1, and let Hh ⊂ H01 (R) be the subspace of piecewise linear continuous functions subordinate to this triangulation. The finite element domain embedding method is a straightforward discretization of (1.2). It determines uh ∈ Hh such that Z   ¯ h dx ∀ v ∈ Hh . (3.6) (∇uh , ∇vh )[L2 (ω)]2 + (∇uh , ∇vh )[L2 (Ω)]2 = fv R

The numerical approximation of the solution of (3.1) offered by the finite element domain embedding method is the restriction of uh on ω. We have Theorem 3.1. Let u0 ∈ H01 (R) be the limit of u that is determined by the domain embedding equation (3.2). Then u0 |ω solves the Neumann problem (3.1). Let uh ∈ Hh be the solution

FINITE ELEMENT DOMAIN EMBEDDING METHODS

13

of the finite element domain embedding equation (3.6). We have the following estimate on the difference between uh |ω and u0 |ω .
√ (3.7) |uh − u0 |1,ω ≤ C  |u0 |1,Ω + C inf |u0 − vh |1,ω + |u0 − vh |1,Ω . vh ∈Hh

Here, C is a constant independent of h, , and f .

Proof. Since uh is the projection of u in the subspace Hh with respect to the inner product defined by the bilinear form in the left hand side of the domain embedding equation (3.2), we have √ √ |u − uh |1,ω + |u − uh |1,Ω ≤ |u − vh |1,ω + |u − vh |1,Ω ∀ vh ∈ Hh . Thus by (3.5)

|u0 − uh |1,ω ≤ ≤ ≤ ≤

|u0 − u |1,ω + |u − uh |1,ω√  |u0|1,Ω + |u − vh |1,ω + |u − vh |1,Ω √ √  0  |u0|1,Ω + |u0 − vh |1,ω + |u − u0 |1,ω + |u − v |u − u0 |1,Ω h |1,Ω + √ 0 3/2 0 0 (2  +  )|u |1,Ω + |u − vh |1,ω + |u − vh |1,Ω ∀ vh ∈ Hh .

The estimate thus follows.



This theorem is valid as long as ω is Lipschitz. Based on this result, we can establish the following error estimate on the finite element domain embedding method for Neumann problem under our assumptions that Γ ∈ C 2 and f ∈ L2 (ω). Under this assumption, by the regularity of elliptic differential equation [12], we have u0 |ω ∈ H 2 (ω) and ku0k2,ω . kf k0,ω . Since u0 |Ω is harmonic and it shares the value with u0 |ω on Γ, we have ku0 k2,Ω . ku0 k1.5,Γ . ku0 k2,ω . kf k0,ω . We can extend u0 |ω to R, see [20], to obtain a function u¯0 ∈ H 2 (R) such that u¯0|ω = u0 |ω and k¯ u0 k2,R . ku0 k2,ω . Similarly, we extend u0 |Ω to R to obtain a 0 2 function u ∈ H (R) such that u0 |Ω = u0 |Ω and ku0 k2,R . ku0 k2,Ω . Note that u¯0 |Ω 6= u0 |Ω , u0 |ω 6= u0 |ω , and u¯0 = u0 = u0 on Γ. Theorem 3.2. Under the assumption that Γ ∈ C 2 and f ∈ L2 (ω), there is a positive constant C that is independent of , h, and f , such that the following error estimate of the finite element domain embedding method for Neumann problems holds.  √  (3.8) |uh − u0 |1,ω ≤ C  +h +  h kf k0,ω . Before proving this theorem, we introduce some notations to describe triangles, vertices, and their relative relations with Γ. Let T be the set of all the open triangular elements of the finite element partition of R. We divide T into three distinctive subsets such that T = Tω ∪ TΓ ∪ TΩ . And TΓ = {τ ∈ T ; τ ∩ Γ 6= ∅}, Tω = {τ ∈ T ; τ ⊂ ω}, and TΩ = {τ ∈ T ; τ ⊂ Ω}.

Among the triangles of Tω , we let Tω0 be the subset of those that share a vertex with some triangle in TΓ . And define TΩ1 as a subset of TΩ similarly. See Figure 3 in which triangles in TΓ are shaded, those in Tω0 are marked by 0, and those in TΩ1 marked by 1. With a slight abuse of notations, we shall also use TΓ to denote the domain formed by the union of all triangles in TΓ , which is the interior of ∪τ ∈TΓ τ¯. And we use T , Tω , TΩ , Tω0 , and TΩ1 in the same way. Thus T = R, and ω ⊂ T ω ∪ TΓ , etc.

14

SHENG ZHANG ω 0

0 0

0 0

0 0

0

0 0

0

0 0

0 1

1 1

0 0

0

1 1

0

0 0

1

1

1

0

0

1 1

0

0

0

0

0

0

0

1 1

0

1

1 1

1

0

1

0

1

1 1

1

1

1

1 1

1

0

0

1 1

1 1

1 1

Ω

Figure 3. A portion of uniform mesh around Γ. Let V be the set of all the vertices in R of the triangulation. We divide V into three distinctive subsets such that V = Vω ∪ VΓ ∪ VΩ , in which VΓ = {ν ∈ V; ν is a vertex of a triangle of TΓ , or ν ∈ Γ}, Vω = {ν ∈ V; ν ∈ ω and ν is not a vertex in VΓ },

VΩ = {ν ∈ V; ν ∈ Ω and ν is not a vertex in VΓ }.

Note that TΓ could be empty, for example when Γ is straight and aligns with the triangulation, in which case VΓ is composed nodes on Γ. We let Vω0 = {ν ∈ Vω ; and ν is a vertex of a triangle in Tω0 },

VΩ1 = {ν ∈ VΩ ; and ν is a vertex of a triangle in TΩ1 }. Proof of Theorem 3.2. What we need to do is constructing √a finite element interpolation vh ∈ Hh to u0 such that |u0 − vh |1,ω . h and |u0 − vh |1,Ω . h. We define vh by requiring vh (ν) = u¯0 (ν) for ν ∈ Vω ∪ VΓ and vh (ν) = u0 (ν) for ν ∈ VΩ . Note that thus defined interpolation interpolates u0 on all the vertices except those in VΓ ∩ Ω. By the standard argument of finite element interpolation, we have |vh − u¯0 |21,τ ≤ Ch2 |¯ u0 |22,τ ∀ τ ∈ Tω ∪ TΓ .

Summing up these estimates, we get Thus

u0|22,R ≤ Ch2 |u0|22,ω ≤ Ch2 kf k20,ω . u0 |22,ω∪TΓ ≤ Ch2 |¯ |vh − u¯0 |21,ω∪TΓ ≤ Ch2 |¯

(3.9)

|vh − u0 |21,ω = |vh − u¯0 |21,ω ≤ Ch2 kf k20,ω .

By the standard argument of finite element interpolation, we have Summing up, we get (3.10)

|vh − u0 |21,τ ≤ Ch2 |u0|22,τ ∀ τ ∈ TΩ \ TΩ1 . |vh − u0 |21,T

1 Ω \T Ω

≤ Ch2 |u0|22,Ω ≤ Ch2 kf k20,ω .

The remaining part of Ω, i.e., TΩ1 ∪ (Ω ∩ TΓ ), is covered by a strip of width O(h), which meets the definition of γδ of Section 2 with δ = O(h). Thus, according to Lemma 2.2, we

FINITE ELEMENT DOMAIN EMBEDDING METHODS

15

have |u0|21,T 1 ∪(Ω∩TΓ ) ≤ Chku0 k22,R ≤ Chkf k20,ω .

(3.11)

Ω

We are left with estimating |vh |21,T 1 ∪(Ω∩TΓ ) . We consider this separately on Ω ∩ TΓ and Ω TΩ1 . We see that |vh |21,Ω∩TΓ ≤ |vh |21,TΓ . Since vh interpolates u¯0 at every vertex of triangles in TΓ , which itself is a cluster of shape regular triangles contained in a strip of with O(h), by Lemma 2.3, we have |vh |21,TΓ ≤ Chk¯ u0 k22,R . Thus |vh |21,Ω∩TΓ ≤ Chk¯ u0 k22,R ≤ Chkf k20,ω .

(3.12)

To estimate |vh |21,T 1 , we introduce an auxiliary interpolation ITΩ1 u0 that interpolates u0 on Ω every vertex of triangles in TΩ1 . Since TΩ1 is composed a cluster of shape regular triangles covered by a strip of width h, by Lemma 2.3, we have |ITΩ1 u0 |21,T 1 ≤ Chku0 k22,R ≤ Chkf k20,ω .

(3.13)

Ω

From the definition we see that vh (ν)−ITΩ1 u0 (ν) = 0 for all ν ∈ VΩ1 . And vh (ν)−ITΩ1 u0 (ν) = u¯0 (ν) − u0 (ν) for all ν ∈ VΓ ∩ Ω. Thus we have X X |vh − ITΩ1 u0 |21,T 1 ≤ C (vh (ν) − ITΩ1 u0 (ν))2 = C (¯ u0 (ν) − u0 (ν))2 . Ω

ν∈VΓ ∩Ω

ν∈VΓ ∩Ω

For each ν ∈ VΓ ∩ Ω, we construct a shape regular triangle τν with ν being one of its u0 − u0 ) be the linear interpolation vertices, and the other two vertices sitting on Γ. Let Iτν (¯ of u¯0 − u0 on τν . Since u¯0 = u0 on Γ, we have u0 − u0 )|21,τν ' (¯ |Iτν (¯ u0 (ν) − u0 (ν))2 . Thus X

ν∈VΓ ∩Ω

(¯ u0 (ν) − u0 (ν))2 ≤ C

X

ν∈VΓ ∩Ω

u0 − u0 )|21,τν . |Iτν (¯

The cluster of triangles {τν ; ν ∈ VΓ ∩ Ω} satisfies the condition of Lemma 2.3 with δ = h. We therefore have X u0 − u0 k22,R . u0 − u0 )|21,τν ≤ Chk¯ |Iτν (¯ ν∈VΓ ∩Ω

From these estimates, we see that |vh − ITΩ1 u0 |21,T 1 ≤ Chk¯ u0 − u0 k22,R ≤ Chkf k20,ω .

(3.14)

Ω

Therefore, by (3.13) and (3.14), u0 − u0 k22,R ) ≤ Chkf k20,ω . |vh |21,T 1 ≤ |ITΩ1 u0 |21,T 1 + |vh − ITΩ1 u0 |21,T 1 ≤ Ch(ku0 k22,R + k¯ Ω

Ω

Ω

Combining this and (3.12), we get (3.15)

|vh |21,T 1 ∪(Ω∩TΓ ) ≤ Chkf k20,ω . Ω

The theorem now follows from (3.9), (3.10), (3.11), and (3.15).



16

SHENG ZHANG

From this theorem we see that when Γ ∈ C 2 and f ∈ L2 (ω), and if we take  = Ch, we have a domain embedding finite element method on uniform mesh that achieves the full accuracy with an error of O(h). Although the triangulation was assumed to be uniform, it is clear that Theorem 3.2 holds for any triangulation that is shape regular, (but not necessarily quasi-uniform,) with a maximum mesh size h. The result of Theorem 3.1 is valid as long as ω is a Lipschitz domain, which does not guarantee the piecewise regularity of u0 . For example, when ω is a convex polygon, then u0 |ω ∈ H 2 (ω) but u0 |Ω can not be expected to have the H 2 regularity. If ω is a convex polygon and Γ aligns with the mesh, one can define vh by simply requiring it interpolate u0 at every node in V. Then one uses the regularity u0 |ω ∈ H 2 (ω) and u0 |Ω ∈ H 1.5(Ω) [15] to show the validity of the estimate in Theorem 3.2. If ω is a non-convex polygon, then u0 |ω ∈ H 2 (ω) may not hold, in which case, a local refinement may be needed at reentrant corners of ω to make |u0 − vh |1,ω . h. This can be done by taking a small rectangle around a corner, which conforms to the uniform mesh, and on which one uses a fine uniform mesh. Thus doing, FFT’s on the global uniform mesh and local uniform fine mesh can be combined in preconditioning the discrete system, in vine of additive multi-level Schwarz method or Chimera method [8]. For Robin boundary condition, the theory is similar to that of Neumann problem, so we will be very brief. We seek function u on ω solving the homogeneous Robin problem (3.16)

−∆u = f on ω,

∂u + ku = 0 on Γ. ∂n

Here k is a smooth, bounded, and positively valued function on Γ. The weak formulation is Z Z (∇u, ∇v)[L2(ω)]2 + kuvds = f vdx ∀ v ∈ H 1 (ω), (3.17) Γ ω 1 u ∈ H (ω). On the continuous level, the domain embedding method determines u ∈ H01 (R) such that Z Z   ¯ (3.18) (∇u , ∇v)[L2(ω)]2 + kuvds + (∇u , ∇v)[L2(Ω)]2 = fvdx ∀ v ∈ H01 (R). Γ

R

This equation also fits in the form of (2.2), in which the B operator has closed range in V if ω is Lipschitz. Therefore, as  → 0, u converges to a limit u0 ∈ H01 (R). The limit is harmonic on Ω and u0 |ω solves the Robin problem (3.17). By Lemma 2.4 we have (3.19)

ku − u0 k1,R ' |u − u0 |1,ω + ku − u0 k0,Γ '  |u0 |1,Ω .

The finite element domain embedding method is a straightforward discretization of (3.18) on a uniform mesh on R. The same results of Theorems 3.1 and 3.2 hold for the Robin problem, and the argument is exactly the same. The Robin problem provides an alternative approach to designing domain embedding for Dirichlet problem, other than the one we analyze below. In the equation (3.18), we can fix  and take k as a parameter. When k → ∞, the solution of (3.18) then converges to a function (at a rate much lower than 1/k [21]) whose restriction on ω solves the homogeneous Dirichlet problem. This method is essentially that of [3]. But the above theory of domain embedding method for Robin problem does not extend to this approach for Dirichlet problem.

FINITE ELEMENT DOMAIN EMBEDDING METHODS

17

3.2. Dirichlet boundary condition. We seek a function u on ω solving the homogeneous Dirichlet problem whose weak formulation is Z (∇u, ∇v)[L2(ω)]2 = f vdx ∀ v ∈ H01 (ω), (3.20) ω u ∈ H01 (ω). Under the assumption that ω is simply connected, on the continuous level, the domain embedding method determines u ∈ H01 (R) such that Z  −1  f¯vdx ∀ v ∈ H01 (R). (3.21) (∇u , ∇v)[L2(ω)]2 +  (∇u , ∇v)[L2(Ω)]2 = R

This is a penalty formulation. We extend u to a function u0 on R such that u0 |ω = u and u0 |Ω = 0. It is easy to see that for all v ∈ H01 (R)  Z 0  0  0 ¯ fvdx − (∇u , ∇v)[L2(ω)]2 . (3.22) (∇(u − u ), ∇v)[L2 (ω)]2 + (∇(u − u ), ∇v)[L2 (Ω)]2 =  R

This equation fits in the form of (2.2) by properly defining the operators and spaces. We see that ker B = {v ∈ H01 (R); v = 0 on Ω}

and, more importantly, the right hand side of (3.22), as a functional, annihilates this subspace. Thus there is a unique u1 ∈ H01 (R) such that ∆u1 = 0 on ω (since u1 ∈ (ker B)⊥ H) and Z  1 0 ¯ (3.23) (∇u , ∇v)[L2(Ω)]2 = f vdx − (∇u , ∇v)[L2 (ω)]2 ∀ v ∈ H01 (R). R

In strong form, this function satisfies (3.24)

−∆u1 = 0 in Ω,

∂u1 ∂u0 = on Γ, ∂n− ∂n

and u1 = 0 on ∂R.

Here n is the unit outward normal to Γ viewed as boundary of ω, and n− is the outward normal of Ω, which is opposite to n. When Ω is a Lipschitz domain, the operator B arising from fitting (3.22) into (2.2) has closed range. By Lemma 2.4 we have (3.25)

ku − u0 −  u1 k1,R ' |u − u0 −  u1 |1,Ω ' 2 |u1|1,ω .

From (3.24) we see that |u1 |1,ω = 0, which is equivalent to |u1 |1,Ω = 0, if and only if u0 satisfies the homogeneous Neumann condition on Γ, or equivalently, u0 ∈ H 2 (R). In this case u ≡ u0, otherwise, u converges to u0 at the sharp rate of . Remark. If ω is not simply connected, then Ω could have isolated components whose union is, say, Ω0 . In this case, the limit u0 of u does not satisfy the homogeneous Dirichlet condition on Γ, and the domain embedding method (3.21) fails. A remedy is adding, for example, the term −1 (u , v)L2 (Ω0 ) or −1 (u , v)L2 (Γ) (or something else as long as ker B can be identified with H01 (ω)) to the left hand side of the domain embedding equation (3.21). Then the method would be convergent and the convergence rate would be retained.

18

SHENG ZHANG

The finite element domain embedding method is a straightforward discretization of (3.21) in the finite element space Hh ⊂ H01 (R). In this subsection, as for the Neumann and Robin problems, we assume the triangulation is uniform. The numerical method determines uh ∈ Hh such that Z −1   (3.26) (∇uh , ∇vh )[L2 (ω)]2 +  (∇uh , ∇vh )[L2 (Ω)]2 = f¯vh dx ∀ v ∈ Hh . R

We have the following estimate on the error of this method. Theorem 3.3. Let u0 ∈ H01(R) be the limit of u determined by the domain embedding equation (3.21). Then u0 |ω solves the Dirichlet problem (3.20) and u0 |Ω = 0. Let uh ∈ Hh be the solution of the finite element domain embedding equation (3.26). We have
√ 1 √ |u − vh1 |1,ω + |u1 − vh1 |1,Ω (3.27) kuh − u0 k1,ω ≤ C  |u1|1,R + C  1inf vh ∈Hh   1 0 0 0 + C 0inf |u − vh |1,ω + √ |vh |1,Ω .  vh ∈Hh Here u1 is defined by (3.23) and C is a constant independent of h, , and f . Proof. Since uh is the projection of u in the subspace Hh with respect to the inner product defined by the bilinear form in (3.21), we have 1 1 |u − uh |1,ω + √ |u − uh |1,Ω ≤ |u − vh |1,ω + √ |u − vh |1,Ω ∀ vh ∈ Hh .   Thus, by using (3.25), we see that ku0 − uh k1,ω ≤ C|u0 − uh |1,R ≤ C (|u0 − u |1,R + |u − uh |1,R ) 1 ≤ C( |u1|1,R + |u − vh |1,ω + √ |u − vh |1,Ω )  1  0 1 ≤ C( |u |1,R + |u − u −  u |1,ω + |u0 +  u1 − vh |1,ω 1 + √ (|u − u0 −  u1 |1,Ω + |  u1 − vh |1,Ω ) ∀ vh ∈ Hh .  We write any vh ∈ Hh as vh = vh0 +  vh1 and use (3.25) again to obtain √ ku0 − uh k1,ω ≤ C[( + 1.5 + 2 )|u1|1,R +  |u1 − vh1 |1,ω + |u1 − vh1 |1,Ω

1 + |u0 − vh0 |1,ω + √ |vh0 |1,Ω ]. 

The desired estimate then follows.



This result is valid as long as Ω is a Lipschitz domain. In the following, we assume that ω ⊂⊂ R, Γ ∈ C 2 , and f ∈ L2 (ω). Under this assumption, we have the regularity that u0 |ω ∈ H 2 (ω) ∩ H01 (ω). Thus ∂u0 /∂n ∈ H 1/2 (Γ). This in turn means that ∂u1 /∂n− ∈ H 1/2 (Γ), cf., (3.24). Therefore, u1 |Ω ∈ H 2 (Ω). Since u1 |ω is harmonic and shares values with u1 |Ω on Γ, so u1 |ω ∈ H 2 (ω). Furthermore, we have the estimates ku0 k2,ω ≤ Ckf k0,ω , ku1k2,ω ≤ Ckf k0,ω , and ku1 k2,Ω ≤ Ckf k0,ω . We shall also use the fact that u0 |ω has a H 2 extension onto R [20], denoted by u¯0 , such that k¯ u0 k2,R ≤ Cku0 k2,ω ≤ Ckf k0,ω .

FINITE ELEMENT DOMAIN EMBEDDING METHODS

19

Theorem 3.4. Under assumption that ω ⊂⊂ R, Γ ∈ C 2 , and f ∈ L2 (ω), for Dirichlet problem the following error estimate of the finite element domain embedding method based on uniform mesh holds. √ √ √ ku0 − uh k1,ω ≤ C( +  h + h  + h)kf k0,ω .

Here, C is a constant independent of , h, and f .

Proof. We need to estimate the three terms in the right hand side of (3.27). The first one is trivial. For the other two, we need to construct finite element interpolation vh1 to u1 and vh0 to u0 , respectively. The construction of vh1 aims to estimate the infimum
√ 1 1 1 1 |u − v | + |u − v | . inf 1,ω 1,Ω h h 1 vh ∈Hh

We observe that this is in the same form as the second term in the right hand side of the estimate for Neumann problem, c.f., (3.7), except that the roles of Ω and ω are switched. We use exactly the same technique of proving Theorem 3.2 to estimate this term. We define vh1 ∈ Hh such that vh1 (ν) = u1 (ν) for all ν ∈ VΩ ∪ VΓ and vh1 (ν) = u1 (ν) for all ν ∈ Vω . Here u1 ∈ H 2 (R) is an extension of u1 |Ω to R such that u1 |Ω = u1 |Ω . With such defined interpolation, we have √ |u1 − vh1 |1,Ω ≤ Chkf k0,ω and |u1 − vh1 |1,ω ≤ C hkf k0,ω . Thus

(3.28)

inf 1

vh ∈Hh

√

√
|u1 − vh1 |1,ω + |u1 − vh1 |1,Ω ≤ C(  h + h)kf k0,ω .

We now turn to the last term in (3.27). I.e.,   1 0 0 0 (3.29) inf |u − vh |1,ω + √ |vh |1,Ω . 0 ∈H  vh h The best choice of the interpolation to estimate this one is taking vh0 ∈ Hh such that vh0 (ν) = u0 (ν) ∀ ν ∈ Vω ,

vh0 (ν) = 0 ∀ ν ∈ VΓ ∪ VΩ .

(See the remark below for the reason why this is the best.) This choice makes |vh0 |1,Ω = 0. To estimate the other term in (3.29), we observe that (3.30)

|u0 − vh0 |21,ω = |u0 − vh0 |21,Tω \Tω0 + |u0 − vh0 |21,Tω0 + |u0|21,ω∩TΓ .

Note that vh0 interpolates u0 at all the vertices of triangles in Tω \ Tω0 . And on each of the triangle τ ∈ Tω \ Tω0 , we have |u0 − vh0 |21,τ ≤ Ch2 |u0 |22,τ . Summing up, we get (3.31)

|u0 − vh0 |21,Tω \Tω0 ≤ Ch2 |u0|22,Tω \Tω0 ≤ Ch2 |u0 |22,ω .

Since Tω0 is covered by a strip of width O(h) that satisfies the condition of Lemma 2.2, so we have (3.32)

|u0|21,Tω0 ≤ Chk¯ u0 k22,R .

Here u¯0 is an H 2 extension of u0 |ω to R such that k¯ u0k2,R ≤ Cku0k2,ω . Note that a vertex 0 = u0 (ν) for all of any triangle in Tω0 either belongs to Vω0 or belongs to ω ∩ VΓ . Since P vh (ν) 0 0 2 0 0 ν ∈ Vω and vh (ν) = 0 for all ν ∈ ω ∩ VΓ , we have |vh |1,Tω0 ≤ C ν∈Vω0 (u (ν))2 . For each ν ∈ Vω0 , we construct a shape regular triangle τν with ν being a vertex, and the other two

20

SHENG ZHANG

vertices sitting on Γ. Let Iτν u¯0 be the linear interpolation of u¯0 on this triangle. We have (u0 (ν))2 ' |Iτν u¯0 |21,τν . Thus X |Iτν u¯0 |21,τν ≤ Chk¯ u0 k22,R . (3.33) |vh0 |21,Tω0 ≤ C ν∈Vω0

In the last step, we used Lemma 2.3. Finally, according to Lemma 2.2, we have (3.34)

|u0 |21,ω∩TΓ ≤ Chk¯ u0 k22,R .

The estimates (3.32), (3.33), and (3.34), together lead to |u0 − vh0 |21,Tω0 + |u0 |21,ω∩TΓ ≤ Chk¯ u0 k22,R . √ This and (3.31) show that |u0 − vh0 |1,ω ≤ C hk¯ u0k2,R , which together with (3.27) and (3.28) complete the proof.  √ Remark on sharpness of the estimate on (3.29). The estimated upper bound C h on (3.29) is rather big. This estimate was established by taking the interpolation vh0 to u0 . One of the reasons for the estimate to be so big is that the chosen vh0 does not offer any approximation to u0 on ω ∩ TΓ , where vh0 is identically equal to zero. Even so, the estimate is sharp. To see this, we consider an example of homogeneous Dirichlet boundary value problem in which ω is the unit circle, f = 1, and the solution is u0 |ω = (1 − x2 − y 2)/2. We embed ω in a larger rectangle R, on which we introduce a uniform triangulation of mesh-size h. The unit circle inevitably cuts √ through Ch triangles.0 The area of ω ∩ TΓ is of the magnitude h. And we have |u0 |1,ω∩TΓ ' h. The interpolation vh defined in the √ proof is zero on ω ∩ TΓ , which makes the √ error bound to be at least C h. But the estimate C h is sharp. One may try to construct a “better” piecewise linear approximation to u0 to make |u0 −vh0 |1,ω∩TΓ smaller. However, if the 0 |1,ω∩TΓ . vh0 were to offer any approximation to u0 , one would at least have that |vh0 |1,ω∩TΓ ' |u √ 0 0 0 0 ' |vh |1,ω∩TΓ ' |u |1,ω∩TΓ ' √ h. The Since vh is piecewise linear on TΓ , we have |v√ h |1,Ω∩T √Γ bound for (3.29) thus obtained would exceed h/ , which is much bigger than h. This argument shows that the estimate of Theorem 3.4 on the error of the finite element domain embedding method for Dirichlet problem is sharp. The optimal convergence rate √ √ O( h) can be achieved by taking the penalty parameter as  = O( h). Smaller  does not help in improving the accuracy but not hurt the√accuracy either. From a computational point of view, the best value of  should be  = C h. Since this accuracy is rather poor, it seems desirable to find ways to enhance it. (3.35)

3.3. Dirichlet boundary condition with adjusted mesh. One way to reduce the magnitude of (3.29) is adjusting the uniform mesh around Γ by moving the near-by nodes onto Γ and reconnecting some mesh lines to form a polygonal interpolation to Γ in the mesh. See Figure 1. The adjusted mesh must fulfill the following requirements: 1) The total number of nodes keeps unchanged. 2) No mesh-line segment has one end in ω and the other in Ω. 3) The adjusted triangulation is shape regular and quasi-uniform. Fulfilling these requirements ensures, on one hand, that the discrete Laplacian on the adjusted mesh is spectrally equivalent to the one on the uniform triangulation thus retains the effectiveness of preconditioning the discrete system by fast Poisson solvers. On the other hand, the magnitude of (3.29) is significantly reduced, and thus that the accuracy of the finite element domain embedding

FINITE ELEMENT DOMAIN EMBEDDING METHODS

21

method is accordingly enhanced. The B¨orgers’s algorithm [6] exactly carries out such adjustment. It is interesting to note that the algorithm was proposed in a context of solving Neumann problem by domain embedding methods [7]. Our above analysis shows that the mesh adjustment is actually not necessary for Neumann problems, while it has a significant effect for Dirichlet problem. It turns out that effect of this mesh adjustment is different for whether ω is convex. For convex ω, we prove that the full accuracy of O(h) is achieved by the domain embedding finite element method on such adjusted mesh if one takes  = O(h). Smaller  does not hurt the accuracy. If ω is not convex, then the accuracy can only be improved to O(h2/3 ) by taking  = Ch2/3 . But in this case, smaller  could diminish the improvement and sets the accuracy back to that of uniform mesh. This is a case in which the balance between  and h is delicate. We keep the notations of triangles and vertices introduced for uniform triangulation. For example, TΓ still comprises those triangles that intersect Γ, which are shaded in Figure 4. In addition, we shall use P to denote the set of all triangles enclosed by the polygonal interpolation of Γ. In consistence to our earlier convention, P also stands for the enclosed polygonal domain. ω

ω 0

0

0

0 0

0 h2

1

0 0

1

0

1 1

0 1

1

1 1

1 Ω

Convex ω

1 Ω

Non-convex ω

Figure 4. A portion of adjusted mesh around Γ. Theorem 3.5. We assume that ω ⊂⊂ R, Γ ∈ C 2 , and f ∈ L2 (ω). The triangulation T of R is obtained by adjusting a uniform mesh around Γ in the way described above. The finite element domain embedding method (3.26) is a straightforward discretization of the penalty formulation (3.21) on the adjusted mesh. Then if ω is convex, we have √ √ (3.36) ku0 − uh k1,ω ≤ C( +  h + h  + h)kf k0,ω .

If ω is not convex, we have

√ √ h ku0 − uh k1,ω ≤ C( +  h + h  + h + √ )kf k0,ω .  Here, C is a constant independent of , h, and f . (3.37)

Proof. The proof is based on the estimate (3.27) that was established in Theorem 3.3, which is still valid even when T is not uniform. There is nothing new in estimating the first two terms in the right hand side of (3.27). I.e., we have  |u1|1,R ≤ C  kf k0,R and

22

SHENG ZHANG

√ √ √  inf vh1 ∈Hh ( |u1 − vh1 |1,ω + |u1 − vh1 |1,Ω ) ≤ C( h + h )kf k0,R . We only need to analyze the third term   1 0 0 0 (3.38) inf |u − vh |1,ω + √ |vh |1,Ω . 0 ∈H  vh h √

We let vh0 ∈ Hh be such a piecewise linear function that interpolates u0 for all the vertices ν ∈ V. When ω is convex, this interpolation is consistent with that defined in the proof of Theorem 3.4. It is, however, different if ω is not convex, since the latter was required to be zero on the shaded TΓ , see Figure 4. When ω is convex, we have P ⊂ ω. And |u0 − vh0 |21,ω = |u0|21,ω\P + |u0 − vh0 |21,P .

Since ω \ P is covered by a strip γδ with δ = Ch2 , see the left figure in Figure 4, by u0 k22,R . Here u¯0 ∈ H 2 (R) is a H 2 extension of u0 |ω such Lemma 2.2, we have |u0 |21,ω\P ≤ Ch2 k¯ that k¯ u0 k2,R ≤ Cku0 k2,ω . Since vh0 interpolates u0 at the vertices of every triangle τ ∈ P, by the standard argument, we see |u0 − vh0 |21,P ≤ Ch2 |u0 |22,P . Therefore, |u0 − vh0 |1,ω ≤ Chkf k0,ω . Since vh0 = 0 on Ω, we get 1 |u0 − vh0 |1,ω + √ |vh0 |1,Ω ≤ Chkf k0,ω .  The estimate (3.36) then follows from this and Theorem 3.3. When ω is not convex, we have |u0 − vh0 |21,ω = |u0 |21,ω∩P c + |u0 − vh0 |21,ω∩P . Here P c = R \ P. Since ω ∩P c is contained in a strip of width O(h2 ), by Lemma 2.2, we have u0 −vh0 |21,P . Since vh0 interpolates u¯0 on u0 k22,R . We see that |u0 −vh0 |21,ω∩P ≤ |¯ |u0 |21,ω∩P c ≤ Ch2 k¯ the vertices of every triangle in P, by the standard argument of finite element interpolation, u0 k22,P . Therefore, we have |¯ u0 − vh0 |21,P ≤ Ch2 k¯ |u0 − vh0 |1,ω ≤ Chkf k0,ω .

(3.39)

Now, we consider |vh0 |1,Ω . Since vh0 is zero on Ω ∩ P c , we have |vh0 |21,Ω = |vh0 |21,Ω∩P . Let PΩ = {τ ∈ P; τ ∩ Ω 6= ∅}, which is a cluster of triangles contained in a strip of width O(h) Γ. On every τ ∈ PΩ , vh0 interpolates u¯0 . According to Lemma 2.3, we have P around 0 2 u0 k22,R . For each τ ∈ PΩ , the portion τ ∩ Ω is either contained in a τ ∈PΩ |vh |1,τ ≤ Chk¯ trapezoid of height O(h2 ) by one side of τ , or a smaller triangle of diameter O(h2 ) at a vertex of τ . In any case, we have |vh0 |21,τ ∩Ω ≤ Ch|vh0 |21,τ . We thus have X X |vh0 |21,Ω∩P = |vh0 |21,τ ∩Ω ≤ Ch |vh0 |21,τ ≤ Ch2 k¯ u0 k22,R . τ ∈PΩ

τ ∈PΩ

u0 k2,R . This together with (3.39) shows that Therefore, |vh0 |1,Ω ≤ Chk¯   1 0 h 0 0 √ |vh |1,Ω ≤ C(h + √ )k¯ inf |u − v | + u0 k2,R . 1,ω h 0 ∈H   vh h The proof of (3.37) is completed.



FINITE ELEMENT DOMAIN EMBEDDING METHODS

23

The estimates of this theorem are sharp. The sharpness argument is similar to that follows Theorem 3.4. From (3.36) we see that when ω is convex, then the finite element domain embedding method (3.26) with the adjusted mesh achieves the full accuracy O(h) in the H 1 norm by taking  = O(h). From a computational point view, the best value then is  = Ch. But smaller  wont hurt the accuracy. From (3.37) we see that if ω is not convex, then the full accuracy of O(h) can not be achieved by the mesh adjustment. In this case, there is an optimal value for . It is  = Ch2/3 , which leads to an overall accuracy of order h2/3 in the H 1 norm. Unlike Neumann problem with uniform mesh or Dirichlet problem on convex domain with adjusted mesh, for this case, a smaller  thwarts the accuracy of the finite element domain embedding. It could reduce the order back to h1/2 that is the order of accuracy achieved by the uniform mesh. References [1] R.A. Adams, Sobolev spaces, Academic Press, New York, 1975. [2] D.N. Arnold, R. Falk. A uniformly accurate finite element method for the Reissner–Mindlin plate, SIAM J. Numer. Anal., 26:1276-1290, 1989. [3] I. Babuska, The finite element method with penalty, Math. Comp., 27:221-228, 1973. [4] N.S.Bakhvalov, A.V.Knyazev, Preconditioned iterative methods in a subspace for linear algebraic equations with large jumps in the coefficients, in Domain decomposition methods in science and engineering, D. Keyes and J.Xu eds., Contemporary Mathematics, Vol 180, AMS, Providence, 1994. [5] P. Bochev, R.B. Lehoucq, On the finite element solution of the pure Neumann problem, SIAM Review, 47:50-66, 2005. [6] C. B¨orgers, A triangulation algorithm for fast elliptic solvers based on domain embedding, SIAM J. Numer. Anal., 27:1187-1196, 1990. [7] C. B¨orgers, O.B. Widlund, Finite element domain imbedding methods, SIAM J. Numer. Anal., 27:962977, 1990. [8] F. Brezzi, J.-L. Lions, O. Pironneau, The Chimera method for a model problem, in Numerical Mathematics and Advanced Applications -ENUMATH 2001, F. Brezzi, A. Buffa, S. Corsaro, A. Murli , Eds., 817-826, Springer, 2003. [9] P.G. Ciarlet, Finite element method for elliptic equations, North-Holland, 1978. [10] S. Del Pino, O. Pironneau, A fictitious domain based general PDE solver, in Numerical methods for scientific computing: variational problems and applications, Y. Kuznetsov, P. Neittanmaki, O. Pironneau eds., CIMNE, Barcelona, 2003. [11] K.O. Friedrichs, H.B. Keller, A finite difference scheme for generalized Neumann problems, in Numerical solutions of partial differential equations, J.H. Bramble ed., Academic Press, New York and London, 1966. [12] D. Gilbarg, N.S. Trudinger, Elliptic partial differential equations of second order, Springer-Verlag, Berlin, 1983. [13] R. Glowinski, T.W. Pan, Error estimate for fictitious domain/penalty/finite element methods, Calcolo, 29:125-141, 1992. [14] R. Glowinski, T.W. Pan, R. O. Wells Jr., X. Zhou, Wavelet and finite element solutions for the Neumann problem using fictitious domains, J. Comp. Phys., 126:40-51, 1996. [15] P. Grisvard, Singularities in boundary value problems, Masson, Springer-Verlag, 1992. [16] A. Knyazev, O. Widlund, Lavrentiev regularization + Ritz approximation = uniform finite element error estimates for differential equations with rough coefficients, Math. Comp., 72:17-40, 2003. [17] Yu. A. Kuznetsov, Domain decomposition and fictitious domain methods with distributed Lagrange multipliers, in Domain decomposition methods in science and engineering, N. Debit, M. Garbey, R. Hoppe, D. Keyes, Y. Kuznetsov, J. Periaux eds., CIMNC, Barcelona, 2002.

24

SHENG ZHANG

[18] R. Lazarov, S. Lu, S.V. Pereverzev, On the balancing principle for some problems of numerical analysis, Report 2005-25, RICAM, Austrian Academy of Sciences, 2005. [19] S. Pereverzev, E. Schock, On the adaptive selection of the parameter in regularization of ill-posed problems, SIAM J. Numer. Anal., 43:2060-2076, 2005. [20] E.M. Stein, Singular integrals and differentiability properties of functions, Princeton University Press, 1970. [21] S. Zhang, Equivalence estimates for a class of singular perturbation problems, C. R. Acad. Sci. Paris, Ser. I, 342:285-288, 2006. [22] S. Zhang, A domain embedding method for mixed boundary value problems, C. R. Acad. Sci. Paris, Ser. I, 343:287-290, 2006. [23] S. Zhang, Sharp convergence rate of domain embedding methods for various boundary conditions, submitted, available at www.math.wayne.edu/~sheng/embedding.pdf.