Coupling non-conforming discretizations of PDEs by spectral
approximation of the Lagrange multiplier space
Simone Deparisa , Luca Pegolottia,∗
a Institute

of Mathematics, École Polytechnique Fédérale de Lausanne, Station 8, EPFL, CH–1015 Lausanne, Switzerland

arXiv:1802.07601v1 [math.NA] 21 Feb 2018

Abstract
This work focuses on the development of a non-conforming domain decomposition method for the approximation of PDEs based on weakly imposed transmission conditions: the continuity of the global
solution is enforced by a discrete number of Lagrange multipliers defined over the interfaces of adjacent
subdomains. The method falls into the class of primal hybrid methods, which also include the well-known
mortar method. Differently from the mortar method, we discretize the space of basis functions on the
interface by spectral approximation independently of the discretization of the two adjacent domains;
one of the possible choices is to approximate the interface variational space by Fourier basis functions.
As we show in the numerical simulations, our approach is well-suited for the solution of problems with
non-conforming meshes or with finite element basis functions with different polynomial degrees in each
subdomain. Another application of the method that still needs to be investigated is the coupling of solutions obtained from otherwise incompatible methods, such as the finite element method, the spectral
element method or isogeometric analysis.
Keywords: Partial Differential Equations, Non-conforming method, Domain decomposition

1. Introduction
In numerical analysis, domain decomposition methods are techniques for the splitting of Partial Differential Equations (PDEs) into smaller and coupled problems defined over subsets of the original domain.
The splitting may be motivated by physical reasons, for instance when the subdomains are characterized
by different governing equations (e.g. in fluid-structure-interaction problems [1]) or by discretization needs,
should it be required to employ specific methods (e.g. finite element method or spectral element method)
or specific polynomial degrees in certain regions of the domain [2]. Moreover, domain decomposition
methods have become particularly important for the solution of large scale problems on multiprocessors
or clusters, as they allow the mapping of the subproblems on separate cores [3].
Domain decomposition methods are typically based either on iterative or direct procedures [4]. In
the first class of techniques the continuity on the interfaces of the solution, of its normal derivatives or
combinations of the two are strongly imposed. Typically, these methods require solving the problems
defined on the subdomains separately multiple times while imposing artificial boundary conditions based
on the solutions at the previous iteration. The type of boundary conditions employed on each subdomain
is a peculiarity of each algorithm, so that the literature on the topic commonly refers to the DirichletDirichlet algorithm, the Dirichlet-Neumann algorithm, and so on; see e.g [2] for details. These strategies
allow reducing the size of the linear systems to be solved and, most importantly, to compute the solution
on each subdomain in parallel.
In this paper, we present an approach belonging to the class of direct procedures in which the continuity
conditions (often called transmission conditions) are weakly imposed through the use of suitable Lagrange
multipliers. Our method is applied to PDEs written in primal hybrid formulation, and for this reason
it shares some of the features of the well-known mortar method [5, 6]. This was originally proposed to
solve PDEs by combining spectral elements and finite elements, or by combining finite element spaces
with different polynomial degrees, in non overlapping portions of the domain [7]. Since then, the mortar
method has become the non-conforming method of choice in many areas of computational science and
∗ Corresponding

author
Email address: luca.pegolotti@epfl.ch (Luca Pegolotti)

Preprint submitted to arXiv

March 7, 2018

engineering, for example in contact mechanics [8], solid mechanics [9], fluid mechanics [10] and fluidstructure interaction problems [11]; see also [12, 13, 14]. The implementation of the mortar method is not
straightforward, as the algorithm is based on L2 -projections of the traces of functional spaces defined on
a group subdomains – the masters – onto the interfaces of the adjacent ones – the slaves. INTERNODES
[15, 16], a recently developed method for the treatment of non-conforming meshes, overcomes this issue
by treating the transmission conditions with the interpolation of basis functions of the master domains
onto the interfaces of the slaves.
As in the mortar method, our approach is based on the idea that the global problem can be subdivided
into a set of smaller problems coupled with weak conditions relying on basis functions defined on the
interfaces. In the mortar method, such basis functions are obtained from the trace space of the adjacent
slave domains. This choice is convenient from the analysis standpoint but makes the implementation of
the method cumbersome. Another drawback is that the final solution is dependent on the choice of master
and slave domains. The originality of our method is to consider basis functions on the interfaces which
are completely independent of the discretization of the neighboring domains: in this paper, we employ
spectral basis functions (specifically, Fourier basis functions). This comes with the advantage of obtaining
a solution which is indifferent to the choice of master and slave domains. Moreover, the accuracy of the
coupling of solutions at the interfaces is easily tuned by varying the number of basis functions on the
common boundary. Our approach can be interpreted as a specialization of the three-fields method [17],
where the space of the three Lagrange multipliers used to weakly impose the continuity of the solution
is (a priori) independent of the spaces defined on the adjacent domains. As the functional spaces in the
subdomains are mutually independent, our choice of basis functions is well-suited for the coupling of
solutions obtained on non-conforming (at the interfaces) meshes, with finite element spaces with different
polynomial degrees, or with different numerical methods, e.g. finite element method, spectral element
method, or isogeometric analysis [18, 19].
The paper is structured as follows. In Section 2, we present the method on an elliptic problem defined
over a domain partitioned into two regions. Section 3 focuses on the discretization of the weak formulation
derived in Section 2. In Section 4, we briefly compare our method with other non-conforming methods,
namely the mortar method, INTERNODES, and the three-field method, and focus on the similarities
and peculiarities with respect to our approach. In Section 5 we address the matter of the stability
of the method, which is strictly related to the inf-sup condition. In Section 6, the method is used to
solve two-dimensional benchmark problems with finite element discretizations in the subdomains: the
Poisson problem on two subdomains (Section 6.2) and the Navier-Stokes equations on five subdomains
(Section 6.3). Finally, in Section 7 some conclusions are drawn.
1.1. Notation
The notation adopted in this paper is standard and commonly found in the literature; see e.g. [20].
Given a generic open and bounded domain Ω embedded in Rd , we define, for all ϕ, ψ : Ω → R and all
ϕ, ψ : Ω → Rd
Z
Z
(ϕ, ψ)Ω :=

ϕψ dx,

(φ, ψ)Ω :=

Ω

ϕ · ψ dx,
Ω

and consider the following Hilbert spaces
L2 (Ω) := {ϕ : Ω 7→ R : (ϕ, ϕ)Ω < ∞},
H 1 (Ω) := {ϕ ∈ L2 (Ω) : ∇ϕ ∈ [L2 (Ω)]d },
H(div; Ω) := {φ ∈ [L2 (Ω)]d : divφ ∈ L2 (Ω)},
with the associated norms
kϕk2L2 (Ω) := (ϕ, ϕ)Ω ,
kϕk2H 1 (Ω) := (ϕ, ϕ)Ω + (∇ϕ, ∇ϕ)Ω ,
kφk2H(div;Ω) := (φ, φ)Ω + (divφ, divφ)Ω .
Given a measurable set Σ ⊆ ∂Ω (where ∂Ω denotes the boundary of Ω), we also define
HΣ1 (Ω) := {ϕ ∈ H 1 (Ω) : ϕ = 0 on Σ}.

2

Under the assumption of sufficient regularity of Σ, there exists a unique linear and continuous application γΣ : H 1 (Ω) 7→ L2 (Σ) called trace operator [21, 22] such that γΣ ϕ = ϕ|Σ for all ϕ ∈ H 1 (Ω) ∩ C 0 (Ω̄),
having indicated with C 0 (Ω̄) the space of continuous functions over the closure of Ω. The range of such
operator is denoted H 1/2 (Σ) ⊂ L2 (Σ). We recall that
kηkH 1/2 (Σ) :=

inf

ϕ∈H 1 (Ω)
ϕ|Σ =η

kϕkH 1 (Ω)

is a norm for H 1/2 (Σ) [23].
For each linear Hilbert space H, we denote H0 the space of linear and bounded functionals on H, namely
its dual space. In particular, we adopt the notation H −1 (Ω) := (H 1 (Ω))0 and H −1/2 (Σ) := (H 1/2 (Σ))0 .
The action of an element of the dual space ξ ∈ H0 on an element of the Hilbert space ϕ ∈ H is indicated
hξ, ϕiH0 or simply hξ, ϕi whenever ambiguity does not arise. Furthermore, we will simply indicate hξ, ϕiΣ
the duality in H −1/2 (Σ). We note that, for ϕ ∈ H 1 (Ω), we will adopt the abuse of notation hξ, ϕiΣ with
ξ ∈ H −1/2 (Σ) to indicate the duality of ξ with the trace of ϕ on Σ. Even though not mathematically
rigorous, this notation is commonly used in the literature.
2. Theory of primal hybrid methods
In this section, we recall the theory of primal hybrid methods for the solution of Partial Differential
Equations (PDEs). These approaches are based on the primal hybrid principle [24], according to which
the continuity across subdomains is weakened by means of Lagrange multipliers. We refer the reader to
[25, 26, 27] for the theory of primal hybrid methods. Here, we recall the main ideas by following the
presentation in [23]. We also restrict ourselves to only two partitions of the domain; however – as we shall
see in Section 3.1 – the method extends to an arbitrary number of partitions.
We are interested in solving a generic PDE described by a second order elliptic operator on an open
and bounded domain Ω with homogeneous Dirichlet boundary conditions on ∂Ω. Specifically, we assume
that a(ϕ, ψ) for ϕ, ψ ∈ H 1 (Ω) is the bilinear form corresponding to the elliptic operator and f is a given
forcing term; we consider problems whose weak formulation can be written as:
1
(Ω), such that
(W1) given f ∈ H −1 (Ω), find u ∈ H∂Ω
1
∀v ∈ H∂Ω
(Ω).

