H M Verhelst
email: h.m.verhelst@tudelft.nl

P Weinmüller
email: pascal.weinmueller@mtu.de

A Mantzaflaris
email: angelos.mantzaflaris@inria.fr

T Takacs

D Toshniwal
email: d.toshniwal@tudelft.nl

A comparison of smooth basis constructions for isogeometric analysis

Keywords: Isogeometric Analysis, Unstructured Splines, Kirchhoff-Love shell, Biharmonic Equation

In order to perform isogeometric analysis with increased smoothness on complex domains, trimming, variational coupling or unstructured spline methods can be used. The latter two classes of methods require a multi-patch segmentation of the domain, and provide continuous bases along patch interfaces. In the context of shell modeling, variational methods are widely used, whereas the application of unstructured spline methods on shell problems is rather scarce. In this paper, we therefore provide a qualitative and a quantitative comparison of a selection of unstructured spline constructions, in particular the D-Patch, Almost-C 1 , Analysis-Suitable G 1 and the Approximate C 1 constructions. Using this comparison, we aim to provide insight into the selection of methods for practical problems, as well as directions for future research. In the qualitative comparison, the properties of each method are evaluated and compared. In the quantitative comparison, a selection of numerical examples is used to highlight different advantages and disadvantages of each method. In the latter, comparison with weak coupling methods such as Nitsche's method or penalty methods is made as well. In brief, it is concluded that the Approximate C 1 and Analysis-Suitable G 1 converge optimally in the analysis of a bi-harmonic problem, without the need of special refinement procedures. Furthermore, these methods provide accurate stress fields. On the other hand, the Almost-C 1 and D-Patch provide relatively easy construction on complex geometries. The Almost-C 1 method does not have limitations on the valence of boundary vertices, unlike the D-Patch, but is only applicable to biquadratic local bases. Following from these conclusions, future research directions are proposed, for example towards making the Approximate C 1 and Analysis-Suitable G 1 applicable to more complex geometries.

Introduction

Present day engineering disciplines depend on Computer-aided design (CAD) and numerical simulation models for physics for design and analysis. Typically, geometries designed in CAD are converted to meshes for an analysis with numerical techniques like Finite Element Methods (FEMs). Since the geometry description in CAD is based on splines whereas meshes for simulation are based on linear geometry approximations, geometric data is lost during this conversion.

Isogeometric analysis [START_REF] Hughes | Isogeometric analysis: CAD, finite elements, NURBS, exact geometry and mesh refinement[END_REF] is the bridge between CAD and Computer-Aided Engineering (CAE), since it is employing splines as a basis for geometric design and numerical analysis. In practice, an isogeometric analysis and optimization workflow can be seen as depicted in Fig. 1. Starting with a geometry from CAD as well as from material parameters, boundary conditions et cetera from CAE, isogeometric simulations and eventually geometry or topology optimization can be performed. The step connecting the inputs from CAD and CAE is referred to as IGA Setup in Fig. 1. This step takes care of the preparation for the simulation step, including the pre-processing of the geometry, if needed, and the construction of the isogeometric discretization space.

Due to the arbitrary smoothness of spline basis functions, isogeometric analysis has several advantages over conventional finite element methods. For example: (i) the introduction of krefinements, which are proven to provide high accuracy per degree of freedom [START_REF] Bressan | Approximation in FEM, DG and IGA: a theoretical comparison[END_REF][START_REF] Sande | Explicit error estimates for spline approximation of arbitrary smoothness in isogeometric analysis[END_REF]; (ii) high accuracy in eigenvalue problems, e.g. for structural vibrations [START_REF] Cottrell | Isogeometric analysis of structural vibrations[END_REF][START_REF] Hughes | Duality and unified analysis of discrete approximations in structural dynamics and wave propagation: Comparison of p-method finite elements with k-method NURBS[END_REF][START_REF] Hughes | Finite element and NURBS approximations of eigenvalue, boundaryvalue, and initial-value problems[END_REF]; or (iii) geometric exactness in parametric design and interface problems, e.g. applied to the parametric design of prosthetic heart valves [START_REF] Xu | A framework for designing patient-specific bioprosthetic heart valves using immersogeometric fluid-structure interaction analysis[END_REF]. Furthermore, the C 1 -smooth discretization spaces allow to solve equations such as the biharmonic equation, the Cahn-Hilliard equations or the Kirchhoff-Love shell equations without introducing auxiliary variables. However, due to the tensor-product structure of the spline basis, higher-order smoothness can be enforced easily only on domains that allow simple patch partitions (e.g. an L-shape or an annulus), whereas on geometrically and topologically more complicated domains alternative approaches are required to solve equations that require basis functions of higher-order continuity.

For more complicated domains, the IGA setup block in Fig. 1 involves a pre-processing step of either the geometry, the system of equations or the solution space to solve the original system of equations. In Fig. 2, this pre-processing step is subdivided into three options: trimmed domain approaches, unstructured splines and variational coupling methods. Given an initial geometry (cf. Fig. 3a), the trimmed domain approaches alter the tensor-product domain by defining parts of the domain that are physical or non-physical (cf. Fig. 3b). In case of unstructured splines or variational coupling methods, the geometry is decomposed into multiple different patches (cf. Fig. 3c) on which continuity conditions are enforced by constructing a smooth basis (unstructured splines) or by adding extra terms to the system of equations (variational coupling approaches). In Section 2 of this paper, a review of trimmed domain approaches, unstructured splines and variational coupling methods is provided. Examples include immersed methods, degenerate patches and Nitsche's method, respectively. In case of simple geometries (and given the right inputs) the methods are identical.

As shown in Fig. 2, each class of methods has its own characteristics and previous work has provided several comparisons of methods among each other, which are elaborated more in Section 2. In the context of the workflow sketched in Fig. 1, unstructured splines provide a valuable alternative to the other methods, since they are constructed for a fixed topology and hence the computational costs of their construction are not related to changing shapes or moving domains. However, recent developments mainly focused on different unstructured spline methods separately, rather than providing a valuable comparison. In this paper, we therefore provide a qualitative and a quantitative comparison of a selection of unstructured spline constructions.

We consider finite, piece-wise polynomial spline constructions, hence we do not include rational constructions or infinite representations, such as subdivision surfaces. More precisely, we compare examples of (globally) G1 -smooth multi-patch constructions (the Analysis-Suitable G 1 construction of [START_REF] Farahat | Isogeometric analysis with C1-smooth functions over multi-patch surfaces[END_REF] and the Approximate C 1 construction of [START_REF] Weinmüller | Construction of approximate C1 bases for isogeometric analysis on two-patch domains[END_REF]), the D-Patch method of [START_REF] Toshniwal | Smooth cubic spline spaces on unstructured quadrilateral meshes with particular emphasis on extraordinary points: Geometric design and isogeometric analysis considerations[END_REF] and the Almost-C 1 construction of [START_REF] Takacs | Almost-C1 splines: Biquadratic splines on unstructured quadrilateral meshes and their application to fourth order problems[END_REF], motivated in Section 2.3. The selected methods are qualitatively compared based on their properties, and quantitatively based on several different examples with biharmonic and Kirchhoff-Love shell equations. The aim of this paper is to provide a fair comparison 1 of these methods, providing a good overview of the strengths and weaknesses of each method in different cases.

The paper is outlined as follows: in Section 2, a detailed overview of the methods appearing in Fig. 2 is provided. In Section 3 we provide a qualitative analysis of the four constructions that are discussed in this paper, while in Section 4 we provide a quantitative analysis of all methods.

There, we present five benchmark problems solving either a biharmonic or a Kirchhoff-Love equation. These benchmark problems serve different purposes and we compare which method, in which setting, performs best. In Section 5, we conclude this paper based on the findings from the previous sections and we provide directions for future research.

Variational coupling

Figure 2: Inside the IGA Setup block from Fig. 1, three methods are distinguished. Firstly, trimmed domain approaches use trimming curves or surfaces to identify parts of a tensor-product domain as the actual domain. However, since elements can be trimmed poorly, specialized quadrature rules and solver preconditioners are typically needed. Alternatives to trimming are weak coupling or unstructured spline methods. For both classes of methods, a geometry with a given topology needs to be decomposed into multiple sub-domains (i.e. patches) via quadrilateral meshing. Given a quadrilateral mesh, weak methods assemble extra penalty terms into the equation to be solved, or add extra equations to be solved to satisfy continuity constraints. Lastly, unstructured spline constructions can be used to couple multiple domains by constructing a continuous basis. These methods, however can only be used on manifold geometries and conforming meshes. When these requirements are satisfied, unstructured spline pre-processing is required before the unstructured spline construction can take place. The pre-processing is highlighted and will be elaborated on more in Fig. 7 in Section [START_REF] Sande | Explicit error estimates for spline approximation of arbitrary smoothness in isogeometric analysis[END_REF]. Here, the domain is decomposed into several patches Ω i which together define the full domain Ω.

Multi-patch isogeometric analysis: literature review

As discussed in the introduction of this paper, in particular in Fig. 2, three classes of methods for the modelling of complicated domains can be characterised: trimmed domain approaches, variational coupling methods and unstructured splines. The goal of all methods is to achieve a certain level of continuity across the whole analysis domain such that multi-patch isogeometric analysis can be performed for example for the Kirchhoff-Love shell model [START_REF] Kiendl | Isogeometric shell analysis with Kirchhoff-Love elements[END_REF], the biharmonic equation or the Cahn-Hilliard equation [START_REF] Gómez | Isogeometric analysis of the Cahn-Hilliard phase-field model[END_REF]. As shown in Fig. 3, trimmed domain approaches use the fact that parts of tensor-product geometries are trimmed away, using trimming curves to separate regions of interest and regions that should be omitted, see Fig. 3b. Variational coupling approaches and unstructured splines are defined on multi-patch domains, typically following from a segmentation of the original domain, see Fig. 3c. In case of variational coupling methods, the system of equations is enriched with terms that will enforce continuity (typically in a weak sense) between the patches. In case of unstructured spline constructions, a basis is constructed on the multi-patch object, where certain smoothness is enforced strongly. When starting from a trimmed geometry, the step of creating a multi-patch domain decomposition (i.e. untrimming) from an arbitrary geometry with an arbitrary topology is a very important step in the application of weak coupling methods and unstructured spline constructions, as can be seen in the flowchart in Fig. 2. In this paper, however, the topic of untrimming will not be discussed as it is out of scope of our study. Hence, the reader is referred to [START_REF] Massarwi | Untrimming: Precise conversion of trimmed-surfaces to tensor-product surfaces[END_REF][START_REF] Hiemstra | Towards untrimmed NURBS: CAD embedded reparameterization of trimmed B-rep geometry using frame-field guided global parameterization[END_REF] for an overview of these methods.

In this section, an overview of the trimmed domain approaches (Section 2.1), variational coupling methods (Section 2.2) and unstructured splines (Section 2.3) is provided. A fourth method, which will not be discussed in this section, is to introduce auxiliary variables for derivatives of the solution, so that C 1 continuity requirements are reduced to C 0 and standard interface coupling can be used. These so-called mixed formulations are common in conventional FEM, although recent advances have also been made for Kirchhoff-Love plates and shells and the biharmonic eigenvalue problem [START_REF] Rafetseder | A decomposition result for kirchhoff plate bending problems and a new discretization approach[END_REF][START_REF] Rafetseder | A new mixed approach to Kirchhoff-Love shells[END_REF][START_REF] Kosin | A new mixed method for the biharmonic eigenvalue problem[END_REF].

Trimming approaches

Trimming is a technique where so-called trimming curves or surfaces separate parts of tensorproduct spline domains to define a geometry. Trimming is a common technique to represent complex geometries in CAD, and typically geometries consist of multiple trimmed patches with boundary and interface curves trimming the actual patches. We refer to the work [START_REF] Marussig | A Review of Trimming in Isogeometric Analysis: Challenges, Data Exchange and Simulation Aspects[END_REF] for an overview of trimming methods in isogeometric analysis. Generalizing the idea of trimming to techniques where curves or surfaces are used to define the domain of interest as trimmed domain approaches, several approaches have been proposed to perform simulations on complex geometries, including the finite cell method [START_REF] Parvizian | Finite cell method[END_REF][START_REF] Düster | The finite cell method for three-dimensional problems of solid mechanics[END_REF][START_REF] Schillinger | The Finite Cell Method: A Review in the Context of Higher-Order Structural Analysis of CAD and Image-Based Geometric Models[END_REF], Cut-FEM [START_REF] Burman | CutFEM: Discretizing geometry and partial differential equations[END_REF] or immersed methods [START_REF] Hsu | Direct immersogeometric fluid flow analysis using B-rep CAD models[END_REF][START_REF] Kamensky | An immersogeometric variational framework for fluid-structure interaction: Application to bioprosthetic heart valves[END_REF].

