STABLE DISCRETIZATIONS AND IETI-DP SOLVERS FOR THE STOKES SYSTEM IN MULTI-PATCH IGA

. We are interested in a fast solver for the Stokes equations, discretized with multi-patch Isogeometric Analysis. In the last years, several inf-sup stable discretizations for the Stokes problem have been proposed, often the analysis was restricted to single-patch domains. We focus on one of the simplest approaches, the isogeometric Taylor–Hood element. We show how stability results for single-patch domains can be carried over to multi-patch domains. While this is possible, the stability strongly depends on the shape of the geometry. We construct a Dual-Primal Isogeometric Tearing and Interconnecting (IETI-DP) solver that does not suffer from that effect. We give a convergence analysis and provide numerical tests


Introduction
Isogeometric Analysis (IgA) was introduced in [14] as a technique for discretizing partial differential equations (PDEs); see also [8] and references therein. The original idea is to improve the integration of simulation and computer aided design (CAD), compared to the classical finite element (FEM) simulation. This is achieved by representing both the computational domain and the solution of the PDE as linear combination of tensor-product B-splines or non-uniform rational B-splines (NURBS). Simple computational domains can be parameterized using a single geometry mapping. More complicated domains are usually composed of multiple patches, each parameterized with its own geometry mapping. Such domains are called multi-patch domains. We are interested in fast solvers for the Stokes system, discretized using IgA on multi-patch domains.
For the discretization of the Stokes equations, we need inf-sup stable discretizations. Several inf-sup stable elements from the FEM world have been generalized for the IgA framework, like Nédélec, Raviart-Thomas and Taylor-Hood elements, cf. [4,6,9]. These methods have in common that the same grid is used both for the velocity and the pressure; an alternative approach based on different grids for velocity and pressure is the subgrid approach, cf. [4]. In this paper, we focus on the generalized Taylor-Hood element. A stability estimate was proven in [4] for tensor-product B-splines and it was later extended to hierarchical splines in [3]. The analysis provides lower bounds for the inf-sup constant that is independent of the grid size. Numerical obviously cheaper to construct and to evaluate, we propose and analyze it in this paper. To the knowledge of the authors, this preconditioner has not been proposed before.
We give a condition number bound for the Schur complement formulation of the IETI-DP solver, preconditioned with the scaled Dirichlet preconditioner (see Thm. 5.12). This analysis uses many results that have been developed in [22] for the Poisson problem. The analysis is explicit with respect to grid sizes, the patch diameters, the spline degree, and the inf-sup constants for the local problems. Also this seems to be novel, compared to previous results, like [17], which depended on the inf-sup constant for the global problem. Numerical experiments for the proposed method are provided, but we also refer to [25], where alternative choices of the primal degrees of freedom and alternative setups of the scaled Dirichlet preconditioner are numerically tested.
The remainder of this paper is organized as follows. We present the model problem in Section 2. In Section 3, we introduce an inf-sup stable discretization for multi-patch domains and prove the stability. A IETI-DP solver is proposed in Section 4, which is analyzed in the subsequent Section 5. We conclude the main part of the paper with Section 6, where the results from numerical experiments are presented and analyzed. In Section 7, we give some conclusions and give some outlook on further challenges. The Appendix A contains some of the proofs.

The model problem
As model problem, we consider the Stokes equations with homogeneous Dirichlet boundary conditions in two dimensions. In detail, the model problem is as follows. Let Ω ⊂ R 2 be an open and bounded domain with Lipschitz boundary Ω. 2 (Ω) and (Ω) denote the standard Lebesgue and Sobolev spaces on Ω. Moreover, The existence and uniqueness of a solution to problem (1) is known for any domain Ω with Lipschitz boundary; for a proof, see, e.g., [2], for further information also [11,19] and references therein. The analysis is based on Brezzi's theorem [5], where one shows that there are constants 0 < ≤ and 0 < ≤ such that one has coercivity (∇u, ∇u) 2 (Ω) ≥ ‖u‖ 2 inf-sup stability and boundedness In (3), we do not explicitly mention that u ̸ = 0. Also formulas with suprema that follow are to be understood in that way. The only non-trivial condition is the inf-sup stability (3). Coercivity (2) is a direct consequence of Friedrichs' inequality (cf., e.g., [21], Lem. 1.31) and boundedness (4) (with = 1 and = √ = √ 2) follows directly from the Cauchy-Schwarz inequality.

