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

We study two techniques for correcting the geometrical error associated with domain approximation by a polygon. The first was introduced some time ago <cit.> and leads to a nonsymmetric formulation for Poisson's equation. We introduce a new technique that yields a symmetric formulation and has similar performance. We compare both methods on a simple test problem.

## Authors

• 1 publication
• 17 publications
• 1 publication
05/31/2020

### Equal higher order analysis of an unfitted discontinuous Galerkin method for Stokes flow systems

In this work, we analyze an unfitted discontinuous Galerkin discretizati...
03/16/2021

### Isogemetric Analysis and Symmetric Galerkin BEM: a 2D numerical study

Isogeometric approach applied to Boundary Element Methods is an emerging...
03/09/2020

### Variational Time Discretizations of Higher Order and Higher Regularity

We consider a family of variational time discretizations that are genera...
12/16/2019

### Higher order Galerkin-collocation time discretization with Nitsche's method for the Navier-Stokes equations

We propose and study numerically the implicit approximation in time of t...
03/08/2021

### A new Besse-type relaxation scheme for the numerical approximation of the Schrödinger-Poisson system

We introduce a new second order in time Besse-type relaxation scheme for...
07/03/2021

### Corrected trapezoidal rules for singular implicit boundary integrals

We present new higher-order quadratures for a family of boundary integra...
03/04/2022

### Constructing Nitsche's method for variational problems

Nitsche's method is a well-established approach for weak enforcement of ...
##### This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

## 1. Introduction

When a Dirichlet problem on a smooth domain is approximated by a polygon, an error occurs that is suboptimal for quadratic approximation [1, 10, 11]. However, this can be corrected by a modification of the variational form [2]. Here we review this approach and suggest a new one.

Let be a smooth, bounded, two-dimensional domain. Consider the Poisson equation with Dirichlet boundary conditions:

 (1) −Δu=f in Ω,u=g on ∂Ω.

We assume that and are sufficiently smooth that can be extended to be in , where contains a neighborhood of the closure of .

One way to discretize (1) is to approximate the domain by polygons , where the edge lengths of are of order in size. Then conventional finite elements can be employed, with the Dirichlet boundary conditions being approximated by the assumption that on [3], with

appropriately defined. For example, let us suppose for the moment that

and we take as well. In particular, we assume that is triangulated with a quasi-uniform mesh of maximum triangle size , and the boundary vertices of are in . We define where

 Wkh={v∈C(Ωh):v|T∈Pk(T),∀T∈Th}.

Then the standard finite element approximation finds satisfying

 (2) ah(uh,v)=(f,v)L2(Ωh),∀v∈˚Wkh,

where . Here we assume that is extended smoothly outside of .

This approach for

(piecewise linear approximation) leads to the error estimate

 ∥u−uh∥H1(Ωh)≤Ch∥u∥H2(^Ω).

However, when this approach is applied with piecewise quadratic polynomials (), the best possible error estimate is

 (3) ∥u−uh∥H1(Ωh)≤Ch3/2,

which is less than optimal order by a factor of . 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 2.1. 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.

There have been many techniques introduced to circumvent the loss of accuracy with quadratics (and higher-order piecewise polynomials) [11, 6]. However, all of them require some modification of the quadrature for the elements at the boundary.

Here we review an approach by Bramble et al. [2] that solves directly on , but with a modified variational form based on the method of Nitsche [6]. The method [2] has been modified and applied in many ways [4]. However, the method in [2] leads to a non symmetric bilinear form. Given this shortcoming we define a new method that is symmetric and solves the problem on that has similar convergence results. As we will see in the next section, one main idea in [2] is that one uses a Taylor series of the solution near the boundary to define appropriate boundary conditions on . We should mention that this idea has been used recently (see for example [5, 8]).

## 2. The Bramble-Dupont-Thomée approach

The method [2] of Bramble-Dupont-Thomée (BDT) achieves high-order accuracy by modifying Nitsche’s method [6] applied on . We assume that and we do not necessarily assume that the boundary vertices of belong to . The bilinear form used in [2] is

 (4) Nh(u,v)=ah(u,v)−∫∂Ωh∂u∂nvds−∫∂Ωh(u+δ∂u∂n)(∂v∂n−γh−1v)ds

Here, denotes the outward-directed normal to and

 δ(x)=min{s>0:x+sn∈∂Ω}.

Contrast the definition of to the closely related function defined by

 d(x)=min{|x−y|:y∈∂Ω}.

For simplicity the assume that . Then the BDT method will find such that

 Nh(uh,v)=∫Ωhfvdx for all v∈Wkh.

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

