Construction of approximate $C^1$ bases for isogeometric analysis on two-patch domains

In this paper, we develop and study approximately smooth basis constructions for isogeometric analysis over two-patch domains. One key element of isogeometric analysis is that it allows high order smoothness within one patch. However, for representing complex geometries, a multi-patch construction is needed. In this case, a $C^0$-smooth basis is easy to obtain, whereas $C^1$-smooth isogeometric functions require a special construction. Such spaces are of interest when solving numerically fourth-order PDE problems, such as the biharmonic equation and the Kirchhoff-Love plate or shell formulation, using an isogeometric Galerkin method. With the construction of so-called analysis-suitable $G^1$ (in short, AS-$G^1$) parametrizations, as introduced in (Collin, Sangalli, Takacs; CAGD, 2016), it is possible to construct $C^1$ isogeometric spaces which possess optimal approximation properties. These geometries need to satisfy certain constraints along the interfaces and additionally require that the regularity $r$ and degree $p$ of the underlying spline space satisfy $1 \leq r \leq p-2$. The problem is that most complex geometries are not AS-$G^1$ geometries. Therefore, we define basis functions for isogeometric spaces by enforcing approximate $C^1$ conditions following the basis construction from (Kapl, Sangalli, Takacs; CAGD, 2017). For this reason, the defined function spaces are not exactly $C^1$ but only approximately. We study the convergence behavior and define function spaces that converge optimally under $h$-refinement, by locally introducing functions of higher polynomial degree and lower regularity. The convergence rate is optimal in several numerical tests performed on domains with non-trivial interfaces. While an extension to more general multi-patch domains is possible, we restrict ourselves to the two-patch case and focus on the construction over a single interface.


Introduction
Isogeometric Analysis (IGA), which is introduced in [16], is a method for numerical simulation combining Finite Elemente Analysis (FEA) with Computer Aided Design (CAD). Within the IGA framework, the same spline functions are used for the exact representation of the CAD geometry and for the approximation of the FEA solution. IGA offers many advantages over classical (piecewise linear) finite elements by providing a basis of high smoothness and high polynomial degree. It is therefore ideal for solving high order partial differential equations (PDEs) over geometries comprised of a single patch. However, most geometries of interest are not given as a single patch, but are represented by a collection of patches forming a so-called multi-patch domain. Note that, in general, CAD models are composed of trimmed patches, cf. [27], which we do not address here.
In this paper, we assume that the geometry is represented by a two-patch parametrization where the patch parametrizations are matching along the interface. On such C 0 -matching, two-or multi-patch domain, one can construct a C 0 -smooth basis in a rather straighforward way, see e.g. [4,35]. Enforcing higher order smoothness over multi-patch domains is however non-trivial, except in regular regions, as in [7,33]. As a consequence, standard basis constructions, which are only C 0 -smooth over patch interfaces, cannot be used directly for solving high-order PDEs. In the following, we focus on fourth order problems, such as the biharmonic equation or a Kirchhoff-Love plate or shell formulation. There are several different methods to overcome the problem of reduced smoothness.
One way is to use a C 0 multi-patch basis and to solve the fourth order problem using Nitsche's method. This approach is studied e.g. in [2,13,32]. Due to the reduced regularity, additional integral terms are derived over all interfaces and a penalty term is introduced to the problem statement. In this way, the C 1 -error, i.e., the jump of the normal derivative across the interface, is penalized. Thus, using Nitsche's method, one has to derive a more complicated variational formulation, depending on the underlying PDE and discretization space, and one has to find a good choice for the penalty parameter, which is also not always straightforward. Another approach is using the mortar method, see e.g. [14] or [6] for C 0 -coupling. The mortar method is based on constraint minimization where the coupling constraints are enforced using Lagrange multipliers. The correct choice of the discrete Lagrange multiplier space, leading to a stable formulation, is non-trivial. Moreover, the mortar method results in a saddle-point problem of larger size than the original problem.
A different possibility to solve forth order problems over multi-patch domains is to perform strong C 1coupling, where the basis functions are coupled strongly across the interfaces, thus creating a C 1 -smooth space over the multi-patch domain. The first work, which is following the idea of strong C 1 -coupling, is the so-called bending strip method, see [24,25]. The idea was later generalized to construct C 1 bases over multi-patch domains as in [10,17,22,23,28,29]. See also [15] for a summary of related approaches.
We follow the constructions in [18,20], which are based upon the findings in [10,22], where an explicit formula for constructing a C 1 basis at the interface is stated. As developed in [12], the C 1 continuity of an isogeometric function is equivalent to the G 1 geometric continuity of its graph surface. This geometric continuity, cf. [30], yields so-called gluing data for each interface from which one can construct a C 1 basis. However, within the isogeometric framework, this construction is only possible for analysis-suitable G 1 (in short AS-G 1 ) geometries which were characterized in [10]. AS-G 1 is defined by having linear gluing data for each interface. This class of geometries contains for instance bilinear patches. However, for most geometries the gluing data is not linear.
Hence, all approaches based on strong C 1 -coupling across interfaces have a significant problem: they can only be applied to certain geometries. If the geometry is not AS-G 1 , one may locally raise the polynomial degree or reduce the continuity requirements to obtain a sufficiently large space. In [8,9] the authors follow the former strategy, by constructing a C 1 -smooth spline space at the interface using spline functions of a higher polynomial degree. A similar strategy is also proposed in [19]. In this article, we intend to follow the latter strategy by properly reducing the continuity requirements.
In Figure 1 we compare the approximation powers of different example parametrizations. While the geometry depicted in Subfigure 1a is AS-G 1 , the geometry in Subfigure 1b is not. Constructing the C 1smooth space for both geometries and solving the biharmonic equation, we observe the following behaviour as plotted in Subfigure 1c: the discretization using the C 1 -smooth space over the AS-G 1 geometry yields optimal convergence rates and the C 1 -smooth space over the non AS-G 1 geometry does not allow any convergence. This lack of convergence can be circumvented as follows: instead of using the gluing data, we introduce so-called approximated gluing data which is then used to construct the basis functions at the interface. Since the gluing data is now approximated, the resulting space is only approximately C 1 -smooth. However, with the correct choice of approximated gluing data, the optimal convergence rate is restored. This can be observed in Subfigure 1c. In this paper, we focus on two-patch domains and extend the construction of C 1 basis functions from AS-G 1 geometries to general two-patch geometries.
The outline of the paper is as follows. We start with the definition of the model problem, more precisely the biharmonic equation, in Section 2. In Section 3, we recall the definition of B-splines and introduce the notation we use. The description of the geometry mapping is given in Section 4. The main part of the basis construction is explained in Section 5. Here, we define the gluing data and the spaces which are used for solving the biharmonic problem. In Section 6 the properties of the approximate C 1 spaces are discussed.  . In (c) we compare the resulting convergence rates. While the convergence rate is optimal for the AS-G 1 geometry, the C 1 -smooth discretization over the non AS-G 1 geometry does not allow convergence. However, a properly constructed approximate C 1 -smooth discretization over the same geometry again yields optimal rates.
The discrete problem is stated in Section 7 which is used to obtain the numerical results shown in Section 8.