Stable discretizations
In the following, we introduce a conforming discretization of the Stokes equations which again satisfies the conditions of Brezzi's theorem. Certainly, equations (2) and (4) carry directly over to conforming discretizations. The story is different for the inf-sup condition, which has to be verified for the discretized problem as well. For single-patch Isogeometric Analysis, such stable discretizations have been introduced previously. After introducing the representation of the computational domain in Section 3.1 and the standard concepts of isogeometric functions in Section 3.2, we replicate the details of the isogeometric Taylor-Hood element, which we use in our further considerations, in Section 3.3. In Section 3.4, we discuss the extension of these results to multi-patch Isogeometric Analysis and the dependence of the inf-sup constant on the shape of the computational domain. In Section 3.5, we present and discuss numerical results that illustrate the dependence of the stability on the shape of the geometry.

Representation of the geometry
We assume that the computational domain Ω ⊂ R 2 is composed of non-overlapping patches Ω ( ) , i.e., the domains Ω ( ) are open and bounded domains with Lipschitz boundary such that where denotes the closure of the set . We need that that the patches form an admissible decomposition, i.e., that there are no T-junctions. This assumption is necessary to allow a fully matching discretization, which is a prerequisite for an 1conforming discretization. Recently, a IETI solver for the Poisson equation was proposed that allows a decomposition including T-junctions, cf. [23]. That approach uses a discontinuous Galerkin method in order to couple the patches. Since we focus on conforming discretizations, we cannot use an analogous approach.
For any patch index , the set Γ ( ) contains the indices ℓ of patches Ω (ℓ) that share an edge with Ω ( ) . The common vertices of two or more patches -that are not located on the (Dirichlet) boundary -are denoted by 1 , . . . , . For each = 1, . . . , , the set ( ) contains the indices of all patches Ω ( ) such that ∈ Ω ( ) . We assume that the number of patches sharing one vertex is uniformly bounded.
Each patch Ω ( ) is parameterized by a geometry mapping which can be continuously extended to the closure of the parameter domain̂︀ Ω. In IgA, the geometry mapping is typically represented using B-splines or NURBS. As usual, the computational methods do not depend on such a representation. We only assume that the geometry mappings are not too much distorted, i.e., that the following assumption holds.
We need one more assumption in order to analyze the inf-sup stability of the global problem. This is a condition that is specific for the analysis of the Stokes equations. The assumption guarantees that the interfaces are bent by uniformly less than 180 degrees. If Assumption 3.3 holds, each subdivision of Ω into patches that does not satisfy Assumption 3.4 can be converted into a subdivision satisfying this condition by (uniformly) subdividing the patches sufficiently often.

Isogeometric functions
On the parameter domain̂︀ Ω = (0, 1) 2 , we choose a B-spline space, which depends on a freely chosen vector of breakpoints for each patch and each spacial direction ∈ {1, 2}, a freely chosen degree parameter p ∈ N := {1, 2, 3, . . .} and a freely chosen smoothness parameter s ∈ {0, 1, . . . , p − 1}. Based on these vectors of breakpoints, we introduce spline spaces of degree p and smoothness s: where P p is the space of polynomials of degree p. For each such set, we choose the basis that is obtained by the Cox-de Boor formula (cf. [8], Eqs. (2.1) and (2.2)); for the application of the Cox-de Boor formula, one uses a knot vector obtained from the vector of breakpoints by repeating the first and the last breakpoint p + 1 times and by repeating all other breakpoints p − s times. Based on these univariate splines, we introduce the corresponding tensor-product spline space as discretization space on the parameter domain̂︀ Ω, and equip it with the standard tensor-product basis. The function spaces on the physical patches Ω ( ) are defined via the pull-back principle, so we define a space of functions The grid sizê︀ ℎ on the parameter domain and the grid size ℎ on the physical patch are defined bŷ︀ where the definition of the latter is motivated by Assumption 3.3. We assume that the grids are quasi-uniform. Note that Assumption 3.3 allows us to relate the norm of the function on the physical patch and the corresponding function on the parameter domain. There is a constant > 0, only depending on the constant from Assumption 3.3, such that Using a standard Poincaré inequality (cf., e.g., [21], Lem. 1.27), we obtain wherê︀ is the Poincaré constant for the parameter domain̂︀ Ω = (0, 1) 2 . This means that the Poincaré constant for Ω ( ) only depends on and . A completely analogous result for the Friedrichs' inequality is straight forward: For all patches Ω ( ) , where at least one edge is located on the (Dirichlet) boundary, we have using a standard Friedrichs' inequality (cf., e.g., [21], Lem. 1.31) wherê︀ F is the Friedrichs' constant for the parameter domain̂︀ Ω.