a(u, v) = hf, vi

(1)

In the sequel, we will often use the Poisson problem with homogeneous boundary conditions
−∆u = f
u=0

in Ω,
on ∂Ω,

(2)

as representative of this class of problems. In this specific case, a(ϕ, ψ) = (∇ϕ, ∇ψ)Ω .
Let us assume that the domain Ω can be partitioned into two non-overlapping open and bounded
domains, such that Ω = Ω1 ∪ Ω2 and Ω1 ∩ Ω2 = ∅; we denote Γ the interface between the two domains,
i.e. Γ = Ω1 ∩ Ω2 . Our goal is to solve, rather than the global problem W1, two local and coupled problems
defined on the partitions Ωi , such that the global solution can be constructed by combining the solutions
1
of the local problems. To this end, let us introduce the functional spaces X (i) = H∂Ω∩∂Ω
(Ωi ) and
i
X := {ϕ ∈ L2 (Ω) : ϕ|Ωi ∈ X (i) for i = 1, 2},
which is a Hilbert space when endowed with the (broken) norm
kϕk2X

:=

2
X

kϕ|Ωi k2H 1 (Ωi ) .

i=1
1
The space H∂Ω
(Ω) is characterized as a subspace of X under suitable conditions [23] which we will state
in the following Lemma and motivate in its proof. We remark that an analogous result is presented in
[25].

3

Proposition 1. [25, cf. Proposition 2.1.1]
1
e := {ϕ ∈ X :
H∂Ω
(Ω) ≡ V

2
X
hφ · ni , ϕi∂Ωi = 0

∀φ ∈ H(div; Ω)},

i=1

where ni is the outward unit vector normal to ∂Ωi .
Proof. Firstly, let us recall that for any φ ∈ H 1 (div; Ω) and for i = 1, 2, Green’s formula gives [25, Lemma
2.1.1]
Z
(∇ϕ · φ + ϕdivφ) dx = hφ · ni , ϕi∂Ωi

∀ϕ ∈ X .

(3)

Ωi

Eq. (3) implies that, since the left hand side is bounded, the duality of φ · ni ∈ H −1/2 (∂Ωi ) on the trace
of ϕ|Ωi is well-defined.
e ⊂ H 1 (Ω), because for all ϕ ∈ V
e ⊂ L2 (Ω)
Clearly V
∂Ω
Z

2

|∇ϕ| dx =
Ω

2 Z
X

|∇ϕ|2 dx < ∞,

Ωi

i=1

1
and ϕ = 0 on ∂Ω in H 1/2 (∂Ω). Let us show the other inclusion. For every ϕ ∈ H∂Ω
(Ω), it holds that
1
ϕ ∈ L2 (Ω) and ϕ|Ωi ∈ X (i) for i = 1, 2, which implies that H∂Ω
(Ω) ⊂ X . Moreover, by applying Green’s
formula as in Eq. (3), we find for all φ ∈ H(div; Ω)
2
2 Z
X
X
hφ · ni , ϕi∂Ωi =
i=1

i=1

Z
(∇ϕ · φ + ϕdivφ) dx =

Ωi

(∇ϕ · φ + ϕdivφ) dx = hφ · n, ϕi∂Ω = 0,

(4)

Ω

where the last equality comes from the fact that ϕ has null trace on the boundary ∂Ω. Since Eq. (4)
1
e it must be H 1 (Ω) ≡ V.
e
shows that H∂Ω
(Ω) ⊂ V,
∂Ω
Remark 1. If the PDE is equipped with Dirichlet conditions on ∂ΩD and Neumann conditions on ∂ΩN ,
the space X must be defined such that ϕ|Ωi belongs to X (i) for i = 1, 2. In this case it is not sufficient
1
to ask that φ ∈ H(div; Ω) in the definition of V to have the equivalence between H∂Ω
(Ω) and V. In
particular, the space H(div; Ω) must be restricted to functions φ such that hφ · n, ϕi∂Ω = 0 for each
ϕ ∈ V.
P2
The condition i=1 hφ · ni , ϕi∂Ωi = 0 for each φ ∈ H(div; Ω) is global, in the sense that it involves the
trace of ϕ|Ωi on the whole ∂Ωi , even though it essentially constrains the restrictions of ϕ to Ω1 and Ω2
to have the same trace at the common interface Γ. Unfortunately, splitting the dualities into two parts
corresponding to Γ and ∂Ω \ Γ is not allowed, as the restrictions of the traces to portions of ∂Ωi can lead
to unbounded dualities. To overcome this issue, we introduce
1/2

(i)

H00 (Γi ) := {η ∈ H 1/2 (Γi ) : E0 η ∈ H 1/2 (∂Ωi )},
with norm

(i)
kηkH 1/2 (Γi ) := kE0 ηkH 1/2 (∂Ωi ) ,
00

(i)
E0 η

where
is the trivial extension by zero of η to the whole boundary of ∂Ωi and Γi = ∂Ωi \ ∂Ω. In the
1/2
1/2
1/2
following, we will consider H00 (Γ) := H00 (Γ1 ) ∩ H00 (Γ2 ). Let us define the spaces
1/2
X00 := {ϕ ∈ X : [ϕ]Γ ∈ H00 (Γ)},

where [ϕ]Γ denotes by our convention the difference of the traces of ϕ|Ω2 and ϕ|Ω1 on Γ, and
−1/2

Λ := H00

(Γ),

(5)

with norm
kηkΛ := kηkH −1/2 (Γ1 ) + kηkH −1/2 (Γ2 ) .
00

00

Furthermore, we introduce the bilinear form
b(ϕ, ξ) := hξ, [ϕ]Γ iΛ
4

1
for ϕ ∈ X00 and ξ ∈ Λ. It can be easily verified [23] that another characterization of H∂Ω
(Ω) analogous
to that in Proposition 1 is given by
1
H∂Ω
(Ω) ≡ V := {ϕ ∈ X00 : b(ϕ, ξ) = 0 ∀ξ ∈ Λ}.

(6)

1
In the sequel, we will use the letter V to refer to H∂Ω
(Ω).
We are now ready to state the primal hybrid formulation of the original weak formulation W1. We
remark that, whenever applied to functions of X , the bilinear form a(·, ·) is to be intended as the sum of
the bilinear forms restricted to the two subdomains.

(W2) given f ∈ H −1 (Ω), find u ∈ X00 and λ ∈ Λ such that
a(u, v) + b(v, λ) = hf, vi

∀v ∈ X00 ,

b(u, η) = 0

∀η ∈ Λ.

(7)

Proposition 2. If u ∈ V is a solution of W1 and there exists λ ∈ Λ such that
b(v, λ) = hf, vi − a(u, v)

∀v ∈ X00 ,

(8)

then (u, λ) ∈ X00 × Λ is a solution of W2. On the other hand, if (u, λ) ∈ X00 × Λ is a solution of W2,
then u ∈ V and u is a solution of W1.
Proof. Let u ∈ V be a solution of W1, then u ∈ X00 and the second condition in Eq. (7) is satisfied
because of the definition (6). The first condition in Eq. (7) is satisfied when choosing λ ∈ Λ such that
Eq. (8) is verified. Conversely, if (u, λ) ∈ X00 × Λ is a solution for W2, then u ∈ V because of the second
condition in Eq. (7). Moreover, for each v ∈ V, b(v, ξ) = 0 for all ξ ∈ Λ and, in particular, for ξ = λ, and
the first condition in Eq. (7) becomes Eq. (1).
Remark 2. If we consider the Poisson equation (2), then Eq. (8) is verified by taking φ = −∇u and by
choosing λ ∈ Λ such that φ · n1 = λ, n1 being the outward unit vector normal to ∂Ω1 . Indeed, by using
integration by parts we find for all v ∈ X00
b(v, λ) = hφ · n1 , [v]Γ iΛ = −

2 Z
X
i=1

=

Z
2
X
i=1

Ωi

∇u · ni v ds

Γ

!

Z
f v dx −

∇u · ∇v dx

= hf, vi − a(u, v),

Ωi

where we used the fact that n1 = −n2 . Note that, if we defined the jump across the interface of a function
ϕ ∈ X00 as the difference of the traces on Γ of ϕ|Ω1 and ϕ|Ω2 , then φ · n2 = λ. Hence, the Lagrange
multiplier in Eq. (7) plays the role of the normal derivative of u at the interface Γ [28], with the direction
of the normal at the interface being determined by the definition of the jump.
3. Discretization of the primal hybrid formulation
We now consider the discretization of the weak formulation W2. We take two arbitrary finite dimensional functional spaces X h,(1) ⊂ X (1) and X h,(2) ⊂ X (2) spanned by two sets of basis functions
(1)
(1)
(2)
(2)
ϕi ∈ X (1) (with i = 1, . . . , nbf ) and ϕi ∈ X (2) (with i = 1, . . . , nbf ) respectively. We assume that
h,(1)
h,(2)
functions in X
and X
can be trivially extended by zero in the other domain and that such extension belong to X00 . The discrete version of the global space X00 is consequently obtained by considering
(1)
(2)
the space X h ⊂ X00 of dimension dim(X h ) = nbf = nbf + nbf and spanned by the basis functions
(1) n

