Scaled boundary isogeometric analysis with C1 coupling for Kirchhoff plate theory

Although isogeometric analysis exploits smooth B-spline and NURBS basis functions for the definition of discrete function spaces as well as for the geometry representation, the global smoothness in so-called multipatch parametrizations is an issue. Especially, if strong C1 regularity is required, the introduction of function spaces with good convergence properties is not straightforward. However, in 2D there is the special class of analysis-suitable G1 (AS-G1) parametrizations that are suitable for patch coupling. In this contribution we show that the concept of scaled boundary isogeometric analysis fits to the AS-G1 idea and the former is appropriate to define C1-smooth basis functions. The proposed method is applied to Kirchhoff plates and its capability is demonstrated utilizing several numerical examples. Its applicability to non-trivial and trimmed shapes is demonstrated.


Introduction
Isogeometric Analysis (IGA) is a concept which was introduced by Hughes et al. [1] and developed to a widely used and very successful approach within numerical analysis and computational mathematics. It connects the fields of geometry representation in Computer Aided Design (CAD) with the framework of Finite Elements (FE). For both, geometry parametrization and discretization space definition, one utilizes B-splines and NURBS (Non-Uniform Rational B-Splines). In IGA we find different built-in features that allow for interesting applications. For example, at least in the so-called single patch case, smoothness of test and ansatz functions is no issue anymore since the underlying B-spline and NURBS basis functions can be chosen with high global regularity. This is an advantage compared to classical FE ansatz spaces which are mostly continuous or discontinuous mappings. Clearly, also C 1 -smooth FE spaces can be defined without the ideas of IGA. But, for this regularity requirement and in particular if we look for C k basis functions where k > 1, the definitions become complicated. Whereas utilizing splines, a regularity increase can be implemented efficiently. This comes along with the k-refinement ansatz. The latter means a simultaneous degree and regularity elevation step. Such a change of the discretization spaces has no pendant in standard FE theory. The possibility to vary regularities in IGA implies two important aspects. On the one hand, IGA is a suitable discretization approach if one deals with high-order problems. On the other hand, increasing the regularity without changing the degree leads to a significant reduction of degrees of freedom although the approximation behaviour changes only slightly. Especially, as long as p > r, where p stands for the B-spline degree and r for regularity, the observed convergence rates are basically the same. For more information on IGA we recommend [1][2][3]. Unfortunately, if one considers complex geometries, multipatch IGA meshes are necessary, i.e. the computational domain is decomposed in several subdomains that are associated to isogeometric spaces. These patch-wise spaces have to be coupled in order to get high regular basis functions on the whole domain. And this coupling is in general not trivial. Basically, when coupling the patches one faces similar problems like for the coupling of classical FE mesh elements. This means, the smoother the basis functions across patch interfaces, the harder the coupling procedure gets. Even worse, for different geometries locking effects arise, i.e. bad approximation results, if we want strong global C r -smoothness with r ≥ 1; see [4]. Consequently, in literature one can find the concept of weak patch-coupling to define proper multipatch IGA spaces; see [5][6][7]. But in 2D, if one restricts oneself to a special class of parametrizations, the so-called analysis-suitable G 1 parametrizations, the strong C 1 coupling of patches is possible; cf. [4]. Besides, for the latter geometries the problem of C 1 -locking can be avoided. Further we note the framework of polar splines, see e.g. [8], that are useful to handle disk like domains and to define globally smooth basis functions. In this article, we combine the approach and the results from the mentioned reference [4] to show that scaled boundary IGA (SB-IGA) is suitable to introduce C 1 -smooth basis functions. In SB-IGA we assume the domain to be defined by a scaling of its boundary w.r.t. to some given scaling center; see [9][10][11]. At first glance, this requires star-convex domains and introduces a singular parametrization mapping. Nevertheless, it fits to the fact that in CAD the computational domains are often represented via its boundaries and furthermore we explain in the subsequent parts why it is still useful for non-star-convex geometries. To demonstrate the C 1 regularity and capabilities we apply our approach to the problem class of Kirchhoff plates. The related strong PDE formulation is of fourth order and hence in the classical weak formulation second derivatives appear ; see [12,13]. Furthermore, we briefly explain, why the proposed scaled boundary meshes (SB-meshes) are convenient when dealing with trimming. Latter domain modification is often applied in IGA and central if topologically complicated domains appear. We refer to [14] for details regarding trimming. More precisely, the structure of the proposal is the following. We start with a very brief discussion of B-splines and SB-IGA in Sec. 2 to clarify the mathematical notions. We proceed with two sections, Sec. 3 and Sec. 4, dedicated to the coupling problem. We discuss the standard two-patch case as introduced in [4] and in a second section we show the application to SB-meshes. Afterwards, in Sec. 5, we briefly explain generalization possibilities, where we emphasize the application to trimmed geometries. In Sec. 6 we look at several numerical examples in the context of Kirchhoff plate theory. In Sec. 7 we discuss the problem of stability coming along with the singular parametrizations. We close with a short conclusion in Sec. 8.