The advantage of these methods is that the trimmed CAD geometries could directly be used for analysis. However, when only small parts of the physical domain are cut, leading to small cut elements, numerical difficulties can occur in the conditioning of the system, leading to solver instabilities or accuracy problems [START_REF] De Prenter | Stability and conditioning of immersed finite element methods: Analysis and remedies[END_REF]. Therefore, the analysis of complex trimmed geometries via methods like the FCM typically require special quadrature schemes to take into account small cut cells [START_REF] Divi | Error-estimate-based adaptive integration for immersed isogeometric analysis[END_REF] or preconditioners to stabilize the numerical analysis [START_REF] De Prenter | Multigrid solvers for immersed finite element methods and immersed isogeometric analysis[END_REF]. In the context of Kirchhoff-Love shell modelling, isogeometric analysis on trimmed geometries has been performed in several studies [START_REF] Guo | Weak Dirichlet boundary conditions for trimmed thin isogeometric shells[END_REF][START_REF] Guo | Variationally consistent isogeometric analysis of trimmed thin shells at finite deformations, based on the STEP exchange format[END_REF][START_REF] Coradello | Hierarchically refined isogeometric analysis of trimmed shells[END_REF] including some with focus on multi-patch coupling [START_REF] Coradello | Hierarchically refined isogeometric analysis of trimmed shells[END_REF][START_REF] Coradello | Adaptive isogeometric analysis on two-dimensional trimmed domains based on a hierarchical approach[END_REF][START_REF] Coradello | Coupling of non-conforming trimmed isogeometric Kirchhoff-Love shells via a projected super-penalty approach[END_REF][START_REF] Guo | A parameter-free variational coupling approach for trimmed isogeometric thin shells[END_REF][START_REF] Leidinger | Explicit dynamic isogeometric B-Rep analysis of penalty-coupled trimmed NURBS shells[END_REF].

Variational coupling methods

We define variational coupling methods as methods that modify the system of equations to enforce certain continuity across patch interfaces. Examples of these methods are penalty methods, Nitsche's methods, mortar methods or Langrangian penalized methods. In the context of Kirchhoff-Love shell analysis, these weak coupling methods have received a lot of attention in previous studies and an overview is provided by [START_REF] Paul | Isogeometric continuity constraints for multi-patch shells governed by fourth-order deformation and phase field models[END_REF]. Firstly, an in-plane coupling was proposed in [START_REF] Kiendl | Isogeometric analysis and shape optimal design of shell structures[END_REF] together with a method for coupling non-manifold patches using the so-called bending strip method [START_REF] Kiendl | The bending strip method for isogeometric analysis of Kirchhoff-Love shell structures comprised of multiple patches[END_REF]. Later, weak coupling approaches have been developed for multi-patch domains. Here, coupling terms can be added inside the existing variational formulation (referred to as Nitsche's or penalty methods) or imposed by Lagrange multipliers (referred to as mortar methods).

Several works on Nitsche techniques (cf. [START_REF] Nitsche | Über ein Variationsprinzip zur Lösung von Dirichlet-Problemen bei Verwendung von Teilräumen, die keinen Randbedingungen unterworfen sind[END_REF]) for isogeometric analysis have been published

starting from the imposition of boundary conditions [START_REF] Ruess | Weakly enforced essential boundary conditions for NURBS-embedded and trimmed NURBS geometries on the basis of the finite cell method[END_REF], towards multi-patch coupling and the coupling of patches [START_REF] Ruess | Weak coupling for isogeometric analysis of non-matching and trimmed multi-patch geometries[END_REF], later using a non-symmetric parameter-free Nitsche's method [START_REF] Schillinger | The non-symmetric Nitsche method for the parameter-free imposition of weak boundary and coupling conditions in immersed finite elements[END_REF].

Nitsche's methods have been applied to Kirchhoff plates [START_REF] Hu | Skew-symmetric Nitsche's formulation in isogeometric analysis: Dirichlet and symmetry conditions, patch coupling and frictionless contact[END_REF], Kirchhoff-Love shells [START_REF] Guo | Weak Dirichlet boundary conditions for trimmed thin isogeometric shells[END_REF][START_REF] Guo | Nitsche's method for a coupling of isogeometric thin shells and blended shell structures[END_REF][START_REF] Benzaken | Nitsche's method for linear Kirchhoff-Love shells: Formulation, error analysis, and verification[END_REF], hyperelastic 2D elasticity [START_REF] Du | Nitsche's method for non-conforming multipatch coupling in hyperelastic isogeometric analysis[END_REF] and the biharmonic equation [START_REF] Weinmüller | Construction of approximate C1 bases for isogeometric analysis on two-patch domains[END_REF][START_REF] Weinmüller | An approximate C1 multi-patch space for isogeometric analysis with a comparison to Nitsche's method[END_REF] and for modelling local subdomains [START_REF] Bouclier | A Nitsche-based non-intrusive coupling strategy for global/local isogeometric structural analysis[END_REF] for elasticity simulations. The advantages of Nitsche's methods are that the formulation is variationally consistent and requires only mild stabilization, which can be performed automatically, by estimating the stability parameter. However, the involved integral terms are complicated expressions that impose high implementation and assembly efforts. Therefore, coupling approaches using only penalization have been developed [START_REF] Duong | A new rotation-free isogeometric thin shell formulation and a corresponding continuity constraint for patch boundaries[END_REF][START_REF] Breitenberger | Analysis in computer aided design: Nonlinear isogeometric B-Rep analysis of shell structures[END_REF][START_REF] Herrema | A framework for isogeometric-analysis-based optimization of wind turbine blade structures[END_REF][START_REF] Leonetti | A robust penalty coupling of non-matching isogeometric Kirchhoff-Love shell patches in large deformations[END_REF][START_REF] Pasch | A priori penalty factor determination for (trimmed) NURBS-based shells with Dirichlet and coupling constraints in isogeometric analysis[END_REF]. Although several improvements have been made in these works, the main disadvantage of penalty methods is that a suitable penalty parameter has to be chosen. Using the super penalty approach [START_REF] Coradello | Coupling of non-conforming trimmed isogeometric Kirchhoff-Love shells via a projected super-penalty approach[END_REF][START_REF] Coradello | A projected super-penalty method for the C 1 -coupling of multi-patch isogeometric Kirchhoff plates[END_REF], the computation of the penalty parameter can be automated. However, this method has not yet been tested for non-linear shell problems or on 'dirty' geometries. Both Nitsche's and penalty methods can be used to couple geometries that are non-manifold, i.e. geometries that have out-of-plane connections like stiffened structures by penalizing changes in the angle of patches on an interface. Furthermore, the methods can handle interfaces with non-matching parameterizations.

Instead of adding coupling terms in the variational form, as is done in Nitsche's and penalty methods, mortar methods [START_REF] Bernardi | Domain Decomposition by the Mortar Element Method[END_REF] add extra degrees of freedom by introducing Lagrange multipliers which are required to resolve additional coupling conditions. The use of mortar methods to couple non-conforming isogeometric sub-domains was first done by [START_REF] Hesch | Isogeometric analysis and domain decomposition methods[END_REF]. In [START_REF] Dornisch | The weak substitution method -an application of the mortar method for patch coupling in NURBS-based isogeometric analysis[END_REF] the FEA-based approach of [START_REF] Bernardi | Domain Decomposition by the Mortar Element Method[END_REF] was extended for NURBS-based IGA, but the aim was to develop a method for C 0 -coupling for Reissner-Mindlin shells, hence insufficient for isogeometric Kirchhoff-Love shells. A mortar method aiming to establish C 1 coupling is given in [START_REF] Bouclier | Development of a new, more regular, mortar method for the coupling of NURBS subdomains within a NURBS patch: Application to a non-intrusive local enrichment of NURBS patches[END_REF] and a method that provides C n continuity was given by [START_REF] Dittmann | Weak Cn coupling for multipatch isogeometric analysis in solid mechanics[END_REF][START_REF] Dittmann | Crosspoint modification for multi-patch isogeometric analysis[END_REF]. Furthermore, G 1 mortar coupling, referred to as extended mortar coupling, was presented in [START_REF] Schuß | Multi-patch isogeometric analysis for Kirchhoff-Love shell elements[END_REF] for Kirchhoff-Love shells, based on a coupling in least square sense. On the other hand, in [START_REF] Benvenuti | Isogeometric multi-patch $C 1$-mortar coupling for the biharmonic equation[END_REF] a mortar method to enforce C 1 coupling for the biharmonic equation was developed, where the Lagrange multiplier spaces are constructed similarly to [START_REF] Brivadis | Isogeometric mortar methods[END_REF] for C 0 -coupling. An approach to reduce computational costs involved in finding Lagrange multipliers is called dual mortaring [START_REF] Wohlmuth | A Mortar Finite Element Method Using Dual Spaces for the Lagrange Multiplier[END_REF], where Lagrange multipliers are eliminated using a compact dual basis. This approach has been developed for Bezier elements [START_REF] Zou | Isogeometric Bézier dual mortaring: Refineable higher-order spline dual bases and weakly continuous geometry[END_REF] and it has been applied for Kirchhoff-Love shells [START_REF] Miao | Isogeometric Bézier dual mortaring: The Kirchhoff-Love shell problem[END_REF] and a bi-orthogonal spline space has been presented for weak dual mortaring for patch coupling [START_REF] Wunderlich | Biorthogonal splines for optimal weak patchcoupling in isogeometric analysis with applications to finite deformation elasticity[END_REF]. In [START_REF] Horger | A hybrid isogeometric approach on multi-patches with applications to Kirchhoff plates and eigenvalue problems[END_REF], a hybrid method was provided and applied to Kirchhoff plates, which combines mortar methods and penalty methods. Lastly, a comparison of Nitsche, penalty and mortar methods is given by [START_REF] Apostolatos | Weak imposition of constraints for structural membranes in transient geometrically nonlinear isogeometric analysis on multipatch surfaces[END_REF]. For a more complete overview of mortar methods for isogeometric analysis, the reader is referred to [START_REF] Hesch | Frontiers in Mortar Methods for Isogeometric Analysis[END_REF]. In general, mortar methods have the advantage over Nitsche's methods that there are no parameters involved and that the implementation efforts are lower. However, the disadvantage is that a suitable spline space needs to be found for the Lagrange multipliers [START_REF] Benvenuti | Isogeometric multi-patch $C 1$-mortar coupling for the biharmonic equation[END_REF][START_REF] Brivadis | Isogeometric mortar methods[END_REF][START_REF] Dornisch | An isogeometric mortar method for the coupling of multiple NURBS domains with optimal convergence rates[END_REF]. Like Nitsche's and penalty methods, mortar methods can handle non-matching parameterizations and non-manifold interfaces, the latter by similar penalization of interfacing patches.

Unstructured splines

Compared to weak coupling methods, unstructured spline constructions do not alter the system of equations to be solved. Instead, the computational basis is modified such that it satisfies continuity conditions across patch interfaces. Unstructured splines are typically constructed for in-plane (i.e. manifold) interfaces and not on out-of-plane (i.e. non-manifold) interfaces, since the notion of smoothness is uniquely defined only in the former setting. However, unstructured spline constructions for non-manifold interfaces are possible, e.g. as in [START_REF] Cirak | Subdivision shells with exact boundary control and non-manifold geometry[END_REF][START_REF] Ying | Nonmanifold subdivision[END_REF][START_REF] Moulaeifard | Subdivide and Conquer: Adapting Non-Manifold Subdivision Surfaces to Surface-Based Representation and Reconstruction of Complex Geological Structures[END_REF] in the context of subdivision. Furthermore, unstructured spline constructions are typically constructed on conforming interfaces, i.e. interfaces with matching meshes, but, as long as the patch parameterizations are matching, this can be overcome by taking the knot vector union of the interface patches. However, the advantage of unstructured spline constructions is that as soon as the basis is constructed for a certain untrimmed geometry, there are no additional costs involved other than evaluation costs for changing shapes, which make unstructured spline bases suitable for shape optimization problems. In case of topology changes or large changes of the shape, however, the mesh topology of the unstructured spline space has to be changed as well. Unlike weak methods, which are typically based on the introduction of penalties (e.g. in terms of energy), unstructured spline constructions are typically provided as generic geometric methods that are applicable to any equation that requires C 1 coupling across multi-patch interfaces. With the advance of isogeometric analysis, the interest in parametrically C 1 and geometrically G 1 splines has grown. An overview of smooth multi-patch discretizations for isogeometric analysis can be found in [START_REF] Hughes | Smooth multi-patch discretizations in Isogeometric Analysis[END_REF],