(1)

(2) n

(2)

bf
bf
bf
{ϕi }ni=1
= {ϕi }i=1
∪ {ϕi }i=1
.
Pnbf
ui ϕi . In the numerical applications in Section 6,
The solution can be then approximated as u ≈ uh = i=1
we will consider standard finite element Lagrangian basis functions built over suitable triangulations T h,(1)
and T h,(2) of Ω1 and Ω2 respectively for the discretization of X (1) and X (2) ; we will always assume that
such triangulations meet standard regularity requirements [20], but we do not require the conformity of

5

the global mesh T h = T h,(1) ∪T h,(2) . We define conforming meshes those meshes for which the intersection
of two elements is either null, a vertex or a whole edge; in non-conforming meshes, on the contrary, two
elements can also share portions of their edges. The discretization parameter h is generic and defines
a family of discretized spaces; when using finite elements, for example, h refers to the maximum edge
length of an element – often called mesh size – in the triangulations of Ω1 and Ω2 . More generally, h
could be also considered a characteristic of the single subdomain, since – as we already mentioned –
the discretizations in Ω1 and Ω2 are independent one of the other and could be obtained from different
discretization methods (e.g. finite elements for Ω1 and isogeometric analysis for Ω2 ).
δ
Our proposition is to discretize
Λ is
PnΓ Λ as Λ by using a set of basis functions ξi ∈ Λ, such that λ ∈h,(1)
δ
,
approximated as λ ≈ λ = i=1 λi ξi . We remark that we characterize the refinement levels for X
X h,(2) and Λδ with different discretization parameters h and δ: this is to indicate that the discretization
of Λ is indeed independent of the discretization on Ω1 and Ω2 . For instance, in the two-dimensional case,
a suitable choice would consist of choosing as ξi the basis functions associated to the low-frequencies of
the Fourier basis defined on the common interface Γ, and the accuracy of the discretization of Λδ can be
increased independently of h by adding Fourier basis functions to the set ξi . In the numerical simulations
of Section 6 we will follow this approach. Alternative possibilities for the discretization of the Lagrange
multiplier space include other spectral basis functions, such as e.g. Legendre or Chebyshev polynomials.
The discrete space for the approximation of V is then defined as
V h,δ := {ϕh ∈ X h : b(ϕh , ξ δ ) = 0 ∀ξ δ ∈ Λδ },
Remark 3. V h,δ is not a subspace of V. As a matter of fact, if Λδ is not equal to Λ, then there may exist
ξ ∈ Λ, ξ 6∈ Λδ such that b(ϕh , ξ) 6= 0 for some ϕh ∈ V h,δ , and therefore ϕh 6∈ V. If we replaced V by
V h,δ in W1, we would obtain a non-conforming method, i.e. a numerical method in which the discretized
search space is not contained into the continuous search space. The generalized version of Cea’s lemma for
this family of methods is Strang’s second lemma [29], which states that the solution uh of the discretized
version of W1 satisfies
!
|a(u, wh ) − hf, wh i|
h
h
,
(9)
ku − u kV h,δ ≤ C
inf ku − v kV h,δ + sup
kwh kV h,δ
v h ∈V h,δ
wh ∈V h,δ
where C > 0 and k · kV h,δ is a norm for V h,δ . Note that the consistency error – i.e. the second term of the
right hand side in Eq. (9) – is identically zero for each wh ∈ V h,δ if V h,δ ⊂ V because u is a solution of
W1.
The discretization of W2 is simply obtained by replacing the continuous functional spaces with their
discrete counterparts, namely:
(W3) given f ∈ H −1 (Ω), find uh ∈ X h and λδ ∈ Λδ such that
a(uh , v h ) + b(v h , λδ ) = hf, v h i

∀v h ∈ X h ,

b(uh , η δ ) = 0

∀η δ ∈ Λδ .

(10)

By expanding uh and λδ on their respective bases, Eq. (10) can be rewritten in system form as
"
#   
f
A BT u
=
,
(11)
λ
0
B
0
where Aij = a(ϕj , ϕi ), Bij = b(ϕj , ξi ), ui = ui , λ = λi and fi = hf, ϕi i. By arranging the basis functions
ϕn and the degrees of freedom such that all the basis functions corresponding to Ω1 come before those
of Ω2 , system (11) can be written as

   
A1
0 −B1T
u1
f1


A2 B2T  u2  = f2  ,
 0
λ
0
−B1 B2
0
R (1)
R (2)
where (B1 )ij = Γ ϕj ξi ds and (B2 )ij = Γ ϕj ξi ds are coupling matrices. Clearly, B1 and B2 are likely
(1)

to be sparse, as only the basis functions ϕi

(2)

and ϕi

6

not vanishing on Γ lead to non-zero integrals.

(1)

ξi

Ω3

Ω1

(2)

ξi



A1

 0

 0

−B11
0

Γ2
Γ1

Ω2

0
A2
0
B21
−B22

0
0
A3
B31
B32

T
−B11
T
B21
T
B31
0
0


0
T 
−B22

T 
B32


0 
0

Figure 1: On the left, example of a three-way partition of Ω with two interfaces; each interface corresponds to a Lagrange
multiplier space. On the right, corresponding system matrix obtained from the discretization of the spaces.

In this paper, the computation of the coupling matrices is performed by numerically integrating by
Gauss quadrature rules [30] the integrals. Let us consider for instance the case of B1 in the two dimensional
Sn(1)
(1)
(1)
el,Γ
case. The triangulation T h,(1) induces on Γ a partition into nel,Γ elements, i.e. Γ = i=1
Ei . Given a
Gauss quadrature rule of order 2q − 1, the approximation of each term of B1 is computed as
(1)

Z
(B1 )mn =
Γ

(1)

nel,Γ

ϕ(1)
n ξm

ds =

XZ
i=1

Γ

nel,Γ

ϕ(1)
n ξm

ds ≈

q
XX

gq
gq
|det(Ji )|ϕ(1)
n (φi (xj ))ξm (φi (xj ))ωj ,

(12)

i=1 j=1

(1)
(1)
where det(Ji ) is the determinant of the Jacobian of the map φi : Ei → (−1, 1) from Ei to the reference
gq
th
interval (−1, 1), xj is the j Gauss quadrature node in (−1, 1) and ωj is the associated weight. As it
is evident from Eq. (12), in order to compute the approximation of B1 it is sufficient to being able to
(1)
gq
evaluate the product ϕn (φi (xgq
j ))ξm (φi (xj )) at each quadrature node.

3.1. Generalization to multiple subdomains
In the previous sections we decided to limit ourselves to the case where the domain of the PDE is
partitioned into two subdomains. This choice is motivated mainly by the fact that considering the generic
case of multiple subdomains leads necessarily to complexity in the notation. We refer the reader to the
already mentioned references [23, 27] for examples of how the functional spaces we considered in Section 2
could be adapted to the case of multiple subdomains. One aspect that differentiates our approach from
other methods (such as the mortar method) is that, in the discretization process, our method requires
to define a set of basis functions for the Lagrange multiplier space of each interface. These bases can be
chosen independently one from the other.
As an example, Fig. 1 (left) shows a three-way partition of Ω into three domains with two interfaces. Each of the interfaces Γ1 and Γ2 requires the definition of a corresponding space for the Lagrange
multipliers. After the discretization, the matrix of the algebraic system can be written as displayed in
R
(i) (j)
Fig. 1 (right), where the matrices (Bij )mn = Γj ϕn ξm ds discretize the coupling between the ith domain
and the j th interface. We remark that the signs of the coupling matrices are determined by the definition
of the normals at each interface. Notice that, besides the choice of the orientation of the normals at the
interfaces, there is no hierarchy among the subdomains.
4. Relationship with other non-conforming methods
4.1. Relationship with the mortar method
The mortar method can be derived from the same problem written in primal hybrid formulation W2
we considered in Section 2 [28]. Its discretized weak formulation could be rewritten in the form of a
saddle-point problem similar to Eq.(7) in which the space of Lagrange multiplier ΛhM depends on the
discretization of either Ω1 or Ω2 ; see [31, 32]. In particular, the classic mortar method requires assigning
7

