Obtaining higher-order Galerkin accuracy when the boundary is polygonally approximated

1 Introduction

When a Dirichlet problem on a smooth domain is approximated by a polygon, an error occurs that is suboptimal for piecewise quadratic approximations [1], [2], [3], [4]. There are several approaches that have been developed to overcome the sub-optimality. Two common approaches are: (1) methods that use curved elements (see, for example [5], [6], [7], to name a few) and (2) methods that use augmented polygonal approximations. The method we develop here falls within the second category. A non-exhaustive list of such methods includes [8], [9], [10], [11], [12], [13], [14], [15], [16], [17], [18], [19].

To put our contribution in context, we consider a model problem and previous numerical methods for this problem. Let Ω be a piecewise smooth, bounded, two-dimensional domain that is on only one side of its boundary. Consider the Poisson equation with Dirichlet boundary conditions:

We assume that f and g are sufficiently smooth that u can be extended to be in W 2 , q ( Ω ̂ ) , where Ω ̂ contains a neighborhood of the closure of Ω , for some q in the range 1 < q ⩽ ∞ . For a Lipschitz domain, this is possible with q < 4 / 3 .

One way to discretize (1) is to approximate the domain Ω by polygons Ω h , where the edge lengths of ∂ Ω h are of order h in size. Then conventional finite elements can be employed on each Ω h , with the Dirichlet boundary conditions being approximated by the assumption that u h = g ̂ on ∂ Ω h [20], with g ̂ appropriately defined. For example, let us suppose for the moment that g ≡ 0 and we take g ̂ ≡ 0 as well. In particular, we assume that Ω h is triangulated with a nondegenerate mesh T h of maximum triangle size h , and the boundary vertices of Ω h are in ∂ Ω . We define W ̊ h k ≔ H 0 1 ( Ω h ) ∩ W h k , where

and P k is the set of polynomials of total degree ⩽ k . Then the standard finite element approximation finds u h ∈ W ̊ h k satisfying

where a h ( u , v ) ≔ ∫ Ω h ∇ u ⋅ ∇ v d x . Here we assume that f is extended smoothly outside of Ω .

This approach for k = 1 (piecewise linear approximation) leads to the error estimate

However, when this approach is applied with piecewise quadratic polynomials ( k = 2 ) , the best possible error estimate [20] is

which is less than optimal order by a factor of h . The reason of course is that we have made only a piecewise linear approximation of ∂ Ω . Table 1 summarizes some computational experiments for the test problem in Section 3. We see a significant improvement for quadratics over linears, but there is almost no improvement with cubics. Moreover, we will see that a significant improvement using quadratics can be obtained using simple approaches that modify the variational form.

The first method using a polygonal approximate domain to obtain higher order schemes was based on modifying the method of Nitsche [5] and appeared in 1972 [8]. There are two main ingredients in the scheme in ref. [8]:

1. Nitsche’s penalization to both enforce the boundary condition and stabilize the method;

2. Taylor’s expansions to give approximate boundary conditions on Ω h .

The idea is to think of u as the solution of a PDE on the approximate domain, with boundary conditions that need to be approximated.

Methods for any polynomial order are developed in ref. [8], and it is proved that the error estimates are optimal in the energy norm. The bulk of ref. [8] analyzes a non-symmetric extension of Nitsche’s method, but a symmetrized version is defined [8], eq.  (3.12)] for the lowest-order scheme for which the same estimates hold by trivial modifications, and estimates are described in the case that the distance to the boundary is only approximated.

Two years later [8], was extended [15] to cover both 2 and 3 dimensions, variable coefficients, and estimates in H − 1 as well as H 1 and L 2 . The paper [15] uses the symmetric form throughout, it allows both linear and nonlinear approximations of the boundary (that is, the approximating domains can be curved), and it also applies the elliptic theory to nonlinear parabolic equations.

The work of refs. [8], [15] has been extended in many ways [9], [10], [11], [16], [21], [22], [23], [24] over the intervening decades.