and a small overview is provided below. We distinguish between enforcing parametric continuity, i.e., the type of continuity between mesh elements within a regular tensor-product spline patch, and general geometric continuity, cf. [START_REF] Groisser | Matched Gk-constructions always yield Ck-continuous isogeometric elements[END_REF]. In the following, three types of constructions are classified, depending on their continuity on patch interfaces, around vertices and in the patch interior:

• Patch coupling with geometric continuity on patch interfaces and parametric continuity inside patches.

• Patch coupling with parametric continuity everywhere.

• Patch coupling with parametric continuity almost everywhere.

Although other constructions outside of these categories exist, e.g., [START_REF] Majeed | Isogeometric analysis using manifold-based smooth basis functions[END_REF][START_REF] Yuan | Mapped B-spline basis functions for shape design and isogeometric analysis over an arbitrary parameterization[END_REF], our review is restricted to the aforementioned categories since the methods considered in Sections 3 and 4 fall into these categories.

Geometric continuity on patch interfaces and parametric continuity inside patches

This first category of unstructured spline constructions assumes that a fixed C 0 -matching multi-patch parametrization is given. On this multi-patch domain, a C 1 -smooth isogeometric space is constructed. As shown in [START_REF] Groisser | Matched Gk-constructions always yield Ck-continuous isogeometric elements[END_REF], for any isogeometric function the C 1 condition over each interface is equivalent to a G 1 geometric continuity condition of the graph surface corresponding to the function. If the domain is planar and the patches are bilinear, then the C 1 constraints can be resolved and a C 1 spline space was constructed by [START_REF] Kapl | Isogeometric analysis with geometrically continuous functions on twopatch geometries[END_REF] and applied to the isogeometric analysis of the biharmonic equation in [START_REF] Kapl | Isogeometric analysis with geometrically continuous functions on planar multi-patch geometries[END_REF]. It could be shown in [START_REF] Kapl | A family of C 1 quadrilateral finite elements[END_REF] and [START_REF] Grošelj | A super-smooth C1 spline space over planar mixed triangle and quadrilateral meshes[END_REF] that C 1 splines over bilinear quadrilaterals and mixed (bi)linear quadrilateral/triangle meshes possess optimal approximation properties. Furthermore, the work [START_REF] Kapl | Cs-smooth isogeometric spline spaces over planar bilinear multi-patch parameterizations[END_REF] studied the arbitrary C n -smooth spline space for bi-linear multi-patch parameterizations, based on their previously published findings.

Considering general C 0 -matching multi-patch domains, the work of [START_REF] Collin | Analysis-suitable G1 multi-patch parametrizations for C1 isogeometric spaces[END_REF] introduces the class of analysis-suitable G 1 (AS-G 1 ) multi-patch parameterizations which includes bi-linear patches.

This AS-G 1 condition is in general required to obtain optimal approximation properties. The condition implies that the gluing data for G 1 continuity is linear, which is explained in more detail in Section 3.1. While it could be shown in [START_REF] Kapl | Construction of analysis-suitable G1 planar multi-patch parameterizations[END_REF] that all planar multi-patch domains possess AS-G 1 reparametrizations, creating AS-G 1 surface domains is more difficult. Several strategies to achieve this were introduced in [START_REF] Farahat | Isogeometric analysis with C1-smooth functions over multi-patch surfaces[END_REF], thus making C 1 -smooth multi-patch parameterizations applicable to biharmonic equations and isogeometric Kirchhoff-Love shell models [START_REF] Farahat | Isogeometric analysis for multi-patch structured Kirchhoff-Love shells[END_REF]. In the work of [START_REF] Reichle | Smooth multi-patch scaled boundary isogeometric analysis for Kirchhoff-Love shells[END_REF] the construction of [START_REF] Collin | Analysis-suitable G1 multi-patch parametrizations for C1 isogeometric spaces[END_REF] is used to develop a scaled-boundary model for smooth Kirchhoff-Love shells, similar to the approach of [START_REF] Arf | Scaled boundary isogeometric analysis with C 1 coupling for Kirchhoff plate theory[END_REF] for Kirchhoff plates.

Alternatively to constructing an AS-G 1 parameterization, one can relax the smoothness conditions. This was done in [START_REF] Weinmüller | Construction of approximate C1 bases for isogeometric analysis on two-patch domains[END_REF], where the construction of an approximate C 1 (Approx. C 1 ) space is presented. The basis construction is explicit, possesses the same degree-of-freedom structure as an AS-G 1 space, but the C 1 condition is not satisfied exactly but only approximately. It defaults to the AS-G 1 construction when the AS-G 1 requirements are met. In [START_REF] Weinmüller | An approximate C1 multi-patch space for isogeometric analysis with a comparison to Nitsche's method[END_REF] a comparison of the presented space with Nitsche's method was performed, yielding optimal convergence results without the need of a coupling terms. More details on the Approx. C 1 method are provided in Section 3.2.

Parametric continuity everywhere

The starting point for this class of constructions is different from the previous. Here we create smooth splines in a parametric sense between neighboring mesh elements. Such parametric C 1 conditions are easy to resolve, but they lead to singularities at vertices of valencies other than four, so-called extraordinary vertices. This is due to the conflicting coupling conditions on partial derivatives around the EVs, which lead to all partial derivatives to vanish there. Inspired by the Degenerate Patch (D-Patch) approach from [START_REF] Reif | A Refineable Space of Smooth Spline Surfaces of Arbitrary Topological Genus[END_REF], the works of [START_REF] Toshniwal | Smooth cubic spline spaces on unstructured quadrilateral meshes with particular emphasis on extraordinary points: Geometric design and isogeometric analysis considerations[END_REF][START_REF] Nguyen | Refinable C1 spline elements for irregular quad layout[END_REF] provide C 1 smooth spline spaces for multi-patch geometries with parametric smoothness everywhere. On extraordinary vertices (EVs), which is a junction between 3 or 5 or more patches (i.e. valence ν > 2, ν 4), the original D-Patch method shows a singularity of the basis in EVs combined with a reduction of degrees of freedom in this point. An improvement of the D-Patch method was presented in [START_REF] Nguyen | Refinable C1 spline elements for irregular quad layout[END_REF],

by splitting elements around the EVs such that every element is associated to four degrees of freedom. However, this construction does not have non-negativity and is based on PHT splines, which have limited smoothness. A new design and analysis framework for multi-patch geometries was presented in [START_REF] Toshniwal | Smooth cubic spline spaces on unstructured quadrilateral meshes with particular emphasis on extraordinary points: Geometric design and isogeometric analysis considerations[END_REF], based on D-Patches with T-splines for refinement and non-negative splines yielding optimal convergence properties. This was also demonstrated in [START_REF] Casquero | Seamless integration of design and Kirchhoff-Love shell analysis using analysis-suitable unstructured T-splines[END_REF] for isogeometric Kirchhoff-Love shells. In the work [START_REF] Wei | Analysis-suitable unstructured T-splines: Multiple extraordinary points per face[END_REF], it is motivated that this construction can also be used if only one element around the EV is isolated. More details on the D-Patch method are provided in Section 3.3.

Alternatively, subdivision surface based constructions lead to unstructured splines that are parametrically continuous everywhere, cf. [START_REF] Barendrecht | Isogeometric analysis for subdivision surfaces[END_REF][START_REF] Riffnaller-Schiefer | Isogeometric shell analysis with nurbs compatible subdivision surfaces[END_REF][START_REF] Pan | Isogeometric analysis based on extended catmull-clark subdivision[END_REF][START_REF] Barendrecht | Efficient quadrature rules for subdivision surfaces in isogeometric analysis[END_REF][START_REF] Zhang | Subdivision surfaces with isogeometric analysis adapted refinement weights[END_REF]. However, such approaches require an infinite number of polynomial pieces around each EV. Thus, we discard them for our comparison.

Moreover, in general their approximation properties are severely reduced near EVs [START_REF] Takacs | Approximation properties over self-similar meshes of curved finite elements and applications to subdivision based isogeometric analysis[END_REF].

Parametric continuity almost everywhere

As mentioned previously, imposing parametric continuity everywhere leads to singularities at all EVs. Thus, instead of constructing a space with full parametric continuity, spaces with parametric continuity almost everywhere except around the EVs can also be considered. This way, one ends up with regular, smooth rings around EVs which then need to be filled in some way. Such so-called hole-filling techniques are commonplace in geometric modeling and can also be used to construct smooth spaces for isogeometric analysis, cf. [START_REF] Nguyen | C1 finite elements on non-tensor-product 2d and 3d manifolds[END_REF][START_REF] Karčiauskas | Smooth multi-sided blending of biquadratic splines[END_REF][START_REF] Karčiauskas | Generalizing bicubic splines for modeling and IGA with irregular layout[END_REF][START_REF] Karčiauskas | Refinable multi-sided caps for bi-quadratic splines, Vision, modeling, and visualization: 26th international symposium on vision, modeling, and visualization virtual[END_REF][START_REF] Karčiauskas | Multi-sided completion of C2 bi-3 and C1 bi-2 splines: A unifying approach[END_REF][START_REF] Karčiauskas | Least Degree G1-Refinable Multi-Sided Surfaces Suitable For Inclusion Into C1 Bi-2 Splines[END_REF]. We focus here on the simplest possible way of resolving this issue, which is to enforce only C 0 -smoothness near the EVs and G 1 at the EV, namely the Almost-C 1 construction proposed in [START_REF] Takacs | Almost-C1 splines: Biquadratic splines on unstructured quadrilateral meshes and their application to fourth order problems[END_REF]. Similar constructions, which enforce no additional smoothness near EVs were proposed for mixed quadrilateral/triangle meshes in [START_REF]Quadratic splines on quad-tri meshes: Construction and an application to simulations on watertight reconstructions of trimmed surfaces[END_REF] and for arbitrary degree multi-patch B-splines with enhanced smoothness (MPBES) in [START_REF] Buchegger | Adaptively refined multi-patch B-splines with enhanced smoothness[END_REF].

The Almost-C 1 construction we consider here yields piece-wise biquadratic splines which are C 1 in regular regions and which have reduced smoothness around extraordinary vertices, independent of the valence or the location (i.e. interior or boundary EVs). In contrast to that, most commonly used hole-filling approaches yield exactly C 1 -smooth spaces but introduce locally polynomials of higher degree, or require a higher degree to start with, such as the construction presented in [START_REF] Marsala | G1 -Smooth biquintic approximation of Catmull-Clark subdivision surfaces[END_REF], which converts Catmull-Clark subdivision surfaces to G 1 -smooth piece-wise biquintic elements. While exact smoothness is of relevance for geometric modeling, it is not necessary from an analysis point of view.

Qualitative comparison

In the qualitative comparison of this paper, we focus on the properties of different unstructured spline constructions and their implication on the application of these constructions in a workflow as in Fig. 1. More precisely, we comment on the continuity of each construction and their nestedness properties and we aim to provide a set of requirements for the unstructured spline pre-processing block in Fig. 2. Since the qualitative comparison of the considered methods in this paper mostly covers properties of the methods and their implications, mathematical details about the construction or convergence properties are not provided. For more details, the reader is referred to [START_REF] Farahat | Isogeometric analysis with C1-smooth functions over multi-patch surfaces[END_REF] for the Analysis-Suitable G 1 (AS-G 1 ) method, which extends the 2D construction from [START_REF] Kapl | An isogeometric C1 subspace on unstructured multi-patch planar domains[END_REF], to [START_REF] Weinmüller | Construction of approximate C1 bases for isogeometric analysis on two-patch domains[END_REF] for the Approximate C 1 (Approx. C 1 ) method, to [START_REF] Toshniwal | Smooth cubic spline spaces on unstructured quadrilateral meshes with particular emphasis on extraordinary points: Geometric design and isogeometric analysis considerations[END_REF] for the Degenerate Patches (D-Patch) and to [START_REF] Takacs | Almost-C1 splines: Biquadratic splines on unstructured quadrilateral meshes and their application to fourth order problems[END_REF] for the Almost-C 1 method. However, for the qualitative comparison, some key terms are introduced as preliminaries.

