An isogeometric mortar method for the coupling of multiple NURBS domains with optimal convergence rates

We investigate the mortar finite element method for second order elliptic boundary value problems on domains which are decomposed into patches Ωk\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Omega _k$$\end{document} with tensor-product NURBS parameterizations. We follow the methodology of IsoGeometric Analysis (IGA) and choose discrete spaces Xh,k\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$X_{h,k}$$\end{document} on each patch Ωk\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Omega _k$$\end{document} as tensor-product NURBS spaces of the same or higher degree as given by the parameterization. Our work is an extension of Brivadis et al. (Comput Methods Appl Mech Eng 284:292–319, 2015) and highlights several aspects which did not receive full attention before. In particular, by choosing appropriate spaces of polynomial splines as Lagrange multipliers, we obtain a uniform infsup-inequality. Moreover, we provide a new additional condition on the discrete spaces Xh,k\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$X_{h,k}$$\end{document} which is required for obtaining optimal convergence rates of the mortar method. Our numerical examples demonstrate that the optimal rate is lost if this condition is neglected.


Introduction
Our main focus is on the theoretical foundation of an IsoGeometric Analysis (IGA) approach to mixed finite element methods using non-uniform rational B-splines (NURBS) for the finite element discretization. We investigate the mortar finite element method on domains B J. Stöckler joachim.stoeckler@math.tu-dortmund.de 1 Fachgebiet Statik und Dynamik, Brandenburgische Technische Universität Cottbus-Senftenberg, Konrad-Wachsmann-Allee 2, 03046 Cottbus, Germany = K k=1 k ⊂ R 2 which are decomposed into patches k with tensor-product NURBS parameterizations For the original mortar method we refer to [6,7,24]. A comprehensive description of the finite-element discretization using the methodology of IsoGeometric Analysis (IGA) is given in [12]. Prior uses of the mortar method in IGA have also been documented in [2,[15][16][17]. Its centerpieces are discrete spaces X h,k ⊂ H 1 ( k ), which are pushforward of tensor-product NURBS spaces on the parameter domain [0, 1] 2 and which have the same or higher degree given by the parameterization F k . Furthermore, discrete spaces M h,l of Lagrange multipliers are defined for the representation of weak continuity conditions across the interfaces γ l . In [12] the pushforward of polynomial spline spaces of the same or lower degree are proposed.
Our work is an extension of [12] and highlights several aspects which did not receive full attention before.
1. In our numerical experiments based on the method in [12], using parameterizations by quartic NURBS with multiple interior knots and quartic splines as Lagrange multipliers, we observed that the optimal approximation rate claimed in Theorem 6 of [12] is not obtained. In Sect. 7 we show that a simple additional condition on the Lagrange multiplier spaces M h,l is sufficient in order to obtain the optimal approximation rate, see Assumption 2(b) and Proposition 7.5. We give a short explanation here. The typical smoothness assumption u| k ∈ H p k +1 ( k ) leads to optimal L 2 -discretization errors O(h p k +1 k ) in a non-parametric setting, where p k denotes the polynomial degree of the finite element space X h,k , see [6,7,24]. For the parametric setting of IGA, the notion of "bent Sobolev spaces" on patches k was introduced in [3,4], in order to remedy the lack of smoothness of the pushforward of tensor-product NURBS with multiple knots. This leads to optimal rates of the approximation errors where W h,l denotes the trace space of X h,k on the interface γ l . As a third ingredient, the isogeometric mortar method in [12] requires that the approximation error of the normal derivative across γ l by the Lagrange multiplier space M h,l satisfies inf μ∈M h,l ∂u ∂ν l − μ We observe in Sect. 7, that this order is not always achieved, if M h,l satisfies the assumptions of Theorem 6 in [12]. The defect results from the presence of multiple knots in the original (coarse) knot vector for the parameterization F k . By a careful analysis of the trace operator with respect to non-smooth curves, as in the book by Nečas [19], we provide the following result: assume that ξ ∈ (0, 1) is an interior knot of the NURBS parameterization of γ l (i.e. parameterized by F k restricted to a boundary line of [0, 1] 2 ). If ξ has multiplicity 2 or higher, then its multiplicity as a knot for the Lagrange multiplier space M h,l should be increased at least by one, in order to achieve the required approximation error O(h ). The precise formulation is given in terms of an augmented knot vector in Proposition 7.1. 2. The formulation of the mortar method as a saddle point problem involves Lagrange multiplier spaces M h,l for the representation of weak continuity conditions across the interfaces γ l . In order to achieve the optimal approximation order of the mortar method, a good choice for M h,l is the pushforward of the polynomial spline space of the same degree as X h,k , if γ l is an edge of the patch k , and with suitable modifications at both endpoints of γ l , see [12,Section 4.3]. It is a conjecture in [12] that this choice justifies a uniform infsup-condition. In Sect. 4.2 we define a spline space M 1 h,l of the same degree with another type of endpoint modifications, for which we provide the analytical proof of the infsup-inequality in Sect. 6. As a side-effect of the new definition of the space of Lagrange multipliers, we show in Sect. 8 that the sparsity of the mass matrix is increased slightly. 3. In the geometrically conforming case, γ kl = ∂ k ∩ ∂ l is either empty, a vertex or a full edge of both patches. In engineering applications of CAD-tools, the decomposition of can often be geometrically non-conforming and includes Tintersections of the patch boundaries. The setting in [2,12,[15][16][17] allows certain types of T-intersections, but excludes configurations of staircase type, where γ kl is neither an edge of k nor of l . We provide an adaptation of the discrete spaces X h,k and allow full flexibility of designing the multi-patch layout of the geometry.
In order to keep the overhead of notations small, we present our analytical results about the infsup-condition and the a-priori error estimates for a simple class of second order elliptic problems. We proceed to more elaborate models in elasticity in our numerical experiments in Sect. 9. The scope of applications of the IGA mortar method has recently been extended to a class of contact problems in [1,20,22], and the infsupcondition and the a-priori error estimates by [12] provided the theoretical foundation. We believe that our results in Sects. 6 and 7 will be valuable ingredients for further developments in this direction.
We begin with the elliptic problem We assume that α is a uniformly positive definite matrix with entries α i, j ∈ L ∞ ( ) and β ∈ L ∞ ( ) is nonnegative. Then the bilinear form a is coercive. In the weak formulation, we look for u ∈ H 1 0,D ( ) such that We often use Sobolev spaces H s ( ) with smoothness order s ≥ 0 which is not always an integer. The norm and semi-norm are denoted as usual by v s, and |v| s, . We end the introduction by a short outline of our work. We repeat from [12] the geometric description and the general setting of the weak formulation as a saddle point problem in Sects. 2 and 3. Section 4 gives the definitions of the discrete spaces X h,k and the new spaces M 1 h of Lagrange multipliers. Sections 5 and 6 deal with the L 2stability of the mortar projection and the uniform infsup-inequality. Section 7 provides a detailed analysis of the approximation order of the discrete solution. In Sect. 8 we give some information about the implementation, with special emphasis on the mass matrix, and Sect. 9 provides several numerical results for the Poisson problem and for elasticity problems.

NURBS description of the geometry
The geometrical setting is formulated as in [12]. In order to fix the notations, we repeat some of the material in the referenced article. Let 1 , . . . , K ⊂ R 2 define a nonoverlapping decomposition of the domain by curvilinear quadrilaterals (patches), i.e.
The decomposition may be geometrically non-conforming, as we allow T-intersections at the boundaries of the patches. A sketch of such decomposition is given in Fig. 1. Each patch k is parameterized by a homeomorphism F k :ˆ = [0, 1] 2 → k which is bi-Lipschitz. Throughout this article, F k is a bivariate NURBS parameterization F k (ξ 1 , ξ 2 ) = i∈I k C k,iNk,i (ξ 1 , ξ 2 ), ( k × (2) k ⊂ˆ and with positive weights {ω k,i , i ∈ I k }. More specifically,N k,i is defined bŷ where (B p k k,i ) i∈I k is the basis of tensor-product B-splines of polynomial degree p k in both coordinate directions. We assume throughout this article that all knot sequences are open; i.e., the first and last knot have multiplicity p k + 1.
For 1 ≤ j, k ≤ K , k = j, we define the interface γ jk as the interior of the intersection γ jk = ∂ j ∩ ∂ k . We consider only those pairs ( j, k) with non-empty γ jk and renumber the interfaces as γ l , l = 1, . . . , L. For each interface γ l , we choose one of the adjacent patches as master cell m(l) and the other one as slave cell s(l) , so that γ l = ∂ m(l) ∩ ∂ s(l) . This assignment allows an arbitrary choice for each interface, i.e. k can be a master cell for one of its boundary lines and a slave cell for another boundary line. We defineγ l = F −1 s(l) (γ l ), which may be only a portion of a boundary line of [0, 1] 2 , and use the notation F l := F s(l) |γ l for the parameterization of the interface.

Remark 2.1
In a geometrically conforming case, every interface γ l is a full edge of s(l) and m(l) . The pre-image F −1 s(l) (γ l ) is a boundary linê NC1 The interface γ l is a full edge of the slave patch s(l) . As before,γ l is as in (4) and we use the abbreviated formγ l = [0, 1]. This restriction on the choice of the master/slave correspondence was used in [12]. It limits the applicability of the mortar method since staircase patch topologies cannot be treated. Furthermore, an additional error term appears in the a priori estimate, as will be seen in (25). This term is due to the fact that the approximation of the weak solution u by NURBS on the master cell does not interpolate at both endpoints of the interface. NC2 There is no geometrical restriction on the choice of the master/slave correspondence. Then we introduce a C 0 -line as an extension of the ending interface of a T-intersection into the adjacent patch k . This is done by inserting a knot ζ of multiplicity p k (or increasing the multiplicity of an existent knot to p k ) into the relevant knot sequence (1) k or (2) k . It does not change the patch geometry, but affects the initial parameterization F k by using a larger set of NURBS basis functions. The immediate effect on the discrete space X h,k is equivalent to splitting the patch k along the line F k (ζ, t) (or F k (t, ζ )), t ∈ [0, 1], and connecting the control points by shared degrees of freedom. By doing so, we get a setting where only geometrical conforming cases are present. Full flexibility in the choice of the master/slave correspondence is obtained and all forms of multi-patch layouts can be computed. This extension leads to the additional element lines visible in Figs. 12 and 27. It is important to note that the required refinement does not propagate to further patches. No additional error term is introduced in (24), in contrast to the configuration in NC1. After this initial modification has been done for all T-intersections, each interface γ l = F l (γ l ) is a parameterized NURBS curve with an open knot sequence l ⊂γ l = [ξ l,1 , ξ l,2 ] ⊂ [0, 1], and ξ l,1 , ξ l,2 are the parameter values of the endpoints of γ l .
The interface γ l inherits its NURBS parameterization F l :γ l → γ l from the slave cell s(l) . We denote its degree by q l = p s(l) , its knot vector by l ⊂γ l , the polynomial B-splines byB q l l,i , the NURBS weights by w l,i and the univariate NURBS basis functions bŷ As usual, both endpoints of the parameter intervalγ l are knots with multiplicity q l + 1, such that l = {θ l,1 = · · · = θ l,q l +1 < θ l,q l +2 ≤ · · · ≤ θ l,n l < θ l,n l +1 = · · · = θ l,n l +q l +1 } and n l denotes the dimension of the initial NURBS space onγ l before h-refinement.
For later use, we also define the univariate NURBS basis functionsN m l, j of degree p m(l) with respect to the master cell m (l) . The extra superscript m is used here in order to point out the role of the master cell.
The standard inner product on γ l is given by where τ l = |F l | is the length of the tangent vector of γ l . More generally, we also use weighted inner products with positive weight functions ρ l and uniform bounds 0 < c ≤ ρ l ≤ C for all l. We will often choose ρ l =ŵ 2 l τ l • F −1 l , and then obtain The inner products are defined for pairs u, v ∈ L 2 (γ l ), or for u ∈ H 1/2 (γ l ) and v ∈ (H 1/2 (γ l )) . The induced norm is equivalent to the standard L 2 -norm on γ l . We make the following assumptions on the parameterization which are essential for our method.

Remark 2.2 (a)
The methods proposed in this paper can also be used for nonwaterproof geometries, whereby from the mathematical point of view a variational crime is committed. This results in an additional error which does not vanish in the fine limit. However, as shown in [12], the mortar method is robust with respect to these kinds of non-matching interfaces and sufficient accuracy for engineering applications is obtained.
The assumption γ l ∈ C 1,1 is natural in the sense that the tangent τ l and the unit normal ν l are Lipschitz-continuous along γ l . This complies with the practical experience, that vertices on an interface γ l should be treated by splitting the adjacent domains accordingly, thus introducing an additional interface. In fact, this is incorporated in all standard CAD programs. However, modifications by hand could possibly result in a non-smooth interface or such points of C 0 -continuity may already exist in industrial CAD files. Therefore, they should be properly treated by a patch coupling method. If no additional interface is desired, we can keep the patch geometry, but split γ l at the relevant location into two parts γ l,1 , γ l,2 and then follow the same method described in Remark 2.1 for T-intersections of type NC2. Thus our results in this article can be extended to geometries with interfaces containing C 0 -points.

Weak formulation of the saddle point problem
The natural function space used in mixed finite element methods for the Dirichlet problem is the direct product The bilinear form a in (1) is canonically extended to X × X by The continuity and coercivity of this extension of a are obvious. Clearly, X is not a subspace of H 1 0,D ( ). For the formulation of weak continuity conditions on every interface γ l , we define the space of Lagrange multipliers where ν l denotes the outer normal of s(l) along γ l . With the notation [v] l = (v s(l) − v m(l) )| γ l for the jump of v ∈ X across γ l , and with the weighted inner product (5), we define the bilinear form with weight functions ρ l =ŵ the weak solution of (2) is obtained from the solution of the following saddle point problem: This problem has a unique solution (u, λ) ∈ X × M, whose first component u ∈ V satisfies (2). The second component λ ∈ M represents the (weighted) normal flux.
Due to the presence of the weight ρ l in (7), its component λ l ∈ (H 1/2 (γ l )) is The usual stability analysis provides the upper bound The reason why we use a weighted inner product instead of the standard inner product of L 2 (γ l ) will become clear when we consider the numerical computation of the mass matrix in Sect. 8.

Discretization of the saddle point problem using NURBS and B-splines
The discretization of the saddle point problem is obtained by choosing a suitable pair of finite-dimensional spaces X h ⊂ X and M h ⊂ M, where the subscript h represents a vector h = (h k ) k=1,...,K for the spatial resolutions on the patches k . The mixed method for the solution of (7) is formulated as follows: The solution space for u h will be denoted by If the weak solution u of (2) is locally smooth, i.e. u ∈ H 1 0,D ( ) ∩ K k=1 H p k +1 ( k ), we look for error estimates of the form Here we use the mesh-dependent norm on M h as in [11,24] with σ = −1/2. In this article, we follow [12] for the definition of the space X h , but we propose an alternative space M h of Lagrange multipliers.

The space X h defined by NURBS
The initial knot sequences k ⊂ [0, 1] 2 were used for the parameterizations F k : [0, 1] 2 → k of the patches k . (In a geometrically non-conforming case NC2 of Remark 2.1, we assume that the knots ξ l,1 , ξ l,2 are already existent with multiplicity p k in (1) k or (2) k , respectively.) The finite element discretization uses refined knot sequences on each patch independently. Every refined knot sequence h,k = (1) h,k × (2) h,k defines a grid of two open knot sequences in [0, 1]. We use shorthand notations ξ (r ) j = ξ (r ) h,k, j for the knots in coordinate direction r = 1, 2, whenever their association with patch k and h-refinement is clear from the context. The index j has range 1 ≤ j ≤ n (r ) + p k + 1, where n (r ) = n (r ) h,k , r = 1, 2, denotes the dimension of the space of univariate NURBS with knots (r ) h,k . To ensure that all NURBS basis functions are continuous in the patch k , we require that all interior knots in (0, 1) have at most multiplicity p k . The relevant parameters for the h-refinement of patch and the representative mesh-size h k for k is defined as The parameters h (r ) j are useful for the Bernstein-type inequalities in (17). They replace the single stepsizes which are often used for piecewise linear finite elements. The following terminology for univariate knot sequences is a slight generalization of (local) quasi-uniformity as defined in [4, Assumption 2.1].
It is called quasi-uniform of order p, if there exist 0 < c 1 ≤ c 2 such that for all h A family of bivariate knot sequences h = (1) h × Xi (2) h is called shape-regular, if the ratio of the diameter and the smallest edge of every rectangle defined by the distinct knots is uniformly bounded.
The spaceX h,k is spanned by all bivariate NURBSN h,k,i with knots h,k ⊂ [0, 1] 2 . We mention that the vector (ω h,k,i ) of weights in the definition (3) of the NURBS basis functions is obtained from the initial weights (ω k,i ) by "knot insertion"; i.e., the denominatorŵ remains the same before and after the refinement. As a common definition in IGA methods, the space X h,k is the pushforward Our assumptions about the knots guarantee that X h,k ⊂ C( k ), so that X h ⊂ X . Following Brivadis et al. [12], for each interface γ l we define the trace spaces and their counterparts onγ l W h,l = span(N h,l, j ; 1 ≤ j ≤ n h,l ),Ŵ h,l,0 = span(N h,l, j ; 2 ≤ j ≤ n h,l − 1).
Here, the univariate NURBS basis functionsN h,l, j are defined onγ l . The knot sequence is h,l = (r ) h,s(l) , r = 1 or 2, in the geometrically conforming case and in case NC1 of Remark 2.1, whereas h,l is only the part (r ) h,s(l) ∩γ l and the multiplicity of both endpoints is raised to q l +1. We denote by n = n h,l the dimension ofŴ h,l . Note that all linear combinations y = n j=1 c jNh,l, j satisfy the endpoint interpolation conditions y(ξ l,1 ) = c 1 and y(ξ l,2 ) = c n . This explains that W h,l,0 can be defined using the reduced index set 2 ≤ j ≤ n − 1.

The space M h defined by B-splines
The definition of the space M h of Lagrange multipliers is based on the space of polynomial splines of degree q l , Before we turn to the definitions on the physical domain, we provide some background about B-splines on the parameter intervalγ l and about the local spline projector. We drop the indices h, l and mention explicitly, in which way our results depend on the degree q and the knot vector . The fundamental stability result for B-splines [13] states that there is a constant κ > 0, which only depends on the degree, such that for arbitrary coefficients c j ∈ R the following inequality holds: κ is called the condition number of B-splines of degree q. Another well-known result for B-splines is the recurrence relation for derivatives Combined with (15), the Bernstein-type inequality follows easily for allv ∈Ŝ q ( ), where h min = min i (θ i+q − θ i ) > 0 is assumed. An important projector onto the spline space is the local spline projector from [21, Section 4.6]ˆ : where the linear functionals σ j only take values of f in the support ofB q j . Its "local" approximation properties are given next.
the so-called support extension of I i . Then for all 0 ≤ s ≤ q + 1 and f ∈ H s (γ ) , we have and if is locally quasi-uniform of order q, as in Definition 4.1, then The constant C 0 depends only on q, the constant C 1 depends on q and the bounds of local quasi-uniformity in (14).

Remark 4.3
The second result in [4,Lemma 4.3] is for semi-norms of order 1 ≤ r ≤ s and for locally quasi-uniform knot sequences, which means that adjacent non-empty knot intervals have comparable lengths. The proof in [4] uses a Bernstein inequality for splines. The slightly weaker assumption of "local quasi-uniformity of order q" is sufficient for the error estimate in the H 1 -semi-norm.
We are now ready to describe in more detail the choice of two subspacesM t h,l ⊂ S q l ( h,l ), whose pushforward will be our Lagrange multiplier spaces of dimension n − 2, with upper index t = 0 or t = 1 denoting the chosen alternative. Choosing subspaces of co-dimension 2 is typical in order to achieve stable infsup-conditions, as will be explained in detail.
The definition for t = 0 was proposed by Brivadis et al. [12, page 305] and describes the " p/ p setting with boundary modification," where p/ p stands for choosing the same degree forM h,l and the trace space W h,l . As before, we can drop indices h, l without causing confusion and let be an open knot vector onγ = [0, 1]. The spacê contains all splines whose polynomial pieces in the first and last knot intervals have degree at most q − 1. (This is also a common choice for Lagrange multipliers in FEM discretizations.) We summarize the results in [12]: If n ≥ q + 2, then (a)M 0 has dimension n − 2 and a basis with local support ; in particular, ρ j = 0 for j ≥ q + 2 and σ j = 0 for j ≤ n − q − 1, (b)M 0 contains all polynomials of degree q − 1.
The authors in [12] explain that the discrete infsup-condition where c does not depend on the h-discretization, has not been proved yet. They provide numerical evidence for its validity. We next define an alternative spaceM 1 ⊂Ŝ q ( ), for which we can prove the infsup-condition analytically for arbitrary knot vectors, i.e., without any assumptions on quasi-uniformity. It is defined as the orthogonal complement of a 2-dimensional subspaceÊ For the definition of the spaceÊ, we use two "boundary" B-splinesB 2q 1 andB 2q n−q of double degree, which are defined with respect to the given knot sequence ; in particular, the first knot ofB 2q 1 and the last knot ofB 2q n−q have multiplicity q + 1. Therefore, the function value and derivatives up to order q − 1 ofB 2q 1 at ξ = 0 vanish, and the same holds forB 2q n−q at ξ = 1. The basis ofÊ iŝ where the coefficients α j , β j can be computed recursively via (16). (We take a closer look at these coefficients in the proof of Lemma 5.3.) We obtain the following properties of the orthogonal complementM 1 .
where ρ j = 0 for j ≥ 2q + 2 and σ j = 0 for j ≤ n − 2q − 1. Proof Since we assume n ≥ q + 2, there is at least one interior knot and the functionŝ e 1 ,ê 2 satisfy the boundary conditionŝ This proves their linear independence. The basis in (20) is obtained by subtracting the orthogonal projection onÊ from every B-splineB q j , 2 ≤ j ≤ n − 1. Note that the basis functionsê 1 ,ê 2 have support Hence, we obtainμ j =B q j for all 2q + 2 ≤ j ≤ n − 2q − 1, and the remaining basis elements have support in [0, θ 3q+2 ] or [θ n−2q , 1], respectively. This is the result in (a). Being q-th derivatives of functions with compact support, bothê 1 andê 2 are orthogonal to all polynomials of degree q − 1. This implies result (b).
Finally, we can define the subspace of Lagrange multipliers on the physical domain. In order to gain a positive effect on sparsity of the mass matrix, as will be explained in Sect. 8, we use a modified pushforward of the spline spaceM 1 Here,ŵ l is the denominator of the univariate NURBS basis functions onγ l , which is the same for all h-refinements. This gives (9). We follow the standard approach to a priori error estimates (11) as in [24], and adapt the proofs to the IGA setting. We will develop our results under the following assumptions on the h-refinement. Assumption 2 (a) All refinement knot sequences h,k are quasi-uniform of order p k and shape-regular with constants independent of h. (b) If θ ∈γ l is an interior knot for the initial parameterization F l of γ l , and if its multiplicity is κ ≥ 2 (= multiple knot), then the same knot θ appears at least with multiplicity κ + 1 in the refinement knot sequence h,s(l) of the slave cell associated with γ l .

Roadmap to a priori error estimates
The second assumption (b) is new and was not observed to be required in [12]. Its importance will become clear in Theorem 7.4, where we prove error bounds for quasi-interpolation of the normal derivative ∂u k ∂n along γ l . Multiple knots in the initial parameterization F l can have the negative effect, that the normal derivative is not an element of the bent Sobolev space H q l −1 (γ l , h,l ). But this property is needed for error bounds of the consistency error. We will show, however, that this defect is repaired by Assumption 2(b).
The following theorems will be proved in Sect. 7.
in the geometrically conforming case and for type (NC2) of Remark 2.1. Moreover, we obtain for type (NC1) where L 0 ⊂ {1, . . . , L} denotes the subset of indices such that γ l is a proper subset of an edge of m(l) . The constants C, C 1 , C 2 do not depend on the h-refinement.
Theorem 4.6 With the same assumptions as above we have in the geometrically conforming case and type (NC2). The same adaptation as in (25) is needed for type (NC1). The constant C does not depend on the h-refinement.
Before we prove both theorems, we consider the infsup-condition with constant c independent of h, in Sect. 6. The method of proof is well established, see e.g. [9]. It makes use of the uniform L 2 -stability of the mortar projections defined in [5] for all interfaces γ l , Here ρ l is the positive weight function in (6). We treat the mortar projections in Sect.

5.
For most of the results in the following sections, there is some technical overhead caused by the IGA setting and the transition between the physical domain and the parametric domain. This transition is not always trivial, as the necessity of the condition in Assumption 2(b) indicates.

Mortar projection
We investigate the mortar projections (27) Our first observation is a simple connection ofQ with the orthogonal projection P : L 2 (0, 1) →Ŝ q ( ). Lemma 5.1 Letê 1 ,ê 2 be the functions in (19) and assume n ≥ q + 2, i.e. has at least one interior knot. Then we have for all f ∈ L 2 (0, 1) and the orthogonality remains valid if any combination ofê 1 andê 2 is added toP f . By the property (21) of the boundary values ofê 1 ,ê 2 , the function on the right-hand side of (29) has boundary values 0, and thus is an element ofŜ q 0 ( ). Clearly, it is the unique element ofŜ q 0 ( ) with orthogonality in (28).
The proof of the L 2 -boundedness ofQ is much more intricate, as it requires a detailed stability analysis of the basis functionsê 1 ,ê 2 in Lemma 5.3. As a profit, we obtain a precise upper bound which only depends on the degree and does not require any quasi-uniformity of the knots. We recall the definition of κ in (15).
Proof Based on (29), it is sufficient to show that Note that the function values ofP f at both endpoints are also coefficients in the B-spline representation The general stability result for B-splines in (15) implies So it remains to prove that holds for arbitrary constants c 1 , c n ∈ R. This is provided by the second inquality in (31) in the following Lemma.
where κ is the condition number of B-splines of degree q in (15).

Proof
Since is an open knot vector (with respect to the given degree), we have q+1 . We first look atê 1 in (19). We want to find upper bounds for the absolute values |α j | in We use an inductive argument and define intermediate derivatives The recurrence relation (16) starts from one coefficient α (0) 1 = 1 and gives Within this range of indices for j, we always have θ j = 0. Moreover, the right-hand side is the sum of two terms of equal sign (−1) j−1 . The inequalities are proved by induction: they are obvious for j = 1 and any r , and follow for 2 ≤ j ≤ r + 1 from For the last identity, we used the recursive relation (32) with j = 1 and α (r −1) 0 = 0. For r = q, we obtain the coefficients α j ofê 1 and have proved in (33) that where the last identity is well-known in combinatorics. The stability result in (15) is applied toê The same upper bound is obtained for the norm ofê 2 β n √ h n . Therefore, we obtain by Minkowski's inequality and this provides the upper bound in (31). Furthermore, the lower bound in (31) comes for free: by our assumption n ≥ q + 2 and (21), the first and last coefficients of the spline c 1ê1 , respectively. Hence, the stability result (15) gives the lower bound in the lemma.
The mortar projection Q h,l : L 2 (γ l ) → W h,l,0 on the physical domain was defined in (27) and can be described as follows. The open knot vector = h,l inγ l is used in the definition (28) ofQ =Q h,l . Then we have in (23), we obtain from (5) Hence, the operator norm in (30) is only changed by the influence of the initial parameterization F l of γ l . More precisely, if c denotes the constant on the right-hand side of (30), we have and likewise Hence, we obtain the following result. (27) have uniformly bounded operator norm on L 2 (γ l )

Theorem 5.4 The mortar projections Q h,l in
where c is the constant on the right-hand side of (30).

Remark 5.5
The estimate of Theorem 5.4 immediately extends to the mesh-dependent norm (12) Moreover, if the knot sequence h,l is quasi-uniform of order p s(l) , we also obtain the H 1 0 -stability of Q h,l by a similar argument as given in Lemma 1.3 in [24]. Indeed, for a functionf ∈ H 1 0 (γ l ) we can find an approximantv ∈Ŝ q l 0 ( h,l ) with homogeneous boundary conditions such that where C does not depend onf . (This is property (Sc) in [24, p.12].) By the projection propertyQ h,lv =v and a Bernstein inequality (here quasi-uniformity of order p s(l) of the knots is needed), we obtain Furthermore, by a typical interpolation argument, we obtain that A further generalization of this result to locally quasi-uniform knot sequences is not known to us.

Discrete infsup-condition
In this section, we prove the uniform infsup-inequality for the h-refined discrete spaces X h and M 1 h . The method of proof is well known, see e.g. [9], and mainly based on Theorem 5.4. We include some details as far as explicit constants are concerned. The infsup-inequality for the spline spaceŜ q 0 ( ) and the spaceM 1 in Proposition 4.4 on the parameter domain is given with respect to L 2 -norms.

Theorem 6.1 Let be an open knot vector on
where c is the constant on the right-hand side of (30).
Proof The proof follows [9]. Letμ ∈M 1 . By duality, the definition ofQ and the bound in (30), we have This gives the result in (35). As mentioned before, an analogous result for the Lagrange multiplier spaceM 0 (" p/ p setting with boundary modification") in [12] is not known. (ii) Since the dimensions ofM 1 andŜ q 0 ( ) agree, the order of the spaces in the infsupinequality can be switched, i.e.
In fact, the best possible lower bound on the right-hand side of (35) and (36) is the smallest singular value of the "mixed Gramian" ( φ j , ψ k ) j,k=1,...,n−2 with respect to L 2 -orthonormal bases of both spaces. For more details see [14].
In order to obtain the infsup-condition on the physical interfaces γ l , we follow the arguments in [24,Lemma 1.9] with the necessary changes caused by the parameterization. We use the mesh-dependent norm μ −1/2,h in (12) of elements μ ∈ M 1 h and the bilinear form b ρ in (6) with weight functions ρ l which depend only on the parameterizations F s(l) .

Theorem 6.3 Let Assumption 2(a) be satisfied. Then
where the constant C only depends on the initial parameterization.
h . By duality and by writing the weighted inner product φ, μ l ρ l as a standard product ρ l φ, μ l in L 2 (γ l ), we have The definition of the mortar projection in (27) and With the maximizing element φ l ∈ W h,l,0 , which is normalized by h Next we let φ l be the extension by zero to all of ∂ s(l) . Lemma 5.1 in [8] provides us with a function v l ∈ X h , which is zero on all patches k with k = s(l) and which is an extension of φ l to s(l) , such that .
The constant C 2 depends only on the geometry of the patches k , not on their discretization. By Assumption 2(a) and a Bernstein inequality for φ l , there is a constant Since Let r l denote the number of times the patch s(l) is counted as a slave domain. Then the Cauchy-Schwarz inequality and (39) and (40) imply Combined with (38) and (41), this gives This allows us to choose c = (

A priori error estimate
We let (u, λ) ∈ X × M be the solution of (7) and (u h , λ h ) ∈ X h × M 1 h be the solution of (9). Recall from (8) that λ| γ l is the weighted flux The restrictions to k and γ l will be denoted by u k , u h,k , and λ l , λ h,l , respectively. We will assume from now on, that the matrix-valued function α in (1) is for 1 ≤ k ≤ K ; this means that all partial derivatives up to order q − 1 exist and are Lipschitz continuous.
The standard approach to estimates of the a-priori error uses the coercivity of a(v, v) on X h and the identities (7), (9), in order to arrive at (10). Then the approximation error and the consistency error lead to the estimate (cf. [10]) Our results in Sects. 7.2-7.3 provide bounds for E a and E b and yield the proof of Theorem 4.5. The proof of Theorem 4.6 follows in Sect. 7.4. Because our setting differs from [24, Section 1.2], we include a step-by-step description, although some of the arguments are similar. We start with some new results on the approximation order by splines on physical interfaces in Sect. 7.1.

Spline approximation on interfaces
In this part, we provide some new results on the approximation order by splines on physical interfaces γ l . Some difficulties arise from the lack of a "standard" trace theorem for Sobolev spaces on domains with non-smooth boundary. We use the notion of C k,1 -curves for k-times differentiable regular curves with Lipschitz continuous kth derivative. We use results from the book [19] about Sobolev spaces H s (γ ) with non-integer s and non-smooth curves γ . To keep the notations simple, we consider a bounded domain ⊂ R 2 and a boundary curve γ ⊂ ∂ , which is a regular C 1,1 -curve and which has a NURBS parameterization F : The following notations appear frequently: the denominator of all NURBS basis functions isŵ, the length of the tangent is τ = |F |, the unit normal on γ is ν and the weight in (6) is ρ =ŵ 2 τ . Note that non-smooth interfaces γ are typical for the IGA mortar method, as soon as B-splines of degree q with interior knots are used for the parameterization F (rather than just Bernstein polynomials for a polynomial curve). Moreover, knot sequences with multiple knots occur in practical examples as the output of standard CAD surface models (e.g. Rhino, CGAL). If the open knot vector contains an interior knot θ i of multiplicity κ i ≥ 1, then the parameterization F has local smoothness C q−κ i near θ i , and γ is a C k,1 -curve near F(θ i ) with k = q − κ i . This has the following effect on the trace of u ∈ H q+1 ( ): we do not obtain u| γ ∈ H q+1/2 as for smooth curves, but only u| γ ∈ H k min +1 with k min = min i (q − κ i ), as stated in [19,Remark 4.3,p. 85].
One step towards the resolution of this defect was already described by Bazilevs et al. [3], and used in [12]. They introduce the notion of bent Sobolev spaces on the parameter domain. By the result in [4, Lemma 4.1], the (univariate) bent Sobolev space of integer order 0 ≤ s ≤ q + 1 is Its norm is the broken Sobolev norm of order s with respect to the knot intervals of .
. The spline f is constructed in order to "swallow up" all jumps of derivatives D k f , 1 ≤ k ≤ s − 1, at the knots θ i . Our goal in this subsection are optimal approximation rates of the local spline projectorˆ in (18), applied to the pullback of u ∈ H q+1 ( ) and its normal derivative on γ . As a first step, we specify the order of the bent Sobolev spaces for both pullbacks. It is at this point, where multiple knots of play a special role. We define an augmented knot vector + by increasing the multiplicity of θ i by one, if θ i is an interior knot of multiplicity κ i ≥ 2. Later on, Assumption 2(b) on the h-refinement requires that h should be a refinement of + .
In a similar way, we have η ∈ H q−1/2 (I i ) for all i. Here we note that the unit normal ν and the weight ρ on γ are piecewise analytic functions. The smoothness across θ i , however, is reduced by one, because the tangent and normal vectors of γ reduce the smoothness of η. We obtain We conclude that η ∈ H q−1 ((0, 1), + ) is satisfied for the augmented knot vector + .

Remark 7.2
Note that there is no need to enlarge , when has only simple knots, which means that γ is a C q−1,1 -curve. Therefore, the defect of the smoothness of the pullback η with regard to H q−1 ((0, 1), ) could not be observed in the numerical examples of [12], where knot sequences with simple knots were used throughout. In Sect. 9 and Appendix A we provide some numerical examples to demonstrate this defect.
The method in [3,4] If, in addition, h is a refinement of + , then  25] to non-integer s for both pullbacks v and η: there appears an extra factor h 1/2 in the error estimates of the next theorem as compared to the results in [4]. This will provide the optimal order for the approximation error E a in Sect. 7.2 and the consistency error E b in Sect. 7.3. For the corresponding error estimates on the physical interface γ , we define the push-forward of the local spline projector as in [4,Eq. (3.5)] Note that F provides the connection between the parameter domain and the physical interface, whereas the factorŵ is responsible for switching between NURBS basis functions and B-splines. The following error estimates are obtained by performing an analysis of super-convergence similar to the techniques in [23]. We present its proof in Appendix A.

Theorem 7.4 Let h be a quasi-uniform refinement of and h
With the same assumptions as in Proposition and where D denotes the tangential derivative.
If, in addition, h is a refinement of + , then and The constants C do not depend on h.
We mention without proof that the same estimates as in Theorem 7.4 are valid for the modified local spline projector h : H q+1 (0, 1) →Ŝ q ( h ) in [4, Proposition 2.3] which matches the boundary values at ξ = 0 and ξ = 1. This operator will be used in the proof of Theorem 7.6.
A further extension of (49) will provide the main ingredient for finding the upper bound of the consistency error in Sect. 7 as stated in (22), where θ h,n denotes the largest interior knot of h . For a quasiuniform refinement, and if h = max i (θ h,i+q −θ h,i ) is small enough, these intervals are disjoint and, more importantly, contained in the first/last knot interval of the initial knot sequence . Therefore, F| J h,i , i = 1, 2, is analytic. This is an important observation in order to obtain the following result.
Proof Let η h be the orthogonal projection of η =ŵ(λ • F) onto the spline spacê S q ( h ). By (49), we have In order to obtain the orthogonal projection of η ontoM 1 h , we subtract the function 1) dξ.
Here we used the assumption that the supports ofê h,1 andê h,l are disjoint. Orthogonality ofê h,i to all polynomials of degree q − 1 and the Cauchy-Schwarz inequality give The assumption on h allows us to use η ∈ H q−1/2 (J h,i ), by the standard trace theorem. This gives and ψ L 2 (0,1) ≤ |d 1 | + |d 2 | ≤ Ch q−1/2 u H q+1 ( ) .

Approximation error
The method in [24, Section 1.2] shows how to find an optimal bound for the approximation error E a . We consider the geometrically conforming case in Theorem 7.6 and provide the extension to the non-conforming cases (NC1) and (NC2) without proof in Remark 7.7.
Theorem 7.6 Let Assumptions 1 and 2(a) be satisfied, let h k be the meshsize in (13), and assume that every interface γ l is a full edge of s(l) and m(l) . Further assume that the weak solution u of (2) satisfies u ∈ H 1 ( ) and u k = u| k ∈ H p k +1 ( k ) for all 1 ≤ k ≤ K , where p k is the degree of the NURBS parameterization of k . Then where the constant C does not depend on the h-refinement.

Proof
The method described in [24, Lemma 1.4] can be used almost verbally. The constant C in the following estimates may change from step to step, but remains independent of h.
where C does not depend on h k . Moreover, for each interface γ l , 1 ≤ l ≤ L, and k = s(l), the tensor-product approach implies that w h,k | γ l is a local approximant of u k | γ l and interpolates u k at both endpoints of γ l , so (u k −w h,k )| γ l ∈ H 1 0 (γ l ). Moreover, by Theorem 7.4 we obtain the bound on the physical interface

Standard interpolation theory of Banach spaces provides the error estimate in the H
Since we restrict ourselves to the geometrically conforming case here, we also obtain The assumptions on the weak solution u imply that the jump [u] l is zero. Nevertheless, the jump can be non-zero, due to different refinements of the NURBS parameterizations F k and F m(l) . By the interpolation conditions at both endpoints of γ l , we still have [w h ] l ∈ H 1 0 (γ l ) (here, the geometrically conforming case is assumed) and Next we describe the construction of the approximant v h ∈ V h . We use the mortar projection Q h,l and define φ l = Q h,l ([w h ] l ) ∈ W h,l,0 . Let φ l be the extension of φ l to ∂ k by zero. Then we choose v l ∈ X h , which is zero in all patches j = k and satisfies where C only depends on the geometry of k . We obtain from the quasi-uniformity, the H

Finally, we define the function
It is an element of V h , because for every 1 ≤ l ≤ L, so that the definitions (6) and (27) give b ρ (v h , μ) = 0 for every μ ∈ M 1 h . For each k let L k be the set of indices 1 ≤ l ≤ L with s(l) = k . Then we obtain the error bound Note that the cardinality of L k is bounded by a constant which depends only on the global geometry; in our case of 2-D tensor-product patches we have L k ≤ 4. Therefore, summing the squared norms leads to the result.

Remark 7.7
The geometrically non-conforming scenario NC1 of Remark 2.1 is treated as follows. Let L 0 ⊂ {1, . . . , L} be the subset of indices such that γ l is a proper subset of an edge of m(l) . Then we cannot assume, as done in the proof of Theorem 7.6, that w h,m(l) interpolates u at both endpoints of γ l . Hence, the jump [w h ] l in (53) is not in H 1 0 (γ l ). However, the method in [24, p.18] describes a way how we obtain the upper bound The norm of the last term can be reduced by considering the restriction of u to a strip l ⊂ m(l) adjacent to γ l and of width h m(l) . For more details we refer to [24]. Things are much easier in the geometrically non-conforming scenario NC2. Recall from Remark 2.1 that the elements of X h,k need only be continuous across the prolongations of T-intersections from adjacent patches. Therefore, in the first step of the proof of Theorem 7.6, both components w h,s(l) and w h,m(l) can be chosen to interpolate u at both endpoints of γ l . Then we obtain the same result as for the geometrically conforming case.

Consistency error
The upper bound for the consistency error E b is obtained as in [24, Lemma 1.8].
Theorem 7.8 Let Assumptions 1 and 2 and the assumptions in Proposition 7.1 and 7.5 be satisfied for all cells k and interfaces γ l . Further assume that the weak solution u of (2) satisfies u ∈ H 1 ( ) and u k = u| k ∈ H p k +1 ( k ) for all 1 ≤ k ≤ K . Then we have where the constant does not depend on the h-refinement.
Proof For each interface γ l , we define the pullback The Cauchy-Schwarz inequality gives For the first sum in (57), we apply the result of Proposition 7.5 and obtain A bound for the second sum in (57) is obtained in analogy to [11,Lemma 3.5

] as
So we have obtained the upper bound in (56).
We can now combine the upper bounds for the approximation error in (52) or (55) and the consistency error in (56). Hence, by (43) we have proved the result of Theorem 4.5.

Error bound for the Lagrange multiplier
We next prove Theorem 4.6 for the error bound of the Lagrange multiplier It is derived from the infsup-condition and the approximation error in the standard way as shown in [24, p. 26]. We only describe the geometrically conforming case. For the nonconforming situation, the same adaptation as in (55) is needed.

Proof of Theorem 4.6
The same method as in [24, p.25] The result follows from Theorem 4.5 and Proposition 7.5.

Implementation and stability of mass matrices
For our implementation of the mortar finite element method, we follow [5, p. 194] and use the mass matrix M h,l for the substitution of all nodal coefficients of s(l) , which are associated with degrees of freedom in the interior of the interface γ l . We show that the local basis for the Lagrange multiplier space M 1 h,l is uniformly stable with respect to the h-refinement; this means that the condition number for this substitution is uniformly bounded.
Let 1 ≤ l ≤ L, k = s(l), q = p s(l) and n = n h,l be the parameters associated with the interface γ l . We often drop the indices h, l without further notice. First, we develop the expressions for the entries of the mass matrix which are associated with the basis elements ofM 1 h,l in Proposition 4.4. The rows of the mass matrix are indexed by 2 ≤ i ≤ n − 1 according to the enumeration of the basis elements ofM 1 h,l . There are two blocks of columns in the mass matrix. The first block is associated with finite elements on k which are non-zero on γ l , Here, w j is the weight factor in the numerator of the NURBS functionN j =N h,l, j . Both functionsμ 1 i ,B q j are splines inŜ q ( h,l ), so there are fast methods for the computation of these integrals. The analysis in Appendix B shows the following. This property is the same as the "spectral equivalence" in [24, p.13]. Note that this block is banded due to the local support of the basis.
The second block is associated with finite elements on the master cell m(l) which are non-zero on γ l . Note that the parameter intervalγ = F −1 m(l) (γ l ) can be a proper subset of (0, 1) in the geometrically non-conforming case. We denote the univariate NURBS, which are associated with the master cell and are non-zero onγ l , byN m j , 1 ≤ j ≤ñ. Note that a sign factor is introduced by the jump operation. In both geometrically conforming and non-conforming cases, we havẽ In a conforming geometry, we often have F −1 m(l) • F l = id, so this term can be omitted, but it must be included in the integral for a non-conforming geometry. The function w l is the denominator of the initial NURBS basis functions of the parameterization.
As an alternative basis ofM 1 h,l in our current implementation, we use basis functions (η i : 2 ≤ i ≤ n − 1) ofM 1 h,l such that the corresponding entries in the first block are instead of (58), and in the second block instead of (59). The basis will be chosen such that the square submatrix is a diagonal matrix. For this definition, we use the dual B-splines which satisfy the biorthogonality relation Here δ i, j denotes the Kronecker delta. The support of the dual B-splines is the whole intervalγ l . Although the new basis will have global supportγ l , the mass matrix is sparse due to the biorthogonality. and note that all coefficients α j with j > q + 1 and β j with j < n − q are zero.
In particular, the block M η h,l in (61) is diagonal. Proof Part (a) follows from two simple observations. First, the specified splinesη i , 2 ≤ i ≤ n −1, are linearly independent, because B q i only appears in the representation ofη i . Secondly, their orthogonality toê 1 andê 2 is evident. Part (b) follows directly from the biorthogonality of B-splines and dual B-splines.
Based on this result, we can implement the substitution method as follows. We let u h = (u h,k ) 1≤k≤K be an element of X h and fix 1 ≤ l ≤ L. Recall the notation for the jump • We use U for the coefficients associated with the slave cell and U m for the coefficients associated with the master cell, so • Fix 2 ≤ i ≤ n − 1. When we use the biorthogonality (62) and the representation (63) in the weak continuity condition, we obtain Hence, the substitution can be used for all coefficients U i , 2 ≤ i ≤ n − 1.

Remark 8.3 (i)
Note that the nodal coefficients U 1 and U n associated with the two endpoints of γ l in the slave-cell are treated as free parameters. Indeed, as in [5] and for the geometrically conforming case, we define separate control points for each cell k meeting at a vertex of the decomposition. Hence, different interfaces are completely decoupled and there is no interference of the control points on γ l with other cells than s(l) and m(l) . The same decoupling remains valid for the geometrically non-conforming types NC1 and NC2 in Remark 2.1. This is clear for NC1, because both endpoints of γ l are vertices of s(l) and can be treated in the same way as in the geometrically conforming case. For the type NC2, the control points U 1 and U n in (64) ( h,l ). In [16] we described a modification of the mortar method which avoids matrix inversion by the use of approximate dual B-splines. It introduces a second consistency error for the error analysis in Sect. 7. The error analysis will be presented in our forthcoming work.

Numerical examples
We present our numerical results for two examples. In Sect. 9.1 we consider the Poisson equation on the unit square, and Sect. 9.1.1 considers a benchmark problem of linear elasticity, namely an elastic plate with hole. In both cases, the numerical results can be compared to an analytical solution. Thus, the convergence behavior of the proposed mortar formulation can be assessed properly. Several discretizations are tested for each example. All computations are performed using an in-house isogeometric analysis code within Matlab.

Poisson equation solved on the unit square
In this example the Poisson equation − u = f is solved on the unit square = (0, 1) 2 . The manufactured analytical solution for the body load f (x, y) = 13 sin(3x) sin(2y) is u(x, y) = sin(3x) sin(2y). The associated boundary conditions are g = ∇u · ν on Neumann boundaries N with the outer normal vector ν, and u = 0 on the Dirichlet boundary D . Within this work the lower edge (x ∈ (0, 1), y = 0) is chosen as Dirichlet boundary, whereas all other edges are Neumann boundaries. Other choices are possible, but do not yield significantly different results, see the work of Zou et al. [25] where different combinations of boundary conditions are compared. In the following, we study several different possibilities to discretize the domain with NURBS patches. The accuracy of the computations is assessed with the help of the L 2 -error u − u h 0, which is plotted over the maximal element diagonal h. Since f is infinitely smooth, optimal error bounds for discretizations of degree p have the order O(h p+1 ) according to the general theory of finite elements [10].

Two conforming patches with straight interface
The first discretization scheme serves as a reference and uses a decomposition into two rectangular patches which intersect at x = 0.5. The discrete spaces X h,k , k = 1, 2, are tensor-product splines relative to conforming, equally spaced knot sequences and equal degrees 2 ≤ p ≤ 5 in both patches. The canonical parameterizations F 1 (ξ 1 , ξ 2 ) = (ξ 1 /2, ξ 2 ) and F 2 (ξ 1 , ξ 2 ) = ((1 + ξ 1 )/2, ξ 2 ) and the constant NURBS denominator w k = 1, k = 1, 2, lead to a constant weight function ρ = 1/2 for the bilinear form b ρ . Figure 2a shows the resulting error, if strong continuitiy conditions across the interface γ are used; here the two patches are coupled by shared degrees of freedom along γ and no consistency error occurs. Equivalently, the space V h consists of all tensor-product splines of coordinate degree p on the full unit square with Dirichlet boundary conditions for y = 0, whose knot vector (1)  Both proposed methods yield the same accuracy level and the same expected convergence rates as the conforming reference computations (Fig. 2a). Thus, the global error of u h is not affected by using the mortar method instead of a direct connection by shared degrees of freedom.

Two non-conforming patches with curved interface with internal C 1 -continuity
Our next example demonstrates the importance of Assumption 2(b), if an interface has limited smoothness. We choose two NURBS patches with a curved interface with initial degree p ini = 2, knot sequence = {0, 0, 0, 0.5, 1, 1, 1} and control points as listed in Table 1. The point (x, y) = (0.5, 0.5) corresponds to the parameter ξ = 0.5 on the interface, where the curve has only C 1 -smoothness. A non-conforming discretization with an element ratio of 2 : 3 along the interface is obtained by choosing this ratio for the coarsest mesh and then performing uniform subdivision in both patches. The same degrees are chosen in both patches. The maximal element diameter h is alike in both patches; see Fig. 3.
In order to demonstrate the effect of the limited smoothness of γ on the a priori error in Sec. 7, we first perform computations for degrees 2 ≤ p ≤ 5 without reduced    5, 1). We demonstrate in Fig. 4 that the orthogonal projection of ∇u · ν onto the spline space of degree p has an L 2 -error of size O(h 3/2 ), regardless of the degree 2 ≤ p ≤ 5. Therefore, the order of approximation in Proposition 7.5 is not achieved for 3 ≤ p ≤ 5. This results in a defect of the consistency error E b and, finally, in a degraded a-priori error as can be seen in Figs. 5 and 6. The effect for p = 3 is somehow damped for the L 2 -error in Fig. 5, but clearly visible for the L ∞ -error in Fig. 6. On the other hand, when we increase the multiplicty of ζ = 0.5 to p as described in Assumption 2(b) (at least for p ≥ 3), then the optimal convergence rate p + 1 for the a-priori error is recovered for all considered degrees and both cases t = 0 and 1, see Figs. 7 and 8. This clearly shows that the reduction of smoothness at interior knots is needed in order to recover the optimal convergence rates.   Fig. 9 Poisson equation: Discretization scheme with curved interface with initial order p = 3 and different internal continuities (left). Coarsest non-conforming mesh for this scheme with control points for the parameterization of degree 3 (right). The points at the interface with limited internal continuity are labeled by C 2 , C 1 and C 0

Two non-conforming patches with higher order curved interface with different continuities
Our next discretization scheme defines two NURBS patches with a curved interface γ with initial degree p = 3, knot sequence = {0, 0, 0, 0, 0.3, 0.3, 0.5, 0.5, 0.5, 0.7, 1, 1, 1, 1} and control points as listed in Table 2. At the points ξ = 0.3 and ξ = 0.7, the interface curve is only C 1 and C 2 -continuous, respectively. At the point ξ = 0.5, the curve is C 0 -continuous, and has a kink as can be seen in Fig. 9. The points on the interface with reduced internal continuity are marked on the right side of Fig. 9.
With this example we demonstrate the importance of Assumption 2(b) for different orders of smoothness and also include the treatment of points with C 0 -continuity as explained in Remark 2.2(b). A non-conforming discretization with an element ratio of 2 : 3 along the interface is obtained by choosing a coarse decomposition with this ratio and subsequent uniform subdivision, see Fig. 9. The same degree 3 ≤ p ≤ 5 is chosen for both patches. The error of u − u h 0, is shown in Fig. 11a, b for both spaces of Lagrange multipliers (t = 0 or 1) and for different choices of the reduction of smoothness at interior knots of . The mortar method which employs reduced smoothness in the definition of the discrete spaces X h,k according to Remark 2.2 and Assumption 2(b) shows optimal convergence (solid lines in 11(a) and (b)). For comparison, in computations labeled by NR we perform no reduction of smoothness; i.e. the knot vectors h related to the interface have knots 0.3, 0.5, 0.7 with multiplicities p − 1, p, p − 2, respectively. The convergence rate is not optimal for all considered degrees 3 ≤ p ≤ 5. It seems to be limited to O(h 1.5 ), which follows the mathematical reasoning in Appendix A. An intermediate case for the reduction of smoothness is labeled N1, if the C 0 -point at ξ = 0.5 is treated by the method in Remark 2.2, but no reduction of smoothness is performed at 0.3, 0.7; then the convergence rate is optimal only for p = 3, but deficient for higher degrees. Another case is labeled N2, if in addition to N1 we use reduced smoothness at ξ = 0.3, but not at ξ = 0.7; then the optimal convergence rate is obtained and the results cannot be distinguished by the eye from the solid lines. It seems that the defect of the consistency error shown in Fig. 10 is damped, but we have no analytical explanation yet. These results give further numerical evidence for the a priori error estimate in Sect. 7 and at the same time indicate that it might be sufficient to reduce the continuity in C 0 and C 1 -points. However, a reduction at all initial internal knots is not counterproductive. For a better understanding of the size of the consistency error E b in (42), Fig. 10 shows the L 2 -error of the orthogonal projection of ∇u · ν onto the full spline spacê S p ( h ) for 3 ≤ p ≤ 5. The labels NR, N1, and N2 indicate the partial reduction of smoothness as explained before. No label is used for the solid line in Fig. 10b, which shows the results where the reduction of smoothness at both knots 0.3 and 0.7 is performed according to Assumption 2(b) and the point ξ = 0.5 is treated as in Remark 2.2.

Discretization with ten patches and different kinds of intersections
Our final discretization scheme is with ten non-conforming patches, where seven Tintersections and two star-intersections are present. A sketch of the patch layout is given on the left side in Fig. 12. The interfaces are straight lines and the geometry is waterproof, i.e., the parametrizations match along the interfaces. In the convergence studies, every patch is uniformly subdivided by the same refinement factor and equal degrees are used. The presence of T-intersections, where one patch boundary borders two other patches, together with the additional knots in the patches with T-intersections yields non-conforming meshes at most interfaces. In the initial mesh given on the right of Fig. 12, no interior knots are present besides the additional knots at T-intersections,   which are resolved by the method described as NC2 in Remark 2.1. Three sample meshes are given in Fig. 13 for a better depiction of the obtained non-conformity. From a mathematical point of view there is no criterion for an optimal choice of patches to be master or slave. However, experience in numerical simulations in [12,16] suggests to use the patch with more elements as slave patch. The numerical results for this case are given in Fig. 14. We obtain optimal convergence rates for all degrees  Fig. 14a, b. This shows that the proposed formulation is able to handle multiple interfaces with different kinds of intersections properly. We also performed the same computation for an inverted master/slave classification with essentially no changes.

Linear elasticity solved on an elastic plate with hole
In this example the differential equations for linear elasticity are solved for an infinite plate with hole, where uniaxial tension is applied in x → ±∞. The mechanical equilibrium on a domain ⊂ R 2 is given by div σ + f = 0 with the boundary conditions u| D = 0 and σ · n| N = g, where n is the outer normal vector with respect to N . In the chosen planar linear elastic context with isotropic constitutive law, the relations between stresses σ , strains ε and displacements u are given by σ = λ tr (ε) I + 2με and where λ and μ are the Lamé parameters, I is the identity matrix and ∇ is the gradient operator. Accordingly, the saddle point problem (7) changes to: where X and M are vector valued extensions of the spaces used in Eq. (7). The bilinear forms a(u, v) and b ρ (u, μ) are given by respectively. For our computations, we limit the domain to finite size and apply the tractions of the exact solution, which can be found e.g. in [2], as Neumann boundary conditions on the relevant edges. Furthermore, we consider only one quarter due to symmetry and apply the associated symmetry boundary conditions. Thus, we consider a domain = (x, y) ∈ (−4, 0) × (0, 4) : x 2 + y 2 ≥ 1 , see Fig. 15. The known analytical solution of this problem allows performing convergence studies for a complex stress distribution and is thus commonly used for numerical studies, especially in the framework of isogeometric analysis.
In the following, we study three different possibilities to discretize the domain with NURBS patches. In all cases the geometry is modeled exactly. The accuracy of the computations is assessed with the help of the L 2 -error norm σ − σ h 0, which is plotted over the maximal element diagonal h. According to the theory of finite elements [10], the slope in the double logarithmic diagram should be p, which is indicated by the dashed lines in the same diagram. In order to allow a straightforward comparison of the mortar results to the reference case, we use the same axes in all figures for this example. The reference example for a conforming method of degrees 2 ≤ p ≤ 5 is provided in Fig. 17a.

Two patches with curved interface with internal C 1 -continuity and conforming meshes
The first decomposition is with two NURBS patches with a curved interface with initial order p ini = 2, knot sequence = {0, 0, 0, 0.5, 1, 1, 1} and control points as listed in Table 3. At the point ξ = 0.5, the interface γ is only C 1 -continuous. The mesh  is chosen conforming in both patches. A sketch of this discretization scheme along with the coarsest computed mesh is given in Fig. 16. Finer meshes for the convergence analysis are obtained by uniform subdivision. This discretization is used to assess the ability of the mortar method to capture the same approximation rate as conforming finite element methods. A computation where the two patches are coupled by shared degrees of freedom along the interface, is used as reference. These results are labeled by conf. in Fig. 17a. The numerical results for the mortar method with Lagrange multiplier space M t h,l , t = 0 or 1, are given in Fig. 17b, c. While both methods yield a lower accuracy level for coarse discretizations than the conforming method, they seem to catch up in terms of accuracy for fine meshes as almost the same error level as in the reference computations is obtained for small stepsizes h. Thus, the expected convergence rates are obtained, even though the error level is higher for coarse meshes. The mortar method does not improve the error as compared to conforming FEM, but it is clearly competitive. We mention that the use of the mortar method does not primarily lie in computations with conforming meshes, after all.
Another way of comparing the conforming FEM and the proposed mortar methods is obtained by plotting the logarithmic error log( σ v M − σ h v M / σ v M ) of the von Mises stresses σ v M = σ 2 11 + σ 2 22 − σ 11 σ 22 + 3σ 2 12 for a fixed discretization. We use a subdivision factor of 20 and degree p = 5 in Fig. 18. In the stress error plot of the reference computations (Fig. 18a) the largest error is inside the domain, while for both proposed methods (Fig. 18b, c) the largest errors are at the end points of the interface. This peculiarity can be explained as follows: In the conforming case, the number of control points along the interface is equal on both sides of the interface, i.e. n h,l =ñ h,l . By using the endpoint modification for M t h,l , the control points at both endpoints of the slave patch become free parameters. This leads to a situation, where the number n h,l − 2 of slave control points is smaller than the number of master control points, which can result in a substantial increase of the consistency error near both endpoints. As a simple work-around, we subdivide the first and last two rows of elements in patch 2 . By doing so, the mesh turns into a non-conforming mesh, where the number of interface slave control points is larger than the number of interface master control points. This reduces the concentrated error at the endpoints of the interface significantly. A sample mesh is given in Fig. 19a. The corresponding stress error plots are given in Fig. 19b, c.

Two patches with curved interface with internal C 1 -continuity and non-conforming mesh
The second discretization scheme uses the same initial discretization as in Sect. 9.2.1, but now non-conforming meshes are chosen. This is obtained by using a subdivision factor of 2 j + 1 along the interface in patch 2 and j in patch 1 . This yields a ratio of 2 j + 1 : j elements along the interface. In the second parametric direction, the subdivision factor is chosen in a complementary way in order to obtain a similar number of elements in both patches. The patch 2 with the smaller stepsize along the interface is chosen as the slave patch. The obtained meshes for a subdivision factor of j = 5 and j = 10 of this discretization scheme are given in Fig. 20. This discretization is used to assess the ability of the proposed methods to use discretizations with curved NURBS interfaces with limited internal smoothness and strongly non-conforming knot sequences. The error norms of the global stress distri- bution computed using the Lagrange multiplier spaces M t h,l , t = 0, 1, are given in Fig. 21a, b. Both methods yield optimal convergence rates and the error levels are comparable to the conforming case using shared degrees of freedom given in Fig. 17a.
Besides accuracy, also efficiency of the methods is studied. This is done by comparing computational costs for the individual stages of the computations. Furthermore, the stability of the methods is assessed by means of the condition number of the global stiffness matrix. The results of Sect. 9.2.1 for the coupling by shared degrees of freedom are used as reference. These results are labeled by conf.
The computational cost for the formation of the global stiffness matrix is compared in Fig. 22, where the results are shown in CPU seconds on a contemporary dual core notebook with 8 GB of RAM. The peaks in the diagrams are due to background activity of the operating system. The entries are computed by optimal integration as proposed in [18] and directly assembled into the sparse matrix format of Matlab, whereby vectorized assembly is used. The computational cost for conforming meshes given in Fig. 22a is lower than for the mortar methods (Fig. 22b, c), but they range in the same order of magnitude. The difference occurs in the assembly process, when the entries of the slave interface control points are assembled to master interface degrees of freedom according to (65). The excessive growth of computational costs in the fine limit is attributed to limitations of RAM, which was restricted to 8 GB for this study. The computational cost for the formation of the mass matrices as explained in Sect. 8 are given in Fig. 23 for both proposed mortar methods with M t h,l , t = 0, 1. There is no significant difference between both methods. In the conforming case, this cost is saved since a direct connection by shared degrees of freedom is used. Note that the cost is in the same order of magnitude as the cost for the formation of the global stiffness matrix, besides the fact that the coupling matrices require the computation of a line integral only, whereas the global stiffness matrix requires computations of a global surface integral. However, there is some potential for a speed-up of the computation of the entriesm η l,i, j of the mass matrices in (60). First, the computation of the line integrals in (60) uses the time-consuming iterative point inversion algorithm for the mapping from the slave to the master patch, which is required in every integration point. This routine could be written as an external routine in C or Fortran. Secondly, we did not yet implement a fast method for the inversion of the Gram matrix of the B-splines on the slave patch in order to compute the dual B-splines in (63), see also Remark 8.3.
The computational cost for the solution of the global system of equations is given in Fig. 24. There are no significant differences between the computations using shared degrees of freedom (Fig. 24a) and both proposed mortar methods (Fig. 24b, c). The computational cost for the solution grows almost linearly, as can be expected for the used sparse matrix format. It is roughly two orders of magnitude smaller than the cost for the formation of the stiffness matrix.
The influence of the proposed mortar methods on the sparsity of the global stiffness matrix is depicted in Fig. 25, whereby the number of non-zero entries (nz) of the stiffness matrix is given below the diagrams. The sparsity pattern of computations using about 3200 elements are compared between coupling by shared degrees of freedom and the mortar method. The upper left block represents the stiffness matrix  Fig. 25c is about 0.5 % smaller than for t = 0 in Fig. 25b). It can be observed in the very last rows and columns of the stiffness matrix. The nonzero entries for t = 0 result from the interrelation between all interface control points U i in the slave patch with both control points U 1 and U n at the endpoints of the interface, see (65). This interrelation is reduced to only few interface control points near the endpoints for t = 1, because most coefficients α i and β i in (65) are zero. This advantage is only small for 2D-problems, but will be more pronounced for 3Dproblems. Moreover, in our planned extension of the mortar method by the use of an h-dependent bilinear form b h,ρ in our future work, we reduce all blocks of the stiffness matrix to banded form by the application of locally supported "approximate duals" instead of the dual B-splines.
The condition number of the global stiffness matrix is compared in Fig. 26 between computations using shared degrees of freedom and the proposed mortar methods. The results show that the condition number is not perceivably affected by using the proposed mortar methods, which is an important requirement for robust and accurate computations.

