Arab. J. Math.

(2022) 11:407–426
https://doi.org/10.1007/s40065-022-00386-w Arabian Journal of Mathematics

Dina N. Alsheikh · Theodoros Katsaounis

A posteriori error estimators for discontinuous Galerkin

method for diffusion problems, based on the hypercircle
method

Received: 4 April 2022 / Accepted: 28 July 2022 / Published online: 5 September 2022
© The Author(s) 2022

Abstract In this work, we apply the hypercircle method to Discontinuous Galerkin (DG) approximations of
second order diffusion problems featuring inhomogeneous Dirichlet and Neumann boundary conditions. We
focus on the interior penalty discontinuous Galerkin (IPDG) approximations of diffusion problems in primal
variational formulation and produce a Prager–Synge theorem for such DG methods. Using the hypercircle
method, we derive an a posteriori error estimator in terms of an equilibrated flux. The estimator is proven to
be reliable and efficient. Numerical results are presented which illustrate the estimator’s performance.

Mathematics Subject Classification (2000) 35J25 · 65N15 · 65N22 · 65N30 · 65N66 · 65N50

1 Introduction

The hypercircle method, also known as the two-energy principle, was introduced in 1947 by Prager and Synge
[25,26] as a method to estimate the error arising from the approximation of boundary value problems for
diffusive partial differential equations. The method is based on splitting the original problem into two easy-
to-solve problems. These problems are solved separately, a common solution is found which represents the
solution of the original problem in a suitably-defined function space. In this case, the solution to the original
problem is located on or inside a hypercircle in that function space. If the bilinear form of the variational
formulation of the original problem is positive definite, then the hypercircle is bounded [25] and the location
of the solution on a hypercircle provides an upper and lower bound for the solution on any neighbourhood
of any prescribed point in the domain. If we choose the center of the hypercircle as an approximation of the
original solution, then its error in a mean-square sense, is easily known.
However, the hypercircle method was forgotten for many years. It resurfaced again when there was a need for
efficient and reliable a posteriori error estimators for the finite element approximation of diffusion-boundary
value problems (cf. e.g. [2,6,8,27]). The equilibrated a posteriori error estimator based on the hypercircle
method has a big advantage over the standard residual-type a posteriori error estimator that does not contain
any generic constants except for possible data oscillations, which leads to better effectivity indices.
The hypercircle method was designed for conforming finite elements approximations, but it can be used as
well as for nonconforming approaches such as discontinuous Galerkin (DG) methods. Indeed, in [7], a Prager–
Synge theorem has been formulated and used to derive an equilibrated a posteriori error estimator for DG

D. N. Alsheikh (B) ·
Jubail Industrial College, Jubail Industrial City, Saudi Arabia
E-mail: sheikhd@ucj.edu.sa

T. Katsaounis
Univ. of Crete, IACM-FORTH, Heraklion, Greece
E-mail: thodoros.katsaounis@uoc.gr

123
408 Arab. J. Math. (2022) 11:407–426

approximations of the Poisson equation. This is done by constructing equilibrated flux and post-processing the
DG approximation via a related mixed formulation from [4], involving suitably defined numerical fluxes across
the interfaces of the finite element meshes. Related approaches can be found in [1,12,16], which however are
not directly based on a Prager–Synge-like theorem.
DG methods have been used since the 1970s and have been so popular since then, [13,14,22] and the
references therein, for their flexibility in local approximation of the solution of partial differential equations
(PDE’s), flexibility in adaptive local mesh refinement, and for their good stability properties. When it comes to
adaptive methods based on equilibrated a posteriori error estimators, DG methods have a big advantage over
C 0 - conforming finite elements. In C 0 -conforming finite elements, the equilibration is done on patches of
elements (union of elements sharing a common nodal point), as in [8], whereas for interior penalty DG (IPDG)
methods equilibration is performed element-wise which is less expensive to compute [7,12].
In this work, we extend the work done in [7] to cover DG approximations of diffusion problems with
non-zero boundary conditions and high contrast coeffiecients. Further, [7] does not address data oscillations
due to approximations, this is another point that we consider here. For such problems, we apply the hypercircle
method to IPDG methods of any polynomial degree, we formulate a Prager–Synge theorem and derive an a
posteriori error estimator in terms of an equilibrated flux vector. The reliability of the error estimator follows
directly from the Prager–Synge theorem, whereas its efficiency can be established by showing that the estimator
is bounded from above by a residual-type a posteriori error estimator for IPDG methods whose efficiency has
been shown in the literature, (cf., e.g., [5,17,19]).
This paper is organised as follows, a brief introduction to the DG approximation is given followed by
focus on the IPDG methods in Sect. 2. We introduce a posteriori error estimate based on the Prager and Synge
theorem in Sect. 3.1 and we show the post processing of the DG methods to conferring FE we used in Sect.
3.2. We introduce a discussion of possible Data oscillations. In Sect. 3.3 and in Sect. 3.4 we construct an
equilibrated flux. In Sect. 3.5 we show the efficiency of our estimator. Sect. 4.1 consists of an introduction to
the adaptive method we are using and some numerical results are given in Sect. 4.2 to show the performance
of the estimator. our conclusions were discussed in section 5.

2 DG methods for diffusion problems

We will use standard notations from Lebesgue and Sobolev space theory [6,9]. In particular, for a bounded
domain Ω ⊂ Rd for d ≥ 2 and D ⊂ Ω we denote by (·, ·)0,D and  · 0,D the L 2 -inner product and
the associated L 2 -norm. We further refer to H 1 (Ω) as the Sobolev space with norm  · 1,Ω and seminorm
1
| · |1,Ω and to H 2 (Γ ), Γ ⊆ ∂Ω, as the associated trace space, whereas H0,Γ 1 (Ω) stands for the subspace

of H -functions with vanishing trace on Γ . Moreover H (div, Ω) denotes the Hilbert space of vector fields
1
q ∈ L 2 (Ω)2 such that ∇ · q ∈ L 2 (Ω) equipped with the graph norm.
We start by assuming Ω ⊂ R2 to be a bounded polygonal domain consisting of pairwise disjoint polygonal
subdomains Ωk , 1  k  m. We refer to Γ = Γ D ∪ Γ N , Γ D ∩ Γ N = ∅, as the boundary of Ω, consisting of
a Dirichlet part Γ D and a Neumann part Γ N . We consider the diffusion problem
−∇ · (a∇u) = f in Ω, (1a)
u = u D on Γ D , (1b)
nΓ N · a∇u = u N on Γ N . (1c)