Firstly, a quadrilateral mesh (quad mesh) is a mesh of quadrilateral elements, representing a (planar) surface geometry. The quadrilaterals can be represented by tensor B-splines of any degree which can be mapped onto a parametric unit-square. Typically, when the tensor B-spline quadrilaterals have different sizes in different directions or even different refinement levels, assemblies of these patches are typically referred to as multi-patches. An example of a multi-patch is given in Fig. 3c. The conversion of a quad-mesh with many elements to a multi-patch with a smaller number of patches derived from groups of elements can be done using the procedure described in Fig. 4. Here, a half-edge mesh is traversed and elements are collected into groups corresponding to final patches. The vertices of the elements in one group (i.e. patch) form the control net of the bi-linear patch.

Secondly, for parametrically smooth constructions, different classes of vertices are considered. For so-called extraordinary vertices (EVs) , these constructions typically are different. An interior extraordinary vertex (interior EV) is a vertex on a quad mesh on which three or more than four patches meet. The number of patches coming together at a vertex is referred to as the valence, denoted by ν. Furthermore a boundary extraordinary vertex (boundary EV) is a vertex on the boundary of the quad mesh with valence ν ≥ 3. For geometrically smooth constructions, the construction depends on the geometry around the vertex rather than the valence of the vertex.

Hence, for these constructions the notion of EVs is irrelevant.

Lastly, a refinement of a spline space is called nested if the refined spline space is fully contained in the unrefined space. As a consequence, the geometry is exact under element refinement, which is beneficial from an analysis point of view.

Analysis-suitable G 1

The analysis-suitable G 1 (AS-G 1 ) construction is a novel approach in isogeometric analysis that was introduced for planar geometries and surfaces in [START_REF] Collin | Analysis-suitable G1 multi-patch parametrizations for C1 isogeometric spaces[END_REF], but a construction which extends [START_REF] Kapl | An isogeometric C1 subspace on unstructured multi-patch planar domains[END_REF] for planar domains to surfaces is detailed in [START_REF] Farahat | Isogeometric analysis with C1-smooth functions over multi-patch surfaces[END_REF]. This construction ensures that basis functions at interfaces have C approach is based on the concept that G k -smooth surfaces can produce C k -smooth isogeometric functions [START_REF] Groisser | Matched Gk-constructions always yield Ck-continuous isogeometric elements[END_REF]. When dealing with general C 0 -matching multi-patch domains, the so-called AS-G 1 conditions must be satisfied to ensure optimal approximation. If these conditions are met, a C 1 -smooth subspace of the isogeometric space can be constructed, which is sufficiently large.

Such geometries are referred to as analysis-suitable geometries. However, it should be noted that the C 1 -smooth multi-patch isogeometric space generally depends on the geometry, as discussed in [START_REF] Kapl | Dimension and basis construction for analysis-suitable G1 two-patch parameterizations[END_REF]. To overcome this issue, an Argyris-like space was proposed in [START_REF] Kapl | An isogeometric C1 subspace on unstructured multi-patch planar domains[END_REF], which has a dimension that is independent of the geometry.

Given an interface between two patches, the C 1 continuity condition at the interface is defined by a linear combination of tangent vectors and transversal derivatives, which is referred to as gluing data [START_REF] Collin | Analysis-suitable G1 multi-patch parametrizations for C1 isogeometric spaces[END_REF]. The C 1 smooth basis functions at the interface, or more generally at the edge, can be described by the first order Taylor expansion of the trace and the transversal derivative.

It is shown in [START_REF] Collin | Analysis-suitable G1 multi-patch parametrizations for C1 isogeometric spaces[END_REF], that the ideal choice for the space-representation of the trace and transversal derivative is S(p, r -1, h)2 and S(p -1, r -2, h), respectively. These basis functions have local support and are linearly independent, but they depend on the gluing data and, therefore, on the geometry reparameterization itself. To ensure that the basis functions form a C 1 -smooth subspace of the isogeometric space, and to maintain the nestedness of the spline spaces, it is necessary to have gluing data as a linear function which fulfils all analysis-suitable geometries.

For instance, all bi-linear patches meet this requirement. However, if a geometry is not analysissuitable, it can be reparameterized using the technique presented in [START_REF] Kapl | Construction of analysis-suitable G1 planar multi-patch parameterizations[END_REF].

For any vertices in the quad mesh, to describe the C 1 condition is not that straight-forward. In order to keep it general, the vertex basis functions is constructed by the C 2 interpolation using the C 1 basis functions from the corresponding edges. As a consequence, the vertex basis functions also have local support and are linearly independent.

Summarizing, the AS-G 1 construction can be constructed by three different, linearly independent sub-spaces: interior, edge and vertex space. They can be described as follows:

• Interior space: basis functions that have zero values and derivatives on the patch edges and vertices.

• Interface space: basis functions that have vanishing function values up to the second derivatives at the vertices.

• Vertex space: C 2 interpolating functions at the vertex, i.e., basis functions that have nonvanishing C 2 data at the vertex.

The AS-G 1 construction with the interface and vertex constructions as described above are fully C 1 over the whole domain. In addition, the AS-G 1 construction can only be constructed when the degree of the basis is p ≥ 3 and the regularity is reduced as r ≤ p -2.

Figure 5a presents a local region around and EV with valence five with line styles indicating different continuity levels on patch or element boundaries (see the caption of Fig. 5). For the AS-G 1 construction, the continuity at the vertex is C 2 by construction. Furthermore, the continuity at the interior element interfaces is C p-2 due to the restriction on keeping the isoparametric concept. Lastly, since the AS-G 1 construction provides a G 1 surface, the patch interfaces are C 1 by construction [START_REF] Groisser | Matched Gk-constructions always yield Ck-continuous isogeometric elements[END_REF].

In sum, the core ideas behind the AS-G 1 construction are as follows:

• Degree, regularity, continuity

The spline space is fully C 1 , hence suitable to solve fourth-order problems. However, the computation of the space requires analysis-suitability of the parameterization as well as degree p ≥ 3 and regularity r ≤ p -2 for the basis functions.

• Limitations on construction

The space can be constructed on fully unstructured quadrilateral meshes with both interior and boundary extraordinary vertices. The construction of the basis functions is independent of the location or valence of the EVs. However, the analysis-suitability condition imposes a requirement on the geometries on which the construction can be constructed.

Furthermore, the geometry parameterization is not changed.

• Nestedness

The spline spaces are nested.

• Refinement procedure

Refinement procedure is standard (by knot insertion) since the parameterization does not change.

Approximate C 1

The Approximate C 1 construction [START_REF] Weinmüller | Construction of approximate C1 bases for isogeometric analysis on two-patch domains[END_REF] provides, as the name suggests, approximately C 1 continuity on interfaces and vertices, more precisely the construction provides C 1 continuity in the refinement limit. The Approx. C 1 construction shares similarities with the AS-G 1 construction, but the main difference between the construction of the Approx. C 1 and the AS-G 1 spaces is that it relaxes the AS-G 1 condition on the geometry, i.e., it allows geometries with non-linear gluing data. In fact, the exact gluing data are splines with higher polynomial degree and lower regularity or even piece-wise rational. As a consequence, trying to extend the construction for AS-G 1 parameterizations directly to non-AS-G 1 geometries yields complicated basis functions that are challenging to evaluate and integrate accurately. To overcome this issue and obtain a construction with more easily definable basis functions, the gluing data are approximated. However, this approximation means that the C 1 condition is no longer satisfied exactly but only approximately.

By utilizing the approximation of the gluing data, the Approximate C 1 construction incorporates the concept of different spline spaces found in the AS-G 1 construction. In this case, the interior, vertex, and interface basis functions fulfill the same conditions as in the AS-G 1 construction, but the degree and regularity differ between these spaces. Specifically, the sub-spaces for the AS-G 1 construction have p ≥ 3 and r ≤ p -2, while the Approximate C 1 construction employs an interior space with p ≥ 3 and r ≤ p -1, along with vertex and interface spaces that have locally reduced smoothness based on the approximation of the gluing data. Consequently, on the one hand the Approximate C 1 construction restores the potential for maximal smoothness of isogeometric functions in the refinement limit, but the nestedness of the basis is lost. On the other hand, the approximation of the gluing data in the Approximate C 1 construction does not require analysis-suitability for the optimal convergence rate, unlike the AS-G 1 construction. This feature makes the method applicable to more complex geometries. When the Approximate C 1 construction is applied to an analysis-suitable geometry with p ≥ 3 and r ≤ p -2, and the gluing data approximation is exact, the construction becomes equivalent to the AS-G 1 construction. In sum, the core ideas behind the Approx. C 1 construction are as follows:

• Degree, regularity, continuity

The spline space is approximately C 1 and fully C 1 in the limit of refinement. This makes the spline space suitable to solve fourth-order problems. Contrary to the AS-G1 construction, the spline space approximates the gluing data, allowing maximal smoothness in the interior space (r = p -1) for degrees p ≥ 3.

• Limitations on construction

The space can be constructed on fully unstructured quadrilateral meshes with both interior and boundary extraordinary vertices. The construction of the basis functions is independent of the location or valence of the vertices. Contrary to AS-G 1 the analysis-suitability condition is not needed. However, the construction requires a G 1 condition at the interfaces of surfaces.

• Nestedness

The spline spaces are not nested.

• Refinement procedure

Refinement procedure is standard since the parameterization does not change.

D-patch

The relative ease of imposing parametric smoothness for splines has led to the development of degenerate Bezier patches, or D-patches [START_REF] Reif | A Refineable Space of Smooth Spline Surfaces of Arbitrary Topological Genus[END_REF], which can be used to build C 1 smooth splines on unstructured quadrilateral meshes with no boundary extraordinary vertices. The constructions can be formulated for splines of any bi-degree [START_REF] Hughes | Smooth multi-patch discretizations in Isogeometric Analysis[END_REF], and there are no restrictions on their smoothness in the locally-structured regions of the mesh. In the locally-unstructured regions of the mesh (i.e., in a neighbourhood of an extraordinary vertex), the splines are C 1 smooth and first-order degenerate. Note that this degeneracy means that the spline spaces are not necessarily H 2 -conforming, but numerical evidence shows that they can still be used to solve fourth-order problems.

Specifically, imposition of strong C 1 smoothness around an extraordinary vertex requires that the splines vanish up to first order at the extraordinary vertex. This degeneracy trivially implies matching first derivatives at the extraordinary vertex (since all of them vanish) but does not imply C 1 smoothness of the resulting spline functions and the geometries built using them. As

shown in [START_REF] Reif | A Refineable Space of Smooth Spline Surfaces of Arbitrary Topological Genus[END_REF], additional conditions can be imposed upon certain higher-order mixed derivatives to ensure this desired C 1 smoothness. Furthermore, the effect of these additional constraints can be localised to a neighbourhood of the extraordinary vertex by imposing them on a subdivided representation of the splines [START_REF] Nguyen | Refinable C1 spline elements for irregular quad layout[END_REF]. This means that a patch-based representation of C 1 D-patch splines takes functions that are in S(p, r, h/2) on each patch, where almost all basis functions are in S(p, r, h), except a few basis functions supported in a neighbourhood of extraordinary points (the number of basis functions depends on the valence).

The D-patch construction allows for nested refinements of the spline spaces [START_REF] Reif | A Refineable Space of Smooth Spline Surfaces of Arbitrary Topological Genus[END_REF]. If different orders of smoothness are being imposed in locally-structured and locally-unstructured regions of the mesh, then nested refinements produce spline spaces with a higher number of basis functions supported in the vicinity of extraordinary points (the number depends on the refinement-level), see [START_REF] Toshniwal | Smooth cubic spline spaces on unstructured quadrilateral meshes with particular emphasis on extraordinary points: Geometric design and isogeometric analysis considerations[END_REF] for instance. On the other hand, a patch-based approach allows for a simpler implementation by limiting the smoothness across patch interfaces to C 1 ; the smoothness in patch-interiors can still be arbitrarily chosen. However, special care should be taken when using D-patches with nested refinements -the degeneracy of the splines near extraordinary vertices means that, with mesh refinements, the shape regularity of the mesh starts to worsen with refinements and the finite element matrices become very ill-conditioned.