B-splines and SB-IGA
In this section we want to introduce briefly some basic notion and clarify the SB-IGA ansatz. First, we state a short overview of B-spline functions, B-spline spaces, respectively. Exploiting [2,3] for a short exposition, we call an increasing sequence of real numbers Ξ := {ξ 1 ≤ ξ 2 ≤ · · · ≤ ξ n+p+1 } for some p ∈ N knot vector, where we assume 0 = ξ 1 = · · · = ξ p+1 , ξ n+1 = · · · = ξ n+p+1 = 1, and call such knot vectors p-open. Further the multiplicity of the j-th knot is denoted by m(ξ j ). Then the univariate B-spline functions B j,p (·) of degree p corresponding to a given knot vector Ξ is defined recursively by the Cox-DeBoor formula: and if p ∈ N ≥1 we set where one puts 0/0 = 0 to obtain well-definedness. The knot vector Ξ without knot repetitions is denoted by {ψ 1 , . . . , ψ N }. The multivariate extension of the last spline definition is achieved by a tensor product construction. In other words, we set for a given tensor knot vector Ξ := Ξ 1 × · · · × Ξ d , where the Ξ l = {ξ l 1 , . . . , ξ l n l +p l +1 }, l = 1, . . . , d are p l -open, and a given degree vector p := (p 1 , . . . , p d ) for the multivariate case with d as the underlying dimension of the parametric domain Ω = (0, 1) d and I the multi-index set B-splines fulfill several properties and for our purposes the most important ones are: • If for all internal knots the multiplicity satisfies 1 ≤ m(ξ j ) ≤ m ≤ p, then the B-spline basis functions B i,p are globally C p−m -continuous. Therefore we define in this case the regularity integer r := p − m. Obviously, by the product structure, we get splines B i,p which are C r l -smooth w.r.t. the l-th coordinate direction if the internal multiplicities fulfill 1 ≤ m(ξ l j ) ≤ m l ≤ p l , r l := p l − m l , ∀l ∈ 1, . . . , d in the multivariate case. We write in the following r := (r 1 , . . . , r d ), for the regularity vector to indicate the smoothness. In case of r i < 0 we have discontinuous splines w.r.t. the i-th coordinate direction. To later to emphasize the regularity of the splines we introduce a upper index r and write in the following B r i,p , B r i,p , respectively.
• The support of the spline B r i,p is contained in the interval [ξ i , ξ i+p+1 ].
• We have the partition of unity property i B r i,p = 1.
The space spanned by all univariate splines B r i,p corresponding to given knot vector and degree p and global regularity r is denoted by For the multivariate case one can define the spline space as the product space of proper univariate spline spaces. To obtain more flexibility it could be useful to introduce a strictly To define discrete spaces based on splines we require a mapping F : Ω := (0, 1) d → R d which parametrizes the computational domain. The knots stored in the knot vector Ξ, corresponding to the underlying splines, determine a mesh in the parametric domain Ω, namely M : under the mapping F, i.e. M := {F(K) | K ∈ M }, gives us a mesh structure in the physical domain. By inserting knots without changing the parametrization we can refine the mesh, which is the concept of h-refinement; see [1,2]. For a mesh M we can introduce the global mesh size through h := max{h K | K ∈ M }, where for K ∈ M we denote with h K := diam(K) the element size and M is the underlying parametric mesh.
The underlying idea of SB-IGA takes into account that in CAD applications the computational domain is often represented by means of its boundary. As long as the region of interest is starshaped we can choose a scaling center x 0 ∈ R d and the domain is then defined by a scaling of the boundary w.r.t. to x 0 . In the planar case, which is the one we focus on here, and in view of the isogeometric analysis we have some boundary NURBS curve γ(ζ) = i C i N r i,p (ζ), C i ∈ R 2 and define the SB-parametrization of Ω through Fig. 1).
Depending on the situation it might be useful to allow for more flexibility and we replace in this article the prefactor ξ by a degree 1 polynomial of the form q(ξ) := c 1 ξ + c 2 , c 1 > 0, c 2 ≥ 0, i.e. . . . In other words, there are so-called control points C i,j ∈ R 2 associated to the NURBS (ζ, ξ) → N r i,p (ζ) B r j,p (ξ) ∈ N r p ⊗ S r p which define F, namely For reasons of simplification we suppose equal degree and regularity w.r.t. each coordinate direction. Due to the SB ansatz we obtain in the physical domain Ω layers of control points and it is C 1,1 = C 2,1 = · · · = C n1,1 ; cf. Fig. 2.
Remark 1. The structure of the domain and the control points in Fig. 2 are reminiscent of the mentioned polar splines framework; see [8]. In fact, the idea of degenerating an edge of the standard parametric domains is analogous to the SB-IGA ansatz followed in this article. In some sense, one can interpret the polar spline approach as a special case of a scaled boundary representation. But we emphasize that there are also several differences between the mentioned polar spline approach and the content of this publication. First of all, we do not work with periodic spline functions and the computational domain boundary domain can be a non-smooth curve. Secondly, our treatment of the scaling center basis functions differs, namely, we use always the original B-spline functions from the non-degenerate parametric domain, whereas in [8] triangular Bernstein polynomials are considered.
In particular, we do not construct polar spline basis functions. Besides, in the subsequent parts we look at the coupling of different SB-parametrizations as well as at the coupling of different star-shaped blocks. This leads in general to geometries that can not be treated directly with polar splines.
The isogeometric test functions, starting point for discretization methods, are then defined as Figure 3: The boundary can be determined by the concatenation of several curves. In this situation the patch-wise defined discrete function spaces have to be coupled in order to obtain the wanted global smoothness.
the push-forwards of the NURBS, namely If the domain boundary ∂Ω is composed of different curves γ (m) , one defines parametrizations for each curve as written above and we are in the field of multipatch geometries; cf. Fig. 3. To be more precise, for a n-patch geometry we have m=1,...,n and F m is defined analogous to (5). Unfortunately, especially if the curves meet not G 1 but high regularity of the corresponding IGA test functions is required, then the coupling gets involved. IGA spaces in the multipatch framework are straight-forwardly defined as denotes the IGA space corresponding to the m-th patch, to F m , respectively. For all the patch coupling considerations we suppose the next assumption.