Here, we suppose a to be a piecewise constant function with a|Ωk > 0, 1  k  m, and f ∈ L 2 (Ω),
1
u D ∈ H 2 (Γ D ) and u N ∈ L 2 (Γ N ). Further, nΓ N denotes the exterior unit normal on Γ N and we set V =
{u ∈ H 1 (Ω) : u|Γ D = u D } and H D1 (Ω) = {v ∈ H 1 (Ω) : v|Γ D = 0}.
We define then, the weak formulation of (1a)-(1c) as : we seek u ∈ V such that for all v ∈ H D1 (Ω) it holds
  
a∇u · ∇v dx = f v dx + u N v ds. (2)
Ω Ω ΓN

For the IPDG approximation of the boundary value problem (1a)–(1c) we consider a locally quasi-uniform
simplicial triangulation Th of a computational domain which aligns with the subdivision of Ω. We denote by
Eh the set of edges of Th in Ω̄, where Eh is divided into EhI , EhD := Eh ∩ Γ D , EhN := Eh ∩ Γ N , the set of interior

123
Arab. J. Math. (2022) 11:407–426 409

edges in Th , Dirichlet boundary edges, and Neumann edges respectively. We further denote by h K , K ∈ Th
the diameter of K , and by h E , E ∈ Eh the length of edge E. Further, for l ∈ N, we refer to Pl (K ) as the set
of polynomials of degree at most l on K . Due to the local quasi-uniformity of the triangulation, there exist
constants 0 < c < C such that
c hE ≤ hK ≤ C hE, E ∈ ∂K. (3)

For a function vh ∈ L 2 (Ω) with vh | K ∈ Pl (K ), K ∈ Th , and an interior edge E = K + ∩ K − , K ± ∈ Th , we

set vh± := vh | E∩K ± , and define the average and jump across E as usual, according to
1 + 
v + vh− , E ∈ Eh ,
{vh } E := 2 h (4a)
vh | E , E ∈ EhD ∪ EhN ,
 +
vh − vh− , E ∈ Eh ,
[vh ] E := (4b)
vh | E , E ∈ EhD ∪ EhN .
Moreover, we refer to n E , E ∈ Eh , E = K + ∩ K − , as the unit normal vector pointing from K + to K − and to
n E , E ∈ EhD ∪ EhN , as the exterior unit normal vector. We introduce now the DG spaces

Vh = {vh ∈ L 2 (Ω)| vh | K ∈ Pk (K ), K ∈ Th }, (5)

Vh = {qh ∈ L 2 (Ω)2 | qh | K ∈ Pk (K )2 , K ∈ Th }. (6)
DG methods can be derived from a mixed formulation of (1a)-(1c), [4] and amount to the computation of
(u h , ph ) ∈ Vh × Vh such that for all (vh , qh ) ∈ Vh × Vh it holds
⎛ ⎞
   
⎝ a −1 ph · qh dx + u h ∇ · qh dx − ũ∂ K · qh ds ⎠ = 0, (7a)
K ∈T h K K ∂K
⎛ ⎞
  
⎝ (ph · ∇vh − f vh ) dx − p̃∂ K · n∂ K vh ds ⎠ = 0, (7b)
K ∈T h K ∂K

where ũ∂ K and p̃∂ K are suitably defined numerical fluxes.

2.1 The IPDG method

Setting ũ E := ũ∂ K | E , p̃ E := p̃∂ K | E , E ∈ Eh , the IPDG method is obtained by the following specification of
the numerical fluxes
⎧
⎨ {u h } E n E , E ∈ EhI ,
ũ E := (u h − u D )n E , E ∈ EhD , (8a)
⎩
uh nE , E ∈ EhN ,
⎧
⎨ {a∇u h } E − αh −1E [u h ] E n E , E ∈ EhI ,
−1
p̃ E := a∇u h − αh E (u h − u D )n E , E ∈ EhD , (8b)
⎩
u N nE , E ∈ EhN ,
where α > 0 is a penalty parameter. This choice of the numerical fluxes allows us to eliminate the dual variable
ph from (7a),(7b) by choosing ph = a∇u h and qh = a∇vh in (7a). We are thus led to the primal variational
formulation of the IPDG method: we seek u h ∈ Vh such that,
ahI P (u h , vh ) = l(vh ), ∀vh ∈ Vh , (9)

where the bilinear form ahI P (·, ·) : Vh × Vh → R is defined by

 
ah (u h , vh ) :=
IP
a∇u h · ∇vh dx
K ∈T h K

123
410 Arab. J. Math. (2022) 11:407–426

 
− ({n E · a∇u h } E [vh ] E + [u h ] E {n E · a∇vh } E ) ds
E∈EhI ∪EhD E
 
α
+ [u h ] E [vh ] E ds, (10)
hE
E∈EhI ∪EhD E

and the functional l(·) : Vh → R reads as follows:

   α 
l(vh ) := f vh dx + u D vh ds
hE
K ∈T h K E∈EhD E
   
− u D n E · a∇vh ds + u N vh ds. (11)
E∈EhD E E∈EhN E

Remark 2.1 Another choice of the numerical fluxes ũ E and p̃ E , E ∈ Eh , gives rise to a different DG approxima-
tion. We refer to [4] for an overview. In the remaining of this paper, we will focus on the IPDG approximation,
but note that both the derivation of the equilibrated a posteriori error estimator and its analysis in terms of
reliability and efficiency can be done the same way.
The IPDG bilinear form (10) and the functional in (11) are not well-defined for u, v ∈ V , since (n E ·∇u)| E , E ∈
Eh , does not belong to L 2 (E). This can be cured by using a lifting operator L : L 2 (EhI ∪ EhD ) → Vh as in [18],
defined according to
   
L(v) · aqh dx = [v] E n E · {aqh } E ds. (12)
K ∈T h K E∈EhI ∪EhD E

As has been shown e.g. in [18,23], the lifting operator is continuous, so there exists a constant C L ≥ 0
depending only on the shape regularity of the triangulation such that

L(v)20,Ω ≤ C L h −1
E [v] E 0,E , v ∈ L (Eh ∪ Eh ).
2 2 I D
(13)
E∈EhI ∪EhD

Using the lifting operator, we define the bilinear form ãhI P (·, ·) : (V + Vh ) × (V + Vh ) → R and the functional
l˜h (·) by means of
   
