1 Introduction
In this article we study numerical approximations of the Laplace equation,
(1.1)
(1.1a)
-
Δ
u
=
0
in
Ω
,
(1.1b)
u
=
g
on
Γ
,
based on linear finite elements in convex polygonal domains
Ω
⊂
ℝ
2
. Especially, we focus on the case that the regularity of the Dirichlet boundary datum g is restricted such that the solution fails to be in
H
2
(
Ω
)
but only belongs to
H
r
(
Ω
)
with some
r
∈
(
1
,
2
)
.
Due to the lack of regularity, the convergence rate
of the numerical approximations is not optimal regarding the approximation properties of the finite elements.
To be more specific, under the assumption
g
∈
H
3
/
2
(
Γ
)
(see Section 2 for the definition of the spaces
H
s
(
Γ
)
) one expects linear convergence in the
H
1
(
Ω
)
-
and quadratic convergence in the
L
2
(
Ω
)
-norm (at least if using an appropriate discretization of the Dirichlet boundary datum g such as its
L
2
(
Γ
)
-projection) for the finite
element approximations
u
h
of u, this is, the estimate
(1.2)
∥
u
-
u
h
∥
L
2
(
Ω
)
+
h
∥
u
-
u
h
∥
H
1
(
Ω
)
≤
c
h
2
∥
g
∥
H
3
/
2
(
Γ
)
is fulfilled on a sequence of quasi-uniform meshes with
h
>
0
denoting
the mesh parameter. The assumption
g
∈
H
3
/
2
(
Γ
)
is, however, in many applications very
restrictive. For less regular
data g the convergence rate is lower (
r
-
1
in
H
1
(
Ω
)
and r in
L
2
(
Ω
)
) if quasi-uniform meshes
are used for the computation of the approximations. Such situations may occur when the data are established from noisy measurements or when the data are highly oscillatory by nature. Nonsmooth boundary data also arise
in the computation of parametric minimal surfaces inside a boundary
skeleton with kinks and cusps, see also the example in
Section 5. In order to restore the
optimal convergence rate we propose the usage of boundary concentrated meshes, i.e., a local mesh refinement towards
the boundary of the computational domain, leading to a better
resolution of the irregular contributions of the solutions while
exploiting higher interior regularity. By the local refinement the computational complexity is not significantly increased if the grading is moderate enough (see Section 2 for more details).
The idea of using local mesh refinement to counteract
singularities in the solution is a well-known strategy used, e.g., to
resolve corner singularities in polygonal domains
[25, 9, 7],
edge singularities in polyhedral domains
[1, 4, 2] or to
resolve boundary layers in singularly perturbed problems [29, 28, 3]. A mesh refinement towards the whole boundary of the
computational domain, similar to the approach studied in this article,
has recently been applied in [26] for the
approximation of normal derivatives,
which is, e.g., of interest in the context of
Dirichlet boundary control problems.
Furthermore, we want to mention [12] where
adaptive finite elements for problems with inhomogeneous Dirichlet
boundary conditions are considered. There the authors only assume
g
∈
H
1
(
Γ
)
or even less and their adaptive method produces meshes refined to the
boundary, similar to the ones we obtain by a priori refinement. We also note that the idea of boundary concentrated meshes has also been applied in the context of the boundary concentrated finite element method [22, 23]. This method is based on a geometric refinement of the mesh towards the whole boundary such that the element size at the boundary is of order h while in addition the polynomial degree of the ansatz functions is increased towards the interior where the solution of the problem is more regular. Assuming a regularity of
H
1
+
β
(
Ω
)
for the solution, the boundary concentrated finite element method admits an error of order up to
h
β
in
H
1
(
Ω
)
. This is as for a standard discretization based on linear finite elements on quasi-uniform triangulations. However, the computational complexity is optimal (up to logarithmic factors). In contrast, the polynomial degree of the ansatz functions within our approach is fixed and the triangulations are gradually refined towards the whole boundary such that the element size at the boundary is of order
h
α
with some
α
≥
1
depending on the regularity of the Dirichlet boundary datum while in the interior the element size is of order h. Compared to a standard discretization with linear finite elements on quasi-uniform triangulations the computational complexity of our approach asymptotically stays the same if the grading is moderate enough (see again Section 2 for more details). However, optimal convergence rates in
L
2
(
Ω
)
and
H
1
(
Ω
)
can be achieved, which we elucidate next in detail.
As main result of this article we show that in convex polygonal domains
(1.2) remains true
for less regular boundary data
g
∈
H
t
(
Γ
)
with arbitrary
t
∈
(
1
2
,
3
2
]
if
the meshes are appropriately refined towards the boundary Γ.
The meshes we consider are constructed in such a way that the
mesh size is of order
h
1
μ
for elements at the boundary, and,
to guarantee that element diameters are not sharply varying, of order
h
ρ
T
1
-
μ
for elements away
from the boundary, where
ρ
T
denotes the distance of the element to the boundary. The parameter
μ
∈
(
0
,
1
]
defines the acuteness
of the refinement, where
μ
=
1
corresponds to a quasi-uniform family of meshes. We establish sharp bounds for μ, depending on the regularity index t of the boundary data
g
∈
H
t
(
Γ
)
, such that the convergence rates as in (1.2) are retained.
The bounds for μ will also depend on the discretization of the Dirichlet boundary data which is used in the discrete scheme. We distinguish between the
L
2
(
Γ
)
-projection and the nodal interpolant of the Dirichlet boundary data.
If the
L
2
(
Γ
)
-projection is used, we show optimal convergence of the
L
2
(
Ω
)
-error, provided that the mesh is refined according to
μ
≤
1
4
+
t
2
.
If the discrete boundary data are established by the nodal interpolant,
we show that the slightly stronger refinement
μ
≤
t
2
is
required to obtain optimal convergence in the
L
2
(
Ω
)
-norm (note that already for
t
=
3
2
mesh refinement is necessary, see also Remark 1). In either cases optimal convergence of the
H
1
(
Ω
)
-error is guaranteed under the assumption
μ
≤
t
-
1
2
.
The proofs of these results rely on new finite element error estimates
for the Poisson equation in weighted norms where the weights
are certain powers of the regularized distance to the boundary.
In order to focus on the core difficulty, a boundary
datum with reduced regularity, we consider in this article the Laplace
equation, but the results can directly be extended to problems with additional inhomogeneous right-hand side and with further but minor modifications to general second order elliptic
problems.
The article is structured as follows. In Section 2 we
describe the precise problem setting, introduce the required function spaces
as well as the finite element approach used to solve
(1.1). Section 3 is devoted to
regularity results for variational solutions of
(1.1). In particular, new a priori estimates in weighted
norms are proved. The main results of this article, error estimates
for the approximations of (1.1) on boundary-concentrated
meshes, are proved in Section 4. The theoretical
results are confirmed by numerical experiments presented in
Section 5.
2 Notation and Preliminaries
Throughout this article
Ω
⊂
ℝ
2
is a convex polygonal
domain with boundary Γ.
The corners are denoted by
𝒄
j
,
j
=
1
,
…
,
d
, numerated
counter-clockwise. The edge connecting
𝒄
j
and
𝒄
j
+
1
(
𝒄
d
+
1
=
𝒄
1
by convention) is denoted by
Γ
j
and has length
L
j
.
The classical Sobolev spaces are denoted as usual by
H
k
(
Ω
)
,
k
∈
ℕ
0
, which can be extended to a non-integer order
s
>
0
when
equipped with the Sobolev–Slobodetskij norm.
By
∥
⋅
∥
L
2
(
Ω
)
and
(
⋅
,
⋅
)
L
2
(
Ω
)
we
denote the norm and the inner product in the space
L
2
(
Ω
)
=
H
0
(
Ω
)
.
The Sobolev spaces on the boundary can be defined in the same way by local charts. However, as Γ is only of class
C
0
,
1
the standard definition is limited to the case that the regularity index s belongs to the interval
[
0
,
1
]
. Instead we use a piecewise (regarding the edges
Γ
j
) definition with additional compatibility conditions at the corners for the introduction of Sobolev spaces on the boundary. Thereby, we essentially rely on the ideas of [18, Section 1.5.2]. For each edge
Γ
j
,
j
=
1
,
…
,
d
, and for real-valued
s
∈
[
0
,
2
]
we use the classical definition for the Sobolev spaces
H
s
(
Γ
j
)
,
j
=
1
,
…
,
d
. The space
H
s
(
Γ
)
is then the subspace of
∏
j
=
1
d
H
s
(
Γ
j
)
defined by the following conditions: For each
j
=
1
,
…
,
d
and
f
=
(
f
1
,
…
,
f
d
)
∈
∏
j
=
1
d
H
s
(
Γ
j
)
there holds
no extra condition
if
s
∈
[
0
,
1
2
)
,
f
j
(
𝒄
j
+
1
)
=
f
j
+
1
(
𝒄
j
+
1
)
if
s
∈
(
1
2
,
2
]
,
I
j
(
f
)
:=
∫
0
δ
|
f
j
(
x
j
(
t
)
)
-
f
j
+
1
(
x
j
(
t
)
)
|
2
t
d
t
<
∞
if
s
=
1
2
.
In the latter case
δ
>
0
is arbitrarily small and
x
j
(
t
)
is a parametrization of the boundary near
𝒄
j
+
1
chosen in such a way
that for
t
>
0
,
x
j
(
t
)
is the point on
Γ
j
+
1
having distance
t to the corner
𝒄
j
+
1
and for
t
<
0
we arrive at the
point on
Γ
j
having distance
-
t
to
𝒄
j
+
1
.
One can show that
∥
u
∥
H
s
(
Γ
)
=
(
∑
j
=
1
d
∥
u
∥
H
s
(
Γ
j
)
2
)
1
2
if
s
∈
[
0
,
1
2
)
or
s
∈
(
1
2
,
2
]
,
∥
u
∥
H
s
(
Γ
)
=
∑
j
=
1
d
(
∥
u
∥
H
s
(
Γ
j
)
2
+
I
j
(
u
)
2
)
1
2
if
s
=
1
2
are norms in
H
s
(
Γ
)
. Moreover, one has that
H
s
(
Γ
)
corresponds to the natural trace space of
H
s
+
1
/
2
(
Ω
)
for
s
∈
(
0
,
2
)
. In particular, for
s
∈
(
0
,
2
)
the trace operator can be continuously extended to an operator from
H
s
+
1
/
2
(
Ω
)
to
H
s
(
Γ
)
, which is surjective, i.e., each function in
H
s
(
Γ
)
can be extended to a function from
H
s
+
1
/
2
(
Ω
)
, see [18, Theorem 1.5.2.8], [20, Theorem 4.2.7] and [8, Section 2].
The weak formulation of the Laplace problem (1.1) reads:
(2.1)
Find
u
∈
H
1
(
Ω
)
such that
u
≡
g
on
Γ
and
(
∇
u
,
∇
v
)
L
2
(
Ω
)
=
0
for all
v
∈
H
0
1
(
Ω
)
.
Under the assumption
g
∈
H
t
(
Γ
)
with
t
∈
[
1
2
,
3
2
]
the weak formulation is well-defined. Note that one might also consider the case
t
∈
[
0
,
1
2
)
,
but then, equation (1.1) has to be considered in the very weak sense which exceeds the scope of this article. For further reading in this context, we refer to [5, 6, 11, 14]. In the sequel, we will even restrict ourselves to the case
t
∈
(
1
2
,
3
2
]
. This ensures the well-posedness of the nodal interpolant on the boundary, which is introduced at the end of this section.
The numerical approximations to (2.1) are defined as
follows: We consider a family of admissible triangulations
𝒯
h
of
Ω. These triangulations are assumed to be shape-regular, i.e., a
minimal angle condition is satisfied. Moreover, in order to compensate the
singular behavior near the boundary Γ, coming from the irregular
boundary data g, we consider locally refined meshes. To this end, we define
for each
T
∈
𝒯
h
the element diameter
h
T
=
diam
(
T
)
and the distance to
the boundary
ρ
T
=
dist
(
T
,
Γ
)
. The mesh parameter is denoted by
h
=
max
T
∈
𝒯
h
h
T
.
For some refinement parameter
μ
∈
(
0
,
1
]
we
assume that the family of meshes fulfills
(2.2)
(2.2a)
c
¯
h
ρ
T
1
-
μ
≤
h
T
≤
c
¯
h
ρ
T
1
-
μ
if
ρ
T
>
0
,
(2.2b)
c
¯
h
1
μ
≤
h
T
≤
c
¯
h
1
μ
if
ρ
T
=
0
,
with some constants
c
¯
,
c
¯
>
0
that are independent of h.
The case
μ
=
1
corresponds to a quasi-uniform family of meshes. The smaller
μ is, the stronger the mesh is refined, yielding a boundary concentrated mesh.
The refinement has an effect on the number of elements
N
cells
, and thus, on the
computational complexity:
If
μ
>
1
2
, there holds
N
cells
∼
h
-
2
. The computational complexity asymptotically coincides
with the complexity if quasi-uniform meshes are used.
If
μ
=
1
2
, there holds
N
cells
∼
h
-
2
|
ln
h
|
, see [26, Remark 4].
We have a slight increase in the complexity compared to
quasi-uniform meshes.
If
μ
<
1
2
, there holds
N
cells
∼
h
-
1
μ
. The complexity blows up compared to quasi-uniform meshes.
The finite-dimensional function spaces used in our method are
V
h
:=
{
v
h
∈
C
(
Ω
¯
)
:
v
h
|
T
∈
𝒫
1
(
T
)
for all
T
∈
𝒯
h
}
,
V
0
h
=
V
h
∩
H
0
1
(
Ω
)
.
The traces of functions from
V
h
form the finite-dimensional function space
V
h
∂
=
tr
(
V
h
)
.
The finite-element approximations of (2.1),
u
h
∈
V
h
,
are defined by
(2.3)
u
h
≡
P
h
∂
g
on
Γ
,
(
∇
u
h
,
∇
v
h
)
L
2
(
Ω
)
=
0
for all
v
h
∈
V
0
h
.
Here,
P
h
∂
:
H
t
(
Γ
)
→
V
h
∂
is a projection operator onto
V
h
∂
. In the sequel, we will either use the
L
2
(
Γ
)
-projection
Q
h
∂
defined by
(
g
-
Q
h
∂
g
,
v
h
)
L
2
(
Γ
)
=
0
for all
v
h
∈
V
h
∂
,
or the nodal interpolant
I
h
∂
satisfying
[
I
h
∂
g
]
(
a
)
=
g
(
a
)
for all boundary nodes a of
𝒯
h
. Note that
I
h
∂
is only well
defined for
t
∈
(
1
2
,
3
2
]
due to
H
t
(
Γ
)
↪
C
(
Γ
)
.
3 Regularity Results
Before proving regularity results we introduce some technical tools
used in the rest of the article.
First, we define a dyadic decomposition of the computational domain
Ω, cf. [27]. For some
I
∈
ℕ
we define the subsets
Ω
J
,
J
=
0
,
…
,
I
, by
Ω
J
:=
{
x
∈
Ω
:
dist
(
x
,
Γ
)
∈
(
d
J
+
1
,
d
J
)
}
with
d
J
=
2
-
J
for
J
=
0
,
…
,
I
,
d
I
+
1
=
0
.
Without loss of generality we may assume that there is an interior domain
Ω
-
1
:=
Ω
∖
(
⋃
J
=
0
I
Ω
¯
J
)
which is not empty. The outer-most subset is assumed to have width
d
I
=
c
I
h
1
μ
with some sufficiently large but mesh-independent
number
c
I
>
0
specified later.
Note that this choice implies that the index
I is mesh-dependent, more precisely,
I
∼
|
ln
h
|
.
These sets form a dyadic decomposition of our computational domain
(3.1)
Ω
=
⋃
J
=
-
1
I
Ω
J
.
Furthermore, we define the patches with the neighboring sets
Ω
J
′
:=
Ω
max
{
-
1
,
J
-
1
}
∪
Ω
¯
J
∪
Ω
min
{
J
+
1
,
I
}
,
Ω
J
′′
:=
Ω
max
{
0
,
J
-
1
}
′
∪
Ω
¯
J
∪
Ω
min
{
J
+
1
,
I
-
1
}
′
.
In the sequel we will use regularity results in weighted norms involving as weight functions the distance function ρ to the boundary or its regularized version σ defined by
ρ
(
x
)
=
dist
(
x
,
Γ
)
and
σ
(
x
)
=
d
I
+
ρ
(
x
)
,
Note that in the regularized distance function σ we use the
mesh-dependent constant
d
I
from (3.1).
The weight functions are related to the dyadic decomposition in the
sense that
d
J
∼
σ
(
x
)
∼
ρ
(
x
)
for all
x
∈
Ω
J
,
J
=
-
1
,
…
,
I
-
1
,
d
I
∼
σ
(
x
)
.
We first investigate regularity results for the Poisson equation
(3.2)
(3.2a)
-
Δ
w
=
f
in
Ω
,
(3.2b)
w
=
0
on
Γ
,
which will serve as an auxiliary problem used for a duality argument
in the proof of Theorem 2.
Lemma 1.
For any
f
∈
L
2
(
Ω
)
the solution
w
∈
H
0
1
(
Ω
)
of
(3.3)
(
∇
w
,
∇
v
)
L
2
(
Ω
)
=
(
f
,
v
)
L
2
(
Ω
)
for all
v
∈
H
0
1
(
Ω
)
fulfills for
β
∈
[
1
,
3
2
)
the estimate
∥
σ
β
∇
2
w
∥
L
2
(
Ω
)
≤
c
∥
σ
β
f
∥
L
2
(
Ω
)
.
Proof.
In [26, proof of Lemma 2] one finds the local estimate
∥
∇
2
w
∥
L
2
(
Ω
J
)
≤
c
(
d
J
-
1
∥
∇
w
∥
L
2
(
Ω
J
′
)
+
∥
f
∥
L
2
(
Ω
J
′
)
)
.
Using
σ
∼
d
J
in
Ω
J
′
yields after summation
(3.4)
∥
σ
β
∇
2
w
∥
L
2
(
Ω
)
≤
c
(
∥
σ
β
-
1
∇
w
∥
L
2
(
Ω
)
+
∥
σ
β
f
∥
L
2
(
Ω
)
)
.
To bound the term
∥
σ
β
-
1
∇
w
∥
L
2
(
Ω
)
, we use the
product rule and obtain
∥
σ
β
-
1
∇
w
∥
L
2
(
Ω
)
2
=
(
σ
2
(
β
-
1
)
∇
w
,
∇
w
)
L
2
(
Ω
)
(3.5)
=
(
∇
(
σ
2
(
β
-
1
)
w
)
,
∇
w
)
L
2
(
Ω
)
-
2
(
β
-
1
)
(
w
σ
2
β
-
3
∇
σ
,
∇
w
)
L
2
(
Ω
)
.
For the first term on the right-hand side we apply (3.3), the Cauchy–Schwarz inequality and get
(3.6)
(
∇
(
σ
2
(
β
-
1
)
w
)
,
∇
w
)
L
2
(
Ω
)
=
(
σ
2
(
β
-
1
)
w
,
f
)
L
2
(
Ω
)
≤
∥
σ
β
-
2
w
∥
L
2
(
Ω
)
∥
σ
β
f
∥
L
2
(
Ω
)
.
To bound the second term in (3) we proceed as follows: First we observe that for
β
=
1
the term vanishes. For
β
∈
(
1
,
3
2
)
we introduce the sets
Ω
Γ
j
:=
{
𝐱
∈
Ω
:
ρ
(
𝐱
)
=
dist
(
𝐱
,
Γ
j
)
}
and local coordinate systems
(
x
j
,
y
j
)
obtained by affine transformations
F
j
(
𝐱
j
)
=
B
j
𝐱
j
+
𝒄
j
with
certain rotation matrices
B
j
chosen in such a way that
F
j
(
0
,
0
)
=
𝒄
j
and
F
j
(
L
j
,
0
)
=
𝒄
j
+
1
. This construction gives
σ
(
F
j
(
x
j
,
y
j
)
)
=
d
I
+
y
j
if
F
(
x
j
,
y
j
)
∈
Ω
Γ
j
. By
y
¯
j
(
𝐱
)
:
[
0
,
L
j
]
→
(
0
,
∞
)
we denote the function which allows to represent
Ω
J
by means of the local coordinates,
Ω
Γ
j
:=
{
F
j
(
x
j
,
y
j
)
∈
ℝ
2
:
x
j
∈
(
0
,
L
j
)
,
y
j
∈
(
0
,
y
¯
j
(
x
j
)
}
.
The function w after transformation to the new coordinates reads
w
~
j
(
x
j
,
y
j
)
=
w
(
F
j
(
x
j
,
y
j
)
)
.
With this transformation and the integration-by-parts formula we obtain for
β
∈
(
1
,
3
2
)
(
⇒
2
β
-
3
<
0
)
∫
Ω
w
σ
2
β
-
3
∇
σ
⋅
∇
w
d
x
=
∑
j
=
1
d
∫
0
L
j
∫
0
y
¯
j
(
x
)
(
d
I
+
y
j
)
2
β
-
3
w
~
j
(
x
j
,
y
j
)
∂
y
j
w
~
j
(
x
j
,
y
j
)
d
y
j
d
x
j
=
∑
j
=
1
d
∫
0
L
j
∫
0
y
¯
j
(
x
)
1
2
(
d
I
+
y
j
)
2
β
-
3
∂
y
j
(
w
~
j
(
x
j
,
y
j
)
2
)
d
y
j
d
x
j
=
∑
j
=
1
d
∫
0
L
j
(
-
∫
0
y
¯
j
(
x
j
)
2
β
-
3
2
(
d
I
+
y
j
)
2
(
β
-
2
)
w
~
j
(
x
j
,
y
j
)
2
d
y
j
(3.7)
+
1
2
(
d
I
+
y
¯
j
(
x
j
)
)
2
β
-
3
w
~
j
(
x
j
,
y
¯
j
(
x
j
)
)
2
)
d
x
j
≥
0
.
Insertion of (3) and
(3.6) into (3)
and using
β
∈
[
1
,
3
2
)
yields
(3.8)
∥
σ
β
-
1
∇
w
∥
L
2
(
Ω
)
2
≤
c
∥
σ
β
-
2
w
∥
L
2
(
Ω
)
∥
σ
β
f
∥
L
2
(
Ω
)
.
It remains to show an estimate for the term
∥
σ
β
-
2
w
∥
L
2
(
Ω
)
. The integration-by-parts formula
gives
∥
σ
β
-
2
w
∥
L
2
(
Ω
)
2
=
∑
j
=
1
d
∫
0
L
j
∫
0
y
¯
j
(
x
)
(
d
I
+
y
j
)
2
(
β
-
2
)
w
~
(
x
j
,
y
j
)
2
d
y
j
d
x
j
=
∑
j
=
1
d
∫
0
L
i
(
-
∫
0
y
¯
j
(
x
)
2
2
β
-
3
(
d
I
+
y
j
)
2
β
-
3
w
~
j
(
x
j
,
y
j
)
∂
y
j
w
~
j
(
x
j
,
y
j
)
d
y
j
+
1
2
β
-
3
(
d
I
+
y
¯
j
(
x
j
)
)
2
β
-
3
w
~
j
(
x
j
,
y
¯
j
(
x
j
)
)
2
)
d
x
j
≤
c
∥
σ
β
-
2
w
∥
L
2
(
Ω
)
∥
σ
β
-
1
∇
w
∥
L
2
(
Ω
)
.
Note that the boundary integral term is non-positive due to
2
β
-
3
<
0
and
y
¯
j
(
x
j
)
≥
0
for
x
j
∈
(
0
,
L
j
)
and can be bounded by zero. Division by
∥
σ
β
-
2
w
∥
L
2
(
Ω
)
yields the weighted Poincaré
inequality
(3.9)
∥
σ
β
-
2
w
∥
L
2
(
Ω
)
≤
c
∥
σ
β
-
1
∇
w
∥
L
2
(
Ω
)
.
Finally, we insert
(3.9) into (3.8)
and obtain
∥
σ
β
-
1
∇
w
∥
L
2
(
Ω
)
2
≤
c
∥
σ
β
f
∥
L
2
(
Ω
)
2
.
Together with (3.4) we infer the assertion.
∎
Lemma 2.
Let
g
∈
H
t
(
Γ
)
with some
t
∈
(
1
2
,
3
2
]
and
let
u
∈
H
1
(
Ω
)
be the solution
of
u
|
Γ
≡
g
,
(
∇
u
,
∇
v
)
L
2
(
Ω
)
=
0
for all
v
∈
H
0
1
(
Ω
)
.
Then there holds
u
∈
H
t
+
1
/
2
(
Ω
)
and
ρ
3
2
-
t
∇
2
u
∈
L
2
(
Ω
)
. Furthermore,
∥
u
∥
H
t
+
1
/
2
(
Ω
)
+
∥
ρ
3
2
-
t
∇
2
u
∥
L
2
(
Ω
)
≤
c
∥
g
∥
H
t
(
Γ
)
.
Proof.
The regularity
u
∈
H
t
+
1
/
2
(
Γ
)
follows from the fact
that g possesses an extension
G
∈
H
t
+
1
/
2
(
Ω
)
to the domain, see the explanations in Section 2,
and a standard shift theorem in convex domains, see, e.g., [13], applied to
-
Δ
u
0
=
Δ
G
∈
H
t
-
3
/
2
(
Ω
)
,
u
0
|
Γ
=
0
.
Combining these results yields
u
=
G
+
u
0
∈
H
t
+
1
/
2
(
Ω
)
.
To prove the weighted regularity we additionally take into account
∥
ρ
3
2
-
t
∇
2
u
∥
L
2
(
Ω
)
≤
c
∥
u
∥
H
t
+
1
/
2
(
Ω
)
,
proved in [19, Lemma 2.3].
∎
4 Error Estimates
The aim of this section is to prove error estimates for the
approximations (2.3). We consider sequences of meshes
refined according to (2.2). Here, we make use of
the dyadic decomposition (3.1) again. Note
that such a family of meshes is locally quasi-uniform within
the sets
Ω
J
, in the sense that for each
J
=
-
1
,
…
,
I
there holds
h
T
∼
h
d
J
1
-
μ
for all
T
∩
Ω
J
≠
∅
.
This fact directly allows to prove the following error estimates for the nodal
interpolant
I
h
:
C
(
Ω
¯
)
→
V
h
.
Lemma 3.
Under the assumption
u
∈
H
2
(
Ω
J
′
)
there holds
∥
u
-
I
h
u
∥
L
2
(
Ω
J
)
+
h
d
J
1
-
μ
∥
∇
(
u
-
I
h
u
)
∥
L
2
(
Ω
J
)
≤
c
h
2
d
J
2
(
1
-
μ
)
∥
∇
2
u
∥
H
2
(
Ω
J
′
)
.
The fact that meshes are locally quasi-uniform in each
Ω
J
as well as the relation
σ
(
x
)
∼
d
J
for
x
∈
Ω
J
allows to directly infer global estimates involving the weight
function σ:
Lemma 4.
For each
u
∈
H
2
(
Ω
)
the following
interpolation error estimate is fulfilled:
∥
σ
2
(
μ
-
1
)
(
u
-
I
h
u
)
∥
L
2
(
Ω
)
+
h
∥
σ
μ
-
1
∇
(
u
-
I
h
u
)
∥
L
2
(
Ω
)
≤
c
h
2
∥
∇
2
u
∥
L
2
(
Ω
)
.
As a first result, we prove a weighted
L
2
(
Ω
)
-norm error
estimate for the approximation of the Poisson equation
(3.2)
used later in a duality argument to prove the main result.
Theorem 1.
Let
β
∈
[
1
,
3
2
)
and
w
∈
H
0
1
(
Ω
)
∩
H
2
(
Ω
)
be
arbitrary.
Let
w
h
∈
V
0
h
be its Ritz projection defined by
(
∇
v
h
,
∇
(
w
-
w
h
)
)
L
2
(
Ω
)
=
0
for all
v
h
∈
V
0
h
.
For the refinement parameter
μ
=
2
-
β
∈
(
1
2
,
1
]
there holds the estimate
∥
σ
-
β
(
w
-
w
h
)
∥
L
2
(
Ω
)
≤
c
h
∥
∇
2
w
∥
L
2
(
Ω
)
.
Proof.
The proof relies on a duality argument. To this end, we introduce the
functions
φ
=
σ
-
β
(
w
-
w
h
)
and
z
∈
H
0
1
(
Ω
)
∩
H
2
(
Ω
)
solving
-
Δ
z
=
σ
-
β
φ
in
Ω
,
z
=
0
on
Γ
.
Together with the dyadic decomposition (3.1) we
arrive at
∥
σ
-
β
(
w
-
w
h
)
∥
L
2
(
Ω
)
2
=
(
σ
-
β
φ
,
w
-
w
h
)
L
2
(
Ω
)
=
(
∇
z
,
∇
(
w
-
w
h
)
)
L
2
(
Ω
)
=
(
∇
(
z
-
I
h
z
)
,
∇
(
w
-
w
h
)
)
L
2
(
Ω
)
≤
∑
J
=
-
1
I
∥
∇
(
w
-
w
h
)
∥
L
2
(
Ω
J
)
∥
∇
(
z
-
I
h
z
)
∥
L
2
(
Ω
J
)
(4.1)
≤
∑
J
=
-
1
I
(
∥
∇
(
w
-
I
h
w
)
∥
L
2
(
Ω
J
′
)
+
d
J
-
1
∥
w
-
w
h
∥
L
2
(
Ω
J
′
)
)
∥
∇
(
z
-
I
h
z
)
∥
L
2
(
Ω
J
)
.
The last step follows from the local error estimate [15, Theorem
3.4]. The terms involving an interpolation error in
(4) can be treated with
Lemma 3. Using also
h
=
d
I
μ
c
I
-
μ
yields
∥
σ
-
β
(
w
-
w
h
)
∥
L
2
(
Ω
)
2
≤
c
∑
J
=
-
1
I
(
h
∥
∇
2
w
∥
L
2
(
Ω
J
′′
)
+
c
I
-
μ
d
J
μ
-
2
∥
w
-
w
h
∥
L
2
(
Ω
J
′
)
)
d
J
2
-
μ
∥
∇
2
z
∥
L
2
(
Ω
J
′
)
.
We get together with
σ
(
x
)
∼
d
J
for
x
∈
Ω
J
and the discrete Cauchy–Schwarz inequality
(4.2)
∥
σ
-
β
(
w
-
w
h
)
∥
L
2
(
Ω
)
2
≤
c
(
h
∥
∇
2
w
∥
L
2
(
Ω
)
+
c
I
-
μ
∥
σ
-
β
(
w
-
w
h
)
∥
L
2
(
Ω
)
)
∥
σ
β
∇
2
z
∥
L
2
(
Ω
)
.
The a priori estimate from Lemma 1 gives
∥
σ
β
∇
2
z
∥
L
2
(
Ω
)
≤
c
∥
φ
∥
L
2
(
Ω
)
=
c
∥
σ
-
β
(
w
-
w
h
)
∥
L
2
(
Ω
)
. Dividing the inequality by
∥
σ
-
β
(
w
-
w
h
)
∥
L
2
(
Ω
)
and choosing
c
I
>
0
sufficiently large such that
c
c
I
-
μ
<
1
2
allows to apply a kick-back argument and (4.2) becomes
∥
σ
-
β
(
w
-
w
h
)
∥
L
2
(
Ω
)
≤
c
h
∥
∇
2
w
∥
L
2
(
Ω
)
.
∎
The previous result is also valid for
β
=
3
2
,
μ
=
1
2
, but with an
additional logarithmic factor
|
ln
h
|
in the
right-hand side, see
[26, Theorem 7].
Furthermore, we need the following interpolation error estimate in a weighted norm:
Lemma 5.
Let
u
∈
H
t
+
1
/
2
(
Ω
)
and
ρ
3
2
-
t
∇
2
u
∈
L
2
(
Ω
)
for any
t
∈
(
1
2
,
3
2
]
.
Then there holds the estimate
∥
σ
1
-
μ
∇
(
u
-
I
h
u
)
∥
L
2
(
Ω
)
≤
c
h
(
∥
u
∥
H
t
+
1
/
2
(
Ω
)
+
∥
ρ
3
2
-
t
∇
2
u
∥
L
2
(
Ω
)
)
,
provided that
μ
≤
1
/
4
+
t
/
2
.
Proof.
We define the set
S
h
:=
⋃
{
T
∈
𝒯
h
:
ρ
T
=
0
}
.
For elements
ρ
T
>
0
we can exploit the interior regularity
u
∈
H
2
(
T
)
, the relation
max
x
∈
T
σ
(
x
)
≤
c
ρ
T
and the mesh criterion
h
T
≤
c
h
ρ
T
1
-
μ
and deduce
∥
σ
1
-
μ
∇
(
u
-
I
h
u
)
∥
L
2
(
Ω
∖
S
h
)
≤
c
h
∥
ρ
3
2
-
t
∇
2
u
∥
L
2
(
Ω
)
.
For elements
T
∈
𝒯
h
with
ρ
T
=
0
we use a different
approach. The relation
max
x
∈
T
σ
(
x
)
≤
c
h
1
μ
and standard interpolation error estimates
exploiting regularity in fractional-order Sobolev spaces, see [16], give
∥
σ
1
-
μ
∇
(
u
-
I
h
u
)
∥
L
2
(
S
h
)
≤
c
h
1
-
μ
μ
h
t
-
1
/
2
μ
|
u
|
H
t
+
1
/
2
(
S
h
)
.
Insertion of
μ
≤
1
4
+
t
2
yields the desired estimate.
∎
Now we are in a position to prove the first main result of this article.
Theorem 2 (A Priori Estimate in
L
2
(
Ω
)
).
Let
g
∈
H
t
(
Γ
)
with
t
∈
(
1
2
,
2
]
. Denote by
u
∈
H
1
(
Ω
)
the solution of (2.1)
and by
u
h
∈
V
h
its finite element approximations solving (2.3). Assume that the refinement parameter fulfills
μ
≤
1
4
+
t
2
if
P
h
∂
=
Q
h
∂
and
μ
≤
t
2
if
P
h
∂
=
I
h
∂
.
Then there holds the error estimate
∥
u
-
u
h
∥
L
2
(
Ω
)
≤
c
h
2
∥
g
∥
H
t
(
Γ
)
.
Proof.
We use a duality argument and introduce the function
w
∈
H
0
1
(
Ω
)
solving
-
Δ
w
=
u
-
u
h
in
Ω
,
w
=
0
on
Γ
.
This implies
∥
u
-
u
h
∥
L
2
(
Ω
)
2
=
(
u
-
u
h
,
-
Δ
w
)
L
2
(
Ω
)
(4.3)
=
(
∇
(
u
-
u
h
)
,
∇
w
)
L
2
(
Ω
)
+
(
g
-
P
h
∂
g
,
∂
n
w
)
L
2
(
Γ
)
.
We consider the first term on the right-hand side of
(4). We denote by
w
h
∈
V
0
h
the
function satisfying
(
∇
v
h
,
∇
(
w
-
w
h
)
)
L
2
(
Ω
)
=
0
for all
v
h
∈
V
0
h
. Employing the Galerkin orthogonality and
(
E
h
P
h
∂
g
-
u
h
)
,
(
I
h
u
-
E
h
I
h
∂
g
)
∈
V
0
h
,
see also
[10, Theorem 6.1], with an
arbitrary discrete extension operator
E
h
:
V
h
∂
→
V
h
yields
(
∇
(
u
-
u
h
)
,
∇
w
)
L
2
(
Ω
)
=
(
(
∇
(
u
-
I
h
u
+
E
h
(
I
h
∂
g
-
P
h
∂
g
)
)
,
∇
(
w
-
w
h
)
)
L
2
(
Ω
)
(4.4)
≤
(
∥
σ
1
-
μ
∇
(
u
-
I
h
u
)
∥
L
2
(
Ω
)
+
∥
σ
1
-
μ
∇
E
h
(
I
h
∂
g
-
P
h
∂
g
)
∥
L
2
(
Ω
)
)
∥
σ
μ
-
1
∇
(
w
-
w
h
)
∥
L
2
(
Ω
)
.
For the term depending on the interpolation error we apply
Lemma 5 and the a priori
estimate from Lemma 2 to conclude
(4.5)
∥
σ
1
-
μ
∇
(
u
-
I
h
u
)
∥
L
2
(
Ω
)
≤
c
h
∥
g
∥
H
t
(
Γ
)
.
We study the term
∥
σ
1
-
μ
E
h
(
I
h
∂
g
-
P
h
∂
g
)
∥
L
2
(
Ω
)
arising in (4).
The term vanishes in case of
P
h
∂
=
I
h
∂
.
Otherwise, if
P
h
∂
=
Q
h
∂
, we choose
E
h
to be the zero-extension satisfying
∥
∇
E
h
v
h
∥
L
2
(
Ω
)
≤
c
h
-
1
2
μ
∥
v
h
∥
L
2
(
Γ
)
,
see the last equation in [24, proof of Lemma 3.3].
Taking into account
σ
(
x
)
≤
c
h
1
μ
for
x
∈
T
with
ρ
T
=
0
we deduce the estimate
∥
σ
1
-
μ
∇
E
h
(
I
h
∂
g
-
P
h
∂
g
)
∥
L
2
(
Ω
)
≤
c
h
1
-
μ
μ
h
-
1
2
μ
∥
I
h
∂
g
-
Q
h
∂
g
∥
L
2
(
Γ
)
=
c
h
1
-
2
μ
2
μ
∥
Q
h
∂
(
I
h
∂
g
-
g
)
∥
L
2
(
Γ
)
(4.6)
≤
c
h
2
t
+
1
-
2
μ
2
μ
∥
g
∥
H
t
(
Γ
)
≤
c
h
∥
g
∥
H
t
(
Γ
)
.
The last two steps follow from the stability of the
L
2
(
Γ
)
-projection
Q
h
and standard estimates for the
nodal interpolant
I
h
(note that
t
>
1
2
), as well as
the refinement criterion
μ
≤
1
4
+
t
2
.
Estimates (4),
(4.5) and (4) give
(4.7)
(
∇
(
u
-
u
h
)
,
∇
w
)
L
2
(
Ω
)
≤
c
h
∥
g
∥
H
t
(
Γ
)
∥
σ
μ
-
1
∇
(
w
-
w
h
)
∥
L
2
(
Ω
)
.
To bound the latter term on the right-hand side
we apply the local
H
1
(
Ω
)
-norm error
estimate from [15, Theorem 3.4] on the
dyadic decomposition from (3.1) to obtain
d
J
μ
-
1
∥
∇
(
w
-
w
h
)
∥
L
2
(
Ω
J
)
≤
c
(
d
J
μ
-
1
∥
∇
(
w
-
I
h
w
)
∥
L
2
(
Ω
J
′
)
+
d
J
μ
-
2
∥
w
-
w
h
∥
L
2
(
Ω
J
′
)
)
and with
σ
(
x
)
∼
d
J
for
x
∈
Ω
J
we get after summation over
J
=
-
1
,
…
,
I
∥
σ
μ
-
1
∇
(
w
-
w
h
)
∥
L
2
(
Ω
)
≤
c
(
∥
σ
μ
-
1
∇
(
w
-
I
h
w
)
∥
L
2
(
Ω
)
+
∥
σ
μ
-
2
(
w
-
w
h
)
∥
L
2
(
Ω
)
)
.
We insert the interpolation error estimate from
Lemma 4 and the estimate for
σ
μ
-
2
(
w
-
w
h
)
from Theorem 1 to
obtain
(4.8)
∥
σ
μ
-
1
∇
(
w
-
w
h
)
∥
L
2
(
Ω
)
≤
c
h
∥
∇
2
w
∥
L
2
(
Ω
)
≤
c
h
∥
u
-
u
h
∥
L
2
(
Ω
)
.
Insertion of (4.8) into
(4.7) yields
(4.9)
(
∇
(
u
-
u
h
)
,
∇
w
)
L
2
(
Ω
)
≤
c
h
2
∥
g
∥
H
t
(
Γ
)
∥
u
-
u
h
∥
L
2
(
Ω
)
We consider the second term on the right-hand
side of (4). In case of
P
h
∂
=
Q
h
∂
we
use the orthogonality of
Q
h
∂
and standard estimates to
get
(
g
-
Q
h
∂
g
,
∂
n
w
)
L
2
(
Γ
)
=
(
g
-
Q
h
∂
g
,
∂
n
w
-
Q
h
∂
(
∂
n
w
)
)
L
2
(
Γ
)
≤
c
h
t
μ
∥
g
∥
H
t
(
Γ
)
h
1
2
μ
∥
∂
n
w
∥
H
1
/
2
(
Γ
)
(4.10)
≤
c
h
2
∥
g
∥
H
t
(
Ω
)
∥
u
-
u
h
∥
L
2
(
Ω
)
.
In the last step we used again
μ
≤
1
4
+
t
2
and the standard
trace and a priori estimate
∥
∂
n
w
∥
H
1
/
2
(
Γ
)
≤
c
∥
w
∥
H
2
(
Ω
)
≤
c
∥
u
-
u
h
∥
L
2
(
Γ
)
.
In case of
P
h
∂
=
I
h
∂
we instead estimate as follows:
(
g
-
I
h
∂
g
,
∂
n
w
)
L
2
(
Γ
)
≤
c
h
t
μ
∥
g
∥
H
t
(
Γ
)
∥
∂
n
w
∥
H
1
/
2
(
Γ
)
(4.11)
≤
c
h
2
∥
g
∥
H
t
(
Γ
)
∥
u
-
u
h
∥
L
2
(
Ω
)
,
where
μ
≤
t
2
was used in the latter step.
The previous estimates (4.10) and
(4), together with
(4.9) are inserted into
(4) and the result is proved.
∎
Remark 1.
When the boundary data are established with the nodal interpolant
P
h
∂
=
I
h
∂
,
a stronger refinement seems to be required. This is due to estimate
(4) which is not sharp. In numerical experiments
one observes that also a weaker refinement is sufficient. If there additionally holds
g
∈
W
s
,
1
(
Γ
)
with some
s
∈
(
1
2
,
2
]
, then, instead of
(4), we estimate as follows. We use the
embedding
H
1
/
2
(
Γ
)
↪
L
q
(
Γ
)
which is valid for
arbitrary
q
<
∞
and
W
s
,
1
(
Γ
)
↪
W
s
-
η
,
1
+
ε
(
Γ
)
for
ε
>
0
and
η
=
ε
1
+
ε
.
Then we obtain
(
g
-
I
h
∂
g
,
∂
n
w
)
L
2
(
Γ
)
≤
∥
g
-
I
h
∂
g
∥
L
1
+
ε
(
Γ
)
∥
∂
n
w
∥
L
(
1
+
ε
)
/
ε
(
Γ
)
≤
c
h
s
-
η
μ
∥
g
∥
W
s
-
η
,
1
+
ε
(
Γ
)
∥
∂
n
w
∥
H
1
/
2
(
Ω
)
≤
c
h
2
∥
g
∥
W
s
,
1
(
Γ
)
∥
u
-
u
h
∥
L
2
(
Ω
)
.
The last step holds if
μ
<
s
2
for sufficiently small ε.
In addition, the assumption
g
∈
W
s
,
1
(
Γ
)
also implies
g
∈
H
t
(
Γ
)
with
s
=
t
+
1
2
due to an embedding. Thus, the grading condition
μ
<
s
2
=
1
4
+
t
2
is sufficient, which is weaker compared to
μ
<
t
2
.
Theorem 3 (A Priori Estimate in
H
1
(
Ω
)
).
Let
g
∈
H
t
(
Γ
)
with
t
∈
(
1
2
,
3
2
]
. Denote by
u
∈
H
1
(
Ω
)
the solution of (2.1)
and by
u
h
∈
V
h
its finite element approximations solving (2.3). Assume that the
refinement parameter fulfills
μ
≤
t
-
1
2
. Then there holds the
error estimate
∥
∇
(
u
-
u
h
)
∥
L
2
(
Ω
)
≤
c
h
∥
g
∥
H
t
(
Γ
)
.
Proof.
Let
E
h
:
V
h
∂
→
V
h
be the zero-extension. The
Galerkin orthogonality and the property
I
h
u
-
u
h
-
E
h
(
P
h
∂
g
-
I
h
∂
g
)
∈
V
0
h
then give
∥
∇
(
u
-
u
h
)
∥
L
2
(
Ω
)
2
=
(
∇
(
u
-
u
h
)
,
∇
(
u
-
I
h
u
+
E
h
(
P
h
∂
g
-
I
h
∂
g
)
)
)
L
2
(
Ω
)
(4.12)
≤
∥
∇
(
u
-
u
h
)
∥
L
2
(
Ω
)
(
∥
∇
(
u
-
I
h
u
)
∥
L
2
(
Ω
)
+
h
-
1
2
μ
∥
P
h
∂
g
-
I
h
∂
g
∥
L
2
(
Γ
)
)
,
where the latter term was obtained by an application of the estimate
∥
∇
E
h
v
h
∥
L
2
(
Ω
)
≤
c
h
-
1
2
μ
∥
v
h
∥
L
2
(
Γ
)
, used already in
the proof of Theorem 2.
To bound the interpolation error, we proceed as in Lemma 5, i.e., we distinguish between elements touching the boundary belonging to
S
h
=
⋃
{
T
∈
𝒯
h
:
ρ
T
=
0
}
and elements in the interior
Ω
∖
S
h
. For elements in the interior we exploit the higher interior regularity
u
∈
H
2
(
T
)
and the mesh criterion
h
T
≤
c
h
ρ
T
1
-
μ
to deduce
∥
∇
(
u
-
I
h
u
)
∥
L
2
(
Ω
∖
S
h
)
≤
c
h
∥
ρ
1
-
μ
∇
2
u
∥
L
2
(
Ω
)
.
For elements
T
⊂
S
h
we use standard interpolation error estimates exploiting regularity in fractional-order Sobolev spaces, see [16]. This gives
∥
∇
(
u
-
I
h
u
)
∥
L
2
(
S
h
)
≤
c
h
t
-
1
/
2
μ
|
u
|
H
t
+
1
/
2
(
S
h
)
.
Combining the previous results yields together with
Lemma 2 and
μ
≤
t
-
1
2
,
∥
∇
(
u
-
I
h
u
)
∥
L
2
(
Ω
)
≤
c
h
(
∥
ρ
3
2
-
t
∇
2
u
∥
L
2
(
Ω
)
+
∥
u
∥
H
t
+
1
/
2
(
Ω
)
)
≤
c
h
∥
g
∥
H
t
(
Γ
)
.
For the term involving
P
h
∂
g
-
I
h
∂
g
we proceed analogous to
(4). For
P
h
∂
=
I
h
∂
the
term vanishes and in case of
P
h
∂
=
Q
h
∂
we use
h
-
1
2
μ
∥
P
h
∂
g
-
I
h
∂
g
∥
L
2
(
Γ
)
≤
c
h
2
t
-
1
2
μ
∥
g
∥
H
t
(
Γ
)
≤
c
h
∥
g
∥
H
t
(
Γ
)
.
We insert the previous two estimates into
(4), divide by
∥
∇
(
u
-
u
h
)
∥
L
2
(
Ω
)
and arrive at the assertion.
∎
Remark 2.
The grading conditions to guarantee first order convergence in
H
1
(
Ω
)
and second order convergence in
L
2
(
Ω
)
, respectively, are illustrated depending on the regularity index t in Figure 1.
5 Numerical Experiments
5.1 A problem with singularities in the boundary data
In a first experiment we study problem (1.1) on the
unit square domain
Ω
:=
(
0
,
1
)
2
and boundary data
(5.1)
g
(
x
,
y
)
=
|
x
-
4
7
|
1
2
+
|
y
-
2
3
|
1
2
which has singularities in the points
(
4
7
,
0
)
,
(
4
7
,
1
)
,
(
0
,
2
3
)
,
(
1
,
2
3
)
and there holds
g
∈
H
t
(
Γ
)
for
t
<
1
only.
In our C++ implementation based on the finite element library MooNMD
[21] we compute the approximation error in the
L
2
(
Ω
)
- and
H
1
(
Ω
)
-norms on a sequence of
hierarchical meshes obtained by a newest-vertex-bisection
algorithm, see Figure 2 (a). The meshes satisfy the refinement criterion
(2.2) and we choose different refinement
parameters to check whether the rates are optimal and the refinement
conditions are sharp. To compute the error norms, we prolongate the
numerical solution
u
h
to the finest mesh
𝒯
h
min
and,
instead of the true errors
∥
u
-
u
h
∥
, we compute
∥
u
h
min
-
u
h
∥
, which are
sufficiently good approximations if
h
≫
h
min
.
The approximate solution is illustrated in Figure 2 (b).
The computational results are summarized in Figure 3.
The results confirm our predictions in the following sense.
To obtain optimal convergence for the given boundary data g
projected with the
L
2
(
Γ
)
-projection
P
h
∂
=
Q
h
∂
,
Theorem 2 requires the refinement parameter μ to be
chosen according to
μ
≤
1
4
+
t
2
. This is
indeed confirmed by the error plot in
Figure 3 (a). Note that the error curves
corresponding to
μ
=
0.75
are parallel
to the hypotenuse of the slope triangle representing a convergence
rate of two. If the nodal interpolant
P
h
∂
=
I
h
∂
is used to establish
g
h
, the refinement criterion
μ
≤
t
2
from
Theorem 2 seems to be too pessimistic. However,
the function g from (5.1) belongs even to
W
s
,
1
(
Γ
)
for any
s
<
3
2
. As explained in
Remark 1 this regularity and the refinement
criterion
μ
<
s
2
, so in our case
μ
<
3
4
, guarantees optimal
convergence as well. This is also confirmed in the numerical
experiment.
The convergence behavior is obviously different in the
H
1
(
Ω
)
-norm.
According to Theorem 3 a stronger mesh refinement is
necessary to obtain the best-possible convergence order of one. This can be
seen in Figure 3 (b). The convergence rates of
∥
u
-
u
h
∥
H
1
(
Ω
)
for the choices
μ
=
0.75
and
μ
=
0.6
are obviously between 0.5 and the optimal rate 1. Only for
μ
=
0.5
the error curve is roughly parallel to the hypotenuse of
the slope triangle representing the order one.
Figure 2
Mesh and solution for the numerical experiment from Section 5.1.
Figure 3
Computational results for the benchmark problem from
Section 5.1.
Plot shows error propagation
u
-
u
h
in the
L
2
(
Ω
)
- and
H
1
(
Ω
)
-norm for
different refinement parameters μ. The slope triangles indicate the
convergence rates. Solid lines represent the computational error for the case
that
P
h
∂
=
Q
h
∂
(
L
2
(
Γ
)
-projection) is used to project the boundary data
and the dashed lines for
P
h
∂
=
I
h
∂
(nodal interpolant).
5.2 Noisy Boundary Data
In a second experiment we consider the Laplace equation (1.1)
on
Ω
=
(
0
,
1
)
2
with boundary data
g
(
x
,
y
)
=
sin
(
2
π
x
)
+
sin
(
2
π
y
)
+
ε
(
x
,
y
)
.
Here,
ε
∈
V
h
noise
∂
,
h
noise
=
2
⋅
10
-
3
, is a discrete Gaussian noise satisfying
ε
(
𝒂
)
∼
𝒩
(
0
,
0.05
)
for each node
𝒂
of an equidistant boundary mesh
𝒯
h
noise
∂
. The numerical solution of
this problem is depicted in Figure 4.
As ε is piecewise linear and globally continuous, there holds
g
∈
H
t
(
Γ
)
for arbitrary
t
<
3
2
. We expect that the finite element
approximations
u
h
of the exact solution u converge almost linearly on
quasi-uniform meshes. However,
this is an asymptotic result and the optimal convergence is only achieved for sufficiently fine meshes
at the boundary so that the noise can be resolved accurately enough.
In this example the refinement to the boundary has the positive
effect that the region where the approximations converge, so even
for larger mesh parameters h, is reached earlier. This is exactly reflected in
Figure 5 (b) where the error
u
-
u
h
in the
H
1
(
Ω
)
-norm is plotted. For
μ
=
1
,
0.75
,
0.6
and 0.5 the
approximations achieve linear convergence for
h
<
4
⋅
10
-
3
,
1
⋅
10
-
2
,
2
⋅
10
-
2
and
3
⋅
10
-
2
, respectively.
The choice of the projection operator
P
h
∂
makes a difference for the
L
2
(
Ω
)
-norm error on coarse meshes. In Figure 5 (a)
we observe that for meshes, which do not yet resolve the noise accurately
enough, the results seem to be better when the
L
2
(
Γ
)
-projection is
used to establish the discrete boundary data.
Figure 5
Computational results for the benchmark problem from
Section 5.2.
Plot shows error propagation
u
-
u
h
in the
L
2
(
Ω
)
- and
H
1
(
Ω
)
-norm for
different refinement parameters μ. The slope triangles indicate the
convergence rates. Solid lines represent the computational error for the case
that
P
h
∂
=
Q
h
∂
(
L
2
(
Γ
)
-projection) is used to project the boundary data
and the dashed lines for
P
h
∂
=
I
h
∂
(nodal interpolant).
6 Conclusion and Outlook
The main results of this article are the error estimates for the linear
finite element approximation of the Laplace equation, see
Theorems 2 and 3. In
particular, we proved a relation between the
regularity index t of the boundary datum
g
∈
H
t
(
Γ
)
and the
required strength of a boundary refinement to guarantee optimal
convergence rates, i.e., rate one in the
H
1
(
Ω
)
- and rate two in the
L
2
(
Ω
)
-norm. Worth to mention is, that
different refinement strengths are required, see
Figure 1, which is also confirmed in the numerical
experiments.
It turns out that using the
L
2
(
Γ
)
-projection to project the boundary data into the finite
element space yields slightly better results than the nodal
interpolant from the theoretical of view. To prove a similar result for the
nodal interpolant we have to assume additional regularity for the boundary
datum, more precisely,
g
∈
W
s
,
1
(
Γ
)
with
s
∈
(
1
,
2
]
, see
Remark 1.
In a further experiment we have seen that boundary concentrated meshes
also yield better computational results for noisy and oscillatory data.
Our theory does not cover irregular boundary data
g
∈
H
t
(
Γ
)
with
t
∈
[
0
,
1
2
)
. In this situation the weak
formulation (2.1) is not even well defined and
equation (1.1) has to be understood in a very weak
sense. Error estimates for a different finite element discretizations
have been derived in [5, 6], see also
[17, 11, 14] for
related results. Boundary concentrated
meshes may also be applied in this situation to get the full order of convergence.
Also the extension to non-convex domains is an interesting issue. In
this case, rigorous modifications are required. First, the
H
2
(
Ω
)
-regularity for the dual solution used in
Theorem 1 is not given any more. Instead
regularity in weighted Sobolev spaces, with weights related to the
singular corners, have to be employed. Optimal convergence can in this
case only be guaranteed, when an additional local refinement towards
the corners is used in the case that the corner singularities entering
the solution are less regular than the boundary datum g.
However, the proofs become much more technical in this case.
The extension to three-dimensional domains is also possible. However, the proof of the main result requires rigorous modifications. In addition, although optimal convergence with respect to the mesh size can be achieved for arbitrary
t
∈
(
1
2
,
3
2
)
, the main problem is the cardinality of the resulting meshes and hence the computational complexity. To be more precise, the relation
N
cells
∼
h
-
3
, which holds for quasi-uniform meshes in 3D, is violated already for refinement parameters
μ
≤
2
3
, see [7, Remark 3.1]. However, according to Theorem 2 for instance, such a refinement is already necessary for
t
≤
5
6
as in this case
μ
≤
1
4
+
t
2
implies
μ
≤
2
3
. This disadvantage may be overcome if local refinement towards critical parts of the boundary is used instead of a refinement towards the whole boundary.
Moreover, it could be interesting to apply the ideas of this article in the context of the boundary concentrated finite element method. For geometric meshes with element size
h
1
μ
at the whole boundary, where μ satisfies
μ
≥
1
2
and the corresponding conditions from Theorem 2 or Theorem 3, respectively, this might result in a finite element method which not only allows for optimal finite element error estimates with respect to h but also has optimal computational complexity (up to logarithmic factors). However, the proof of sharp error estimates, especially in the
L
2
(
Ω
)
-norm, is an open problem.
Acknowledgements
We would like to take the opportunity to thank our PhD supervisor
Professor Thomas Apel for his assistance during our time as PhD students
at the Universität der Bundeswehr München.
References
[1]
T. Apel and B. Heinrich,
Mesh refinement and windowing near edges for some elliptic problem,
SIAM J. Numer. Anal. 31 (1994), no. 3, 695–708.
10.1137/0731037Search in Google Scholar
[2]
T. Apel, A. L. Lombardi and M. Winkler,
Anisotropic mesh refinement in polyhedral domains: error estimates with data in
L
2
(
Ω
)
,
ESAIM Math. Model. Numer. Anal. 48 (2014), no. 4, 1117–1145.
10.1051/m2an/2013134Search in Google Scholar
[3]
T. Apel and G. Lube,
Anisotropic mesh refinement for a singularly perturbed reaction diffusion model problem,
Appl. Numer. Math. 26 (1998), no. 4, 415–433.
10.1016/S0168-9274(97)00106-2Search in Google Scholar
[4]
T. Apel and S. Nicaise,
Elliptic problems in domains with edges: Anisotropic regularity and anisotropic finite element meshes,
Partial Differential Equations and Functional Analysis,
Progr. Nonlinear Differential Equations Appl. 22,
Birkhäuser, Boston (1996), 18–34.
10.1007/978-1-4612-2436-5_2Search in Google Scholar
[5]
T. Apel, S. Nicaise and J. Pfefferer,
Discretization of the Poisson equation with non-smooth data and emphasis on non-convex domains,
Numer. Methods Partial Differential Equations 32 (2016), no. 5, 1433–1454.
10.1002/num.22057Search in Google Scholar
[6]
T. Apel, S. Nicaise and J. Pfefferer,
Adapted numerical methods for the Poisson equation with
L
2
boundary data in nonconvex domains,
SIAM J. Numer. Anal. 55 (2017), no. 4, 1937–1957.
10.1137/16M1062077Search in Google Scholar
[7]
T. Apel, A.-M. Sändig and J. R. Whiteman,
Graded mesh refinement and error estimates for finite element solutions of elliptic boundary value problems in non-smooth domains,
Math. Methods Appl. Sci. 19 (1996), no. 1, 63–85.
10.1002/(SICI)1099-1476(19960110)19:1<63::AID-MMA764>3.0.CO;2-SSearch in Google Scholar
[8]
D. N. Arnold, L. R. Scott and M. Vogelius,
Regular inversion of the divergence operator with Dirichlet boundary conditions on a polygon,
Ann. Sc. Norm. Super. Pisa Cl. Sci. (4) 15 (1988), no. 2, 169–192.
Search in Google Scholar
[9]
I. Babuška, R. B. Kellogg and J. Pitkäranta,
Direct and inverse error estimates for finite elements with mesh refinements,
Numer. Math. 33 (1979), no. 4, 447–471.
10.1007/BF01399326Search in Google Scholar
[10]
S. Bartels, C. Carstensen and G. Dolzmann,
Inhomogeneous Dirichlet conditions in a priori and a posteriori finite element error analysis,
Numer. Math. 99 (2004), no. 1, 1–24.
10.1007/s00211-004-0548-3Search in Google Scholar
[11]
M. Berggren,
Approximations of very weak solutions to boundary-value problems,
SIAM J. Numer. Anal. 42 (2004), no. 2, 860–877.
10.1137/S0036142903382048Search in Google Scholar
[12]
S. Bertoluzza, E. Burman and C. He,
An a posteriori error estimate of the outer normal derivative using dual weights,
SIAM J. Numer. Anal. 60 (2022), no. 1, 475–501.
10.1137/20M1358219Search in Google Scholar
[13]
M. Dauge,
Elliptic Boundary Value Problems on Corner Domains,
Lecture Notes in Math. 1341,
Springer, Berlin, 1988.
10.1007/BFb0086682Search in Google Scholar
[14]
K. Deckelnick, A. Günther and M. Hinze,
Finite element approximation of Dirichlet boundary control for elliptic PDEs on two- and three-dimensional curved domains,
SIAM J. Control Optim. 48 (2009), no. 4, 2798–2819.
10.1137/080735369Search in Google Scholar
[15]
A. Demlow, J. Guzmán and A. H. Schatz,
Local energy estimates for the finite element method on sharply varying grids,
Math. Comp. 80 (2011), no. 273, 1–9.
10.1090/S0025-5718-2010-02353-1Search in Google Scholar
[16]
T. Dupont and R. Scott,
Polynomial approximation of functions in Sobolev spaces,
Math. Comp. 34 (1980), no. 150, 441–463.
10.1090/S0025-5718-1980-0559195-7Search in Google Scholar
[17]
D. A. French and J. T. King,
Approximation of an elliptic control problem by the finite element method,
Numer. Funct. Anal. Optim. 12 (1991), no. 3–4, 299–314.
10.1080/01630569108816430Search in Google Scholar
[18]
P. Grisvard,
Elliptic Problems in Nonsmooth Domains,
Monogr. Stud. Math. 24,
Pitman, Boston, 1985.
Search in Google Scholar
[19]
T. Horger, J. M. Melenk and B. Wohlmuth,
On optimal
L
2
- and surface flux convergence in FEM,
Comput. Vis. Sci. 16 (2013), no. 5, 231–246.
10.1007/s00791-015-0237-zSearch in Google Scholar
[20]
G. C. Hsiao and W. L. Wendland,
Boundary Integral Equations,
Appl. Math. Sci. 164,
Springer, Berlin, 2008.
10.1007/978-3-540-68545-6Search in Google Scholar
[21]
V. John and G. Matthies,
MooNMD—a program package based on mapped finite element methods,
Comput. Vis. Sci. 6 (2004), no. 2–3, 163–169.
10.1007/s00791-003-0120-1Search in Google Scholar
[22]
B. N. Khoromskij and J. M. Melenk,
An efficient direct solver for the boundary concentrated FEM in 2D,
Computing 69 (2002), no. 2, 91–117.
10.1007/s00607-002-1452-2Search in Google Scholar
[23]
B. N. Khoromskij and J. M. Melenk,
Boundary concentrated finite element methods,
SIAM J. Numer. Anal. 41 (2003), no. 1, 1–36.
10.1137/S0036142901391852Search in Google Scholar
[24]
S. May, R. Rannacher and B. Vexler,
Error analysis for a finite element approximation of elliptic Dirichlet boundary control problems,
SIAM J. Control Optim. 51 (2013), no. 3, 2585–2611.
10.1137/080735734Search in Google Scholar
[25]
L. A. Oganesjan and L. A. Ruhovec,
Variational-difference schemes for second order linear elliptic equations in a two-dimensional region with a piecewise-smooth boundary,
Zh. Vychisl. Mat. Mat. Fiz. 8 (1968), 97–114.
10.1016/0041-5553(68)90008-6Search in Google Scholar
[26]
J. Pfefferer and M. Winkler,
Finite element error estimates for normal derivatives on boundary concentrated meshes,
SIAM J. Numer. Anal. 57 (2019), no. 5, 2043–2073.
10.1137/18M1181341Search in Google Scholar
[27]
A. H. Schatz and L. B. Wahlbin,
Interior maximum norm estimates for finite element methods,
Math. Comp. 31 (1977), no. 138, 414–442.
10.1090/S0025-5718-1977-0431753-XSearch in Google Scholar
[28]
A. H. Schatz and L. B. Wahlbin,
On the finite element method for singularly perturbed reaction-diffusion problems in two and one dimensions,
Math. Comp. 40 (1983), no. 161, 47–89.
10.1090/S0025-5718-1983-0679434-4Search in Google Scholar
[29]
G. I. Shishkin,
Approximation of solutions of singularly perturbed boundary value problems with a corner boundary layer,
Zh. Vychisl. Mat. Mat. Fiz. 27 (1987), no. 9, 1360–1374, 1438.
10.1016/0041-5553(87)90044-9Search in Google Scholar