to Ω1 or to Ω2 the role of master and slave domains. The basis functions of ΛhM are chosen as the trace
of the basis functions defined over the triangulation of the slave domain which do not vanish on Γ; the
polynomial order of the basis functions on the extrema of Γ is usually decreased by one.
With respect to the mortar method, we believe that the main advantages of our approach are the
following
1. the accuracy of the coupling can be increased or decreased by varying the number of Lagrange
multipliers at the interfaces independently of the discretization in the subdomains;
2. the solution is independent of the partition of the subdomains into master and slaves;
3. the computation of the coupling matrices does not require projections between meshes, which makes
the implementation of the method easier.
One drawback of our method is that the Lagrange multiplier space has to be rich enough to provide the
necessary accuracy, but coarse enough to satisfy the inf-sup condition as described in Section 5.
4.2. Relationship with INTERNODES
The INTERNODES (INTERpolation for NOnconforming DEcompositionS) method [15, 16] is based
on an interpolation approach, rather than the L2 -projection approach which characterizes the mortar
method. Given each interface, the two adjacent subdomains are given the role of master and slave domains.
Similarly to the mortar method, the traces of the (finite element or spectral element) basis functions
defined over the meshes of the master and slave domains are used to enforce the continuity of the solution
and the normal stresses. More precisely, two interpolation operators – or intergrid operators – are defined:
the interpolation operator from the master to the slave domain is used to ensure the continuity of the
solution, while the interpolation operator from the slave to the master domain enforces the continuity of
the normal fluxes. INTERNODES has been proven to retain the optimal convergence properties of the
mortar method. For more information about the method and its analysis, we refer the reader to [33].
Being INTERNODES closely related to the mortar method, compared to the former our approach
offers the same advantages we presented in Section 4.1 except for 3. Indeed, INTERNODES has the big
advantage of being simple to implement and allowing for small geometric non-conformity. We believe that
our method is as simple to implement as INTERNODES and that it can be extended to non-conforming
geometries with the help of localized Rescaled Radial Basis Interpolation [16, 34]. One complexity of
INTERNODES comes from the special treatment of integrals at the intersection of the interface Γ with
portions of the boundary where non-homogeneous Neumann conditions are imposed. The method proposed here does not need such special treatment.
4.3. Relationship with the three-field method
The three-field method was originally proposed in [17] and analyzed in [35]. Compared to the mortar
method, it has had significantly less impact on the domain decomposition community.
The multidomain extension of the weak formulation W1 by the three fields method reads [7]:
(W4) for i = 1, 2, find u(i) ∈ X (i) , σ (i) ∈ H −1/2 (Γ) and ψ ∈ H 1/2 (Γ) such that
a(u(1) , v (1) ) − hσ (1) , v (1) iH −1/2 (Γ) = hf, v (1) i

∀v (1) ∈ X (1) ,

hη (1) , ψ − u(1) iH −1/2 (Γ) = 0

∀η (1) ∈ H −1/2 (Γ),

hσ (1) + σ (2) , ρiH −1/2 (Γ) = 0

∀ρ ∈ H 1/2 (Γ),

hη

(2)

,ψ − u

(2)

a(u

(2)

,v

(2)

iH −1/2 (Γ) = 0

) − hσ

(2)

,v

(2)

iH −1/2 (Γ) = hf, v

(2)

i

∀η

(2)

∀v

(2)

∈H

−1/2

∈X

(2)

(13)
(Γ),

.

It can be proven ([7, Proposition 1.7.1]) that if u is the solution of W1 and u(i) , σ (i) , ψ are solutions of
W4, then u(i) = u|Ωi , σ (i) = (∇L u · ni )Γ (where ∇L u · ni indicates the conormal derivative of u with
respect to the normal vector ni ), and ψ = u|Γ .
The weak formulation W2 we derived in Section 2 can be interpreted as a particular case of W4.
1/2
Indeed, let us firstly restrict the space H 1/2 (Γ) to its embedded subset H00 (Γ) and let us consider
−1/2
the particular case in which λ = σ (1) = −σ (2) ∈ Λ = H00 (Γ): then, the third equation in Eq. (13)

8

1/2

is automatically satisfied for all choices of ρ ∈ H00 (Γ). Moreover, subtracting the second and fourth
equations evaluated at the same η (1) = η (2) = η ∈ Λ yields
hη, ψ − u(1) iΛ − hη, ψ − u(2) iΛ = hη, u(2) − u(1) iΛ .

(14)

1/2

Obviously, this duality is well defined only if u(2) − u(1) belongs to H00 (Γ). We, therefore, set u ∈ X00
such that u(1) = u|Ω1 and u(2) = u|Ω2 ; Eq. (14) can be then rewritten as b(u, η) = 0 for all η ∈ Λ, i.e. the
second equation in Eq. (7). The first equation in Eq. (7) is found by adding the first and last equations in
(13) tested for all v ∈ X00 such that v (1) = v|Ω1 , v (2) = v|Ω2 ; observe that also in this case it is necessary
to restrict the search space for v to X00 , in order to ensure the well-posedness of b(v, λ) = hλ, v (2) − v (1) iΛ .
Although W2 and W4 are equivalent, their discretizations are not. Indeed, in the three-field method,
it is necessary to define the discretizations of the variational spaces of σ (1) , σ (2) and ψ. In contrast, when
discretizing W2, the third equation of W4 is not approximated but solved exactly and the second and
fourth are merged into a single equation. As we have shown, setting σ (1) = −σ (2) is efficient because
it allows to automatically satisfy the third equation in Eq. (13), thus reducing the number of variables.
Our approach limits to one the number of spaces to be discretized for each interface, thus allowing better
control of the stability of the method.
5. Inf-sup condition of the discretized problem
Problems W2 and W3 are saddle-point problems [21]. As such, their well-posedness depends on the
Ladyschenskaja-Babuška-Brezzi inf-sup condition [36], which sets the requirements for the uniqueness of
the solution as well as the stability of the sequence of problems depending on the discretization parameters
(e.g. the mesh size h or the number of basis functions on the interface nΓ ). We refer the reader to [36] and
[37] for a comprehensive description of the inf-sup condition from the functional and algebraic point of
view respectively. In this section, we specifically address the well-posedness of W3, and we limit ourselves
to recall that if the space Λ is characterized as in (5), the continuous problem W2 has a unique solution
[23].
Before stating the main stability result for W3, we recall that we characterize a(·, ·) and b(·, ·) as
continuous if there exist κa > 0 and κb > 0 such that a(ϕ, ψ) ≤ κa kϕkX kψkX for every ϕ, ψ ∈ X and
b(ϕ, ξ) ≤ κb kϕkX kξkΛ for every ϕ ∈ X , ξ ∈ Λ.
The following theorem prescribes the conditions for the well-posedness of W3.
Theorem 1. [37, Theorem 3.2] Assume that a(·, ·) and b(·, ·) are continuous with constants ka > 0 and
kb > 0, and that there exist α > 0 and β > 0 such that X h , V h,δ and Λδ satisfy the conditions
inf

v h ∈V h,δ

inf

a(v h , v h )
≥ α,
kv h k2X
sup

η δ ∈Λδ v h ∈X h

b(v h , η δ )
≥ β.
kv h kX kη δ kΛ

(15)

Then W3 has a unique solution. Moreover, there exists a constant C ≥ 0, depending only on κa , κb , α
and β, such that


h
δ
h
δ
ku − u kX + kλ − λ kΛ ≤ C
inf ku − v kX + inf kλ − η kΛ ,
(16)
v h ∈X h

η δ ∈Λδ

where (u, λ) is the solution of W2.
5.1. Numerical computation of the inf-sup constant
The inf-sup condition (15) is satisfied whenever Λδ is sufficiently “small” compared to X h . In the applications in Section 6 we ensure that β exists by numerically computing an approximation βe with the ap(1)
(2)
(1)
(2)
proach presented in [38], which we briefly summarize here. Let us suppose that XX ∈ R(n +n )×(n +n )
nΓ ×nΓ
h 2
δ 2
and XΛ ∈ R
are norm matrices such that kv kY = (XY v, v) and kη kΛ = (XΛ η, η) for every
v ∈ X h and every η ∈ Λδ . In the previous expressions, we denoted (·, ·) the standard scalar product in

9

Rm (m = n(1) + n(2) or m = nΓ ) and v and η the vectors of degrees of freedom of v h and η δ . Then, we
have
b(v h , η δ )
(Bv, η)
= inf sup
η6=0 v6=0 (XX v, v)1/2 (XΛ η, η)1/2
η δ ∈Λδ v h ∈X h kv h kX kη δ kΛ


−1/2 T
w,
X
B
η
X
1
= inf
sup
η6=0 (XΛ η, η)1/2
(w,
w)1/2
1/2
w=XX v6=0

1/2

1/2
−1/2
−1/2
XX B T η, XX B T η
BXX−1 B T η, η
= inf
= inf
.
η6=0
η6=0
(XΛ η, η)1/2
(XΛ η, η)1/2

βe = inf

sup

Introducing now the following generalized eigenvalue problem
"
# 

 
0 0
v
XX B T v
= −σ
,
η
0 XΛ η
B
0

(17)

and recognizing that we have

BXX−1 B T η = σXΛ η

⇒

σ=

BXX−1 B T η, η



(XΛ η, η)

,

(18)

we conclude that βe can be computed as the square root of the minimum eigenvalue of Eq. (17), i.e.
√
βe = σmin . For an application of this strategy, we refer the reader to the results presented in Fig. 5.
5.2. Convergence result for saddle-point problems
We close this Section by focusing on the convergence of problem W3. The following theorem gives a
sharper bound than Eq. (16) to the estimate of the approximation error.
Theorem 2. [21, Theorem 16.6] Let the assumptions of Theorem 1 be satisfied. Then the solution (u, λ)
of W2 and the solution (uh , λδ ) of W3 satisfy the following error estimates