Model problem
In this paper, we focus on the biharmonic equation. Let Ω be a bounded open subset of R 2 with a sufficiently smooth (piecewise Lipschitz) boundary ∂Ω and a given source function f . We consider the fourth order problem with the boundary conditions ϕ = g 0 on ∂Ω and (2) where g 0 and g 1 are given. One can use different boundary conditions such as ϕ = g 0 and ∂ n ϕ = g 1 on ∂Ω, but in this paper we focus on the boundary conditions stated in (2) and (3). We assume that all functions f , g 0 and g 1 are sufficiently smooth, i.e., f ∈ H −2 (Ω), g 0 ∈ H −3/2 (∂Ω), g 1 ∈ H −1/2 (∂Ω). Note that (2) is enforced as an essential boundary condition, which can be eliminated. Hence, we assume from now on that the problem is homogeneous. Let The weak formulation of the problem (1)-(3) is the following.
where the bilinear form is defined as a(ϕ, ψ) := Ω ∆ϕ ∆ψ dx and the right hand side as where ∂ n is the normal derivative at the boundary.
We solve Problem 1 using an isogeometric discretization, cf. [16]. As it is common in IGA, we assume that the domain Ω is parametrized with B-spline patches. A discretization space can then be defined on Ω based on the same B-spline space as the geometry parametrization. Thus, we recall the definition of B-splines in the next section. A more detailed introduction to IGA can be found, e.g., in [4,11].

B-spline spaces
In this section, a brief overview of B-splines is stated. Given positive integers p, r and n, and an (uniform) mesh, with mesh size h = 1/n, we define the (univariate) spline space S(p, r, h) of degree p and regularity r, with r < p, as The open knot vector Ξ := (ξ 1 , ..., ξ N +p+1 ) with N = p + 1 + (p − r)(n − 1) satisfies We mention here that, for simplicity, the regularity is assumed to be the same at all interior knots and consequently all the knots have the same multiplicity. Given the knot vector Ξ and a polynomial degree p one can define the B-spline functions denoted as b i , 1 ≤ i ≤ N , using the Cox-de Boor recursion, see [31]. We have A conforming discretization of Problem 1 requires H 2 -regularity of the discretization space. To this end, we assume that the underlying spline space is in H 2 , which is equivalent to C 1 -smoothness.
The definitions can be extended to the two-dimensional case by means of a tensor-product structure. Let h = (h 1 , h 2 ) be the pair of (uniform) mesh-sizes and (Ξ 1 , Ξ 2 ) be the two knot vectors, one for each direction. Additionally, we define the bivariate B-spline functions as b i := b 1,i1 ⊗ b 2,i2 where any univariate B-spline function has the degree p s and the regularity r s , s ∈ {1, 2}. The tensor-product spline space S(p, r, h) is spanned by the bivariate B-spline functions, yielding where p = (p 1 , p 2 ) and r = (r 1 , r 2 ).