In sum, the core ideas behind the D-patch spline construction are the following:

• Degree, regularity, continuity

The spline space is fully C 1 . In general, the degeneracy of derivatives means that the spaces are H 2 -nonconforming, however numerical evidence supports their use in solving fourth-order problems. The construction can be formulated for splines of any degree and the smoothness away from extraordinary vertices can be chosen arbitrarily.

• Limitations on construction

The space can be constructed on unstructured quadrilateral meshes with no boundary extraordinary vertices.

• Nestedness

The spline spaces can be refined in a nested manner, however the resulting mesh have poor shape regularity and the corresponding finite element matrices may be very ill-conditioned.

• Refinement procedure

Refinement procedures can be derived from standard B-spline knot insertion.

Almost C 1

Almost-C 1 splines are defined on a general, conforming quadrilateral mesh. They are piecewise biquadratic and possess mixed smoothness, i.e., they are C 1 in regular regions, while the smoothness near extraordinary vertices, i.e., vertices with valence different from four, is reduced.

To be precise, they are C 1 smooth at all vertices (including extraordinary vertices) and across all edges except for the ones emanating from an extraordinary vertex. Moreover, while they are defined to be biquadratic on all regular elements, they are piece-wise biquadratic splines (with one inner knot in each direction) on all elements that are neighboring an extraordinary vertex. Details can be found in [START_REF] Takacs | Almost-C1 splines: Biquadratic splines on unstructured quadrilateral meshes and their application to fourth order problems[END_REF]. As a consequence, a patch based representation of Almost-C 1 splines takes functions that are in S(2, 1, h/2) on each patch, where almost all basis functions are in S(2, 1, h), except a few basis functions supported in a 1-ring neighbourhood of extraordinary points (the number depends on the valence).

A central feature of Almost-C 1 splines is the mixed smoothness imposition described above.

In particular, this choice of mixed smoothness only depends on the current refinement level of the mesh. That is, standard C 1 -smoothness is enforced across all edges at the current refinement level except the ones that are incident upon extraordinary vertices, where only C 0 smoothness is enforced. Additionally, these smoothness conditions are combined with G 1 smoothness imposition at each extraordinary vertex. This means that almost-C 1 splines do not yield nested spaces when refining. As a result, the refinement process essentially amounts to a projection of coarse Almost-C 1 splines onto the refined Almost-C 1 spline space. This projection can be chosen in many different ways and can have a significant impact on the limit surface description as well as isogeometric simulations using these spaces. In [START_REF] Takacs | Almost-C1 splines: Biquadratic splines on unstructured quadrilateral meshes and their application to fourth order problems[END_REF] a smoothing and refinement procedure is proposed that results in a C 1 -smooth limit surface for sufficiently regular input data.

Let us briefly summarize the refinement procedure here. We assume that we are given a quad mesh and associate a control point with each face of the mesh. The initial smoothing step guarantees that all control points associated to the one ring around an extraordinary vertex are coplanar.

Having given such an initial control point grid, we then refine the geometry using explicit subdivision rules as specified in [START_REF] Takacs | Almost-C1 splines: Biquadratic splines on unstructured quadrilateral meshes and their application to fourth order problems[END_REF][START_REF]Quadratic splines on quad-tri meshes: Construction and an application to simulations on watertight reconstructions of trimmed surfaces[END_REF]. The rules are the same as for quadratic tensor-product B-splines in regular regions and maintain the coplanarity near extraordinary vertices.

In sum, the core ideas behind the Almost-C 1 spline construction are the following:

• Degree, regularity, continuity

The spline space locally reproduces biquadratic polynomials and it is sufficiently smooth to be able to solve fourth order problems.

• Limitations on construction

The splines can be constructed on fully unstructured quadrilateral meshes, in particular, those that contain both interior and boundary extraordinary vertices.

• Nestedness

Since the spaces are not nested, the convergence behavior of Almost-C 1 splines depends on how the geometry parameterization is refined.

• Refinement procedure

An initial geometry and a refinement procedure can be constructed in such a way, that the limit geometry parameterization is normal continuous everywhere.

Thus, the concept introduced in [START_REF] Takacs | Almost-C1 splines: Biquadratic splines on unstructured quadrilateral meshes and their application to fourth order problems[END_REF] is quite flexible, since the initial smoothing procedure and the refinement procedure are not unique and can be tailored to the needs coming from geometric modeling, e.g., one may want to reproduce Doo-Sabin subdivision surfaces, thus having to modify the subdivision rule for refinement accordingly. The spline space that is introduced on each refinement level can be seen as a simple hole-filling construction, which is sufficient for numerical analysis.

Conclusions

In this section, a summary of the construction and the properties of the analysis-suitable G 1

(AS-G 1 ), the approximate C 1 (Approx. C 1 ), the degenerate patches (D-patch) and the Almost-C 1 methods have been provided, referring to the relevant publications for the mathematical details.

For each method, comments have been provided on the degree, regularity and continuity of the space, on the limitations of the construction in terms of the quadrilateral mesh, on nestedness for refinement and on the refinement procedure itself. In addition, Figs. 5 and6 provides detailed information on the local continuity of the constructions around an extraordinary vertex.

The aim of the qualitative analysis of the methods in this paper is to provide a comparison of a set of properties and requirements of each method and their implications on their applicability. While the subsections presented before provide a brief description of the properties of the methods and the reason behind these properties and requirements, Table 1 provides a side-byside comparison of each method based on the subsections before. In particular, the table lists the (i,ii) requirements on degree and regularity for the constructions, iii geometrical or topological limitations if applicable, (iv,v) the continuity of the constructed bases in the interior and on the interfaces and element boundaries and vi nestedness of the constructed basis.

Following from Table 1, the requirements for construction of the unstructured spline bases are summarized in Fig. 7 as pre-processing conditions that have to be satisfied for each unstructured spline construction in the process depicted in Fig. 2. The degree and regularity conditions (cf. i,ii in Table 1) must be satisfied for each construction, e.g. by performing projections on suitable spline spaces or by knot insertion routines. Furthermore, the geometric or topological limitations (cf. iii in Table 1) impose additional constraints that the geometry must satisfy.

Table 1: Summary of the requirements for the construction and the properties for each of the considered bases. The construction requirements include the degree and regularity of the basis used for construction as well as geometrical or topological properties of the input geometry. The properties include the continuity on interfaces, vertices and in the interior of the unstructured spline construction, as well as the nestedness property.

Requirements

AS-G 1 Approx. C 1 D-Patch Almost-C 1 (i) Degree p ≥ 3 p ≥ 3 p ≥ 3 p = 2 (ii) Regularity r ≤ p -2 r ≤ p -1 r ≤ p -1 r = 1 (v) Geometrical / topological limita- tions Analysis- suitability G 2 continuity BEVs: ν ≤ 3, C 1 continuity C 1 continuity Properties AS-G 1 Approx. C 1 D-Patch Almost-C 1 (iii) Interface & Vertex Continuity C 1 C 1 in the limit C 1 C 1 in the limit (iv) Interior conti- nuity C p-2 C p-1 C p-1 C 1 (vi) Nestedness Yes No Yes No p ≥ 3 r ≤ p -1 p ≥ 3 r ≤ p -2 p ≥ 2 r ≤ p -1 p = 2 r ≤ p -1 From Unstructured Spline Constraints G 2 geometry Analysis suitability C 1 geometry & BEV ν ≤ 3 C 1 geometry
To Unstructured Spline Construction D -P a t c h A p p r o x . The first row represents requirements on degree p and regularity r. If not satisfied, the geometry can be projected onto a space that satisfies the requirement, or degree elevation or reduction steps can be performed together with refinement operations. The second row depicts the requirements on the geometry parameterization; these blocks can be satisfied by changing the geometry.

C 1 A S -G 1 A l m o s t C 1

Quantitative comparison

In this section a quantitative comparison between the methods provided in Section 3 is provided. In addition, variational coupling methods are compared if applicable. The quantitative comparison is composed of various benchmark problems, each providing a different conclusion with respect to the methods considered:

Biharmonic problem on a planar domain (Section 4.1) The first example entails solving the biharmonic problem on a planar domain. The goal of this example is to assess the convergence properties of all considered unstructured spline constructions, hence the problem will be solved on a simple analysis-suitable geometry without EVs on the boundary, such that every method from Section 3 can be applied and compared to the manufactured solution.

Linear Kirchhoff-Love shell analysis on a surface (Section 4.2) The second example entails solving the Kirchhoff-Love shell equation on curved domains. The goal of this example is to demonstrate the performance of the unstructured spline construction for simple shell problems. Therefore, comparison will be made to single-patch results and penalty coupling from [START_REF] Herrema | A framework for isogeometric-analysis-based optimization of wind turbine blade structures[END_REF].

Spectral analysis on a planar domain (Section 4.3) In the third example, spectral analysis of a plate equation is performed. The goal of this example is to assess the spectral properties of the unstructured spline methods compared to a variational approach and a single patch, since the spectral properties of highly continuous bases have been demonstrated to be superior over non-smooth bases [START_REF] Cottrell | Studies of refinement and continuity in isogeometric structural analysis[END_REF].

Modal analysis of a complex geometry (Section 4.4) In the fourth example, a modal analysis is performed on a complex geometry extracted from a quad-mesh. The goal of this example is to demonstrate the applicability and performance of the unstructured spline methods on a large-scale, more complicated geometry.

Stress analysis in a curved shell (Section 4.5) Lastly, the fifth example involves the analysis of stress fields in shells. The goal of this example is to assess the performance of unstructured spline constructions and a penalty method when it comes to stress reconstruction in shells.

For the Kirchhoff-Love shell, the stresses are obtained by taking gradients of the deformed geometry, hence of the solution. This means that for C 1 bases, stresses are C 0 . This might be unfavourable in engineering applications where local stress fields are of importance, e.g. fatigue analysis.

In all examples except the complex geometry in Section 4.4, the domain decomposition from out the benchmark problems. In Table 2 the combinations of p and r and the methods that are compared for these bases are provided. For the biharmonic problem and the spectral analysis (Sections 4.1 and 4.3) Nitsche's method is used for comparison, see [START_REF] Weinmüller | An approximate C1 multi-patch space for isogeometric analysis with a comparison to Nitsche's method[END_REF] for more details. When solving the Kirchhoff-Love shell equations, the penalty method is used for comparison, see [START_REF] Herrema | A framework for isogeometric-analysis-based optimization of wind turbine blade structures[END_REF] Figure 8: Multi-patch decomposition of a simple domain into six patches. The domain has two EVs in the interior (valence 3 and 5) and no boundary EVs.

of the geometry. In case of mesh convergence results, refinements can be performed in different ways. On the one hand, the original geometry can be refined and a new construction with a new geometry approximation can be performed. On the other hand, the geometry resulting from the construction in the first refinement level can be refined in a nested way, such that the geometry does not change after the first mesh. In the quantitative comparison, all refinements are performed in a nested way, unless specified otherwise.

Biharmonic equation on a planar domain

The first benchmark entails the biharmonic equation on a planar domain. The purpose of this example is to assess the convergence properties of the unstructured spline methods described in Section 3. We mainly follow the structure of [START_REF] Weinmüller | An approximate C1 multi-patch space for isogeometric analysis with a comparison to Nitsche's method[END_REF]. The biharmonic equation is solved on a unit square Ω = [0, 1] 2 with the patch segmentation from Fig. 8. The biharmonic equation is defined by

∆ 2 φ = f. (1)
In the present example, convergence is analysed with respect to a manufactured solution φ(x 1 , x 2 ) = (cos(4πx 1 ) -1)(cos(4πx 2 ) -1), [START_REF] Bressan | Approximation in FEM, DG and IGA: a theoretical comparison[END_REF] such that the right-hand-side function becomes:

f (x 1 , x 2 ) = 256π 4 (4 cos(4πx) cos(4πy) -cos(4πx) -cos(4πy)) (3)
Furthermore, on all boundaries of the domain, the manufactured solution and its derivatives are imposed as Dirichlet and Neumann boundary conditions, respectively:

φ = φ(x 1 , x 2 ) ∂ n φ = ∂ n φ      on Γ, (4)
where Γ = ∂Ω, n is the unit outward normal vector on Γ. The biharmonic equation from Eq. ( 1)

with boundary conditions Eq. ( 4) can be discretized by obtaining the weak formulation, see [START_REF] Weinmüller | An approximate C1 multi-patch space for isogeometric analysis with a comparison to Nitsche's method[END_REF],

inserting Eq. ( 4) and by defining an approximation of the solution φ as φ h . Furthermore, a weak coupling can be established through Nitsche's method. For the mathematical details behind the discretization of the biharmonic equation and optionally adding Nitsche interface coupling terms, we refer to [START_REF] Weinmüller | Construction of approximate C1 bases for isogeometric analysis on two-patch domains[END_REF][START_REF] Weinmüller | An approximate C1 multi-patch space for isogeometric analysis with a comparison to Nitsche's method[END_REF]. For the D-Patch and Almost C 1 constructions, the geometry is smoothed upon construction. The geometry used for evaluation of the weak formulation is constructed by using an L 2 -projection of the geometry from the coarsest space which is projected onto the smooth basis of each refinement level. For the D-Patch, the non-negative smoothness matrix for vertex smoothing is used. Although this matrix produces non-nested meshes, it provides the highest rates of convergence. Furthermore, the factor β (cf. [10, sec. 5.1]) is chosen as β = 0.4 as used by [START_REF] Toshniwal | Smooth cubic spline spaces on unstructured quadrilateral meshes with particular emphasis on extraordinary points: Geometric design and isogeometric analysis considerations[END_REF], or β = 1.2, and halved in each refinement level.

To evaluate the unstructured spline constructions from Section 3, the numerical approximation φ h is compared to the manufactured solution φ in the L 2 -, H 1 -and H 2 -norms on the multipatch segmentation from Fig. 8. The bi-linear segmentation is refined and degree elevated until the desired degree p and regularity r from Table 2 are obtained. In addition, a Nitsche coupling of the patches is employed for comparison.

The results for the comparison are presented in Fig. an analysis-suitable geometry with regularity p -2, the approximate gluing data becomes exact, hence the same as in the AS-G 1 construction. The D-patch in this case shows better convergence of the L 2 , H 1 and H 2 errors for β = 1.2 than for β = 0.4. Though, for both choices of β, the convergence is sub-optimal, as was observed in the work by [START_REF] Casquero | Seamless integration of design and Kirchhoff-Love shell analysis using analysis-suitable unstructured T-splines[END_REF]. Furthermore, the L 2 -norm increases at the last point of the D-Patch results, due to ill-conditioning of the system of equations.

Lastly, for degree p = 3 and regularity r = 2 the Approx. C 1 , D-patch and Nitsche's method are compared. The observations are as for the p = 3, r = 1 case.

10 -3 10 -2 10 -1 10 0 10 -9 10 -5 10 -1

10 3 1 1 1 2 1 2 Element size h L 2 , H 1 , H 2 -error p = 2, r = 1 10 -3 10 -2 10 -1 10 0 1 2 1 3 1 4 Element size h p = 3, r = 1 10 -3 10 -2 10 -1 10 0 1 2 1 3 1 4 Element size h p = 3, r = 2 L 2 H 1 H 2 AS-G 1 Approx-C 1 D-Patch, β = 0.4 D-Patch, β = 1.2 Almost-C 1 Nitsche
Figure 9: Errors for the AS-G 1 , Approx. C 1 , D-Patch and Almost-C 1 construction for the biharmonic problem on the domain in Fig. 8. The L 2 , H 1 and H 2 errors with respect to the analytical solution are plotted with different line styles in the top row. Furthermore, all results are plotted against the element size h and the expected convergence rates are given by the triangles.

Overall, the results show expected convergence behaviour for all considered spline constructions compared to theoretical results and compared to a Nitsche coupling method. However, the D-Patch method does not converge for very fine meshes, due to ill-conditioning of the system matrix.

Linear Kirchhoff-Love shell analysis on a surface

We solve the linear Kirchhoff-Love shell equations on two geometries to demonstrate the convergence behaviour of the methods on curved surfaces. To this end, two benchmark examples are considered. Firstly, a hyperbolic paraboloid surface is constructed with shape, inspired by [START_REF] Farahat | Isogeometric analysis with C1-smooth functions over multi-patch surfaces[END_REF]:

r(ξ 1 , ξ 2 ) = ξ 1 ξ 2 ξ 2 1 -ξ 2 2 (5)
The left-side of the hyperbolic paraboloid is clamped (u = 0) and the other sides are free. Furthermore, a distributed load with magnitude 8000t is applied with t the thickness, see Fig. 10.

Secondly, an elliptic paraboloid shaped-domain is modelled, with equation

r(ξ 1 , ξ 2 ) = ξ 1 ξ 2 1 -2 ξ 2 1 + ξ 2 2 (6)
For this shape, a point load with magnitude 10 8 t is applied in the middle of the domain. The corners of the domain are only fixed in vertical z direction to allow sliding in the xy-plane. One corner is fixed in all directions to create a well-posed problem. For both hyperbolic paraboloid (Fig. 10) and elliptic paraboloid (Fig. 11) the multi-patch segmentation from Fig. 8

2 ) = ξ 1 ξ 2 ξ 2 1 -ξ 2 2 , ξ 1 , ξ 2 ∈ [-1/2.
2 ) = ξ 1 ξ 2 1 -2 ξ 2 1 + ξ 2 2 , ξ 1 , ξ 2 ∈ [-1/2.1/2].
On the corners of the domain, the vertical displacements are set to zero u z = 0 and one corner is fixed in-plane as well. Furthermore, a point load with magnitude P = 10 8 t is applied in the middle of the geometry. The results of both analyses are given in Figs. 12 and13. Here, different unstructured spline constructions are tested on patch-bases with different degrees and regularities, as reported in Table 2. For each combination of degree p and regularity r, the energy norm

Kirchhoff material with

W h int = 1 2 u ⊤ h K h u h is
plotted against the number of degrees of freedom, with u h the discrete displacement vector and K h the discrete linear stiffness matrix. From the results in Figs. 12 and 13, a few observations can be made. Firstly, the Approx. C 1 and AS-G 1 methods show slow convergence on the hyperbolic paraboloid geometry, while the convergence on the elliptic paraboloid geometry is similar to the single-patch convergence. The slow convergence for the hyperbolic paraboloid shell is also observed in [START_REF] Farahat | Isogeometric analysis with C1-smooth functions over multi-patch surfaces[END_REF]. Since the results of the same constructions on the elliptic paraboloid geometries do not show slower convergence, the slow convergence is hypothetically a result of the double curvature with different signs of the shell. Secondly, the D-Patch and Approx. C 1 show comparable convergence to the penalty method on both geometries, which is slightly slower than the convergence of the single-patch results. This is explained by the fact that the degrees of freedom are more optimally allocated for the single-patch parameterization. Lastly, the results obtained by the penalty method for different penalty parameters α show convergence with a rate similar to the D-Patch and Almost-C 1 methods for penalty parameters α ∈ {1, 10}. For α = 100 the penalty method is still converging to the same solution, but convergence starts after a few refinement steps.

Spectral analysis on a planar domain

In this example, the spectral properties of the unstructured spline constructions on a multipatch domain are considered. From [START_REF] Cottrell | Isogeometric analysis of structural vibrations[END_REF] it is known that isogeometric analysis has the advan- Figure 12: Bending energy norm W h int = 1 2 u ⊤ h K h u h for the hyperboloid geometry from Fig. 10 with a patch segmentation as in Fig. 8. The results are presented for different combinations of the degree p and regularity r for all unstructured spline constructions. In addition, the results for a penalty method with parameter α ∈ {1, 10, 100} are provided for comparison. Figure 13: Bending energy norm W h int = 1 2 u ⊤ h K h u h for the paraboloid geometry from Fig. 11 with a patch segmentation as in Fig. 8. The results are presented for different combinations of the degree p and regularity r for all unstructured spline constructions. In addition, the results for a penalty method with parameter α ∈ {1, 10, 100} are provided for comparison.

tage over C 0 Finite Element Analysis with respect to spectra for eigenvalue problems. Smooth isogeometric discretization provide converging spectra with spline degree p, whereas the spectra obtained by C 0 FEA diverge with p and typically have optical branches. Similarly, when patches with C 0 continuity are considered, optical branches are introduced and the accuracy of the spectral approximation decreases [START_REF] Nguyen | A variational approach based on perturbed eigenvalue analysis for improving spectral properties of isogeometric multipatch discretizations[END_REF]. In this benchmark problem, we compare the basis constructions from Table 2 on Fig. 8 on their spectral properties. For Nitsche's method, we use different values for the coupling parameter to assess its influence on the spectrum.

For the problem at hand, we consider a unit-square domain with parametric lay-out from Fig. 8 for simplicity. We consider modal analysis using the plate equation. The stiffness operator of the free vibration plate equation is similar to the biharmonic equation from Eq. ( 1), and the inertia is included on the right-hand-side:

D∆ 2 w = -ρt ∂ 2 w ∂τ 2 (7)
Assuming that w(x, y, τ) is harmonic, i.e. w(x, y, τ) = ŵ(x, y) exp{iωτ} with ω a frequency, the equation simplifies to

D∆ 2 ŵ = ω 2 ∂ 2 ŵ ∂τ 2 . (8)
Here, D = Et 8) is a generalized eigenvalue problem with eigenpairs (ω i , v i ) where ω i is the i th eigenfrequency and v i the i th mode shape. The mode shape for a simply supported unit plate with n × m half-waves is given by v nm (x, y) = sin (nπx) sin (mπy) [START_REF] Weinmüller | Construction of approximate C1 bases for isogeometric analysis on two-patch domains[END_REF] with corresponding eigenfrequency

ω nm = (n 2 + m 2 )π 2 D ρt . (10)
In addition, the numerical solution to Eq. ( 8) is obtained by solving the following generalized eigenvalue problem

D Ω ∆w∆φ dΩ = ω 2 ρt Ω uφ dΩ (11)
With φ a test function, see Section 4.1. In further representation of the solutions, we employ the index i such that ω i < ω i+1 and we use the subscript h for numerically obtained solutions.

Figure 14 presents the spectra for different degrees, regularities and for different methods.

Here, the vertical axis represents the ratio of the numerically obtained eigenfrequency over the analytical eigenfrequency with index i, thus ω h,i /ω i . The horizontal axis represents the fraction of the eigenfrequency index i over the total number of eigenmodes. The total number of eigenmodes is equal to the number of degrees of freedom in the system. The results are presented for the degrees and regularities as in Table 2.

Firstly, the p = 2, r = 1 plot shows that Nitsche's method oscillates for all considered values of the penalty parameter. Furthermore, in the part where it is not oscillating, the ratio ω i.h /ω i is higher than for the D-patch and Almost-C 1 method. Additionally, the D-patch and Almost-C 1 methods show a significant difference with respect to the single patch result, which is due to the non-Cartesian multi-patch segmentation of Fig. 8 and the fact that the analytical mode shapes are Cartesian. For the p = 3, r = 1 and p = 3, r = 2 bases similar conclusions can be drawn. Although for the p = 3, r = 1 case the Approx. C 1 method seems worse than the D-patch method, the opposite is true for p = 3, r = 2. Hence, it can be concluded that no method outperforms another, but that all unstructured spline constructions perform better than Nitsche's method.

Modal analysis of a complex geometry

The next example for the quantitative analysis in this paper involves the modal analysis on a larger-scale complex geometry, depicted in Fig. 15a. The goal of this example is to show the usability of the considered constructions on an off-the-shelf industrial geometry. The geometry is represented as a mesh consisting of 15895 vertices, 31086 edges and 62172 faces. This geometry is converted to bi-linear patches using the procedure discussed in Fig. 4 in Section 3.

The interface and boundary curves of the patches are given in Fig. 15b and the final multi-patch object is given in Fig. 15c. The latter has 3 EVs of valence 3, 10 EVs of valence 5 and 16 bEVs.