Stable discretizations for the single-patch case
As discretization space for the single-patch case, we use the isogeometric Taylor-Hood element, as proposed in [4]. It uses the same grid for all velocity components and for the pressure and can be defined based on any underlying spline degree parameter p ∈ N and any underlying smoothness s ∈ {0, . . . , p − 1}.
The idea of the isogeometric Taylor-Hood element is to use splines of degree p + 1 and smoothness s, which vanish on the (Dirichlet) boundary, for the velocity and splines of degree p and smoothness s with vanishing mean value for the pressure. Our approach is to use these spaces for each of the patches Ω ( ) , however we have to modify the spaces accordingly. So, Dirichlet boundary conditions are not to be imposed on Ω ( ) , but only on Γ As basis for the spacê︀ V ( ) , we choose the basis functions of the standard tensor-product B-spline basis that vanish on̂︀ Γ ( ) D . Their images under the geometry function G form the basis for the space V ( ) . The bases for︀ ( ) and ( ) are defined analogously. In [4], it was shown that the isogeometric Taylor-Hood element is inf-sup stable if the grid is sufficiently fine. Certainly, this only holds if we have boundary conditions on all of Ω ( ) and the averaging condition locally, i.e., we have where the inf-sup constant is independent of the grid size ℎ , but it depends on G and the constant from Assumption 3.5. Since the discretization is conforming, coercivity (2) and boundedness (4) are also satisfied for the discretion problem.
Remark 3.6. Extensive numerical experiments indicate that the constant is independent of the spline degree p if at least one inner knot is used in each parametric direction. However, at the time of writing, no such proof is known to the authors.

Stable discretization in the multi-patch case
In this section, we introduce the global function spaces V ⊂ [ 1 0 (Ω)] 2 and ⊂ 2 0 (Ω). Since we set up a conforming discretization, we need that the space V is continuous. To be able to set up a continuous global function space, we need that the discretization is fully matching, i.e., that the following assumption holds.
Assumption 3.7. For every interface Γ ( ,ℓ) between two patches, the following statement holds true. For any basis function in the basis for V ( ) having support on Γ ( ,ℓ) , there is exactly one basis function in the basis for V (ℓ) such that they agree on the interface Γ ( ,ℓ) .
This assumption holds if the spline degree, the vector of breakpoints and the geometry mapping agree on all common interfaces. For each of the matching basis functions in Assumption 3.7, we set the corresponding coefficients to have the same value. In this way, we obtain an 1 -conforming discretization space. Note, this is not done for the pressure space since it only needs to be 2 -conforming. We can now state the overall discretization space. For the velocity, we use The discretized Stokes problem reads as follows. Find (u, ) ∈ V × such that Now, we prove an inf-sup stability result for the multi-patch case, which uses the inf-sup stability of the continuous problem, i.e., (3), and the inf-sup stability result for the single-patch case, i.e., (8).
Before we can prove the main inf-sup result, we need some auxiliary results. Note that is the direct sum of the space of function with zero average on each patch, and the space of patchwise constant functions. First, we state the existence of a Fortin operator.
Lemma 3.8. There exists an operator Π F : and where F ≥ 1 is a constant that only depends on the constants from the Assumptions 3.2-3.4.
The proof of this lemma is given in the Appendix A. We now show an inf-sup estimate for the pressure space of patchwise constants. Lemma 3.9. We have the inf-sup estimate where is as in (3) and F is as in Lemma 3.8.
Proof. Let 1 ∈ 1 be arbitrary but fixed. From the continuous inf-sup condition (3), it follows that there exists By setting u := Π F v and using Lemma 3.8 and (13), we get which finishes the proof.
Using this inf-sup result and the patchwise inf-sup result (8), we can show a global discrete inf-sup result.
Theorem 3.10. Let V × be the generalized Taylor-Hood space as defined in this Section. We have the inf-sup result where is as in (3), as in (8), as in (4) and F as in Lemma 3.8.
Before we prove this theorem, we give some remarks.
Remark 3.11. ℎ is independent of the grid sizes ℎ since the local inf-sup constants and the Fortin constant F are independent of ℎ and the inf-sup constant for the continuous problem is inherently independent of the discretization. Robustness in the spline degree is obtained if the local inf-sup constants are robust in the spline degree, which is an unproven conjecture for the isogeometric Taylor-Hood element, see Remark 3.6.
Remark 3.12. We observe that the local inf-sup constants and the Fortin constant F are independent of the shape of the overall domain Ω (both only depend on the parameterization and thus also of the shape of the individual patches and the maximum number of patches meeting in one vertex), so any shape-dependence observed in numerical results is due to the shape-dependence of , the inf-sup constant of the continuous problem and thus inherent to the Stokes equations themselves.
Remark 3.13. Moreover, we observe that Theorem 3.10 is not restricted to the isogeometric Taylor-Hood element. So, the proofs can be applied to any conforming, locally inf-sup stable pair of discretization space, where the velocity space is fully matching at the interfaces and where all the bi-quadratic functions are in the underlying spacê︀ V ( ) .
The space 2 0 (Ω) has to be replaced by the space 2 (Ω) since the pressure can be uniquely determined in this framework. Also for this setting, the existence of a Fortin operator (Lem. 3.8) can be shown. The construction of the operator (as presented in the Appendix A) is the same as for pure Dirichlet conditions with the sole difference that also correction functions ( ,ℓ) associated to the edges that are located on the Neumann boundary have to be taken into account. Having Lemma 3.8, all other results follow naturally.
Proof of Theorem 3.10. Let ∈ be arbitrary but fixed and let = 0 + 1 with 0 ∈ 0 and 1 ∈ 1 . From (8), we know that there are non-zero functions The suprema in the inf-sup conditions are scaling invariant, so we can restrict ourselves to Note that u vanishes on the interfaces between the patches. Summing up, we obtain By applying integration-by-parts patchwise and using ∇ 1 = 0 and u| Ω ( ) = 0, we obtain that By combining this with (15), we obtain Lemma 3.9 gives together with boundedness (4) Using the triangle inequality, we have further By adding ( + −1 F ) times inequality (16) and̂︀ times inequality (17), we finally get where we make use of F ≥ 1 and ,̂︀ ≤ .