κb
κa
inf ku − v∗h kX +
inf kλ − η δ kΛ ,
(19)
ku − uh kX ≤ 1 +
α v∗h ∈V h,δ
α ηδ ∈Λδ




κa
κa
κb
κa κb
kλ − λδ kΛ ≤
1+
inf ku − v∗h kX + 1 +
+
inf kλ − η δ kΛ .
β
α v∗h ∈V h,δ
β
αβ ηδ ∈Λδ
Moreover, the following error estimate holds
inf

v∗h ∈V h,δ

ku − v∗h kX ≤


1+

κb
β



inf ku − v h kX .

v h ∈X h

(20)

Theorem 2 shows that, whenever the space of Lagrange multipliers is rich enough (namely the second
term in Eq. (19) becomes negligible compared the first one), the approximation of u is essentially bounded
by the best approximation error on X . However, this richness may lower the inf-sup constant β and
therefore loose the approximation (20). It is therefore important to find the correct balance.
We remark that in Eq. (9) we have that increasing the size of the Lagrange multipliers space is
equivalent to lowering the size of V h,δ and, consequently, the supremum in its right hand side. The
two error estimates in (9) and (19) are therefore two equivalent ways of expressing the fact that, if the
continuity over the interface Γ is enforced strongly enough, the error converges to zero as the error due
to the spatial discretization in X h . As we show in the next section, we are then able to recover the usual
convergence orders for u with respect to the mesh size h when using the finite element method.

10

Non-orthonormal basis functions

Orthonormal basis functions

condition number

1019

h = 1/20
h = 1/28
h = 1/40
h = 1/56
h = 1/80

10

10
1014

108

109

106
104

104
1

7

13

19

25

31

1

nΓ

7

13

19

25

31

nΓ

Figure 2: Condition number of the discretized matrix in Eq. (11) for the Poisson problem on two subdomains (see Section 6.2)
with conforming meshes vs number of basis functions for the Lagrange multiplier space, with different refinement levels of
the (conforming) meshes. On the left, we consider non-orthonormal “half” Fourier basis functions (21), while on the right
we consider their orthonormalization.

6. Numerical results
In this section, we focus on the performance of the method presented in Section 2 on two-dimensional
problems defined over the unit square. The numerical simulations we present are performed with a set of
Matlab scripts which can be freely downloaded1 .
For all the simulations, we employ standard piecewise polynomial Lagrangian basis functions defined
over suitable triangulations in the subdomains. Regarding the choice of basis functions for Λδ , we already
anticipated in Section 2 that in this paper we investigate the possibility of using low-frequency Fourier
basis functions built on the interface Γ.
6.1. Choice of basis functions for the Lagrange multipliers
Given an interface with length L, we consider ξ1 = 1 and, for i = 1, . . . , nω
ξ2i (s) = sin(ωi πs),

ξ2i+1 (s) = cos(ωi πs),

(21)

where s is the arc length of the interface Γ, ωi = i/L, and nω is the number of considered frequencies; it
Γ
holds that nΓ = 2nω +1. With this definition, the set {ξi }ni=1
forms an orthogonal basis with respect to the
2
L (0, 2L) scalar product. We choose to employ such basis – instead of the standard Fourier basis orthogonal (or orthonormal) with respect to the L2 (0, L) scalar product – because, by considering basis functions
with periodicity L, we would impose an unnecessary periodicity constraint, in particular, the equality of
the functions in Λδ and their derivatives at the extrema of Γ. As a result, we empirically observed that
by employing the standard L2 (0, L) orthonormal Fourier basis functions the optimal convergence of the
finite element method is retrieved for larger values of nΓ compared to the choice in Eq. (21). However,
utilizing non-orthonormal basis functions (21) has a dramatic influence on the condition number of the
resulting linear system, which has exponential growth with the increasing number of basis functions on
the interface; see Fig. 2 (left).
In order to retain the convergence order attained by using the Fourier modes in Eq. (21) and, at
the same time, control the condition number of the system, we propose an orthonormalization strategy
based on the Gram-Schmidt algorithm or, equivalently, on the QR decomposition [39]. Even though
the coefficients of the orthonormal basis generated by (21) with these algorithms could be analytically
derived, their exact expression quickly becomes complex with nΓ becoming large. With our approach, we
aim at obtaining an approximation of such coefficients relying on a fine sampling of the basis functions
on the interval (0, L). We remark that, in addition to allowing to effortlessly compute a large number of
orthonormal basis functions, our approach has the advantage to be general enough to be applied to any
set of non-orthonormal basis functions.
Γ
s
Let {ξi }ni=1
be the set of non-orthonormal basis functions defined on Γ. Moreover, let {xi }ni=1
be
distinct sample points distributed over the interval (0, L), where L still denotes the length of the interface.
1 https://github.com/lucapegolotti/coupling_scripts

11

1

1.5

y

1
0.5

0.5
0

0

0

0.5
x

10

0.5
x

10

0.5
x

1

1

−0.5

∂uh,(1) /∂x on Γ
∂uh,(2) /∂x on Γ
λδ
∂uex /∂x on Γ

0
−1
−2
−3

0

0.2

0.4

0.6
y

0.8

10

0.2

0.4

0.6

0.8

10

0.2

0.4

y

0.6

0.8

1

y

Figure 3: Contour lines of the solution (top row) and derivative of the solution at the interface Γ (bottom row) when nΓ = 1
(left column), nΓ = 3 Fourier modes (middle column) and nΓ = 5 Fourier modes (right column) are used to characterize the
space Λδ . The red and blue solid lines in the bottom plots represent the partial derivatives with respect to x – computed at
the interface Γ – of the numerical solutions in Ω1 and Ω2 respectively.
PnΓThe dash-dot green line is computed by reconstructing
λi ξi . The results are obtained on a mesh conforming
λδ as a linear combination of the Fourier basis functions, i.e. λδ = i=1
at the interface, with quadratic Lagrangian polynomials on both domains and with mesh size h = 1/20.

s
We now introduce the functions {κi }ni=1
, which we identify with the set of standard LagrangianR piecewise
linear basis functions centered on each sample point xi , and the associated mass matrix Mij = Γ κi κj dx.
Let V = [v1 , v2 , . . . , vnΓ ] ∈ Rns ×nΓ be the matrix of the evaluations of the basis functions on the sample
points, namely Vij = ξj (xi ). We remark that, for each i = 1, 2, . . . , nΓ , we have kξi k2L2 (0,L) ≈ vTi M vi . Since
M is a positive-definite matrix, it admits a unique Cholesky decomposition and there exists C ∈ Rns ×ns
such that C T C = M . Let us now consider the unit matrix Q ∈ Rns ×nΓ and the upper triangular matrix
R ∈ RnΓ ×nΓ such that the truncated QR decomposition of CV reads

CV = QR.
By construction, we have
(C −1 Q)T M C −1 Q = QT C −T C T CC −1 Q = I,
thus, the columns of C −1 Q represent evaluations at the sample points of functions orthonormal on (0, L)
with respect to the L2 product. The matrix R performs the change of variable from the frame of reference
of the new orthonormal basis functions to the frame of reference of the non-orthonormal basis functions.
If the sampling is sufficiently fine, we speculate that the elements of the matrix R−1 well approximate
the coefficients which are computed by applying the Gram-Schmidt algorithm to the continuous nonΓ
orthonormal basis functions {ξi }ni=1
and, in particular, that
ξiGS =

nΓ
X

−1
ξj Rji
.

(22)

j=1

From a practical perspective, the matrix R−1 is suitable to compute the coupling matrix B GS with
respect to the orthonormal Fourier basis functions, knowing the coupling matrix computed without
12

orthonormalization B. Indeed, we have
GS
Bij
=

Z
Γ

ξiGS ϕj ds =

nΓ Z
X
k=1

Γ


nΓ
X
−1
−1
ξk ϕj ds Rki
=
Bkj Rki
,
k=1

GS
or equivalently Bij
= R−T B. Therefore, the condition number of the system can be controlled by
multiplying the coupling matrices B by the matrix R−1 ; observe that, being R an upper triangular
matrix, the application of its inverse is performed with negligible cost. The matrix R−1 depends only
on the choice of the non-orthonormal basis functions and can be then computed a priori. We remark
that, with this approach, the orthonormal basis functions are never explicitly computed. Moreover, since
the discrete space is exactly the same, the approximation properties and the convergence orders are not
changed. Fig. 2 (right) shows that, after the orthonormalization of the Fourier basis functions (21) by the
algorithm we presented, the system is more stable and the condition number increases with the number
of Fourier basis functions nΓ dependently on the refinement level of the mesh h.

