# BDDC preconditioners for divergence free virtual element discretizations of
the Stokes equations

Tommaso Bevilacqua, Simone Scacchi  Tommaso Bevilacqua 111Dipartimento di
matematica, Universitá degli studi di Milano, Via Saldini 50, 20133 Milano,
Italy; e-mail<EMAIL_ADDRESS>Simone Scacchi 222Dipartimento di
matematica, Universitá degli studi di Milano, Via Saldini 50, 20133 Milano,
Italy; e-mail<EMAIL_ADDRESS>
# BDDC preconditioners for divergence free virtual element discretizations of
the Stokes equations

Tommaso Bevilacqua, Simone Scacchi  Tommaso Bevilacqua 111Dipartimento di
matematica, Universitá degli studi di Milano, Via Saldini 50, 20133 Milano,
Italy; e-mail<EMAIL_ADDRESS>Simone Scacchi 222Dipartimento di
matematica, Universitá degli studi di Milano, Via Saldini 50, 20133 Milano,
Italy; e-mail<EMAIL_ADDRESS>
###### Abstract

The Virtual Element Method (VEM) is a new family of numerical methods for the
approximation of partial differential equations, where the geometry of the
polytopal mesh elements can be very general. The aim of this article is to
extend the balancing domain decomposition by constraints (BDDC) preconditioner
to the solution of the saddle-point linear system arising from a VEM
discretization of the two-dimensional Stokes equations. Under suitable
hypotesis on the choice of the primal unknowns, the preconditioned linear
system results symmetric and positive definite, thus the preconditioned
conjugate gradient method can be used for its solution. We provide a
theoretical convergence analysis estimating the condition number of the
preconditioned linear system. Several numerical experiments validate the
theoretical estimates, showing the scalability and quasi-optimality of the
method proposed. Moreover, the solver exhibits a robust behavior with respect
to the shape of the polygonal mesh elements. We also show that a faster
convergence could be achieved with an easy to implement coarse space, slightly
larger than the minimal one covered by the theory.

Keywords: Virtual element method, divergence free discretization, saddle-point
linear system, domain decomposition preconditioner.

##  1 Introduction

The balancing domain decomposition by constraints (BDDC) preconditioner is an
iterative substructuring method for the solution of partial differential
equations (PDEs), that belongs to the class of nonoverlapping domain
decomposition algorithms [24, 25]. BDDC, first introduced in [14] for elliptic
problems, represents an evolution of the balancing Neumann-Neumann
preconditioner [25]. We also remark that BDDC presents several features in
common with the dual-primal finite element tearing and interconnecting (FETI-
DP) algorithm. In particular, the BDDC and FETI-DP operators share almost the
same eigenvalues [22, 8], thus they exhibit analogous convergence properties.
Both BDDC and FETI-DP have been successfully developed for finite and spectral
element discretizations of several physical problems governed by PDEs, see
e.g. [19, 28, 15, 23]. In particular, regarding the Stokes equations, they
have been studied in [21, 20]. In recent years, BDDC and FETI-DP algorithms
have been also extended to various innovative discretizations techniques for
PDEs, such as Mortar discretizations [18], discontinuous Galerkin methods [16,
11], isogeometric analysis [17, 27], weak Galerkin methods [26] and virtual
element methods [4, 5].

The Virtual Element Method (VEM), introduced in the pioneering paper [2],
represents a generalization of the finite element method (FEM), that can
easily handle general polytopal meshes. The core idea behind VEM is to use
approximated discrete bilinear forms, whose computation requires only the
integration of polynomials on the element boundary and interior. The resulting
discrete solution is conforming and the accuracy guaranteed by such discrete
bilinear forms turns to be sufficient to achieve the correct order of
convergence. The advantage of these methods is that they can be applied on a
wide choice of general polygonal meshes without the need to integrate complex
non-polynomial functions on the elements, keeping an high degree of accuracy.

In the VEM literature only a few studies have focused on the construction and
analysis of preconditioners for VEM approximations of PDEs; see [1, 9, 10,
12]). BDDC for VEM discretizations of scalar elliptic problems have been first
introduced in [4, 5] and then extended to mixed formulations of scalar
elliptic equations in [13]. To our knowledge, the development of effective
non-overlapping domain decomposition preconditioners for VEM discretizations
of the Stokes equations is still an open problem.

The novelty of the present study is to develop a BDDC preconditioner for the
divergence free VEM discretization of the two-dimensional Stokes equations
introduced in [3]. Our algorithm represents an extension to VEM of the BDDC
preconditioner proposed in [21] for FEM discretizations of the Stokes
equations with discontinuous pressure spaces. We prove a convergence rate
estimate of the preconditioned system, independent of the number of subdomains
and polylogarithmic with respect to the ratio $H/h$, where $H$ denotes the
subdomain size and $h$ the mesh size. Such an estimate yields the scalability
and quasi-optimality of the resulting algorithm. Several numerical tests
confirm the theoretical estimate and show the robustness of the solver with
respect to different polygonal meshes.

The paper is organized as follows: in Section 2 we introduce the continuous
problem and its variational formulation; in Section 3 we describe the VEM
discretization; in Section 4 we introduce the domain decomposition tecnique
and the BDDC preconditioner; in Sections 5 and 6 we describe the theoretical
aspects, while in Section 7 we report several numerical results; finally in
Section 8 we draw the conclusions.

##  2 Continuous problem

Let $\Omega\subseteq\mathbb{R}^{2}$, with $\Gamma=\partial\Omega$, and
consider the stationary Stokes problem on $\Omega$ with homogeneous Dirichet
boundary conditions:

$\begin{cases}\text{Find }(\mathbf{u},p)\text{ such that}\\\
-\nu\mathbf{\Delta}\mathbf{u}-\nabla p=\mathbf{f}\qquad&\text{ in
}\mathit{\Omega}\\\ \text{div }\mathbf{u}=0\qquad&\text{ in
}\mathit{\Omega}\\\ \mathbf{u}=0\qquad&\text{ on }\mathit{\Gamma},\end{cases}$
(1)

where $\mathbf{u}$ and $p$ are the velocity and the pressure fields,
respectively. Furthermore $\mathbf{\Delta}$, div and $\mathbf{\nabla}$ denote
the vector Laplacian, the divergence and the gradient operators. Finally,
$\mathbf{f}$ represents the external force, while $\nu>0$ is the viscosity.

Let us consider the spaces:

$\displaystyle\mathbf{V}:=[H^{1}_{0}(\mathit{\Omega})]^{2},\qquad
Q:=L^{2}_{0}(\mathit{\Omega})=\bigg{\\{}q\in L^{2}(\mathit{\Omega})\quad
s.t.\quad\int_{\mathit{\Omega}}q\text{ d}\Omega=0\bigg{\\}}$ (2)

with norms:

$\displaystyle\|\mathbf{v}\|_{1}:=\|\mathbf{v}\|_{[H^{1}(\mathit{\Omega})]^{2}},\quad\|q\|_{Q}:=\|q\|_{L^{2}(\mathit{\Omega})}.$
(3)

We assume $\mathbf{f}\in[H^{-1}(\mathit{\Omega})]^{2}$, and $\nu\in
L^{\infty}(\mathit{\Omega})$ uniformly positive in $\mathit{\Omega}$. Let the
bilinear forms
$a(\cdot,\cdot):\mathbf{V}\times\mathbf{V}\rightarrow\mathbb{R}$ and
$b:\mathbf{V}\times Q\rightarrow\mathbb{R}$ be defined as:

$a(\mathbf{u},\mathbf{v}):=\int_{\Omega}\nu\mathbf{\nabla u}:\mathbf{\nabla
v}\text{ d}\Omega\qquad\text{for all }\mathbf{u},\mathbf{v}\in\mathbf{V}$ (4)
$b(\mathbf{v},q):=\int_{\Omega}\text{div }\mathbf{v}q\text{
d}\Omega\qquad\text{for all }\mathbf{u}\in\mathbf{V},q\in Q.$ (5)

Then a standard variational formulation of problem (1) reads:

$\begin{cases}\text{find }(\mathbf{u},p)\in\mathbf{V}\times Q\text{ such
that}\\\
a(\mathbf{u},\mathbf{v})+b(\mathbf{v},p)=(\mathbf{f},\mathbf{v})&\text{for all
}\mathbf{v}\in\mathbf{V},\\\ b(\mathbf{u},q)=0&\text{for all }q\in
Q,\end{cases}$ (6)

where

$\displaystyle(\mathbf{f},\mathbf{v}):=\int_{\Omega}\mathbf{f}\cdot\mathbf{v}\text{
d}\Omega.$

It is well-known that:

  * •

$a(\cdot,\cdot)$ and $b(\cdot,\cdot)$ are continuous, i.e.

$\begin{array}[]{llll}|a(\mathbf{u},\mathbf{v})|&\leq&\|a\|\|\mathbf{u}\|_{1}\|\mathbf{v}\|_{1}&\quad\text{for
all }\mathbf{u},\mathbf{v}\in\mathbf{V}\vspace{0.2cm},\\\
|b(\mathbf{v},q)|&\leq&\|b\|\|\mathbf{v}\|_{1}\|q\|_{Q}&\quad\text{for all
}\mathbf{v}\in\mathbf{V}\text{and }q\in Q,\end{array}$

where $\|a\|$ and $\|b\|$ are the usual norm of the two bilinear forms;

  * •

$a(\cdot,\cdot)$ is coercive i.e., there exists a positive constant $\alpha$
such that

$|a(\mathbf{v},\mathbf{v})|\geq\alpha\|\mathbf{v}\|_{1}^{2}\quad\text{for all
}\mathbf{v}\in\mathbf{V};$

  * •

the bilinear form $b(\cdot,\cdot)$ satisfies the inf-sup condition [6], i.e.

$\exists\beta>0\text{ such
that}\quad\sup_{\mathbf{v}\in\mathbf{V},\mathbf{v}\neq\mathbf{0}}\frac{|b(\mathbf{v},q)|}{\|\mathbf{v}\|_{1}}\geq\beta\|q\|_{Q}\quad\text{for
all }q\in Q.$ (7)

Therefore, problem (6) has a unique solution
$(\mathbf{u},p)\in\mathbf{V}\times Q$ such that

$\displaystyle||\mathbf{u}||_{1}+||p||_{Q}\leq
C||\mathbf{f}||_{H^{-1}(\mathit{\Omega})},$ (8)

where the constant $C$ depends only on $\mathit{\Omega}$ and $\nu$; see [6].

##  3 Virtual element discretization

We present here the discretization of problem (1), based on the virtual
element space introduced in [3], that is designed to solve a Stokes-like
problem element-wise. In particular we will use the reduced space presented in
section 5 of [3], that, exploiting the divergence free property of the
solution, allows to save a lot of degrees of freedom especially when the
polynomial degree $k$ is large. We recall here the definition of the local
spaces. Let $\\{{\mathcal{T}}_{h}\\}_{h}$ be a sequence of triangulations of
$\Omega$ into general polygonal elements $K$ with

$\displaystyle h_{K}:=\text{diameter}(K),\quad
h:=\sup_{K\in{\mathcal{T}}_{h}}h_{K}.$

We suppose that, for all $h$, each element $K\in{\mathcal{T}}_{h}$ satisfies
the following assumptions:

  * •

$(\mathbf{A1})$ $K$ is star-shaped with respect to a ball of radius
$\geq\gamma h_{k}$,

  * •

$(\mathbf{A2})$ the distance between any two vertices of K is $\geq ch_{K}$,

  * •

$(\mathbf{A3})$ the triangulation ${\mathcal{T}}_{h}$ is quasi-uniform, i.e.
there exist positive constants $c_{0},\,c_{1}$ such that for any two elements
$K$ and $K^{\prime}$ in ${\mathcal{T}}_{h}$ we have $c_{0}\leq
h_{K}/h_{K^{\prime}}\leq c_{1}$.

where $\gamma$ and $c$ are positive constants.

######  Remark 3.1.

These hypotheses could be weakened as in [2], for example assuming that every
$K$ is a union of a finite (and uniformly bounded) number of star-shaped
domains, each satisfying ($\mathbf{A1}$).

We also assume that the scalar viscosity field $\nu$ is piecewise constant
with respect to the decomposition ${\mathcal{T}}_{h}$, i.e. $\nu$ is constant
on each polygon $K\in{\mathcal{T}}_{h}$.

For $k\in\mathbb{N}$, let us define the spaces:

  * •

$\mathbb{P}_{k}(K)$ the set of polynomials on $K$ of degree $\leq k$,

  * •

$\mathbb{B}_{k}(K):=\\{v\in C^{0}(\partial K)\text{ s.t.
}v_{|e}\in\mathbb{P}_{k}(e)\quad\forall\text{ edge }e\in\partial K\\}$,

  * •

$\mathit{G}_{k}(K):=\nabla(\mathbb{P}_{k+1}(K))\subseteq[\mathbb{P}_{k}(K)]^{2}$,

  * •

$\mathit{G}_{k}(K)^{\perp}\subseteq[\mathbb{P}_{k}(K)]^{2}$ the
$L^{2}$-orthogonal complement to $\mathit{G}_{k}(K)$.

On each element $K\in\mathit{T}_{h}$ we define, for $k\geq 2$, the following
finite dimensional local virtual element spaces:

$\begin{split}\mathbf{\widehat{V}}_{h}^{K}:=\bigg{\\{}\mathbf{v}\in[H^{1}(K)]^{2}\text{
s.t. }\mathbf{v}_{|\partial K}\in\mathbb{B}_{k}(\partial K)]^{2},\\\
\begin{cases}-\nu\mathbf{\Delta}\mathbf{v}-\nabla
s\in\mathit{G}_{k-2}(K)^{\perp},\\\ \text{div
}\mathbf{v}\in\mathbb{P}_{0}(K),\end{cases}\text{for some }s\in
L^{2}(K)\bigg{\\}}\end{split}$ (9)

and

$\displaystyle Q_{h}^{K}:=\mathbb{P}_{0}(K).$ (10)

Now it is possible to introduce suitable sets of degrees of freedom for the
local approximations fields.
Given a function $\mathbf{v}\in\mathbf{\widehat{V}}_{h}^{K}$ we take the
following linear operators $\mathbf{D_{\widehat{V}}}$, split into three
subsets:

  * •

$\mathbf{D_{\widehat{V}}1}$: the values of $\mathbf{v}$ at the vertices of the
polygon $K$,

  * •

$\mathbf{D_{\widehat{V}}2}$: the values of $\mathbf{v}$ at $k-1$ distinct
points of every edge $e\in\partial K$ (for the implementation we will take the
$k-1$ internal points of the $(k+1)$-Gauss-Lobatto quadrature rule in $e$),

  * •

$\mathbf{D_{\widehat{V}}3}$: the moments of the values of $\mathbf{v}$