Numerical exploration of the inf-sup constant
In this subsection, we present numerical experiments that illustrate the dependence of the discrete inf-sup constant ℎ on the grid size, the spline degree and the shape of the geometry. This is done by deriving the condition number of the (negative) Schur complement, preconditioned with the inverse of the mass matrix for the pressure. The relation between this condition number and the inf-sup constant is where ℎ ≈ 1 is the discrete version of the boundedness constant from (4). As computational domain, we consider the Yeti-footprint, consisting of 21 patches, and domains obtained by combining a few of the patches. The Yeti-footprint is depicted in Figure 1; the patches are represented by different colors. The lines within each patch represent the coarsest (ℓ = 0) grid on each patch. We obtain finer grids by performing ℓ uniform refinement steps (ℓ = 0, 1, 2, 3) and test for various choices of the spline degree parameters p (p = 1, 2, . . . , 5).
The computed condition numbers are displayed in Table 1. The left table shows the results for the single-patch domain Ω (1) and the right table shows the results for the full Yeti-footprint Ω (1) ∪ · · · ∪ Ω (21) . Due to the size the full Yeti-footprint, we could not calculate the condition numbers for ℓ = 3, so these are left out. As we see from the tables, the constants depend neither on the grid size (this is predicted by the theory), nor on the spline degree (cf. Rem. 3.6). We notice that the condition numbers are significantly larger for the full domain, which is a property of the geometry (cf. Rem. 3.12).
To explore the dependence on the shape of the domain further, we calculate condition numbers for partial Yeti-footprints Ω (1) ∪ · · · ∪ Ω ( ) with = 1, 2, 3, 4, 5, 21. The corresponding condition numbers are computed for ℓ = 2 and p = 2 and presented in Table 2. We observe that the condition numbers increase steadily from