Assumption 1 (Regular patch coupling). sdakjhjkh
• In each patch we use NURBS and B-splines with the same degree p and regularity r.
• The control points at interfaces match, meaning the control points of meeting patches coincide at the respective interface.
Thus it is justified to write for the set of parametric basis functions in the m-th patch for proper n (m) 1 , n 2 ∈ N >1 . The main part of the article is dedicated to the C 1 coupling of such SB-IGA patches, i.e. face spaces of the form where the singularity of the F i at x 0 for c 2 = 0 requires a special attention.

Classical planar two-patch coupling
For reasons of simplification, we restrict ourselves until Sect. 5 to the situation of two-patch geometries. The aspects of coupling can be straightforwardly generalized to three or more patches. First, before we come back to SB-IGA we turn towards the simple case of a classical two-patch parametrization, meaning there are no singular points for the F m . For such a situation we explain the needed conditions for a C 1 coupling and exploit the notion of analysis-suitable G 1 planar parametrizations. Here we mainly focus on the results and the framework of [4]. The setting is now the following. Assume that we have two mappings F (S) , S ∈ {L, R} corresponding to the left and right side of a two-patch situation as displayed in Fig. 4. Figure 4: A planar non-degenerate two-patch domain.
One notes the colors in the mentioned Fig. 4 that should determine the parametrization orientation used here and in the following. Further we assume that the F (S) are diffeomorphisms and write Γ = F −1 (Γ) := {(0, ξ) |ξ ∈ [0, 1]} for the interface. And, using the notation from the previous section, we can write F 1 := F (L) (· − 1, ·), F 2 := F (R) , i.e. there is a shift in the parametric domain. In view of isogeometric analysis we suppose two functions g (L) ∈ C 1 (Ω (L) ), g (R) ∈ C 1 (Ω (R) ) that can be combined to a continuous mapping g : Ω → R defined by g |Ω (S) = g (S) and we want to know which properties of the g (S) lead to a C 1 -regular g. The latter is the case if and only if the graph G := {(F(ζ), g(F(ζ))) | ζ ∈Ω} as a surface in 3D has well-defined tangent planes along the patch interface. This condition can be formulated by means of auxiliary coefficient functions α (S) , β as written down in the next lemma which uses Proposition 2 and Definition 2 of [4].
Proof. This follows directly by Definition 2 and Proposition 2 in [4].
Note that in (8) the last term of the sum on the left side does not change if R is replaced by L, due to the assumed continuity of F and g. In the IGA context, one considers cases like ĝ (L) (· − 1, ·) ∈ N r,r p,p andĝ (R) ∈ N r,r p,p . More precisely, the standard two-patch isogeometric spaces without coupling conditions are defined in the parametric domain through Although there is a clear criterion for C 1 regularity, the actual calculation of test functions meeting the conditions within the scope of numerical methods is in general not trivial. Especially if test functions are defined separately for each patch, like for the isogeometric spaces, and a suitable global C 1 -smooth linear combination of the test functions is sought one can observe bad approximations properties and a loss of convergence in different situations. This is a reason why in [4] a class a special parametrizations are introduced leading to optimal convergence of the isogeometric spline spaces under h-refinement, namely the so-called analysis-suitable G 1 parametrizations.
Then the parametrization F is called analysis-suitable G 1 .
As already mentioned, the function spaces of interest are C 1 -regular spline spaces Proof. Compare Theorem 1 in [4].
An important point is the appearance of C 1 locking if we choose p = r + 1 even if F is AS-G 1 as shown by Theorem 2 in [4]. After the consideration of the regular two-patch case in 2D we face now the coupling problem for SB-parametrizations.

Planar SB-IGA with C 1 coupling
Here we show that planar SB-IGA parametrizations are quasi analysis-suitable, except at the scaling center. This is important to obtain good convergence properties for the C 1 -coupled test functions. The problem with the singular point is addressed in the subsequent second subsection.

SB-parametrization as quasi AS-G 1
Analogous to the classical planar two-patch case we first look at two scaled boundary patches as displayed in Fig. 5. This means, we have constants c 1 > 0, c 2 ≥ 0 with two boundary curves γ (S) and a common scaling center x 0 and the patch parametrizations are of the form x0 Figure 5: Two-patch situation for SB-IGA.

Assumption 2.
Here and in the rest of the article we assume that the boundary curves γ (S) are parametrized in a strong C 1 sense and it is