Moreover, the material parameters specified for a steel material. That is, the density of the material is ρ = 7850 • 10 -6 [tonnes/mm3 ], the shell thickness is t = 10 [mm], the Young's modulus is

E = 210 • 10 3 [MPa]
and the Poisson's ratio is ν = 0.3 [-]. All the sides of the geometry are kept free, meaning that the modal analysis results will consist of six modes with zero eigenfrequencies: the rigid body modes. In the sequel, we list the results for deformation modes only.

After the creation of the linear multi-patch object, h-, p-and k-refinement steps can be performed to construct a multi-basis corresponding to the patch lay-out on which unstructured splines can be constructed. For the Almost-C 1 and D-Patch constructions, the bases are constructed by refining and elevating the initial linear basis up to the desired degree and regularity, after which the the Almost-C 1 and D-Patch basis and geometry are computed. An Almost-C 1 geometry is provided in Fig. 15c.

The AS-G 1 construction requires an analysis-suitable geometry, which can be constructed following [START_REF] Farahat | Isogeometric analysis with C1-smooth functions over multi-patch surfaces[END_REF], which is based on the planar construction developed in [START_REF] Kapl | Construction of analysis-suitable G1 planar multi-patch parameterizations[END_REF]. However, the geometry from Fig. 15b is only C 0 -smooth due to the original linear mesh it is constructed from.

An algorithm to automatically pre-process the geometry to obtain an analysis-suitable G 1 surface is not yet developed. The algorithm from [START_REF] Farahat | Isogeometric analysis with C1-smooth functions over multi-patch surfaces[END_REF] requires AS-G 1 gluing data, which cannot be prescribed directly on a C 0 surface. If the surface is not pre-processed to be AS-G 1 , no suitable gluing data can be found and the basis construction is not applicable. Although the Approx. C 1 construction does not require an analysis-suitable re-parameterization, it does require G 1 smoothness at the interfaces. If this condition is not satisfied there exists no C 1 construction that can be approximated by this method. For both methods, the required pre-processing efforts are nontrivial or not demonstrated on industrial geometries, and therefore left out of the scope of this paper 3 .

Furthermore, penalty methods have been used in the context of modal analysis on a 27 patch composite wind-turbine blade in [START_REF] Herrema | A framework for isogeometric-analysis-based optimization of wind turbine blade structures[END_REF], where the variation of the element size of interface elements seems rather small. In the present paper, an attempt was made to apply the penalty method on the geometry in Fig. 15c, but unidentifiable vibration modes were obtained, possibly because of the large variation of element lengths across the interfaces of the domain, challenging the determination of a suitable penalty parameter α.

Table 3 presents the eigenfrequencies for the first four deformation modes of the car side panel for the D-Patch and the Almost-C 1 constructions with degree p = 2 and regularity r = 1 Table 3: Eigenfrequencies of the Almost-C 1 and D-Patch constructions for the car geometry in Fig. 15. The results of an ABAQUS FEA simulation using the S4R element are provided as a reference. The mode-shapes are plotted in Fig. 16. 3. From these results, it can be observed that the Almost-C 1 and D-Patch methods provide eigenfrequencies in the same range and that the eigenfrequencies are mostly converging in the second digit.

Method

Moreover, the eigenfrequencies of the D-Patch and Almost-C 1 methods for coarse meshes and p = 2, r = 1 already provide reasonable estimates compared to higher degrees and refinements.

On the other hand, the results obtained using an ABAQUS S4R element show convergence in the second digit, and slightly lower frequencies than the IGA results, possibly because the FEM uses a different geometry approximation. Overall, it can be concluded from this benchmark problem that the Almost-C 1 and D-Patch are more robust for industrial and large scale geometries, that are represented by at least C 0 -conforming quadrilateral meshes, compared to the Approx. C 1 and AS-G 1 methods due to the pre-processing efforts required by the latter. Furthermore, these methods are parameter-free, making them robust also with respect to penalty methods.

Stress analysis in a curved shell

An interesting application for smooth unstructured spline construction is for the use of thin shell analysis for engineering applications. Not only displacements (Section 4.2) or vibrations (Section 4.4) are of interest, but also stress evaluations, for example for fatigue analysis. In the last example, we demonstrate the performance of all methods in Table 2 on the evaluation of stresses in a curved Kirchhoff-Love shell. Since the Kirchhoff-Love shell formulation is displacement-based, the displacements are C 1 continuous across patch interfaces for C 1 constructions. The stresses, however, are based on the gradients of the displacements, hence their continuity theoretically is C 0 for a perfect C 1 coupling. In this example, we elaborate on the Von Mises membrane stress field resulting from the 6-patch elliptic paraboloid from Fig. 13. The stress fields are plotted for bases with degree and regularity from Table 2 and additionally for a basis with p = 4, r = 2. Note that the regularity r of these bases is the regularity in the patch interior.

In Fig. 16, the stress fields for the elliptic paraboloid example from Fig. 13 are provided.

From these results, it can immediately be seen that the stress field for a single patch parameterization with basis p = 2, r = 1 exposes the elements of the basis because of the C 0 continuity across elements. Similar effects are seen for the D-patch, Almost-C 1 and the penalty method.

Increasing the degree of the basis while keeping the regularity the same results in a p = 3, r = 1 basis. The element continuity is still C 0 for the stresses, but the higher continuity of the basis within the element results in a slightly improved stress field, as can be seen from the single patch, the D-patch and penalty methods. The Approx. C 1 and AS-G 1 methods in addition show a better stress field around the EVs compared to the D-patch with only small wiggles in the inner contour.

Increasing the smoothness by going to p = 3, r = 2 shows that the Approx. C 1 method predicts the stress field very well over the whole domain but with the wiggles in the inner contour, and that the D-patch suffers from the singularity at the EVs. Lastly, the p = 4, r = 2 plots show that the wiggles in the inner contour are eliminated for the Approx. C 1 and the AS-G 1 methods and that the artifacts of the D-patch around the EV are still there but to a lesser extent. Finally, the results of the penalty method in Fig. 16 show it is able to provide an accurate representation of the stress fields. As seen from Fig. 13, penalty factors α = 1 and α = 10 provide good convergence in the bending energy norm. Indeed, the stress fields for the fixed 64 × 64 element meshes in Fig. 16 confirm that for these penalty factors the stress fields accurately represent the single patch stress fields, despite small artifacts around the EVs for α = 1. For a higher penalty factor of α = 100, the stress fields following from the penalty method are not guaranteed to be accurate, showing the downside of this method.

Overall, the stress analysis for multiple combinations and regularities shows that the Almost-C 1 method is generally unfavourable since it is only applicable for p = 2, r = 1 hence C 0 stress fields, suffering from a lack of continuity over the whole domain. This also makes the D-patch as applicable as the Approx. C 1 method in terms of degree and regularity combinations. Comparing the D-patch with the Approx. C 1 and the AS-G 1 methods, it is shown that the D-patch suffers from the singularity in the EVs when reconstructing stresses, whereas the other two methods are able to recover the stress fields without problems. Moreover, this example has also shown the advantage of smooth unstructured spline constructions for stress analyses, since their continuity across (almost) all of the domain is ensured, contrary to the penalty method. Lastly, this example shows the advantage of IGA in general over lower-order methods like FEA, since the higherdegree bases (e.g. p = 4, r = 2) provides smooth stress fields compared to lower-degree bases (p = 2, r = 1).

Conclusions

In this section a quantitative comparison of the AS-G 1 , the Approx. C 1 , the D-Patch and the Almost-C 1 constructions is provided. The methods have been assessed on different aspects: i) convergence of the biharmonic equation (Section 4.1); ii) convergence of the linear Kirchhoff-Love shell (Section 4.2); iii) eigenvalue spectrum approximation (Section 4.3); iv) application to a large-scale complex geometry (Section 4.4) and; v) the reconstruction of stress fields (Section 4.5). From these analyses, the following conclusions can be drawn: • All methods converge in a theoretical setting to the same solution for the biharmonic equation (Sections 4.1 and 4.2). However, the convergence behaviour of the D-Patch method is sub-optimal and affected by conditioning issues for large meshes. Furthermore, the Approx. C 1 and AS-G 1 methods give worse convergence compared to other methods for the hyperbolic paraboloid shell but good convergence rates for the elliptic paraboloid shell example.

p = 2, r = 1 p = 3, r = 1 p = 3, r = 2 p = 4, r = 2 D-Patch Almost-C 1 Approx. C 1 AS-G 1 σ V M [MPa]
• From a spectral analysis on the biharmonic equation Section 4.3 it can be concluded that there is no best unstructured spline construction. Depending on the degree and regularity, small difference in the eigenvalue spectra are observed between the methods. Comparing with Nitsche's method, however, it is concluded that the unstructured spline constructions considered in this paper perform consistently better. This is also confirmed by the applied modal analysis on the car geometry Section 4.4, where penalty method fails to find accurate eigenfrequencies, possibly because of an unsuitable penalty parameter.

• From the applied modal analysis on a complex geometry, it can also be concluded that the Almost-C 1 and D-Patch constructions are more straight-forward to apply to a complex geometry extracted from a mesh. This is due to the fact that the Approx. C 1 and AS-G 1 constructions require, respectively, a G 2 geometry and an analysis-suitable geometry, which are both not trivial to construct from an originally C 0 -continuous mesh. Instead, the D-Patch and Almost-C 1 constructions require a C 1 geometry, which is easier to construct in general.

• From the stress fields presented in Section 4.5 following from the analysis in Section 4.2, it can be concluded that the AS-G 1 and Approx. C 1 methods provide excellent stress fields.

The D-Patch also provides good stress fields, but inaccuracies are found around the EVs, possibly because of the singularity close to the EV. The Almost-C 1 method is considered inaccurate for stress analysis because of a lack of higher-degree generalizations. Lastly, comparison with penalty methods shows that the unstructured spline constructions generally provide a robust parameter-free approach for coupling, whereas the penalty method requires careful selection of the penalty parameter.

Overall, our finding suggest that the Almost-C 1 and D-Patch are generally easier to construct, but for certain problems they have limited accuracy. On the other hand, the AS-G 1 or Approx. C 1 discretisations require more pre-processing efforts, but provide optimal convergence, hence accuracy. This, however, depends on the input geometry: generic quad-meshes might require more pre-processing efforts than C 1 -matching parameterisations. Lastly, the results provided in this section have shown that strong coupling methods have certain advantages over weak methods, and therefore provide an interesting alternative.

Conclusions and future work

In this paper, we provide a qualitative and quantitative comparison of unstructured spline constructions for smooth multi-patches in isogeometric analysis. The general advantage of unstructured spline constructions over trimming or variational coupling methods is that they are parameter-free, do not require specialized solvers and are typically once constructed in a shape optimization workflow. The goal of this paper is to compare the analysis-suitable G 1 (AS-G 1 ) the approximate C 1 (Approx. C 1 ), the degenerate patches (D-Patch) and the Almost-C 1 constructions with respect to qualitative aspects (i.e. constraints for application) and quantitative aspects (i.e. numerical performance).

From the qualitative analysis, it followed that each method required a different set of constraints to be satisfied before the constructions could be applied, see Fig. 7 and Table 1. Degree and regularity constraints can be satisfied by knot insertion routines or re-fitting, which are relatively straight-forward. The constraint on analysis-suitability for the AS-G 1 and the constraint on G 2 continuity for the Approx. C 1 method require dedicated reparametrisation routines, such as the one presented by [START_REF] Kapl | Construction of analysis-suitable G1 planar multi-patch parameterizations[END_REF]. The fact that D-Patches are restricted to geometries without boundary extraordinary vertices requires redefinition of the quadrilateral mesh. Lastly, the fact that the Almost-C 1 method is only defined for bi-quadratic bases (p = 2) restricts the inter-element continuity to C 1 through the whole domain. Depending on the application and the availability of existing routines in software, different unstructured spline constructions are favourable, depending on the geometric flexibility or desired degree and regularity .

From the quantitative analysis, some conclusions can be drawn on the considered unstructured spline constructions and between unstructured spline constructions compared to variational methods such as Nitsche's method or a penalty method. From the analysis, it was in general observed that depending on the problem type, the different methods have their advantages and disadvantages. Firstly, simple biharmonic equations (see Section 4.1) and linear shells (see Section 4.2) provided good results for all methods. However, the AS-G 1 and Approx. C 1 methods showed slow convergence for the double-curved shell and the D-patch suffered from illconditioning for fine meshes. The Almost-C 1 provided good results in general, however it is only applicable on bi-quadratic splines. Secondly, all methods showed superiority over Nitsche's method for the computation of an eigenvalue spectrum for plate vibrations (see Section 4.3) and no significant differences between the unstructured spline constructions have been observed.