A IETI-DP solver for the Stokes system
In this section, we outline the setup of the proposed IETI-DP method for solving the discretized Stokes problem (9). For the IETI-DP method, we have to assemble the variational problem locally. So, the still uncoupled problem is to find Here, we use the notation ! = to remind ourselves that the coupling is still missing. By discretizing these bilinear forms using the tensor-product bases, we obtain linear systems Here and in what follows, underlined quantities refer to the coefficient representations of the corresponding functions. We first represent the spaces V ( ) as a direct sum where V ( ) I is spanned by the basis functions which vanish on the interfaces and V ( ) Γ is spanned by the remaining functions, i.e., the functions that are active on the interfaces (which includes the primal degrees of freedom). Assuming a corresponding ordering of the basis functions, we have Analogous to the case of the Poisson problem, the local systems correspond to pure Neumann problems, unless the corresponding patch contributes to the Dirichlet boundary. These systems are not uniquely solvable since the constant velocities are in the null space of ( ) . To ensure that the system matrices of the patch-local problems are non-singular, we introduce primal degrees of freedom, whose continuity across the patches is enforced strongly (continuity conditions): -the function values of the velocity at each of the corners of the patch, i.e., -the integrals of the normal components of the velocity, i.e., Since the constraint (20) is vector-valued, there are actually 2 primal degrees of freedom for each corner. Overall, for patches that do not contribute to the Dirichlet boundary, there are 12 primal degrees of freedom related to the continuity conditions. The matrix ( ) C evaluates the primal degrees of freedom associated to the patch Ω ( ) ; thus the relation Additionally, we introduce primal degrees of freedom for the pressure in order to be able to realize the condition that the average pressure vanishes (averaging conditions). This is done by fixing the average pressure on each patch individually to zero and by allowing patchwise constant pressure functions in the primal problem. The matrix ( ) A evaluates the average pressure on the patch, i.e., Thus, the relation The corresponding Lagrangian multipliers are denoted by The continuity of the velocity between the patches is enforced by the matrices ( ) . The term evaluates to a vector containing the differences of the coefficients of any two matching (cf. Assumption 3.7) basis functions. Here, we do not include the vertex values (see Fig. 2) since these are primal degrees of freedom anyway. Note that the constraint matrices ( ) are redundant to the condition on the integrals of the normal components of the velocity over the edges. The relation guarantees the continuity of the velocity function across the patches.
Moreover, we introduce the primal problem, i.e., the global problem for the primal degrees of freedom. We use a ( ) -orthogonal basis for the primal degrees of freedom. This basis is represented in terms of the basis  We define the system matrix, right-hand side and jump matrix for the primal problem as So far, we have a pressure averaging condition for the patch-local problems and the primal problem allows for patchwise constant pressure modes. So, in order to obtain unique solvability of the problem, we need to add a global condition that guarantees that the average pressure vanishes. So, we augment the primal system and obtain the following primal system where Π ∈ R 1×( Π+ ) is such that the relation Π x Π = 0 guarantees that the average of the pressure vanishes. Correspondingly, we define¯Π := (︀ Π 0 )︀ . So, we are finally able to write down the overall IETI-DP system: ⎛ . . .
For solving this linear system, we take the Schur complement with respect to the Lagrange multipliers . This means that we solve¯= , where¯: The patch-local solutions are then recovered bȳ The final solution is then obtained by distributingx Π to the patches using the matrices Ψ ( ) . We solve the system (23) with a conjugate gradient solver. (We will show in Sect. 5 that¯is indeed positive semidefinite.) The conjugate gradient solver is preconditioned with the scaled Dirichlet preconditioner. Usually, the scaled Dirichlet preconditioner refers to the local solution of Dirichlet problems of the corresponding differential equation, which would mean that we should consider the Stokes equations. However, this is not necessary. Indeed, the local Dirichlet problems are solved to realize the 1/2 -norm. For this purpose, it is sufficient to only consider local Dirichlet problems of the Poisson equation. A recent numerical study, cf. [25], has shown that this approach is not only simpler, but also leads to better convergence behavior. Thus, we define the local Schur complements via Then, the scaled Dirichlet preconditioner is given by where := 2 is set up based on the principle of multiplicity scaling.