Corrections of arbitrary order, involving terms for are studied in [2], but for simplicity we restrict attention to the first-order correction to Nitsche’s method given in (4). The error estimates obtained in [2] are as follows

 |||u−uh|||1≤Chk∥u∥Hk+1(Ω)+Ch7/2∥u∥W2∞(Ω),

where

 |||v|||1:=(ah(v,v)+h−1∫∂Ωhv2ds+h∫∂Ωh(∂v∂n)2ds)1/2.

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

### 2.1. An example of a circle

We consider a numerical example. Consider the case where is a disc of radius centered at the origin, in which case we have . However, it is more difficult to evaluate . We have for , where denotes the outward normal to . We can write , and . Since , we have

 R2=(x⋅t)2+(x⋅n+δ(x))2=|x|2−(x⋅n)2+((x⋅n+δ(x))2.

Then

 δ(x)=±√R2−|x|2+(x⋅n)2−x⋅n.

Note that for , and . Since , we must pick the plus sign, so

 δ(x)=√R2−|x|2+(x⋅n)2−x⋅n.

It is not hard to see that in this case.

This problem is simple to implement using the FEniCS Project code dolfin [7]. We take , , and in the computational experiments described subsequently. Computational results for this example are given in Table 2 where we see optimal order approximation for , improvement for over (suboptimal by a factor ), and no improvement for quintics. These errors are depicted in Figure 2.

## 3. 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 possible to symmetrize the method as we discuss in Section 8. Here we develop a technique that leads to a symmetric system. Moreover, this method does not require the parameter(s) from Nitsche’s method. For Nitsche’s method to succeed, must be chosen appropriately [9].

We first separate to its piecewise linear part and its curvilinear part. We will assume that where is a piecewise linear segment and are and no where linear. We let the end points of to be .

For the method in this section we assume that the vertices of belong to and hence might not be a subdomain of . Hence, we need to define in this case. We assume that for every that is there is a unique smallest number in absolute value such that

 x+δ(x)n(x)∈∂Ω.

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

 (5) Γ±=∪{e∈Eh:±δ|eo>0},

where denotes the interior of . Let . We assume the following.

###### Assumption 1.

We assume that all the vertices of belong to and that each (for ) is a vertex of . Finally, we assume that

 ∂Ωh=Γ0∪Γ.

Our method is based on a Robin type of boundary condition on . In fact, our method will be based on the closely related problem:

 −Δw= f, on Ω, w= 0, on Γ0, w+δ∂w∂n= ^g, on Γ.

Here we define for and not a vertex of . The key here is that, using that vanishes on , for ( not a vertex of ) we have

 (6) u(x)+δ∂u∂n(x)=^g(x)−δ22∂nnu(z),

where lies in the line segment with end points and .

Now we can write the method. We start by defining the finite element space we will use

 Vkh={Wkh:v=0 on Γ0,v(x)=0 for all vertices % of x of ∂Ωh}.

Also define

 Vkh(g)={Wkh:v=gI on Γ0,v(x)=Ig(x) for all% vertices of x of ∂Ωh}.

where is a suitable approximation of and is a piecewise polynomial of degree at most on .

The bilinear form is given by

 bh(u,v):=ah(u,v)+ch(u,v),

where

 ch(u,v)=∫Γδ−1uvds.

Then the method solves:

Find such that

 (7) bh(uh,v)=∫ΩhFv+∫Γδ−1^gvds. for all v∈Vkh.

Here

 F={fonΩ∩ΩhI1fonΩh∖Ω,

where

is the linear interpolant onto

. Note that we can define only knowing on . Alternatively, if we have an analytic representation of we can define as a smooth extension of outside of .

## 4. Error Analysis

### 4.1. Stability Analysis

Unfortunately, the bilinear form is not positive definite. However, we will be able to prove stability of method. In order to do so, we need to decompose the space into its boundary contribution and interior contribution. More precisely, we can write

 Vkh=˚Wkh⊕Bkh,

where . We will define a norm on :

 ∥v∥2a:=ah(v,v)

and a semi-norm

 |v|2c:=∫Γv2|δ|ds.

Note that is in fact a norm on .

The following crucial lemma will allow us to prove stability.

###### Lemma 1.

There exists a constant such that

 (8) ∥v∥a≤c1√h|v|c for all v∈Bkh.
###### Proof.

Let be the collection of edges that are a subset of and let be triangles such that has an edge in . Then, if and using inverse estimates we have

 ∥v∥2a=∑T∈TΓh∥∇v∥2L2(T)≤∑T∈TΓhCh2T∥v∥2L2(T)≤∑e∈EΓhChe∥v∥2L2(e).

The result is complete after we use that for . ∎

We note that may not be well defined for all . Therefore, we need to make an assumption on such that this is not the case.

###### Assumption 2.

We assume that is such that

 (9) |ch(u,v)|<∞∀u,v∈Vkh.

For example, if has a lower bound as follows, then (9) will hold. Suppose that the end points of are and . Then we assume that there exists a constant and a such that

 |x−x0|p|x−x1|p≤c|δ(x)| for all x∈e,

where is independent of . Under these conditions, Assumption 2 holds.

We can now prove the stability result.

###### Theorem 1.

We assume that Assumption 1 and Assumption 2 hold. Suppose that is a bounded linear function on and suppose that solves

 bh(uh,v)=G(v), for all v∈Vkh.

Then, assuming we have

 ∥uh∥a≤2⎛⎜⎝supvh∈˚Wkh|G(vh)|∥vh∥a⎞⎟⎠+113c1√h⎛⎝supvh∈Bkh|G(vh)||vh|c⎞⎠.

and

 |uh|c≤32⎛⎜⎝supvh∈˚Wkh|G(vh)|∥vh∥a⎞⎟⎠+53⎛⎝supvh∈Bkh|G(vh)||vh|c⎞⎠.
###### Proof.

We know we can write where and . Define by

 ϕh=⎧⎨⎩shonΓ+−shonΓ−0onΓ0.

Note that . Now we can estimate .

 |sh|2c=ch(sh,ϕh)=bh(uh,ϕh)−ah(uh,ϕh)=G(ϕh)−ah(uh,ϕh).

Hence, we have

 |sh|2c≤ ⎛⎝supvh∈Bkh|G(vh)||vh|c⎞⎠|ϕh|c+∥uh∥a∥ϕh∥a ≤ ⎛⎝supvh∈Bkh|G(vh)||vh|c⎞⎠|sh|c+c1√h(∥wh∥a+c1√h|sh|c)|sh|c.

Here we used (8) twice. In particular, we used . Assuming we have

 34|sh|2c≤⎛⎝supvh∈Bkh|G(vh)||vh|c⎞⎠|sh|c+12∥wh∥a|sh|c.

Hence,

 (10) |sh|c≤43⎛⎝supvh∈Bkh|G(vh)||vh|c⎞⎠+23∥wh∥a

Next,

 ∥wh∥2a=ah(wh,wh)=ah(uh,wh)−ah(sh,wh)=bh(uh,wh)−ah(sh,wh)=G(wh)−ah(sh,wh).

We therefore have

 ∥wh∥2a≤⎛⎜⎝supvh∈˚Wkh|G(vh)||vh|a⎞⎟⎠∥wh∥a+∥sh∥a∥wh∥a.

Hence, we obtain using (10)

 ∥wh∥a≤ ⎛⎜⎝supvh∈˚Wkh|G(vh)|∥vh∥a⎞⎟⎠+∥sh∥a ≤⎛⎜⎝supvh∈˚Wkh|G(vh)|∥vh∥a⎞⎟⎠+c1√h∥sh∥c ≤⎛⎜⎝supvh∈˚Wkh|G(vh)|∥vh∥a⎞⎟⎠+43c1√h⎛⎝supvh∈Bkh|G(vh)||vh|c⎞⎠+13∥wh∥a

Thus we arrive at

 ∥wh∥a≤32⎛⎜⎝supvh∈˚Wkh|G(vh)|∥vh∥a⎞⎟⎠+2c1√h⎛⎝supvh∈Bkh|G(vh)||vh|c⎞⎠.

From this and (10) we get

 |uh|c=|sh|c≤32⎛⎜⎝supvh∈˚Wkh|G(vh)|∥vh∥a⎞⎟⎠+53⎛⎝supvh∈Bkh|G(vh)||vh|c⎞⎠.

Finally,

 ∥uh∥a≤ ∥wh∥a+∥sh∥a≤∥wh∥a+c1√h∥sh∥c ≤ 2⎛⎜⎝supvh∈˚Wkh|G(vh)|∥vh∥a⎞⎟⎠+113c1√h⎛⎝supvh∈Bkh|G(vh)||vh|c⎞⎠.

We can now prove error estimates after we make an assumption more stringent than Assumption 2.

###### Assumption 3.

Suppose that the end points of are and . Then we assume that there exists a constant such that

 |x−x0||x−x1|≤β|δ(x)| for all x∈e,

where is independent of .

Note that this assumption does not allow and to be tangent on the vertices of . Assumption 3 implies Assumption 2; in particular, the example after Assumption 2 holds with .

###### Theorem 2.

We assume Assumptions 1 and 3 hold. We assume that solves (1) and belongs to where . We assume that where is the Lagrange interpolant of . Let solve (7) and assume that solves (1) then we have

 ∥u−uh∥a≤ Chk∥u∥Hk+1(^Ω)+Chk+1/2∥u∥Wk+1,∞(Γ) +C(h4∥u∥W4,∞(^Ω)+h7/2∥u∥W2,∞(^Ω)).

and

 |u−uh|c≤ Chk∥u∥Hk+1(^Ω)+Chk∥u∥Wk+1,∞(Γ) +C(h4∥u∥W4,∞(^Ω)+h3∥u∥W2,∞(^Ω)).
###### Proof.

We let . Then we see that

 bh(eh,v)=G(v) for all v∈Vkh,

where , and .

Note that using integration by parts we have

 G1(v)= ∫ΩhFv+∫Γδ−1^gvds−∫Ωh(−Δu)vdx−∫Γ(∂u∂n+1δ(u−^g))v = ∫Ωh∖Ω(I1(−Δu)−(−Δu))vdx−∫Γ(∂u∂n+1δ(u−^g))v.

First consider then we have

 |G1(v)|≤h2∥u∥W4,∞(^Ω)∥v∥L1(Ωh∖Ω)

However, we have

 ∥v∥L1(Ωh∖Ω)≤ Ch2∥v∥L∞(Ωh∖Ω) ≤ Ch3∥∇v∥L∞(Ω) ≤ Ch2∥∇v∥L2(Ω)=Ch2∥v∥a.

Therefore, we get

 supv∈˚Wkh|G1(v)||v|a≤Ch4∥u∥W4,∞(^Ω).

Now consider .

 G1(v)= h4∥u∥W4,∞(^Ω)∥v∥L∞(Γ)+∥δ−1/2(δ∂u∂n+u−^g)∥L∞(Γ)∥v∥c ≤ h7/2∥u∥W4,∞(^Ω)∥v∥L2(Γ)+h3∥u∥W2,∞(^Ω)∥v∥c ≤ h9/2∥u∥W4,∞(^Ω)|v|c+h3∥u∥W2,∞(^Ω)∥v∥c.

Here we used (6).

Hence,

 √h(supv∈Bkh|G1(v)||v|c)≤C(h7/2∥u∥W2,∞(^Ω)+h5∥u∥W4,∞(^Ω)).

Now lets consider . If we let then

 G2(v)=ah(u−uI,v)≤∥u−uI∥a∥v∥a

Hence,

 supv∈˚Wkh|G2(v)|∥v∥a≤Chk∥u∥Hk(^Ω).

Now let we then have

 G2(v)=∥u−uI∥a∥v∥a+|u−uI|c|v|c≤c1√h∥u−uI∥a∥v∥c+|u−uI|c|v|c.

Let , with end points and . Then, we have . Hence, using Assumption 3 we get

 (u−uI)2(x)|δ(x)|≤Cβ∥∂t(u−uI)∥L∞(e).

Thus,

 ∫e(u−uI)2|δ|ds≤Cβ|e|∥∂t(u−uI)∥2L∞(e).

We then obtain the following estimate, after summing over all edges ,

 |u−uI|2c≤C∥∂t(u−uI)∥2L∞(Γ).

We get the following inequality after using approximation properties of the Lagrange interpolant:

 |u−uI|c≤Chk∥u∥Wk+1,∞(Γ).

Therefore, we have

 √hsupv∈Bkh|G2(v)||v|c≤Chk+1/2(∥u∥Wk+1,∞(Γ)+∥u∥Hk+1(^Ω)).

Combining the above results we get

 supv∈˚Wkh|G(vh)|∥vh∥a≤C(hk∥u∥Hk+1(^Ω)+h4∥u∥W4,∞(^Ω)).
 √hsupv∈Bkh|G(v)||v|c≤ C(h7/2∥u∥W2,∞(^Ω)+h5∥u∥W4,∞(^Ω)) +Chk+1/2(∥u∥Wk+1,∞(Γ)+∥u∥Hk+1(^Ω)).

The result now follows from Theorem 1.

## 5. Implementation

One feature of Nitsche’s method, that is preserved with BDT, is that one uses the full space of piecewise polynomials without restriction at the boundary. The modification of