The two-patch geometry
We assume that the domain Ω is given as the union of two non-overlapping subdomains, i.e., we have open subdomains Ω (S) , for S ∈ {L, R}, such that with a single interface Γ which is defined as Here Ω (S) denotes the closure of Ω (S) . In this paper, we always consider the notation S ∈ {L, R} where L denotes the left patch and R the right patch. Furthermore, each Ω (S) is a spline patch with the geometry where S (S) = S(p (S) ,r (S) ,ĥ (S) ) is a tensor-product spline space as defined in Section 3 and Ω = [0, 1] 2 . We assume that the mappings are regular, i.e., Moreover, we assume that the patch interface is along an entire edge of both patches. Without loss of generality, on each patch the interface is parametrized by (u, v) ∈ {0} × (0, 1), which can be achieved by a simple reparametrization (a combination of translation, rotation and symmetry). Furthermore, we assume that the patch parametrizations agree along the interface, summarized in the following.
Assumption 2 (C 0 -conformity at the interface). The parametrizations of the two patches meet C 0 along the interface, i.e., For simplicity, we assume S As a consequence, the left and the right patch share the same tangential derivative along the interface. Hence, the unit tangent vector is given by where τ (v) = t(v) . We denote the outward pointing unit normal vector to Ω (S) by n (S) . Along the interface Γ, one can compute the normal vector n = n (L) = −n (R) , which satisfies the following proposition.
Proposition 1. Given the geometry mapping F (S) , the normal vector can be expressed as Proof. Since n ⊥ t 0 , the vector ∂ u F (S) can be uniquely described as a linear combination of n and t 0 , i.e., where λ and µ are the two unknown factors. Using the vector projection of ∂ u F (S) onto t 0 gives us the first unknown Since n and t 0 are unitary vectors, we have Then, the desired result, including the sign of µ, follows directly from the definition of n = n (L) = −n (R) . Figure 2 gives an overview of the domain setting. Figure 2: Example of the general setting for the two-patch parametrization.

The isogeometric discretization
In this section, we define isogeometric functions over two-patch domains and discuss their continuity conditions. As developed in [10,18,20,22], the class of analysis-suitable G 1 (in short AS-G 1 ) geometries allows optimal approximation. In [20], the Argyris isogeometric space A is introduced as the direct sum of the single patch-interior, edge and vertex components. Its name is derived from the fact that the vertex space is obtained from interpolating C 2 -data at every vertex, similar to the Argyris finite element. Moreover, the edge space can be split in degrees of freedom for function values as well as normal derivative (or general crossing derivative) values along the interface.
However, the construction for the space A is only possible for certain geometries, i.e., for AS-G 1 parametrizations, which are discussed in Section 6, Definition 2. For general geometries, a different approach to construct an (approximate) C 1 isogeometric space is introduced in this section and discussed in more detail in Section 6. Since only two-patch domains are considered, a slight modification of the structure of the space is performed: there is no need to define separate vertex spaces, thus the space A is split into the patch-interior spaces and the interface space containing all functions that have non-vanishing trace or crossing derivative at the entire interface, including the vertices.

Spline spaces of mixed regularity
In Definition 1 we introduce spline spaces of mixed regularity. Following the definition in (5), the uniform spline space S(p, r, h) has n = 1/h polynomial segments with the distinct inner knots {h, 2h, ..., (n − 1)h}.
Note that the notation S(p, (r,r), (h,ĥ)) simplifies to S(p, r, h) if the geometry space S(p,r,ĥ) has no inner knots, i.e.,ĥ = 1, or if r =r. In Figure 3 an example is depicted.
Corresponding to the spaces S (S) containing the geometry mappings, we define the discretization spaces as S (L) and S (R) , with For the sake of simplicity, we only consider discretizations where the spline spaces of both patches are matching at the interface. This leads to the following restriction on the discrete space.
Assumption 3 (Matching two-patch discretization). We assume that the discrete spline spaces are matching at the interface, i.e., we have S 2 := S(p 2 , (r 2 ,r 2 ), (h 2 ,ĥ 2 )) = S In the following we define the C 0 isogeometric space, followed by the definition of the C 1 isogeometric space. Section 5.4 introduces the gluing data which is needed for the construction of the approximate C 1 basis, which is described in Section 5.5.

The space of C 0 isogeometric functions
The space of C 0 isogeometric functions on Ω is given as Following standard FEM notation, we denote the discrete space with a subscript h. Here h represents the mesh size of the spline spaces S (S) , for S ∈ {L, R}. The mesh size of V 0 h in physical space is always of the same order as h. As the spline spaces are matching at the interface, one can easily construct a basis for the C 0 isogeometric space since for each function with non-vanishing trace on one side there exists exactly one function on the other side having the same trace.

The space of C 1 isogeometric functions
The space of C 1 isogeometric functions on Ω is given as One can describe the C 1 continuity of a function at the interface by studying the geometric continuity of its graph surface. The graph surface Σ ⊂ Ω × R of an isogeometric function ϕ : Ω → R consists of the two graph surface patches Considering only regularly parametrized patches, one can see that the C 1 continuity of an isogeometric function at the interface Γ is equivalent to the G 1 geometric continuity of its graph parametrization, i.e., there exists for each point at the interface a well-defined tangent plane to the graph surface. The tangent plane is well-defined if and only if the graph surfaces fullfill for all v ∈ [0, 1] This condition is known as G 1 (geometric) continuity, cf. [30]. Figure 4 illustrates the G 1 continuity. Proposition 2. An isogeometric function ϕ ∈ V 0 h belongs to V 1 h if and only if its graph surface Σ is geometrically continuous of order 1, in short G 1 , at the interface Γ.
For a discussion and generalizations of the equivalence above, see [10,12,22]. Figure 4: A visualization of G 1 continuity of the graph surface. To obtain G 1 continuity, the four vectors, two of them being equal, are coplanar (and span the tangent plane) for each point along the interface.