Moreover, we assume that the boundary curves are chosen in such a way that for every
And for c 2 = 0 the parametrization is AS-G 1 except at the scaling center x 0 . In other words, the condition in Definition 1 is fulfilled.
Proof. We show the assertion for the interesting case c 2 = 0. If c 2 > 0 similar proof steps can be used. First let δ > 0 and ∂ ζ γ (R) (0) ∦ ∂ ζ γ (L) (0). Obviously, F |Ω δ is globally continuous. In view of (11) and Assumption 2 there are d 1 , By the orientation of the parametrization as elucidated in Fig. 6 it has to be d 1 = 0, d 2 = 0 and d 1 d 2 < 0. Namely, for the set It is easy to check that it has to be c > 0 otherwise we Figure 6: The orientation of the parametrization leads to d1 d2 < 0. Lemma 2 and Lemma 3 suggest that C 1 coupling for SB-IGA planar parametrizations leads to adequate approximation results of the underlying isogeometric spaces. In the paper [4] an important property of the AS-G 1 geometries is the large enough interface spaces of traces and transversal directional derivatives, see (27) and (28) in [4]. More precisely, for the SB-IGA two-patch case we look at the space where d is the transversal vector at the interface defined through (23) in [4].
Note that we can think of constant α (S) due to the proof of Lemma 3. Now we want to explain, why for SB-IGA with NURBS boundary curves the corresponding interface space still contains enough elements, at least away from the scaling center. The problem with the singularity is then faced in the subsequent section.
Proof. We follow the proof idea from [4]. First, let c 2 = 0 and assume a singular parametrization. By the properties of the NURBS, it is straightforward to see that there are functionsĉ (S) Obviously, the composed mappingĝ, i.e.ĝ | Ω (S) =ĝ (S) is continuous and except of a shift in the parameter ζ they are feasible NURBS basis functions in N r p ⊗ S r p . Moreover, one hasĝ (S) (ζ, ξ) ∈ O(ξ 2 ) and in view of the considerations in Sec. 4 . Basically, an analogous argumentation yields the second part of the assertion.
In view of the Theorem 1 in [4], the space V M,1 h seems appropriate for approximations away from the scaling center. From now on, we concentrate on the more interesting case with a singular parametrization, i.e. we assume that we have the scaling factor q(ξ) = ξ; cf. (5).

Approximation in the scaling center
Clearly, a special consideration of the behavior near the scaling center is needed. But a simple calculation shows that only such splines may cause problems at the scaling center which have nonvanishing values or derivatives in points (ζ, 0). For this purpose assume w.l.o.g.
. Now we check that the push-forward φ :=φ • F −1 m has a well-defined value and derivatives in x 0 . The former is obvious and we concentrate on the derivatives. The chain rule yields where γ = (γ 1 , γ 2 ). If we consider the derivatives away from the singular point, we get by assumption that the matrix on the right-hand side above is invertible, i.e.
Hence, it is By linearity and the definition of the B-spline basis functions we can suppose w.
For this case we study the derivatives when ξ → 0. With (13) one sees But this implies directly that φ(x, y) has a well-defined x-derivative in the scaling center, namely This means, it is justified to remove in each patch initially before patch coupling all the parametric basis functions N r i,p · B r j,p with j ≤ 2. But clearly, to preserve the approximation ability of SB-IGA test functions, we have to introduce new basis functions in the physical domain that determine the function value and derivatives at the scaling center. In the planar case, three additional test functions are sufficient, where we exploit the isoparametric paradigm to define preliminary test functions φ i,sc ∈ C 0 (Ω) with Latter requirements can be easily satisfied if we use the entries of the geometry control points as coefficients for the parametric pendantsφ i,sc . To be more precise, if we have Note that in the two-patch case we have i,sc . For example, the geometry in Fig. 2 (a) with control points in Fig. 2 (b) leads to the three scaling center functions in Fig. 7 below.
For a C 1 spline j c j B r j,p (ξ) the values in the first mesh interval are completely determined by the terms with j ≤ p + 1. Thus using the partition of unity property we get directly that φ (m) 1]} and hence φ 1,sc = 1, φ 2,sc = x, φ 3,sc = y in a neighborhood of x 0 . In other words, we can choose the latter three functions for the determination of values and derivatives at x 0 , i.e. we add them to the set of basis functions used for the coupling step. Note that the continuity of the global parametrization F implies the continuity of the composed φ i,sc . After we have seen how the scaling center issue can be handled, we have to face now the coupling conditions. This is content of the next section.