Recently, the related shifted-boundary method has been developed [17], [18] and analyzed [25], [26]. It uses the two main ingredients (i.e., Nitsche’s penalization and Taylor’s expansions) to develop methods for Poisson’s problem, fluid flow problems and advection–diffusion problems. The method seems to be stable when the boundary of Ω h is within O ( h ) of the boundary of Ω . Cheung et al. [12] develop a method using average Taylor expansions, and they prove stability in the case the boundary of Ω h is within O ( h 3 / 2 ) of the boundary of Ω . Finally, a series of papers [13], [14], [19], develops hybridizable discontinuous Galerkin (HDG) methods using polygonal approximations. HDG has its own stabilization mechanism which they use. A distinctive feature is that they approximate both u and ∇ u independently. The latter has been named the “transfer path method” (TPM) [27] and been developed not only for HDG methods, but also for mixed finite element methods [28], [29]. Therefore, it does not require a stabilization parameter.

In this paper, we introduce a new method that uses a first order Taylor expansion which naturally leads to a Robin-type boundary problem on Ω h . However, we do not use Nitsche’s stabilization, and thus we avoid the need to choose a penalty parameter. Another parameter-free method is developed in refs. [30], [31], [32].

The Robin-type method studied here is naturally symmetric, but it can be not positive definite. Despite this, we are able to prove optimal estimates for piecewise quadratics and cubics. The analysis is quite different from what is used for methods that are rooted in Nitsche’s method, such as refs. [8], [15] and its descendents. It is not a simple perturbation of the standard arguments. For this reason, we restrict attention to the model problem (1) in order to make our analysis as transparent as possible. We assume for simplicity that the vertices of the boundary Ω h belong to boundary of Ω (i.e., a fitted mesh). We give numerical results that show that our estimates are sharp.

2 The Bramble–Dupont–Thomée approach

We start by reviewing the method [8] of Bramble–Dupont–Thomée (BDT) which achieves high-order accuracy by modifying Nitsche’s method [5] applied on Ω h . We assume that Ω h ⊂ Ω and we do not necessarily assume that the boundary vertices of Ω h belong to ∂ Ω . The bilinear form used in ref. [8] is

Here, n denotes the outward-directed normal to ∂ Ω h and

Contrast the definition of δ to the closely related function d defined by

see Figure 1.

Figure 1:

Definitions of (A) δ and (B) d .

For simplicity the assume that g = 0 . Then the BDT method will find u h ∈ W h k such that

If δ were 0, this would be Nitsche’s method on Ω h .

Corrections of arbitrary order, involving terms δ ℓ ∂ ℓ u ∂ n ℓ for ℓ > 1 are studied in ref. [8], but for simplicity we restrict attention to the first-order correction to Nitsche’s method given in (5). The error estimates obtained in ref. [8] are as follows

where

Thus using the variational form (5) leads to an approximation that is optimal-order with quadratics and cubics and is only suboptimal for quartics by a factor of h .

3 An example: the circle

We consider a numerical example. Consider the case where Ω is a disc of radius R centered at the origin, in which case we have d ( x ) = R − | x | . However, it is more difficult to evaluate δ ( x ) . We assume that the vertices of ∂ Ω h lie on ∂ Ω . We have x + δ ( x ) n ∈ ∂ Ω for x ∈ ∂ Ω h , where n denotes the outward normal to Ω h . Let t denote the unit tangent vector to Ω h such that n × t points up. We can write x = ( x ⋅ n ) n + ( x ⋅ t ) t , and ( x ⋅ t ) 2 = | x | 2 − ( x ⋅ n ) 2 . Since | x + δ ( x ) n | = R , we have

Then

Note that for x ∈ ∂ Ω h , | x | ⩽ R and x ⋅ n > 0 . Since δ ( x ) ⩾ 0 , we must pick the plus sign, so

Note that this formula does not depend on the placement of the boundary vertices. If x i is a boundary vertex, then | x i | = R , and δ ( x i ) = 0 .

It is not hard to see that d − δ = O ( h 4 ) in this case.

This problem is simple to implement using the FEniCS Project code dolfin [33]. We take R = 1 , u ( x , y ) = 1 − ( x 2 + y 2 ) 3 , and f = 36 ( x 2 + y 2 ) 2 in the computational experiments described subsequently. Computational results for this example are given in Table 2, where we see optimal order approximation for k ⩽ 3 , improvement for k = 4 over k = 3 (suboptimal by a factor h − 1 / 2 ), and no improvement for quintics. These errors are depicted in Figure 2.

Figure 2:

Errors u h − u I in (A) L 2 ( Ω h ) and (B) H 1 ( Ω h ) as a function of the maximum mesh size for the BDT method with γ = 100 . The asterisks indicate data for (A) k = 4 and (B) k = 5 . The interpolant u I is defined in (16).

4 A new method based on a Robin-type approach

One issue with the BDT method is that the resulting linear system is not symmetric, although it is easily symmetrized as we discuss in Section 11. Here we develop a technique that leads directly to a symmetric system. More importantly, this method does not require the parameter(s) from Nitsche’s method. For Nitsche’s method to succeed, γ must be chosen appropriately [34]. Naturally, there is a price to pay for having a parameter-free method. We will have to make some restrictive assumptions about the domain boundary that would not be required for the success of BDT.

We first separate ∂ Ω into its piecewise linear part and its curvilinear part. We will assume that

where Γ 0 is a finite union of piecewise linear segments and the S i ’s are C 2 and nowhere linear. The rectilinear part Γ 0 of the boundary may be empty, as is the cases for our numerical examples. Denote by

the intersection points (if any) where S i and S i + 1 meet. We make the following assumption on the domain.

In Figure 3, we show a domain with two S i ’s and exactly one point y i , where the two circles meet tangentially, the point (2, 1). The domain may be written as

where D ( x , y ) denotes the unit disc centered at ( x , y ) . The dashed lines indicate this construction via constructive solid geometry.

Figure 3:

A P-shaped domain (solid lines) with one point y i . The dashed lines indicate the definition (9) via constructive solid geometry.

It may be that different S i are disconnected and have no end points, as in our example in Section 9.2. So the set of y i ’s could be empty ( ϰ = 0 ) .

For the method in this section we assume that the vertices of Ω h belong to ∂ Ω , and hence Ω h might not be a subdomain of Ω . Thus, we need to define δ in this case. We assume that for every x ∈ ∂ Ω h there is a unique smallest number δ ( x ) in absolute value such that

For x ∈ Γ 0 , we have δ ( x ) = 0 provided that Γ 0 ⊂ ∂ Ω h .

We assume that the approximate domain boundary ∂ Ω h can be decomposed into three parts, as follows. Let E h be the edges of ∂ Ω h . Define

where e o denotes the interior of e . Let Γ = Γ + ∪ Γ − . We make the following assumption on the boundary approximations.

This assumption means that δ cannot change sign in the interior of an edge.

Note that Γ depends on h , but we omit a subscript to simplify notation.

Our method resembles a Robin-type of boundary condition on Γ . It is similar to the closely related problem:

Here we define

for x ∈ Γ . The key here is that, using that u = g on ∂ Ω , for x ∈ Γ ( x not a vertex of ∂ Ω h ) we have

Here E u denotes an extension of u outside of Ω to allow for the possibility that Ω h ⊄ Ω . The function S defined by

is piecewise smooth for smooth u and will play a significant role in the study of quadrature error in Section 6.3. We postpone to Section 13 a discussion of the smoothness of S . But suffice it to say that

and S ̂ is piecewise smooth.

Now we can define the method. Assume that Ω h is triangulated by a nondegenerate mesh, and we denote by h Γ the maximum edge length on the boundary and by h Ω the maximum mesh size over the entire domain. By nondegenerate, we mean the following. For each triangle T in a mesh T h , let ρ min T be the radius of the largest circle lying inside T and let ρ max T be the radius of the smallest circle containing T . Define

A nondegenerate mesh family is one for which ρ T h ⩾ ρ > 0 for each mesh in the family. In the following, many important constants will depend on ρ T h , but we will suppress mentioning the dependence to simplify the discussion.

Let { ψ j } be the usual Lagrange nodal basis for W h k , and define the interpolant

for any continuous function u on the closure of Ω h .

We start by defining one finite element space we will use:

where W h k is defined in (2). Also define

where g I ∈ C ( ∂ Ω h ) is a suitable approximation of g and is a piecewise polynomial of degree at most k on ∂ Ω h . For example, we can take g I to be the interpolant, defined in (16), of g ̂ defined in (11). However, note that the definition (18) does not require values of g I except on Γ 0 and at vertices of ∂ Ω h where we do know these values.

The bilinear form for the new method is given by

where

The form c ̃ h ( ⋅ , ⋅ ) depends on h because Γ depends on h . Then the method solves: Find u h ∈ V ̊ h k ( g ) such that, for all v ∈ V ̊ h k ,