ãh (u, v) :=
IP
a∇u · ∇v dx − (L(u) · a∇v
K ∈T h K K ∈T h K
 α 
+ a∇u · L(v)) dx + [u] E [v] E ds, (14a)
hE
E∈Eh E
   α 
l˜h (v) := f v dx + u D v ds
hE
K ∈T h K E∈Eh
D
E
   
− L(u D ) · a∇v dx + u N v ds. (14b)
K ∈T h K E∈EhN E

It is easy to show that the following holds

ãhI P (u h , vh ) = ahI P (u h , vh ), l˜h (vh ) = lh (vh ), ∀u h , vh ∈ Vh .

123
Arab. J. Math. (2022) 11:407–426 411

2.2 Coercivity and continuity of the IPDG bilinear form

Introducing the mesh-dependent IPDG norm on V + Vh ,

 1  α
v21,h,Ω := a 2 ∇v20,K + [v] E 20,E , v ∈ V + Vh .
hE (15)
K ∈T h E∈EhI ∪EhD

It is not difficult to show, cf., [19,20], that for sufficiently large α > 0 the IPDG norm and the mesh-dependent
energy norm are equivalent, i.e., there exists a positive constant γ such that

ãhI P (v, v) ≥ γ v21,h,Ω , v ∈ V + Vh . (16a)

On the other hand, there exists a constant Γ > γ such that for any α > 0

ãhI P (v, w) ≤ Γ v1,h,Ω w1,h,Ω , v, w ∈ V + Vh . (16b)

In particular, it follows from (16a) and (16b) that for a sufficiently large α > 0, the bilinear form ahI P (·, ·)
is continuous and coercive, thus the IPDG problem (9) admits a unique solution u h ∈ Vh . It has been shown
before [24] that to ensure the coercivity of the bilinear form and the continuity of the IPDG approximation,
the penalty parameter has to be chosen according to α = amax O((k + 1)2 ).

3 A posteriori error estimate based on the Prager–Synge theorem

Theorem 3.1 Prager–Synge theorem for diffusion Problems, two-energy principle Let p ∈ H (div, Ω) and
v ∈ V . Further, let u ∈ V be the solution of (2) and suppose that p satisfies the equilibrium condition

−∇ · p = f in L 2 (Ω), (17)

and the boundary condition

nΓ N · p = u N in L 2 (Γ N ). (18)

Then it holds
1   1 1  
a 2 ∇v − a −1 p 20,Ω = a 2 ∇ (v − u) 20,Ω + a 2 ∇u − a −1 p 20,Ω . (19)

Proof Observing (2), (17),(18), and the fact that v − u ∈ H0,Γ

1
D
(Ω), partial integration yields

 
(a∇u − p, ∇(v − u))0,Ω = (a∇u, ∇(v − u))0,Ω + (∇ · p, v − u)0,Ω − nΓ N · p, v − u 0,Γ
  N

= u N − nΓ N · p, v − u 0,Γ = 0.
N

This orthogonality relation implies

1 1 1 1
a 2 ∇(v − u) + a 2 (∇u − a −1 p)20,Ω = a 2 ∇(v − u)20,Ω + a 2 (∇u − a −1 p)20,Ω .

which proves (19). 



123
412 Arab. J. Math. (2022) 11:407–426

3.1 A posteriori error estimator

For K ∈ Th let f K be the L 2 -projection of f onto Pk−1 (K ), and for E ∈ EhD resp., E ∈ EhN let u D | E resp. u N | E
be the L 2 -projection of u D resp. u N onto Pk (E). We define f h ∈ L 2 (Ω) and u D,h ∈ L 2 (Γ D ), u N ,h ∈ L 2 (Γ N )
by f h | K = f K , K ∈ Th , and u D,h | E = u D | E , E ∈ EhD , u N ,h | E = u N | E , E ∈ EhN . The Prager–Synge theorem
eq
(3.1) can be applied to the IPDG approximation involving an equilibrated flux vector ph , if we replace f, u D ,
and u N by f h , u D,h , and u N ,h , respectively. In particular, we denote by û h ∈ Vh the IPDG approximation of
the weak solution û ∈ H 1 (Ω) of the diffusion problem

−∇ · a∇ û = f h , in Ω,
û = u D,h , on Γ D ,
nΓ N · a∇ û = u N ,h , on Γ N .

Note that in this case, the associated mixed formulation reads

⎛ ⎞
   
⎝ a −1 ph · qh dx + u h ∇ · qh dx − û∂ K · qh ds ⎠ = 0, (20a)
K ∈T h K K ∂K
⎛ ⎞
  
⎝ (ph · ∇vh − f vh )dx − p̃∂ K · n∂ K vh ds ⎠ = 0. (20b)
K ∈T h K ∂K

where the numerical fluxes û∂ K and p̂∂ K are given by

⎧
⎨ {ûh } E n E , E ∈ EhI ,
û E := (û h − u D,h )n E , E ∈ EhD , (21a)
⎩
û h n E , E ∈ EhN ,
⎧
⎨ {a∇ û h } E − αh −1
E [û h ] E n E , E ∈ EhI ,
p̂ E := a∇ û h − αh −1 (û h − u D,h )n E , E ∈ EhD , (21b)
⎩ E
u N ,h n E , E ∈ EhN .

3.2 Postprocessing to conforming finite element spaces

Denoting by Vhc := Vh ∩ C(Ω) ⊂ V the C 0 -conforming Lagrange finite element space associated with Vh ,
the function v in (19) can be chosen as v = û conf
h ∈ Vhc , where û conf
h is obtained by postprocessing the IPDG
approximation û h ∈ Vh . Let Nh be a set of Lagrange nodal points for the elements in Vh and let κi be the
L

number of elements that share the nodal point xi ∈ NhL . We note that κi = 1, if xi belongs to the interior of
an element K ∈ Th or if xi is situated in the interior of an edge E ∈ EhD ∪ EhN , while κi > 1 if xi ∈ NhL ∩ Eh .
The function û conf
h is defined by its nodal values

1 
h (x i ) :=
û conf û h | K (xi ). (22)
κi
K ∈Th ,xi ∈K

The following estimate has been provided by Theorem (2.2) in [19]:

 
∇(û h − û conf
h )0,K 
2
h −1
E [û h ]0,E .
2
(23)
K ∈T h E∈Eh

