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

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.


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 nonpolynomial 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.

Continuous problem
Let Ω ⊆ R 2 , with Γ = ∂Ω, and consider the stationary Stokes problem on Ω with homogeneous Dirichet boundary conditions: where u and p are the velocity and the pressure fields, respectively.Furthermore ∆, div and ∇ denote the vector Laplacian, the divergence and the gradient operators.Finally, f represents the external force, while ν > 0 is the viscosity.
Let us consider the spaces: with norms: We assume f ∈ [H −1 (Ω)] 2 , and ν ∈ L ∞ (Ω) uniformly positive in Ω.Let the bilinear forms a(•, •) : V × V → R and b : V × Q → R be defined as: b(v, q) := Ω div vq dΩ for all u ∈ V, q ∈ Q. ( Then a standard variational formulation of problem (1) reads: where It is well-known that: where a and b are the usual norm of the two bilinear forms; • a(•, •) is coercive i.e., there exists a positive constant α such that • the bilinear form b(•, •) satisfies the inf-sup condition [6], i.e.
∃β > 0 such that sup Therefore, problem (6) has a unique solution (u, p) ∈ V × Q such that where the constant C depends only on Ω and ν; see [6].

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 {T h } h be a sequence of triangulations of Ω into general polygonal elements K with We suppose that, for all h, each element K ∈ T h satisfies the following assumptions: • (A1) K is star-shaped with respect to a ball of radius ≥ γh k , • (A2) the distance between any two vertices of K is ≥ ch K , • (A3) the triangulation T h is quasi-uniform, i.e. there exist positive constants c 0 , c 1 such that for any two elements K and where γ 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 (A1).
We also assume that the scalar viscosity field ν is piecewise constant with respect to the decomposition T h , i.e. ν is constant on each polygon K ∈ T h .
For k ∈ N, let us define the spaces: On each element K ∈ T h we define, for k ≥ 2, the following finite dimensional local virtual element spaces: and Now it is possible to introduce suitable sets of degrees of freedom for the local approximations fields.Given a function v ∈ V K h we take the following linear operators D V , split into three subsets: • D V 1: the values of v at the vertices of the polygon K, • D V 2: the values of v at k − 1 distinct points of every edge e ∈ ∂K (for the implementation we will take the k − 1 internal points of the (k + 1)-Gauss-Lobatto quadrature rule in e), • D V 3: the moments of the values of v Furthermore, for the local pressure, given q ∈ Q K h , we consider the linear operators Since D V and D Q are unisolvent respectively of V K h and Q K h , we can define the global virtual element spaces: and with obvious associated sets of global degrees of freedom.

Discrete problem
Referring to [3], we can now state the discrete virtual element problem By construction the discrete bilinear form a h (•, •) is stable (uniformly) with respect to the V norm and also obviously the bilinear form b(•, •).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 Ω 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 V h and Q h defined in (11) and ( 12), there exists a positive β, independent of h, such that: A consequence of the previous proposition is the following statement.

verifying the estimate
We have also a convergence result Theorem 3.2.Let (u, p) ∈ V × Q be the solution of problem (6) and (u h , p h ) ∈ V h × Q h be the solution of problem (13).Then it holds and 4 Construction of the BDDC preconditioner In this section, we first divide the domain Ω 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.

Domain decomposition
We decompose the domain Ω into N non-overlapping subdomains Ω 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 Γ = (∪∂Ω i ) \ ∂Ω; we define also Γ i = ∂Ω i ∩ Γ as the interface of an individual subdomain Ω i .
According to [5], where more details could be found, we recall two requirements on the subdomain partition: • (S1) Each subdomain Ω i is the union of polygonal elements of the triangulation T h and the number of polygons forming an individual subdomain is uniformly bounded; • (S2) If a face of a subdomain intersects ∂Ω, then the measure of this set is comparable to that of ∂Ω i .Similarly, if an edge of a subdomain intersects ∂Ω, the length of this intersection is bounded from below in terms of the diameter of ∂Ω i .
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.

Decomposition of the virtual element spaces
The discrete variational problem (13) can be written, in matrix form, as the following saddle-point linear system: where the matrices A and B are associated with the discrete bilinear forms a h (•, •) and b(•, •).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 V × Q instead of V h × 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 (D V 1 and D V 2) and interior dofs (D V 3 and D V 4).Following the notations introduced in [21], we decompose the discrete velocity and pressure space V and Q into: V I and Q I are direct sums of subdomain interior velocity spaces V (i) I , and subdomain interior pressure spaces Q (i) I , respectively, i.e., The elements of V I have support in the subdomain Ω i and vanish on its interface Γ i , while the elements of Q (i) I are restrictions of elements in Q to Ω i .V Γ is the space of the traces on Γ of functions in V and Q 0 is the subspace of Q with constant values q (i) 0 in the subdomain Ω i .We denote the space of interface velocity variables of the subdomain Ω i by V (i) Γ , and the associated product space by Γ ; generally functions in V Γ are discontinuous across the interface.
Γ is the operator which maps functions in the continuous interface velocity space V Γ to their subdomain components in the space V (i) Γ .We denote the direct sum of the R With the decomposition of the solution space given in (19), the global saddle-point problem (18) can be written as: find Remark 4.1.Here the lower left block of ( 21) is zero because the bilinear form b(u I , q 0 ) vanishes for any v I ∈ V I and q 0 ∈ 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., Γ .Correspondingly, the right-hand side vector f I consists of subdomain vectors f (i) I , and f Γ is assembled from the subdomain components f (i) Γ ; we denote the spaces of the right-hand side vectors f I and f Γ by F I and F Γ 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: 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 (u I ) and (u I ) in the system (22).To do so, we solve two independent Dirichlet problems: Then, substituting the solutions of ( 25) and (26) in we obtain the global interface saddle-point problem: where the right-hand side g ∈ F Γ × F 0 is given by We note that S is assembled from the subdomain Stokes Schur complements S (i) , which are defined by: given Denoting by S Γ the direct sum of the S Γ , then S Γ is given by and then we set Finally we see from (30), that the action of S (i) on a vector can be evaluated by solving a Dirichlet problem on the subdomain Ω i as in ( 25) and ( 26), so it is not necessary to assemble the matrix S because only its action is required.
In the next section we introduce a BDDC preconditioner for problem (28), where the operator of the preconditioned problem is symmetric and positive definite, so we will use the PCG method to solve it.

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 Γ indicates dofs living on the interface, Π and ∆ are instead used to distinguish dofs of Γ 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 • refers to a continuous space, the • 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 V Γ , V Π 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 V ∆ is the direct sum of the subdomain dual interface velocity spaces ∆ , 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 V Γ 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 (i) Γ is the operator which maps a function in the space V Γ to its component in ∆ as the operator which maps the space V Γ to its dual component in the space ∆ , and it is a map from V Γ into V Γ .The relationships among the previous spaces and operators are summarized in the following diagram: In order to define certain scaling operators, which will be used in the definition of the BDDC preconditioner, see (39) , we introduce a positive scaling factor δ † i (x) for the nodes on the interface Γ i of each subdomain Ω i .For the type of problem we will use in the numerical experiment (incompressible Stokes problems), we simply define the δ † i (x) as the pseudoinverse counting functions, so: 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 ∆ , only one for row, by the corresponding scaling factor δ † i (x).We construct also the scaled operator R D,Γ as the direct sum of R Γ,Π and R (i) D,∆ .After the change of basis, the interface velocity Schur complement S Γ is defined on the partially assembled interface velocity space V Γ by: given Defining by R Γ the operator that maps the space V Γ into the product space V Γ associated with the set of subdomains, we observe that S Γ can be obtained from the Schur complements S (i) Γ by assembling only the primal interface velocity part, i.e. as As we saw before (30) the global interface Schur operator S Γ is obtanied by fully assembling the S Γ across the subdomain interface, therefore it can be also obtained from S Γ by further assembling the dual interface velocity part, S Γ = R T Γ S Γ R Γ .So we need to define an operator B 0Γ , which maps the partially assembled interface velocity space V Γ into F 0 , the space of right hand sides corresponding to Q 0 , and it is obtained from B 0Γ by assembling the dual interface velocity part on the subdomain interfaces, i.e.
we can write S, the operator of the global interface problem (28), as The preconditioner for solving the global saddle-point problem ( 28) is where we have defined and so we have the BDDC preconditioned problem: find What we need in our implementation is to determine the action S −1 q for any given q = (q Γ , q 0 ) ∈ F Γ × F 0 , so we have to solve the linear system Given the definition of S Γ in (5.3), we have that solving (5.10) is equivalent to solve where Π .Now using a block factorization we obtain where R ∆,i maps ∆ , the set of right hand sides corresponding to V ∆ .The matrix S CC , relatively to the primal constraints, has to be completely assembled in this way where we have defined the maps from Finally we define the matrix where R Π0 is the map between the space F Γ × F 0 and F Π × F 0 .

Theoretical estimates
We now present an estimate for the eigenvalues of the preconditioned operator M −1 S, following the theory developed in [21] and adapting it to our VEM formulation.We can do so because the space V Γ 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 (i) Γ , defined in (28), are symmetric and positive definite.
Proof.We know from (30) that the Schur complement related to the velocity is defined by: given w By the coercivity of a(•, •) we know that the matrices are symmetric and positive definite and so the left two by two upper block of the lefthand-side of (48) has the same number of negative eigenvalues of the all matrix.Now, the left-hand-side matrices of (48) are congruent to: 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 and b (i) the restrictions to subdomain Ω i of the bilinear forms a, a h and b, respectively.Then, we introduce the |.| S (i) Γ and |.| S Γ seminorms defined by and a norm and a seminorm on the space with consequently the norm .1/2,Γ and seminorm |.| 1/2,Γ defined on the space V Γ by There exist positive constant c 1 and c 2 , independent of H, h and the shape of subdomains, such that where β 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 v Γ ∈ V Γ , we define the operators The above condition uniquely defines S and T .Now given v Γ ∈ V Γ , let v H Γ ∈ V be the discrete harmonic extension of v Γ , i.e. the unique function in V which equals v Γ on Γ and satisfies ∀i = 1, ..., N : By the stability of the discrete harmonic extension and the stability of the discrete bilinear form a h [3], we have on each subdomain: 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 (v Γ ), S(v Γ )) = 0, we have: Applying ( 14) on the subdomains, we have, for some c > 0: Applying Cauchy-Schwarz to the first term in (54) and using (55), we have: Finally, from (53) and summing on the subdomains 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 : and then: The operators S Γ and S Γ , given in ( 31) and (36), are both symmetric and positive definite, because of the Dirichlet boundary conditions on ∂Ω and provided that sufficiently many primal constraints are chosen.We can then define the S Γ and S Γ norms on the spaces V Γ and V Γ by We then define two spaces, whose utility is that, restricted to such spaces, the interface problem operators S of ( 28) and S of ( 42) are positive semi-definite.As in [21], we give the following: Definition 5.1.Given the discrete spaces V Γ and V Γ , we define the two subspaces Lemma 5.3.The interface operator S of (28), restricted to the subspace V Γ,B × Q 0 is positive semi-definite.The same is true for S of (36 We define the S and S seminorms on the benign subspaces Now we define an average operator E D = R R T D , which maps V Γ × Q 0 , with generally discontinuous interface velocities, to elements with continuous interface velocities in the same space.For any where , provides the average of the interface velocities across the interface Γ. Recalling that we can split 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.
where n is the outward normal of ∂Ω i .We can equivalently write Assumption 2. There exists a positive constant C, which is independent of H, h and the number of subdomains, such that With these two assumptions, we have the following results (proof of 5.4 in [21]): 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 where β is the inf-sup stability constant of (14).
We have from the definition of the S-seminorm, that where the last inequality follows from Lemma 5.2.We have, from Assumption 2 and Lemma 5.2 Consequently we have We have the following lemma (proof in [21]): Lemma 5.6.Any vector of the form u = (0, p 0 ) ∈ V Γ,B × Q 0 is an eigenvector of the preconditioner operator M −1 S with eigenvalue equal to 1.
Theorem 5.1.Let Assumptions 1 and 2 hold.The preconditioned operator M −1 S is then symmetric, positive definite with respect to the bilinear form •, • S on the benign space V Γ,B × Q 0 .Its minimum eigenvalue is 1 and its maximum eigenvalue is bounded by Here, C is a constant which is independent of H, h and the number of subdomains, and β is the inf-sup stability constant defined in (14).
Proof.We know from Lemma 5.6, that any vector of the form u = (0, p 0 ) ∈ V Γ,B × Q 0 is an eigenvector of the preconditioned operator M −1 S with an eigenvalue equal to 1.It is sufficient to find lower and upper bounds of the quotient M −1 Su, u S / u, u S , for any u = (u Γ , p 0 ) ∈ V Γ,B × Q 0 , where u Γ is non zero and therefore u, u S > 0.
Lower bound: We have from the fact that From the Cauchy-Schwartz inequality and the fact that S = R T S R, we find that Therefore from (65) and (66), Since, we obtain, from equations (67) and (68), that u, u S ≤ u, M −1 Su 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.
R and by using Lemma 5.5, we have Therefore from equation (68), we have Using the Cauchy-Schwarz inequality and equation (70), we have This gives and the upper bound of the theorem.

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 Γ ij , which is shared by a pair of subdomains Ω i and Ω j , we impose for a fixed selection of the normal n ij of Γ 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 Γ ij , i.e.
By the definition of the average operator E D,∆ we have that the average interface velocity is ∆ ) on each edge and hence 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 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 v Γ ∈ V Γ we have: with C positive constant independent of H, h and the number of subdomains.

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 Ω = [0, 1] × [0, 1], applying homogeneous Dirichlet boundary conditions on the whole ∂Ω.We choose the load term f by imposing that the analytical solution is u(x, y) = − sin(πx) sin(πx) sin(2πy) sin(πy) sin(πy) sin(2πx) , p(x, y) = sin(πx) − sin(πy).
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 .
(a) QUAD meshes.In the following tables, we report the number of iterations to solve the global interface saddle-point problem (28) 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).
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 subdo-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.
(a) QUAD meshes.main 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.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.

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.

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.

Figure 1 :
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.

Figure 2 :
Figure 2: Examples of the different type of meshes we consider in our numerical experiments.
The primal space consists of the subdomain vertices and two basis functions per subdomain edge.The number of subdomains is fixed to N = 16.The primal space consists of the subdomain vertices and two basis functions per subdomain edge.The ratio H/h is fixed to 4.
The primal space consists of the subdomain vertices and one basis function per subdomain edge.The number of subdomains is fixed to N = 16.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 :
Figure 3: Plots of the spectral condition number (κ 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.

Figure 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.

Figure 5 :
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.

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.

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.