Enforcing the coupling conditions
We shortly mention an adaption of the approach of [4] to compute the actual C 1 -regular test functions. We want to emphasize that there are also other ways to enforce the coupling conditions, e.g. a least squares ansatz. Nevertheless, in view of clarity we restrict ourselves to one coupling procedure. Furthermore we explain the steps for a two-patch situation like Fig. 5.
We apply similar steps as in [4] to our context which can be summarized as follows. Initially, we start with the set of all the non-coupled basis functions Here we extend an uncoupled basis function toΩ by setting it to zero on the remaining patch. Further, below we writeB k for the set of parametric basis functions after the k-th step of the coupling procedure. The hat "· " is used if we consider the functions in the parametric domainΩ and we drop it if we mean the corresponding obvious push-forwards.
1. Remove all basis functions corresponding to N r i,p (· + 1) · B r j,p , N r i,p · B r j,p with j ≤ 2. We have now a new basisB 1 .
Latter step is done in order to remove the problematic test functions near the scaling center.
2. The remaining test functions inB 1 are coupled continuously which leads to a modified basis The C 0 coupling is easily achieved due to Assumption 1. This means we once more reduce again the number of basis functions. This means if one has a weak linear formulation of some problem, the assembly of matrices of the form (A) is a proper bilinear form and B 6 = {φ 6 1 , φ 6 2 , . . . , φ 6 N6 }, can be computed from the uncoupled system matrix ( Applying the above steps we can calculate C 1 -regular basis functions in the two-patch case. In case of a SB-parametrization which consists of more than two patches the coupling is done for each interface according to the mentioned approach. However, the definition of the additional scaling center test functions only has to be done once.
Remark 3. Later in the numerics part, we use Gauss quadrature rules to compute the matrix M J and other appearing integrals. In particular, we do not need to evaluate the basis functions at the singular point. Thus, the computation of the matrix M J is well-defined.

Remarks on generalizations
In the first part, the C 1 -coupling was explained for the two-patch case, but as already mentioned the approach can be generalized to situations with more patches. Hereto one enforces the C 1 coupling at each interface, but in case of c 2 = 0, i.e. classical SB-IGA, the scaling center test functions φ i,sc are defined once for the complete multipatch domain.

Non-star-shaped domains
Although we are now able to handle various geometries we are still limited to star-shaped domains. But frequently in applications the computational domain is not star-shaped. Nevertheless, in a special situation the C 1 -coupling can be generalized without loosing the (quasi) AS-G 1 structure. Namely, lets assume a decomposition of the domain Ω into star-shaped subdomains Ω m , where the interfaces between the different subdomains are straight lines; see for example Fig. 8. Then we know how to couple the patches within each subdomain and the new interfaces between the star-shaped subdomains again fit to the AS-G 1 framework since the two elements corresponding to that interface have a w.l.o.g. bilinear parametrization (see Fig. 8). And due to Proposition 3 in [4] we have that bilinear multipatch parametrizations are AS-G 1 . This property becomes evident in the numerical examples shown in the last part of the article.  Figure 8: A non-star-shaped domain which is divided into star-shaped subdomains that have noncurved interfaces. Such multipatch structures are still suitable for a C 1 coupling. If the interface between SB-parametrizations with two different scaling centers is a straight line, then w.l.o.g. the patches meet as two bilinear patches and the parametrization is quasi AS-G 1 , i.e. AS-G 1 except at the singular points.

Application to trimmed domains
Trimming, i.e. the cut off of domain parts utilizing trimming curves or surfaces, is a fundamental operation within CAD and is applied frequently. For more information and a detailed study we refer to [14]. But as explained in latter reference the implementation of general trimming procedures is quite complicated and for different approaches there arise several difficulties. Especially the stable integration over trimmed geometries might be an issue. SB-IGA with its boundary representation is suitable for the consideration of trimming and we want to explain here briefly, how we can handle C 1 -coupling if we incorporate trimming. The basic idea is to interpret the trimmed domain as a new computational domain which is defined by appropriate new boundary curves and to apply the coupling approach from above. First, we explain the procedure by means of a simple example and concerning more complicated situations we add some remarks later. Let the boundary curves γ (m) , m = 1, . . . , n of the untrimmed domain Ω be given; see the blue square boundary in Fig. 9 (a). Further, let a trimming curve γ T be given such that there are two intersections with the boundary ∂Ω, e.g. γ (1) (ζ (1) ) = γ T (s 1 ) and γ (2) (ζ (2) ) = γ T (s 2 ) for proper ζ (k) ∈ [0, 1]. Then, we have to check which curve segments belong to the boundary of the wanted trimmed domain Ω T . For the example in Fig. 9 (a) we assume that the curve segments γ T ,γ (1) ,γ (2) and γ (3) , γ (4) define the new boundary of Ω T . If we parameterize now the modified boundary partsγ T ,γ (1) ,γ (2) utilizing standard NURBS or Bsplines, i.e. e.g.γ T ∈ (N r p ) 2 we are again in the situation of classical untrimmed SB-IGA provided that we have a suitable scaling center for the trimmed domain; see Remark 4 below. The exact boundary representation of the trimmed domain is always possible which can be seen as follows. Let γ : [0, 1] → R 2 be a NURBS curve in (N r p ) 2 and 0 ≤ ζ (1) < ζ (2) ≤ 1. Inserting knots at ζ (1) , ζ (2) s.t. the multiplicity of the knots ζ (1) , ζ (2) is p + 1, we obtain NURBS with discontinuities at the mentioned two knots that represent the original curve γ; cf. [2]. Consequently, each new NURBS basis function is non-zero only in one of the three intervals (1) , ζ (2) ], I 3 = [ζ (2) , 0], see Fig. 10.
Hence, choosing the appropriate NURBS and after a simple re-scaling it is easy to see that we can describe each of the three curve segments γ |Ii with NURBS curves such that they start and end in points in {γ(0), γ(ζ (1) ), γ(ζ (2) ), γ(1)}. In other words, if the intersection parameter values in the trimming example above are known, we just have to apply some knot insertion steps to find exact parametrizationsγ (m) : [0, 1] → ∂Ω T of the trimmed domain boundary curves.
Clearly, the trimming example in Fig. 9 is very simple and in general situations some issues might occur. However, most problems can be handled as the next remarks indicate.
Remark 4. On the one hand, if we have the boundary curves of the trimmed domain, we need to compute a suitable scaling center, but the trimmed domain may loose its star-shape structure. Nevertheless, the trimmed domain can be partitioned into smaller star-shaped regions. We propose a naive ansatz, namely we divide the trimmed domain via straight cut lines into smaller star-shaped blocks in order to maintain the analysis-suitability. In more detail, we repeat trimming steps with straight lines but incorporate all the subdomains arising from the trimming with that straight lines. This simple partition approach is illustrated in Fig. 11. There we have the square domain with a trimming curve that leads to a non-star-shaped new domain. But the application of an division step with a proper straight line leads to a two star-domains for which scaling centers can be chosen, see Fig. 11 (b).
Remark 5. On the other hand, there could be more than 2 intersection points between trimming curve and original boundary. Then, the overall trimming still can be done in an iterative manner, i.e. we cut off successively simple regions. Also the situation of a trimming curve completely defined in the interior of the domain is manageable. One partitions the trimmed domain via straight lines into star-shaped blocks and handles each of this subdomains with SB-IGA. We refer to Fig. 12 for an illustration.
Remark 6. We assume always that the SB-mesh boundary is given by NURBS curves, i.e. curves of the form γ : [0, 1] → R 2 . Then, experiments not shown here indicate that if one uses exact boundary representations of the trimmed domain based on knot insertion together with the iso-parametric paradigm one might obtain SB meshes which are not optimal. To be more precise, one might get meshes with very thin elements which are not optimal in the sense of condition. Maybe if we only approximate the boundary curves with NURBS, e.g. with a L 2 projection ansatz, we might relax latter issue. γ T γ (1) γ (2) γ (3) γ (4) γ T γ (2)γ      (c) Possible SB-mesh after the trimming. Figure 12: Here we see a trimming curve that would lead to a non-star domain. We can divide the trimmed domain into two star-shaped parts using a straight cut line.