where E f is an extension of f outside of Ω , not necessarily the same extension used for u . For the moment, we assume that we can compute c ̃ h ( g ̂ , v ) exactly, based on the formula (11). But this is not a practical method because it requires us to work with a space V ̊ h k whose functions are forced to vanish at prescribed points. However, it is easier to analyze this limited method and from this we learn the key issues for more complicated (and practical) versions.

Integration by parts gives

Define the extension error

For the two-circle problem in Section 9.2, we have analytical expressions for u and f , and using these for the extensions, we get − Δ E u − E f = 0 , so the term E can typically be ignored. Combining (12), (21), and (22), we get

We can estimate the second term via

where Ω Δ Ω h = Ω \ Ω h ∪ Ω h \ Ω . Here, the normal direction is the normal to Ω h . We postpone the proof of (25) to Section 12. Our method thus appears at first to be a standard variational crime with small error. However, the stability of the method is not standard.

The role of the form c ̃ h ( ⋅ , ⋅ ) appears to be of two parts. Due to the singularity, it forces adoption of the boundary conditions. But in addition, it forms the required correction in view of (25).

4.2 The linking lemma

Due to the singularity of δ − 1 , c ̃ h ( u , v ) may not be well defined for all u , v ∈ V ̊ h k . Therefore, we first show that this is not the case if we make Assumption 1.

Lemma 1.

Under Assumption 1,

(35) | c ̃ h ( u , v ) | < ∞ ∀ u , v ∈ V ̊ h k .

Proof.

Assumption 1 implies that δ is non-zero in the interior of each edge. If δ ′ ≠ 0 at the end of each edge, then δ − 1 v is bounded for all v ∈ V h k . The condition δ ′ ≠ 0 is obvious for edges whose vertices are in the interior of each segment S i , since the curvature of S i in the interior would be bounded away from zero, but at a boundary, the segment could be flat. But examining Figure 4 shows that δ ′ ≠ 0 there as well. The tangent to S i at the left vertex has zero slope, but the chord connecting that vertex and the next has a positive slope, due to the strict monotonicity of the boundary arc. The difference of slopes is δ ′ . □

Figure 4:

Bound for δ for the example in (38).

Assumption 1 does not hold for a domain of the form

(36) Ω f = ( x , y ) : | x | < 1 , f ( x ) < y < 1

if f is defined by

(37) f ( x ) = 0 , x ⩽ 0 , sin ( π / x ) e − 1 / x , x > 0 .

In particular, if we take boundary mesh points at (0,0) and ( 1 / n , 0 ) , then the corresponding edge e with these vertices in the polygonal approximation will be on the x -axis. Thus δ = f in this interval and c h ( v , v ) < ∞ only if v vanishes in that interval. Although we must rule out such domains for the new method, the BDT method would not be deterred by such boundaries.

On the other hand, if f is monotone, Assumption 1 does hold for Ω f . For f defined by

(38) f ( x ) = 0 , x ⩽ 0 , e − 1 / x , x > 0 ,

the slope of δ in [ 0,1 / n ] is e − n at the origin. Therefore Assumption 1 holds for reasonable domains.

Lemma 2.

Under Assumption 1, there exists a constant c 1 > 0 such that

(39) ‖ v ‖ a ⩽ c 1 h | v | c ∀ v ∈ B h k .

Proof.

Let E h Γ be the collection of edges that are a subset of Γ and let T h Γ be triangles T such that T has an edge in E h Γ . Then, if v ∈ B h k and using inverse estimates [20] we have

(40) ‖ v ‖ a 2 = ∑ T ∈ T h Γ ‖ ∇ v ‖ L 2 ( T ) 2 ⩽ ∑ T ∈ T h Γ C h T 2 ‖ v ‖ L 2 ( T ) 2 ⩽ ∑ e ∈ E h Γ C h e ‖ v ‖ L 2 ( e ) 2 ⩽ ∑ e ∈ E h Γ C h e max x ∈ e | δ ( x ) | ‖ | δ | − 1 / 2 v ‖ L 2 ( e ) 2 .

The result is complete after we use that

(41) max x ∈ e | δ ( x ) | ⩽ C h e 2

for all e ∈ E h Γ , which holds since δ is essentially the error in a piecewise linear approximation of the boundary. □