Condition number analysis for the IETI solver
The convergence rates of the conjugate gradient solver are estimated based on the condition number of the preconditioned system sD¯. Following the framework introduced in [18], we rewrite the whole problem equivalently as a formulation only living to the skeleton. First, we define the skeleton formulation associated to the main saddle point matrix ( ) via Note that the inverse is well-defined due the constraint on the average pressure. Next, we define corresponding function spaces. Let W := W (1) × · · · × W ( ) with be the skeleton space and ℋ : represents the functions that satisfy the primal degrees of freedom homogeneously. The spacẽ︁ W, in which all the approximate solutions live in, is given bỹ︁ where w = (w (1) , · · · , w ( ) ). The first two lines in definition of this space refer to continuity constraints on the velocity, both on the corner values of the velocities and the integrals of the normal components of the velocity on all of the edges. Since the original formulation of the discretized Stokes problem (9) uses a continuous discretization space for the velocity, these continuity conditions are obviously satisfied by the solution. However, the definition also introduces a third class of constraints of the form ∫︁ i.e., that the inflow equals the outflow on each patch. Let u = (u (1) , . . . u ( ) ) be the exact solution with skeleton representation w ( ) := u ( ) | Ω( ) and let Be the characteristic functions. By integration by parts and since the functions˜Ω( ) ( ) : are contained in , we know from (9) that ∫︁ which vanishes due to the Dirichlet boundary conditions. This shows that (27) holds for the solution. The spacẽ︁ W Π is the orthogonal complement tõ︁ W Δ iñ︁ W, i.e., we definẽ︁ , the corresponding underlined symbol, here w ( ) , denotes the representation of the corresponding function with respect to the basis for the space W ( ) . So, functions in the all of these spaces are represented with respect to the same basis.
For the analysis, we introduce the following lemma that allows to write expressions involving inverses of matrices as suprema.
First we note the equivalence of ( ) and ( ) .
Lemma 5.2. We have for all skeleton functions w ( ) ∈ W ( ) , where ( ) is as defined in (24) and ( ) is as defined in (26).
Proof. Recalling (21), we have that where ( ) is the mass matrix that is obtained from discretizing (·, ·) 2 (Ω ( ) ) with the basis functions in the basis for ( ) and 1 = (1, · · · , 1) ⊤ is a vector of ones of the corresponding size. Certainly, we have Using the definition of ( ) and by expanding the products, we obtain IΓ is as defined in (24) and We first observe that Next, we observe that is the overall divergence and ℋ is the matrix representation of the discrete harmonic extension. This means that we have using Lemma 5.1 Using the Gauss rule and the choice ( ) = 1, we obtain and thus the desired bound from below.
Proof. Using the definition of¯( ) , the block structure of¯( ) , the fact that we can reorder the entries in¯( ) , block Gaussian elimination and (26), we obtain Since w ( ) ∈̃︁ W Proof. Since ( ) A = (0, · · · , 0, 1, 0, · · · , 0), where the non-zero coefficient is in the -th column, we immediately obtain from the definition (22) (32) Since the matrix¯( ) is non-singular, the system has a unique solution. Choose which is possible since ( ) C evaluates the primal degrees of freedom and the integrals of the normal components of the velocity variable on each edge are primal degrees of freedom. So, ( ) is just such that sum over the edges. Using the Gauss rule, we further have which shows Analogously to Lemma 5.3, we show that the operator¯Π corresponds to taking the maximum iñ︁ W Π . Before we give a proof, we introduce some useful notation by collecting local contributions to global matrices and vectors.
Notation 5.5. The matrix is a block row matrix containing the corresponding patch-local contributions, like := (︀ (1) · · · ( ) )︀ . The matrices Ψ ΓC and Ψ ΓA are block column matrices containing the corresponding patch-local contributions. The matrices , , Γ and are block diagonal matrices containing the corresponding patch-local contributions. The vectors, like w, are the corresponding block vectors containing the patch-local contributions.
Lemma 5.6. The matrix¯Π is symmetric and positive semidefinite and satisfies  (33) and (26), we have Using Lemma 5.4, we have where is a block-diagonal matrix containing where W C := {w C : Π = 0 ⇒ ⊤ ⊤ Ψ ΓC w C = 0}. Define w := Ψ ΓC w C and let w = (w (1) , · · · , w ( ) ) ∈̃︁ W be the function associated to the coefficient vector w. Observe that the condition Since functions iñ︁ W satisfy homogeneous Dirichlet boundary conditions and since they satisfy a continuity condition on the averages of the normal components of the velocity, we also have By combining these results, we obtain that ∫︀ Ω ( ) w ( ) · n ( ) d = 0 for all = 1, . . . , . Since w is in the image space of the primal basis functions, we also know that w satisfies the remaining conditions in the definition of︁ W and that it is orthogonal tõ︁ W Δ . This shows w ∈̃︁ W Π . The reverse direction, i.e., that for each w ∈̃︁ W Π , there is some w C ∈ W C with w := Ψ ΓC w C is straight forward. So, we obtain which is what we wanted to show. This representation immediately shows that the symmetric matrix¯Π is positive semidefinite.
From the Lemmas 5.3 and 5.6, we immediately obtain that the matrix¯is symmetric positive semidefinite and that the following result holds.
Lemma 5.7. The identity holds for all .
Proof. Using orthogonality, the Cauchy-Schwarz inequality and Lemmas 5.3 and 5.6, we have The Cauchy-Schwarz inequality is satisfied with equality if the corresponding terms are equal, i.e., Let W * ⊂̃︁ W be the subset of functions that satisfy these conditions. Due to scaling invariance and the fact that the Cauchy-Schwarz inequality is satisfied with equality and W * ⊂̃︁ W, we have which finishes the proof.
This lemma is standard, for a proof, see, e.g., Lemma 4.16 of [22]. Since we exclude the corners, this statement is standard, see, e.g., [18].
holds, where only depends on the constants from the Assumptions 3.3 and 3.5.
Proof. Within this proof, we write if there is a constant that only depends on the constants from the Assumptions 3.3 and 3.5 such that ≤ . Let v ∈̃︁ W and w ∈ W with coefficient representations such that w = −1 ⊤ v. Using Theorem 4.2 of [22] (which depends on the constants from Assumptions 3.3 and 3.5), we obtain where we apply Theorem 4.2 of [22] to both components of the velocity variable separately. By applying Lemma 4.15 of [22] (which depends on the constants from Assumptions 3.3 and 3.5) to both velocity components separately, we further obtain where Λ := 1 + log p + max =1,..., log ℎ and | | ∞ 0 (Γ ( ,ℓ) ) := inf ∈R ‖ − ‖ ∞ (Γ ( ,ℓ) ) . Using Lemma 5.8 and the triangle inequality, we further obtain )︂ .
By applying Theorem 4.2 and Lemma 4.14 of [22] (which depend on the constants from Assumptions 3.3 and 3.5) again to both velocity components, we arrive at )︂ .
Using a Poincaré inequality (6), we have from which the desired result follows.
If standard Krylov space methods are applied to the singular matrix , preconditioned with a non-singular preconditioner sD , all iterations live in the corresponding factor space. The convergence behavior is dictated by the essential condition number of sD , cf. Remark 23 of [18]. The essential condition number for a positive semidefinite matrix is the ratio between the largest eigenvalue and the smallest positive eigenvalue.
Theorem 5.12. Provided that the IETI-DP solver is set up as outlined in the previous section, the condition number of the preconditioned system satisfies where is as in (4), is as in (8) and is a constant that only depends on the constants from the Assumptions 3.3 and 3.5.
Proof. Within this proof, we write if there is a constant that only depends on the constants from the Assumptions 3.3 and 3.5 such that ≤ . For an upper bound, we use Lemma 5.2, (37), Lemma 5.11 and the fact that has a full rank, where is a vector space of the corresponding dimension and 2 := p(1 + log p + max log ℎ ) 2 . This provides an upper bound for the eigenvalues of −1 sD¯.
Next, we estimate the smallest non-zero eigenvalue. We consider the following generalized eigenvalue problem: = are eigenvectors with eigenvalue 0. Since all eigenvectors with non-zero eigenvalue are −1 sD -orthogonal to the eigenvectors with eigenvalue 0, these eigenvectors have to be in 1 . Using 1 ∈ 1 ∖{0} and Lemma 5.10, we have √︁ Using Lemma 5.9, we know that v := −1 ⊤ w ∈̃︁ W, so using Lemma 5.2 we further obtain which provides a lower bound for the positive eigenvalues. Having these eigenvalue bounds, we immediately obtain the desired bound on the essential condition number.
With the observations of Remark 3.14, the results for the IETI solver can be generalized to mixed boundary conditions.

