Smooth multi-patch scaled boundary isogeometric analysis for Kirchhoff–Love shells

In this work, a linear Kirchhoff–Love shell formulation in the framework of scaled boundary isogeometric analysis is presented that aims to provide a simple approach to trimming for NURBS-based shell analysis. To obtain a global C1\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$C^1$$\end{document}-regular test function space for the shell discretization, an inter-patch coupling is applied with adjusted basis functions in the vicinity of the scaling center to ensure the approximation ability. Doing so, the scaled boundary geometries are related to the concept of analysis-suitable G1\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$G^1$$\end{document} parametrizations. This yields a coupling of patch boundaries in a strong sense that is restricted to G1\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$G^1$$\end{document}-smooth surfaces. The proposed approach is advantageous to trimmed geometries due to the incorporation of the trimming curve in the boundary representation that provides an exact representation in the planar domain. Since the coupling of star-shaped blocks is feasible, a mesh generation also for complex domains is implementable. The potential of the approach is demonstrated by several problems of untrimmed and trimmed geometries of Kirchhoff–Love shell analysis evaluated against error norms and displacements. Lastly, the applicability is highlighted in the analysis of a violin structure including arbitrarily shaped patches.


Introduction
In modern engineering applications such as automotive engineering and aerospace engineering, computer-aided design (CAD) is a key method to design structural components.The mathematical underlying of these drawings usually consists of B-splines and non-uniform rational B-splines (NURBS) and provides that by using isogeometric analysis (IGA) [1] for the numerical analysis, the model for the design and the analysis are the same.A major issue of numerical analyses is the treatment of trimmed geometries whereas parts of the initial geometry are cut out [2].This implies, that a numerical analysis in common frameworks is not conductible anymore in a straightforward manner.Several approaches are presented in the literature for treating trimming such as [3][4][5] mostly considering an approximation of the trimming curve and a limitation to quadrilateral sections.
Isogeometric analysis is especially suitable for Kirchhoff-Love shells [6] since the formulation is derived on the change of the normal vector which requires well-defined second order derivatives of the basis functions.This is naturally fulfilled within single patches of quadratic or higher order.However, across patch boundaries, IGA patches are not naturally C p−1 continuous.This implies that the Kirchhoff-Love shell formulation is not applicable for multi-patch structures in IGA.As complex shell structures often consist of multi-patch geometries, C 1 continuity needs to be enforced at patch boundaries.Common approaches to ensure higher continuity on multiple patches are Nitsche methods [7,8], mortar methods [9,10], and penalty methods [11][12][13][14], among others [15,16].However, all of these methods are coupling techniques in a weak sense that can only represent the isogeometric discretization space approximately C 1 smooth.To ensure exact C 1 smooth multi-patch geometries, a special class of G 1 continuous multi-patch surfaces, namely analysis-suitable G 1 surfaces [17] can be utilized in the two-dimensional domain.An approach of analysis-suitable G 1 spaces that will be utilized in this work, is presented in [18] for Kirchhoff-Love shells with restriction to quadrilateral patches.
A promising approach to take the exact boundary of geometries more into account is the scaled boundary isogeometric analysis (SB-IGA) [19,20].It combines the advantages of the scaled boundary finite element method [21,22] and the isogeometric analysis providing patches of arbitrarily shaped boundaries.Since the discretization of SB-IGA patches does not only consist of multipatches for trimmed patches but also in a single IGA-patch domain, a coupling approach is inherently necessary to obtain C 1 -continuity on the SB-IGA domain.
This contribution presents a Kirchhoff-Love shell in the framework of scaled boundary isogeometric analysis with analysis-suitable G 1 surfaces that are especially suitable for complex, trimmed geometries by incorporating the trimming curve into the boundary representation of the geometry in a simple way.Moreover, to preserve the approximation ability of the SB-IGA test functions, special scaling center test functions have to be introduced.A step that can also be observed if smooth polar splines are considered [23].The proposed approach is tested on its general performance by untrimmed shell structures and the potential is later outlined by trimmed shell structures in the linear case.The advantages are summarized as follows.
• A linear Kirchhoff-Love shell formulation is presented in isogeometric boundary representation.
• A patch coupling of smooth patches across boundary edges is applied utilizing the analysissuitable G 1 concept including modifications in the vicinity of the scaling center.
• Advantages of the approach to discretized star-shaped domains of an arbitrary number of boundaries are outlined.
• Examples of discretizations for complex topologies are presented and evaluated.
The work is structured as follows.In Section 2, the linear formulation of the Kirchhoff-Love theory is provided.Section 3 presents the fundamental concept of SB-IGA and the notation herein.Further, Section 4 discusses the choice of basis functions for C 1 continuity in the planar domain and the application to Kirchhoff-Love shells with the incorporation of trimming.The proposed approach is tested in numerical examples in Section 5. Finally, the contribution is concluded in Section 6 including an outlook to further research.