The C 1 condition across the interface Γ
If the graph surface is G 1 continuous, then there exist functions One can uniquely determine the functions α (L) , α (R) and β up to a common function γ : where t is defined as in (7). Furthermore, there exist non-unique functions is satisfied for all v ∈ [0, 1]. One possible choice for the functions is 8 cf. [10,Proposition 1]. The functions α (S) and β or more generally, the functions α (S) and β (S) are called gluing data. The first two lines of (10) are equivalent to Proposition 2 can be reformulated with the help of the gluing data.
, assuming γ(v) ≡ 1, the functions α (S) and β fulfill α (S) ∈ S(2p 2 ). The functions β (L) and β (R) are in general piecewise rational functions with regularityr In order to obtain an optimal convergence rate for the gluing data, see Proposition 5, we need the following smoothness condition for the gluing data.
Assumption 4. We assume that the gluing data satisfies Note thatr (S) 2 ≥ 2 is a sufficient condition for Assumption 4.

Construction of an approximate C 1 basis
A basis construction for AS-G 1 two-patch geometries was developed in [18]. In the following, we provide a variation of that approach, which extends the construction to general geometries by relaxing the smoothness condition. Instead of constructing the C 1 -smooth isogeometric space exactly, we define a basis of isogeometric functions which are only approximately C 1 -smooth. We call the resulting space V 1 h the approximate C 1 isogeometric space on Ω. It is defined as where the interface space A Γ is the space of functions which have non-vanishing traces or derivatives at the interface and A (S) • are the patch-interior spaces, which have support on Ω (S) and have vanishing value and normal derivative on Γ, hence, they satisfy A where A Γ,+ spans certain traces along the interface Γ and A Γ,− has vanishing trace and (approximately) spans certain normal derivatives along Γ. We use the notation · to signify that the spaces are only approximately C 1 and in general V 1 h ⊂ V 1 h . Details on the behaviour of the normal derivative across the interface are discussed in Section 6. In Figure 5 the construction of the spaces is illustrated.

The patch-interior spaces
The patch-interior spaces are spanned by those isogeometric basis functions which have vanishing function values and derivatives at the interface Γ, that is,  • . The eliminated basis functions are replaced with approximate C 1 basis functions which span the space A Γ , see (c). Here, the basis functions, which are represented with stars, span certain traces along the interface and those, which are visualized as diamonds, span certain (approximate) normal derivatives along the interface. The corresponding spaces are denoted with A Γ,+ and A Γ,− , respectively.

The approximated gluing data
In order to construct the space A Γ we introduce an approximation of the gluing data. The approximated gluing data is taken from the spline space S( p, ( r,r 2 −1), (h 2 ,ĥ 2 )), where we prescribe the polynomial degree p ≥ 2 and the regularity r ≥ 1 in advance. Then α (S) , β (S) are computed by a projection operator P h onto S( p, ( r,r 2 − 1), (h 2 ,ĥ 2 )) with where α (S) is defined in (11) and β (S) in (13). The functions α (S) and β (S) are called the approximated gluing data. They do not fulfill (10) exactly, but approximately. For a suitable projection operator, we have . If the gluing data satisfies α (S) , β (S) ∈ S( p, ( r,r 2 − 1), (h 2 ,ĥ 2 )), then (18) is exactly zero. One can also set p = r ≥ 0, if the gluing data is a polynomial function of degree p, i.e., α (S) , β (S) ∈ P p . In this case, the requirement p ≥ 2 can be dropped. Using the approximated gluing data we construct basis functions along the interface.

The interface space
We define a basis following the approach presented in [18, Section 5.2], where we replace the gluing data with the approximated gluing data. Let {b + j } 1≤j≤N+ be the basis for the spline space and let {b − j } 1≤j≤N− be the basis for the spline space where N ± are the dimensions of the corresponding spaces. The optimal choice for the degrees p + and p − and the regularities r + and r − will be discussed later. The interface space with approximate C 1 continuity is given as with We have by definition where r * 2 = min( r, r + − 1, r − ). Depending on p * 2 , r * 2 andr 2 , the interface space A Γ is clearly not necessarily a subspace of V 0 h and therefore does not yield an isoparametric discretization. Note that the isogeometric concept is violated by using a spline space of lower regularity and (in general) higher degree near the interface. A brief overview of the steps for constructing the interface spaces is shown in Figure 6. Figure 6: The steps for constructing the approximate C 1 space V 1 h : the dashed arrows illustrate the construction steps and the solid line corresponds to the projection step.