eq eq
A flux vector ph ∈ Vh is called an equilibrated flux vector, if ph ∈ H (div, Ω) and if it satisfies the equilibrium
equation
eq
−∇ · ph = f h in each K ∈ Th , (24a)

123
Arab. J. Math. (2022) 11:407–426 413

and the boundary condition

eq
nΓ N · ph = u N ,h ineachE ∈ EhN . (24b)
As we shall see below, the Prager–Synge theorem (3.1) gives rise to an equilibrated a posteriori error estimator
⎛ ⎞1 ⎛ ⎞1 ⎛ ⎞1 ⎛ ⎞1
2 2 2 2
  ⎜  eq 2 ⎟ ⎜  eq 2 ⎟
:= ⎝ (η K ,1 )2 ⎠ +2⎝ (η K ,2 )2 ⎠
eq eq eq
ηh +⎝ (η E,1 ) ⎠ + ⎝ (η E,2 ) ⎠ (25)
K ∈T h K ∈T h E∈EhI E∈EhD

eq eq
where the element-related terms η K ,i , 1  i  2, and edge-related terms η E,i , 1  i  2, are given as follows:
eq eq
η K ,1 := ph − a∇ û h 0,K , K ∈ Th , (26a)
eq 1
η K ,2 := a ∇(û h − û conf
2
h )0,K , K ∈ Th , (26b)
eq − 21
η E,1 := αh E [û h ] E 0,E , E ∈ EhI , (26c)
eq − 21
η E,2 := αh E u D,h − û h 0,E , E ∈ EhD . (26d)

3.3 Data oscillations

We provide the following auxiliary result concerning the data oscillation due to the approximation of f, u D ,
and u N by f h , u D,h , and u N ,h respectively.
Lemma 3.2 Let z h ∈ Vh be the IPDG approximation of the weak solution z ∈ H 1 (Ω) of the diffusion problem
−∇ · a∇z = f − f h , in Ω, (27a)
z = u D − u D,h , on Γ D , (27b)
nΓ N · a∇z = u N − u N ,h , on Γ N . (27c)
Then it holds
 
z h 21,h,Ω  h 2K  f − f h 20,K + h −1
E u D − u D,h 0,E
2

K ∈T h E∈EhD

+ h E u N − u N ,h 20,E . (28)
E∈EhN

Proof For vh ∈ Vh and p0 ∈ L 2 (Ω) with p0 | K ∈ P0 (K ), K ∈ Th , we have

 
ahI P (z h , vh ) = ( f − f h , vh − p0 )0,K − (u D − u D,h , n E · a∇vh )0,E
K ∈T h E∈EhD
 
+ αh −1
E (u D − u D,h , vh − p0 )0,E + (u N − u N ,h , vh − p0 )0,E ,
E∈EhD E∈EhN

where we have used that ( f − f h , p0 )0,K = 0, K ∈ Th , (u D − u D,h , p0 )0,E = 0, E ∈ EhD , and (u N −

u N ,h , p0 )0,E = 0, E ∈ EhN . For vh = z h , it follows from (16a) that
 
z h 21,h,Ω   f − f h 0,K z h − p0 0,K + u D − u D,h 0,E n E · a∇z h 0,E
K ∈T h E∈EhD
 
+ αh −1
E u D − u D,h 0,E z h − p0 0,E + u N − u N ,h 0,E z h − p0 0,E . (29)
E∈EhD E∈EhN

123
414 Arab. J. Math. (2022) 11:407–426

The Poincaré-Friedrichs inequality asserts

⎛ ⎞

z h − p0 0,K  ⎝∇z h 0,K + | (z h − p0 )dx|⎠ .
K

−1

Hence, if we choose p0 = |K | z h dx, we obtain
K

z h − p0 0,K  h K ∇z h 0,K . (30)

Moreover, the trace inequalities for functions v K ∈ Pk (K ) (cf.,e.g., [28])

−1 −1
v K 0,∂ K  h K 2 v K 0,K , n E · a∇v K 0,∂ K  h K 2 a∇v K 0,K

imply that for E ∈ EhD , resp.E ∈ EhN and E ⊂ ∂ K for some K ∈ Th it holds

−1
z h − p0 0,E ≤ z h − p0 0,∂ K  h K 2 z h − p0 0,K , (31a)
−1
n E · a∇z h 0,E ≤ n∂ K · a∇z h 0,∂ K  h K 2 a∇z h 0,K . (31b)

In view of (30), (31a), (31b), and using (3) we deduce from (29)
  −1
z h 21,h,Ω ≤ h K  f − f h 0,K ∇z h 0,K + h E 2 u D − u D,h 0,E ∇z h 0,K
K ∈T h E∈EhD
 −1  1
+ αh E 2 u D − u D,h 0,E ∇z h 0,K + h E2 u N − u N ,h 0,E ∇z h 0,K . (32)
E∈EhD E∈EhN

The assertion (28) then follows from (32) by a suitable application of Young’s inequality. 

In view of the previous result, the data oscillations will be denoted by

⎛ ⎞1 ⎛ ⎞1 ⎛ ⎞1
2 2
2
 ⎜ ⎟ ⎜ ⎟
osch := ⎝ osc2K ( f )⎠ +⎝ osc2E (u D )⎠ + ⎝ osc2E (u N )⎠ , (33)
K ∈T h E∈EhD E∈EhN

where osc K ( f ), osc E (u D ), and osc E (u N ) are given by

osc K ( f ) := h K  f − f h 0,K , K ∈ Th ,
− 21
osc E (u D ) := h E u D − u D,h | E 0,E , E ∈ EhD ,
1
osc E (u N ) := h E u N − u N ,h | E 0,E ,
2
E ∈ EhN .

Theorem 3.3 Let u ∈ V be the weak solution of (1a)-(1c), let û h ∈ Vh be the solution of the IPDG approxi-
eq
mation (9) with f, u D , and u N replaced by f h , u D,h , u N ,h . Moreover, let ph ∈ H (div, Ω) be an equilibrated
eq
flux vector and û conf
h ∈ Vhc be the extension of û h to Vhc as defined by (22). Finally, let ηh be the equilibrated a
posteriori error estimator given by (25) and let osch be the data oscillations (33). Then there exists a constant
C > 0 independent of h such that
eq
u − û h 1,h,Ω ≤ ηh + C osch . (34)

123
Arab. J. Math. (2022) 11:407–426 415

Proof Let u h ∈ Vh be the solution of (9). Then we have

⎛ ⎞1 ⎛ ⎞1
2 2
 