Numerical results
In this section, we present numerical results that illustrate the efficiency of the proposed IETI-DP solver. In Section 6.1, we present results for domains that have been previously considered in IgA and which are fully covered by the presented theory. In Section 6.2, we present a more physical test example using boundary conditions that go beyond the model problem considered for the theory. We solve this problem on two computational domains: a B-spline approximation of a quarter annulus with 64 patches and the Yeti-footprint, where we split the patches uniformly such that the domain has 84 patches (rather than 21), see Figure 3. On each patch, we obtain a discretization space by performing ℓ = 1, 2, . . . uniform refinement steps. The problem is discretized as outlined in Section 3 and solved using the IETI-DP method proposed in Section 4. The conjugate gradient solver is started using a random initial guess and stopped when  the Euclidean norm of the residual vector is reduced by a factor of = 10 −6 compared to the Euclidean norm of the initial residual vector. While performing the conjugate gradient method, we also estimate the condition number based on the underlying Lanczos iteration. The local linear systems are solved using a sparse LU-solver. All experiments have been implemented using the G+Smo library 1 and have been performed on the Radon1 cluster 2 in Linz.

Results for quarter annulus and Yeti-footprint
We present the iteration counts and the condition numbers in Tables 3 (Quarter annulus) and 4 (Yetifootprint). We see that the iteration counts and the condition numbers grow about linearly in ℓ, which corresponds to a growth like log ℎ , which is slower than (log ℎ ) 2 , predicted by the theory. The growth in the spline degree parameter p seems to be linear or sublinear. Here, the theory does not tell the complete story since we do not actually know how the inf-sup constant depends on the spline degree. For both domains, the condition numbers and the iteration counts are very satisfactory, keeping in mind the results from Section 3.5.