Numerical examples
In this section, the proposed method is investigated regarding its performance by means of several numerical experiments. There are different application examples and methods where C 1 regularity of test and ansatz functions is required for the numerical calculations. Here we want to study the application in structural mechanics of Kirchhoff plate theory for which we can exploit the increased inter-patch regularity ansatz from above. The approach is based on the assumptions made by Kirchhoff who states that "planes perpendicular to the mid surface will remain plane and perpendicular to the deformed mid surface".

Kirchhoff plate theory
Kirchhoff plates determine a two-dimensional fourth-order boundary value problem governed by the bi-Laplace operator. The description of the formulation below follows the notation and derivation presented in [17] and [18]. We consider domains Ω ⊂ R 2 that have a sufficient smooth boundary ∂Ω = Γ such that the unit normal vector n is well-defined. Further, the boundary is partitioned into parts for the specification of the energetically conjugate deflections and shear forces Γ = Γ u ∪ Γ Q and rotations and bending moments Γ = Γ φ ∪ Γ M . Further, we suppose that Γ u ∩ Γ Q = ∅ and Γ φ ∩ Γ M = ∅, respectively. The strong form of the Kirchhoff plate formulation is stated as with ∇(·) as the gradient operator, ∆(·) the Laplace operator and Ψ(·) as the third order differential operator Moreover, u is the deflection of the plate, g the load per unit area and u Γ , φ Γ , M Γ and Q Γ the prescribed deflections, rotations, bending moments and shear forces acting on the boundary. D denotes the bending stiffness consisting of the thickness t, the Poisson ratio ν and the Young modulus E. As t, ν and E are assumed to be constant over the domain, yielding an isotropic, homogeneous material, the bending stiffness is defined as The classical primal weak form for the computation of approximate solutions has the form for a suitable test function space V h that depends on the boundary conditions. Moreover, a(u h , v h ) is a bilinear form defined as Furthermore, the linear functional F is denoted as We remark that in the numerical tests no external shear forces Q Γ and bending moments M Γ are applied except for the L-bracket 6.5. Moreover, the boundary conditions for u Γ (18) and φ Γ (19) need to be enforced strongly, while the conditions of Q Γ (20) and M Γ (21) are natural boundary conditions. Provided proper boundary conditions and suitable V h ⊂ H 2 (Ω), the conditions of the Lax-Milgram theorem are satisfied and the weak form has a unique solution u h ∈ V h . In the examples below, the discretization of the solution field is performed by the previously defined SB-IGA test functions including the coupling approach at the patch boundaries, which ensure C 1continuity within the whole domain that is required for the plate formulation. In other words, the discrete spaces are of the form V h ⊂ V M,1 h ; see (10), since boundary conditions are taken into account. In the following examples, we mean mesh size h if the underlying parametric meshes for each patch have 1/h equidistant subdivisions with respect to both parametric coordinate directions.