6.2. The Poisson problem
Let us consider the global Poisson problem (2) on the domain Ω = (0, 1)2 , where we take f such that
uex = 100xy(1 − x)(1 − y) sin(1/3 − xy 2 ) is the exact solution. We divide Ω into Ω1 = (0, 0.5) × (0, 1) and
Ω2 = (0.5, 1) × (0, 1).
We numerically solve the problem on Ω1 and Ω2 by employing structured triangular conforming and
non-conforming meshes with varying mesh size h. The conforming meshes are obtained by subdividing the
domain in the x- and y-direction in the same number of elements. On the other hand, the non-conforming
meshes are built by taking in the y-direction of Ω2 N + 1 elements, N being the number of elements in
the y-direction in Ω1 as well as the total number of elements in the x-direction.
Fig. 3 shows how the solutions on Ω1 and Ω2 obtained with a conforming mesh with N = 20 elements
in each direction change with respect to the number of basis functions on the interface. The results are
obtained with quadratic Lagrangian polynomials in both subdomains. From the contour lines plots in the
top row, it appears that the two solutions match quite accurately at the interface with 5 Fourier basis
functions (nω = 2). In the second row of Fig. 3, we plot the approximation by finite differences of the
derivative of the solution with respect to x in the two domains, which is equal to the normal derivative
of uh,(1) and to the opposite of the normal derivative on uh,(2) on Γ respectively. Observe that, as we
already highlighted in Remark 2, the Lagrange multiplier λδ takes the role of the normal derivative of uh
on Γ.
Let us address the convergence of the global solution to the exact one with respect both to the mesh
size h and the number of basis functions on the interface nΓ . To this end, we consider meshes with
total number of elements in the x-direction N = 20, 28, 40, 56, 80, 114, 160 and we solve the problems
by employing quadratic Lagrangian basis functions in both subdomains. Fig. 4 (top row) depicts the
decaying of the error in X -norm (the broken norm) with respect to h, as well as the convergence of the
error obtained by solving the problem on a single mesh of Ω (in black dashed line). When employing
both conforming and non-conforming meshes, the error is optimal – in the sense that we recover the
theoretical order of convergence h2 of quadratic finite elements for the H 1 -error – when nΓ is large
enough, e.g. nΓ ≥ 13. If nΓ is too small, on the contrary, the solution is unable to converge to the exact
solution with h and reaches a stagnation point. We remark that this result is perfectly consistent with
Strang’s second lemma (9) and with the stability result in Theorem 2: whenever the space of Lagrange
multiplier is rich enough (which is equivalent to requiring that V h,δ be a good approximation of V), the
best approximation error of the interpolation is recovered.
Remark 4. In our numerical simulations with non-conforming meshes, we observed that instabilities
arise when using coarse meshes and low-order quadrature rules for the computations of the approximate
integrals of B1 and B2 in Eq. (12). Fig. 4 (bottom row, left) shows that, when using for example 2 Gauss
quadrature nodes, the error increases with nΓ when h = 1/20, 1/28, 1/40. By increasing the order of the
quadrature rule and choosing 4 Gauss quadrature nodes this issue is completely fixed; see Fig. 4 (bottom
row, right). The plots in Fig. 4 (right column) are obtained from the same data. We did not encounter
stability problems when using conforming meshes, even with low-order quadrature rules.
Fig. 5 shows the variation of the estimate of the inf-sup constant βe – computed as the square root
of the minimum eigenvalue of the generalized eigenvalue problem (17), as described in Section 5.2 –
when the number of basis functions on the interface changes; the estimate refers to the simulation of
13

Non-conforming mesh

Conforming mesh
10

0

nΓ = 1
nΓ = 3
nΓ = 5
nΓ = 7
nΓ = 9
nΓ = 13
Global solution

error

10−1
10−2

10−4

h2

h2

10−3

6.3 · 10−3

1.8 · 10−2
h

6.3 · 10−3

5 · 10−2

2 Gauss quadrature nodes
10

1.8 · 10−2
h

5 · 10−2

4 Gauss quadrature nodes

0

h
h
h
h
h
h
h

error

10−1
10−2

=
=
=
=
=
=
=

1/20
1/28
1/40
1/56
1/80
1/114
1/160

10−3
10−4

1

7

13

19

25

31

1

7

13

nΓ

19

25

31

nΓ

Figure 4: The top row shows the convergence of global error with respect to the mesh size h and number of basis functions
on the interface nΓ with conforming and non-conforming meshes. The colored lines show the errors in the X -norm, i.e. the
broken norm. The black dashed line represents the H 1 -error obtained by solving the global problem on a conforming mesh,
which can be regarded as the union of two conforming meshes on Ω1 and Ω2 . The bottom row shows the convergence of
the error on the non-conforming meshes with respect to nΓ with varying mesh sizes; in particular, the discretization of the
terms corresponding to the bilinear form b(·, ·) in Eq. (12) has been performed with 2 and 4 Gauss quadrature nodes on
the left and on the right plots respectively. Note that the two plots at the right both refer to the same data obtained with
a 4-node quadrature rule.

h = 1/20
h = 1/28
h = 1/40
h = 1/56
h = 1/80
h = 1/114
h = 1/160

1.5

βe

1

0.5

0

1

7

13

19

25

31

NΓ
Figure 5: Decaying of the inf-sup constant βe with respect to nΓ computed on the Poisson problem on conforming meshes,
using quadratic polynomial basis functions on both subdomains and the orthonormal basis functions ξiGS on the interface.
The constant is approximated as the square root of the minimum eigenvalue of Eq. (18).

14

Ω1

Ω2

100

H 1 -error

10−1

nΓ = 1
nΓ = 3
nΓ = 5
nΓ = 7
nΓ = 9
nΓ = 13
Global solution

h

10−2

h2

10−3
10−4

6.3 · 10−3

1.8 · 10−2
h

5 · 10−2

6.3 · 10−3

1.8 · 10−2
h

5 · 10−2

Figure 6: Convergence of the H 1 -error in Ω1 (left) and Ω2 (right) when linear and quadratic Lagrangian basis functions are
employed in the two domains respectively. The meshes are not conforming at the interface. The black dashed lines refer to
the global H 1 -error obtained with linear (on the left) and quadratic (on the right) elements on the whole Ω.

the Poisson equations with conforming meshes and quadratic polynomial basis functions. Due to the
−1/2
difficulties in computing the H00 - norm for the Lagrange multiplier, we replaced the estimate given by
Eq. (17) with a surrogate where the space X h is substituted with the space spanned by the traces on Γ
of the finite element basis functions ϕi ; the L2 -norm is used both for such space and Λδ . For the result
in Fig. 5, we employed the orthonormal basis functions ξiGS computed as in Eq. (22), so that XΛ = I;
therefore, from Eq. (18) it follows that βe is simply found as the square root of the minimum eigenvalue
of B GS XX−1 (B GS )T . In Fig. 5, each curve presents a plateau phase in which the inf-sup constant stays
approximately constant at around βe ≈ 1.41 with the increment of nΓ . The amplitude of such plateau
phase increases when h becomes smaller. Indeed, we observe that βe starts decreasing for smaller values
of nΓ when the meshes are coarser and that, conversely, for finer meshes the inf-sup constant varies
relatively little in the range nΓ ∈ (1, 31). We remark that, combined with the condition number shown
in Fig. 2 (right), this result ensures that for each refinement level, we are able to obtain the optimal
convergence of the finite element method when the basis functions are orthonormal. Indeed, refining the
mesh has the effect of both increasing the range of stability of the linear system – see Fig. 2 (right) – and
increasing the number of basis functions at the interface that can be employed without reaching the fast
decaying region of βe in Fig. 5. With regard to this last point, we recall that it is important to prevent the
inf-sup constant to become too small because it appears at the denominator of the constant multiplying
the best approximation errors on u and on λ in the error estimates of Theorem 2.
We focus now on the solution of the problem when employing non-conforming meshes, linear Lagrangian basis functions in Ω1 and quadratic Lagrangian basis functions in Ω2 . Fig. 6 shows that the
H 1 -error scales in the two subdomains as the best approximation error of the local (to the subdomain)
basis: we recover first order convergence in Ω1 and second order convergence in Ω2 . We remark that the
convergence of the global error in the broken norm is determined by the rate in Ω1 – being the error in
such subdomain much larger than that in Ω2 – and it is of first order. In Fig. 6 we also show with black
dashed lines the global H 1 error obtained when solving the problem with linear (in the left plot) and
quadratic (in the right plot) basis functions on the whole Ω. As expected, the accuracy obtained with
mixed polynomial degrees lies between the accuracies achieved while using only linear and only quadratic
basis functions.
6.3. The Navier Stokes equations
In this section, we test the flexibility of our method by solving the Navier-Stokes equations on Ω =
(0, 1) × (0, 1)
−µ∆u + (u · ∇)u + ∇p = f in Ω,
divu = 0

in Ω,

u = g on ΓD ,

(23)

σ(u, p)n = h on ΓN ,
where u and p are velocity and pressure respectively, µ ∈ R is the viscosity, ΓD and ΓN are portions of
the boundary such that ΓD ∩ ΓN = ∂Ω and ΓD ∪ ΓN = ∅, f is a given forcing term, g and h are the
15

1

1

⌦1

0.7

⌦5
4

3

⌦4

0.15

0.35

1

⌦3

2

⌦2

ku − uh kH 1 (Ω) + kp − ph kL2 (Ω)

10−2

10−3

h2

10−4

10−5
7.8 · 10−3

0.15

2.3 · 10−2

6.25 · 10−2

