An approximate $C^1$ multi-patch space for isogeometric analysis with a comparison to Nitsche's method

We present an approximately $C^1$-smooth multi-patch spline construction which can be used in isogeometric analysis (IGA). A key property of IGA is that it is simple to achieve high order smoothness within a single patch. To represent more complex geometries one often uses a multi-patch construction. In this case, the global continuity for the basis functions is in general only $C^0$. Therefore, to obtain $C^1$-smooth isogeometric functions, a special construction for the basis is needed. Such spaces are of interest when solving numerically fourth-order problems, such as the biharmonic equation or Kirchhoff-Love plate/shell formulations, using an isogeometric Galerkin method. Isogeometric spaces that are globally $C^1$ over multi-patch domains can be constructed as in (Collin, Sangalli, Takacs; CAGD, 2016) and (Kapl, Sangalli, Takacs; CAGD, 2019). The constructions require so-called analysis-suitable $G^1$ parametrizations. To allow $C^1$ spaces over more general multi-patch parametrizations, we need to increase the polynomial degree and relax the $C^1$ conditions. We adopt the approximate $C^1$ construction for two-patch domains, as developed in (Weinm\"uller, Takacs; CMAME, 2021), and extend it to more general multi-patch domains. We employ the construction for a biharmonic model problem and compare the results with Nitsche's method. We compare both methods over complex multi-patch domains with non-trivial interfaces. The numerical tests indicate that the proposed construction converges optimally under $h$-refinement, comparable to the solution using Nitsche's method. In contrast to weakly imposing coupling conditions, the approximate $C^1$ construction is explicit and no additional terms need to be introduced to stabilize the method. Thus, the new proposed method can be used more easily as no parameters need to be estimated.


Introduction
Computer Aided Design (CAD) is used to create digital models of geometric objects.In a CAD model the object is usually described with the help of two-dimensional curves as well as three-dimensional surfaces and volumes.These free-form curves and surfaces/volumes can be described by means of splines, i.e., piecewise Email addresses: pascal.weinmueller@jku.at(Pascal Weinmüller), thomas.takacs@ricam.oeaw.ac.at (Thomas Takacs) Preprint submitted to Elsevier February 10, 2022 arXiv:2202.04516v1[math.NA] 9 Feb 2022 polynomial functions.Such CAD models can be used in many applications, e.g., for simulations based on the Finite Element Method (FEM) or Isogeometric Analysis (IGA), which is considered in this paper.IGA, as introduced in [16], uses the same spline functions that are used to construct the geometry also for the discretization spaces for the computation of numerical simulations.One advantage of IGA over classical higher-order FEM is that it provides basis functions with high smoothness and high polynomial degree, making it ideal for solving high order partial differential equations (PDEs) over single patches.In most applications, however, the geometries can usually not be described with one patch.Multi-patch domains composed of a collection of several patches, or related concepts, are needed.
In this paper, we restrict ourselves to planar multi-patch domains where the patch parametrizations are matching along the interfaces.While constructing C 0 smooth basis functions is quite straightforward, see, e.g., [10,38], imposing higher smoothness in a multi-patch setting is non-trivial.This makes solving higher order equations a more challenging task.In the following, we focus on fourth order problems, such as the biharmonic equation or a Kirchhoff-Love plate or shell formulation.
Two basic ways to get around this problem are to impose the C 1 smoothness weakly or strongly, that is, to adjust the variational problem or to construct special basis functions with higher smoothness, respectively.Following the first approach, one way to solve fourth order equations while keeping discontinuous, patch-wise basis functions is to employ a discontinuous Galerkin (dG) discretization as studied in [29,30].dG methods approximate the solution with patch-wise defined functions, which are discontinuous across patch interfaces.As a consequence, the variational formulation contains additional integral terms.See also [1,14,33].In this paper we use a Nitsche formulation for a C 0 smooth multi-patch discretization.In that case, a stability term, penalizing the jump of the normal derivative, is added to ensure the coercivity of the bilinear form.Hence, the corresponding stability parameter must be chosen sufficiently large.On the other hand, a too large stability parameter penalizes the jump of the normal derivative too much, which leads to locking of the solution.Therefore, the optimal range for the stability parameter must be determined, which we examine in more detail in this paper.
Another option to circumvent C 1 smoothness is to use a mixed/hybrid formulation as in [34,36,37,43].Here the PDE is reformulated in such a way, by introducing an extra field, that the resulting mixed formulation is of lower order.The obtained formulation has a saddle-point structure.To solve this problem efficiently, one needs a suitable preconditioner.Furthermore, it has a larger system to solve in comparison to the original problem.
At last, C 1 smoothness can be enforced weakly by using the mortar method, see, e.g., [5,15] for C 0coupling, where the coupling constraints are enforced using Lagrange multipliers, also resulting in a saddlepoint problem.In addition to the challenges of the saddle-point structure, there is the difficulty of finding a suitable discrete Lagrange multiplier space such that the resulting formulation is stable.
There are several different strategies when enforcing smoothness over multi-patch domains by strong C 1 coupling, i.e., constructing C 1 bases along the interfaces.One first attempt is the so-called bending strip method, see [26,27].More general formulations imposing geometric continuity over multi-patch domains are later developed in [9,13,18,24], while the related approaches presented in [25,31,32] follow a more local construction of geometrically continuous splines.See also [17] for a summary of related approaches.A significant problem of all strong C 1 coupling methods are the limitations they pose on the underlying geometries.It was shown in [9] that a standard isogeometric multi-patch discretization possesses optimal approximation properties only if the parametrization of the domain is a so-called analysis-suitable G 1 multipatch parametrization.Even though many geometries can be reparametrized, cf.[20], this is a significant restriction on the geometry.
It is nonetheless possible to use explicit C 1 constructions over more general multi-patch parametrizations by increasing the polynomial degree locally, as in [7,8], and/or by relaxing the smoothness conditions, as in [41,42].In [42] basis functions with higher polynomial degree and lower regularity are introduced locally at the interface between two patches, which are not exactly C 1 at the interface, but nevertheless yield optimal convergence rates in numerical tests.
In this paper, we focus on two aspects.One is to extend the approximate C 1 method introduced in [42] to multi-patch domains.Thus, a construction for basis functions at vertices needs to be developed.The second is to compare the approximate C 1 method with Nitsche's method, for which we derive conditions on the stability parameter.Thereby we see that the error obtained with the approximate C 1 method agrees with the error from Nitsche's method.
The outline of the paper is as follows: We start with basic notations for B-splines and multi-patch geometries needed for the paper, and introduce the C 0 isogeometric multi-patch space.Then the model problem, more precisely, the biharmonic problem is stated and we define two weak formulations, the standard and the Nitsche formulation.In Section 4 we give the C 1 smoothness conditions at the interfaces.Next, we present the basis constructions for the approximate C 1 method for the multi-patch geometries and introduce the C 1 space used for the approximate C 1 method.In Section 6 we analyse how the parameter in Nitsche's method must choosen for the method to be stable, which is numerically shown in Section 7.There we also compare Nitsche's method with the approximate C 1 method on several examples.

Preliminaries
In this section, we give a brief overview of B-splines and present some notation concerning the multipatch geometry.We start with the introduction of B-splines and their spaces.After that we present the definition of multi-patch geometries together with the underlying multi-patch topology, which is needed for the basis construction.The section concludes with the description of the C 0 isogeometric space.

B-splines
Given positive integers p, r and n, with r ă p, and a (uniform) mesh with mesh size h " 1{n, the open knot vector Ξ :" pξ 1 , ..., ξ N `p`1 q with N " p `1 `pp ´rqpn ´1q satisfies Ξ " p 0, . . ., 0 loomoon pp`1q´times q. (1) For simplicity, we assume that the regularity is the same at all interior knots and consequently all the knots have the same multiplicity.In general, the knot vector does not need to be uniform and the knot multiplicity may be different for different knots, see also [42].Since we assume uniform knot multiplicity p ´r, the inter-element continuity is defined by r, i.e., we have C r -smoothness at each knot.As usual, the B-spline basis functions b i , with i " 1, ..., N , can be constructed using the Cox-de Boor recursion, see [35].
We have Spp, r, hq " spantb i , i " 1, ..., N u which is the (univariate) spline space of degree p, regularity r and mesh size h.The spline space consists of functions which are piece-wise polynomials and C r in r0, 1s, more precisely, Spp, r, hq :" tw P C r pr0, 1sq : w| pih,pi`1qhq P P p for all i " 0, . . ., n ´1u. ( The univariate B-splines can be extended to the two-dimensional case by means of a tensor product structure.The multivariate spline space is defined in the parametric domain p Ω " r0, 1s 2 by Spp, r, hq where p " pp 1 , p 2 q, r " pr 1 , r 2 q and i " pi 1 , i 2 q denote the corresponding parameter pairs.The tensorproduct spline basis function is the product of two univariate basis functions, i.e., b i pu, vq " b i1 puqb i2 pvq where u, v P r0, 1s.We assume throughout the paper that p 1 " p 2 " p, r 1 " r 2 " r and h 1 " h 2 " h.

Multi-patch geometry
Let Ω be a bounded open subset of R 2 with a sufficiently smooth boundary BΩ.Moreover, let Ω be given through a multi-patch segmentation consisting of non-overlapping patches Ω pkq , k P M P " t1, ..., Ku, where K ą 0 is the total number of patches, i.e., Ω " ď with Ω pkq X Ω plq " H for all k ‰ l.Moreover, we assume that no hanging nodes exist.Each patch Ω pkq is a spline patch with the geometry mapping F pkq P pS pkq q 2 with where S pkq " Spp pkq , r pkq , h pkq q is a tensor-product spline space of patch Ω pkq .For simplicity, we assume that the spaces in all patches are the same, i.e., p pkq " p, r pkq " r and h pkq " h for all k P M P .The construction can also be used for different spline spaces with different degrees, regularities and mesh sizes as long as the interfaces are (partially) matching.We assume that all patch parametrizations F pkq are regular, i.e., det ∇F pkq pu, vq ě c ą 0, @ pu, vq P r0, 1s 2 , and the closure of Ω pkq possesses no self-overlaps, i.e., @pu, vq, pu 1 , v 1 q P p Ω : pu, vq ‰ pu 1 , v 1 q ñ F pkq pu, vq ‰ F pkq pu 1 , v 1 q.
In this paper, we introduce a local and a global notation for the mesh objects (i.e.edges and vertices) of the multi-patch.While the local notation describes the mesh objects of a single patch, the global notation concerns the relation of edges and vertices on the whole geometry.We introduce the pair pk, sq as a short notation for the index of the edge E pkq s .Then for the interface we introduce the index pair κ " pk, lq, if there exist pk, s k q, pl, s l q, with E pkq s k " E plq s l .Furthermore, we define the set M I which collects the ordered pairs of patch-indices for all interfaces and the set M E which collects the pairs pk, sq, each corresponding to a boundary edge.The sets are defined as M I " tκ " pk, lq | Ds k , s l P t1, . . ., 4u, s.t.E pkq s k " E plq s l and k ă lu, M E " tσ " pk, sq | Ds P t1, . . ., 4u, s.t.E pkq s P BΩu.
In addition, we denote the interface of κ " pk, lq P M I with the notation I κ just as we denote the boundary edges of σ " pk, sq P M E with the notation E σ .We assume that for two adjacent patches Ω pkq and Ω plq the patch parametrizations agree along the interface I κ , summarized in the following.
Assumption 1 (C 0 -conformity at the interfaces).The parametrizations of the two adjacent patches k and l meet C 0 along the interface I κ , with κ " pk, lq, i.e., there exists an Euclidean motion R κ : p E pkq s k Ñ p E plq s l (a mapping which is a combination of rotation, translation and reflection), such that F pkq pu, vq " F plq pR κ pu, vqq @ pu, vq P p E pkq s k .
Further, there exist three different types of vertices: (a) corner vertices, (b) interface-boundary vertices or (c) inner vertices.We have

R BΩ.
Similar to the edges, we introduce the set M V which collects all (unique) vertices V ι .Let ν be the valence of vertex V ι and Ω pk1q , ..., Ω pkν q be the patches around the vertex V ι .Then we have and the set is defined with ordered tuples of patch-indices: s kν and k ν ą ... ą k 1 and Epk 1 , s 1 q R tpk 1 , s k1 q, . . ., pk ν , s kν qu s.t.
Figure 1 gives an example of the nomenclature of the topology for a multi-patch parametrization.
Figure 1: The multi-patch notation: on the right side the global notation is defined and on the left side the local notation of patch k 1 where we also have the mapping F pk 1 q from the parameter domain p Ω " r0, 1s 2 to the patch Ω pk 1 q .We start with the local notation: each patch has four edges and four vertices.As an example, the local notation of patch k 1 is shown in the figure.Its local edges are indicated by E pk 1 q s and its vertices by V pk 1 q s , where s P t1, ..., 4u.In the global setting, the given geometry is constructed with three patches, indicated with Ω pk 1 q , Ω pk 2 q and Ω pk 3 q .Furthermore the geometry has three interface edges, labeled with Iκ, κ P t1, ..., 3u, six boundary edges, denoted as Eσ, σ P t1, ..., 6u and seven vertices, Vι, ι P t1, ..., 7u.From the topology of the geometry, we see that E pk 1 q 1 and E pk 1 q 2 are interface edges and E pk 1 q 3 and E pk 1 q 4 are boundary edges.Furthermore, the vertex V pk 1 q 2 is an interior vertex, the vertices V pk 1 q 1 and V pk 1 q 3 are interface-boundary vertices and the vertex V pk 1 q 4 a corner vertex.

C 0 isogeometric multi-patch spaces
In this section we define isogeometric function spaces over multi-patch domains.An isogeometric function ϕ : Ω Ñ R is defined such that ϕ ˝Fpkq P S pkq .
For this definition to be consistent at the patch-interfaces, we assume global C 0 -smoothness, resulting in the following definition.The C 0 isogeometric multi-patch space on Ω is given as X h " tϕ P C 0 pΩq such that w pkq " ϕ ˝Fpkq P S pkq @ k P M P u.
Since the patch parametrizations are meeting C 0 at the interface, see Assumption 1, and the knot vectors are assumed to be uniform, see (1), the meshes are conforming along the interfaces.Thus we achieve C 0smoothness at the interfaces if the corresponding spline coefficients at both sides of the interface are equal.Therefore a basis for the isogeometric C 0 multi-patch space X h can be constructed easily, cf.[2] for a more detailed description.

The model problem and its weak formulations
As a model problem we consider the biharmonic equation.We first set up the equation with two different combinations of boundary conditions.Then we present the weak formulation in the continuous setting, as well as in a discontinuous setting, following Nitsche's method.The section concludes with the two discrete formulations, a Nitsche discretization having discontinuous and a strong discretization having approximately continuous derivatives across interfaces.

The model problem
Let Ω be a multi-patch geometry as in Section 2.2 and let f : Ω Ñ R be a given source funtion on Ω.The biharmonic equation is given as with the boundary conditions where Γ N Y Γ L " BΩ, Γ N X Γ L " H and n is the unit outward normal vector on BΩ.The functions g 0 , g 1 and g 2 are all assumed to be sufficiently smooth.The boundary conditions concerning the function value and normal derivative can be imposed as essential boundary conditions, hence the problem can be homogenized, see also Section 5.3.So we assume from now on, that the problem is already homogeneous, that is, g 0 " g 1 " 0.

The standard weak formulation
We introduce the space where for any m P N and for any open, sufficiently smooth domain D the space H m pDq is the standard Sobolev space over D Ă R 2 with the standard scalar product pϕ, ψq H m pDq :" We denote by p¨, ¨qH 0 pDq " p¨, ¨qL 2 pDq the scalar product of the standard Lebesgue space L 2 pDq with the norm ¨ L 2 pDq .The weak formulation of problem ( 4)-( 6) is the following.
where B n is the normal derivative at the boundary.

The Nitsche formulation
Since the geometry is given by a collection of subdomains Ω " Ť k Ω pkq , a discontinuous Galerkin-type approach is a natural alternative to a fully conforming discretization.Therefore, we introduce for each m P N the broken Sobolev space H m pΩq :" tϕ P H 1 pΩq | ϕ| Ω pkq P H m pΩ pkq qu with the norm and inner product

¸1{2
and pϕ, ψq H m pΩq :" Moreover, we introduce the discretization space Let t˝u κ " 1 2 p˝p lq `˝pkq q| Iκ denote the average and ˝ κ " p˝p lq ´˝pkq q| Iκ denote the jump across the interface I κ .Since V 0 " H 2 pΩq X X 0 we have ϕ κ " 0 for any function ϕ P V 0 and for all interfaces I κ .
We follow Nitsche's method, that is, we reformulate the model problem following a symmetric interior penalty Galerkin approach.Thus, we consider the following problem, cf. in [29,30].
where η κ ą 0 is a prescribed stability parameter assigned to each interface I κ .
Conditions on the stability parameters η κ are discussed in Section 6.

The discrete formulations
In this subsection, we formulate two discrete problems.We have introduced the isogeometric space X h Ă H 2 pΩq in (3).Thus, it yields a suitable discretization space X h,0 " X h X X 0 for Problem 2. However, the space V h,0 " X h X V 0 " X h X C 1 pΩq is in general too restrictive, cf.[9].Hence, we introduce a different space r V h,0 ‰ V h,0 to discretize Problem 1.The construction of this space is described in Section 5.While the first space X h,0 fulfills the required conformity X h,0 Ă X 0 by definition, the second space does not in general fulfill the conformity relation r V h,0 Ę V 0 .The reason for this is that the space is spanned with basis functions that are not C 1 at the interfaces, but only approximately C 1 .However, in the limit the jump of the normal derivative across the interface vanishes by construction, see [42].Therefore, for the discretization we treat the functions from r V h,0 as if they were C 1 at the interfaces.As a consequence, the additional terms vanish in the variational formulation and no interface integrals need to be calculated.Furthermore, in some special cases, we achieve exact C 1 smoothness at the interface and we have r V h,0 Ă H 2 pΩq which is described in more detail in Remark 3.
Discretizing Problem 1 using r V h,0 , we obtain the following discrete problem.
Discretizing Problem 2 using the space X h,0 , we obtain the following discrete problem.
As pointed out before, the space V h,0 " X h X V 0 Ă H 2 pΩq is not a suitable discretization space for Problem 3, since its approximation power is in general drastically reduced.In Section 7 we compare with numerical experiments the solutions of Problem 4 and Problem 3.

Normal derivatives and C 1 smoothness conditions at interfaces
In order to solve fourth order problems on multi-patch domains, we need to give a description of the normal derivative of an isogeometric function across an interface.This is necessary both for the definition of the bilinear forms p¨, ¨qB h and p¨, ¨qC h as well as for the definition of the isogeomtric space r V h,0 , which we develop in detail in Section 5. Let us focus first on one edge of a patch Ω pkq , without loss of generality we consider the edge with u " 0, that is, " tF pkq p0, vq : v P r0, 1su.
We now define the tangential derivative along the edge to be tpvq :" B v F pkq p0, vq, and the unit tangent vector t 0 pvq :" tpvq τ pvq , where τ pvq " tpvq .We define the outer normal vector of BΩ pkq to be n k , with β pκ,kq pvq α pκ,kq pvq , with α pκ,kq pvq " det `Bu F pkq p0, vq, tpvq ˘, following the proof of [9, Proposition 2].Given a function ϕ : Ω Ñ R, with ϕ ˝Fpkq " f pkq , the normal derivative of ϕ along the edge can be described by " ´1 α pκ,kq pvq ´Bu f pkq p0, vq ´βpκ,kq pvqB v f pkq p0, vq ¯, since Here ∇ x denotes the gradient in physical space and ∇ denotes the gradient in pu, vq-coordinates, in particular, Hence, a patch-wise defined function ϕ is C 1 -smooth along the interface I κ between Ω pkq and Ω plq , iff Here we consider the restriction to E pkq s k to be the limit from the side Ω pkq , whereas the restriction to E plq s l denotes the limit from Ω plq .
In general, the isogeometric function is C 1 , if its graph parametrization is G 1 , which is determined by the relation For a further discussion on the statement above, see, e.g., [13,24].Due to this relation, the functions α pκ,kq , α pκ,lq and β pκ,kq , β pκ,lq are called gluing data.When constructing basis functions related to boundary edges, we introduce artificial gluing data, by setting α pk,sq " 1 and β pk,sq " 0 iff the edge E pkq s is a boundary edge.

The construction of the discrete space for the approximate C 1 method
In this section, we explain the basis construction for the space r V h,0 used in Problem 3. To do this, we first describe the basis functions in the local setting, i.e., patch-wise, and then define the global basis functions by gluing the functions together at the interfaces and at the vertices.
In the local setting, we introduce three different types of subspaces: the patch interior, the edge and the vertex space.Each patch Ω pkq can be described by these spaces, more precisely, each patch can be divided into nine subspaces: one patch interior, four edge and four vertex spaces, corresponding to the topology of the patch geometry.Each subspace (interior, edge and vertex space) is spanned by different basis functions.Let Ω pkq be the patch in which we have the interior space A pkq ˝, the four edge spaces A pkq E,s and the four vertex spaces A pkq V,s , s " 1, ..., 4. In Subsection 5.1 we describe the construction in the parameter domain.Figure 2 shows an overview of the local patch-wise separation.
Figure 2: The nine subspaces of patch k: one interior, four edge and four vertex spaces.In the right figure we give an example for the dofs.Thus, we choose the knot vector Ξ " p0, 0, 0, 0, 1{4, 1{2, 3{4, 1, 1, 1, 1q for both directions to obtain the number of dofs.Follow the construction in Subsection 5.1, we have nine interior dofs, three for each edge and six for each vertex.Similar to the element notation in FEM, the arrows at the edges represent the dofs for the normal derivative at the edge and the circles at the vertices the dofs for the first and second derivatives at the vertices.
To obtain the global basis functions, we need to match the basis functions at the interfaces and vertices to ensure the C 0 and the approximate C 1 continuity.This procedure is explained in Subsection 5.2.The resulting spaces are defined as A Iκ and A Vι .Further, we denote the boundary space as A Bσ for each boundary edge.As the result, we have

The patch-local subspaces
In this subsection we explain how the pull-backs x and A pkq V,s are constructed.We visualize them with the help of an example as shown in Figure 2.There we choose the knot vector Ξ " p0, 0, 0, 0, 1{4, 1{2, 3{4, 1, 1, 1, 1q in both directions.Before we can explain the construction in detail, we need to introduce the concept of the approximated gluing data in the following subsection.

The approximated gluing data
As shown in [22], a basis of the (exactly) C 1 -smooth isogeometric space can be constructed from the gluing data, which appear linearly in the formula similar to (16).The gluing data α pκ,kq and β pκ,kq defined in (11) are in general rather complex.While the functions α pκ,kq are splines from Sp2p ´1, r ´1, hq, the functions β pκ,kq are even piecewise rational functions with regularity r ´1.Thus, this results in the pull-back of the isogeometric function being a non-trivial rational function.To obtain "nicer" basis functions, that is, piecewise polynomials of a controlable, bounded degree, we introduce the approximated gluing data as spline functions, which are computed by a projection into Spr p, r r, hq with the operator P h , that is, r α pκ,kq :" P h pα pκ,kq q and r β pκ,kq :" P h pβ pκ,kq q.
For the numerical experiments, we fix the spline parameters of the approximated gluing data to r p " p ´1 and r r " r p ´1 to obtain optimal convergence rates, in accordance with [42].One can use a higher polynomial degree and/or lower regularity to approximate the gluing data to improve the approximation of the C 1 continuity (or even in some special cases to obtain an exactly C 1 -smooth space), but this does not improve the results in the numerical experiments, see [42].Similar to (10), we can express the approximate normal vector r n k for the patch Ω pkq on the edge E pkq 4 (u " 0) as where the functions in the linear combination are given as r a are defined equivalently.The choice of the spaces S `and S ´is derived from [19,Corollary 7] for AS-G 1 geometries, which also turned out to be the ideal choice to obtain optimal convergence rates in the numerical tests on general geometries, cf.[42].We have by construction ´, specifies the directional derivative in the direction of the approximate normal vector r n k .Hence, we have for The term B Recall that r n k « n k " ´nl « ´r n l .Here we consider the restriction to E pkq s k to be the limit from the side Ω pkq , whereas the restriction to E plq s l denotes the limit from Ω plq .
In other words, we impose an exact coupling of the approximate normal derivatives B r n k and B r n l .To obtain exact C 1 smoothness we need B r n k " B n k , which is achieved by r n k " ´r n l " n.For further discussions, see Remark 3. In the case of an boundary edge, we replace the approximate gluing data with α pk,sq " 1 and β pk,sq " 0.

The patch interior basis functions
The patch interior space is defined as where I pkq ˝" tpi 1 , i 2 q P Z 2 : 3 ď i 1 ď N ´2, 3 ď i 2 ď N ´2u and tb pkq j u are the basis functions of the tensor-product B-spline space S pkq " Spp, r, hq of dimension N ˆN , with p ě 2 and p ´1 ě r ě 1.In Figure 3 an example is shown.and normal derivatives at all edges of Ω pkq , that is, Moreover, the patch interior spline space satisfies x A pkq ˝Ă C 1 pr0, 1s 2 q.

The edge basis functions
Without loss of generality, let the edge E pkq s be such that it corresponds to u " 0, i.e., s " 4. Any edge can be rotated and translated in such a way that it coincides with this configuration.
We define the space x A pkq E,s as the span of those basis functions that have non-vanishing trace or approximate normal derivative along the interface and vanishing value, gradient and Hessian at both endpoints of the interface.More precisely, we have x A pkq E,s,`" spantf pkq s rb j , 0s : 4 ď j ď N `´4u and x A pkq E,s,´" spantf pkq s r0, b j s : 3 ď j ď N ´´3u, where f pkq s r¨, ¨s is defined as in ( 16), tb ì u and tb j u, N `and N ´are the bases and dimensions of the spaces S `and S ´, respectively.In Figure 4 we give an example of the edge space.

The vertex basis functions
For simplicity of the notation, we assume that the vertex is at x " V pkq 1 " F pkq p0, 0q.Then we collect all the edge basis functions, which have non-vanishing C 2 -data on one of the two adjacent edges, i.e., which fulfill one of the following conditions More precisely, we define three sets of basis functions B b.e." tf To construct the vertex space we perform C 2 interpolation at the vertex for all three sets of functions.To do so, we prescribe C 2 -data Φ " tΦ i1,i2 u i1,i2 , with 1 ď i 1 , i 2 ď 3 and i 1 `i2 ď 4, in physical space and interpolate the pull-backs using the three spaces defined above.We then add the first two interpolations and subtract the third.The resulting functions are denoted by g pkq 1 rΦ i1,i2 s.We refer to Appendix A for details of the construction.We have by construction g pkq s rΦ i1,i2 s P Spp `r p ´1, mintr r, r, p ´2u, hq b Spp `r p ´1, mintr r, r, p ´2u, hq.
The space for the vertex V pkq s is defined as where I V " tpi 1 , i 2 q P Z 2 : 1 ď i 1 , i 2 ď 3, and i 1 `i2 ď 4u.Thus, the dimension of the space x A pkq V,s is always six. Figure 5 shows an example for the vertex space.

Construction of the global space
In this subsection we first describe the coupling conditions.Then, the coupling conditions are used to connect the local (patch-wise) spaces to define the global space r V h .Considering one interface I κ between patches Ω pkq and Ω plq , we assume for all isogeometric functions ϕ h P r V h that Moreover, we assume that for each vertex V ι the functions ϕ h P r V h are C 2 -smooth at the vertex, that is, the limit of the function value, gradient and Hessian is the same on all patches sharing the vertex V ι .
Remark 2. We recall the estimate from [42, Theorem 1], yielding bound where ϕ h is defined satisfying ( 22)-( 23).Here C ą 0 depends on the geometry, but not on the mesh size.Therefore, the C 1 error depends on the choice of the polynomial degree for the approximate gluing data.We also observe that higher polynomial degree and/or lower regularity for the approximated gluing data does not lead to better results in the numerical experiments, see [42].
where the patch interior spaces are defined as A pkq ˝" tϕ h P C 0 pΩq : ϕ h ˝Fpkq P x A pkq ˝and ϕ h ˝Fplq " 0 for all l ‰ ku, the interface spaces as ˝Fpmq " 0 for all m R tk, lu and ϕ h satisfies (23) , for all i " 1, . . ., ν ϕ h ˝Fplq " 0 for all l R tk 1 , . . ., k ν u and / -, and the boundary edge spaces as A basis for the global space can be derived immediately from the local bases.For patch interior and  s kν rΦ pkν q i1,i2 s are coupled for each index pair pi 1 , i 2 q P I V .In Figure 6 an example visualizing edge and vertex basis functions is given.Since the structure of the construction is similar to the AS-G 1 construction in, for example, [19,22] we obtain linear independent basis functions for the space r V h as stated in the following lemma.
Lemma 4. The space r V h is the direct sum of its subspaces A pkq ˝, A Iκ , A Vι and A Bσ .Moreover, the coupling described above yields a basis for each of the subspaces, which in turn yields a global basis.Remark 3. We now briefly discuss the C 1 smoothness of the space r V h which is discussed in more details in [42,Subsection 6.3].In some cases, the C 1 condition in (23) at the interface I κ is exact, i.e., r n k " ´r n l .
A sufficient condition would be r n k " n k " ´nl " ´r n l .This is the case, when the projection of the gluing data is exact, i.e., r α pSq pvq " α pSq pvq and r β pSq pvq " β pSq pvq for S P tk, lu.The condition holds for all v P r0, 1s if the gluing data satisfies α pkq , β pkq P Spr p, r r, hq.Then the condition in ( 23) is actually an exact C 1 condition and the jump in (24) vanishes.As a result we then have A Iκ Ď C 1 pΩq.

Recall the boundary conditions
We assume that each set Γ N and Γ L is the union of boundary edges of patches, i.e., the boundary conditions can change only at vertices of the multi-patch domain.
The boundary condition ∆ϕ " g 2 is naturally enforced in the equation on the right hand side.The other two boundary conditions are enforced strongly by encorporating them into the space.To do so we need to define functions that satisfy the boundary conditions for general g 0 and g 1 , as well as functions that have homogeneous boundary conditions spanning the space r V h,0 " tϕ h P r V h : ϕ " 0 on BΩ and B n ϕ " 0 on Γ N u.
We collect all functions which are used to approximate g 0 and g 1 in the space r V h,BΩ Ă r V h .We then have Similar to the global space r V h , we split the space r V h,BΩ into separate contributions where M V,BΩ denotes the indices of all boundary vertices.Let B σ " E pkq s be a boundary edge.For B σ Ă Γ N we set A Bσ,BΩ " A Bσ , collecting all functions from Let V ι be the vertex at the boundary and the corresponding vertex space is constructed by To obtain the correct subspace A Vι,BΩ , we compute the kernel of the space, evaluated with the value, that is with λ i P R, i P t1, ..., 6u : For the boundary space, we use the functions which span Note that the dimension of the kernel and of the boundary space depends on the considered boundary conditions and on the geometry.
Remark 4. To obtain an exact kernel at the vertices, an interpolation of the approximate gluing data at the boundary points is required.However, the interpolation can be omitted.In this case, the kernel must be calculated to a h-dependent tolerance.

Existence and uniqueness of the solution using Nitsche's method
In this section, the optimal choice of the stability parameter yielding coercivity and boundedness of the bilinear form of Problem 4 is derived.Thus, the existence and uniqueness of the solution of Problem 4 is shown.We prove the coercivity and boundedness of the form in the following dG-norm ϕ 2 X h " pϕ, ϕq X h where pϕ, ψq X h :" p∆ϕ, ∆ψq H 0 pΩq `pϕ, ψq C h .
We start with a description of the discrete space for Nitsche's method.Let the space X pkq h be the standard tensor-product spline space for the patch k.Then the space X h is the collection of the patch spaces X pkq h and, in addition, the dofs along the interfaces are matching to ensure C 0 smoothness, see, e.g., [10,38].For the Problem 4, we set X h,0 " X h X X 0 to fulfill the boundary conditions.In order to prove the coercivity and boundedness, we need to bound the average at the interface which is bounded as follows: Lemma 5. Let ϕ h P X h .For any I κ with κ P M I , we have where c κ phq ą 0 is a constant, which depends on the mesh size h and on the patch parametrizations, but not on the function ϕ h .
Proof.Using the definition of the average, we have with the triangle inequality Since ϕ h is from a finite dimensional space, the supremum sup exists and we denote it by c k .Similarly, we obtain an upper bound c l on Ω plq .Thus, the proof is complete with c κ " 1 2 pc k `cl q, which does, in general, depend on h.
Assumption 2. For any I κ with κ P M I , the constant c κ phq from Lemma 5 satisfies where c κ is an h-independent constant.
In practice, the constant c κ phq in Lemma 5 can be computed by a generalized eigenvalue problem, following the same steps as in [11], for example.For η κ sufficiently large, Nitsche's formulation as given in Problem 4 is coercive and bounded which is stated in the following Theorem.
Theorem 1.For all ϕ h , ψ h P X h,0 , let η κ be such that η κ ą 16 h c κ phq, then pϕ h , ϕ h q A h ě µ ϕ h Proof of Theorem 1.By definition, we have With the help of Young's inequality with δ κ ą 0 and Lemma 5 we have we have for the coercivity thus the bilinear form is coercive, if for all κ.This means that the stability parameter η κ must satisfy η κ ą 16 h c κ phq and, under Assumption 2, that η κ ą 16 c κ .
For the boundedness we obtain from the Cauchy-Schwarz inequality, the triangle inequality, Lemma 5 and the fact that s P t1, ..., 4u for all ϕ h , ψ h P X h,0 .Using the estimate above and the Cauchy-Schwarz inequality, we obtain boundedness of the bilinear form If moreover Assumption 2 is satisfied, then the constant is independent of h.
Using Theorem 1, we can apply the Lax-Milgram theorem, e.g., as in [12], and consequently have existence and uniqueness of the solution to Problem 4. The existence and uniqueness of the solution of Problem 3 is straightforward.

Numerical experiments
In this section we perform numerical experiments on four C 0 multi-patch geometries -denoted by Example I-IV.On a biharmonic model problem we compare the two methods considered in this work, i.e., the approximate C 1 discretization and Nitsche's method.More precisely, we solve on each geometry Problem 3 using the discrete space r V 1 h,0 and Problem 4 using the discrete space X h,0 .For simplicity, let p " p As boundary conditions in ( 5)-( 6) we consider Γ L " H for Example I and II, and Γ N " H for Example III and IV.
Let ϕ h be the discrete solution of either Problem 3 or Problem 4. Since the exact solution is smooth, we expect in both cases optimal convergence rates in the mesh size h, i.e., ϕ ´ϕh H 2 pΩq " Oph p´1 q. (27) In Subsection 7.1, the four geometries are presented.A comparison of the errors for different polynomial degrees is given in Subsection 7.2.In Subsection 7.3 we study the influence of the solution using Nitsche's method on the choice of the stability parameter.Finally, in Subsection 7.4 we conclude with a comparison of the (exact) C 1 -smooth discretization on the reparametrized AS-G 1 geometry, following the approach presented in [20], with the approximate C 1 method.All tests are implemented within the open-source C++ library G+Smo, cf.[28].

The model geometries
We consider four geometries, the first two describe the same domain, that is, the unit square, but with different parametrizations, the last two are more application-oriented geometries.The four geometries are shown in Subfigures 7a-7d.Example I consists of six bilinear patches and the gluing data is consequently linear.It follows that we have exact C 1 smoothness at the interfaces, see Remark 3. In contrast to Example I, Example II is made of bicubic patches and has curved interfaces.The gluing data of this geometry is not linear and therefore, the discete space is only approximate C 1 .The same is true for Example III and Example IV: both geometries have non-linear gluing data and hence the spaces are not exact C 1 .Example III describes a turtle with bicubic patches, while Example IV is a part of a picture of a car and is inspired by examples from [6,20].All patches in Examples I-IV are Bézier patches.The corresponding exact solutions are depicted in Figures 7e-7h

Convergence analysis
We compare the convergence rates of the errors measured in the L 2 -, H 1 -and H 2 -norms for the polynomial degrees p P t3, 4, 5u and use the maximum regularity r " p ´1.In all four examples we compute the error using the approximate C 1 method and Nitsche's method, represented by a dashed and solid line, respectively.To obatin the stability term in Nitsche's method, we solve Lemma 5 with the eigenvalue problem at a fixed h 0 and choose where c κ is the largest eigenvalue.
In Figure 9, we plot the results for both methods and see that the H 1 -and H 2 -errors differ only slightly.In the plots, these lines almost overlap and converge optimally with the rate stated in (27).In the L 2 -norm, the error for the approximate C 1 method in Examples I and II is almost the same as the error for Nitsche's method, while in Examples III and IV the error for Nitsche's method is slightly smaller.One reason could be, that the approximate C 1 basis is slightly more restrictive near boundary vertices.Nevertheless, we see in the examples that both methods solve the biharmonic equation optimally with similar error values.
For Example II, we additionally compare the errors of the approximate C 1 method and Nitsche's method with a single patch parametrization on the same domain (unit square).There we plot the results as a function of the number of dofs, see Figure 8.As expected, the single patch parametrization needs the least number of dofs for a given error.The approximated method and Nitsche's method are almost the same: the approximate C 1 method needs a slightly smaller number of degrees of freedom (dofs) than Nitsche's method for a fixed mesh-size.However, it can be seen that the gap between the single patch and the other two methods, with smaller mesh-size, becomes smaller.

Dependence on the stability parameter in Nitsche's method
In this subsection we examine the dependence of the error using Nitsche's method on the stability parameter.We assume that the parameter is chosen globally, i.e., the parameter is the same for each interface.Therefore, the bilinear form p¨, ¨qC h changes to pϕ, ψq C h " η h In our study we vary the stability parameter over a range from 4h 0 ¨10 ´3 to 4h 0 ¨10 4 and compute the errors for a fixed mesh size h 0 .Figure 3 shows the errors in the L 2 -, H 1 -and H 2 -norm for mesh size h 0 " 1{2 6 and different polynomial degrees.In comparison, we plot the error obtained with the approximate C 1 method with a dashed line, which is independent of the stability parameter.From the numerical results it is evident that Nitsche's method does not converge properly for large values of the stability parameter.The reason for this is that a too large parameter leads to an over-penalization of the jump of the normal derivative, which, in general, leads to locking of the solution.On the other hand, a too small stability parameter leads to instability.Also, we see a 'spike' occuring for a value of η close to the constant c κ phq in Lemma 5. A similar behavior with the occurence of 'spikes' is also observed in [11] where the authors use Nitsche's method for imposing Dirichlet boundary conditions.7.4.Comparing an AS-G 1 reparametrization with the approximate C 1 discretization In the last example, we follow the reparametrization strategy as in [20].Here, a non AS-G 1 geometry is reparametrized into an AS-G 1 geometry.The consequence is that although the geometry and the discretization spaces can be constructed with exactly C 1 -smooth basis functions as in [20,22], the geometry may change as a result of the reparametrization.We choose the geometry of Example IV and reparametrized the non-AS-G 1 geometry.Figure 11 shows the results.There, the differences between the non-AS-G 1 geometry and the AS-G 1 geometry are shown represented in black (original) and red lines (reparametrized).For both methods we chose the maximum possible regularity.That is, we select r " p ´1 for the approximate C 1 method and r " p ´2 for the AS-G 1 discretization -which is a necessary restriction derived from the construction, cf.[9].If r " p ´1 then the C 1 -smooth subspace X h X C 1 pΩq reduces in general to global polynomials when restricted to one interface, which results in locking of the numerical solution.In Table 1, for fixed mesh size h " 1{2 5 , the number of dofs and the errors are given.There we see that both methods yield very similar errors.However, the major difference between both methods is the number of dofs, which can be explained by the different regularities: since the AS-G 1 discretization requires r " p ´2, the approximate C 1 method with r " p ´1 needs much fewer dofs to obtain the same error levels.Table 1: The dofs and errors for the AS-G 1 geometry and the non-AS-G 1 geometry.

Conclusion and future work
We extend the basic construction from [42] to general multi-patch domains.Therefore, we introduce a construction for basis functions around vertices using interpolation of functions that are approximately C 1 -smooth across interfaces.Three different kinds of spaces are created in the construction: the patch interior spaces, the edge (interface and boundary) spaces and the corner spaces (both for boundary vertices and inner vertices), which are derived from the topology of the multi-patch geometry.This creates spaces that locally possess higher polynomial degrees and lower regularity, with the exception of the patch interior space, which is a standard isogeometric space.As a result we get non-nested spaces.In contrast to discretization spaces over AS-G 1 parametrizations, as in [22], which require r ď p ´2, the approximate C 1 method also allows us to choose spline spaces of maximum regularity r " p ´1.
Moreover, we compare the approximate C 1 method with Nitsche's method.In the numerical experiments we see that both methods converge optimally and the error values are almost the same.While one has to determined a suitable stability parameter for Nitsche's method, no such tuning is needed for the approximate C 1 method.Thus, to summarize, the approximate C 1 method provides an explicit and simple to implement alternative to weak (H 2 -nonconforming) and exact (H 2 -conforming) methods to solve fourth order problems, exemplified on a biharmonic model problem.The advantages of the approximate C 1 method are that the method can be applied on any C 0 -conforming multi-patch parametrization and does not depend on any non-trivial parameter choices.
In the future we want to study several aspects of the method, such as convergence and stability properties and extensions that result in nested spaces.This would allow an adaptive construction with THB-splines, following the work as in [3,4].Moreover, we want to extend the construction to C 0 -non-conforming (nonmatching) interfaces and to surface domains.In such a context, the approximate C 1 method could be a viable option to discretize Kirchhoff-Love shell problems.Another possible direction of research is the extension to volumetric domains, where C 1 -smooth discretizations, in general, yield suboptimal convergence rates, cf.[23].Since the approximate C 1 method has no interface integrals, it would be interesting to combine the method with a multigrid solver, see [39,40].

2 .
C 1 expansion along one edge For the construction of the basis functions, we use a Taylor expansion of the trace and of the transversal derivative as stated in [19, Proposition 5].For simplicity, we consider again the edge p E pkq 4 on Ω pkq , as in Section 4. Then the functions are defined for all u, v P r0, 1s by f pkq 4 rb `, b ´spu, vq " b `pvq pb 1 puq `b2 puqq `´r α pκ,kq pvqb ´pvq `r β pκ,kq pvqpb `q1 pvq ¯h p b 2 puq, (16) where b `P S `" Spp, p ´1, hq and b ´P S ´" Spp ´1, p ´2, hq.Representations f pkq s rb `, b ´s for the functions along other edges p E pkq s

f pkq 4
rb `, b ´s P S 1 pp, r, hq b S 2 pp `r p ´1, mintr r, r, p ´2u, hq.While the first variable in the function, i.e., b `, describes the trace at the edge, the second variable, i.e., b

r n k ϕLemma 1 .
h is an approximation of the normal derivative at the edge due to the definition of the approximated gluing data, see[42, Proposition 4].For the interface I κ , if the pull-back of ϕ h is given by f pkq s k rb `, b ´s on Ω pkq and by f plq s l rb `, b ´s on Ω plq , then using (15) and (16) gives us B r n k ϕ h | E pkq s k " ´Br n l ϕ h | E plq s l .

Lemma 2 .
Let b pkq j P x A pkq ˝, then the isogeometric function ϕ h | Ω pkq " b pkq j ˝pF pkq q ´1 has vanishing traces The dofs of the patch interior space.(b) The basis functions of the patch interior space.

Figure 3 :
Figure 3: The patch interior basis: While the red dots represent the dofs that are eliminated, that is, the corresponding basis functions have non-zero function values or non-zero normal derivatives at the patch boundary, the black dots represent the dofs for the patch interior space.

Lemma 3 .Figure 4 :
Figure 4: We consider here the edge with u " 0. For each edge we obtain seven trace basis functions depicted with a square and six normal derivative basis functions represented by a pentagon.Then the edge basis functions which have influence up to order two at the vertices are eliminated (marked in red).The remaining three basis functions span the edge space and are shown in Subfigures 4b-4d.

Figure 5 :
Figure 5: As an example the basis functions at the vertex u " v " 0 are shown in Subfigures 5b-5g.Those basis functions are computed by the C 2 -interpolation at the vertex.

Figure 6 :
Figure 6: Examples of global basis functions at the interface (marked with a thick line) and vertex (marked with a dot).Only the relevant patches are plotted.
Ex. IV: exact solution

Figure 7 :
Figure 7: The geometries which are used for the numerical results and their exact solutions.

Figure 8 :
Figure 8: The errors versus the dofs.

C 1 Figure 9 :
Figure 9: Convergence rate with different polynomial degrees.

Figure 11 :
Figure 11: Non-AS-G 1 geometry (black) vs. AS-G 1 reparametrization (red) Starting with the local setting, each single patch Ω pkq has four edges E for I κ