$\displaystyle\int_{K}\mathbf{v}\cdot\mathbf{g}^{\perp}_{k-2}\text{
d}K\qquad\text{for all
}\mathbf{g}^{\perp}_{k-2}\in\mathbf{G}_{k-2}(K)^{\perp}$

Furthermore, for the local pressure, given $q\in Q^{K}_{h}$, we consider the
linear operators $\mathbf{D_{Q}}$:

  * •

$\mathbf{D_{Q}}$: the moment

$\displaystyle\int_{K}q\text{ d}K.$

Since $\mathbf{D_{\widehat{V}}}$ and $\mathbf{D_{Q}}$ are unisolvent
respectively of $\mathbf{\widehat{V}}_{h}^{K}$ and $Q_{h}^{K}$, we can define
the global virtual element spaces:

$\displaystyle\mathbf{\widehat{V}}_{h}:=\\{\mathbf{v}\in[H^{1}_{0}(\mathit{\Omega})]^{2}\quad\text{s.t.}\quad\mathbf{v}_{|K}\in\mathbf{\widehat{V}}^{K}_{h}\quad\text{for
all }K\in{\mathcal{T}}_{h}\\}$ (11)

and

$\displaystyle Q_{h}:=\\{q\in L^{2}_{0}(\mathit{\Omega})\quad\text{s.t.}\quad
q_{|K}\in Q^{K}_{h}\quad\text{for all }K\in{\mathcal{T}}_{h}\\},$ (12)

with obvious associated sets of global degrees of freedom.

###  3.1 Discrete problem

Referring to [3], we can now state the discrete virtual element problem

$\begin{cases}\text{find }(\mathbf{u}_{h},p_{h})\text{ such that}\\\
a_{h}(\mathbf{u}_{h},\mathbf{v}_{h})+b(\mathbf{v}_{h},p_{h})=(\mathbf{f}_{h},\mathbf{v}_{h})&\quad\text{for
all }\mathbf{v}_{h}\in\mathbf{\widehat{V}}_{h}\\\
b(\mathbf{u}_{h},q_{h})=0&\quad\text{for all }q_{h}\in Q_{h}\end{cases}$ (13)

By construction the discrete bilinear form $a_{h}(\cdot,\cdot)$ is stable
(uniformly) with respect to the $\mathbf{V}$ norm and also obviously the
bilinear form $b(\cdot,\cdot)$. Therefore, to prove the existence and
uniqueness of the solution of the problem (13) is necessary only a suitable
inf-sup condition. For our work, we will only need this condition for the
subdomains in which $\Omega$ will be divided into. In this way the local
subdomains problem, as weel as the global one, will be well posed. The proof
of the following inf-sup condition could be found in [3].

######  Proposition 3.1.

Given the discrete spaces $\mathbf{\widehat{V}}_{h}$ and $Q_{h}$ defined in
(11) and (12), there exists a positive $\tilde{\beta}$, independent of h, such
that:

$\sup_{\mathbf{v}_{h}\in\mathbf{\widehat{V}}_{h},\mathbf{v}_{h}\neq\mathbf{0}}\frac{|b(\mathbf{v}_{h},q_{h})|}{\|\mathbf{v}_{h}\|_{1}}\geq\tilde{\beta}\|q_{h}\|_{Q}\quad\text{for
all }q_{h}\in Q_{h}.$ (14)

A consequence of the previous proposition is the following statement.

######  Theorem 3.1.

Problem (13) has a unique solution
$(\mathbf{u}_{h},p_{h})\in\mathbf{\widehat{V}}_{h}\times Q_{h}$, verifying the
estimate

$\displaystyle\|\mathbf{u}_{h}\|_{1}+\|p_{h}\|_{Q}\leq
C\|\mathbf{f}\|_{0}\text{.}$ (15)

We have also a convergence result

######  Theorem 3.2.

Let $(\mathbf{u},p)\in\mathbf{V}\times Q$ be the solution of problem (6) and
$(\mathbf{u}_{h},p_{h})\in\mathbf{\widehat{V}}_{h}\times Q_{h}$ be the
solution of problem (13). Then it holds

$\displaystyle\|\mathbf{u}-\mathbf{u}_{h}\|_{1}\leq
Ch^{k}(|\mathbf{f}|_{k-1}+|\mathbf{u}|_{k+1})$ (16)

and

$\displaystyle\|p-p_{h}\|_{Q}\leq
Ch^{k}(|\mathbf{f}|_{k-1}+|\mathbf{u}|_{k+1}+|p|_{k})\text{.}$ (17)

##  4 Construction of the BDDC preconditioner

In this section, we first divide the domain $\Omega$ into subdomains and
introduce appropriate function spaces, paragraph 4.1. Then, in paragraph 4.2,
we show how the global interface saddle-point problem takes form and then in
4.3 we define the BDDC preconditioner that allows us to use a preconditioned
conjugate gradient method (PCG) for its solution.

###  4.1 Domain decomposition

We decompose the domain $\Omega$ into $N$ non-overlapping subdomains
$\Omega_{i},i=1,2,...N$, of characteristic diameter $H$. Each subdomain is a
union of shape regular elements and the nodes on the boundaries of neighboring
subdomain match across the interface
$\Gamma=(\cup\partial\Omega_{i})\setminus\partial\Omega$; we define also
$\Gamma_{i}=\partial\Omega_{i}\cap\Gamma$ as the interface of an individual
subdomain $\Omega_{i}$. According to [5], where more details could be found,
we recall two requirements on the subdomain partition:

  * •

(S1) Each subdomain $\Omega_{i}$ is the union of polygonal elements of the
triangulation ${\mathcal{T}}_{h}$ and the number of polygons forming an
individual subdomain is uniformly bounded;

  * •

(S2) If a face of a subdomain intersects $\partial\Omega$, then the measure of
this set is comparable to that of $\partial\Omega_{i}$. Similarly, if an edge
of a subdomain intersects $\partial\Omega$, the length of this intersection is
bounded from below in terms of the diameter of $\partial\Omega_{i}$.

Figure 1: Interface of the subdomains (excluding the nodes on the boundary):
red circles indicate the vertices of the subdomains, whereas black circles
indicate the remainder interface nodes.

Restricting to the two-dimensional case, although the theory of iterative
substructuring ([25] Section 4.2) does not cover the general cases where the
boundary of a subdomain is not a straight line (as we have in our
implementation, since we use general polygonal meshes), we can anyway define
vertices and interface relatively easily. We say that a node $x$ belongs to
the interface of a subdomain if it belongs to at least two subdomains, while a
node $x$ is a vertex of a subdomain if it belongs to more than two subdomains
(Figure 1). This is the rule that we used in the implementation to split our
mesh in the different subdomains.

###  4.2 Decomposition of the virtual element spaces

The discrete variational problem (13) can be written, in matrix form, as the
following saddle-point linear system:

$\left[\begin{array}[]{cc}A&B^{T}\\\ B&0\\\
\end{array}\right]\left[\begin{array}[]{c}\mathbf{u}\\\ p\\\
\end{array}\right]=\left[\begin{array}[]{c}\mathbf{f}\\\ 0\\\
\end{array}\right]$ (18)

where the matrices A and B are associated with the discrete bilinear forms
$a_{h}(\cdot,\cdot)$ and $b(\cdot,\cdot)$. In the remainder of the paper, we
omit the underscore $h$ since we will always refer to the finite dimensional
space and so we write $\mathbf{\widehat{V}}\times Q$ instead of
$\mathbf{\widehat{V}}_{h}\times Q_{h}$, only for sake of simplifying the
notation. Referring to the notations of the previous section, we naturally
split the degrees of freedom (dofs) of the velocity components into boundary
dofs ($\mathbf{D_{V}1}$ and $\mathbf{D_{V}2}$) and interior dofs
($\mathbf{D_{V}3}$ and $\mathbf{D_{V}4}$). Following the notations introduced
in [21], we decompose the discrete velocity and pressure space
$\mathbf{\widehat{V}}$ and $Q$ into:

$\mathbf{\widehat{V}}=\mathbf{V}_{I}\bigoplus\mathbf{\widehat{V}}_{\Gamma}\text{,}\quad
Q=Q_{I}\bigoplus Q_{0}\text{.}$ (19)

$\mathbf{V}_{I}$ and $Q_{I}$ are direct sums of subdomain interior velocity
spaces $\mathbf{V}_{I}^{(i)}$, and subdomain interior pressure spaces
$Q_{I}^{(i)}$, respectively, i.e.,

$\mathbf{V}_{I}=\bigoplus_{i=1}^{N}\mathbf{V}_{I}^{(i)}\text{,}\quad
Q_{I}=\bigoplus_{i=1}^{N}Q_{I}^{(i)}\text{.}$ (20)

The elements of $\mathbf{V}_{I}^{(i)}$ have support in the subdomain
$\Omega_{i}$ and vanish on its interface $\Gamma_{i}$, while the elements of
$Q_{I}^{(i)}$ are restrictions of elements in Q to $\Omega_{i}$.
$\mathbf{\widehat{V}}_{\Gamma}$ is the space of the traces on $\Gamma$ of
functions in $\mathbf{\widehat{V}}$ and $Q_{0}$ is the subspace of $Q$ with
constant values $q_{0}^{(i)}$ in the subdomain $\Omega_{i}$. We denote the
space of interface velocity variables of the subdomain $\Omega_{i}$ by
$\mathbf{V}_{\Gamma}^{(i)}$, and the associated product space by
$\mathbf{V}_{\Gamma}=\prod_{i=1}^{N}\mathbf{V}_{\Gamma}^{(i)}$; generally
functions in $\mathbf{V}_{\Gamma}$ are discontinuous across the interface.
$R_{\Gamma}^{(i)}:\mathbf{\widehat{V}}_{\Gamma}\rightarrow\mathbf{V}_{\Gamma}^{(i)}$
is the operator which maps functions in the continuous interface velocity
space $\mathbf{\widehat{V}}_{\Gamma}$ to their subdomain components in the
space $\mathbf{V}_{\Gamma}^{(i)}$. We denote the direct sum of the
$R_{\Gamma}^{(i)}$ with $R_{\Gamma}$.
With the decomposition of the solution space given in (19), the global saddle-
point problem (18) can be written as: find
$(\mathbf{u}_{I},p_{I},\mathbf{u}_{\Gamma},p_{0})\in(\mathbf{V}_{I},Q_{I},\mathbf{\widehat{V}}_{\Gamma},Q_{0})$,
such that:

$\left[\begin{array}[]{cccc}A_{II}&B_{II}^{T}&\widehat{A}_{\Gamma I}^{T}&0\\\
B_{II}&0&\widehat{B}_{I\Gamma}&0\\\ \widehat{A}_{\Gamma
I}&\widehat{B}_{I\Gamma}^{T}&\widehat{A}_{\Gamma\Gamma}&\widehat{B}_{0\Gamma}^{T}\\\
0&0&\widehat{B}_{0\Gamma}^{T}&0\\\
\end{array}\right]\left[\begin{array}[]{c}\mathbf{u}_{I}\\\ p_{I}\\\
\mathbf{u}_{\Gamma}\\\ p_{0}\\\
\end{array}\right]=\left[\begin{array}[]{c}\mathbf{f}_{I}\\\ 0\\\
\mathbf{f}_{\Gamma}\\\ 0\\\ \end{array}\right].$ (21)

######  Remark 4.1.

Here the lower left block of (21) is zero because the bilinear form
$b(\mathbf{u}_{I},q_{0})$ vanishes for any $\mathbf{v}_{I}\in\mathbf{V}_{I}$
and $q_{0}\in Q_{0}$. To keep this property, when the change of basis for the
pressure space is applied, it is important to take care of the fact that the
shape and dimension of the elements is different.

The blocks related to the continuous interface velocity are assembled from the
corresponding subdomain submatrices, e.g.,
$\widehat{A}_{\Gamma\Gamma}=\sum_{i=1}^{N}{R_{\Gamma}^{(i)}}^{T}\widehat{A}_{\Gamma\Gamma}^{(i)}R_{\Gamma}^{(i)}$
and
$\widehat{B}_{0\Gamma}=\sum_{i=1}^{N}\widehat{B}_{0\Gamma}^{(i)}R_{\Gamma}^{(i)}$.
Correspondingly, the right-hand side vector $\mathbf{f}_{I}$ consists of
subdomain vectors $\mathbf{f}_{I}^{(i)}$, and $\mathbf{f}_{\Gamma}$ is
assembled from the subdomain components $\mathbf{f}_{\Gamma}^{(i)}$; we denote
the spaces of the right-hand side vectors $\mathbf{f}_{I}$ and
$\mathbf{f}_{\Gamma}$ by $\mathbf{F}_{I}$ and $\mathbf{F}_{\Gamma}$
respectively.

By employing a symmetric permutation, the leading two by two blocks in the
coefficient matrix can be rewritten as a block diagonal matrix with blocks
corresponding to independent subdomain problems. We show here how such a
matrix takes form in the simplest case of two subdomains:

$\left[\begin{array}[]{cccccc}A_{II}^{(1)}&{B_{II}^{(1)}}^{T}&0&0&{A_{\Gamma
I}^{(1)}}^{T}&0\\\ B_{II}^{(1)}&0&0&0&B_{I\Gamma}^{(1)}&0\\\
0&0&A_{II}^{(2)}&{B_{II}^{(2)}}^{T}&{A_{\Gamma I}^{(2)}}^{T}&0\\\
0&0&B_{II}^{(2)}&0&B_{I\Gamma}^{(2)}&0\\\ A_{\Gamma
I}^{(1)}&{B_{I\Gamma}^{(1)}}^{T}&A_{\Gamma
I}^{(2)}&{B_{I\Gamma}^{(2)}}^{T}&\widehat{A}_{\Gamma\Gamma}&{\widehat{B}_{0\Gamma}}^{T}\\\
0&0&0&0&\widehat{B}_{0\Gamma}&0\end{array}\right]\left[\begin{array}[]{c}\mathbf{u}_{I}^{(1)}\\\
p_{I}^{(1)}\\\ \mathbf{u}_{I}^{(2)}\\\ p_{I}^{(2)}\\\ \mathbf{u}_{\Gamma}\\\
p_{0}\\\ \end{array}\right]=\left[\begin{array}[]{c}\mathbf{f}_{I}^{(1)}\\\
0\\\ \mathbf{f}_{I}^{(2)}\\\ 0\\\ \mathbf{f}_{\Gamma}\\\ 0\\\
\end{array}\right].$ (22)

In the rest of this section the computations are always performed in the case
of two subdomains. The extension to the general case with more subdomains is
natural, but the computations are clearly more involved.
We proceed eliminating, by static condensation, the independent subdomain
variables $(\mathbf{u}_{I}^{(1)},p_{I}^{(1)})$ and
$(\mathbf{u}_{I}^{(2)},p_{I}^{(2)})$ in the system (22). To do so, we solve
two independent Dirichlet problems:

$\left[\begin{array}[]{cc}A_{II}^{(1)}&{B_{II}^{(1)}}^{T}\\\
B_{II}^{(1)}&0\end{array}\right]\left[\begin{array}[]{c}\mathbf{u}_{I}^{(1)}\\\
p_{I}^{(1)}\\\ \end{array}\right]+\left[\begin{array}[]{cc}{A_{\Gamma
I}^{(1)}}^{T}&0\\\
B_{I\Gamma}^{(1)}&0\end{array}\right]\left[\begin{array}[]{c}\mathbf{u}_{\Gamma}\\\
p_{0}\\\ \end{array}\right]=\left[\begin{array}[]{c}\mathbf{F}_{I}^{(1)}\\\
0\\\ \end{array}\right]\text{,}$ (23)
$\displaystyle\left[\begin{array}[]{cc}A_{II}^{(2)}&{B_{II}^{(2)}}^{T}\\\
B_{II}^{(2)}&0\end{array}\right]\left[\begin{array}[]{c}\mathbf{u}_{I}^{(2)}\\\
p_{I}^{(2)}\\\ \end{array}\right]+\left[\begin{array}[]{cc}{A_{\Gamma
I}^{(2)}}^{T}&0\\\
B_{I\Gamma}^{(2)}&0\end{array}\right]\left[\begin{array}[]{c}\mathbf{u}_{\Gamma}\\\
p_{0}\\\ \end{array}\right]=\left[\begin{array}[]{c}\mathbf{F}_{I}^{(2)}\\\
0\\\ \end{array}\right]\text{,}$ (34)

thus

$\left[\begin{array}[]{c}\mathbf{u}_{I}^{(1)}\\\ p_{I}^{(1)}\\\
\end{array}\right]=\left[\begin{array}[]{cc}A_{II}^{(1)}&{B_{II}^{(1)}}^{T}\\\
B_{II}^{(1)}&0\end{array}\right]^{-1}\left(\left[\begin{array}[]{c}\mathbf{F}_{I}^{(1)}\\\
0\\\ \end{array}\right]-\left[\begin{array}[]{cc}{A_{\Gamma I}^{(1)}}^{T}&0\\\
B_{I\Gamma}^{(1)}&0\end{array}\right]\left[\begin{array}[]{c}\mathbf{u}_{\Gamma}\\\
p_{0}\\\ \end{array}\right]\right)\text{,}$ (35)
$\left[\begin{array}[]{c}\mathbf{u}_{I}^{(2)}\\\ p_{I}^{(2)}\\\
\end{array}\right]=\left[\begin{array}[]{cc}A_{II}^{(2)}&{B_{II}^{(2)}}^{T}\\\
B_{II}^{(2)}&0\end{array}\right]^{-1}\left(\left[\begin{array}[]{c}\mathbf{F}_{I}^{(2)}\\\
0\\\ \end{array}\right]-\left[\begin{array}[]{cc}{A_{\Gamma I}^{(2)}}^{T}&0\\\
B_{I\Gamma}^{(2)}&0\end{array}\right]\left[\begin{array}[]{c}\mathbf{u}_{\Gamma}\\\
p_{0}\\\ \end{array}\right]\right)\text{,}$ (36)

Then, substituting the solutions of (35) and (36) in

$\left[\begin{array}[]{cc}A_{\Gamma I}^{(1)}&{B_{I\Gamma}^{(1)}}^{T}\\\
0&0\end{array}\right]\left[\begin{array}[]{c}\mathbf{u}_{I}^{(1)}\\\
p_{I}^{(1)}\\\ \end{array}\right]+\left[\begin{array}[]{cc}A_{\Gamma
I}^{(2)}&{B_{I\Gamma}^{(2)}}^{T}\\\
0&0\end{array}\right]\left[\begin{array}[]{c}\mathbf{u}_{I}^{(2)}\\\
p_{I}^{(2)}\\\
\end{array}\right]+\left[\begin{array}[]{cc}\widehat{A}_{\Gamma\Gamma}&{\widehat{B}_{0\Gamma}}^{T}\\\
\widehat{B}_{0\Gamma}&0\end{array}\right]\left[\begin{array}[]{c}\mathbf{u}_{\Gamma}\\\
p_{0}\\\ \end{array}\right]=\left[\begin{array}[]{c}\mathbf{F}_{\Gamma}\\\
0\\\ \end{array}\right]$ (37)

we obtain the global interface saddle-point problem:

$\widehat{S}\text{
}\widehat{u}=\left[\begin{array}[]{cc}\widehat{S}_{\Gamma}&{\widehat{B}_{0\Gamma}}^{T}\\\
\widehat{B}_{0\Gamma}&0\\\
\end{array}\right]\left[\begin{array}[]{c}\mathbf{u}_{\Gamma}\\\ p_{0}\\\
\end{array}\right]=\left[\begin{array}[]{c}\mathbf{g}_{\Gamma}\\\ 0\\\
\end{array}\right]=\widehat{\mathbf{g}}\text{,}$ (38)

where the right-hand side $\widehat{\mathbf{g}}\in\mathbf{F}_{\Gamma}\times
F_{0}$ is given by

$\displaystyle\widehat{\mathbf{g}}=\sum_{i=1}^{2}{R_{\Gamma}^{(i)}}^{T}\bigg{\\{}\left[\begin{array}[]{c}\mathbf{f}_{\Gamma}^{(i)}\\\
0\\\ \end{array}\right]-\left[\begin{array}[]{cc}A_{\Gamma
I}^{(i)}&{B_{I\Gamma}^{(i)}}^{T}\\\
0&0\end{array}\right]\left[\begin{array}[]{cc}A_{II}^{(i)}&{B_{II}^{(i)}}^{T}\\\
B_{II}^{(i)}&0\end{array}\right]^{-1}\left[\begin{array}[]{c}\mathbf{f}_{I}^{(i)}\\\
0\\\ \end{array}\right]\bigg{\\}}.$ (47)

We note that $\widehat{S}$ is assembled from the subdomain Stokes Schur
complements $S^{(i)}$, which are defined by: given
${\mathbf{w}}^{(i)}=\mathbf{w}_{\Gamma}^{(i)}\times
q_{0}^{(i)}\in\mathbf{V}_{\Gamma}^{(i)}\times Q_{0}^{(i)}$, determine
$S^{(i)}{\mathbf{w}}^{(i)}\in\mathbf{F}_{\Gamma}^{(i)}\times F_{0}^{(i)}$ such
that

$\left[\begin{array}[]{cccc}A_{II}^{(i)}&{B_{II}^{(i)}}^{T}&{A_{\Gamma
I}^{(i)}}^{T}&0\\\ B_{II}^{(i)}&0&B_{I\Gamma}^{(i)}&0\\\ A_{\Gamma
I}^{(i)}&{B_{I\Gamma}^{(i)}}^{T}&{A}_{\Gamma\Gamma}^{(i)}&{{B}_{0\Gamma}^{(i)}}^{T}\\\
0&0&{B}_{0\Gamma}^{(i)}&0\end{array}\right]\left[\begin{array}[]{c}\mathbf{w}_{I}^{(i)}\\\
q_{I}^{(i)}\\\ \mathbf{w}_{\Gamma}^{(i)}\\\
q_{0}^{(i)}\end{array}\right]=\left[\begin{array}[]{c}\mathbf{0}\\\ 0\\\
S^{(i)}{\mathbf{w}}^{(i)}\end{array}\right]\text{.}$ (48)

Denoting by $S_{\Gamma}$ the direct sum of the $S^{(i)}_{\Gamma}$, then
$\widehat{S}_{\Gamma}$ is given by

$\widehat{S}_{\Gamma}=R_{\Gamma}^{T}S_{\Gamma}R_{\Gamma}=\sum_{i=1}^{2}{R^{(i)}_{\Gamma}}^{T}S^{(i)}_{\Gamma}R^{(i)}_{\Gamma},$
(49)

and then we set

$R=\left[\begin{array}[]{cc}R_{\Gamma}&0\\\ 0&I\end{array}\right]\text{,}\quad
R^{(i)}=\left[\begin{array}[]{cc}R_{\Gamma}^{(i)}&0\\\
0&I\end{array}\right]\text{.}$ (50)

Finally we see from (48), that the action of $S^{(i)}$ on a vector can be
evaluated by solving a Dirichlet problem on the subdomain $\Omega_{i}$ as in
(35) and (36), so it is not necessary to assemble the matrix $\widehat{S}$
because only its action is required.
In the next section we introduce a BDDC preconditioner for problem (38), where
the operator of the preconditioned problem is symmetric and positive definite,
so we will use the PCG method to solve it.

###  4.3 BDDC preconditioner

We now present the BDDC preconditioner, first designed in [21] for finite
element discretizations of the Stokes equations, that we will extend to the
VEM discretization introduced in the previous sections. This preconditioner is
very similar to FETI-DP, but there is a main difference between them: while in
a FETI-DP algorithm the continuity of the solution will not be fully satisfied
until the algorithm has converged, in the BDDC one full continuity is restored
at the end of each iteration step, by using an average operator.
Before entering into the definition of the function space used to construct
the BDDC preconditioner, we briefly justify the choice of our notation. The
subscript $\Gamma$ indicates dofs living on the interface, $\Pi$ and $\Delta$
are instead used to distinguish dofs of $\Gamma$ that belong to the primal and
dual spaces, respectively, defined here below. Two other subscripts are used:
$C$ indicates an operator referred to the coarse space and $D$ is instead used
to highlight that an operator has been rescaled by suitable scaling functions,
defined later. The hat $\widehat{\cdot}$ refers to a continuous space, the
$\widetilde{\cdot}$ means that the space is continuous on primal interface
dofs and discontinuous on the dual ones and finally no hat is used for the
product of local spaces, which is discontinuous at all interface dofs.
As a first step, we introduce a partially assembled interface velocity space
$\mathbf{\widetilde{V}}_{\Gamma}$,

$\mathbf{\widetilde{V}}_{\Gamma}=\mathbf{\widehat{V}}_{\Pi}\bigoplus\mathbf{V}_{\Delta}=\mathbf{\widehat{V}}_{\Pi}\bigoplus\big{(}\prod_{i=1}^{N}\mathbf{V}_{\Delta}^{(i)}\big{)}.$
(51)

$\mathbf{\widehat{V}}_{\Pi}$ is the continuous coarse level primal interface
velocity space which typically is spanned by subdomain vertex nodal basis
functions, and/or by interface edge basis functions with constant values, or
with values of weight functions, on these edge. These basis functions
correspond to the primal interface velocity continuity constraints, which will
be discussed later. We will always assume that the basis has been changed so
that each primal basis function corresponds to an explicit degree of freedom.
In other words, we will have explicit primal unknowns corresponding to the
primal continuity constraints on edges. The primal degrees of freedom are
shared by neighboring subdomains. The complimentary space
$\mathbf{V}_{\Delta}$ is the direct sum of the subdomain dual interface
velocity spaces $\mathbf{V}_{\Delta}^{(i)}$ , which correspond to the
remaining interface velocity degrees of freedom and are spanned by basis
functions which vanish at the primal degrees of freedom. Thus, an element in
the space $\mathbf{\widetilde{V}}_{\Gamma}$ has a continuous primal velocity
and typically a discontinuous dual velocity component.
We now introduce several restriction, extension, and scaling operators between
a variety of spaces. As in [21], $R_{\Gamma}^{(i)}$ is the operator which maps
a function in the space $\mathbf{\widehat{V}}_{\Gamma}$ to its component in
$\mathbf{V}_{\Gamma}^{(i)}$. We define $R_{\Delta}^{(i)}$ as the operator
which maps the space $\mathbf{\widehat{V}}_{\Gamma}$ to its dual component in
the space $\mathbf{V}_{\Delta}^{(i)}$. $R_{\Gamma\Pi}$ is the restriction
operator from the space $\mathbf{\widehat{V}}_{\Gamma}$ to its subspace
$\mathbf{\widehat{V}}_{\Pi}$; $R_{\Pi}^{(i)}$ is the operator which maps
$\mathbf{\widehat{V}}_{\Pi}$ into its $\Gamma_{i}$-component.
$\widetilde{R}_{\Gamma}$ is the direct sum of $R_{\Gamma\Pi}$ and the
$R_{\Delta}^{(i)}$ , and it is a map from $\mathbf{\widehat{V}}_{\Gamma}$ into
$\mathbf{\widetilde{V}}_{\Gamma}$.
The relationships among the previous spaces and operators are summarized in
the following diagram:

${\mathbf{\widetilde{V}}_{\Gamma}}$${\mathbf{V}_{\Gamma}^{(i)}}$${\mathbf{\widehat{V}}_{\Gamma}}$${\mathbf{V}_{\Delta}^{(i)}}$${\mathbf{\widehat{V}}_{\Pi}}$${\mathbf{\widehat{V}}_{\Pi}^{(i)}}$$\scriptstyle{R_{\Gamma}^{(i)}}$$\scriptstyle{R_{\Delta}^{(i)}}$$\scriptstyle{R_{\Pi\Gamma}}$$\scriptstyle{\widetilde{R}_{\Gamma}}$$\scriptstyle{R_{\Pi}^{(i)}}$

In order to define certain scaling operators, which will be used in the
definition of the BDDC preconditioner, see (80) , we introduce a positive
scaling factor $\delta^{\dagger}_{i}(x)$ for the nodes on the interface
$\Gamma_{i}$ of each subdomain $\Omega_{i}$. For the type of problem we will
use in the numerical experiment (incompressible Stokes problems), we simply
define the $\delta^{\dagger}_{i}(x)$ as the pseudoinverse counting functions,
so:

$\delta^{\dagger}_{i}(x):=1/card(I_{x}),\quad x\in\Gamma_{i}$ (52)

where $I_{x}$ is the set of indices of subdomains which have $x$ on their
boundaries and $card(I_{x})$ is the number of these subdomains. Now we can
define the scaled restriction operators $R^{(i)}_{D,\Delta}$, simply
multiplying each non-zero element of $R^{(i)}_{\Delta}$, only one for row, by
the corresponding scaling factor $\delta^{\dagger}_{i}(x)$. We construct also
the scaled operator $\widetilde{R}_{D,\Gamma}$ as the direct sum of
$R_{\Gamma,\Pi}$ and $R^{(i)}_{D,\Delta}$. After the change of basis, the
interface velocity Schur complement $\widetilde{S}_{\Gamma}$ is defined on the
partially assembled interface velocity space $\widetilde{\mathbf{V}}_{\Gamma}$
by: given $\mathbf{v}_{\Gamma}\in\widetilde{\mathbf{V}}_{\Gamma}$,
$\widetilde{S}_{\Gamma}\mathbf{v}_{\Gamma}\in\mathbf{\widetilde{F}}_{\Gamma}$
satisfies