Discretization with ten patches and different kinds of intersections
In order to assess the ability of our method to deal with different types of geometry, we use a discretization scheme with ten non-conforming patches with four T-intersections and three star-intersections. A sketch of the patch layout is given on the left side in Fig. 27. The interfaces are straight and have infinite internal continuity, so Assumptions can be neglected. The geometry is waterproof, as the parametrizations along the interfaces are matching. On the right side in Fig. 27 the coarsest initial mesh is drawn, where the prolongation of all T-intersections as C 0 -continuous lines is already included (type NC2 in Remark 2.1). In the convergence studies, every patch is refined using a number of a · j + b elements, where the values of a and b are given in Fig. 28a for each parametric direction within each patch. The refinement is performed in a way that the lengths of the element spans are as similar as possible in the knot vector under consideration of the prescribed element boundaries which arise due to the prolongation of ending interfaces at T-intersections. The order of the basis functions is chosen uniformly within the whole domain. Three sample meshes are given in Fig. 28 for a better depiction of the obtained non-conformity. The error norm of the global stress distributions is given in Fig. 29. Both proposed methods (Fig. 29a, b) yield optimal convergence rates. The error levels are slightly higher than in the conforming case (Fig. 17a), where shared degrees of freedom are used. It is to be noted that in the conforming case the difference between the individual