Smooth solution on a square plate
At first, the plate formulation is checked on its general performance. A square plate of Ω = [−0.5, 0.5] 2 is subjected to a smoothly distributed source function g that is chosen such that the exact solution is u = cos(πx) 2 cos(πy) 2 . A plate of thickness t ≈ 0.1063 is considered of elastic material with Young's modulus E = 10 4 and Poisson's ratio ν = 0.0, such that the flexural stiffness D = 1. We have clamped boundary conditions and the domain is discretized with four SB-IGA patches. The scaling center is placed with an offset of the plate center, namely x of f = −0.15 and y of f = 0.1, to demonstrate its approximation behaviour for non-symmetric meshes. Fig. 13 exemplary shows the mesh for h = 1/4 and the corresponding deformation plot.  The numerical solution is evaluated with respect to the H 2 seminorm and L 2 norm that are computed for third, fourth and fifth order basis functions, where the error norms are defined as For the corresponding convergence rates, we refer to  The convergence rate indicate for both error estimates optimal convergence rates O(h p−1 ) for the H 2 seminorm and O(h p+1 ) for the L 2 norm. Especially we do not see a C 1 locking effect.

Point load on a square plate
The next example shows the capability of the proposed formulation to consider point loads even in the scaling center. Again, the square plate defined as Ω = [−0.5, 0.5] 2 is considered and the boundaries ∂Ω are simply supported. For the material parameters we choose E = 10 6 and ν = 0 and the thickness is chosen such that D = 1. The point load is defined as F = 1. There is an analytical solution for the displacement of the plate center under a point load, namely using [19] one gets: where L is the length of plate. The reference solution u ref ≈ u ex for the deflection in the center under the load application is obtained from the series above by taking sufficient terms into account. Fig. 15 shows a mesh exemplary and the corresponding deformation plot. Besides we look at the relative errors calculated by |1 − u/u ref |, where the convergence study is shown in Fig. 16. Here, the example points out that the formulation is capable to determine results properly even for loading applied in the singular point. We think that the deviation from the best approximation rates is caused by the lack of regularity for the point load, which is not a function anymore, strictly speaking. In the mentioned figure we display also the bending moment m 11 = D ∂ xx u + ν ∂ yy u, a quantity for which the C 1 regularity is crucial.

Perforated circular disk
Having tested the method for the plate formulation in a simply square domain, trimmed geometries are now incorporated. Therefore, a simply supported disk of radius R = 1 and t = 0.02, E = 10 7 , ν = 0.3 is trimmed by four holes of diameter d hole = 0.1, see Fig. 17 (a). The center of the holes are placed on the x and y axes with a distance of 0.4 from the origin. The disk is loaded by a constant surface load g = 1. The reference solution is obtained from [12] as the maximal deflection in the center u ref = 0.008950, which is a numerically converged solution determined using the commercial software Abaqus with the triangular, linear shell element S3. Since the SB-IGA formulation requires star shaped patches, it is not possible to determine the deflection by a single scaling center but several scaling centers are required; see Sec. 5.1. In the following, a mesh of 25 scaling centers is evaluated (see Fig. 17 (b)) with varying numbers of elements and orders on each patch.  Fig. 17. Furthermore, the convergence is evaluated in 18 (b) for the orders of p = 3, 1/h = 2, 3, 4, 5 and p = 4, 1/h = 2, 3, 4. For a constant surface load nearly exact results are obtained by at least cubic basis functions. This is fulfilled for both orders. A comparison to an automatic meshing approach of triangular Bézier spline elements presented in [12] shows that the incorporation of trimming in the structural description of the model as in the context of SB-IGA has advantages over mesh finding algorithms as the meshing algorithm does not converge towards the reference solution even for very fine meshes.

L-bracket
In the second last example, an L-bracket is investigated analogously to [20] and [21] with the geometry shown in Fig. 19. The bracket has underlying parameters E = 200 · 10 9 , ν = 0 and t = 0.01. The boundary conditions are clamped conditions on the left holes (upper and lower hole) as well as a constant line load of f = 100 on the upper right edge marked in blue. Further, the partitioning of the mesh created with 20 scaling centers is given in Fig. 19 (b). We emphasize that meshes with less scaling centers could be applied, however, the model is not chosen to demonstrate the convergence, but the general application to complex models. p=4 [12] (b) u/u ref and r = 1 is shown. On the right side, the convergence studies of the displacement in the middle of the perforated disk with orders of p = 3, p = 4 (r = 1) is compared to [12].

A clamped violin
In this example, we demonstrate the applicability of the coupling procedure in the context of more complicated geometries involving non-trivial trimmed parts. In more detail, we look at a KL-plate which has the shape given in Fig. 22 (a). This geometry is inspired by the violin example in [22], however, here the setting as well as the geometry differs. We note that the boundaries of Fig. 22 (a) are represented by degree 3 B-splines. A decomposition of the domain into star-shaped blocks with straight interfaces is possible and for the computations we use the mesh structure given in Fig.  22 (b). We consider the case of a uniform load function f = −0.01 and as boundary conditions we require the clamping of the inner boundaries that define the f-holes. All the other boundaries can move freely. Further we choose the material parameters E = 10 5 , ν = 0.1, t = 0.2. Using Bsplines of degree 3 and regularity 1 for the determination of the coupled basis functions we obtain as deformation the result in Fig. 23. Although a reference solution is missing, we regard the numerical deformations as reasonable and the coupling ansatz seems to work for this advanced setting.