$\displaystyle\left[\begin{array}[]{ccccc}A_{II}^{(i)}&B_{II}^{{(i)}^{T}}&A_{\Delta
I}^{{(i)}^{T}}&&\widetilde{A}_{\Pi I}^{{(i)}^{T}}\\\
B_{II}^{(i)}&0&B_{I\Delta}^{(i)}&&\widetilde{B}_{I\Pi}^{(i)}\\\ A_{\Delta
I}^{(i)}&B_{I\Delta}^{{(i)}^{T}}&A_{\Delta\Delta}^{(i)}&&\widetilde{A}_{\Pi\Delta}^{{(i)}^{T}}\\\
&&&\ddots&\vdots\\\ \widetilde{A}_{\Pi
I}^{(i)}&\widetilde{B}_{I\Pi}^{{(i)}^{T}}&\widetilde{A}_{\Pi\Delta}^{(i)}&\ldots&\widetilde{A}_{\Pi\Pi}\end{array}\right]\left[\begin{array}[]{c}\mathbf{v}_{I}^{(i)}\\\
p_{I}^{(i)}\\\ \mathbf{v}_{\Delta}^{(i)}\\\ \vdots\\\
\mathbf{v}_{\Pi}^{(i)}\\\
\end{array}\right]=\left[\begin{array}[]{c}\mathbf{0}\\\ 0\\\
(\widetilde{S}_{\Gamma}\mathbf{v}_{\Gamma})_{\Delta}^{(i)}\\\ \vdots\\\
(\widetilde{S}_{\Gamma}\mathbf{v}_{\Gamma})_{\Pi}^{(i)}\\\ \end{array}\right]$
(68)

Here
$\widetilde{A}_{\Pi\Pi}=\sum_{i=1}^{N}{R_{\Pi}^{(i)}}^{T}A_{\Pi\Pi}^{(i)}R_{\Pi}^{(i)}$,
$\widetilde{A}_{\Pi I}^{(i)}={R_{\Pi}^{(i)}}^{T}A_{\Pi I}^{(i)}$,
$\widetilde{A}_{\Pi\Delta}^{(i)}={R_{\Pi}^{(i)}}^{T}A_{\Pi\Delta}^{(i)}$ and
$\widetilde{B}_{I\Pi}^{(i)}=B_{I\Pi}^{(i)}R_{\Pi}^{(i)}$.

Defining by $\overline{R}_{\Gamma}$ the operator that maps the space
$\widetilde{\mathbf{V}}_{\Gamma}$ into the product space $\mathbf{V}_{\Gamma}$
associated with the set of subdomains, we observe that
$\widetilde{S}_{\Gamma}$ can be obtained from the Schur complements
$S^{(i)}_{\Gamma}$ by assembling only the primal interface velocity part, i.e.
as

$\widetilde{S}_{\Gamma}={\overline{R}_{\Gamma}}^{T}S_{\Gamma}\overline{R}_{\Gamma}.$
(69)

As we saw before (48) the global interface Schur operator
$\widehat{S}_{\Gamma}$ is obtanied by fully assembling the $S_{\Gamma}^{(i)}$
across the subdomain interface, therefore it can be also obtained from
$\widetilde{S}_{\Gamma}$ by further assembling the dual interface velocity
part,
$\widehat{S}_{\Gamma}={\widetilde{R}_{\Gamma}}^{T}\widetilde{S}_{\Gamma}\widetilde{R}_{\Gamma}$.
So we need to define an operator $\widetilde{B}_{0\Gamma}$, which maps the
partially assembled interface velocity space $\widetilde{V}_{\Gamma}$ into
$F_{0}$, the space of right hand sides corresponding to $Q_{0}$, and it is
obtained from $\widetilde{B}_{0\Gamma}$ by assembling the dual interface
velocity part on the subdomain interfaces, i.e.
$\widehat{B}_{\Gamma}=\widetilde{B}_{\Gamma}\widetilde{R}_{\Gamma}$.
Introducing

$\displaystyle\widetilde{R}=\left[\begin{array}[]{cc}\widetilde{R}_{\Gamma}&0\\\
0&I\end{array}\right]\text{,}\quad\widetilde{S}=\left[\begin{array}[]{cc}\widetilde{S}_{\Gamma}&\widetilde{B}_{0\Gamma}^{T}\\\
\widetilde{B}_{0\Gamma}&0\end{array}\right]\text{,}$ (74)

we can write $\widehat{S}$, the operator of the global interface problem (38),
as

$\displaystyle\widehat{S}=\left[\begin{array}[]{cc}\widehat{S}_{\Gamma}&{\widehat{B}_{0\Gamma}}^{T}\\\
\widehat{B}_{0\Gamma}&0\end{array}\right]=\left[\begin{array}[]{cc}{\widetilde{R}_{\Gamma}}^{T}\widetilde{S}_{\Gamma}\widetilde{R}_{\Gamma}&{\widetilde{R}_{\Gamma}}^{T}{\widetilde{B}_{0\Gamma}}^{T}\\\
\widetilde{B}_{0\Gamma}\widetilde{R}_{\Gamma}&0\end{array}\right]=\widetilde{R}^{T}\widetilde{S}\widetilde{R}\text{.}$
(79)

The preconditioner for solving the global saddle-point problem (38) is

$M^{-1}=\widetilde{R}_{D}^{T}\widetilde{S}^{-1}\widetilde{R}_{D},$ (80)

where we have defined

$\displaystyle\widetilde{R}_{D}:=\left[\begin{array}[]{cc}\widetilde{R}_{D,\Gamma}&0\\\
0&I\end{array}\right]\text{,}$ (83)

and so we have the BDDC preconditioned problem: find
$(\mathbf{u_{\Gamma}},p_{0})\in\mathbf{\widehat{V}}_{\Gamma}\times Q_{0}$,
such that

$\displaystyle\widetilde{R}_{D}^{T}\widetilde{S}^{-1}\widetilde{R}_{D}\widehat{S}\left[\begin{array}[]{c}\mathbf{u}_{\Gamma}\\\
p_{0}\\\
\end{array}\right]=\widetilde{R}_{D}^{T}\widetilde{S}^{-1}\widetilde{R}_{D}\left[\begin{array}[]{c}\mathbf{g}_{\Gamma}\\\
0\\\ \end{array}\right]\text{.}$ (88)

What we need in our implementation is to determine the action
$\widetilde{S}^{-1}\mathbf{q}$ for any given
$\mathbf{q}=(\mathbf{q_{\Gamma}},q_{0})\in\widetilde{\mathbf{F}}_{\Gamma}\times
F_{0}$, so we have to solve the linear system

$\left[\begin{array}[]{cc}\widetilde{S}_{\Gamma}&\widetilde{B}_{0\Gamma}^{T}\\\
\widetilde{B}_{0\Gamma}&0\\\
\end{array}\right]\left[\begin{array}[]{c}\mathbf{u}_{\Gamma}\\\ p_{0}\\\
\end{array}\right]=\left[\begin{array}[]{c}\mathbf{q}_{\Gamma}\\\ q_{0}\\\
\end{array}\right]\text{.}$ (89)

Given the definition of $\widetilde{S}_{\Gamma}$ in (5.3), we have that
solving (5.10) is equivalent to solve

$\displaystyle\left[\begin{array}[]{cccccc}A_{II}^{(i)}&B_{II}^{{(i)}^{T}}&A_{\Delta
I}^{{(i)}^{T}}&&\widetilde{A}_{\Pi I}^{{(i)}^{T}}&0\\\
B_{II}^{(i)}&0&B_{I\Delta}^{(i)}&&\widetilde{B}_{I\Pi}^{(i)}&0\\\ A_{\Delta
I}^{(i)}&B_{I\Delta}^{{(i)}^{T}}&A_{\Delta\Delta}^{(i)}&&\widetilde{A}_{\Pi\Delta}^{{(i)}^{T}}&B_{0\Delta}^{{(i)}^{T}}\\\
&&&\ddots&\vdots&\\\ \widetilde{A}_{\Pi
I}^{(i)}&\widetilde{B}_{I\Pi}^{{(i)}^{T}}&\widetilde{A}_{\Pi\Delta}^{(i)}&\ldots&\widetilde{A}_{\Pi\Pi}&\widetilde{B}_{0\Pi}^{T}\\\
0&0&B_{0\Delta}^{(i)}&&\widetilde{B}_{0\Pi}&0\\\
\end{array}\right]\left[\begin{array}[]{c}\mathbf{u}_{I}^{(i)}\\\
p_{I}^{(i)}\\\ \mathbf{u}_{\Delta}^{(i)}\\\ \vdots\\\ \mathbf{u}_{\Pi}\\\
p_{0}\\\ \end{array}\right]=\left[\begin{array}[]{c}\mathbf{0}\\\ 0\\\
\mathbf{q}_{\Delta}^{(i)}\\\ \vdots\\\ \mathbf{q}_{\Pi}\\\ q_{0}\\\
\end{array}\right]$ (108)

where $\widetilde{B}_{0\Pi}=\sum_{i=1}^{N}B_{0\Pi}^{(i)}R_{\Pi}^{(i)}$. Now
using a block factorization we obtain

$\displaystyle\widetilde{S}^{-1}=\sum_{i=1}^{N}\left[\begin{array}[]{ccc}0&0&R_{\Delta,i}^{T}\end{array}\right]\left[\begin{array}[]{ccc}A_{II}^{(i)}&B_{II}^{{(i)}^{T}}&A_{\Delta
I}^{{(i)}^{T}}\\\ B_{II}^{(i)}&0&B_{I\Delta}^{(i)}\\\ A_{\Delta
I}^{(i)}&B_{I\Delta}^{{(i)}^{T}}&A_{\Delta\Delta}^{(i)}\\\
\end{array}\right]^{-1}\left[\begin{array}[]{c}0\\\ 0\\\
R_{\Delta,i}\end{array}\right]+\Phi S_{CC}^{-1}\Phi^{T}\text{,}$ (116)

where $R_{\Delta,i}$ maps $\widetilde{\mathbf{F}}_{\Gamma}\times F_{0}$ into
$\mathbf{F}_{\Delta}^{(i)}$, the set of right hand sides corresponding to
$\mathbf{V}_{\Delta}^{(i)}$. The matrix $S_{CC}$, relatively to the primal
constraints, has to be completely assembled in this way

$\displaystyle\begin{split}S_{CC}=\sum_{i=1}^{N}R_{C}^{{(i)}^{T}}\bigg{\\{}\left[\begin{array}[]{cc}A_{\Pi\Pi}^{(i)}&B_{0\Pi}^{{(i)}^{T}}\\\
B_{0\Pi}^{(i)}&0\end{array}\right]-\left[\begin{array}[]{ccc}A_{\Pi
I}^{(i)}&B_{I\Pi}^{{(i)}^{T}}&A_{\Pi\Delta}^{(i)}\\\
0&0&B_{0\Delta}^{(i)}\end{array}\right]\\\
\left[\begin{array}[]{ccc}A_{II}^{(i)}&B_{II}^{{(i)}^{T}}&A_{\Delta
I}^{{(i)}^{T}}\\\ B_{II}^{(i)}&0&B_{I\Delta}^{(i)}\\\ A_{\Delta
I}^{(i)}&B_{I\Delta}^{{(i)}^{T}}&A_{\Delta\Delta}^{(i)}\\\
\end{array}\right]^{-1}\left[\begin{array}[]{cc}A_{\Pi I}^{{(i)}^{T}}&0\\\
B_{I\Pi}^{(i)}&0\\\ A_{\Pi\Delta}^{{(i)}^{T}}&B_{0\Delta}^{{(i)}^{T}}\\\
\end{array}\right]\bigg{\\}}R_{C}^{(i)}\text{,}\end{split}$ (117)

where we have defined

$\displaystyle R_{C}^{(i)}:=\left[\begin{array}[]{cc}R_{\Pi}^{(i)}&0\\\
0&I\end{array}\right]\text{,}$ (120)

the maps from $\widehat{\mathbf{V}}_{\Pi}\times Q_{0}$ to
${\mathbf{V}}_{\Pi}^{(i)}\times Q_{0}$. Finally we define the matrix

$\displaystyle\Phi=R_{\Pi
0}^{T}-\sum_{i=1}^{N}\left[\begin{array}[]{ccc}0&0&R_{\Delta,i}^{T}\end{array}\right]\left[\begin{array}[]{ccc}A_{II}^{(i)}&B_{II}^{{(i)}^{T}}&A_{\Delta
I}^{{(i)}^{T}}\\\ B_{II}^{(i)}&0&B_{I\Delta}^{(i)}\\\ A_{\Delta
I}^{(i)}&B_{I\Delta}^{{(i)}^{T}}&A_{\Delta\Delta}^{(i)}\\\
\end{array}\right]^{-1}\left[\begin{array}[]{cc}A_{\Pi I}^{{(i)}^{T}}&0\\\
B_{I\Pi}^{(i)}&0\\\ A_{\Pi\Delta}^{{(i)}^{T}}&B_{0\Delta}^{{(i)}^{T}}\\\
\end{array}\right]R_{C}^{(i)}\text{,}$ (128)

where $R_{\Pi 0}$ is the map between the space
$\widetilde{\mathbf{F}}_{\Gamma}\times F_{0}$ and
$\widehat{\mathbf{F}}_{\Pi}\times F_{0}$.

##  5 Theoretical estimates

We now present an estimate for the eigenvalues of the preconditioned operator
$M^{-1}\widehat{S}$, following the theory developed in [21] and adapting it to
our VEM formulation. We can do so because the space $\mathbf{V}_{\Gamma}$
coincides with the analogous space that would be obtained applying the same
procedure with the FEM. This substantially allow us to carry over the theory
formulated for the FEM to the VEM, except for the second assumption we will
see later. For this result, a different proof is necessary and we follow [4],
where a proof independent of the tassellation is given.

We have, as a consequence of a result on the inertia of Schur complements, the
following:

######  Lemma 5.1.

The subdomain Schur complements $S_{\Gamma}^{(i)}$, defined in (38), are
symmetric and positive definite.

###### Proof.

We know from (48) that the Schur complement related to the velocity is defined
by: given $\mathbf{w}_{\Gamma}^{(i)}\in\mathbf{V}_{\Gamma}^{(i)}$, determine
$S_{\Gamma}^{(i)}{\mathbf{w}_{\Gamma}}^{(i)}\in\mathbf{F}_{\Gamma}^{(i)}$ such
that