h
Figure 7: On the left, domain decomposition and computational mesh of Ω = (0, 1)×(0, 1), namely the domain considered in
the numerical simulation of the Navier-Stokes equations. Each triangulation T h is characterized by the following mesh sizes:
h in Ω1 , h/2 in Ω2 , h in Ω3 , h/2 in Ω4 and h in Ω5 . On the right, convergence of the error ku − uh kH 1 (Ω) + kp − ph kL2 (Ω)
against the exact solution (24) with respect to h. The number of basis functions at each interface is the same for each
(1)
(2)
(3)
(4)
evaluation of the error: we set nΓ = 22, nΓ = 18, nΓ = 18 and nΓ = 14.

Dirichlet and Neumann boundary data respectively, and
σ(u, p) = µ∇u − pI
is the stress tensor. The domain is partitioned into five subdomains Ωi with i = 1, . . . , 5 and divided by
four interfaces Γi with i = 1, . . . , 4, as shown in Fig. (7) (left). We define a family of non-conforming triangulations T h characterized by the mesh size h, i.e. the maximum edge length over Ω , which corresponds
to the mesh size in Ω1 ; Ω3 and Ω5 are characterized by approximately the same mesh size h, whereas
Ω2 and Ω4 have mesh size h/2. The choice of using meshes refined in the bottom left and bottom right
subdomains is motivated by the numerical simulation presented in Section 6.3.2, in which the solution of
the Navier-Stokes equations presents recirculation zones localized in these regions of Ω. Differently from
what done in Section 6.2, we only focus on non-conforming meshes and we exclusively employ the inf-sup
stable Taylor-Hood [40] elements with quadratic Lagrangian basis functions for the velocity and linear
Lagrangian basis functions for the pressure. The non-linear discretized system is numerically solved by
Newton’s method.
Remark 5. When applied to the Navier-Stokes equations in two dimensions, the method requires assigning
to each interface two sets of basis functions discretizing the two components of the normal stress. To see
why this is the case, consider the situation in which Ω is subdivided into Ω1 and Ω2 ; let us denote as
always the interface of the two partitions Γ. Multiplying the momentum equation by a test function
v ∈ [HΓ1D (Ω)]2 and integrating by parts on Ω1 leads to
Z
Z
Z
Z
Z
µ
∇u : ∇v dx −
p∇ · v dx − σ(u, p)n · v dx =
f · v dx +
h · v dx.
Ω1

Ω1

Γ

Ω1

∂Ω1 ∩ΓN

The integral on Γ is the coupling term. Each of the two components of the normal stress σ(u, p)n must
be discretized by a set of basis functions. In this paper, we choose for simplicity to use the same set for
the two components of the normal stress.
We consider again Fourier basis functions for the approximation of the normal stresses. Since, as
explained in Remark 5, we need two Lagrange multipliers for representing each normal stress, the number
(i)
of basis functions on Γi is found as nΓ = 2(2nω(i) + 1), where nω(i) is the number of frequencies used on
th
the i interface.

16

6.3.1. Numerical convergence against the exact solution
We consider Eq. (23) with µ = 1, ΓN = {1} × (0, 1), ΓD = ∂Ω \ ΓN , and f , g, and h chosen such that


1
sin(yπ)
u=
(24)
,
p = − x2
exp(x)
2
is the exact solution.
(1)
(2)
(3)
(4)
On the interfaces we set nΓ = 22, nΓ = 18, nΓ = 18 and nΓ = 14. The number of basis functions
on the interfaces is chosen such that the optimal convergence of finite elements is retrieved. Fig. 7 (right)
shows that the following classical error estimate for Taylor-Hood elements
ku − uh kH 1 (Ω) + kp − ph kL2 (Ω) ≤ Ch2

(25)

holds. The norms in (25) should be interpreted as broken norms.
6.3.2. Lid-driven cavity problem
We now focus on the numerical approximation of the classic lid-driven cavity problem [41, 42] with
Reynolds number Re = 500. Specifically, we consider Eq. (23) with µ = 1, ΓD = ∂Ω, g = [U, 0]T with
(1)
U = 500 on (0, 1) × {1} and U = 0 on the rest of the boundary, and f = [0, 0]T . We consider nΓ = 42,
(2)
(3)
(4)
nΓ = 18, nΓ = 18 and nΓ = 14; as in the numerical simulation presented in Section 6.3.1, the number
of basis functions on the interfaces is chosen such that the error on the Lagrange multipliers can be
considered negligible if compared with the finite element error. We remark that, being the solution of the
problem at hand considerably more difficult to capture accurately than the exact solution (24) – because
it features steep gradients and higher Reynolds numbers – it became necessary to increase the number of
basis functions on Γ1 in order to obtain optimal convergence. Since we consider only Dirichlet boundary
conditions, the problem is not well-posed as the pressure is unique up to a constant. We deal with this
issue by fixing the degree of freedom of the pressure in the bottom left corner to zero.
Fig. 8 shows the streamlines obtained by solving the problem with non-conforming meshes belonging
to the family T h characterized by h = 1/16, h = 1/32 and h = 1/128. Furthermore, the streamlines
of a fine solution computed with uniform h = 1/300 are displayed for reference: these are qualitatively
similar to the ones corresponding to Re = 500 reported in e.g. [43], and we, therefore, assume that the
fine solution well approximates the exact solution of the problem. As shown in Fig. 8 (bottom right),
using a smaller mesh size in the regions of the two smaller eddies in the lower part of the domain allowed
us to obtain satisfactory approximations of those secondary recirculation zones even with the coarsest
mesh size h = 1/16. As expected, the differences in the streamlines among the different refinement levels
are more evident in Ω1 , were for each mesh belonging to T h the largest elements are located. Specifically,
we notice that, while the position of the primary eddy is approximated with good accuracy even for
h = 1/16, in the peripheral regions of the domain only the streamlines corresponding to h = 1/128 are
almost indistinguishable from the ones of the exact solution; see Fig. 8 (top right) and Fig. 8 (bottom
left). Table 1 provides the number of finite element degrees of freedom corresponding to each refinement
level and the error in the approximation of the center of each eddy: this quantity is computed as the
Euclidean distance of the points where the velocity field attains minimum velocity (in magnitude) in the
coarse solutions and in the fine solution. The errors become smaller with h. We remark that, for coarse
meshes, the approximation error of the two smaller eddies in the lower part of Ω is one order of magnitude
lower than that of the central eddy; this confirms that employing the non-conforming meshes in T h , which
are characterized by a smaller element size in Ω2 and Ω4 , leads to a satisfactory approximation of the
secondary recirculation regions even for large values of h. It is worth noting that the number of degrees
of freedom reserved to the discretization of the Lagrange multipliers is constant for all the meshes (as
it depends solely on the number of basis functions at each interface) and it is equal to 92: this quantity
is much smaller than the number of degrees of freedom for velocity and pressure for each h. Hence,
the coupling of the finite element spaces is performed by introducing a negligible number of additional
variables.
7. Conclusions
We presented a non-conforming domain decomposition method for non-overlapping subdomains. At
the continuous level, our method and the mortar method are based on the same weak formulations in
which the continuity constraints over the primal (the solution) and the dual (the stresses) variables are
17

1

1

h = 1/16
h = 1/32
h = 1/128
Fine solution

0.8
0.95

y

y

0.6
0.4
0.9
0.2
0

0

0.2

0.4

0.6

0.8

0.85

1

0

5 · 10−2

0.1

x

0.15

x

1

0.4

0.3

y

y

0.95
0.2

0.9
0.1

0.85
0.85

0.9

0.95

1

0
0.6

0.7

0.8

x

0.9

1

x

Figure 8: Streamlines of the lid-driven cavity problem computed with meshes belonging to the family T h , with h = 1/16,
h = 1/32 and h = 1/128. The streamlines of a fine reference solution obtained on a uniform mesh with h = 1/300 are
displayed in black dashed lines. The plots on the top right and bottom row show details over regions of Ω.

dofs

error

h

velocity

pressure

system size

E1

E2

E3

1/16
1/32
1/64
1/128

3’214
12’162
46’982
184’178

445
1’602
6’031
23’334

3’751
13’856
53’105
207’604

2.9e–3
2.4e–4
1.5e–5
8.3e–7

6.5e–4
5.0e–5
4.5e–6
1.6e–6

8.3e–4
8.9e–5
8.4e–6
1.3e–6

Table 1: Degrees of freedom (dofs) for the meshes considered in the lid-driven cavity problem and approximation error of the
center of the three eddies (E1: primary central eddy, E2: bottom right eddy, E3: bottom left eddy). The system size is the
sum of the degrees of freedom of velocity and pressure and the number of degrees of freedom for the Lagrange multipliers,
the latter being constant and equal to 92 for all the meshes. The center of the eddy is numerically found as the point where
the minimum of the magnitude of the velocity field is reached, and the error is computed as the Euclidean distance of such
approximation with the center of the eddy of the fine solution obtained over a uniform mesh with h = 1/300, corresponding
to 722’402 degrees of freedom for the velocity and 90’601 degrees of freedom for the pressure. The coordinates of the centers
of E1, E2 and E3 for the fine solution are: [0.545907,0.593810], [0.879935,0.121555] and [0.059761,0.053354] respectively.

18