Thirdly, the D-Patch and Almost-C 1 showed straight-forward applicability on the problem of a complex geometry (see Section 4.4), whereas the analysis-suitability requirement of the AS-G 1 method and the smoothness requirement of the Approx. C 1 method are non-trivial to satisfy on off-the-shelf industrial geometries.

For the penalty method, no suitable penalty parameter was found, and probably optimal penalty parameters should be chosen per interface rather than globally. Lastly, the AS-G 1 and Approx. C 1 methods provided superior results for stress reconstruction, where the D-Patch suffered around the EVs due to its singular parameterization and the Almost-C 1 method provided bad results due to a lack of higher degrees.

In conclusion, both comparisons provide an overview of the applicability of the methods with respect to the requirements needed to construct them, on the notions of nestedness and in general on the performance of the methods. Overall, it can be concluded from both analyses that among the compared methods, there is no general best construction. More precisely, the quantitative analysis shows that different methods perform differently in different applications,

given that they can be constructed. Furthermore, with the backgrounds and properties provided in the qualitative analysis section, we hope that the present paper provides valuable insights for application of the considered methods to multi-patch problems.

From

ΩFigure 3 :

 3 Figure3: Given an initial geometry Ω (a), trimming (b) uses the curves of the boundary of the original geometry to define the interior domain Ω int and the exterior domain Ω ext . An alternative approach for modelling the domain is to use domain segmentation (c). Here, the domain is decomposed into several patches Ω i which together define the full domain Ω.

Figure 4 :

 4 Figure 4: Procedure to find a multi-patch segmentation from a given mesh. The original mesh in (a) has 46 vertices, 81 edges and 45 faces and the final multi-patch (c) has 20 patches.

Figure

  Figure5bpresents a local region around and EV with valence five with line styles indicating different continuity levels on patch or element boundaries (see the caption of Fig.5). For the Approx. C 1 construction on a fully smooth basis (p ≥ 3 and r = p -1), the interior basis recovers full smoothness on element boundaries, hence C p-1 continuity. In the shaded region around the interfaces and the EV, the continuity is locally reduced by construction of the locally reduced continuous subspace and the approximation of the gluing data. Similar to the AS-G 1 construction, the continuity on the EV is C 2 by construction and the element boundaries are C 1 approximately.

1 Figure 5 :

 15 Figure 5: Schematic representation of the continuity across element boundaries and patch interfaces for the (a) AS-G 1 construction, (b) Approx. C 1 constructions. Thin lines indicate element boundaries and thick lines indicate patch interfaces. Solid lines represent C p-1 continuity, dashed lines represent C p-2 continuity, thick dashed lines represent C 1 interfaces and loosely dashed lines represent approximate C 1 interfaces. A double lined circle represent a C 2 continuous vertex, a filled circle represent a singular vertex and a white filled circled with a single line represents a C 1 continuous vertex. The gray shaded area for the Approx. C 1 represents local reduced continuity.

1 Figure 6 :

 16 Figure 6: Schematic representation of the continuity across element boundaries and patch interfaces for the (a) D-Patch and (b) Almost-C 1 . Line styles are as in Fig. 5.

Figure 7 :

 7 Figure 7: Inside the unstructured spline pre-processing block from Fig. 2. The unstructured spline requirements are depicted in diamond-shaped blocks for methods AS-G 1 , Approx. C 1 , D-Patch and Almost C 1 .The first row represents requirements on degree p and regularity r. If not satisfied, the geometry can be projected onto a space that satisfies the requirement, or degree elevation or reduction steps can be performed together with refinement operations. The second row depicts the requirements on the geometry parameterization; these blocks can be satisfied by changing the geometry.

Fig. 8

 8 Fig. 8 is used to decompose a simple domain into a domain with extraordinary vertices in the interior. Domains with EVs on the boundary are left out of scope, since the D-patch construction would change the outer boundaries of the domain, hence the comparison would involve a significantly different geometry. Since different methods have different constraints on the degree and regularity of the basis, different combinations of the degree p and regularity r are tested through-

9 .

 9 For degree p = 2 and regularity r = 1 the Almost-C 1 , D-patch and Nitsche coupling methods are compared. As expected, the results show consistency between the Almost-C 1 , D-patch and Nitsche's method with expected convergence. The results also show a slight dependency on the factor β for the D-Patch. For degree p = 3 and regularity r = 1, the Approx. C 1 , AS-G 1 , D-patch and Nitsche's method can be compared. The results of the Approx. C 1 and AS-G 1 are exactly the same, since the original geometry is analysis-suitable and contains only bi-linear patches. Then applying the Approx. C 1 to

Figure 10 :

 10 Figure 10: Hyperbolic paraboloid shell geometry with coordinates r(ξ 1 , ξ 2 ) = ξ 1 ξ 2 ξ 2 1 -ξ 2 2 , ξ 1 , ξ 2 ∈ [-1/2.1/2].The left-edge of the hyperbolic paraboloid is clamped, i.e. the displacements and rotations are zero (u = 0 and

Figure 11 :

 11 Figure 10: Hyperbolic paraboloid shell geometry with coordinates r(ξ 1 , ξ 2 ) = ξ 1 ξ 2 ξ 2 1 -ξ 2 2 , ξ 1 , ξ 2 ∈ [-1/2.1/2].The left-edge of the hyperbolic paraboloid is clamped, i.e. the displacements and rotations are zero (u = 0 and ∂uz ∂x = 0).

  Young's modulus E = 200 [GPa] and Poisson's ratio ν = 0.3 [-]. The refinement procedure as described in Section 4.1 is used for the D-Patch and Almost-C 1 constructions.

  Penalty α = 1 Penalty α = 10 Penalty α = 100 D-Patch Approx. C 1 AS-G 1 Almost-C 1

3 /( 12 ( 1 -

 3121 ν 2 )) is the flexural rigidity of the plate with E = 10 5 [Pa] the Young's modulus of the plate, t = 10 -2 [m] the thickness and ν = 0.2 [-] the Poisson's ratio. Furthermore, ρ = 10 5 [kg] is the material density. Equation (

  Nitsche, α = 10 5 Approx. C1 Nitsche, α = 10 6

Figure 14 :

 14 Figure14: Eigenvalue spectra for the biharmonic eigenvalue problem on the domain from Fig.8. The horizontal axes depict the eigenvalue index i over the total number of eigenvalues N. The vertical axes represent the numerical eigenvalue ω i,h over the analytical eigenvalue ω i , both with index i. The results are plotted for different combinations of the degree p and regularity r of the basis. The results for a Nitsche method are given for different penalty parameters α.

  (a) Original quad mesh with 15895 vertices, 31086 edges and 62172 faces. (b) Interface (green) and boundary (red) curves. (c) Final multi-patch segmentation with 307 patches.

Figure 15 :

 15 Figure 15: Geometry of the side panel of a car. The original mesh (a) is traced with the procedure from Fig. 4, yielding a set of boundary and interface curves (b). From these curves, the multi-patch segmentation (c) for isogeometric analysis is constructed following Fig. 4c.

  for the Almost-C 1 construction and with (p, r) = (2, 1), (p, r) = (3, 1) and (p, r) = (3, 2) for the D-Patch. Figure16provides the corresponding mode shapes on the D-Patch geometry with p = 3, r = 2 and the mode shapes have been qualitatively matched to construct Table

Figure 16 :

 16 Figure 16: Out-of-plane deformations of the first four vibration modes of the side of the car from Fig. 15. The results on the left represent the results obtained by the D-Patch construction and the results on the right represent results obtained using ABAQUS (10mm). The mode shapes are all deformation modes warped by the deformation vector and plotted over the undeformed (transparent) geometry.

Figure 16 : 2 α

 162 Figure 16: (Caption on next page).

Figure 16 :

 16 Figure 16: Von Mises membrane stress fields for the single patch, unstructured splines and penalty-coupled multipatch paraboloid from Section 4.2 and Fig. 13 with 64 × 64 elements per patch. The results are provided for different combinations of degree p and regularity r. The color bar represents the stress and the contours are plotted for stress levels σ V M ∈ {10 5 , 10 6 , 10 7 } [MPa].

  1 continuity, while basis functions at vertices have C 2 continuity. The

(a) A simple mesh with boundary edges in black and interior edges in gray. The boundary extraordinary vertices (bEVs), i.e. the vertices on a boundary with valence ν ≥ 3 are denoted by a black circle and the interior extraordinary vertices (iEVs), i.e. interior vertices with valence ν ≥ 4, ν 4 are denoted by gray circles. (b) Illustration of the interface tracing procedure. From each EV all outgoing edges are traced as illustrated until another EV or a boundary is hit. (c) Result of interface tracing from all the EVs. Every patch is now bounded by a a set of boundary and traced interface curves. All patch corners are corners where a traced interface and/or a boundary edge form a corner. Along the hole, different patches are indicated with different shades of gray. In the part bottomright of the hole, every face forms a patch, since all traced curves denoted by colors intersect with other traced curves.

We believe that a comparison like the one presented in this paper is never fully unbiased, since the authors have contributed to different methods in previous publications and do not represent the entire research community.

The notation S(p = (p, p), r = (r, r), h = (h, h)) indicates a two-dimensional spline space with p as the polynomial degree, r as the regularity and h as the mesh size in both directions.

Table2: Degree p and regularity r constraints for each considered method from Section 3, see Table1.p = 2, r = 1 p = 3, r = 1 p = 3, r = 2 D-patch ⋆ ⋆ ⋆ Almost-C 1 ⋆ Approx. C 1 ⋆ ⋆ AS-G 1 ⋆ Nitsche/Penalty ⋆ ⋆ ⋆for more details. In all examples, Dirichlet boundary conditions are applied at the control points and clamped boundary conditions are applied weakly as in[START_REF] Herrema | A framework for isogeometric-analysis-based optimization of wind turbine blade structures[END_REF]. All results are obtained using the Geometry + Simulation modules[START_REF] Jüttler | Geometry + Simulation Modules: Implementing Isogeometric Analysis[END_REF][START_REF] Mantzaflaris | An overview of geometry plus simulation modules[END_REF] and will be published in a separate publication.As discussed in Section 3, the D-patch and Almost-C 1 constructions involve a pre-smoothing

In the case of a different starting point for this benchmark, such as a smooth mesh composed of higher-order quadrilateral elements, e.g. derived from a subdivision surface, instead of a bi-linear mesh, the pre-processing efforts required for the AS-G 1 and Approx. C 1 will be different.

Acknowledgements HMV is grateful to Delft University of Technology, specifically the faculty of Electrical Engineering, Mathematics and Computer Science (EEMCS) and the faculty of Mechanical, Maritime and Materials Engineering (3mE) for the financial support to conduct this research. AM acknowledges support from EU's Horizon 2020 research under the Marie Sklodowska-Curie grant 860843 and DT is grateful to ANSYS Inc. for their financial support. Furthermore, the authors thank Wei Jun Wong from Delft University of Technology for delivering the ABAQUS reference results, and the community of the Geometry + Simulation Modules (a.k.a. gismo) for their continuous effort on improving the code.

In addition, the comparisons in the present paper provide directions for the improvements of the considered methods. For the AS-G 1 and Approx. C 1 methods, restrictions on geometry and parameterization are a bottleneck in the industrial applications. Therefore, it is recommended to expand the applicability of these methods by developing dedicated geometric pre-processing routines. For the D-Patch construction, the limitation of the construction of the basis near ν > 3 boundary EVs calls for the development of routines to eliminate these EVs in quadrilateral multipatches, as discussed in the qualitative comparison. Furthermore, the example of the bi-harmonic equation has shown that the D-Patch can suffer from ill-conditioned system, hence development of pre-conditioners for D-Patch constructions is advised. Lastly, although the Almost-C 1 resolves the downsides of the D-Patch construction, its restriction on the degree of the spline-space is a major disadvantage when plotting stress fields in shell analysis. Therefore, for the Almost-C 1 construction it is recommended to explore expansion to higher degrees.