Kirchhoff-Love shell formulation
The Kirchhoff-Love shell formulation is recalled in this section, based on the derivations presented in [6,8,24].The formulation is derived considering a single patch but can be extended to scaled boundary multi-patch structures as outlined in the following sections.Neglecting cross-sectional shear stresses, the Kirchhoff-Love shell formulation states the assumption that the initial director vectors to the shell center surface remain perpendicular during deformation.The description of the shell is reduced from the domain of the shell body to the mid-surface Ω ⊂ R 3 utilizing a convective covariant space Ω ⊂ R 2 .Utilizing the unit normal vector on the mid-surface, the initial shell structure reduced on the mid-surface is described as where a is the position vector in the initial configuration, a 3 is the unit normal vector on each point on the surface, and θ 3 the thickness coordinate ranging from −t/2 to +t/2.Further, we use the Greek indices α, β, λ, µ ∈ {1, 2} and Latin indices i, j ∈ {1, 2, 3} and drop the dependence of θ α .To obtain the normal vector on the surface, the partial derivatives on the mid-surface are utilized to construct the unit normal vectors with In the following, we will indicate differentiation w.r.t. the corresponding coordinate (•) ,α = ∂ ∂θ α .Further, we define the L 2 -norm of the normal vector a 3 as J. Utilizing the partial derivatives, we introduce the covariant metric coefficients constructed by ( Due to the orthogonality condition, we construct the contravariant basis system a α • a β = δ α β , that yields a relation of the covariant and the contravariant metric system by the inverse operator {a αβ } = {a αβ } −1 , where {•} denotes the matrix components.In the discrete setting of the boundary conditions of the boundary value problem, we denote the boundary of the body Γ = ∂Ω with Dirichlet and Neumann boundary conditions Γ D and Γ N with the assumption that Ω is a smooth surface with well-defined derivatives of the curvature.The boundary conditions of Γ D and Γ N are non-overlapping Γ = Γ D ∪ Γ N .Further, we decompose the boundary conditions into since the prescribed displacements ũ and traction forces Q, as well as the prescribed rotations β and traction moments M at the boundary are of energetically conjugate nature.A similar consideration is done at the corners χ ⊂ Γ with a split into Neumann boundary conditions χ N ∩ Γ N and Dirichlet boundary conditions χ D ∩ Γ D .Having defined the boundary, we choose some suitable discrete space V h ∈ H 2 (Ω) in the sense of finite element methods and depending on the boundary conditions of the problem.Then, the discrete weak form seeks for Where the bilinear form a and the linear form F are expanded as where g is an applied body force, S a prescribed twisting moment at each corner of χ N , and β n the normal rotation [8].In this work, only body forces and traction forces are considered.Additionally, for a complete weak form, the kinematics and the constitutive relation are pending.The quantities ε, κ, n and m denote the membrane strains and bending strains depending on the test function v h and the energetically conjugate membrane forces and bending moments depending on the discretized solution field u h .According to the Kirchhoff kinematical assumptions, transverse shear strains vanish and membrane and bending strains remain.Consequently, the decomposition of the linearized Green-Lagrange strain tensor reads For a rigorous derivation of the linearization, see [8].According to [25] the linearized strain components are and Besides the kinematics, the stress resultants n and m are derived by applying a suitable constitutive relation.Thus, the second Piola-Kirchhoff stress tensor is analytically integrated through the thickness and reads as a decomposition into the stress tensors' components as [13] n and with the materials fourth-order tensor D. We want to remark, that the analytical integration herein assumes a constant thickness of the shell.Applying Hooke's law as a linear constitutive model, the tensor is given as [13] and Herein the parameters ν and E are the specific material parameters of Poisson's ratio and Young's modulus, respectively.It is remarkable, that due to the occurrence of second order derivatives of the basis functions in the bending components, C 1 continuity is required for the computation.This is fulfilled naturally for SB-IGA patches across the elements within each patch, however, at the patch boundaries, only C 0 continuity is given.Thus, C 1 continuity is enforced across patch boundaries in the manner of scaled boundary isogeometric analysis.

B-splines and SB-IGA
In this section, we introduce the notation herein and explain briefly some basic notions in the context of the SB-IGA ansatz.For more details regarding isogeometric analysis, we refer to [26,27].
We start with the definition of B-spline functions and B-spline spaces, respectively.An increasing sequence of real numbers Ξ := {ξ 1 ≤ ξ 2 ≤ • • • ≤ ξ n+p+1 } for some p ∈ N is called knot vector, where we assume 0 , 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 Note one sets 0/0 = 0 to obtain well-definedness.The knot vector Ξ without knot repetitions is denoted by {ψ 1 , . . ., ψ N }.
The extension of the last spline definition to the multivariate case is achieved by a tensor product construction.In other words, we set for a given tensor knot vector 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 global 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 lsmooth 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.To emphasize later the regularity of the splines, we introduce an upper index r and write in the following B r i,p , B r i,p respectively.Here r := (r 1 , . . ., r d ).
• For univariate splines B r i,p , p ≥ 1, r ≥ 0 we have • The support of the spline B r i,p is part of the interval [ξ i , ξ 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 Later, to have more flexibility it could be useful to introduce a strictly positive weight function , the NURBS spaces N r p := 1 W S r p , respectively.Such weighted spline functions are needed for conic section parametrizations.For the multivariate case we just define the needed spaces as product spaces, e.g.
provided proper univariate spaces.
To define discrete spaces in the computational domain Ω we require a parametrization mapping F : Ω := (0, 1) d → Ω ⊂ R d .The knots stored in the knot vector Ξ, corresponding to the underlying splines, determine a mesh in the parametric domain Ω, namely as the knot vector Ξ without knot repetitions.The image of this mesh under the mapping F, i.e.M := {F(K) | K ∈ M }, gives us a mesh structure in the physical domain.Obviously, by inserting knots without changing the parametrization we can refine the mesh, which is known as h-refinement; see [1,26].For a mesh M we can define 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 concept of SB-IGA fits the fact that in CAD applications the computational domain is often described by means of its boundary.Furthermore, as discussed later the boundary representation is particularly suitable if trimming is considered.As long as the region of interest is star-convex we can choose a scaling center z 0 ∈ R d and the domain is then defined by a scaling of the boundary w.r.t. to z 0 .In this article, we restrict ourselves to planar SB domains, and in view of isogeometric analysis we have some boundary NURBS curve and define the SB-parametrization of some Ω through Remark 1.In the considerations below it is allowed to replace the prefactor ξ by a general degree 1 polynomial c 1 ξ + c 2 , c 1 > 0, c 2 ≥ 0, i.e. no difficulties arise.This might be useful in some situations.
By the linearity w.r.t. the second parameter ξ we can assume for In particular, the weight function depends only on ζ.We suppose that there are control points For simplification, we assume 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 The isogeometric spaces utilized for discretization methods are defined as the push-forwards of the NURBS, namely In particular, we suppose the existence of an inverse mapping on the interior of Ω.If the domain boundary ∂Ω is composed of different curves γ (k) , one defines parametrizations for each curve as written above and we get a multi-patch geometry; see Fig. 3. To be more precise, for a n-patch γ (1)  γ ( 2) and F m are SB parametrizations.IGA spaces in the multi-patch framework are straight-forwardly defined as where denotes the IGA space corresponding to the m-th patch, to F m , respectively.If the boundary curves do not meet G 1 but high global regularity is necessary, then the coupling is an issue.For all the patch coupling considerations we suppose the next assumption.

Assumption 1 (Regular patch coupling). sdakjhjkh
• In each patch we use the NURBS and B-splines with the same degree p and regularity r > 0.
Furthermore, the F m are C 1 in the interior of each patch with C 1 -regular inverse.
• The control points at interfaces match, i.e. the control points of the meeting patches coincide along the interface.
Thus, it is justified to write for the set of parametric basis functions in the m-patch Below we look at the C 1 coupling of such SB-IGA patches, i.e. face spaces of the form where the singularity of the F m at z 0 requires a special attention.

From planar SB-IGA to KL shells
The appearance of second derivatives in the variational formulation (4b) for the Kirchhoff-Love shell model requires a proper choice of the test and ansatz spaces for the computation of the shell displacements.An important aspect is the representation of the shell via a middle surface which reduces the shell formulation to a 2D problem w.r.t. the parametric coordinates.We exploit this fact by using C 1 -regular mappings in the parametric domain to express the shell configuration.
R To be more precise, the basic idea for the KL-shell within SB-IGA framework is as follows.We express the initial shell configuration as well as the deformation of the shell utilizing the coupled test functions defined on the parametric domain Ω, where we suppose Ω = ∪ m Ωm to be given as a multi-patch SB-IGA geometry.On the patches Ωm in turn we can introduce basis functions by means of the standard push-forwards.This means our shell middle-surface can be described by for suitable control points P (m) i,j ∈ R 3 ; see Fig. 4. The advantage of such an approach is clear, namely coupling in the planar domain Ω is easier than a strong C 1 coupling of NURBS surface patches in 3D.And the restriction to the special class of SB domains Ω gets useful if trimmed shells are considered.Assumption 2. The initial shell configuration is C 1 smooth in a strong sense.Further, we assume that the shell deformations can be expressed with C 1 mappings.In particular, we do not consider shells with kinks or deformations that lead to non-smooth shell configurations.
Summarising it is enough to introduce globally C 1 basis functions on Ω to obtain test and ansatz functions for the Kirchhoff-Love shell model.The coupling of scalar SB-IGA spaces in planar domains is part of the next section.A componentwise generalization to compute shell displacements is clear and will not be addressed separately.

Planar SB-IGA with C 1 coupling
We explain the procedure of how to get C 1 basis functions in the two-patch situation as shown in Fig. 5 since a generalization to more patches is straightforward.According to the notation from the previous section and in view of isogeometric analysis we define the uncoupled basis functions on Ω = Ω1 ∪ Ω2 Above we extend the basis functions to the whole Ω by setting them to zero on the remaining patches.To obtain the necessary global C 1 smoothness we change the basis B according to the subsequent steps.

Remove basis functions near the scaling center
First, we remove in each patch all parametric and thus physical basis functions which are associated to N r i,p • B r j,p with j < 3.This means we consider the modified basis As a consequence, we see directly that the pull-backs of the remaining basis functions are elements of C 0 ( Ω), i.e. we have a well-defined value 0 in the singular point z 0 .Further, it can be easily verified that the pullbacks onto the parametric domain define mappings C 1 ( Ω).Indeed, an application of the chain rule implies vanishing derivatives of the physical basis functions in the scaling center.For this purpose assume w.l.o.g.
Below we use the abbreviation and we have d(ζ) = 0 due to the assumed invertibility of F. Hence, it is By linearity and the definition of the B-spline basis functions and the assumption For this case we study the derivatives when ξ → 0. With ( 16) one sees But this implies directly that φ(θ 1 , θ 2 ) has a well-defined θ 1 -derivative in the scaling center, namely ∂ θ 1 φ = 0 in z 0 .Analogously one gets the well-defined derivative ∂ θ 2 φ(z 0 ) = 0.With this first basis modification step we avoid the problematic part near the scaling center.

Add scaling center basis functions
To preserve the approximation ability of SB-IGA test functions, we certainly have to introduce new test 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 iso-parametric 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 And then the φ i,sc are determined on Ω via By the properties of B-splines we see directly φ 3,sc = θ 2 in a neighborhood of z 0 .In other words, we can choose the latter three functions for the determination of values and derivatives at z 0 , i.e. we add them to the set of basis functions used for the coupling step.Hence now we have the new set of basis functions The three new test functions are continuous since we assume that the patches match continuously.And we note that the φ i,sc are only defined once for a scaling center.For example, in Fig. 6, the scaling center test functions for a three-patch case are shown.

The coupling step
For the functions in B , we now apply the needed coupling step to obtain the global C 1 regularity.This coupling is in turn composed of two simple substeps which correspond to the procedure used in [17].Namely, we first couple the basis functions from B in a C 0 manner and observe that this is easily achieved due to the regular patch coupling Assumption 1.So we have again a new basis B = B ∩ C 0 ( Ω).Looking at the interface Γ between two patches, e.g.like in Fig. 5, we have for φ ∈ B a well-defined directional derivative in the direction of the interface.Thus it is enough to find linear combinations φ of basis functions in B that have continuous directional derivatives ∇φ • n Γ , where n Γ is a fixed normal vector to the interface.Consequently, we get the C 1 coupled basis functions by the null space of the derivative jump matrix M Γ , where with [[g]] standing for the jump value of the g across the interface and above we assume B := {φ 1 , φ 2 , . . .}. Straightforwardly this null space computation is adapted if several interfaces are involved.
Finally, we have with the basis of the null space of M Γ our C 1 smooth basis functions, i.e. our desired coupled basis B 1 := B ∩ C 1 ( Ω).We denote the space spanned by the coupled basis functions as Remark 2. The steps to get smooth basis functions above can be implemented directly.Even if the computation of the jump matrix for the different matrices seems problematic, we get welldefined integral values since the evaluations are only done in quadrature points away from the singular point.Besides, we note that all the steps from above reduce to a suitable transformation matrix T ∈ R N ×N1 expressing the coupled basis functions w.r.t. the original basis.In particular After we showed the possibility to obtain C 1 mappings in the planar SB context, an obvious question arises.Do the coupled functions yield a reasonable approximation behavior if we utilize them in numerical applications?It is well-known that a strong C 1 coupling might lead to C 1 locking, i.e. the loss or worsening of approximation abilities.Admittedly, we do not give here a strict proof, of whether the C 1 coupling for planar SB domains ensures optimal convergence orders.However, we want to outline why we expect the latter.Numerical tests confirm this conjecture.Following the results from [17], we get optimal convergence orders in the planar multi-patch case, if we restrict ourselves to so-called analysis-suitable G 1 (AS-G 1 )parametrizations and B-spline basis functions with p > r + 1 > 1.The mentioned parametrizations class is defined for the standard two-patch case, cf.Fig. 7, as follows.
Definition 1.In view of Fig. 7 we set The global parametrization F, i.e.F| Ω (S) = F (S) , S ∈ {L, R}, is analysis-suitable G 1 if there are polynomial functions We call then the complete multi-patch parametrization analysis-suitable G 1 , if the latter condition is fulfilled for each interface.With the next lemma we relate the introduced AS-G 1 geometries to the planar SB ansatz.
Figure 7: A standard planar two-patch geometry.
Proof.This can be shown by simple calculations.But it is also clear by the observation that the SBparametrizations can be approximated in a neighbourhood of an interface by bilinear parametrizations.This can can be exemplified if w.l.o.g.F(ζ, ξ) = ξγ(ζ), where the interface corresponds to the points with ζ = 0. Then it is But since bilinear parametrizations are AS-G 1 as shown in [17], the assumption is clear.
Hence, if we are away from the singular point we get with the SB ansatz AS-G 1 parametrizations.And for the approximation near the scaling center we added the scaling center basis functions φ i,sc .This and Theorem 1 in [17] suggest the feasibility of SB-IGA for C 1 coupling.

Generalization to non-star shaped and trimmed geometries
In this section, we shortly explain why the coupling from above can be generalized to more complicated geometries.Up to now, the computational domain is considered to be star-shaped.However, assuming we can partition our domain into several star-shaped blocks, we can do the coupling as explained previously block-wise.The coupling between the different blocks can then be realized analogously and easily if the block interfaces are given by straight lines.Moreover, two patches from two blocks that meet at a straight interface meet w.l.o.g. as two bilinear parametrizations, see Fig. 8.But as already mentioned, the bilinear parametrizations are AS-G 1 meaning we should not see The mentioned natural approach to applying a preliminary partition step to handle more complicated geometries is convenient when considering trimmed planar domains.The trimming operation, i.e. the cut away of domain parts is a standard operation within the IGA community and standing to reason for the description of different geometries.Especially, when dealing with perforated domains trimming curves are adequate to determine the computational domain.The concept of trimming can be formally incorporated within planar SB-IGA.Namely, let us suppose an untrimmed domain Ω is defined as SB multi-patch domain.If a trimming curve γ T defines which parts are omitted, we first compute the intersections between the trimming curve and the boundary curves of the Γ z 0 z 1 Figure 8: A non-star-shaped domain which is divided into star-shaped subdomains that have noncurved interfaces.Such multi-patch structures are still suitable for a C 1 coupling.If the interface between SB-param.with two different scaling centers is a straight line, then w.l.o.g. the patches meet as two bilinear patches and the param.is quasi AS-G 1 , i.e.AS-G 1 except at the singular points.untrimmed domain.First, we think of a non-interior curve, i.e. there are at least two intersections.After possible knot insertions, it is possible to represent all relevant boundary segments by new NURBS curves exactly.If the trimmed domain is not-star shaped we first partition the new domain into star-convex blocks and then choose suitable scaling centers.Otherwise, we choose directly a feasible scaling center.If there is an interior trimming curve, we start with the partition step and then we detect the relevant boundary curves that are again inputs for the blocks-wise SB-parametrizations.Since the situation with an interior trimming curve is more interesting for us and the only one appearing in the numerics section later, we illustrate the proceeding for the latter case in the subsequent Fig. 9.

Algorithm 1 SB parametrization and trimming
Require: Trimming curve γ T ; Untrimmed domain boundary curves γ (i)  1: Compute intersection parameters if there are no intersections then Interior trimming curve

3:
Divide trimmed domain into star-shaped blocks

Go to line 3 9: end if
After we discussed how we can obtain C 1 basis functions on planar SB multi-patch domains, we exploit them below to compute approximations to the Kirchhoff-Love shell model.

Numerical Examples
In the following, the presented formulation is checked on its consistency at first by a simple test as square flat shell.In the further, the example of the Scordelis-Lo roof, known from the so-called 'shell obstacle course' is evaluated [28][29][30].Moreover, the double curved hyperbolic paraboloid [31] is evaluated as it considers negative Gaussian curvature and is a highly challenging example.Besides, the example of the Scordelis-Lo roof and further examples are evaluated in terms of trimmed structures as they can be found in literature [3,5,32] and are especially suitable for the boundary representation.Conducting an h-refinement for all examples, the models are validated by error norms or the displacements at specific points on reference solutions obtained analytically or from the literature.We denote that we further use t as the thickness, E as the Young's modulus, ν as the Poisson's ratio, and h as mesh size for the underlying mesh with 1/h equidistant subdivisions with respect to both parametric coordinate directions.The approach has been implemented using MATLAB 2022 [33] in combination with the open source package GeoPDEs [34].It is specifically designed to solve partial differential equations in the context of isogeometric analysis.We use the model from above, namely basis function in the space V M,1 h for each component.But we note that due to the specific boundary conditions we have to reduce the number of test functions at specific places.

Geometry approximation
The coupling of the basis functions for the Kirchhoff-Love shell is done in the parametric domain Ω and R denotes the initial shell configuration parametrization.Sometimes the latter is not directly available in the form (14).Then, due to the special structure of SB-IGA, it is not clear how to choose in general the control points in order to represent the initial shell configuration in an exact manner, to obtain R in the sense of (14), respectively.Thus we exploit the coupled C 1 -smooth basis functions in the parametric domain to approximate the initial shell shape for the examples below.In this sense, we follow the isogeometric paradigm.To be more precise, if we state below the geometry mapping R, we use in the actual computations an approximated R approx ≈ R. Doing so, we use a L 2 -projection onto the test function space to get R approx .Although this is a naive geometry approximation, it gives us even for coarse meshes reasonable results.Furthermore, several geometries like the hyperbolic paraboloid shape can be still represented exactly, except of rounding errors; see Fig. 10 below.To be more precise, the Euclidean norm difference |R approx (θ 1 , θ 2 ) − R(θ 1 , θ 2 )| is plotted, where the underlying meshes are displayed in Fig. 25 and Fig. 19.Even for relatively coarse meshes we see only small differences and the geometry errors can be assumed to be small.Furthermore, in principle, it is possible to use different NURBS for the initial geometry in order to have exact geometries.But, this requires more implementation effort due to the special form of the SB-IGA basis functions in view of quadrature rules.And on the other hand, we would leave in some sense the idea of iso-geometric analysis.

Smooth solution on a square shell
The first example considering the shell formulation is a flat shell of the square domain Ω = [0, 2] 2 under a smoothly distributed load.It is dedicated to investigating the convergence rates of the approach.The parametric domain coincides with the physical domain in this example.The properties are E = 10 6 , ν = 0.1 and correspondingly t = 12(1 − ν 2 )/E 1/3 .The boundary ∂Ω is fixed and a load g is applied as g z = 4π 4 sin(πx) sin(πy), such that the analytical solution of the deformation is u z,ref = sin(πx) sin(πy).The scaling center is placed in the middle of the structure and the mesh is evaluated for third, fourth and fifth order of basis functions.An exemplary mesh for h = 1/4, p = 3, r = 1 for each patch and a corresponding deformation plot is shown in Fig. 11.
The convergence rates indicate for both error estimations optimal convergence rate O(h p−1 ) for the H 2 seminorm and O(h p+1 ) for the L 2 norm.The deviations for p = 5 might be caused by large condition numbers together with rounding errors.

Scordelis-Lo roof
The Scordelis-Lo roof is a well-known model described by [35].It is a cylinder cutout of 80 • with dimensions of radius R = 25 and length L = 50.The thickness is t = 0.25 and the material parameters are E = 4.32 • 10 8 and ν = 0.0.The structure is supported by rigid diaphragms at the curved edges and free at the straight edges.The Scordelis-Lo roof is investigated in three setups as they are the untrimmed structure, trimmed with an elliptic hole under gravity load, and trimmed with four holes under a single load.The untrimmed example is investigated to compare the boundary parametrization to standard IGA parametrizations.Further, the trimmed roofs are compared to numerical results obtained from the literature.The exact parametrization function is chosen to with L Ω,θ 1 and L Ω,θ 2 being the corresponding side length of Ω in each parametric direction.

Untrimmed Scordelis-Lo roof
At first, the untrimmed roof under a gravity load of g z = 90 per unit area is investigated and the reference solution is given as the vertical displacement at the midpoint of the straight, free edge as u ref = 0.3024 [29], [28].Since the reference considers shear deformation which is not represented in Kirchhoff-Love shells, the reference solution herein refers to [6] with u ref = 0.3006 in terms of the vertical deflection at the same point of interest.The parametric mesh, defined on the domain Ω = [−0.5,0.5], [−1, 1], and the corresponding initial shell configuration are shown in Fig. 13.The parametric mesh considers an offset of the scaling center with θ 1 of f = 0.2 and θ 2 of f = 0.1 to show the generality of the formulation.The results of basis functions of order p = 3, p = 4, and p = 5 are shown in Fig. 14 on the left side.Even though the scaling center is chosen in non-centered discretization, the results converge in a smooth way towards the reference solution.A comparison to a standard IGA computation of the Scordelis-Lo roof shows that the convergence of IGA-basis functions is at better results for fewer degrees of freedom.This is a natural consequence of the SB-IGA approach since four patches are coupled towards one scaling center in this example.This points out, that the SB-IGA formulation is not as computationally efficient as IGA for standard single-patch structures.However, the example shows a similar convergence rate for the IGA and the SB-IGA approach which validates the proposed method.
obtained by overkill solution of the method for a Kirchhoff-Love shell proposed in [5]   The results of the vertical deformation at the investigated point of the model with respect to the reference solution in Fig. 16 show good agreement.Accurate solutions are obtained even from the initial mesh.A slight underestimation is visible for the converged however, the difference is less than 0.5% for n = 6 and p = 5.On the right side, the deformation plot for p = 5 and h = 1/6 is figured.The results of the model can be compared to deformation plots presented in [5, Figure 11 (a)], which proposes an IGA shell with adaptive refinement.Fig. 18 shows the deformed structure with the displacement in z-direction for the SB-IGA approach of p = 3, r = 1 and h = 0.125.The results point out, that the approach derived herein is capable of computing the deformation for trimmed domains of shells.The deformation plot has its maximum and minimum at u z,min ≈ −3.3 • 10 −1 and u z,max ≈ 3.7 • 10 −1 , respectively.These results are also obtained in [5].Good agreement is visible on the whole domain including the area of the load application and the free edges.

Hyperbolic paraboloid
The model of the hyperbolic paraboloid is a saddle structure described in [31].Besides its complexity due to the double-curved structure and the negative curvature, the problem is compared to the ASG 1 approach for standard IGA multi-patches.The structure is clamped one-sided and loaded by a uniformly distributed load of g z = −8000 • t.The properties are length L = 1, Young's modulus E = 2 • 10 11 , Poisson's ratio ν = 0.3 and thickness of t = 1/1000 and t = 1/100.The exact parametrization function of both hyperbolic paraboloid shells is chosen as The reference solution is obtained from [31] as the vertical deflection at the middle of the free, curved edge with u z,t=1/1000 = −0.0063941and u z,t=1/100 = −0.93137• 10 −4 .The partition of the parametric mesh of domain Ω = [−0.5,0.5] 2 and the initial shell configuration is shown in Fig. 19 with the clamped boundary highlighted in blue and the point of interest highlighted in red.The scaling center is placed with an offset in θ 2 -direction as θ 2 of f = −0.1.This yields the possibility to have a higher mesh density towards the Dirichlet boundary conditions naturally enforced by moving the scaling center.Further, the proposed approach is compared to the analysis-suitable G 1 approach for standard IGA patches.Fig. 21 shows a comparison of the results for the order of p = 3 and p = 4.It is important to note, that the approach herein has a regularity of r = 1 for various orders, while the compared approach has a regularity of r = p − 2. The comparison shows, that the Kirchhoff-Love shell in boundary representation seems to converge faster towards the reference solution than the approximate solutions from [18] utilizing non-degenerate IGA multi-patches.

Flat shell with multiple holes
To take the trimming more into account, a flat shell structure with an asymmetrical arrangement is evaluated following [32].The model is presented in Fig. 22 with the parametric mesh and the corresponding deformation.The quadratic shell of side length L = 5 and thickness t = 0.005 is trimmed by 16 circular holes of radius r = 0.025 and clamped at all four sides.The material parameters are E = 8.736 • 10 7 and ν = 0.3.The structure is loaded by a vertical load of g z = sin (πx) • sin (πy).To compare the results to reference solutions from the literature [32], a normalization of the vertical deformation is done as with u z as the computed deflection in z-direction, D the flexural stiffness and L the side length of the shell.Fig. 23 shows the results of the proposed approach for the deformation u z and to compare the results to the reference solution from the literature, the results are evaluated in terms of a normalization.We refer to the reference solution of an SB-FEM plate formulation with scaling in the thickness direction and semi-analytic solution procedure presented in [32,Figure 29 (a)].Good agreement is visible for the proposed approach considering maximal normalized deflection and the contour lines.

Violin
In the last example, we utilize the boundary representation in the most sophisticated extent by evaluating a sketch of a violin following the idea of [3] with modifications for the boundary conditions and the parametrization functions.The violin in this contribution is clamped at the whole outer boundary and trimmed by two so-called "F-holes".It is loaded by a constant vertical dead load of g z = 0.05 per unit area.The parameters are t = 0.25, E = 10 5 and ν = 0.1.The underlying parametrization function of the violin is R(θ 1 , θ 2 ) = θ 1 θ 2 2 exp −0.0025 θ 1 2 • 2 exp −0.01 The parametric mesh and the corresponding initial shell configuration are shown in Fig. 24.The advantages of the scaled boundary approach are especially visible due to the various shapes of the SB-patches including triangles, quadrilaterals, and pentagons.Further, the outer boundary curve of the violin contains vertices that are easily included in the computational domain.Even though there is no reference solution, the results seem reasonable as the deformation is symmetric in the y-axis, and the boundary conditions are fulfilled.

Conclusions
In this contribution, a Kirchhoff-Love shell was derived in the framework of scaled boundary isogeometric analysis.To ensure C 1 -continuity across patches, a coupling of the scaled boundary geometries was applied that fits the concept of analysis-suitable G 1 parametrizations from [17] and a special treatment of the basis functions in the vicinity of the scaling center was implemented.The benefits of the boundary representation technique are pointed out by the incorporation of trimming curves as part of the geometry.The shell formulation is at first tested in general against optimal convergence rates of the H 2 seminorm and the L 2 norm.Further, a comparison to the Kirchhoff-Love shell proposed by [6] and in the terms of multi-patch structures [18] is done.The application of trimming is presented in several examples and good agreement with examples from the literature is obtained.The method is especially powerful when it comes to multi-patch geometries that cannot easily be described by a single IGA patch which is outlined by an example of a violin structure with trimmed F-holes.
In further work, an exact representation of the shells in the physical domain is pending.A procedure to adjust the planar domain and the corresponding parametrization according to the physical domain will improve the applicability of the approach.Besides, the formulation of the Kirchhoff-Love shell may be extended to geometric and material nonlinearity to show the formulations' power not only for the linear case but to make it applicable to more advanced shell problems including buckling or plastic deformations.

Figure 1 : 1 ΩFigure 2 :
Figure 1: Within SB-IGA a boundary description is used, where a well-chosen scaling center is required.

Figure 3 :
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 desired global smoothness.

Figure 4 :
Figure 4: Our shell middle surface is parametrized by some mapping R. We assume in next sections that the domain Ω of this R is given as a SB multi-patch domain.

Figure 6 :
Figure 6: An illustration of the scaling center basis functions for a planar 3-patch example.Underlying (patch-wise) degree and regularity are p = 3, r = 1.
(a) Original domain with trimming curve γT (b) Trimmed domain with cut lines.(c)Possible SB-mesh after the trimming.

Figure 9 :
Figure 9: 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.

Figure 10 :
Figure 10: Here we see the difference between the exact parametrization of the Scordelis-Lo roof with holes and the approximated parametrization mapping utilizing the coupled basis functions.To be more precise, the Euclidean norm difference |R approx (θ 1 , θ 2 ) − R(θ 1 , θ 2 )| is plotted, where the underlying meshes are displayed in Fig.25and Fig.19.Even for relatively coarse meshes we see only small differences and the geometry errors can be assumed to be small.
with order p = 5 and 62.456 active elements.The parametric mesh of domain Ω = [−4, 4] 2 and the initial shell configuration is shown in Fig. 15.The point of interest is marked in red in the figure.Ω (a) Parametric mesh R L = 5 0 θ (b) Initial shell configuration

Figure 15 :
Figure 15: Parametric mesh of the Scordelis-Lo roof with hole and the initial shell configuration of h = 1/2.
Normalized displacement uz/u ref u z (b) Deformation plot of component uz

Figure 16 :
Figure 16: On the left side, the convergence study on of the Scordelis-Lo roof with one hole under gravity load is figured.On the right side a plot of the deformation component u z is shown for a mesh of h = 1/6 and p = 5, r = 1.

Figure 21 :
Figure21: Comparison of the results from the proposed approach with the standard analysissuitable G 1 for Kirchhoff-Love shells proposed in[18] for a thickness of t = 1/100 (left) and an exemplary deformation plot of the overall deformation u for p = 3 and h = 1/24 (right).

Figure 22 :
Figure 22: Parametric representation of the flat shell with multiple holes and an exemplary deformation plot with p = 3 and h = 1/3.

Figure 23 :
Figure 23: Left, the deformation u z is shown for the flat shall with multiple holes.Right, the normalized vertical displacements u *z of the proposed method look similar to results by an SBFEM formulation using normal scaling presented in[32].

Figure 24 :
Figure 24: Parametric representation and initial shell configuration of the violin.On the left, the parametric mesh for h = 1/2 is shown.On the right, the corresponding initial shell configuration is shown.The violin is clamped at the outer boundary in all directions including rotations.

Figure 25 :
Figure 25: Deformation plot of the proposed approach for a violin with p = 3 and h = 1/2.
d are p l -open, and a given degree vector p := (p 1 , . . ., p d ) for the multivariate case