enforced via Lagrange multipliers. As we described in the paper, our choice of discretizing the space of
Lagrange multipliers independently of the spatial discretization in the subdomains offers the advantage of
a straight-forward implementation of the method and the possibility of tuning the accuracy of the coupling
as required by the application; we limited ourselves to considering Fourier basis functions defined over
the interface. However, the saddle-point nature of the problem poses constraints over the richness of the
discretized space for the Lagrange multipliers compared to the degrees of freedom of the primal variable:
we empirically verified that the inf-sup constant can be controlled dependently on the mesh size, in the
sense that finer meshes allow considering larger number of Fourier basis functions, without violating the
inf-sup stability. In the numerical experiments, we showed that the optimal convergence of the finite
element method was recovered for the Poisson problem; this was confirmed both when using conforming
and non-conforming meshes, and when using different polynomial degrees in the subdomains. In the
last part of the paper, we showed that the method can be easily extended to the case of non-elliptic
equations, such as the Navier-Stokes equations, and to cases of partitions of the domain into multiple
subdomains. We were able to recover the optimal convergence rate of finite elements also for the NavierStokes equations by considering a number of basis functions on the interfaces considerably lower than
the number of degrees of freedom of the discretized subdomains. Moreover, we focused on a possible
practical application of the method, i.e. the use of non-conforming structured meshes in the lid-driven
cavity problem for capturing the secondary recirculation regions. We showed that considering smaller
mesh sizes in correspondence of the secondary eddies ensures satisfactory results in terms of accuracy of
the streamlines of the vortex rings.
Acknowledgments
The authors are grateful to Prof. Annalisa Buffa and Prof. Alfio Quarteroni for the fruitful discussions
and their advice on the topics presented in this paper. The research of the authors is supported by the
Swiss National Foundation (SNF), project No. 140184.
References
[1] S. Deparis, M. Discacciati, A. Quarteroni, A domain decomposition framework for fluid-structure
interaction problems, Computational Fluid Dynamics 2004 (2006) 41–58.
[2] A. Toselli, O. B. Widlund, Domain decomposition methods: algorithms and theory, volume 34 of
Springer Series in Computational Mathematics, Springer, 2005.
[3] M. Israeli, L. Vozovoi, A. Averbuch, Domain decomposition methods for solving parabolic pdes on
multiprocessors, Applied numerical mathematics 12 (1993) 193–212.
[4] R. Becker, P. Hansbo, R. Stenberg, A finite element method for domain decomposition with nonmatching grids, ESAIM: Mathematical Modelling and Numerical Analysis 37 (2003) 209–225.
[5] C. Bernardi, A new nonconforming approach to domain decomposition: the mortar element method,
Nonlinear partial equations and their applications (1989).
[6] C. Bernardi, Y. Maday, F. Rapetti, Basics and some applications of the mortar element method,
GAMM-Mitteilungen 28 (2005) 97–123.
[7] A. Quarteroni, A. Valli, Domain decomposition methods for partial differential equations, Numerical
Mathematics and Scientific Computation, Oxford University Press, 1999.
[8] M. A. Puso, T. A. Laursen, A mortar segment-to-segment contact method for large deformation
solid mechanics, Computer methods in applied mechanics and engineering 193 (2004) 601–629.
[9] M. A. Puso, A 3d mortar method for solid mechanics, International Journal for Numerical Methods
in Engineering 59 (2004) 315–336.
[10] A. Ehrl, A. Popp, V. Gravemeier, W. Wall, A dual mortar approach for mesh tying within a
variational multiscale method for incompressible flow, International Journal for Numerical Methods
in Fluids 76 (2014) 1–27.

19

[11] T. Klöppel, A. Popp, U. Küttler, W. A. Wall, Fluid–structure interaction for non-conforming interfaces based on a dual mortar formulation, Computer Methods in Applied Mechanics and Engineering
200 (2011) 3111–3126.
[12] A. Popp, W. Wall, Dual mortar methods for computational contact mechanics–overview and recent
developments, GAMM-Mitteilungen 37 (2014) 66–84.
[13] C. Hesch, A. Gil, A. A. Carreño, J. Bonet, P. Betsch, A mortar approach for fluid–structure interaction problems: Immersed strategies for deformable and rigid bodies, Computer Methods in Applied
Mechanics and Engineering 278 (2014) 853–882.
[14] F. B. Belgacem, L. K. Chilton, P. Seshaiyer, The hp-mortar finite-element method for the mixed
elasticity and stokes problems, Computers & Mathematics with Applications 46 (2003) 35–55.
[15] S. Deparis, D. Forti, P. Gervasio, A. Quarteroni, INTERNODES: an accurate interpolation-based
method for coupling the Galerkin solutions of PDEs on subdomains featuring non-conforming interfaces, Computers & Fluids 141 (2016) 22–41.
[16] D. Forti, Parallel algorithms for the solution of large-scale fluid-structure interaction problems in
hemodynamics (2016).
[17] F. Brezzi, L. D. Marini, A three-field domain decomposition method, Contemporary Mathematics
157 (1994) 27–34.
[18] J. A. Cottrell, T. J. Hughes, Y. Bazilevs, Isogeometric analysis: toward integration of CAD and
FEA, John Wiley & Sons, 2009.
[19] T. J. Hughes, J. A. Cottrell, Y. Bazilevs, Isogeometric analysis: CAD, finite elements, NURBS, exact
geometry and mesh refinement, Computer methods in applied mechanics and engineering 194 (2005)
4135–4195.
[20] A. Quarteroni, A. Valli, Numerical approximation of partial differential equations, volume 23,
Springer Science & Business Media, 2008.
[21] A. Quarteroni, Numerical models for differential problems, volume 8, Springer-Verlag, 2014.
[22] S. Salsa, Partial differential equations in action: from modelling to theory, volume 99, Springer, 2016.
[23] D. Braess, W. Dahmen, C. Wieners, A multigrid algorithm for the mortar finite element method,
SIAM Journal on Numerical Analysis 37 (1999) 48–69.
[24] P.-A. Raviart, J. Thomas, Primal hybrid finite element methods for 2nd order elliptic equations,
Mathematics of computation 31 (1977) 391–413.
[25] D. Boffi, F. Brezzi, M. Fortin, et al., Mixed finite element methods and applications, volume 44,
Springer, 2013.
[26] S. K. Acharya, A. Patel, Primal hybrid method for parabolic problems, Applied Numerical Mathematics 108 (2016) 102–115.
[27] F. B. Belgacem, The mortar finite element method with Lagrange multipliers, Numerische Mathematik 84 (1999) 173–197.
[28] B. I. Wohlmuth, A mortar finite element method using dual spaces for the Lagrange multiplier,
SIAM journal on numerical analysis 38 (2000) 989–1012.
[29] P. G. Ciarlet, The finite element method for elliptic problems, SIAM, 2002.
[30] A. Quarteroni, R. Sacco, F. Saleri, Numerical mathematics, volume 37, Springer Science & Business
Media, 2010.
[31] P. Seshaiyer, M. Suri, Convergence results for non-conforming hp methods: The mortar finite element
method, Contemporary Mathematics 218 (1998) 453–459.

20

[32] P. Seshaiyer, Stability and convergence of nonconforming hp finite-element methods, Computers &
Mathematics with Applications 46 (2003) 165–182.
[33] P. Gervasio, A. Quarteroni, Analysis of the internodes method for non-conforming discretizations of
elliptic equations, MATHICSE report (2016).
[34] S. Deparis, D. Forti, A. Quarteroni, A rescaled localized radial basis function interpolation on noncartesian and nonconforming grids, SIAM Journal on Scientific Computing 36 (2014) A2745–A2762.
[35] F. Brezzi, D. Marini, Error estimates for the three-field formulation with bubble stabilization,
Mathematics of computation 70 (2001) 911–934.
[36] F. Brezzi, On the existence, uniqueness and approximation of saddle-point problems arising from
lagrangian multipliers, Revue française d’automatique, informatique, recherche opérationnelle. Analyse numérique 8 (1974) 129–151.
[37] F. Brezzi, K.-J. Bathe, A discourse on the stability conditions for mixed finite element formulations,
Computer methods in applied mechanics and engineering 82 (1990) 27–57.
[38] F. Ballarin, A. Manzoni, A. Quarteroni, G. Rozza, Supremizer stabilization of POD–Galerkin approximation of parametrized steady incompressible Navier–Stokes equations, International Journal
for Numerical Methods in Engineering 102 (2015) 1136–1161.
[39] A. Ruhe, Numerical aspects of Gram-Schmidt orthogonalization of vectors, Linear algebra and its
applications 52 (1983) 591–601.
[40] P. Hood, C. Taylor, Navier-stokes equations using mixed interpolation, Finite element methods in
flow problems (1974) 121–132.
[41] J. D. Bozeman, C. Dalton, Numerical study of viscous flow in a cavity, Journal of Computational
Physics 12 (1973) 348–363.
[42] U. Ghia, K. N. Ghia, C. T. Shin, High-Re solutions for incompressible flow using the Navier-Stokes
equations and a multigrid method, Journal of computational physics 48 (1982) 387–411.
[43] X. Shi, J. Khodadadi, Laminar fluid flow and heat transfer in a lid-driven cavity due to a thin fin,
Journal of Heat Transfer 124 (2002) 1056–1063.

21