Fig. 28
Elastic plate with hole: Sample meshes for the discretization scheme with ten non-conforming patches and different kinds of intersections. The applied refinement rule is given in a). All finer meshes are obtained by choosing a factor j ∈ N element diameters h is smaller than in the case with ten patches, and thus naturally a lower error level is produced. The computational costs for the formation of the global stiffness matrix ranges in the same order of magnitude as for the conforming case. The difference occurs in the assembly process and grows with the ratio of interface degrees of freedom to domain degrees of freedom, see Sect. 9.2.2. The same statements as in Sect. 9.2.2 can be made about the computational cost for the formation of the mass matrix and for the global solution.
The influence of the proposed mortar methods on the sparsity of the global stiffness matrix is studied in Fig. 30. The sparsity pattern of computations using 3141 elements are compared between both proposed mortar methods. The number of non-zero entries for t = 1 is about 0.8 % smaller than for t = 0. It can be observed that the proposed approach for t = 1 creates fewer interrelations between patches than for t = 0. Furthermore, the interrelation within patches is less pronounced, see the banded structure of the block in the center of the diagrams: For t = 0, there is an interrelation between the control points at both ends of the interface (visualized by the square around the banded structure). For t = 1, there is no such interrelation. The condition number of the global stiffness matrix is compared in Fig. 31 between computations using shared degrees of freedom and the proposed mortar methods. Apart from very coarse meshes, the behavior of the condition number is very similar, both in magnitude and in slope. No negative impact of the proposed coupling method on the condition number can be detected. This shows that robust and accurate computations are possible with the proposed mortar method.
The difficult part of the proof is the split w = w 1 + w 2 , with global smoothness of w 1 ∈ H 1/2 (0, 1) and local smoothness of w 2,i := w 2 | I i ∈ H 1 (I i ), Then the imbedding H 1 (I i ) ⊂ C(I i ) allows us to define the jumps The spline has degree q and simple knots, so it is in H q (0, 1). It is defined to carry all discontinuities of w 2 , such that we obtain global smoothness D q v 2 − w 2 ∈ H 1 (0, 1) by local smoothness in H 1 (I i ) and continuity across ζ i . We now provide more details on this splitting of w = D q (v − v 1 ). Consider the expansion of D q (v − v 1 ) by means of Leibniz' and Faa di Bruno's formulas. By collecting all terms with partial derivatives of u of order q, we obtain w 1 :=ŵ where we let F 1 , F 2 denote the coordinate functions of F : [0, 1] → R 2 . By Assumption 1(b), all factorsŵ(F 1 ) r 1 (F 2 ) r 2 are in C 0,1 [0, 1]. Combined with D (r 1 ,r 2 ) u ∈ H 1 ( ), the trace theorem in [19,Theorem 5.5,p. 95] gives w 1 ∈ H 1/2 (0, 1), |w 1 | H 1/2 (0,1) ≤ C r 1 +r 2 =q |(D (r 1 ,r 2 ) u) • F| H 1/2 (0,1) All other terms of D q (v − v 1 ) contain partial derivatives of u of order less than or equal to q − 1 and derivatives ofŵ, F 1 , F 2 up to order q. Therefore, by u ∈ H q+1 ( ) and analyticity ofŵ and F on every I i , we have Consequently, the jumps d i in (66) are bounded by The numbers d i define the spline v 2 in (67) which is an element of H q (0, 1).
Finally, the last piece in our splitting is the function In order to show global smoothness v 3 ∈ H q+1/2 (0, 1), we observe that The first term w 1 is in H 1/2 (0, 1) and its H 1/2 -seminorm has the bound in (69). The second part w 2 − D q v 2 has global smoothness H 1 (0, 1) as explained in (67). Its H 1/2 -seminorm is bounded by Because D q v 2 is constant on every interval I i , its seminorm vanishes. By (70) and the (local) imbedding H 1 (I i ) → H 1/2 (I i ), we obtain (The constant C from (70) has grown by a factor N .) Therefore, we have shown the global smoothness v 3 ∈ H q+1/2 (0, 1) and obtain the upper bound for its seminorm |v 3 | H q+1/2 (0,1) ≤ C u H q+1 ( ) .
Now we obtain the error estimate (47) for any quasi-uniform refinement h of . Since v 1 + v 2 ∈Ŝ q ( ) can be ignored, the local spline projector in Proposition 4.2 provides the approximation 1) ≤ Ch q+1/2− j u H q+1 ( ) for j = 0 and j = 1.
(b) A similar proof can be given for in (49). We describe the required adaptations and leave further details aside. Especially the last step of the error estimate forˆ h requires more attention, because we use a spline η 2 outside of the spaceŜ q ( + ), which is not reproduced byˆ h .
We start with the same steps as before. The trace theorem provides the local smoothness η ∈ H q−1/2 (I i ) for every 1 ≤ i ≤ N − 1. The projector in (44) is now chosen to define a spline η 1 ∈Ŝ q ( + ), such that η − η 1 ∈ H q−1 (0, 1) and D q−1 η 1 | I i = 0 for all i.
This completes the proof of (49). The upper bounds in (48) and (50) are directly obtained from the definition of h in (46). Indeed, for the L 2 -bound in (48), we use and a similar expression if u| γ is replaced by λ in (50). For the tangential derivative in (48) we use and obtain the desired upper bound from (47).

Appendix B
We prove the uniform stability of the block M h,l of the mass matrix in Proposition 8.1. For this purpose, we let be an arbitrary open knot vector in [0, 1] and n be the dimension ofŜ q ( ). We define the normalized basis functions Note that M h,l is obtained from M 0 by diagonal scaling. For a quasi-uniform knot vector h,l , the condition numbers of both matrices are comparable. Because the matrices are banded, with bandwidth independent of , local quasi-uniformity is sufficient to make the condition numbers of M h,l and M 0 comparable.