Optimal choice of the spline parameters
In order to be refineable spline spaces, the polynomial degrees and regularities of S + and S − need to satisfy The functions f (S) (j,±) as in (21) have to satisfy f (S) (j,±) ∈ C 1 ( Ω), in accordance with Assumption 1. Thus, the regularity r * 2 needs to satisfy r * 2 = min( r, r + − 1, r − ) ≥ 1. To be able to reproduce traces and (approximate) normal derivatives of optimal order, that is, of degree p 2 and p 2 − 1, respectively, the degrees for S + and S − need to satisfy p + = p 2 and p − = p 2 − 1, respectively. Note that a higher polynomial degree for S + and S − will not improve the global approximation properties.
To summarize, we obtain r ≥ 1 as well as From now on, we choose the maximum regularity for the spaces S + and S − in order to achieve the smallest number of degrees of freedom, i.e., r + = p 2 − 1 and r − = p 2 − 2. We conclude from these restrictions that the degree p 2 needs to satsify p 2 ≥ 3.
Assumption 5 (Minimum polynomial degree at the interface). We assume that the polynomial degree for the discrete space at the interface fulfills p 2 ≥ 3.
In Section 6, we discuss the role of the approximated gluing data. Theorem 1 and Conjecture 1 and the numerical experiments show how to choose the degree p and regularity r of the approximated gluing data in order to get optimal convergence rates. To obtain a sufficiently smooth spline approximation of the gluing data, we need the following.
Assumption 6 (Requirement on the approximation of the gluing data). If the gluing data are not polynomial functions, we assume that the approximated gluing data is computed from S( p, ( r,r 2 − 1), (h 2 ,ĥ 2 )), with 1 ≤ r ≤ p − 1.
If the gluing data are polynomials of low degree, they can be reproduced exactly using the space P p = S( p, ∞, 1), whereas if they are polynomials of high degree, they can be approximated using splines from S( p, p−1, h 2 ). Throughout the following section we assume that the gluing data are not polynomial functions.

Properties of the approximate C 1 space
In this section, the properties of the approximated gluing data as well of the approximate C 1 space V 1 h are studied. The properties of the projector P h , which is used to define the approximated gluing data, are described in Subsection 6.1. In Subsection 6.2, the boundedness of the jump of the normal derivative at the interface is proven. We can show that the convergence rate of this error depends only on the polynomial degree p of the approximated gluing data. In Subsection 6.3, we study the special case when the jump of the normal derivative vanishes and we introduced the AS-G 1 case. We conclude this section with a brief discussion of the two possible cases at the boundary, see Subsection 6.4.
We assume that 1 ≤ r ≤ p − 1 and that the operator satisfies the following properties: • it preserves splines, i.e., • it is L ∞ -stable, i.e., there exists a constant C > 0, such that • it interpolates at the boundary, i.e., and • it satisfies the following estimate: there exists a constant C > 0, such that for all 0 ≤ j ≤ 1/ĥ 2 − 1 we have whereω j = [jĥ 2 , (j + 1)ĥ 2 ] and s ( p+1) is the derivative of s of order p + 1.
For each subintervalω j there exists a local projector satisfying (23) and (25), cf. [34,Theorem 6.25]. A modification of that projector, interpolating function values at the global boundary and derivatives up to orderr 2 − 1 at all inner knots jĥ 2 , yields a global projector satisfying (22) and (24). Such a construction is similar to the one presented in [5,Proposition 3.2].
In practice, the required properties can be relaxed, since the projector is used only to prove a bound as in Proposition 5. Note that (24) is required to simplify the imposition of boundary conditions, as described in Section 6.4. The condition may be dropped, as discussed in Remark 2. A desirable property is that applying the operator to a piecewise rational function should be computationally cheap. In the following, the influence of the parameters p and r on the normal jump of the interface space is studied.

Estimating the jump of the normal derivative of the approximate C 1 space
Since the space V 1 h is not exactly C 1 -smooth but only approximately, we want to estimate the jump of the normal derivative across the interface. Let n be the normal vector to the interface, as in (9), and let ϕ at x ∈ Γ, defined as a limit. Similarly, the tangential derivative along the interface is expressed by Note that the tangential derivative is continuous and therefore the limit is well-defined and does not depend on the side. The jump of the normal derivative is defined as It satisfies the following bound.
Theorem 1. Let h 2 be small enough and the gluing data α (S) , β (S) ∈ Cr 2−1 ([0, 1]) . Then we have for all where C > 0 depends on the geometry, but not on the mesh size. For certain configurations the jump may vanish. This is characterized in Section 6.3.
Note that h 2 needs to be sufficiently small, such that for all v ∈ [0, 1] and for S ∈ {L, R}. Such a bound on h 2 always exists, since α (S) converges to α (S) pointwise and 1/α (S) is bounded from above. Before we can state the proof of Theorem 1, some preliminary estimates are needed. The jump of the normal derivative satisfies the following pointwise representation.

Proposition 4.
We have for ϕ h ∈ A Γ and for all x ∈ Γ, with x = F (L) (0, v) and v ∈ [0, 1], that Due to the symmetry of the construction, a similar statement is valid with switched sides.
Proof. Using the definition of ϕ h we have f . Furthermore, one can describe ∂ u F (S) (0, v) as a linear combination of t 0 and n, like in the proof of Proposition 1, and it follows By construction, cf. (20), f for some functions G 1 and G 2 , which are independent of the side S ∈ {L, R}. Hence, we get and and consequently independent of the side S ∈ {L, R}. We obtain In the following we replace ∂ n ϕ h − ∂ n ϕ h , we get a representation for the right patch. Using (26) we obtain Multiplying both sides with α (R) τ /α (R) yields which concludes the proof.
Due to the approximation using the projector P h , the following estimates for the factors E where the constant C depends on p, and on the geometry mappings F (L) and F (R) , but not on the mesh size.
Proof. Throughout the proof, all norms are to be considered L ∞ -norms. We have and similarly The terms τ 2 and 1 α (R) are bounded by definition and depend only on the geometry mapping F (R) . Due to (23) the term α (R) is bounded from above by C α (R) , which in turn depends only on p and on F (R) . Estimate (25) yields where the constants depend only on p and on F (S) . What remains to be shown is an estimate from above for 1/ α (L) . Due to the regularity of patch F (L) , we have which is bounded from above by a constant that depends only on Such an ε exists if h 2 is sufficently small. This concludes the proof.
Based on the projection P h which interpolates the boundary, see (24), it follows that α (S) (v) = α (S) (v) and β (S) (v) = β (S) (v) and consequently E (S) 1}. In addition, it may happen that for certain points v ∈ (0, 1) the approximated gluing data α (S) (v) and/or β (S) (v) match with the gluing data α (S) (v) and/or β (S) (v). In this case, the corresponding factor also vanishes at v. We can now proof Theorem 1.