≤⎝ − u h )20,K ⎠ +⎝ ∇(u h − û h )20,K ⎠
1 1
u − û h 1,h,Ω a ∇(u2 a 2

K ∈T h K ∈T h
⎛ ⎞1 ⎛ ⎞1 ⎛ ⎞1
2 2 2

⎜  eq 2 ⎟ ⎜  eq 2 ⎟ ⎜ ⎟
+⎝ (η E,1 ) ⎠ + ⎝ (η E,2 ) ⎠ + ⎝ (osc E (u D ))2 ⎠ . (35)
E∈EhI E∈EhD E∈EhD

Since z h := u h − û h is the IPDG approximation of (27a)–(27c), the second term on the right-hand side in (35)
can be estimated from above by (28) and thus gives rise to the data oscillation term in (34). For the first term
on the right-hand side in (35) we obtain
⎛ ⎞1 ⎛ ⎞1
2 2
 
⎝ a 2 ∇(u − u h )20,K ⎠ ≤ ⎝ 2 ⎠
1 1
a 2 ∇(u − û conf
h )0,K
K ∈T h K ∈T h
⎛ ⎞1 ⎛ ⎞1
2 2
 
+⎝ 2 ⎠
≤⎝ 2 ⎠
1 1
a 2 ∇(u h − û conf
h )0,K a ∇(u
2 − û conf
h )0,K
K ∈T h K ∈T h
⎛ ⎞1 ⎛ ⎞1
2 2
 
+⎝ 2 ⎠
+⎝ ∇(u h − û h )20,K ⎠
1 1
a 2 ∇(û h − û conf
h )0,K a 2 .
K ∈T h K ∈T h
(36)
Again, the last term on the right-hand side can be estimated from above by (28). The Prager–Synge theorem
eq
3.1 with v = û conf
h and p = ph gives
⎛ ⎞1
2

≤⎝ − a −1 ph )20,K ⎠
1 1 eq
a ∇(u
2 − û conf
h )0,Ω a 2 (∇ û conf
h
K ∈T h
⎛ ⎞1 ⎛ ⎞1
2 2
 
≤⎝ 2 ⎠
+⎝ (∇ û h − a −1 ph )20,K ⎠
1 1 eq
a 2 ∇(û h − û conf
h )0,K a 2 (37)
K ∈T h K ∈T h

Using (36) and (37) in (35) allows to conclude (34). 


Remark 3.4 We emphasize that the constant C in front the data oscillations is the only generic constant in the
reliability estimate (34).

3.4 Construction of an equilibrated flux

eq
We consider now the construction of an equilibrated flux ph ∈ BDMk (Ω, Th ) where and BDMk (K ) stands
for the classical Brezzi-Douglas-Marini element on K , [10]. Denoting by λiK , 1 ≤ i ≤ 3, the barycentric
coordinates of K ∈ Th , we refer to b K := λ1K λ2K λ3K as the element bubble function which vanishes on ∂ K .
For qh | K ∈ BDMk (K ), we define its degrees of freedom(DOF) according to

n E · qh pk ds, pk ∈ Pk (E), E ∈ ∂ K , (38a)
E

qh · ∇ pk−1 dx, pk−1 ∈ Pk−1 (K ), (38b)
K

123
416 Arab. J. Math. (2022) 11:407–426


qh · curl(b K pk−2 )dx, pk−2 ∈ Pk−2 (K ). (38c)
K

Each vector field qh ∈ BDMk (K ) is uniquely determined by (38a)–(38c), [11]. An immediate consequence
is the following result:
Lemma 3.5 For each qh ∈ BDMk (K ) it holds
⎛ ⎧ ⎫
  ⎨ ⎬
q2K (x)dx ≤ c ⎝h K (q K · n∂ K )2 ds + h 2K max (q K · ∇ p)2 dx; p ∈ Pk−1 , max| p(x)| ≤ 1
⎩ x∈K ⎭
K ∂K K
⎧ ⎫⎞
⎨ ⎬
+h 2K max (q K · curl(b K p))2 dx; p ∈ Pk−2 , max| p(x)| ≤ 1 ⎠ . (39)
⎩ x∈K ⎭
K

Proof The proof of (39) follows by a standard scaling argument. 


eq
Denoting by p̂h ∈ Vh the flux associated with û h , the flux ∈ BDMk (Ω, Th ) is specified element-wise
ph
according to
 
eq
n E · ph pk ds = n E · p̂h pk ds, pk ∈ Pk (E), E ∈ ∂ K , (40a)
E E
 
eq
ph · ∇ pk−1 dx = p̂h · ∇ pk−1 dx, pk−1 ∈ Pk−1 (K ), (40b)
K K
 
eq
ph · curl(b K pk−2 )dx = a∇ û h · curl(b K pk−2 )dx, pk−2 ∈ Pk−2 (K ). (40c)
K K

eq
Lemma 3.6 The flux ph ∈ BDMk (Ω, Th ) given by (40a)-(40c) is an equilibrated flux, i.e., it satisfies
(24a),(24b).
Proof It follows from Gauss’s theorem and (38b) that for pk−1 ∈ Pk−1 (K ) it holds
  
eq eq eq
− ∇ · ph pk−1 dx = ph · ∇ pk−1 dx − n∂ K · ph pk−1 ds
K K ∂K
  
= p̂h · ∇ pk−1 dx − n∂ K · p̂∂ K pk−1 ds = f h pk−1 dx. (41)
K ∂K K

eq
Since ∇ · ph and f h belong to Pk−1 (K ), from (41) we deduce (24a). For E ∈ EhN , we choose pk−1 ∈ Pk−1 (K )
eq
with supp( pk−1 ) = E. By the definition of the numerical flux p̂∂ K and the fact that both n E · ph | E and u N ,h | E
belong to Pk−1 (E), (24b) also follows from (41). 


3.5 Efficiency of the equilibrated error estimator

For the symmetric IPDG approximation of second-order diffusion problems on simplicial triangulations,
residual-type error estimators have been derived and analysed by many authors, [5,17,19,20]. These type of
estimators are of the form,

(ηhr es )2 = (ηrKes )2 , (42)
K ∈T h

123
Arab. J. Math. (2022) 11:407–426 417

where the element residuals ηrKes are given by

  
(ηrKes )2 := h 2K  f h + ∇ · a∇ û h 20,K + h E [n E · a∇ û h ] E 20,E + α 2 h −1
E [û h ] E 0,E
2

E∈EhI
 
+ α 2 h −1
E u D,h − û h 0,E +
2
h E u N ,h − n E · a∇ û h 20,E .
E∈EhD E∈EhN

The efficiency estimate from [20] shows

ηhr es  u − û h 1,h,Ω . (43)

eq
The efficiency of the equilibrated a posteriori error estimator ηh follows from (43) and the following result.
eq
Lemma 3.7 Let ηh be given by (25) and let ηhr es be the residual-type a posteriori error estimator (42). Then
there holds

eq
ηh  ηhr es . (44)

eq eq
Proof To estimate η K ,1 , K ∈ Th , from above, we apply Lemma3.5 to ph − a∇ û h . From (21b) and (40a) it
follows that

eq
n E · (ph − a∇ û h ) = n E · (p̂h − a∇ û h ) =
⎧1
⎨ 2 [n E · a∇ û h ] E − αh −1
E [û h ] E , E ∈ EhI , (45)
−1
−αh E (û h − u D,h ), E ∈ EhD
⎩
u N ,h − n E · a∇ û h , E ∈ EhN

On the other hand, due to (24a) we have

eq
∇ · (ph − a∇ û h ) = f h + ∇ · a∇ û h in each K ∈ Th (46)

Finally, in view of (40c) for each K ∈ Th it holds


eq
(ph − a∇ û h ) · curl(b K pk−2 )dx = 0, pk−2 ∈ Pk−2 (K ). (47)
K

Using (45)-(47), from Lemma3.5 we obtain

eq
(η K ,1 )2  (ηrKes )2 . (48)

Moreover, using (23) we find


h −1
eq
(η K ,2 )2  E [û h ] E 0,E .
2
(49)
E∈Eh

In view of (26a)-(26d), the assertion follows from (48) and (49). 



123
418 Arab. J. Math. (2022) 11:407–426

4 Adaptive refinement, numerical results

4.1 Adaptive refinement

The a posteriori error estimation’s main purpose is to indicate the elements where the error is large. Then we
use this information to locally refine those elements, and repeat the finite element computation. The adaptive
approach consists of successive cycles of the steps

SOLVE ⇒ ESTIMATE ⇒ MARK ⇒ REFINE.

In the step SOLVE, we compute the solution û h ∈ Vh of the IPDG approximation. In the second step ESTI-
eq eq eq
MATE, we compute the local components η K ,i and η E,i , 1 ≤ i ≤ 2, of the error estimator ηh (cf. (26a)-(26d)).
To come up with a purely element-based marking strategy, we set
 
eq

2
eq
 eq
ηK := η K ,i + η E,i , K ∈ Th ,
i=1 E∈∂ K

and use standard Dörfler marking [15] in step MARK: Given some constant 0 < θ ≤ 1, we choose a set
M ⊂ Th of elements K ∈ Th such that
eq
 eq
θ ηh ≤ ηK . (50)
K ∈M

In the final step REFINE, marked elements K ∈ M are subdivided by newest vertex bisection. The performance
of the equilibrated a posteriori error estimator will be illustrated in the sequel by numerical examples.

4.2 Numerical examples

We present now numerical examples which validate the optimality and robustness of the equilibrated a posteriori
error estimator based on the hypercircle method. We begin with verifying the optimal behaviour of the a
posteriori estimator by computing experimental order of convergence by means of manufactured solutions.
The efficiency of the a posteriori estimator is evaluated for the classical L-shaped domain problem and in
the case of highly discontinuous coefficients. The numerical results were obtained using the computational
framework FEniCS, [3,21] and the references within.

4.2.1 Validation, optimality of the estimator

We first validate the optimality of the a posteriori estimator by computing convergence rates in the following
case. We consider (1a)–(1c) in Ω = (−8, 8)2 such that the boundary consists of the Dirichlet boundary
Γ D = ∂Ω. For simplicity, we take constant coefficient a(x) ≡ 1, ∀x ∈ Ω and we select f, u D , and u N such
that
1
u 1 (x1 , x2 ) = 10e− 2 (x1 +x2 ) (x12 + x22 ),
2 2

is the exact solution in Ω.

Table 1 The errors and EOC’s for u 1 in L 2 and IPDG norms

Level Cell diam. Num. of cells L 2 -Error EOC (L 2 ) u − û h 1,h,Ω EOC (1, h) ηh EOC (ηh ) EI
0 2 256 4.79e−1 – 3.05 – 13.38 – 4.38
1 1 1024 2.86e−2 4.07 4.03e−1 2.92 2.06 2.70 5.12
2 0.5 4096 1.85e−3 3.95 5.33e−2 2.92 2.76e−1 2.90 5.17
The optimality of the estimator and the effectivity index are also shown

123
Arab. J. Math. (2022) 11:407–426 419

Fig. 1 The mesh and numerical solution contour lines for u 1 . Cubic elements were used for the approximation after each cycle
of refinement of the mesh with θ = 1.0 in the Dörfler marking

For the numerical approximation we use cubic elements for u 1 and we start the realizations with an initial
grid and then perform a series of uniform refinements, calculating the error in the L 2 and IPDG norms. The
experimental order of convergence(EOC) is computed as
   
E1 h1
E OC = log log ,
E2 h2
where E 1 , E 2 and h 1 , h 2 are the errors and maximum cell sizes, respectively, for two different realizations. The
numerical results for u 1 are shown in Table 1. For u 1 the expected convergence rates of 4 and 3 are observed
for the L 2 -norm and I P DG norm respectively. The optimality of the a posteriori estimator ηh is also depicted

123
420 Arab. J. Math. (2022) 11:407–426

Table 2 Adaptive refinement with θ = 0.5 in the Dörfler marking for u 1

Level Num. of cells u − û h 1,h,Ω ηh EI

0 256 3.05e+0 1.34e+1 4.38
1 310 7.52e−1 3.80e+0 5.05
2 544 1.75e−1 8.68e−1 4.96
3 1094 4.92e−2 2.50e−1 5.08
4 2196 1.58e−2 8.15e−2 5.15
5 3872 6.18e−3 3.13e−2 5.06
6 7344 2.36e−3 1.22e−2 5.16
The error in IPDG norm and the corresponding effectivity index are shown

Fig. 2 The mesh and numerical solution contour lines for u 1 . Cubic elements were used for the approximation after each cycle
of refinement of the mesh with θ = 0.5 in the Dörfler marking

123
Arab. J. Math. (2022) 11:407–426 421

ηh
in Table 1. The effectivity index, E I = which is also shown in the last column of the table,
u − û h 1,h,Ω
tends to a value of ∼ 5.
Finally, in Fig. 1 the approximate solution is depicted for the three levels of uniform mesh refinement.

4.2.2 Adaptivity and estimator performance

Next we investigate the performance of the a posteriori estimator in terms of its ability to detect regions of the
solutions where high resolution is required to maintain the overall accuracy of the method.
We consider the aforementioned example u 1 . We performed adaptive refinement of the mesh with θ = 0.5
in the Dörfler marking, and calculated the error in the IPDG norm and the value of the estimator. The results
are shown in Table 2. As in the uniform refinement, the effectivity index stabilizes around the same value
E I ∼ 5. In Fig. 2 the approximate solution is depicted for the three levels of mesh refinement.

4.2.3 Discontinuous coefficients and adaptivity

Next, we study the robustness of the estimator in the case where the coefficient a in (1a)–(1c) is discontinuous in
the domain. Let Ω = (0, 1)2 be such that Ω = ∪i=1 4 Ω , where Ω = {(x , x )|0 < x < 0.5, 0 < x < 0.5},
i 1 1 2 1 2
Ω2 = {(x1 , x2 )|0.5 < x1 < 1, 0 < x2 < 0.5} Ω3 = {(x1 , x2 )|0.5 < x1 < 1, 0.5 < x2 < 1}, and
Ω4 = {(x1 , x2 )|0 < x1 < 0.5, 0.5 < x2 < 1}. The discontinuity of a is along the two lines x = 0.5 and
y = 0.5, that is we define: a1 = a|Ω1 = 100, a2 = a|Ω2 = 1, a3 = a|Ω3 = 100, a4 = a|Ω4 = 1. The
domain Ω with its subdivision and the coefficient a are shown in Fig. 3. The boundary consists of the Dirichlet
boundary Γ D = ∂Ω. We also choose f, u D , and u N such that u 2 is the exact solution of (1a)–(1c).
u 2 (x1 , x2 ) = sin(π x1 )sin(π x2 ).
We use linear elements and we performed several levels of adaptive refinement of the mesh with θ = 0.5
in the Dörfler marking, and calculated the error in the IPDG norm, the value of the estimator and the effectivity
index. The results are shown in Table 3. The effectivity index reaches a stable value of E I ∼ 19. In Fig. 4 the
approximate solution is depicted for the three levels of mesh refinement.
In the second example with discontinuous coefficients we consider again Ω = (0, 1)2 such that Ω =
∪i=1 Ω i , where Ω1 = {(x1 , x2 )|x1 < x2 < 1−x1 }, Ω2 = {(x1 , x2 )|x2 < x1 < 1−x2 }, Ω3 = {(x1 , x2 )|1−x1 <
4

Fig. 3 Subdivision of Ω and distribution of a

Table 3 Discontinuous coefficient a for u 2 : the error in IPDG norm, the value of the estimator, and the effective index E I over
several levels of adaptive refinement of the mesh with θ = 0.5 in the Dörfler marking

Level Num. of cells u − û h 1,h,Ω ηh EI

0 256 3.68 102.16 27.78
1 438 3.09 78.63 25.47
2 860 2.77 52.23 18.87
3 1766 2.11 38.80 18.39
4 3470 1.43 25.62 17.87
5 7112 1.04 19.02 18.33
6 13,402 7.34e−1 13.66 18.61

123
422 Arab. J. Math. (2022) 11:407–426

Fig. 4 Discontinuous coefficient a for u 2 : adaptive mesh refinement and the corresponding numerical solution u h is shown here
for linear elements after each cycle of adaptive refinement of the mesh with θ = 0.5 in the Dörfler marking

Fig. 5 Subdivision of Ω and distribution of a

123
Arab. J. Math. (2022) 11:407–426 423

Table 4 Discontinuous coefficient a for u 3 : the error in IPDG norm, the value of the estimator, and the effective index E I over
several levels of adaptive refinement of the mesh with θ = 0.5 in the Dörfler marking

Level Num. of cells u − û h 1,h,Ω ηh EI

0 256 2.10 28.64 13.67
1 361 2.06 17.42 8.46
2 725 1.98 13.65 6.90
3 1569 1.58 12.57 7.95
4 3053 1.13 10.53 9.36
5 5923 0.75 9.07 12.16
6 11,365 0.51 7.44 14.64

Fig. 6 Discontinuous coefficient a for u 3 : adaptive mesh refinement and the corresponding numerical solution u h is shown here
for linear elements after each cycle of adaptive refinement of the mesh with θ = 0.5 in the Dörfler marking

123
424 Arab. J. Math. (2022) 11:407–426

Fig. 7 u 2 : the adaptive mesh refinement is shown here for k = 2 (right) and k = 3 (left) after 12 adaptive refinement cycles, with
θ = 0.5 in the Dörfler marking

Fig. 8 u 3 : the adaptive mesh refinement is shown here for k = 2 (right) and k = 3 (left) after 12 adaptive refinement cycles, with
θ = 0.5 in the Dörfler marking

x2 < x1 }, and Ω4 = {(x1 , x2 )|1 − x2 < x1 < x2 }. In this case the discontinuity of a is along the two main
diagonals of the domain, cf. Fig. 5.
We choose a1 = a|Ω1 = 1, a2 = a|Ω2 = 100, a3 = a|Ω3 = 1, a4 = a|Ω4 = 100. The boundary Γ = ∂Ω
consists of the Dirichlet boundary Γ D = ∂Ω. We also choose f, u D , and u N such that
u 3 (x1 , x2 ) = 64x12 (1.0 − x1 )2 x22 (1.0 − x2 )2 ,
is the exact solution of (1a)-(1c).
We use linear elements and we performed several levels of adaptive refinement of the mesh with θ = 0.5
in the Dörfler marking, and calculated the error in the IPDG norm, the value of the estimator and the effectivity
index. The results are shown in Table 4. The effectivity index reaches a stable value of E I ∼ 18. In Fig. 6 the
approximate solution is depicted for the three levels of mesh refinement.
The last two examples show how the estimator helps in marking and refining the cells around the interfaces
between the domains where the value of a varies, cf. Figs. 3 and 5. The effectivity indices shows clearly the
efficiency and reliability of the a posteriori error estimator. The refinement improves with the increase in the
polynomial degree k as shown below in Figs. 7 and 8.

123
Arab. J. Math. (2022) 11:407–426 425

5 Conclusions

Based on the hypercircle method also known as the Prager–Synge theorem, we have developed and analyzed
an equilibrated a posteriori error estimator for IPDG approximations of second order diffusion problems. The
error estimator relies on the construction of an equilibrated flux and has been shown to be both efficient and
reliable. Numerical results for some selected examples confirm the theoretical findings.

Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use,
sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original
author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other
third party material in this article are included in the article’s Creative Commons licence, unless indicated otherwise in a credit
line to the material. If material is not included in the article’s Creative Commons licence and your intended use is not permitted
by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To
view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.

Funding No funding was received to assist with the preparation of this manuscript.

Declarations

Conflict of interest The authors have no relevant financial or non-financial interests to disclose.

Data availability statement Not applicable.

References

1. Ainsworth, M.: A posteriori error estimation for discontinuous Galerkin finite element approximation. SIAM J. Numer. Anal.
45(4), 1777–1798 (2007)
2. Ainsworth, M.; Oden, J.T.: A Posteriori Error Estimation in Finite Element Analysis. Wiley, New York (2000)
3. Alnæs, M.S.; Blechta, J.; Hake, J.; Johansson, A.; Kehlet, B.; Logg, A.; Richardson, C.; Ring, J.; Rognes, M.E.; Wells, G.N.:
The fenics project version 1.5. Arch. Numer. Softw. 3(100),9–13 (2015)
4. Arnold, D.N.; Brezzi, F.; Cockburn, B.; Marini, L.D.: Unified analysis of discontinuous Galerkin methods for elliptic
problems. SIAM J. Numer. Anal. 39(5), 1749–1779 (2002)
5. Bonito, A.; Nochetto, R.H.: Quasi-optimal convergence rate of an adaptive discontinuous Galerkin method. SIAM J. Numer.
Anal. 48(2), 734–771 (2010)
6. Braess, D.: Finite Elements: Theory, Fast Solvers, and Applications in Solid Mechanics. Cambridge University Press,
Cambridge (2007)
7. Braess, D.; Fraunholz, T.; Hoppe, R.H.W.: An equilibrated a posteriori error estimator for the interior penalty discontinuous
Galerkin method. SIAM J. Numer. Anal. 52(4), 2121–2136 (2014)
8. Braess, D.; Hoppe, R.H.; Schöberl, J.: A posteriori estimators for obstacle problems by the hypercircle method. Comput.
Vis. Sci. 11(4–6), 351–362 (2008)
9. Brenner, S.C.; Scott, R.: The Mathematical Theory of Finite Element Methods, 3rd edn Springer, Berlin (2008)
10. Brezzi, F.; Douglas, A.; Marini, D.: Two families of mixed finite elements for second order elliptic problems. Numer. Math.
47(1), 217–235 (1985)
11. Brezzi, F.; Fortin, M.: Mixed and Hybrid Finite Element Methods. Springer, Berlin (1991)
12. Cochez-Dhondt, S.; Nicaise, S.: Equilibrated error estimators for discontinuous Galerkin methods. Numer. Methods Partial
Differ. Equ. 24(5), 1236–1252 (2008)
13. Cockburn, B.: Discontinuous Galerkin Methods for Computational Fluid Dynamics, pp. 1–63. Wiley, New York (2017)
14. DiPietro, D.A.; Ern, A.: Mathematical Aspects of Discontinuous Galerkin Methods. Springer, Berlin (2012)
15. Dörfler, W.: A convergent adaptive algorithm for Poisson’s equation. SIAM J. Numer. Anal. 33(3), 1106–1124 (1996)
16. Ern, A.; Stephansen, A.F.; Vohralík, M.: Guaranteed and robust discontinuous Galerkin a posteriori error estimates for
convection–diffusion–reaction problems. J. Comput. Appl. Math. 234(1), 114–130 (2010)
17. Hoppe, R.H.; Kanschat, G.; Warburton, T.: Convergence analysis of an adaptive interior penalty discontinuous Galerkin
method. SIAM J. Numer. Anal. 47(1), 534–550 (2008)
18. Houston, P., Schötzau, D., Wihler, T.P. Mixed hp-Discontinuous Galerkin Finite Element Methods for the Stokes Problem in
Polygons. In: Feistauer, M., Dolejší, V., Knobloch, P., Najzar, K. (eds) Numerical Mathematics and Advanced Applications.,
493 –501 (2004) Springer, Berlin, Heidelberg.
19. Karakashian, O.A.; Pascal, F.: A posteriori error estimates for a discontinuous Galerkin approximation of second-order
elliptic problems. SIAM J. Numer. Anal. 41(6), 2374–2399 (2003)
20. Karakashian, O.A.; Pascal, F.: Convergence of adaptive discontinuous Galerkin approximations of second-order elliptic
problems. SIAM J. Numer. Anal. 45(2), 641–665 (2007)
21. Logg, A.; Mardal, K.-A.; Wells, G.N.; et al.: Automated Solution of Differential Equations by the Finite Element Method.
Springer, Berlin (2012)
22. Rivière, B.: Discontinuous Galerkin Methods for Solving Elliptic and Parabolic Equations: Theory and Implementation.
Society for Industrial and Applied Mathematics, Philadelphia (2008)

123
426 Arab. J. Math. (2022) 11:407–426

23. Schötzau, D.; Schwab, C.; Toselli, A.: Mixed hp-dgfem for incompressible flows. SIAM J. Numer. Anal. 40, 2171–2194
(2002)
24. Shahbazi, K.: An explicit expression for the penalty parameter of the interior penalty method. J. Comput. Phys. 205(2),
401–407 (2005)
25. Synge, J.: The method of the hypercircle in function-space for boundary-value problems. Proc. R. Soc. Lond. A Math. Phys.
Eng. Sci. 191, 447–467 (1947)
26. Synge, J.L.: The Hypercircle in Mathematical Physics. Cambridge University Press, Cambridge (1957)
27. Vejchodskỳ, T.: Local a posteriori error estimator based on the hypercircle method. In: Neittaanmäki, P., et al. (eds.) Pro-
ceedings ECCOMAS 2004. University of Jyvaskylä, Jyvaskylä (2004)
28. Warburton, T.; Hesthaven, J.S.: On the constants in hp-finite element trace inverse inequalities. Comput. Methods Appl.
Mech. Eng. 192(25), 2765–2773 (2003)

Publisher’s Note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional
affiliations.

123