Flow through a rectangle with an obstacle
In this section, we consider a stationary flow through a rectangle with a circular hole. This domain consists 11 patches, see Figure 4. The first four patches, which are adjacent to the circle are parameterized using NURBS. The remaining patches are parameterized using a standard affine mapping, represented as tensor-product Bspline mappings. Even though some patches are parameterized NURBS, we use tensor-product B-splines to set up the discrete function spaces. Starting from a coarsest level with no interior knots, we perform ℓ = 1, 2, 3, 4, 5 uniform refinements to obtain the grid used for the simulation.
The differential equation is set up as follows. We have a zero source term, f = 0, and the boundary conditions are chosen as follows. On {30}×(−2, 2), we choose an outlet boundary condition, that is, we use the homogeneous Neumann condition ∇u · n + n = 0.
This problem does not coincide with our model problem (1) as we have a Neumann boundary condition. Using the Neumann condition, the average of the pressure is uniquely solvable in 2 (Ω). So, we omit the condition on the average pressure. We solve the resulting system as proposed in Section 4. The only difference is that we now do not average the overall pressure, that is, we no longer enforce Π x Π = 0. In Table 5, we present iteration counts and estimated condition numbers for various refinement levels ℓ and spline degrees p. These numbers behave similarly to those presented in Tables 3 and 4. In Figure 5, we present a reconstruction of the solution for ℓ = 4 and p = 2. Note that we changed the sign of the pressure such that it is positive.

Conclusions and outlook
In this paper, we have shown how bounds on the inf-sup constant, previously derived for single-patch isogeometric analysis, can be carried over to the multi-patch case. Moreover, we have proposed and analyzed a IETI-DP solver, where the provided condition number bounds only depending on the inf-sup constants for the patch-local problems and not on the inf-sup constant for the global problem. This means that only the local geometries, and not the global geometry effect the convergence rates.
There are several challenges for future research, like the extension to three dimensions, which would require further work on discrete harmonic extensions for B-splines in this case. Moreover, extensions to the saddle point formulation of (linearized) elasticity problems would be of interest. There, besides the inf-sup constant, also the constant from the Korn inequality would pose challenges.

Appendix A.
Before we give a proof of Lemma 3.8, we construct an operator Π̃︀ that is similar to the Fortin operator. This construction only requires a few basis functions per patch and is best understood as a variation of well-known techniques in the area of finite elements, where the patches Ω ( ) play the role of elements.
Lemma A.1. There exists an operator Π̃︀ : and hold, wherẽ︀ F > 0 only depends on the constants from Assumptions 3.3 and 3.4.
The estimate (5) finishes the proof.
Note that the operator from Lemma A.1 fails to meet the conditions for a Fortin operator due to the additional term of the form −2 ‖u‖ 2 2 (Ω ( ) ) . For the proof of the desired error bound, we need an approximation error estimate that only needs to decrease with the patch size , not with the grid size ℎ . For the construction of such an approximation error estimate, we use a Scott-Zhang operator, which relies on Poincaré estimates. In the following, we verify that the Poincaré constant for the subdomains ( ) from (A.10) only depends on the constants from Assumptions 3.2 and 3.3. where the constant̃︀ only depends on , which is well-bounded due to Assumption 3.2. The statement (A.11) can be transferred to ( ) by applying the same arguments as for (6). The constants only depend on the Jacobians of G and T , which are bounded due to Assumptions 3.2 and 3.3. This finishes the proof.
We now use (A.12) and (A.13) to show (10). It remains to show (11). We have = 0 for all 1 ∈ 1 by applying (A.2) to w, which concludes the proof.
Proof of Lemma 5.1. Since the statements (28) and (29) can be found in the literature, cf. Chapter II, S 1.1 of [11], we only show (30). Let ∈ 2 be such that it maximizes sup ∈ 2 ( , ) ℓ 2 ‖ ‖ . Observe that is only defined up to scaling. So, we introduce the constraint ⊤ = 1. The first order optimality system for the minimizer then reads as follows: where ∈ R and ∈ R 2 are the Lagrange multipliers. The specific form in the third line in (A.14) is obtained by the fact that we may only consider derivatives in the feasible directions. By multiplying the first line in (A.14) from left with ⊤ , we obtain using the second line in (A.14) that ⊤ ⊤ + + ⊤ ⊤ = 0. Since = 0, we know from the third line in (A.14) that ⊤ ⊤ = 0. This shows = − ⊤ ⊤ = − ⊤ . The third line in (A.14) is satisfied if and only if + ⊤ = 0 for some . From this line, the first line in (A.14) and the fourth line in (A.14), we obtain where we make use of ‖ ‖ 2 = 1 and = − ⊤ . This finishes the proof.