$\left[\begin{array}[]{ccc}A_{II}^{(i)}&{B_{II}^{(i)}}^{T}&{A_{\Gamma
I}^{(i)}}^{T}\\\ B_{II}^{(i)}&0&B_{I\Gamma}^{(i)}\\\ A_{\Gamma
I}^{(i)}&{B_{I\Gamma}^{(i)}}^{T}&{A}_{\Gamma\Gamma}^{(i)}\\\
\end{array}\right]\left[\begin{array}[]{c}\mathbf{w}_{I}^{(i)}\\\
q_{I}^{(i)}\\\
\mathbf{w}_{\Gamma}^{(i)}\end{array}\right]=\left[\begin{array}[]{c}\mathbf{0}\\\
0\\\ S_{\Gamma}^{(i)}{\mathbf{w}_{\Gamma}}^{(i)}\end{array}\right]\text{.}$
(129)

By the coercivity of $a(\cdot,\cdot)$ we know that the matrices

$\displaystyle\left[\begin{array}[]{cc}A_{II}^{(i)}&{A_{\Gamma
I}^{(i)}}^{T}\\\ A_{\Gamma I}^{(i)}&{A}_{\Gamma\Gamma}^{(i)}\\\
\end{array}\right]$

are symmetric and positive definite and so the left two by two upper block of
the left-hand-side of (129) has the same number of negative eigenvalues of the
all matrix. Now, the left-hand-side matrices of (129) are congruent to:

$\displaystyle\left[\begin{array}[]{cccc}A_{II}^{(i)}&{B_{II}^{(i)}}^{T}&0\\\
B_{II}^{(i)}&0&0\\\ 0&0&S_{\Gamma}^{(i)}\\\ \end{array}\right]$

and so, by the Sylvester’s law of inertia the velocity Schur complements are
positive definite. ∎

In the following, we denote by $a^{(i)}$, $a_{h}^{(i)}$ and $b^{(i)}$ the
restrictions to subdomain $\Omega_{i}$ of the bilinear forms $a$, $a_{h}$ and
$b$, respectively. Then, we introduce the $|.|_{S_{\Gamma}^{(i)}}$ and
$|.|_{S_{\Gamma}}$ seminorms defined by

$\displaystyle|\mathbf{v}_{\Gamma}^{(i)}|_{S_{\Gamma}^{(i)}}^{2}={\mathbf{v}_{\Gamma}^{(i)}}^{T}S_{\Gamma}^{(i)}\mathbf{v}_{\Gamma}^{(i)}\text{,}\quad|\mathbf{v}_{\Gamma}|_{S_{\Gamma}}^{2}={\mathbf{v}_{\Gamma}}^{T}S_{\Gamma}\mathbf{v}_{\Gamma}=\sum_{i=1}^{N}|\mathbf{v}_{\Gamma}^{(i)}|_{S_{\Gamma}^{(i)}}^{2},$
(130)

and a norm and a seminorm on the space $\mathbf{V}_{\Gamma}^{(i)}$

$\|\mathbf{v}_{\Gamma}^{(i)}\|_{1/2,\Gamma_{i}}^{2}=\|\mathbf{v}_{\Gamma}^{(i)}\|_{[H^{1/2}(\partial\Omega_{i})]^{2}}^{2}\text{,}\qquad|\mathbf{v}_{\Gamma}^{(i)}|_{1/2,\Gamma_{i}}^{2}=|\mathbf{v}_{\Gamma}^{(i)}|_{[H^{1/2}(\partial\Omega_{i})]^{2}}^{2}\text{,}$
(131)

with consequently the norm $\|.\|_{1/2,\Gamma}$ and seminorm
$|.|_{1/2,\Gamma}$ defined on the space $\mathbf{V}_{\Gamma}$ by
$\|\mathbf{v}_{\Gamma}\|_{1/2,\Gamma}^{2}=\sum_{i=1}^{N}\|\mathbf{v}_{\Gamma}^{(i)}\|_{1/2,\Gamma_{i}}^{2}$
and
$|\mathbf{v}_{\Gamma}|_{1/2,\Gamma}^{2}=\sum_{i=1}^{N}|\mathbf{v}_{\Gamma}^{(i)}|_{1/2,\Gamma_{i}}^{2}$.

######  Lemma 5.2.

There exist positive constant $c_{1}$ and $c_{2}$, independent of $H$, $h$ and
the shape of subdomains, such that

$c_{1}\tilde{\beta}^{2}|\mathbf{v}_{\Gamma}|_{S_{\Gamma}}^{2}\leq|\mathbf{v}_{\Gamma}|_{1/2,\Gamma}^{2}\leq
c_{2}|\mathbf{v}_{\Gamma}|_{S_{\Gamma}}^{2}\quad\forall\mathbf{v}_{\Gamma}\in\mathbf{V}_{\Gamma}\text{,}$

where $\tilde{\beta}$ is the inf-sup stability constant defined in (14).

###### Proof.

This proof follows substantially the result presented in Bramble and Pasciak
([7] Theorem 4.1) where a proof for FEM is provided. Given
$\mathbf{v}_{\Gamma}\in\mathbf{V}_{\Gamma}$, we define the operators
$T:\mathbf{V}_{\Gamma}\rightarrow\mathbf{\widehat{V}}$ and
$S:\mathbf{V}_{\Gamma}\rightarrow Q$ satisfying $\forall i=1,...,N$:

$\displaystyle\begin{split}&(1)\quad S(\mathbf{v}_{\Gamma})|_{\Omega_{i}}\in
Q_{I}^{(i)},\\\ &(2)\quad
T(\mathbf{v}_{\Gamma})|_{\Gamma}=\mathbf{v}_{\Gamma},\\\ &(3)\quad
a_{h}^{(i)}(T(\mathbf{v}_{\Gamma}),\mathbf{v})+b^{(i)}(\mathbf{v},S(\mathbf{v}_{\Gamma})=\mathbf{0}\quad\text{for
all}\quad\mathbf{v}\in\mathbf{V}_{I}^{(i)},\\\ &(4)\quad
b^{(i)}(T(\mathbf{v}_{\Gamma}),q)=0\quad\text{for all}\quad q\in
Q_{I}^{(i)}.\\\ \end{split}$ (132)

The above condition uniquely defines $S$ and $T$. Now given
$\mathbf{v}_{\Gamma}\in\mathbf{V}_{\Gamma}$, let
$\mathbf{v}_{\Gamma}^{\mathcal{H}}\in\mathbf{\widehat{V}}$ be the discrete
harmonic extension of $\mathbf{v}_{\Gamma}$, i.e. the unique function in
$\mathbf{\widehat{V}}$ which equals $\mathbf{v}_{\Gamma}$ on $\Gamma$ and
satisfies $\forall i=1,...,N$:

$\displaystyle
a^{(i)}(\mathbf{v}_{\Gamma}^{\mathcal{H}},\mathbf{v})=0\quad\text{for
all}\quad\mathbf{v}\in\mathbf{V}_{I}^{(i)}.$ (133)

By the stability of the discrete harmonic extension and the stability of the
discrete bilinear form $a_{h}$ [3], we have on each subdomain:

$\displaystyle
a_{h}^{(i)}(\mathbf{v}_{\Gamma}^{\mathcal{H}},\mathbf{v}_{\Gamma}^{\mathcal{H}})\leq
c_{3}a^{(i)}(\mathbf{v}_{\Gamma}^{\mathcal{H}},\mathbf{v}_{\Gamma}^{\mathcal{H}})\leq
c_{3}|\mathbf{v}^{(i)}_{\Gamma}|^{2}_{1/2,\Gamma_{i}},$ (134)

where $c_{3}$ is a positive constant independent of $h$, $H$ and the number of
subdomains $N$. Now, by definition of $S$ and $T$, and since
$b^{(i)}(T(\mathbf{v}_{\Gamma}),S(\mathbf{v}_{\Gamma}))=0$, we have:

$\displaystyle
a_{h}^{(i)}(T(\mathbf{v}_{\Gamma}),T(\mathbf{v}_{\Gamma}))=a_{h}^{(i)}(T(\mathbf{v}_{\Gamma}),\mathbf{v}_{\Gamma}^{\mathcal{H}})+b^{(i)}(\mathbf{v}_{\Gamma}^{\mathcal{H}},S(\mathbf{v}_{\Gamma})).$
(135)

Applying (14) on the subdomains, we have, for some $c>0$:

$\displaystyle\begin{split}||S(\mathbf{v}_{\Gamma})||^{2}_{Q_{i}}\leq\tilde{\beta}^{-2}\sup_{\mathbf{w}\in\mathbf{V}_{I}^{(i)}}\frac{b^{(i)}(\mathbf{w},S(\mathbf{v}_{\Gamma}))^{2}}{\|\mathbf{w}\|_{1}^{2}}\leq
c\tilde{\beta}^{-2}\sup_{\mathbf{w}\in\mathbf{V}_{I}^{(i)}}\frac{b^{(i)}(\mathbf{w},S(\mathbf{v}_{\Gamma}))^{2}}{a_{h}^{(i)}(\mathbf{w},\mathbf{w})}\\\
=c\tilde{\beta}^{-2}\sup_{\mathbf{w}\in\mathbf{V}_{I}^{(i)}}\frac{a_{h}^{(i)}(T(\mathbf{v}_{\Gamma}),\mathbf{w})^{2}}{a_{h}^{(i)}(\mathbf{w},\mathbf{w})}\leq
c\tilde{\beta}^{-2}a_{h}^{(i)}(T(\mathbf{v}_{\Gamma}),T(\mathbf{v}_{\Gamma})),\end{split}$
(136)

Applying Cauchy-Schwarz to the first term in (135) and using (136), we have:

$\displaystyle\begin{split}|\mathbf{v}^{(i)}_{\Gamma}|^{2}_{S_{\Gamma}^{(i)}}=a_{h}^{(i)}(T(\mathbf{v}_{\Gamma}),T(\mathbf{v}_{\Gamma}))\leq
a_{h}^{(i)}(T(\mathbf{v}_{\Gamma}),T(\mathbf{v}_{\Gamma}))^{1/2}a_{h}^{(i)}(\mathbf{v}_{\Gamma}^{\mathcal{H}},\mathbf{v}_{\Gamma}^{\mathcal{H}})^{1/2}\\\
+c|\mathbf{v}_{\Gamma}^{\mathcal{H}}|_{H^{1}(\Omega_{i})}||S(\mathbf{v}_{\Gamma})||_{Q_{i}}\leq
a_{h}(T(\mathbf{v}_{\Gamma}),T(\mathbf{v}_{\Gamma}))^{1/2}a_{h}(\mathbf{v}_{\Gamma}^{\mathcal{H}},\mathbf{v}_{\Gamma}^{\mathcal{H}})^{1/2}+\\\
c\tilde{\beta}^{-1}a_{h}(\mathbf{v}_{\Gamma}^{\mathcal{H}},\mathbf{v}_{\Gamma}^{\mathcal{H}})^{1/2}a_{h}^{(i)}(T(\mathbf{v}_{\Gamma}),T(\mathbf{v}_{\Gamma}))^{1/2}\end{split}$

and then:

$c\tilde{\beta}^{2}|\mathbf{v}^{(i)}_{\Gamma}|^{2}_{S_{\Gamma}^{(i)}}=c\tilde{\beta}^{2}a_{h}^{(i)}(T(\mathbf{v}_{\Gamma}),T(\mathbf{v}_{\Gamma}))\leq
a_{h}^{(i)}(\mathbf{v}_{\Gamma}^{\mathcal{H}},\mathbf{v}_{\Gamma}^{\mathcal{H}}).$
(137)

Finally, from (134) and summing on the subdomains

$\displaystyle
c\tilde{\beta}^{2}|\mathbf{v}_{\Gamma}|_{S_{\Gamma}}^{2}=c\tilde{\beta}^{2}a_{h}(T(\mathbf{v}_{\Gamma}),T(\mathbf{v}_{\Gamma}))\leq
a_{h}(\mathbf{v}_{\Gamma}^{\mathcal{H}},\mathbf{v}_{\Gamma}^{\mathcal{H}})\leq
c_{3}|\mathbf{v}_{\Gamma}|_{1/2,\Gamma}^{2},$

which yields the first inequality of the thesis with $c_{1}=c/c_{3}$.

For the second inequality we have, by definition of the discrete harmonic
extension and again the stability of the discrete bilinear form $a_{h}$:

$|\mathbf{v}_{\Gamma}|_{1/2,\Gamma_{i}}^{2}\leq
a^{(i)}(\mathbf{v}_{\Gamma}^{\mathcal{H}},\mathbf{v}_{\Gamma}^{\mathcal{H}})\leq
a^{(i)}(T(\mathbf{v}_{\Gamma}),T(\mathbf{v}_{\Gamma}))\leq
c_{2}a_{h}^{(i)}(T(\mathbf{v}_{\Gamma}),T(\mathbf{v}_{\Gamma}))$ (138)

and then:

$|\mathbf{v}_{\Gamma}|_{1/2,\Gamma}^{2}\leq
c_{2}|\mathbf{v}_{\Gamma}|_{S_{\Gamma}}^{2}\text{.}$ (139)

∎

The operators $\widehat{S}_{\Gamma}$ and $\widetilde{S}_{\Gamma}$, given in
(49) and (69), are both symmetric and positive definite, because of the
Dirichlet boundary conditions on $\partial\Omega$ and provided that
sufficiently many primal constraints are chosen. We can then define the
$\widehat{S}_{\Gamma}$ and $\widetilde{S}_{\Gamma}$ norms on the spaces
$\widehat{V}_{\Gamma}$ and $\widetilde{V}_{\Gamma}$ by

$\displaystyle\begin{split}\|\mathbf{v}_{\Gamma}\|_{\widehat{S}_{\Gamma}}^{2}=\mathbf{v}_{\Gamma}^{T}R_{\Gamma}^{T}S_{\Gamma}R_{\Gamma}\mathbf{v}_{\Gamma}=|R_{\Gamma}\mathbf{v}_{\Gamma}|_{S_{\Gamma}}^{2}\quad\forall\mathbf{v}_{\Gamma}\in\widehat{\mathbf{V}}_{\Gamma}\text{,}\\\
\|\mathbf{v}_{\Gamma}\|_{\widetilde{S}_{\Gamma}}^{2}=\mathbf{v}_{\Gamma}^{T}\overline{R}_{\Gamma}^{T}S_{\Gamma}\overline{R}_{\Gamma}\mathbf{v}_{\Gamma}=|\overline{R}_{\Gamma}\mathbf{v}_{\Gamma}|_{S_{\Gamma}}^{2}\quad\forall\mathbf{v}_{\Gamma}\in\widetilde{\mathbf{V}}_{\Gamma}\text{.}\end{split}$

We then define two spaces, whose utility is that, restricted to such spaces,
the interface problem operators $\widehat{S}$ of (38) and $\widetilde{S}$ of
(89) are positive semi-definite. As in [21], we give the following:

######  Definition 5.1.

Given the discrete spaces $\widehat{\mathbf{V}}_{\Gamma}$ and
$\widetilde{\mathbf{V}}_{\Gamma}$, we define the two subspaces

$\displaystyle\begin{split}\widehat{\mathbf{V}}_{\Gamma,B}=\\{\mathbf{v}_{\Gamma}\in\widehat{\mathbf{V}}_{\Gamma}|\widehat{B}_{0\Gamma}\mathbf{v}_{\Gamma}=0\\},\\\
\widetilde{\mathbf{V}}_{\Gamma,B}=\\{\mathbf{v}_{\Gamma}\in\widetilde{\mathbf{V}}_{\Gamma}|\widetilde{B}_{0\Gamma}\mathbf{v}_{\Gamma}=0\\}.\end{split}$

We call $\widehat{\mathbf{V}}_{\Gamma,B}\times Q_{0}$ and
$\widetilde{\mathbf{V}}_{\Gamma,B}\times Q_{0}$ the benign subspaces of
$\widehat{\mathbf{V}}_{\Gamma}\times Q_{0}$ and
$\widetilde{\mathbf{V}}_{\Gamma}\times Q_{0}$.

######  Lemma 5.3.

The interface operator $\widehat{S}$ of (38), restricted to the subspace
$\widehat{\mathbf{V}}_{\Gamma,B}\times Q_{0}$ is positive semi-definite. The
same is true for $\widetilde{S}$ of (69) restricted to
$\widetilde{\mathbf{V}}_{\Gamma,B}\times Q_{0}$.

We define the $\widehat{S}$ and $\widetilde{S}$ seminorms on the benign
subspaces

$\displaystyle\begin{split}|\mathbf{v}|_{\widehat{S}}^{2}=\mathbf{v}^{T}\widehat{S}\mathbf{v}=\|\mathbf{v}_{\Gamma}\|_{\widehat{S}_{\Gamma}}\quad\forall\mathbf{v}=(\mathbf{v}_{\Gamma},q_{0})\in\widehat{\mathbf{V}}_{\Gamma,B}\times
Q_{0}\text{,}\\\
|\mathbf{v}|_{\widetilde{S}}^{2}=\mathbf{v}^{T}\widetilde{S}\mathbf{v}=\|\mathbf{v}_{\Gamma}\|_{\widetilde{S}_{\Gamma}}\quad\forall\mathbf{v}=(\mathbf{v}_{\Gamma},q_{0})\in\widetilde{\mathbf{V}}_{\Gamma,B}\times
Q_{0}\text{.}\end{split}$

Now we define an average operator $E_{D}=\widetilde{R}\widetilde{R}_{D}^{T}$,
which maps $\mathbf{\widetilde{V}}_{\Gamma}\times Q_{0}$, with generally
discontinuous interface velocities, to elements with continuous interface
velocities in the same space. For any
$\mathbf{v}=(\mathbf{v}_{\Gamma},q_{0})\in\widetilde{\mathbf{V}}_{\Gamma,B}\times
Q_{0}$,

$\displaystyle E_{D}=\left[\begin{array}[]{c}\mathbf{v}_{\Gamma}\\\
q_{0}\end{array}\right]=\left[\begin{array}[]{cc}\widetilde{R}_{\Gamma}&0\\\
0&I\end{array}\right]\left[\begin{array}[]{cc}\widetilde{R}_{D,\Gamma}&0\\\
0&I\end{array}\right]\left[\begin{array}[]{c}\mathbf{v}_{\Gamma}\\\
q_{0}\end{array}\right]=\left[\begin{array}[]{c}E_{D,\Gamma}\mathbf{v}_{\Gamma}\\\
q_{0}\end{array}\right]$ (150)

where $E_{D}=\widetilde{R}\widetilde{R}_{D,\Gamma}^{T}$, provides the average
of the interface velocities across the interface $\Gamma$. Recalling that we
can split $\mathbf{v}=\mathbf{v}_{\Pi}\oplus\mathbf{v}_{\Delta}$, we have
$E_{D}\mathbf{v}=\mathbf{v}_{\Pi}\oplus E_{D,\Delta}\mathbf{v}_{\Delta}$,
where $E_{D,\Delta}\mathbf{v}_{\Delta}$ is the dual part of the averaged
vector. As in the FEM case (see [21]) we need two assumptions to proceed in
the discussion, these will be satisfied when a reasonable choice of the primal
constraints will be done.

######  Assumption 1.

For any $\mathbf{v}_{\Delta}\in\mathbf{V}_{\Delta}$,
$\int_{\partial\Omega_{i}}\mathbf{v}_{\Delta}^{(i)}\cdot\mathbf{n}=0$ and
$\int_{\partial\Omega_{i}}(E_{D,\Delta}\mathbf{v}_{\Delta})^{(i)}\cdot\mathbf{n}=0$,
where $\mathbf{n}$ is the outward normal of $\partial\Omega_{i}$. We can
equivalently write $B_{0\Delta}^{(i)}\mathbf{v}_{\Delta}^{(i)}=0$ and
$B_{0\Delta}^{(i)}(E_{D,\Delta}\mathbf{v}_{\Delta})^{(i)}=0$

######  Assumption 2.

There exists a positive constant C, which is independent of $H$, $h$ and the
number of subdomains, such that

$\displaystyle|\bar{R}_{\Gamma}(E_{D,\Gamma}\mathbf{v}_{\Gamma})|_{1/2,\Gamma}\leq
C\bigg{(}1+\log\left(\frac{Hk^{2}}{h}\right)\bigg{)}|\bar{R}_{\Gamma}\mathbf{v}_{\Gamma}|_{1/2,\Gamma}\text{,}\quad\forall\mathbf{v}_{\Gamma}\in\mathbf{V}_{\Gamma}.$

With these two assumptions, we have the following results (proof of 5.4 in
[21]):

######  Lemma 5.4.

Let Assumption 1 hold. Then
$\widetilde{R}_{D}^{T}\mathbf{v}\in\widehat{\mathbf{V}}_{\Gamma,B}\times
Q_{0}$, for any $\mathbf{v}\in\widetilde{\mathbf{V}}_{\Gamma,B}\times Q_{0}$.

######  Lemma 5.5.

Let Assumptions 1 and 2 hold. Then there exists a positive constant C, which
is independent of $H$, $h$ and the number of subdomains, such that

$\displaystyle|E_{D}\mathbf{v}|_{\widetilde{S}}\leq
C\frac{1}{\tilde{\beta}}\bigg{(}1+\log\left(\frac{Hk^{2}}{h}\right)\bigg{)}|\mathbf{v}|_{\widetilde{S}},\quad\forall\mathbf{v}=(\mathbf{v}_{\Gamma},q_{0})\in\widetilde{\mathbf{V}}_{\Gamma,B}\times
Q_{0}\text{,}$

where $\tilde{\beta}$ is the inf-sup stability constant of (14).

###### Proof.

Given any
$\mathbf{v}=(\mathbf{v}_{\Gamma},q_{0})\in\widetilde{\mathbf{V}}_{\Gamma,B}\times
Q_{0}$, we know, from Lemma 5.4, that
$\widetilde{R}_{D}^{T}\mathbf{v}\in\widehat{\mathbf{V}}_{\Gamma,B}\times
Q_{0}$. Therefore,
$E_{D}\mathbf{v}=\widetilde{R}\widetilde{R}_{D}^{T}\mathbf{v}\in\widetilde{\mathbf{V}}_{\Gamma,B}\times
Q_{0}$. We have from the definition of the $\widetilde{S}$-seminorm, that

$|E_{D}\mathbf{v}|^{2}_{\widetilde{S}}=\|E_{D,\Gamma}\mathbf{v}_{\Gamma}\|^{2}_{\widetilde{S}_{\Gamma}}=|\bar{R}_{\Gamma}(E_{D,\Gamma}\mathbf{v}_{\Gamma})|^{2}_{S_{\Gamma}}\leq
C\frac{1}{\tilde{\beta}^{2}}|\bar{R}_{\Gamma}(E_{D,\Gamma}\mathbf{v}_{\Gamma})|^{2}_{1/2,\Gamma}\text{,}$
(151)

where the last inequality follows from Lemma 5.2. We have, from Assumption 2
and Lemma 5.2

$\begin{split}|\bar{R}_{\Gamma}(E_{D,\Gamma}\mathbf{v}_{\Gamma})|^{2}_{1/2,\Gamma}\leq
C\bigg{(}1+\log\left(\frac{Hk^{2}}{h}\right)\bigg{)}^{2}|\bar{R}_{\Gamma}\mathbf{v}_{\Gamma}|^{2}_{1/2,\Gamma}\\\
\leq
C\bigg{(}1+\log\left(\frac{Hk^{2}}{h}\right)\bigg{)}^{2}|\bar{R}_{\Gamma}\mathbf{v}_{\Gamma}|^{2}_{S_{\Gamma}}\leq
C\bigg{(}1+\log\left(\frac{Hk^{2}}{h}\right)\bigg{)}^{2}\|\mathbf{v}_{\Gamma}\|^{2}_{\widetilde{S}_{\Gamma}}\textit{.}\end{split}$
(152)

Consequently we have

$|E_{D}\mathbf{v}|^{2}_{\widetilde{S}}\leq
C\frac{1}{\tilde{\beta}^{2}}\bigg{(}1+\log\left(\frac{Hk^{2}}{h}\right)\bigg{)}^{2}\|\mathbf{v}_{\Gamma}\|^{2}_{\widetilde{S}_{\Gamma}}=C\frac{1}{\tilde{\beta}^{2}}\bigg{(}1+\log\left(\frac{Hk^{2}}{h}\right)\bigg{)}^{2}|\mathbf{v}|^{2}_{\widetilde{S}}\textit{.}$
(153)

∎

We have the following lemma (proof in [21]):

######  Lemma 5.6.

Any vector of the form
$\mathbf{u}=(\mathbf{0},p_{0})\in\widehat{\mathbf{V}}_{\Gamma,B}\times Q_{0}$
is an eigenvector of the preconditioner operator $M^{-1}\widehat{S}$ with
eigenvalue equal to 1.

######  Theorem 5.1.

Let Assumptions 1 and 2 hold. The preconditioned operator $M^{-1}\widehat{S}$
is then symmetric, positive definite with respect to the bilinear form
$\langle\cdot,\cdot\rangle_{\widehat{S}}$ on the benign space
$\widehat{\mathbf{V}}_{\Gamma,B}\times Q_{0}$. Its minimum eigenvalue is 1 and
its maximum eigenvalue is bounded by

$\displaystyle
C\frac{1}{\tilde{\beta}^{2}}\bigg{(}1+\log\left(\frac{Hk^{2}}{h}\right)\bigg{)}^{2}\text{.}$
(154)

Here, $C$ is a constant which is independent of $H$, $h$ and the number of
subdomains, and $\tilde{\beta}$ is the inf-sup stability constant defined in
(14).

###### Proof.

We know from Lemma 5.6, that any vector of the form
$\mathbf{u}=(\mathbf{0},p_{0})\in\widehat{\mathbf{V}}_{\Gamma,B}\times Q_{0}$
is an eigenvector of the preconditioned operator $M^{-1}\widehat{S}$ with an
eigenvalue equal to 1. It is sufficient to find lower and upper bounds of the
quotient
$\big{\langle}M^{-1}\widehat{S}\mathbf{u},\mathbf{u}\big{\rangle}_{\widehat{S}}/\langle\mathbf{u},\mathbf{u}\rangle_{\widehat{S}}$,
for any
$\mathbf{u}=(\mathbf{u}_{\Gamma},p_{0})\in\widehat{\mathbf{V}}_{\Gamma,B}\times
Q_{0}$, where $\mathbf{u}_{\Gamma}$ is non zero and therefore
$\langle\mathbf{u},\mathbf{u}\rangle_{\widehat{S}}>0$.
Lower bound: Given $\mathbf{u}\in\widehat{\mathbf{V}}_{\Gamma,B}\times Q_{0}$,
let

$\mathbf{v}=\widetilde{S}^{-1}\widetilde{R}_{D}\widetilde{S}\mathbf{u}\in\widetilde{\mathbf{V}}_{\Gamma,B}\times
Q_{0}\text{.}$ (155)

We have from the fact that
$\widetilde{R}^{T}\widetilde{R}_{D}=\widetilde{R}^{T}_{D}\widetilde{R}=I$,

$\langle\mathbf{u},\mathbf{u}\rangle_{\widehat{S}}=\mathbf{u}^{T}\widehat{S}\widetilde{R}^{T}_{D}\widetilde{R}\mathbf{u}=\mathbf{u}^{T}\widehat{S}\widetilde{R}^{T}_{D}\widetilde{S}^{-1}\widetilde{S}\widetilde{R}\mathbf{u}=\langle\mathbf{v},\widetilde{R}\mathbf{u}\rangle_{\widetilde{S}}.$
(156)

From the Cauchy-Schwartz inequality and the fact that
$\widehat{S}=\widetilde{R}^{T}\widetilde{S}\widetilde{R}$, we find that

$\langle\mathbf{v},\widetilde{R}\mathbf{u}\rangle_{\widetilde{S}}\leq\langle\mathbf{v},\mathbf{v}\rangle^{1/2}_{\widetilde{S}}\langle\widetilde{R}\mathbf{u},\widetilde{R}\mathbf{u}\rangle^{1/2}_{\widetilde{S}}=\langle\mathbf{v},\mathbf{v}\rangle^{1/2}_{\widetilde{S}}\langle\mathbf{u},\mathbf{u}\rangle^{1/2}_{\widetilde{S}}\text{.}$
(157)

Therefore from (156) and (157),

$\langle\mathbf{u},\mathbf{u}\rangle_{\widetilde{S}}\leq\langle\mathbf{v},\mathbf{v}\rangle_{\widetilde{S}}\text{.}$
(158)

Since,

$\langle\mathbf{v},\mathbf{v}\rangle_{\widetilde{S}}=\mathbf{u}^{T}\widehat{S}\widetilde{R}^{T}_{D}\widetilde{S}^{-1}\widetilde{S}\widetilde{S}^{-1}\widetilde{R}_{D}\widehat{S}\mathbf{u}=\big{\langle}\mathbf{u},\widetilde{R}^{T}_{D}\widetilde{S}^{-1}\widetilde{R}_{D}\widehat{S}\mathbf{u}\rangle_{\widehat{S}}=\big{\langle}\mathbf{u},M^{-1}\widehat{S}\mathbf{u}\big{\rangle}_{\widehat{S}}\text{,}$
(159)

we obtain, from equations (158) and (159), that
$\langle\mathbf{u},\mathbf{u}\rangle_{\widetilde{S}}\leq\big{\langle}\mathbf{u},M^{-1}\widehat{S}\mathbf{u}\big{\rangle}_{\widehat{S}}$,
which gives a lower bound of 1 for the eigenvalues. Then from Lemma 5.6, we
know that 1 is the minimum eigenvalue of the preconditioned operator.
Upper bound: Given $\mathbf{u}\in\widehat{\mathbf{V}}_{\Gamma,B}\times Q_{0}$,
take $\mathbf{v}\in\widetilde{\mathbf{V}}_{\Gamma,B}\times Q_{0}$ as in (155).
We have, $\widetilde{R}^{T}_{D}\mathbf{v}=M^{-1}\widehat{S}\mathbf{u}$. Since
$\widehat{S}=\widetilde{R}^{T}\widetilde{S}\widetilde{R}$ and by using Lemma
5.5, we have

$\begin{split}\langle
M^{-1}\widehat{S}\mathbf{u},M^{-1}\widehat{S}\mathbf{u}\rangle_{\widehat{S}}=\langle\widetilde{R}^{T}_{D}\mathbf{v},\widetilde{R}^{T}_{D}\mathbf{v}\rangle_{\widehat{S}}=\langle\widetilde{R}\widetilde{R}^{T}_{D}\mathbf{v},\widetilde{R}\widetilde{R}^{T}_{D}\mathbf{v}\rangle_{\widetilde{S}}\\\
=|E_{D}\mathbf{v}|^{2}_{\widetilde{S}}\leq
C^{2}\frac{1}{\tilde{\beta}^{2}}\bigg{(}1+\log\left(\frac{Hk^{2}}{h}\right)\bigg{)}^{2}|\mathbf{v}|^{2}_{\widetilde{S}}\text{.}\end{split}$
(160)

Therefore from equation (159), we have

$\langle
M^{-1}\widehat{S}\mathbf{u},M^{-1}\widehat{S}\mathbf{u}\rangle_{\widehat{S}}\leq
C^{2}\frac{1}{\tilde{\beta}^{2}}\bigg{(}1+\log\left(\frac{Hk^{2}}{h}\right)\bigg{)}^{2}\langle\mathbf{u},M^{-1}\widehat{S}\mathbf{u}\rangle_{\widehat{S}}\text{.}$
(161)

Using the Cauchy-Schwarz inequality and equation (161), we have

$\begin{split}\langle\mathbf{u},M^{-1}\widehat{S}\mathbf{u}\rangle_{\widehat{S}}\leq\langle
M^{-1}\widehat{S}\mathbf{u},M^{-1}\widehat{S}\mathbf{u}\rangle^{1/2}_{\widehat{S}}\langle\mathbf{u},\mathbf{u}\rangle^{1/2}_{\widehat{S}}\\\
\leq
C\frac{1}{\tilde{\beta}}\bigg{(}1+\log\left(\frac{Hk^{2}}{h}\right)\bigg{)}\langle\mathbf{u},\mathbf{u}\rangle^{1/2}_{\widehat{S}}\langle\mathbf{u},M^{-1}\widehat{S}\mathbf{u}\rangle^{1/2}_{\widehat{S}}\text{.}\end{split}$
(162)

This gives

$\langle\mathbf{u},M^{-1}\widehat{S}\mathbf{u}\rangle_{\widehat{S}}\leq
C\frac{1}{\tilde{\beta}^{2}}\bigg{(}1+\log\left(\frac{Hk^{2}}{h}\right)\bigg{)}^{2}\langle\mathbf{u},\mathbf{u}\rangle_{\widehat{S}}\text{,}$
(163)

and the upper bound of the theorem. ∎

##  6 Satisfying the Assumptions

To satisfy the assumptions 1 and 2 in of the previous section we have to
choose properly the primal constraints for the interface velocity space. In
particular to satisfy the Assumption 1, it is not sufficient to choose as
primal constraints the subdomain vertices, i.e. make both components of the
velocity continuous at those nodes, but some extra edge constraints are
necessary. To do so, for each interface edge $\Gamma_{ij}$, which is shared by
a pair of subdomains $\Omega_{i}$ and $\Omega_{j}$, we impose

$\displaystyle\int_{\Gamma_{i}j}\mathbf{v}_{\Gamma}^{(i)}\cdot\mathbf{n}_{ij}=\int_{\Gamma_{i}j}\mathbf{v}_{\Gamma}^{(j)}\cdot\mathbf{n}_{ij}$
(164)

for a fixed selection of the normal $\mathbf{n}_{ij}$ of $\Gamma_{i,h}$.
Proceeding the discussion with this first choice, after changing the
variables, the dual interface velocity component will vanish at the subdomain
vertices and its normal component will have a weighted zero average over each
$\Gamma_{ij}$, i.e.

$\displaystyle\int_{\Gamma_{i}}\mathbf{v_{\Gamma}}^{(i)}\cdot\mathbf{n}_{ij}=\int_{\Gamma_{j}}\mathbf{v_{\Gamma}}^{(j)}\cdot\mathbf{n}_{ij}=0.$

By the definition of the average operator $E_{D,\Delta}$ we have that the
average interface velocity is
$E_{D,\Delta}\mathbf{v}_{\Delta}=\frac{1}{2}(\mathbf{v}_{\Delta}^{(i)}+\mathbf{v}_{\Delta}^{(j)})$
on each edge and hence

$\displaystyle\int_{\Gamma_{i}j}(E_{D,\Delta}\mathbf{v}_{\Delta})^{(i)}\cdot\mathbf{n}_{ij}=0\text{.}$
(165)

In our codes we also choose a strong condition, we decide to require that the
integral of both velocity components have common values across each interface
edge

$\displaystyle\int_{\Gamma_{i}j}\mathbf{v}_{\Gamma}^{(i)}\big{|}_{x}=\int_{\Gamma_{j}}\mathbf{v}_{\Gamma}^{(j)}\big{|}_{x}\text{,}\quad\int_{\Gamma_{i}j}\mathbf{v}_{\Gamma}^{(i)}\big{|}_{y}=\int_{\Gamma_{j}}\mathbf{v}_{\Gamma}^{(j)}\big{|}_{y}\text{.}$
(166)

In this way, we clearly satisfy Assumption 1. The advantage of this condition
is that it is easiest to implement and, enlarging a little the coarse space,
it yields a faster convergence. Assumption 2 is also satisfied, requiring only
vertices as primal constraints, and it derives directly from the following
lemma, proved in [4]:

######  Lemma 6.1.

For all $\mathbf{v}_{\Gamma}\in\widetilde{\mathbf{V}}_{\Gamma}$ we have:

$\displaystyle|E_{D,\Gamma}\mathbf{v}_{\Gamma}|_{1/2,\Gamma}\leq
C\left(1+\log\left(\frac{Hk^{2}}{h}\right)\right)|\mathbf{v}_{\Gamma}|_{1/2,\Gamma}$
(167)

with $C$ positive constant independent of $H$, $h$ and the number of
subdomains.

##  7 Numerical Results

In this section, we provide some numerical tests to study the behavior of the
BDDC preconditioner with respect to the mesh size $h$, the number of
subdomains $N$ and the shape of the polygonal mesh elements. We solve the
Stokes equations on the unit square domain $\Omega=[0,1]\times[0,1]$, applying
homogeneous Dirichlet boundary conditions on the whole $\partial\Omega$. We
choose the load term $\mathbf{f}$ by imposing that the analytical solution is

$\displaystyle\mathbf{u}(x,y)=\left(\begin{array}[]{c}-\sin(\pi x)\sin(\pi
x)\sin(2\pi y)\\\ \sin(\pi y)\sin(\pi y)\sin(2\pi x)\end{array}\right),\quad
p(x,y)=\sin(\pi x)-\sin(\pi y)\text{.}$

(a) QUAD mesh. (b) HEXA mesh.

(c) TRI mesh. (d) CVT mesh.

Figure 2: Examples of the different type of meshes we consider in our
numerical experiments. Table 1: GMRES iteration counts to solve the interface
saddle-point problem without preconditioner varying the mesh size $h$ and the
number of subdomains $N=1/H^{2}$.

1/H 1/h |  | 8
---
it
| 16
---
it
| 32
---
it
| 64
---
it
| 128
---
it
2 | 7 | 13 | 18 | 27 | 40
4 | x | 40 | 57 | 79 | 99
8 | x | x | 140 | 203 | 282
16 | x | x | x | 295 | 437
32 | x | x | x | x | 605

(a) QUAD meshes.

1/H 1/h |  | 8
---
it
| 16
---
it
| 32
---
it
| 64
---
it
| 128
---
it
2 | 18 | 23 | 31 | 42 | 56
4 | x | 91 | 131 | 163 | 199
8 | x | x | 221 | 295 | 434
16 | x | x | x | 486 | 583
16 | x | x | x | x | 976

(b) HEXA meshes.

1/H 1/h |  | 8
---
it
| 16
---
it
| 32
---
it
| 64
---
it
| 128
---
it
2 | 12 | 17 | 25 | 38 | 52
4 | x | 42 | 63 | 82 | 105
8 | x | x | 127 | 190 | 286
16 | x | x | x | 380 | 569
32 | x | x | x | x | 1130

(c) TRI meshes.

1/H 1/h |  | 8
---
it
| 16
---
it
| 32
---
it
| 64
---
it
| 128
---
it
2 | 31 | 44 | 55 | 79 | 101
4 | x | 89 | 104 | 134 | 172
8 | x | x | 267 | 307 | 336
16 | x | x | x | 885 | 924
32 | x | x | x | x | 2522

(d) CVT meshes.

Table 2: PCG iteration counts to solve the interface saddle-point problem with
BDDC preconditioner, varying the mesh size $h$ and the number of subdomains
$N=1/H^{2}$. The primal coarse space is spanned only by the subdomain
vertices.

1/H 1/h |  | 8
---
it
| 16
---
it
| 32
---
it
| 64
---
it
| 128
---
it
2 | 7 | 7 | 7 | 7 | 8
4 | x | 12 | 13 | 15 | 17
8 | x | x | 22 | 23 | 29
16 | x | x | x | 22 | 26
32 | x | x | x | x | 22

(a) QUAD meshes.

1/H 1/h |  | 8
---
it
| 16
---
it
| 32
---
it
| 64
---
it
| 128
---
it
2 | 8 | 9 | 9 | 9 | 9
4 | x | 18 | 19 | 21 | 22
8 | x | x | 26 | 30 | 33
16 | x | x | x | 28 | 30
32 | x | x | x | x | 29

(b) HEXA meshes.

1/H 1/h |  | 8
---
it
| 16
---
it
| 32
---
it
| 64
---
it
| 128
---
it
2 | 8 | 8 | 8 | 9 | 9
4 | x | 15 | 17 | 19 | 20
8 | x | x | 18 | 24 | 27
16 | x | x | x | 20 | 25
32 | x | x | x | x | 20

(c) TRI meshes.

1/H 1/h |  | 8
---
it
| 16
---
it
| 32
---
it
| 64
---
it
| 128
---
it
2 | 16 | 17 | 17 | 17 | 17
4 | x | 26 | 30 | 32 | 33
8 | x | x | 36 | 41 | 41
16 | x | x | x | 50 | 50
32 | x | x | x | x | 51

(d) CVT meshes.

Table 3: PCG iteration counts to solve the interface problem with BDDC
preconditioner, varying the mesh size $h$ and the number of subdomains
$N=1/H^{2}$. The primal coarse space is spanned by the subdomain vertices and
only one basis function per subdomain edge.

1/H 1/h |  | $\kappa_{2}$
---
| 8
---
it
| $\kappa_{2}$
---
| 16
---
it
| $\kappa_{2}$
---
| 32
---
it
| $\kappa_{2}$
---
| 64
---
it
| $\kappa_{2}$
---
| 128
---
it
2 | 1,82 | 7 | 2,11 | 7 | 2,37 | 7 | 2,80 | 7 | 3,24 | 8
4 | x | x | 4,40 | 9 | 5,75 | 10 | 7,20 | 11 | 8,76 | 12
8 | x | x | x | x | 5,78 | 13 | 7,81 | 15 | 10,09 | 16
16 | x | x | x | x | x | x | 6,16 | 16 | 8,46 | 19
32 | x | x | x | x | x | x | x | x | 6,29 | 16

(a) QUAD meshes.

1/H 1/h |  | $\kappa_{2}$
---
| 8
---
it
| $\kappa_{2}$
---
| 16
---
it
| $\kappa_{2}$
---
| 32
---
it
| $\kappa_{2}$
---
| 64
---
it
| $\kappa_{2}$
---
| 128
---
it
2 | 3,35 | 9 | 4,45 | 9 | 5,49 | 9 | 6,64 | 9 | 6,95 | 9
4 | x | x | 5,32 | 13 | 6,86 | 14 | 8,45 | 15 | 10,17 | 16
8 | x | x | x | x | 7,34 | 17 | 10,12 | 19 | 12,37 | 20
16 | x | x | x | x | x | x | 7,97 | 18 | 11,07 | 22
32 | x | x | x | x | x | x | x | x | 8,20 | 18

(b) HEXA meshes.

1/H 1/h |  | $\kappa_{2}$
---
| 8
---
it
| $\kappa_{2}$
---
| 16
---
it
| $\kappa_{2}$
---
| 32
---
it
| $\kappa_{2}$
---
| 64
---
it
| $\kappa_{2}$
---
| 128
---
it
2 | 2,73 | 8 | 3,51 | 8 | 4,28 | 8 | 5,12 | 9 | 6,32 | 9
4 | x | x | 4,01 | 11 | 5,20 | 12 | 6,54 | 13 | 7,98 | 14
8 | x | x | x | x | 5,01 | 15 | 6,83 | 16 | 8,93 | 18
16 | x | x | x | x | x | x | 5,25 | 15 | 7,28 | 17
32 | x | x | x | x | x | x | x | x | 5,32 | 15

(c) TRI meshes.

1/H 1/h |  | $\kappa_{2}$
---
| 8
---
it
| $\kappa_{2}$
---
| 16
---
it
| $\kappa_{2}$
---
| 32
---
it
| $\kappa_{2}$
---
| 64
---
it
| $\kappa_{2}$
---
| 128
---
it
2 | 5,82 | 14 | 6,97 | 15 | 8,16 | 16 | 9,32 | 16 | 10,23 | 16
4 | x | x | 10,20 | 20 | 13,87 | 21 | 15,98 | 22 | 17,22 | 23
8 | x | x | x | x | 22,24 | 27 | 21,43 | 28 | 23,13 | 28
16 | x | x | x | x | x | x | 30,12 | 34 | 28,34 | 33
32 | x | x | x | x | x | x | x | x | 30,28 | 35

(d) CVT meshes.

Table 4: PCG iteration counts to solve the interface saddle-point problem with
BDDC preconditioner, varying the mesh size $h$ and the number of subdomains
$N=1/H^{2}$. The primal coarse space is spanned by the subdomain vertices and
two basis functions per subdomain edge.

1/H 1/h |  | $\kappa_{2}$
---
| 8
---
it
| $\kappa_{2}$
---
| 16
---
it
| $\kappa_{2}$
---
| 32
---
it
| $\kappa_{2}$
---
| 64
---
it
| $\kappa_{2}$
---
| 128
---
it
2 | 1,48 | 7 | 1,68 | 7 | 1,90 | 7 | 2,18 | 8 | 2,51 | 7
4 | x | x | 2,80 | 9 | 3,73 | 10 | 4,78 | 10 | 5,93 | 11
8 | x | x | x | x | 2,99 | 10 | 4,05 | 11 | 5,20 | 13
16 | x | x | x | x | x | x | 2,72 | 9 | 3,67 | 10
32 | x | x | x | x | x | x | x | x | 2,64 | 8

(a) QUAD meshes.

1/H 1/h |  | $\kappa_{2}$
---
| 8
---
it
| $\kappa_{2}$
---
| 16
---
it
| $\kappa_{2}$
---
| 32
---
it
| $\kappa_{2}$
---
| 64
---
it
| $\kappa_{2}$
---
| 128
---
it
2 | 3,33 | 9 | 4,29 | 10 | 5,36 | 10 | 6,68 | 10 | 8.08 | 11
4 | x | x | 4,21 | 12 | 5,29 | 13 | 6,58 | 15 | 7,90 | 15
8 | x | x | x | x | 4,59 | 13 | 5,90 | 14 | 6,65 | 15
16 | x | x | x | x | x | x | 4,79 | 14 | 5,12 | 13
32 | x | x | x | x | x | x | x | x | 4,35 | 12

(b) HEXA meshes.

1/H 1/h |  | $\kappa_{2}$
---
| 8
---
it
| $\kappa_{2}$
---
| 16
---
it
| $\kappa_{2}$
---
| 32
---
it
| $\kappa_{2}$
---
| 64
---
it
| $\kappa_{2}$
---
| 128
---
it
2 | 2,49 | 9 | 3,41 | 9 | 4,38 | 10 | 5,40 | 10 | 6,55 | 10
4 | x | x | 2,96 | 10 | 3,85 | 11 | 5,17 | 12 | 6,36 | 14
8 | x | x | x | x | 3,26 | 10 | 4,26 | 12 | 5,33 | 13
16 | x | x | x | x | x | x | 3,44 | 9 | 4,42 | 11
32 | x | x | x | x | x | x | x | x | 3,49 | 8

(c) TRI meshes.

1/H 1/h |  | $\kappa_{2}$
---
| 8
---
it
| $\kappa_{2}$
---
| 16
---
it
| $\kappa_{2}$
---
| 32
---
it
| $\kappa_{2}$
---
| 64
---
it
| $\kappa_{2}$
---
| 128
---
it
2 | 4,35 | 13 | 5,27 | 14 | 6,63 | 15 | 7,31 | 16 | 8,28 | 16
4 | x | x | 5,20 | 15 | 10,22 | 18 | 13,00 | 19 | 15,63 | 20
8 | x | x | x | x | 9,03 | 20 | 17,52 | 23 | 18,41 | 22
16 | x | x | x | x | x | x | 12,29 | 21 | 19,21 | 24
32 | x | x | x | x | x | x | x | x | 14,43 | 23

(d) CVT meshes.

In the following tables, we report the number of iterations to solve the
global interface saddle-point problem (38) with the non-preconditioned GMRES
method or the PCG method, accelerated by BDDC. Where possible, we estimate the
extreme eigenvalues using the Lanczos trick. Both in case of PCG and GMRES, we
set the tolerance for the relative residual error to $10^{-6}$. Note that in
the tables we marked with an "x" the numerical tests that we do not have
performed because they are not significant.
Our tests have been executed on different types of polygonal meshes and using
the VEM discretization with degree $k=2$ with the divergence free approach,
that means having polynomials of degree 2 on the boundary of each element for
the velocity and piecewise constant functions for the pressure. We underline
the fact that we would have obtained the same behavior, both in terms of
number of iterations and spectral condition number number, also in the case of
neglecting the divergence free property, because the interface problem and the
preconditioner are exactly the same due to the decomposition technique used in
(19) and (20).

The polygonal meshes considered are quadrilateral (QUAD), hexagonal (HEXA),
triangular (TRI) and Voronoi (CVT) (Figure 2).

Table 1 reports the number of iterations to solve the interface saddle-point
problem with the non-preconditioned GMRES. As expected, we observe that the
iteration counts grow when the number of subdomains increases and the mesh
size decreases.

Table 2 reports the number of iterations to solve the interface saddle-point
problem with PCG, preconditioned by BDDC, considering as primal constraints
only the subdomain vertices. In this case the solver appears to be scalable,
since, moving along the diagonals of the table, the iterations remain bounded
when the number of subdomains increase, and quasi-optimal, since, moving along
the rows of the table, the growth of iterations seems logarithmic. The results
also show that the solver suffers more on the Voronoi meshes than on the
others. We recall that with this choice of primal constraints the assumption 1
is not satisfied, therefore the preconditioned system is not positive definite
and we are not able to give an estimate on the eigenvalues.

Table 3 reports the spectral condition number of the preconditioned system and
the iteration counts to solve the interface problem with the PCG method,
preconditioned by BDDC, where the primal constraints are the subdomain
vertices and one basis function for each subdomain edge. In this case both the
assumptions are satisfied, therefore the system is symmetric and positive
definite and we are able to give an estimate of the eigenvalues. The results
confirm the theoretical estimates, since both the condition number and the
number of iterations are independent of number of subdomains (scalability) and
exhibit a logarithmic growth with respect to the ratio $H/h$ (quasi-
optimality).

Table 4 reports the spectral condition number of the preconditioned system and
the iteration counts to solve the interface problem with the PCG method,
preconditioned by BDDC, where the primal constraints are the subdomain
vertices and two basis functions for each subdomain edge. In this case the
system is again symmetric and positive definite, thus we are able to give an
estimate of the eigenvalues. Both the condition number and the iteration
counts exhibit a scalable and quasi-optimal behavior as before, but in this
case the convergence is faster since the coarse problem is slightly larger.

We recall that our code is implemented in Matlab and the tests were performed
in serial, therefore we do not provide an analysis on the time of
computations.

(a) The primal space consists of the subdomain vertices and two basis
functions per subdomain edge. The number of subdomains is fixed to $N=16$. (b)
The primal space consists of the subdomain vertices and two basis functions
per subdomain edge. The ratio $H/h$ is fixed to 4.

(c) The primal space consists of the subdomain vertices and one basis function
per subdomain edge. The number of subdomains is fixed to $N=16$. (d) The
primal space consists of the subdomain vertices and one basis function per
subdomain edge. The ratio $H/h$ is fixed to 4.

Figure 3: Plots of the spectral condition number ($\kappa_{2}$) of the BDDC
preconditioned linear system as a function of the ratio $H/h$ (left) and of
the number of subdomains $N$ (right) for different types of mesh and choices
of the primal space.

(a) QUAD meshes, 16 subdomains. (b) QUAD meshes, fixed local size (H/h=4).

(c) HEXA meshes, 16 subdomains. (d) HEXA meshes, fixed local size (H/h=4).

Figure 4: Plots of the PCG iteration counts of the BDDC preconditioner for
different choices of primal space on quadrilateral (QUAD) and hexagonal (HEXA)
meshes.

(a) TRI meshes, 16 subdomains. (b) TRI meshes, fixed local size (H/h=4).

(c) CVT meshes, 16 subdomains. (d) CVT meshes, fixed local size (H/h=4).

Figure 5: Plots of the PCG iteration counts of the BDDC preconditioner for
different choices of primal space on triangular (TRI) and Voronoi (CVT)
meshes.

In Figure 3, we plot the spectral condition number of the BDDC preconditioner
with the two different choices of primal constraints that satisfy the
assumptions. The left column displays an optimality test, fixing at 16 the
number of subdomains and increasing the ratio $H/h$. We observe the
logarithmic growth of the condition number. The right column displays a weak
scalability test, fixing the ratio $H/h=4$ and increasing the number of
subdomains. In this case we see that the condition number remains bounded when
the number of subdomains increases. We observe a worse behavior for the
Voronoi meshes due to the fact that the boundary of the subdomains are quite
irregular.

Finally, in Figures 4 and 5, we plot the PCG iteration counts of the BDDC
preconditioner for different choices of primal constraints and meshes. The
left column reports an optimality test with 16 subdomains and we observe that
the logarithmic growth is respected, with a smaller number of iterations when
the coarse space is enriched. The right column displays the number of PCG
iterations for a fixed local problem size ($H/h=4$) and we observe that the
number of iterations remains bounded when the number of subdomains increase,
again with a smaller number of iterations for richer primal spaces.

##  8 Conclusions

In this work, we have analyzed BDDC preconditioners to solve the saddle-point
linear system deriving from a divergence free VEM discretization of the steady
two-dimensional Stokes equations. The numerical tests have validated the
convergence estimates, showing the scalability and quasi-optimality of the
algorithm, under appropriate choices of the primal coarse space. We have also
obtained a better behavior and a faster convergence of the method for an
enriched primal space, easy to implement.

##  9 Acknowledgements

The Authors are grateful to INdAM-GNCS for the support, we are also grateful
to Giuseppe Vacca who provided us the initial code for the VEM discretization
of the Stokes system.

## Data Availability Statement

The numerical experiments have been performed using an in-house matlab code
available upon request to the authors.

## Conflict of interest

The authors declare that they have no conflict of interest.

## References

  * [1] P. F. Antonietti, L. Mascotto, and M. Verani.  A multigrid algorithm for the p-version of the virtual element method.  ESAIM: Math. Model. Numer. Anal., 52(1):337–364, 2018.
  * [2] L. Beirão da Veiga, F. Brezzi, A. Cangiani, G. Manzini, L. D. Marini, and A. Russo.  Basic principles of virtual element methods.  Math. Mod. Meth. Appl. Sci., 23(1):199–214, 2013.
  * [3] L. Beirão da Veiga, C. Lovadina, and G. Vacca.  Divergence free virtual elements for the stokes problem on polygonal meshes.  ESAIM: Math. Mod. Numer. Anal., 51(2):509–535, 2017.
  * [4] S. Bertoluzza, M. Pennacchio, and D. Prada.  Bddc and feti-dp for the virtual element method.  Calcolo, 54(4):1565–1593, 2017.
  * [5] S. Bertoluzza, M. Pennacchio, and D. Prada.  FETI-DP for the three dimensional virtual element method.  SIAM J. Numer. Anal., 58(3):1556–1591, 2020.
  * [6] D. Boffi, F. Brezzi, and M. Fortin.  Mixed finite element methods and applications, volume 44.  Springer, 2013.
  * [7] J. H. Bramble and J. E. Pasciak.  A domain decomposition technique for stokes problems.  Appl. Numer. Math., 6(4):251–261, 1990.
  * [8] S. C. Brenner and L.-Y. Sung.  Bddc and feti-dp without matrices or vectors.  Comput. Meth. Appl. Mech. Eng., 196(8):1429–1435, 2007.
  * [9] J. G. Calvo.  On the approximation of a virtual coarse space for domain decomposition methods in two dimensions.  Math. Mod. Meth. Appl. Sci., 28(07):1267–1289, 2018.
  * [10] J. G. Calvo.  An overlapping schwarz method for virtual element discretizations in two dimensions.  Comput. Math. Appl., 77(4):1163–1177, 2019.
  * [11] C. Canuto, L. F. Pavarino, and A. B. Pieri.  BDDC preconditioners for continuous and discontinuous galerkin methods using spectral/hp elements with variable local polynomial degree.  IMA J. Numer. Anal., 34(3):879–903, 2014.
  * [12] F. Dassi and S. Scacchi.  Parallel block preconditioners for three-dimensional virtual element discretizations of saddle-point problems.  Comput. Meth. Appl. Mech. Eng., 372:113424, 2020.
  * [13] F. Dassi, S. Zampini, and S. Scacchi.  Robust and scalable adaptive BDDC preconditioners for virtual element discretizations of elliptic partial differential equations in mixed form.  Comput. Meth. Appl. Mech. Eng., 391:114620, 2022.
  * [14] C. R. Dohrmann.  A preconditioner for substructuring based on constrained energy minimization.  SIAM J. Sci. Comput., 25(1):246–258, 2003.
  * [15] C. R. Dohrmann and O. B. Widlund.  A BDDC algorithm with deluxe scaling for three-dimensional H (curl) problems.  Comm. Pure Appl. Math., 69(4):745–770, 2016.
  * [16] M. Dryja, J. Galvis, and M. Sarkis.  BDDC methods for discontinuous galerkin discretization of elliptic problems.  J. Complex., 23(4-6):715–739, 2007.
  * [17] C. Hofer.  Analysis of discontinuous galerkin dual-primal isogeometric tearing and interconnecting methods.  Math. Mod. Meth. Appl. Sci., 28(1):131–158, 2017.
  * [18] H. H. Kim, M. Dryja, and O. B. Widlund.  A BDDC method for mortar discretizations using a transformation of basis.  SIAM J. Numer. Anal., 47(1):136–157, 2009.
  * [19] A. Klawonn and O. B. Widlund.  Dual-primal FETI methods for linear elasticity.  Comm. Pure Appl. Math., 59(11):1523–1572, 2006.
  * [20] J. Li and X. Tu.  A nonoverlapping domain decomposition method for incompressible stokes equations with continuous pressures.  SIAM J. Numer. Anal., 51(2):1235–1253, 2013.
  * [21] J. Li and O. Widlund.  Bddc algorithms for incompressible stokes equations.  SIAM J. Numer. Anal., 44(6):2432–2455, 2006.
  * [22] J. Li and O. B. Widlund.  Feti-dp, bddc, and block cholesky methods.  Int. J. Numer. Meth. Eng., 66(2):250–271, 2006.
  * [23] D.-S. Oh, O. B. Widlund, S. Zampini, and C. R. Dohrmann.  BDDC algorithms with deluxe scaling and adaptive selection of primal constraints for Raviart–Thomas vector fields.  Math. Comput., 87(310):659–692, 2017.
  * [24] A. Smith, P. Bjørstad, and W. Gropp.  Domain Decomposition: Parallel Multilevel Methods for Elliptic Partial Differential Equations.  Cambridge University Press, 2004.
  * [25] A. Toselli and O. B. Widlund.  Domain decomposition methods-algorithms and theory, volume 34.  Springer Science & Business Media, 2006.
  * [26] X. Tu and B. Wang.  A BDDC algorithm for the stokes problem with weak galerkin discretizations.  Comput. Math. Appl., 76(2):377–392, 2018.
  * [27] O. B. Widlund, S. Zampini, S. Scacchi, and L. F. Pavarino.  Block FETI–DP/BDDC preconditioners for mixed isogeometric discretizations of three-dimensional almost incompressible elasticity.  Math. Comp., 90(330):1773–1797, 2021.
  * [28] S. Zampini.  Dual-primal methods for the cardiac bidomain model.  Mathematical Models and Methods in Applied Sciences, 24(04):667–696, 2014.