Proof of Theorem 1. From Proposition 4 and 5 we obtain
where C depends only on p and on the geometry parametrizations, but not on h or ϕ h . For reasons of symmetry, the same bound is valid for ϕ (R) h . Using a standard Sobolev trace inequality, cf. [1], yields which concludes the proof.

A special case: vanishing jumps From Proposition 4 it can be concluded that the jump of the normal derivative is zero at
for some c = 0. Obviously, a sufficient condition is α (S) (v) = α (S) (v) and β (S) (v) = β (S) (v) for S ∈ {L, R}.
Proposition 6. Let α (S) and β (S) be in the spline space S( p, ( r,r), (h 2 ,ĥ 2 )) and (27) be satisfied. Then the interface space A Γ is C 1 -smooth and, consequently, Moreover, when the gluing data are linear polynomials, i.e., α (S) , β (S) ∈ P 1 , then the interface basis functions are by definition in the space S(p (21). Those geometries are also known as analysis-suitable G 1 geometries, as introduced in [10]. We repeat the definition here.
Note that the AS-G 1 condition requires the existence of linear gluing data. However, we define α (L) , α (R) through the formulas in (11), with γ ≡ 1, and β (L) , β (R) through (13). Hence, even though the gluing data we compute is not linear, there might exist linear gluing data for a different choice of γ or a different splitting of β = α (L) β (R) − α (R) β (L) . Bilinear patches are AS-G 1 and they yield linear gluing data for γ ≡ 1. Thus, we have the following for bilinear patches.
Proposition 7. If the geometry is AS-G 1 with γ ≡ 1, e.g., if both patches are bilinear, then V 1 h ⊆ V 1 h . Hence, when the parametrization is piecewise bilinear, then the isogeometric concept remains and the jump vanishes. It is shown numerically in [10,18,20], that the convergence rates for general AS-G 1 geometries are optimal when solving fourth order problems. Moreover, optimal approximation error bounds were proven in [21] for bilinear patches.

Functions with vanishing trace at the domain boundary
In this section, we characterize the space of interior functions V 1 h,0 = V 1 h ∩ H 1 0 (Ω) of functions with vanishing trace on the domain boundary. We denote its complement, used to impose non-homogeneous boundary conditions, by V 1 h,∂Ω . By definition, cf. in (20), one can see that the interface basis functions which are, in general, not vanishing at the boundary corresponding tov = 0 that is Since this equation needs to be satisfied for all u, we can cancel out the factor 1,2 (u) and get Consequently, Hence, we obtain a non-trivial solution (λ 1 , λ 2 , λ 3 ) = (0, 0, 0) if and only if The kernel is then given by (λ 1 , λ 2 , λ 3 ) = (0, c, −c ψ), c ∈ R. Since the projector P h in (17) interpolates at the boundary, we have β(0) = β(0). To summarize, we need to distinguish two cases: 1. If β(0) = 0, which is equivalent to By properly modifying the functions β (L) and β (R) , i.e., replacing β (S) by for S ∈ {L, R}, we achieve ψ = 0. This simplifies the definition of the kernel, which is then spanned by B * (2,+) = B (2,+) . Thus, if the boundary is smooth atv = 0, the function B * (2,+) belongs to V 1 h,0 , whereas the functions B (1,+) and B (1,−) belong to V 1 h,∂Ω . If the boundary is not smooth, all functions belong to V 1 h,∂Ω . A similar modification can be achieved if the boundary is smooth atv = 1, that is, if β(1) = 0, or if the boundary is smooth on both ends of the interface. In Figure 7 an example of the two cases is depicted. There, the boundary is smooth at the lower end of the interface and non-smooth at the upper end.
Remark 2. Note that the condition in (24), requiring the projector P h to interpolate at the boundary, may be dropped. However, in that case the approximate gluing data β for some constant C that depends on the degree and on the exact geometry, similar to Proposition 5. Thus, one has to compute an approximate kernel up to an h-dependent tolerance.

The two-patch formulation and discretization
In this section, the discrete variational problem is stated. We consider a sequence of approximate Note that the spaces V 1 h are in general not nested, only their underlying C 0 spaces. The behaviour of the approximate solution as the mesh size h goes to zero is studied numerically for different choices of polynomial degrees in Section 8.
7.2. The approximate C 1 isogeometric discretization We consider the discrete space V 1 h,0 . Under the assumptions summarized in Figure 8, we have V 1 h,0 ⊂ X 0 . We then solve the following discrete problem. where Furthermore, we have, by definition, Note that this discrete problem is not an exact discretization of Problem 2, since we omit the jump term, which vanishes in the limit by construction. The speed of convergence then depends on the bounds on the jump term as in Theorem 1, which depend on the approximation of the gluing data. We expect that the two-patch model problem with approximate C 1 -smoothness at the interface satisfies the following a-priori error estimate. Conjecture 1. Let the assumptions of Theorem 1 be satisfied and let p (S) 1 ≥ 2, p 2 ≥ 3 and r 1 , r 2 , r ≥ 1.
We set q = min(p ∈ P p and q = min(p , otherwise. Let ϕ be the solution of Problem 1, with ϕ ∈ H 5/2+ε (Ω) and ϕ| Ω (S) ∈ H q+2 (Ω (S) ), and let ϕ h be the solution of Problem 3, then we have All requirements for constructing the interface space to obtain an optimal convergence rate are summarized in Figure 8.

Numerical experiments
In the following we perform numerical experiments on four two-patch geometries. In those geometries the patches meet C 0 at the interface in accordance with Assumption 2. On each geometry we solve Problem 3 with the exact solution u(x, y) = (cos(4πx) − 1)(cos(4πy) − 1), using the approximate C 1 space V 1 h described in Section 5.5.
Let for simplicity p = p Two-patch geometry: Approx. gluing data: Figure 8: Summary of the requirements for the spaces. It starts with the given C 0 -matching, two-patch geometry where the patches meet C 0 at the interface. The gluing data needs to be at least C 1 . A sufficient condition for that isr 2 ≥ 2. Then the discrete space is constructed as in Definition 1. Furthermore, p 2 ≥ 3 is needed for constructing the interface space. For simplicity, we only assume matching spaces at the interface. Next, the space for the approximated gluing data is defined. Here, one needs p to be sufficiently large and r ≥ 1 to obtain an optimal convergence rate.
The results are reported in Subsection 8.2. In Subsection 8.3 we present the convergence rates of the error measured in L 2 -, H 1 -and H 2 -norms for various polynomial degrees p and p ≥ 3, see Assumption 5. The observed rates are consistent with Conjecture 1. In Subsection 8.4, we conclude the numerical tests with comparisons of varying spline regularity r. All tests are implemented within the open-source C++ library G+Smo, cf. [26].

The geometries for the numerical tests
The numerical tests are performed on four two-patch geometries shown in Subfigures 9a-9d. Example I is an AS-G 1 geometry as described in Definition 2, whereas Examples II, III and IV are non AS-G 1 geometries. The corresponding exact solutions are shown in Subfigures 9e-9h. The first two Examples I and II describe the same domain, but have different parametrizations. While Example I, which is composed of bilinear patches, has a straight interface, Example II is composed of bicubic patches and has a curved interface. Example III and IV are a geometries which both have a corner at one end of the interface. There the boundary space is modified as described in Subsection 6.4. Furthermore, Example III depicts a quarter of a plate with a circular hole, where the circular arc is approximated by a cubic B-spline curve. The Examples I, II and III are constructed with no internal knots, i.e. the geometries are constructed with Bézier patches. Therefore, one can setr 1 =r 2 = ∞ and the gluing data is C ∞ -smooth. In contrast, Example IV is a B-spline geometry with degreep 2 = 3 and regularityr 2 = 2 ergo the gluing data is only C 1 -smooth. For each geometry the gluing data is computed as defined in (11) and in (13). One can see that the gluing data for Example I is linear, see Subfigure 9i, while for Examples II-IV the gluing data is not linear, see Subfigures 9j-9l.

Convergence of the jump of the normal derivative
In this section, we provide convergence results for the jump of the normal derivative of the discrete solution at the interface, that is, we compute where ϕ h is the solution of Problem 3. In all four examples we fixed the polynomial degree to p = 3 and the regularity to r = 1. The approximate gluing data is computed for splines of degree p ∈ {1, 2, 3, 4} and with maximum regularity r = p − 1. The results are shown in Figure 10. Since the first geometry is AS-G 1 , the normal jump is (numerically) zero (see Subfigure 10a) and thus the discrete solution is C 1 -smooth at the interface, up to tolerance. In the other Examples II-IV the geometry is not AS-G 1 . The observed convergence rate p+1 of the normal jump is consistent with Theorem 1, as can be seen in Subfigures 10b, 10c and 10d.

Dependence of convergence rates on approximation of gluing data
We compare the convergence rates of the error measured in the L 2 -, H 1 -and H 2 -norms for varying polynomial degree p of the spline space and varying polynomial degree p of the approximate gluing data. We expect that the approximate gluing data must be of degree p ≥ p − 2, such that the convergence rates are optimal, as stated in Conjecture 1.
Therefore, we plot the L 2 -, H 1 -and H 2 -errors for Examples I to IV for polynomial degrees p ∈ {3, 4, 5}. The approximate gluing data is constructed with the lowest polynomial degree p = max(p−2, 2) and highest possible regularity r = p − 1 to achieve the smallest number of degrees of freedom for the approximate C 1 construction. The results are summarized in Figure 11.
For the AS-G 1 geometry in Example I, we obtain optimal convergence rates for all polynomial degrees p. The rate is independent of the degree of the approximate gluing data, which follows from the fact that the gluing data is linear and the normal jump vanishes for all p ≥ 1. Hence, according to Proposition 7 one obtains an exactly C 1 -smooth space. The errors are plotted in Subfigures 11a, 11e and 11i for degrees p = 3, p = 4 and p = 5, respectively.
While all errors converge optimally for Example I this is not the case for Examples II, III and IV, which we study in the following. One can see that for those non AS-G 1 geometries the rates are not optimal for any degree p if the gluing data is approximated with ( p, r) = (1, 0). The reason for this is that the space A Γ of interface functions is not conforming, as it is not C 1 in a neighborhood of the interface (see (21)). Hence, the space violates Assumption 1 which leads to non-optimal convergence rates.
Note that for p = 3 the convergence rates are optimal if the gluing data is approximated with p ≥ 2, see Subfigures 11b, 11c and 11d. Increasing the degree for the approximate gluing data does not reduce the errors significantly as the curves overlap.
A similar behaviour can be observed for polynomial degree p = 4, as can be seen in Subfigures 11f, 11g and 11h. Optimal rates are obtained for p ≥ 2.
For polynomial degree p = 5 approximating the gluing data with splines with p = 2 is not enough. One requires at least p = 3 to reach an optimal error rate as shown in Subfigures 11j, 11g and 11h. In Subfigure 11k, one cannot observe a difference between p = 2 and p = 3. A possible reason is that the normal jump for Example IV is significant smaller than the H 2 -error. We expect that as the mesh becomes more refined the consistency error from the jump of the normal derivative will dominate the H 2 -error and the convergence rate will deteriorate, similar to the behaviour for Examples II and IV.

Different regularity r
In the following we focus on different regularities and therefore on different numbers of degrees of freedom. For each polynomial degree p ∈ {3, 4, 5} the gluing data is approximated by splines with p = max(p − 2, 2) and r = p−1. The results are shown in Figure 12. One can see that in all examples the errors for r ≤ p−2 are quite similar with respect to the number of degrees of freedom, whereas the errors with r = p − 1 are almost the same with less than half of the degrees of freedom. This shows that using splines of maximum regularity, a comparable error can be achieved with significantly fewer degrees of freedom. Hence, the underlying linear system is much smaller and the computation time can be reduced significantly. Note that for the standard AS-G 1 basis construction as developed in [18] the regularity is bounded by r ≤ p − 2 globally to obtain the nested, isogeometric spaces V 1 h . For the construction we propose here, reduced regularity is only needed for the interface space A Γ , whereas the patch-interior spaces A (S) • , for S ∈ {L, R}, can be constructed with r = p − 1. As a consequence the approximately C 1 -smooth spaces V 1 h are not nested, even though the underlying C 0 -smooth spaces V 0 h are.

Conclusion and future work
In this paper, we constructed and studied approximately C 1 -smooth spaces V 1 h over general two-patch domains which can be used for solving fourth order problems. Following the approach in [18], the basis construction for the space is simple and combines standard patch-wise basis functions with specific interface functions. The construction in [18], which is based on [10,22], describes a basis for analysis-suitable G 1 parametrizations. In that case, the parametrizations need to satisfy the AS-G 1 condition which requires the gluing data to be linear functions. This is a severe restriction, as most (generic) spline parametrizations are not AS-G 1 . Thus, in general, a reparametrization as developed in [17] is necessary.
Our approach relaxes this AS-G 1 condition on the geometry and is applicable on most two-patch domains. The only requirement is, that the gluing data is C 1 , which is always satisfied if the patch parametrizations are at least C 2 . Generic spline parametrizations yield gluing data which is either piecewise polynomial of high degree or rational. Instead of reparametrizing the domain to obtain linear gluing data, we base our construction on an approximation of the given gluing data. Therefore, we only obtain approximate C 1 -smoothness at the interface. In other words, we get a jump of the normal derivative at the interface. The space V 1 h is given as the direct sum of the subspaces A (L) • , which are the patch-interior spaces, and A Γ , which is the interface space. The patch interior spaces are standard isogeometric spaces which

Mesh size h
Error (a) Ex. I with p = 3, r = 1.

Mesh size h
Error (j) Ex. II with p = 5, r = 1.                have vanishing function value and vanishing gradient at the interface. The interface space is composed of functions that span traces as well as functions that span (approximate) normal derivatives at the interface. Since the construction of the interface space is based on approximated, nonlinear gluing data, the space A Γ is of higher polynomial degree and lower regularity locally near the interface, as stated in (21). As a consequence the approximately C 1 -smooth spaces { V 1 h } h are not nested, even though the underlying family of C 0 -smooth spaces {V 0 h } h is refined by knot insertion and therefore nested. The advantage with using the approximated gluing data compared to [18] is that we allow geometries that are not necessarily AS-G 1 geometries. Furthermore, we show that by allowing non-nested spaces, splines of maximum regularity can be used away from the interface. In contrast, the standard construction over AS-G 1 parametrizations requires r ≤ p − 2.
In the future, we want to extend the construction of the approximately C 1 -smooth spaces to multi-patch domains. A possible approach is to introduce vertex spaces by interpolation, as in [20], thus enforcing C 2 super-smoothness at all vertices. Another challenge is to obtain a construction that results in nested spaces. This may be done by avoiding spline spaces of locally higher degree and reduced regularity. Moreover, a construction with uniform degree p everywhere allows the use of a standard Gaussian quadrature rule in all elements, whereas the construction proposed in this paper requires a quadrature rule of higher order in all elements neighboring the interface. Furthermore, due to the non-standard structure of the space near the interface, a complete numerical analysis of the proposed approach is beyond the scope of this paper. An extension of the construction to surfaces is also of practical relevance, since for many applications, a given surface geometry is only smooth up to some prescribed tolerance. Hence, the C 1 -smoothness of functions defined on such surfaces may also be imposed only approximately.