Remarks on stability
The examples from above confirm that we can apply the C 1 -coupled basis functions in the context of the Kirchhoff plate. In particular, it is still possible to handle the high-order problem although we have a singular parametrization mapping. Certainly, this singularity comes along with some difficulties. The main bottleneck is the occurrence of large condition numbers for the underlying linear systems, especially if the mesh sizes gets small. The latter causes problems in the sense of stability of the approximate solution and below we discuss this briefly. First, we bring to mind that the definition of the C 1 basis functions implies the boundedeness of the H 2 -norm, i.e. V M,1 h ⊂ H 2 (Ω) if we have B-spline ansatz spaces and B-spline parametrizations. By construction, we have the smoothness of each basis function in the interior of each mesh element and as already mentioned the global continuity of the first derivatives. Then, if we look at 'Case I' from in Sec. 4. of [23], we get with Theorem 4.4. from latter article that the coupled SB-IGA test functions are in H 2 (Ω m ) for each patch. To see this, we observe that we use for the coupling method only such B-splines B r i,p · B r j,p with j > 2 or the additional scaling center test functions, which are obviously smooth in a whole neighborhood of the scaling center. But, as piece-wise H 2 mappings which are globally C 1 we also get the H 2 property in the whole domain. Since we use a quadrature formula with discrete evaluation points away from the scaling center, we obtain well-defined integral values and one does not see directly the singularity within the computations of system matrix entries. However, we can observe large second derivative entries caused by the mesh degeneration. This comes along with instability problems for the discrete problem in the case of fine meshes. In other words, the condition number increases strongly. This instability issue can be seen in figures Fig. 25 (a)-(b) and Fig. 26 (a). Here we have to note that we use for all examples the standard mldivide() function from MATLAB [15] to solve appearing linear systems. A naive idea to relax the problem of evaluations near the singular point is to reduce the number of basis functions near the scaling center. An easy way without much additional effort, is to combine the basis functions of two meshes. One coarse mesh, for the basis functions with support near the scaling center and one refined mesh, for the areas away from the singular point. In more detail, one could proceed as follows. We have given two SB-meshes for some domain. One fine mesh with subdivisions w.r.t. both parametric coordinates (see Fig. 24 (a)) and secondly we have the analogous SB-mesh of the same domain, but which is only refined in the radial direction, compare Fig. 24. Clearly, the coarse mesh is not suitable to highly accurate approximations. However, the mesh elements adjacent to the scaling center have diameter O(h) and the mesh from Fig. 24 (b) is enough for an approximation in the neighborhood of x 0 . Hence, in view of Fig. 24(a) and Fig. 24 (b) we consider only those basis functions in the fine mesh which have a support outside A which stands for the union of mesh elements adjacent to the singular point. This means we consider Consequently, our modified (uncoupled) SB-IGA test function space is After this we can go through similar steps like in the sections before to define C 1 -regular global basis functions. This idea can straightforwardly generalized to more complex geometries. Although such a combination of two meshes is certainly not the only method to remove problematic basis functions near the singularity, it is easy to implement and it is helpful to demonstrate the stabilizing influence of a such mesh coarsening near the center. Therefore, we repeat below the convergence test from Sec. 6.2 as well as the example of the plate with point load from Sec. 6.3, but now with the above explained two-mesh ansatz. If we look at the results in Fig. 14 and 16, we see an obviously more stable decay behavior compared to the original unstabilized computations.
These first experiments with modified spaces reveal the importance of a good mesh choice, at least when dealing with high order problems. For sure, different refinement strategies should be discussed, too, maybe also in combination with some adaptive scheme. And the application of the naive two mesh ansatz in the context of more complicated test cases is advisable. But, it is not our goal of this article to study the stability issue in detail. This seems a reasonable object of investigations with another publication. We want to conclude this subsection with the following summarizing words. The SB-IGA space with C 1 coupled basis functions suffers from bad conditions numbers, especially if we look at high order problems combined with fine meshes. Nevertheless, there is hope to alleviate this problem. We think that removing basis functions near the scaling center leads to increased stability in the computations without reducing the overall convergence behaviour significantly.

Conclusion
In this contribution, an approach of the incorporation of SB-IGA patches and C 1 coupling in terms of analysis suitable G 1 parametrization was presented. This implies a special consideration of the basis functions at the scaling center. The method was tested in the context of the Kirchhoff plate formulation. It is especially suitable for trimmed models since the boundary representation can easily replace the existing boundary by the trimming curve, no matter if the trimming curve is entirely inside the domain or partially outside. The proposed approach was tested in various ways such as the L 2 norm and the H 2 seminorm against the optimal convergence rates, the bending moments and the deflection at specific points on the accuracy per degree. The stability issues arising at the scaling center for fine meshes are discussed and a stabilization remedy is presented that yields improved results. Moreover, the results were compared to other approaches considering meshing for unfavorably shaped problems and to results from other coupling approaches in isogeometric analysis.
In conclusion, we presented the feasibility of SB-IGA with C 1 coupling in terms of Kirchhoff plates that is especially powerful for sophisticated structures.