domain: info math

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

A Mantzaflaris
email: angelos.mantzaflaris@inria.fr

M Möller
email: m.moller@tudelft.nl

J H Den Besten
email: henk.denbesten@tudelft.nl

Goal-Adaptive Meshing of Isogeometric Kirchhoff-Love Shells

Keywords: Isogeometric Analysis, Buckling, Dual-Weighted Residual Method, Adaptive Meshing, Kirchhoff-Love Shell 1 Solve Estimate Mark Refine Transfer Advance Until converged Addition for solution stepping Adaptive meshing

Mesh adaptivity is a technique to provide detail in numerical solutions without the need to refine the mesh over the whole domain. Mesh adaptivity in isogeometric analysis can be driven by Truncated Hierarchical B-splines (THB-splines) which add degrees of freedom locally based on finer B-spline bases. Labeling of elements for refinement is typically done using residual-based error estimators. In this paper, an adaptive meshing workflow for isogeometric Kirchhoff-Love shell analysis is developed. This framework includes THB-splines, mesh admissibility for combined refinement and coarsening and the Dual-Weighted Residual (DWR) method for computing element-wise error contributions. The DWR can be used in several structural analysis problems, allowing the user to specify a goal quantity of interest which is used to mark elements and refine the mesh. This goal functional can involve, for example, displacements, stresses, eigenfrequencies etc. The proposed framework is evaluated through a set of different benchmark problems, including modal analysis, buckling analysis and non-linear snap-through and bifurcation problems, showing high accuracy of the DWR estimator and efficient allocation of degrees of freedom for advanced shell computations.

Introduction

The idea behind isogeometric analysis (IGA) [START_REF] Hughes | Isogeometric analysis: CAD, finite elements, NURBS, exact geometry and mesh refinement[END_REF] is to bridge the gap between computer aided design (CAD) and finite element analysis (FEA). By employing B-splines or Non-uniform rational B-splines (NURBS) or ordinary B-splines as the basis for FEA, IGA not only provides geometrically exact analysis, but the high smoothness of the spline bases also provides high accuracy per degree of freedom [START_REF] Sande | Explicit error estimates for spline approximation of arbitrary smoothness in isogeometric analysis[END_REF]. The close link with conventional engineering fields such as automotive, offshore, aircraft or civil engineering makes structural analysis with isogeometric analysis a particular field of interest. Besides the performance of the different isogeometric element formulations for thin Kirchhoff-Love shells [START_REF] Kiendl | Isogeometric shell analysis with Kirchhoff-Love elements[END_REF][START_REF] Kiendl | Isogeometric Kirchhoff-Love shell formulations for general hyperelastic materials[END_REF][START_REF] Verhelst | Stretch-Based Hyperelastic Material Formulations for Isogeometric Kirchhoff-Love Shells with Application to Wrinkling[END_REF][START_REF] Alaydin | An updated Lagrangian framework for Isogeometric Kirchhoff-Love thin-shell analysis[END_REF], moderately thick Reissner-Mindlin shells [START_REF] Benson | A large deformation, rotation-free, isogeometric shell[END_REF][START_REF] Hu | Isogeometric analysis of thin Reissner-Mindlin shells: locking phenomena and B-bar method[END_REF][START_REF] Kiendl | Isogeometric collocation methods for the Reissner-Mindlin plate problem[END_REF][START_REF] Sobota | Implicit dynamic analysis using an isogeometric Reissner-Mindlin shell formulation[END_REF][START_REF] Benson | Blended isogeometric shells[END_REF] or thicker solid-like shells [START_REF] Hosseini | An isogeometric solidlike shell element for nonlinear analysis[END_REF][START_REF] Leonetti | An efficient isogeometric solidshell formulation for geometrically nonlinear analysis of elastic shells[END_REF] in static and dynamic simulations, conventional engineering disciplines also rely on accurate modal and (post-)buckling simulations. In addition, the ability to handle complex (multipatch) CAD geometries via trimming [START_REF] Coradello | Hierarchically refined isogeometric analysis of trimmed 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] Leidinger | Explicit dynamic isogeometric B-Rep analysis of penalty-coupled trimmed NURBS shells[END_REF] or patch coupling methods [START_REF] Herrema | Penalty coupling of non-matching isogeometric Kirchhoff-Love shell patches with application to composite wind turbine blades[END_REF][START_REF] Leonetti | A robust penalty coupling of non-matching isogeometric Kirchhoff-Love shell patches in large deformations[END_REF][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][START_REF] Guo | Nitsche's method for a coupling of isogeometric thin shells and blended shell structures[END_REF][START_REF] Coradello | Coupling of non-conforming trimmed isogeometric Kirchhoff-Love shells via a projected super-penalty approach[END_REF] improves the applicability of IGA in structural engineering. For problems with a large number of degrees of freedom or problems with a large number of load/time steps, mesh adaptivity can play a key role in providing efficient simulations for industrial applications.

A loop in an adaptive isogeometric method (AIGM) consists of the steps solve Partial Differential Equation (PDE) at hand, estimate element-wise error contributions, mark elements for refinement, refine (coarsen) marked elements [START_REF] Buffa | Adaptive isogeometric methods with hierarchical splines: Error estimator and convergence[END_REF], see Fig. 1. This process can be repeated in an iterative manner (e.g. for static, buckling or modal analysis), until satisfactory accuracy is achieved or it can be applied (iteratively) within a time/load stepping procedure. A broad overview of the mathematical foundations of AIGMs is given in [START_REF] Buffa | Mathematical Foundations of Adaptive Isogeometric Analysis[END_REF]. In previous works, AIGMs are developed for different applications (solve), using different estimation strategies, marking strategies and often for mesh refinement (coarsening), with a few also providing coarsening strategies [START_REF] D'angella | Efficient extraction of hierarchical B-Splines for local refinement and coarsening of Isogeometric Analysis[END_REF][START_REF] Carraturo | Suitably graded THBspline refinement and coarsening: Towards an adaptive isogeometric analysis of additive manufacturing processes[END_REF][START_REF] Lorenzo | Hierarchically refined and coarsened splines for moving interface problems, with particular application to phase-field models of prostate tumor growth[END_REF][START_REF] Hennig | Projection and transfer operators in adaptive isogeometric analysis with hierarchical B-splines[END_REF][START_REF] Garau | Algorithms for the implementation of adaptive isogeometric methods using hierarchical B-splines[END_REF].

Solve The solve block contains the partial differential equation (PDE) at hand. It can be a physics-based problem, e.g. to solve shell [START_REF] Coradello | Hierarchically refined isogeometric analysis of trimmed shells[END_REF][START_REF] Antolin | A hierarchical approach to the a posteriori error estimation of isogeometric Kirchhoff plates and Kirchhoff-Love shells[END_REF][START_REF] Coradello | Hierarchically refined isogeometric analysis of trimmed shells[END_REF], linear elasticity [START_REF] Verhoosel | Imagebased goal-oriented adaptive isogeometric analysis with application to the micromechanical modeling of trabecular bone[END_REF] or freesurface flow [START_REF] Kuru | Goal-adaptive Isogeometric Analysis with hierarchical splines[END_REF] problems. Also, the solve step involve a non-physics PDE, e.g. for mesh generation [START_REF] Hinz | Goal-Oriented Adaptive THB-Spline Schemes for PDE-Based Planar Parameterization[END_REF]. Estimate Determination of element-wise errors is done in the estimate block. In the works [START_REF] Coradello | Hierarchically refined isogeometric analysis of trimmed shells[END_REF][START_REF] Antolin | A hierarchical approach to the a posteriori error estimation of isogeometric Kirchhoff plates and Kirchhoff-Love shells[END_REF][START_REF] Coradello | Hierarchically refined isogeometric analysis of trimmed shells[END_REF], an error estimator based on a residual-like variational problem in the so-called bubble-space was presented for Kirchhoff plates, Kirchhoff-Love shells and trimmed domains. This method has proven a large decrease in CPU time compared to a residual-based error estimator in the strong form, due to its easy parallelisation and the small block-structure of the linear system to solve. As an alternative to this method, error estimation can also be performed in a goal-oriented fashion, by the Dual-Weighted Residual (DWR) method. This method has been applied in the FEA context in various works [START_REF] Rannacher | Adaptive finite element methods in flow computations[END_REF][START_REF] Gedicke | A posteriori error estimators for convection-diffusion eigenvalue problems q[END_REF][START_REF] Hartmann | Error estimation and adaptive mesh refinement for aerodynamic flows[END_REF][START_REF] Hartmann | Adaptive discontinuous Galerkin finite element methods for the compressible Euler equations[END_REF][START_REF] Hartmann | Adaptive discontinuous Galerkin finite element methods for nonlinear hyperbolic conservation laws[END_REF][START_REF] Hartmann | Adaptive discontinuous Galerkin finite element methods for the compressible Euler equations[END_REF][START_REF] Möller | Adaptive mesh refinement for high-resolution finite element schemes[END_REF][START_REF] Cliffe | Adaptive discontinuous galerkin methods for eigenvalue problems arising in incompressible fluid flows[END_REF] and was used in the works [START_REF] Kuru | Goal-adaptive Isogeometric Analysis with hierarchical splines[END_REF][START_REF] Van Der Zee | Isogeometric analysis-based goal-oriented error estimation for free-boundary problems[END_REF] for Poisson and free-boundary problems, in [START_REF] Dedè | B-spline goal-oriented error estimators for geometrically nonlinear rods[END_REF] for a geometrically non-linear rod, in [START_REF] Hinz | Goal-Oriented Adaptive THB-Spline Schemes for PDE-Based Planar Parameterization[END_REF] PDE-based domain parametrisations and in [START_REF] Verhoosel | Imagebased goal-oriented adaptive isogeometric analysis with application to the micromechanical modeling of trabecular bone[END_REF] for micromechanical modelling of trabecular bone. Fig. 1: A typical flowchart for an adaptive meshing routine. The classical solution stepping depicts a process without adaptive meshing. Here, a solution is obtained by the solve and the solution is advanced (e.g. in time or load step) and recomputed. The adaptive meshing step denotes the additional operations for mesh adaptivity and the Addition for solution stepping includes an additional transfer step in case the adaptive meshing method is applied to solution. The Estimate block provides an error estimation with local contributions per element or per degree of freedom (DoF). The Mark block contains a marking rule that marks elements for refinement based on a specific rule. The Refine block transforms the current mesh to a new mesh, where elements are refined and coarsened based on the marking rule. The block Transfer transfers the previous solution to the new mesh, so that it can be used to recompute the present interval on a modified mesh. This recomputation is performed again in the Solve block and follows through the subsequent blocks, until an adaptivity criterion is reached. For example, a criterion that requires the total error in the mesh to be within certain bounds.

Goal-oriented refinement in general provides element-wise error estimates by solving a linear adjoint problem on the current space and an enriched space. Mark As soon as element-wise error contributions are known, elements can be marked for refinement. This marking is mostly done using using the Dörfler marking strategy [START_REF] Dörfler | A Convergent Adaptive Algorithm for Poisson's Equation[END_REF], as in [START_REF] Buffa | Adaptive isogeometric methods with hierarchical splines: Error estimator and convergence[END_REF][START_REF] Carraturo | Suitably graded THBspline refinement and coarsening: Towards an adaptive isogeometric analysis of additive manufacturing processes[END_REF][START_REF] Gedicke | A posteriori error estimators for convection-diffusion eigenvalue problems q[END_REF], which involves marking the largest elements until their sum exceeds a certain percentage of the total error. An alternative is to mark the elements with an error higher than a threshold (an absolute threshold based on the maximum error) [START_REF] Antolin | A hierarchical approach to the a posteriori error estimation of isogeometric Kirchhoff plates and Kirchhoff-Love shells[END_REF][START_REF] Giannelli | THB-splines: An effective mathematical technology for adaptive refinement in geometric design and isogeometric analysis[END_REF] or based on a relative threshold taking a fixed percentage of the total number of cells for refinement. In [START_REF] Giannelli | THB-splines: An effective mathematical technology for adaptive refinement in geometric design and isogeometric analysis[END_REF] the latter two strategies are discussed. Refine Local refinement for adaptive meshing in isogeometric analysis is enabled by the use of Hierarchical B-splines (HB-splines) [START_REF] Vuong | A hierarchical approach to adaptive local refinement in isogeometric analysis[END_REF], Truncated Hierarchical B-splines (THB-splines) [START_REF] Giannelli | THB-splines: An effective mathematical technology for adaptive refinement in geometric design and isogeometric analysis[END_REF][START_REF] Giannelli | THB-splines: The truncated basis for hierarchical splines[END_REF] or T-splines [START_REF] Bazilevs | Isogeometric analysis using Tsplines[END_REF] amongst other spline constructions, which are reviewed in [START_REF] Hennig | Bézier extraction and adaptive refinement of truncated hierarchical NURBS[END_REF]. HB-splines provides a nested, linear dependent space that violates the partition of unity property. To preserve the latter, THB-splines have been introduced in [START_REF] Giannelli | THB-splines: The truncated basis for hierarchical splines[END_REF]. For (T)HB-splines suitable grading to generate admissible meshes should be taken into account in order to guarantee a bounded error [START_REF] Buffa | Adaptive isogeometric methods with hierarchical splines: Error estimator and convergence[END_REF], for which algorithms have been presented in [START_REF] Bracco | Refinement algorithms for adaptive isogeometric methods with hierarchical splines[END_REF]. In the work of [START_REF] Hennig | Adaptive mesh refinement strategies in isogeometric analysis-A computational comparison[END_REF], a distinction is made between greedy and safe refinement, the former being a refinement of cells with a 1-level difference with adjacent cells and the latter being a refinement complying with the refinement neighborhoods defined in [START_REF] Buffa | Adaptive isogeometric methods with hierarchical splines: Error estimator and convergence[END_REF]. Besides for adaptive meshing for solving PDEs [START_REF] D'angella | Efficient extraction of hierarchical B-Splines for local refinement and coarsening of Isogeometric Analysis[END_REF][START_REF] Carraturo | Suitably graded THBspline refinement and coarsening: Towards an adaptive isogeometric analysis of additive manufacturing processes[END_REF][START_REF] Lorenzo | Hierarchically refined and coarsened splines for moving interface problems, with particular application to phase-field models of prostate tumor growth[END_REF][START_REF] Hennig | Projection and transfer operators in adaptive isogeometric analysis with hierarchical B-splines[END_REF][START_REF] Garau | Algorithms for the implementation of adaptive isogeometric methods using hierarchical B-splines[END_REF][START_REF] Hinz | Goal-Oriented Adaptive THB-Spline Schemes for PDE-Based Planar Parameterization[END_REF], THB-splines have also been succesfully applied in the context of fitting [START_REF] Bracco | THB-Spline Approximations for Turbine Blade Design with Local B-Spline Approximations[END_REF][START_REF] Kiss | Adaptive CAD model (re-)construction with THB-splines[END_REF]. Transfer The transfer from previous time/load-steps onto a new mesh can be done using different methodologies. In the work of [START_REF] Hennig | Projection and transfer operators in adaptive isogeometric analysis with hierarchical B-splines[END_REF] different least-squares approaches are provided. Furthermore, quasi-interpolation [START_REF] Speleers | Effortless quasi-interpolation in hierarchical spaces[END_REF][START_REF] Giust | Local (t)HB-spline projectors via restricted hierarchical spline fitting[END_REF] is a technique that can be used to transfer solutions between hierarchical meshes.

In this paper, we employ goal-oriented adaptive refinement for isogeometric thin shell analysis to facilitate THB-adaptive meshing for a variation of structural analysis problems. The contribution of the paper is threefold. Firstly, we use the Dual-Weighted Residual (DWR) method to derive novel error estimators for structural shell analysis, given goal functionals based on displacements, (principal) stresses and strains, forces and moments. These goal functionals are applied to isogeometric Kirchhoff-Love shells -since it provides a natural separation of bending and membrane terms -but the general framework and the goal functionals presented are easily adapted for other shell formulations. Secondly, we employ the eigenvalue DWR from [START_REF] Rannacher | Adaptive finite element methods in flow computations[END_REF][START_REF] Lathouwers | Spatially adaptive eigenvalue estimation for the SN equations on unstructured triangular meshes[END_REF] for error estimations for modal and buckling analyses. Lastly, the goal functionals are used to drive an adaptive meshing strategy with suitable grading and efficient transfer of solutions by quasi-interpolation method on hierarchical spline spaces [START_REF] Speleers | Effortless quasi-interpolation in hierarchical spaces[END_REF][START_REF] Giust | Local (t)HB-spline projectors via restricted hierarchical spline fitting[END_REF]. This adaptive meshing strategy is applied to non-linear shell analysis with focus on buckling problems with snap-through and bifurcation instabilitiesbeing new applications in the realm of adaptive meshing research for shells.

The paper is structured as follows. In Section 2, the isogeometric Kirchhoff-Love shell analysis proposed by [START_REF] Kiendl | Isogeometric shell analysis with Kirchhoff-Love elements[END_REF] is briefly revised and some basic concepts for structural analysis computations are given. In Section 3, the Dual-Weighted Residual (DWR) method is provided for the isogeometric Kirchhoff-Love shell using the membrane and flexural strain split. However, it can be used for general elasticity problems. Moreover, the section provides the DWR method for eigenvalue problems to compute error estimators for modal and buckling analyses. Thereafter, Section 4 provides the details for adaptivity for isogeometric analysis. This includes the concept of Truncated-Hierarchical B-splines (THB-splines) and admissible refinement. Furthermore, the mark and transfer operations are described. In Section 5, a summary of the preceeding sections is provided by means of a global algorithm for the AIGM for structural analysis computations with load-stepping. In Section 6 the present work is evaluated on numerical benchmark problems, ranging from linear problems with analytical solutions to non-linear shell problems. Finally, Section 7 provides conclusions and and outlook based on this work.

Isogeometric Kirchhoff-Love Shell Analysis

In this section, we provide a brief background on the Kirchhoff-Love shell formulation. For more details on this formulation, we refer to [3-5, 57, 58].

Shell Kinematics

Since Kirchhoff-Love shells satisfy the Kirchhoff Hypothesis [START_REF] Reddy | An Introduction to Nonlinear Finite Element Analysis: With Applications to Heat Transfer[END_REF], the coordinates x of any parametric point θ = (θ 1 , θ 2 , θ 3 ) in the shell surface can be represented by the surface position r(θ 1 , θ 2 ) and contribution in normal direction θ 3 a 3 as

x(θ) = r(θ 1 , θ 2 ) + θ 3 a 3 , (1)
Given the covariant basis of the surface r, defined by a α , α = 1, 2 and the orthogonal unit normal a 3 , the covariant basis of x is defined as follows:

g α = x ,α = a α + θ 3 a 3,α , g 3 = x ,3 = a 3 . (2)
Given the second fundamental form b αβ = a 3 • a α,β = -a 3,β • a α , the metric metric coefficients defined as

g αβ = g α • g β = a αβ -2θ 3 b αβ , (3)
and the contravariant basis vectors g α can simply be obtained by g α = g αβ g β . The undeformed configuration r and the deformed configuration r of the surface are related by r = r + u. From the defintion of the deformation gradient F = g i ⊗ gi , the deformation tensor C can be obtained:

C = F F = g i • g j gi ⊗ gj = g ij gi ⊗ gj . (4)
Note that the deformation tensor is defined in the contravariant undeformed basis gi ⊗ gj . For Kirchhoff-Love shells, it is known that g α3 = g 3α = 0, hence this implies C α3 = C 3α = 0. Since g 33 = 1, which implies C 33 to be unity and meaning that the thickness remains constant under deformation. As a result, the Green-Lagrange strain tensor E = E αβ gα ⊗g β and its decomposition to membrane and bending contributions (ε and κ, respectively) is [START_REF] Kiendl | Isogeometric shell analysis with Kirchhoff-Love elements[END_REF][START_REF] Kiendl | Isogeometric Kirchhoff-Love shell formulations for general hyperelastic materials[END_REF]:

E αβ = 1 2 (g αβ -gαβ ) = 1 2 (a αβ -åαβ ) -2θ 3 b αβ -bαβ = ε αβ + θ 3 κ αβ .
(5)

Constitutive relation

The constitutive relations for the Kirchhoff-Love shell relate the Green-Lagrange strain tensor E to the second Piola-Kirchhoff stress tensor S. For linear elastic materials, this is achieved by:

S αβ = C αβγδ E γδ (6)
where C = C αβγδ gi ⊗ gj ⊗ gk ⊗ gl is the material tensor, which takes for linear materials the form C αβγδ = 4 λµ λ+2µ gαβ gγδ + 2µ gαδ gβγ + gαγ gβδ [START_REF] Goyal | Isogeometric Shell Discretizations for Flexible Multibody Dynamics[END_REF]. For nonlinear hyperelastic constitutive relations, the stress and material tensors are derived from the 3D constitutive relations for (in)compressible materials and due to through-thickness integration, the shell normal force and bending moment tensors n = n αβ gα ⊗g β and m = m αβ gα ⊗g β , respectively, are defined as

n αβ (u) = [-t/2,t/2] S αβ (u) dθ 3 , m αβ (u) = [-t/2,t/2] θ 3 S αβ (u) dθ 3 . (7)
For more details on hyperelastic material models, the reader is referred to [START_REF] Kiendl | Isogeometric Kirchhoff-Love shell formulations for general hyperelastic materials[END_REF][START_REF] Roohbakhshan | Efficient isogeometric thin shell formulations for soft biological materials[END_REF] and specifically for stretch-based ones to [START_REF] Verhelst | Stretch-Based Hyperelastic Material Formulations for Isogeometric Kirchhoff-Love Shells with Application to Wrinkling[END_REF].

Variational Formulation

The shell internal and external equilibrium equations in variational form are derived by the principle of virtual work [START_REF] Kiendl | Isogeometric shell analysis with Kirchhoff-Love elements[END_REF][START_REF] Kiendl | Isogeometric Kirchhoff-Love shell formulations for general hyperelastic materials[END_REF]. The weak formulation follows from the principle of virtual work with virtual displacements φ:

Find u ∈ S s.t.W(u, φ) := δW int -δW ext = Ω n(u) : ε (u, φ) + m(u) : κ (u, φ) dΩ - Ω f (u) • φ dΩ , ∀φ ∈ S (8)
With f (u) the surface load acting on the mid-surface, for the sake of generality defined as a function of the displacements u (e.g. a follower pressure p gives f (u) = pa 3 (u)). Furthermore, ε (u, φ) and κ (u, φ) are the virtual strain components given displacements u and being linear with respect to variation φ, hence W(u, φ) is also linear in its second argument. The coefficients of the variations of the membrane force and bending moment tensors are

(n ) αβ (u, φ) = [-t/2,t/2] C αβγδ (u) dθ 3 ε γδ (u, φ) + [-t/2,t/2] θ 3 C αβγδ (u) dθ 3 κ γδ (u, φ), (m ) αβ (u, φ) = [-t/2,t/2] θ 3 C αβγδ (u) dθ 3 ε γδ (u, φ) + [-t/2,t/2] θ 3 2
C αβγδ (u) dθ 3 κ γδ (u, φ).

(9) Linearizing the virtual work from Eq. ( 8) provides the continuous equivalent of the Jacobian or tangential stiffness matrix for Newton iterations which will be performed to solve the nonlinear weak formulation Eq. ( 8) in a discrete setting [START_REF] Sauer | On the theoretical foundations of thin solid and liquid shells[END_REF]: [START_REF] Sobota | Implicit dynamic analysis using an isogeometric Reissner-Mindlin shell formulation[END_REF] where ε (u, φ, ψ) and κ(u, φ, ψ) are the second variations of the membrane and bending strains and where f is the first variation of the applied force, being nonzero when the force is depending on the displacements u. For the details on these formulations, we refer to previous publications [START_REF] Kiendl | Isogeometric shell analysis with Kirchhoff-Love elements[END_REF][START_REF] Sauer | On the theoretical foundations of thin solid and liquid shells[END_REF][START_REF] Goyal | Isogeometric Shell Discretizations for Flexible Multibody Dynamics[END_REF]. It should be noted that in the undeformed case, u = 0, the internal membrane forces and bending forces, n(u) and m(u), respectively, vanish. As a result, the continuous equivalent for the linear stiffness matrix is:

W (u, φ, ψ) := Ω n (u, ψ) : ε (u, φ) + n(u) : ε (u, φ, ψ) + m (u, ψ) : κ (u, φ) + m : κ (u, φ, ψ) dΩ - Ω f (u, ψ) • φ dΩ
W (φ, ψ) = Ω n (ψ) : ε (φ) + m (ψ) : κ (φ) dΩ - Ω f (ψ) • φ dΩ (11)
Where the • denotes tensors and functions on the undeformed geometry, i.e. with u = 0.

In our implementation, the tangent stiffness matrix is computed using appropriate Gauss-Lengendre quadrature for each element in the hierarchical mesh. We note that more efficient numerical integration approaches exist [START_REF] Pan | Efficient matrix assembly in isogeometric analysis with hierarchical b-splines[END_REF][START_REF] Pan | Efficient matrix computation for isogeometric discretizations with hierarchical b-splines in any dimension[END_REF][START_REF] Giannelli | Weighted quadrature for hierarchical b-splines[END_REF] for hierarchical splines that might further reduce the computational cost.

Structural Analysis

In the present paper, we will provide different goal functionals for different structural analysis applications. Therefore, we briefly recall the different structural analysis types. Firstly, in case of static analysis, the problem as in Eq. ( 8) is solved. In case of quasistatic analysis, load and/or displacement steps are performed successively and in each step, a static solve is performed. Typically, one writes Eq. ( 8) for load-control as

Find u ∈ S s.t.W(u, φ, λ) := δW int -λδW ext = Ω n(u) : ε (u, φ) + m(u) : κ (u, φ) dΩ - Ω λf 0 • φ dΩ , ∀φ ∈ S, λ ∈ R ( 12
)
where λ is the load factor scaling the reference load f 0 . Quasi-static simulations can be solved using simple load or displacement controlled schemes, using arc-length continuation such as Riks' method or Crisfield's method [START_REF] Riks | The Application of Newton's Method to the Problem of Elastic Stability[END_REF][START_REF] Crisfield | A Fast Incremental/Iterative Solution Procedure That Handles[END_REF]. When quasi-static analysis is performed for post-buckling analysis, one or multiple bifurcation points are passed by definition. On a bifurcation point, the determinant of the tangential stiffness matrix K is equal to zero, hence this matrix is singular. To cope with instabilities, a priori perturbations can be applied to the geometry, or a procedure for approximating singular points [START_REF] Wriggers | A quadratically convergent procedure for the calculation of stability points in finite element analysis[END_REF] can be used. In our previous work, we provide more details on arc-length continuation for post-buckling analysis without providing a priori perturbations [START_REF] Verhelst | Equilibrium Path Analysis Including Bifurcations with an Arc-Length Method Avoiding A Priori Perturbations[END_REF].

In the case of modal analysis and buckling analysis, a generalized eigenvalue problem needs to be solved. These eigenvalue problems have the general form

find (µ, v) ∈ R × S s.t. A(v, φ) = µB(v, φ) ∀φ ∈ S ( 13
)
Where µ is the provides the eigenfrequency in modal analysis and the critical load factor in buckling analysis and where v denotes the vibration or buckling mode shape. The operators A and B are bi-linear. For buckling analysis, A(v, φ) = W (u L , v, φ) and B(v, φ) = W (0, v, φ) with M with u L the pre-buckling solution given load λ L . For modal analysis, A(v, φ) = W (0, v, φ) and B(v, φ) = M(v, φ) with M the mass operator:

M(v, φ) = Ω ρvφ dΩ ( 14
)
Where ρ is the density function over the surface.

Dual-Weighted Residual method

This section elaborates on the Dual-Weighted Residual (DWR) method, which is used in the Estimate step of Fig. 1. The DWR is a method to compute the a posteriori error of a solution in terms of a given goal functional of interest, by solving a linear dual problem. The DWR provides a global estimate of the error, but given a partition of unity of the spline space, it can be used to provide an error contribution per basis function.

General Framework

The general framework of the dual weighted residual (DWR) is presented by [START_REF] Hartmann | Adaptive discontinuous Galerkin finite element methods for the compressible Euler equations[END_REF]. For the sake of completeness, we provide a brief overview of the DWR here. Consider the following nonlinear problem to solve

find u ∈ S s.t. W(u, φ) = 0 ∀φ ∈ S (15)
Where W(u) is a semi-linear operator, u is the solution and φ is a test function and S is a suitably chosen vector space including u ∈ S. The approximation of u, denoted by u h can be found by solving the discrete counterpart of Eq. ( 15)

find u h ∈ S p h s.t. W(u h , φ h ) = 0 ∀φ ∈ S p h . ( 16
)
Where u h , φ h are the discrete counterparts of u and φ and the space S p h ⊂ S is a function space on the (isogeometric) mesh T p h (Ω) with mesh size h and order p covering the computational domain Ω. The solution to this problem is typically obtained by iteratively solving

find φ h ∈ S p h s.t. W (u h , φ h , ψ h ) = R(u h , ψ h ) ∀ψ h ∈ S p h . (17)
While updating the discrete solution. From Eqs. ( 15) and ( 16), and the fact that φ h ∈ S p h ⊂ S the Galerkin orthogonality is found:

W(u, φ h ) -W(u h , φ h ) = W(u -u h , φ h ) = 0. ( 18
)
Where the last equality follows from the fact that the operator W(•, •) is linear in its second argument. Let us now define a nonlinear and differentiable goal functional L(u) or quantity of interest, such that

∆L(u h ) = L(u) -L(u h ). ( 19
)
And let the residual be defined as

R(u h , φ h ) = -W(u h , φ h ). ( 20
)
If ξ ∈ S p is the solution to the following dual problem with L (u, ξ) being the linearisation of the goal functional L(u)

find ξ ∈ S p s.t. W (u, ζ, ξ) = L (u, ζ) ∀ζ ∈ S p , (21)
with discrete counterpart

find ξ h ∈ S p h s.t. W (u h , ζ h , ξ h ) = L (u h , ζ h ) ∀ζ h ∈ S p h . (22)
Then the error representation formula can be obtained (see Theorem 4.1 of [START_REF] Hartmann | Adaptive discontinuous Galerkin finite element methods for the compressible Euler equations[END_REF]):

∆L(u h ) = R(u h , ξ -ξ h ), (23)
where ξ ∈ S is the exact dual solution and ξ h ∈ S p h is the discrete dual solution. This solution can be found by solving the adjoint problem on an enriched basis, i.e.

find ξh ∈ Sp h s.t. W (u h , ζh , ξh ) = L (u h , ζh ) ∀ ζh ∈ Sp h , (24)
with ξh and ζh the dual solution and test functions on the enriched space Sp h , respectively. A choice for Sp h is to use the same mesh as for S p h , with the same regularity but with a higher degree, i.e. Sp h = S p+1 h . This is easily achieved using spline bases. When using B-Splines, one can repeat all knots of the knot vector an extra time compared to the original basis, such that S p h ⊂ S p+1 h ⊂ S is a nested space.

Finally, using Eq. ( 32) together with the dual solution ξ h ∈ S p h and the enriched dual solution ξh ∈ Sp h , an estimate for the global error with respect to the goal functional L can be obtained. To obtain the local element-wise error estimations r i for element ω i ∈ T p h (Ω), such that

∆L(u h ) = R(u h , ξh -ξ h ) = ωi∈T p h (Ω) r i . (25)
Element-wise integration of Eq. ( 23) is simply performed to obtain r i . However, as discussed in Section 4.3 it can be beneficial to have strictly positive element error contributions for element labeling. One can either take the absolute values of r i or one can integrate the squared norm of the integrant in Eq. ( 23) to ensure positivity of element error contributions. Obviously, the sum of the element errors would not be equal to ∆L.

For Kirchhoff-Love shells specifically, the operator W(u, φ) and its linearisation W (u, φ, ψ) are used to perform the DWR analysis.

Eigenvalue problems

When the problem of interest is an eigenvalue problem, the DWR routine is slightly different. Here, we follow the works [START_REF] Gedicke | A posteriori error estimators for convection-diffusion eigenvalue problems q[END_REF][START_REF] Cliffe | Adaptive discontinuous galerkin methods for eigenvalue problems arising in incompressible fluid flows[END_REF][START_REF] Lathouwers | Spatially adaptive eigenvalue estimation for the SN equations on unstructured triangular meshes[END_REF][START_REF] Becker | An optimal control approach to a posteriori error estimation in finite element methods[END_REF][START_REF] Giani | Benchmark results for testing adaptive finite element eigenvalue procedures[END_REF][START_REF] Bangerth | Adaptive Finite Element Methods for Differential Equations[END_REF] to give a brief overview of the DWR for eigenvalue problems.

Let us consider the following eigenvalue problem

find (µ, v) ∈ R × S s.t. A(v, φ) = µB(v, φ) ∀φ ∈ S. (26)
Here, A and B are bi-linear operators. Typically, discretizing the system gives the following:

find (µ h , v h ) ∈ R × S p h s.t. A(v h , φ h ) = µB(v h , φ h ) ∀φ h ∈ S p h , (27)
where the eigenpairs vh = (µ h , v h ) are the solutions of the eigenvalue problem. For uniqueness of the problem, the discrete eigenvectors v h is normalized by the condition

[56] B(v h , v h ) = 1. (28)
To derive the DWR method for the eigenvalue problem in Eq. ( 26), the functional V(•, •), such that the following problem should be solved

Find v = (µ, v) ∈ R × S s.t. V(v, φ) = µB(v, φ) -A(v, φ) + τ (B(v, v) -1) = 0, ∀ φ = (τ, φ) ∈ R × S, (29)
where the normalisation condition from Eq. ( 28) is enforced weakly. The discrete counterpart of this equation reads:

Find vh = (µ h , v h ) ∈ R × S p h s.t. V(v h , φh ) = µB(v h , φ h ) -A(v h , φ h ) + τ h (B(v h , v h ) -1) = 0, ∀ φh = (τ h , φ h ) ∈ R × S p h . (30)
Furthermore, a goal-function for the eigenvalues is defined as follows:

L(v) = µ = µB(v, v), (31)
giving

∆L( φh ) = µ -µ h (32)
Using the nonlinear functional V and the goal functional L, the same derivations as in Section 3.1 can be followed to find a system of equations to solve the DWR eigenvalue problem. The Gateaux derivative of V, denoted by V is given by:

V (v, φ, ψ) = ηB(v, ψ) + µB(ψ, φ) -A(ψ, φ) + τ (B(v, ψ) + B(ψ, v)), (33)
where the derivatives A (ψ, φ) and B (ψ, φ) are equal to the bi-linear operators A(u, φ) and B(u, φ) themselves. Furthermore, the solution around which the linerisation is performed is denoted by v = (µ, v), the test functions are denoted by φ = (τ, φ) and the trial functions are denoted by ψ = (η, ψ). Furthermore, the linearisation of the goal functional Eq. ( 31) is

L (v, ψ) = ηB(v, v) + µ[B(v, ψ) + B(ψ, v)], (34)
such that the adjoint eigenvalue problem becomes

Find φ = (τ, φ) ∈ R × S s.t. ηB(v, ψ) + µB(ψ, φ) -A(ψ, φ) + τ (B(v, ψ) + B(ψ, v)) = ηB(v, v) + µ[B(v, ψ) + B(ψ, v)], ∀ ψ = (η, ψ) ∈ R × S, (35)
which is equivalent to solving [START_REF] Lathouwers | Spatially adaptive eigenvalue estimation for the SN equations on unstructured triangular meshes[END_REF] Find

φ = (τ, φ) ∈ R × S s.t. µB(ψ, φ) -A(ψ, φ) + η[B(v, ψ) -B(v, v)] + (τ -µ)[B(v, ψ) + B(ψ, v)], ∀ ψ = (η, ψ) ∈ R × S. (36)
Assuming the operator B is self-adjoint, it follows that the adjoint eigenvalue problem is solved by η = µ and when normalizing the adjoint eigenvectors ψ using

B(v, ψ) = 1. ( 37
)
It can be seen that the adjoint solution (η, ψ) is obtained by solving the eigenvalue problem

find (η, ψ) ∈ R × S s.t. A(ψ, φ) = µB(ψ, φ) ∀(τ, φ) ∈ R × S, (38)
using Eqs. ( 15), ( 23) and [START_REF] Kuru | Goal-adaptive Isogeometric Analysis with hierarchical splines[END_REF] with W = V according to Eq. ( 29) and with ψ denoting the dual eigenvector and η the dual eigenvalue, the error estimation according to the DWR method for an eigenvalue problem is

∆L(v h ) = A(v h , ψ -ψ h ) -µ h B(v h , ψ -ψ h ) + (η -η h )(B(v h , v h ) -1), ( 39
)
for vh = (µ h , v h ) ∈ R × S p h , ψh = (η h , ψ h ) ∈ R × S p h and ψ = (η, ψ) ∈ R × S .
The exact adjoint solution ψh is again approximated by solving Eq. ( 38) on an enriched space Sp h ⊂ S, Sp h ⊃ S p h , providing (η h , ψh ) ∈ R × Sp h . In [START_REF] Hartmann | Adaptive discontinuous Galerkin finite element methods for nonlinear hyperbolic conservation laws[END_REF] different choices for constructing Sp h are given, including an h-refinement and a p-refinement. As in the work of [START_REF] Hinz | Goal-Oriented Adaptive THB-Spline Schemes for PDE-Based Planar Parameterization[END_REF], the second approach is used in the present paper, with the same mesh as for S p h , but with a higher-order and with the same regularity, i.e. Sp h = S p+1 h .

As specificed in the end of Section 2.4, the DWR method for modal analysis requires A(v, φ) = W (0, v, φ) and B(v, φ) = M(v, φ). For buckling analysis, A(v, φ) = W (u L , v, φ) and B(v, φ) = W (0, v, φ) with the first operator defined about a prebuckling solution u L .

Goal functionals for Isogeometric Kirchhoff-Love Shells

The remainder of this section focusses on defining the goal functional L(u), see Eq. ( 19), together with its variation L (u) such that the dual problem (Eq. ( 21)) can be solved and the error estimate (Eq. ( 23)) can be computed.

In general, the goal functional can be defined in a point, on a boundary or over the domain:

L(•) = Ω l(•, x) dΩ , Domain-wise, (40)
L(•) = ∂Ω l(•, x) dΓ , Boundary-wise, (41)
L(•) = i∈I l(•, x i ), Point-wise, (42)
where Ω denotes the integration domain, ∂Ω a side of Ω and I a set of indices corresponding to points x i ∈ Ω. Furthermore, l(•, x i ) denotes a goal functional summant or integrant, which has a variation denoted by l (•, φx i ). The variation of L, denoted by L (•, φ, x i ), directly follows from l (•, φx i ) due to linearity of integrals and summation. In addition, we classify two different types of goal functional integrants, resulting in norm-based and component-based goal functionals. In the former case, l is of the form l = A 2 with variation l = 2A • A . For component-based goal functionals, we define l = A • e i with variation l = A • e i . Here, e i is a unit vector in direction i. In Table 1 we provide some goal functional integrants or summants l(•, x i ). Together with their variations l (•, φx i ), these provide L (•, φ, x i ) due to linearity of integrals and summation. The tensor-based goal functionals refer to goal functionals that could be used for any second-order tensor, e.g. the membrane strain tensor ε(u) or the flexural moment tensor m(u).

Table 1: Overview of the goal functionals. Here u h = u(x) is the discrete deformation tensor depending on position coordinate x, C h = C(u h ) is the deformation tensor based on u h and C(φ h = C (u h , φ) is its variation, T(A) is the transformation of a second-order tensor A from the undeformed contravariant basis to the basis spanned by the principal directions. Note that the variation [T(A)] of T(A) is T(A ), since the spectral decomposition of the deformation tensor itself is just a linear change of tensor basis.

Displacement-norm l u (u h ) = u h (x) 2 l u (u h , φ h ) = 2u h • φ h Displacement-component lu i (u h ) = u h (x) • e i l u (u h , φ h ) = e i • φ h Stretch-norm l λ (u h ) = T(C h ) 2 l λ (u h , φ h ) = 2T(C h ) • φ h Stretch-component l λ i (u h ) = T(C h ) • e i l λ (u h , φ h ) = e i • T(C h ) Tensor-norm l A (u h ) = A 2 l A (u h , φ h ) = 2A • φ h Tensor-component l A i (u h ) = A • e i l A (u h , φ h ) = e i • A
4 Coarsening and Refinement using THB Splines

This section elaborates on coarsening and refinement of isogeometric meshes using THB-splines. In particular, this section elaborates on the Mark, Refine and Transfer blocks of Fig. 1. Firstly, Section 4.1 will provide a brief background on THB-splines, which enable the Refine step of the adaptive meshing flowchart. Then, Section 4.2 elaborates on methods for suitable grading for refinement meshes; which is required to provide admissible refinement with (Truncated) Hierarchical B-spline ((T)HB) bases, given labelled elements. Section 4.3 elaborates on the labeling method for the Mark step, given an element-wise error distribution, taking admissibility into account. Lastly, Section 4.4 elaborates on the quasi-interpolation method that is used to Transfer the solution of one solution step to the next. The notations in this section will be closely related to those used in [START_REF] Carraturo | Suitably graded THBspline refinement and coarsening: Towards an adaptive isogeometric analysis of additive manufacturing processes[END_REF][START_REF] Bracco | Refinement algorithms for adaptive isogeometric methods with hierarchical splines[END_REF]

(T)HB-Splines

Refinement of B-spline meshes can be done using (Truncated) Hierarchical B-splines ((T)HB-splines), of which the details can be found in [START_REF] Vuong | A hierarchical approach to adaptive local refinement in isogeometric analysis[END_REF][START_REF] Giannelli | THB-splines: The truncated basis for hierarchical splines[END_REF]. The conceptual idea behind (T)HB-splines is that they are constructed from a sequence of N nested tensor B-spline spaces in different levels l = 0, ..., N -1, denoted by

V 0 ⊂ V 1 ⊂, ..., V N -1
with associated bases B of degree p on a grid G with elements Q. The parametric domains are defined as Ω = Ω 0 ⊇ Ω 1 ⊇ ... ⊇ Ω N -1 = ∅. By defining the set of active cells by

G := {Q ∈ G : Q ⊂ Ω ∧ Q ⊂ Ω +1 }, the hierarchical mesh is defined as Q = {Q ∈ G : = 0, ..., N -1}. In
Q 0 1 β i 0 0.5 1 ξ Q 0 0.5 1 ξ 0 0.5 1 ξ
Fig. 2: Principles of refinement for different spline bases. The top figures represent the basis on level B 0 , optionally with refined basis functions coloured blue. Bottom pictures represent refined bases. left) uniform refinement (hence B 1 ; middle) HB-refinement; right) THB-refinement, with truncated basis functions coloured yellow. Note that the refinement basis functions are from V 1 . The unrefined unique knot vector in all cases is Ξ = {0, 1/8, 2/8, . . . , 7/8, 1} and the degree of the basis is 2. The bases are generated in G+Smo [START_REF] Jüttler | Geometry + Simulation Modules: Implementing Isogeometric Analysis[END_REF]. that for THB-splines a truncation is performed to ensure partition of unity, which is discussed in more detail in [START_REF] Giannelli | THB-splines: The truncated basis for hierarchical splines[END_REF].

Admissible Meshing

The concept of admissible meshing was discussed in [START_REF] Buffa | Adaptive isogeometric methods with hierarchical splines: Error estimator and convergence[END_REF][START_REF] Buffa | Mathematical Foundations of Adaptive Isogeometric Analysis[END_REF][START_REF] Bracco | Refinement algorithms for adaptive isogeometric methods with hierarchical splines[END_REF]]. An admissible mesh of class m is a mesh of which the truncated basis functions belong to at most m successive levels and mesh admissibility ensures that the number of basis functions acting on a mesh elements does not depend on the number of levels in the hierarchy, but on the parameter m. In order to guarantee mesh admissibility for refinement and coarsening operations, refinement and coarsening neighborhoods are defined such that admissible meshes can be constructed recursively, which is discussed in more detail in [START_REF] Buffa | Adaptive isogeometric methods with hierarchical splines: Error estimator and convergence[END_REF][START_REF] Buffa | Mathematical Foundations of Adaptive Isogeometric Analysis[END_REF][START_REF] Bracco | Refinement algorithms for adaptive isogeometric methods with hierarchical splines[END_REF]. In Fig. 3, we illustrate a simple mesh together with the refinement neighborhood of some selected elements. The T-

refinement neighborhood N r (Q, Q, m) of element Q is defined as N THB r (Q, Q, m) = Q ∈ G -m+1 : ∃Q ∈ S(Q, -m + 2), Q ⊆ Q , ( 43
) (a) (b) (c) (d) (e) (f)
Fig. 3: Recursive marking strategy on the marked element of level on the initial mesh represented in (a). As a first step, the support extension of the marked element is obtained (b), from which the parents that are active on level -1 define the Tneighborhood of the marked cell (c). Starting the same procedure on the marked cells of level -1, the support extension can again be obtained (d) with their corresponding parents on level -2, marking the T-neighborhood of the marked elements of level -1 (e). The complete recursive marking from the marked element in (a) is depicted in (f).

where S(Q, k) is the multi-level support extension with respect to level k. [START_REF] Carraturo | Suitably graded THBspline refinement and coarsening: Towards an adaptive isogeometric analysis of additive manufacturing processes[END_REF]. When coarsening element Q , the coarsening neighborhood ensures that the newly activated basis functions are not active on the surrounding basis functions of level + m. In other words, if element Q of level is the element to be coarsened, then the coarsening neighborhood defined by

The coarsening neighborhood N c (Q) of element Q ∈ G is defined by
N c (Q, Q, m) := Q ∈ G +m-1 : ∃Q ∈ G and Q ⊂ P (Q), with Q ∈ S(Q , ) , (44)
must be empty. Here, P (Q) denotes the parent of Q, i.e. the unique cell Q ∈ G -1 such that Q ⊂ Q . Note the small difference with respect to the definintion given by [START_REF] Carraturo | Suitably graded THBspline refinement and coarsening: Towards an adaptive isogeometric analysis of additive manufacturing processes[END_REF], since the coarsening neighborhood in their work is defined for the element Q of level which will be activated, i.e. Q is the parent of Q for which the coarsening neighborhood is defined here. Given the definition in Eq. ( 44) and given a set of elements marked for refinement M r , a coarsening neighborhood checking elements marked for refinement, can be defined: The cell marks the cell of level that is marked for coarsening to its parent and the cells mark cells that are marked for refinement. The ring around the cell marked for coarsening depicts the region that should be checked for the coarsening neighborhood. That is, it defines the region that should not contain cells of level + 1 (for N c ) or cells of level that are marked for refinement (for N c ). The cells for which N c = ∅ are marked in (d) and the cells with N c = ∅ are marked in (e). The final mesh after refinement and coarsening is depicted in (f). The coarsened elements that satisfy N c ∪ N c = ∅ are marked as coarsened elements.

N r c (Q, Q, m, M r ) := Q ∈ G +m-2 : Q ∈ M r , ∃Q ∈ G -1 and Q ⊂ P (Q), (a)
with Q ∈ S(Q , -1) . [START_REF] Giannelli | THB-splines: An effective mathematical technology for adaptive refinement in geometric design and isogeometric analysis[END_REF] In other words, this is the coarsening neighborhood that checks whether for element Q of level to be coarsened there are elements in the marked set M r that will be part of the coarsening neighborhoord as soon as they are refined; thus it uses Eq. ( 44) with -1. This neighborhood ensures that coursening labelling can be performed conforming with teh Dörfler marking strategy and without refining first. This avoids to compute element error contributions on an in-between mesh which has been refined first. Obviously, if another element with the same parent as Q is marked for refinement, no coarsening should take place. An element can be coarsened

if N r c (Q, Q, m, M r ) = ∅. Combining both neighborhoods, an element Q of level can be coarsened if and only if Nc (Q, Q, m, M r ) = N c (Q, Q, m) ∪ N r c (Q, Q, m, M r ) = ∅.
In Fig. 4, the coarsening neighborhood is illustrated for a simple mesh.

Labeling Methods

Let Q be the ordered set of Q such that e k ≥ e k+1 ∀k ∈ Q, where e k denotes the error of element k. Then, the Dörfler marking strategy [START_REF] Dörfler | A Convergent Adaptive Algorithm for Poisson's Equation[END_REF] is defined as the elements Q i ∈ Q, i = 0, ..., k such that their sum is smaller than a fraction ρ r of the total element error e = i e i :

M Dörfler r = Q i ∈ Q : i k=0 e k < ρ r e . (46)
This marking strategy, however, does not take into account the contributions of the elements that are marked because they are part of a refinement neighborhood of a marked element Q i ∈ M r . Therefore, we define the index set I K r as the set of element indices that span contain elements Q k ∈ Q and their refinement neighborhoods:

I K r := k ∈ 1, ..., K : N r (Q k , Q, m) ∪ Q k : Q k ∈ Q , ( 47
)
and we define κ r as the maximum index for which the sum of all elements with indices i in I κr r is smaller than the error tolerance ρ r e:

κ r := arg max i∈I K r e i < ρ r e, (48)
such that the Dörfler marking including refinement neighborhoods is

M r = {Q k ∈ Q : k ∈ I κr r }. (49)
For marking a set of coarsening elements, M c the Dörfler marking procedure can be followed again. The original Dörfler marking strategy would be coarsening the elements Q i ∈ Q such that their total element error is smaller than a fraction of the total element error ρ c e, with coarsening parameter ρ c :

M Dörfler c = Q i ∈ Q : N k=i e k < ρ c e . (50)
Similar to marking for refinement, the marking rule for coarsening can be specified more precisely by including admissible coarsening. In this case, the elements for which Nc (Q, Q, m, M r ) = ∅ holds are added in the sum of marked elements. Therefore, let us define the index set J K that contains all elements Q k ∈ Q for which the admissible coarsening condition holds, starting from the element with the smallest error, i.e. Q N .

I K c := k ∈ 1, ..., K : Nc (Q, Q, m, M r ) = ∅, Q N -k-1 ∈ Q . ( 51
)
Similar to κ r , we define κ c as the maximum index for which the sum of all elements with indices i in I κc c is smaller than the error tolerance ρ c e:

κ c := arg max i∈I K r e i < ρ c e, (52)
such that the Dörfler taking into account coarsening admissibility is defined as

M c = {Q k ∈ Q : k ∈ I κc c }. (53)
An alternative to the Dörfler marking strategy is a strategy where a given fraction of the total number of elements is marked. In that case, the formulations from Eqs. ( 49) and ( 53) would still hold, but in Eqs. ( 47) and ( 51) the indices κ r and κ c are defined by the sum of the marked elements in respectively I K r and I K c .

Whether to mark a set for refinement or coarsening, i.e. to construct M r and M c depends on the global error ∆L following from the DWR and user-defined tolerances for refinement and coarsening. Let tol r be the tolerance for refinement and tol c the tolerance for coarsening, such that M r = ∅ if and only if ∆L > tol r and M c = ∅ if and only if ∆L < tol r . As a consequence, if tol r ≥ tol c , refinement and coarsening are never performed simultaneously. If tol r < tol c a band with bandwidth tol c -tol r is defined, in which refinement and coarsening are performed simultaneously. In the present work, tolerances are defined such that the latter condition is satisfied, and the adaptivity iterations are terminated when ∆L ∈ [tol r , tol c ], i.e.:

     M r = ∅, M c = ∅ if ∆L < tol r , M r = ∅, M c = ∅ if ∆L > tol c , M r = ∅, M c = ∅ if tol r ≥ ∆L ≥ tol c , (54)
given tol r ≤ tol c .

Note that the total element error e and the total estimated error of the system of equations ∆L are not necessarily the same, since the element error measure e k can be defined in different ways. In case of the DWR method, a natural choice is to choose e k as the element-wise integrals of ∆L from Eq. ( 23). However, integrating the squared norm of the integrant from Eq. ( 23) would yield strictly positive element errors, making the ordering of the set of element errors simple.

Quasi-Interpolation

In the discrete setting, solution of the problem u h is represented by the THB-spline basis φ i ∈ S p h together with the solution coefficients α i ∈ R. In case of analyses with multiple solution steps (e.g. dynamic or quasi-static analysis), mesh refinements can be performed after each solution step. As a consequence, the solution at load step k + 1 is defined on another set of basis functions { φi } ∈ S p h with corresponding coefficients ᾱk i compared to the previous solution at step k. In order to transfer the coefficients α k i to the new basis, an interpolation scheme needs to be used.

Interpolation on a spline basis can be a costly part of the simulation. Global interpolation implies that the contributions of all basis functions are taken into account in the interpolation. This requires solving a large dense system. An efficient way of interpolating spline coefficients for hierarchical basis is a so-called quasi-interpolation scheme [START_REF] Speleers | Effortless quasi-interpolation in hierarchical spaces[END_REF][START_REF] Giust | Local (t)HB-spline projectors via restricted hierarchical spline fitting[END_REF]. Here, on each level of the hierarchical basis a quasi-interpolant is constructed. This quasi-interpolant interpolates a given function f over the support of each basis function individually, to find the coefficient related to that basis function. More precisely, given a function f ∈ C(Ω 0 ), the quasi-interpolant for level is defined as

Λ (f ) = N i=1 λ i, (f )B i, , = 0, . . . , n -1, (55)
where the coefficients λ i, are suitable linear functionals on C(Ω 0 ). Across all levels = 0, . . . , n -1, the interpolant for the function becomes:

Λ(f ) = n-1 =0 i∈I ,Ωn λ i, (f )B i, , (56)
where B T i, ,Ωn is a THB spline of level constructed on domain Ω n . For any basis function B i, , the coefficient λ i, are found by locally interpolating the function f onto all active basis functions B j, , j ∈ J in the support of B i, . This gives coefficients λ j, , j ∈ J of which coefficient i gives λ i, .

Algorithmic Overview

In Fig. 1 the adaptive isogeometric method for solution stepping problems has been presented. Based on Sections 2 to 4, a summarized workflow for adaptive isogeometric shell analysis is depicted in Fig. 5 and Algorithm 1.

The Solve block involves solving the nonlinear isogeometric Kirchhoff-Love shell equation from Eq. ( 16). This variational formulation involves geometric and material non-linearities and can potentially also involve load non-linearities. After solving the Kirchhoff-Love shell problem, the discrete solution vector u h is passed to the Estimate block. Here, the DWR method is solved by computing the adjoint problem in the primal space (Eq. ( 21)) and in the enriched space (Eq. ( 24)). Then, the element-wise error estimate can be obtained by integrating Eq. ( 23) element-wise. The element-wise errors e k can be passed to the Mark block, where elements are marked for refinement (Eq. ( 53)) if the total error ∆L is larger than a lower (refinement) tolerance tol r and a coarsening marking (Eq. ( 53)) is performed if the total error is above an upper (coarsening) tolerance tol c . This implies that if tol c < ∆L < tol r , a combined coarsening and refinement step is performed, as described in Eq. [START_REF] Speleers | Effortless quasi-interpolation in hierarchical spaces[END_REF]. In this case, the coarsening marking from Eq. ( 53) is performed given M r . Given the elements marked for refinement and coarsening, collected in M r and M c respectively, the mesh can be Adapted. In order to start the solution interval again, the start point should be Transfer red to the new mesh and the governing equation can be solved again if the error is not in the interval [tol r , tol c ] or if the number of refinement iterations i exceeds the maximum number of refinement iterations, I max . If the total error is in the interval [tol r , tol c ] or if the number of refinement iterations i exceeds the maximum number of refinement Solve Eq. 16

Estimate

Eqs. [START_REF] Coradello | Coupling of non-conforming trimmed isogeometric Kirchhoff-Love shells via a projected super-penalty approach[END_REF] The adaptive meshing iterations are performed within each solution step until the total error ∆L is contained in the interval [tol r , tol c ], tol r < tol c , following the tolerances in Section 4.3. In case of convergence, the solution is advanced, e.g. with an arc-length iteration. Algorithm 1 provides an algorithm corresponding to this flow-chart.

iterations, I max , the solution can be advanced, e.g. using a load-stepping or an arclength method. Thereafter, the governing equations can be Solved again. Note that if I max = 1, no inner iterations for adaptive meshing are performed.

Numerical examples

In this section, several numerical examples are presented. The examples represent different applications of the theory presented in this paper and -without loss of generality -all employ Isogeometric Kirchhoff-Love shells. The first three examples illustrate the performance of the DWR error estimator and the last three example demonstrate the use of this error estimator for adaptive meshing. More precisely, the numerical examples performed in this section, as well as their purpose, are:

Linear static analysis of a square plate (Section 6.1) A simple example of linear Kirchhoff-Love shell theory is presented. In this case, error estimators using the DWR are computed for different goal functionals and verified using the actual error computed from manufactured solutions. The goal of this benchmark problem is to evaluate the accuracy of the error estimators in linear static analysis. Modal analysis of a circular plate (Section 6.2) Since the analytical eigenvalues and eigenmodes are known for this case, the goal of this benchmark problem is to verify the error estimator for a vibration eigenvalue problem, given in Eq. ( 27) with the stiffness and mass operators A(v, φ) = W (0, v, φ) and B(v, φ) = M(v, φ).

Algorithm 1 An algorithmic summary of the goal-adaptive meshing routine employed in the present work. Figure 5 provides a graphical summary of this algorithm.

1: for loadsteps do 2:

while ∆L ∈ [tol l , tol u ] and i < I max do 3:

Compute the primal solution u h (Eq. 16)

4:

x y z W L f z u = 0 E = 10 6 [Pa] ν = 0.3 [-] L = 1 [m] W = 1 [m] t = 0.01 [m]
Fig. 6: Geometry and parameters for the example of a unit-square plate with a distributed vertical load f z given by Eq. ( 57). The displacements are fixed on all edges.

In the following subsections, the short-hand notations L an = L(u an ), L num = L(u num ), ∆L an = L an -L num and ∆L num = R(u num , ξh -ξ h ) (see Eq. ( 23)) are used, given the analytical and numerical solutions u an and u num , respectively. Furthermore, where relevant, the parameters ρ r , ρ c , tol r and tol c (see Eqs. ( 49), ( 53) and ( 54)) are fixed per example. A study on finding optimal values for these parameters is out of the scope of this paper. Lastly, all simulations are performed using the open-source Geometry+Simulation modules [START_REF] Jüttler | Geometry + Simulation Modules: Implementing Isogeometric Analysis[END_REF].

Linear static analysis of a square plate

For the linear shell example, let us consider a flat plate with unit dimensions L = W = 1 [m], a thickness of t = 10 -2 [m] and with material parameters E = 10 6 [Pa], ν = 0.3, which is clamped on all sides, see Fig. 6. A load vector of

f = 2EAt 3 1 -ν 2 x 4 -2x 3 + 3x 2 -2x + y 4 -2y 3 + 3y 2 -2y + 12x 2 y 2 -12x 2 y -12xy 2 + 12xy + 1 3 e z , (57)
is applied, based on the manufactured solution given by

u an = Ax 2 (x -1) 2 y 2 (y -1) 2 e z . (58)
Using this manufactured solution, any goal functional L an can be evaluated. Solving the primal problem for this linear shell example gives u num , which can be used to compute the DWR error estimate of ∆L num .

In Fig. 7, the results for the linear shell problem are given. The title of each column represents the goal-function that is used for error estimation in this column. The top row provides the errors ∆L an and ∆L num with respect to an uniformly refined mesh size. As can be seen in this figure, ∆L num quickly converges to ∆L an for different splines orders p. In addition, the bottom row of Fig. 7 provides the efficiency of the error estimator, given by ∆L num /∆L an . These figures confirm convergence of the DWR estimates to the analytical goal functional errors for all considered goal Fig. 7: Linear static analysis of a clamped plate with a uniformly distributed load according to Eq. ( 57) according to a manufactured solution from Eq. ( 58). The top row provides ∆L against the uniform mesh size h. The bottom row presents the efficiency of the error estimators against the mesh size h. The markers represent error estimates computed via the DWR and the lines represent the exact error, i.e. the error of the numerical solution with respect to the analytical solution. The results are given for goal functionals (from left to right) L = Ω u dΩ (displacement norm), L = Ω λ 2 dΩ (second principal stretch), L = Ω ε(u) dΩ (membrane strain norm) and L = Ω n(u) • e 1 dΩ (first component of membrane force).

functionals. Only for the membrane strain norm goal functional the error estimate for coarse meshes is inaccurate. This can possibly be explained by the in-plane shear strain that cancels out over the whole domain but which does contribute in the norm ε

Concluding, the linear shell benchmarks shows that for different goal functionals the DWR method provides accurate estimation of the error ∆L of the goal functional L starting at relatively small mesh sizes of h < 10 -1 .

Modal analysis of a circular plate

z R u = ∇u • n = 0 E = 10 6 [Pa] ν = 0.3 [-] ρ = 1 [kg/m 3 ] R = 0.5 [m] t = 0.01 [m]
Fig. 8: Geometry and parameters for a vibrating circular plate with a clamped boundary.

solutions for the eigenfrequencies of the circular plate are obtained by

ω n = γ 2 n D/ρt, ( 59
)
where γ n is the n th root of the equation (I m-1 (γR) -m/RI m (γR))J m (γR) -I m (γR)(J m-1 (γR) -m/RJ m (γR)) = 0 following from a separation of variables solution [START_REF] Verhelst | Equilibrium Path Analysis Including Bifurcations with an Arc-Length Method Avoiding A Priori Perturbations[END_REF], R is the radius of the plate and D = Et 3 /(12(1 -ν 2 )) = 9.16 • 10 -8 is the flexural rigidity. As stated in Section 3.2, the goal functional eigenvalue problems is given in Eq. ( 31), and requires the eigenvalue problem in Eq. ( 26) to be solved with linear operators A(v, φ) = W (0, v, φ) and B(v, φ) = M(v, φ), see Eqs. ( 10) and ( 14).

Figure 9 presents the first four eigenmodes in the top row. Furthermore, ∆L an and ∆L num as a function of the element size for uniformly refined domain are given in the middle row for the first four eigenmodes. These plots show that the approximation for the error ∆L num approximates ∆L an . In the bottom row of Fig. 9 the efficiencies also show that the approximation converges to an efficiency equal to 1. However, for the p = 4 line, the efficiency degrades when the 'exact' error obtained by the analytical solution approaches values around 10 -11 , which is attributed to the approximation of the roots γ i and the precision of the eigenvalue solver.

Concluding, the modal analysis benchmark shows that the DWR method provides accurate estimation of the eigenfrequency error for different considered mode shapes.

Linear buckling analysis of a square plate

Similar to modal analysis, DWR error estimation for buckling analysis also relies on the formulations in Section 3.2. The difference with the modal analysis error estimation is that the buckling analysis error estimation involves the solution of a pre-buckling solution and no mass matrix. As an example for buckling analysis, a square simply supported plate is considered, see Fig. 10,with L and with equal loads is given in [START_REF] Jones | Buckling of Bars, Plates, and Shells[END_REF] and reads:

σ m,n c t = Dπ 2 L 2 m 2 + n 2 , (60)
with D = Et 3 /12(1 -ν 2 ) the flexural rigidity of the plate. Using this expression, the first four unique modes are, indexed in ascending order:

σ 1,1 c t = σ 1 c t = 1.808 [N/m], σ 2,1 c t = σ 1,2 c t = σ 2 c t = 4.519 [N/m], σ 2,2 c t = σ 3 c t = 7.230 [N/m], σ 3,1 c t = σ 1,3 c t = σ 4 c t = 9.038 [N/m].
Figure 11 depicts the analytical error ∆L an and the DWR error estimate ∆L num as a function of the mesh size h for uniform refinements. Both errors are normalized with the analytical value of the critical buckling load. As can be seen in this figure,

u z = 0 u z = 0 u x = u z = 0 E = 10 6 [Pa] ν = 0.3 [-] L = 1 [m] W = 1 [m] t = 0.01 [m]
Fig. 10: Geometry and parameters for the plate buckling example. A distributed load of σt is acting on two boundaries in two different directions, and the other boundaries are simply supported and fixed in out-of-plane direction.

the DWR prediction of the error converges with a rate of convergence of 2(p -1) for for all degrees p until it reaches values of around 10 -10 after which the errors stagnate and increase again (in particular for p = 4). This behaviour can be is similar to the behaviour observed in [START_REF] Liu | Balancing truncation and round-off errors in FEM: One-dimensional analysis[END_REF] and can be attributed to the round-off errors as discussed. These errors occur when the number of degrees of freedom is large enough and the machine precision is limited. In case of a buckling simulation, where the non-linear stiffness operator is constructed on an initial solution of a linear simulation, the influence of round-off errors is expected to occur sooner. In addition, it can also be seen that the error computed using the analytical solution stagnates. This is due to the fact that the numerical approximation of the critical buckling load shows small variations depending on the solution to the linear problem that is solved to obtain the tangential stiffness matrix to compute the generalized eigenvalue problem for buckling. For the results presented in Fig. 11, the load σt = 10 -4 [N/m] was used to compute the buckling linearization.

Non-linear analysis of a pinched thin plate

In a next example, we consider a square membrane subject to a point-load in the middle and with corners fixed in all directions, see Fig.

c t = σ 1 c t = 1.808 [N/m], σ 2,1 c t = σ 1,2 c t = σ 2 c t = 4.519 [N/m], σ 2,2 c t = σ 3 c t = 7.230 [N/m], σ 3,1 c t = σ 1,3 c t = σ 4 c t = 9.038 [N/m].
case are based on displacements as well as on principal stresses:

L(u) = Ω u dΩ , L(u) = Ω σ p • e y dΩ . (61)
Figure 13 presents the deformed membrane for the last step of the adaptive simulation. Furthermore, Fig. 14 presents the estimated error ∆L given the goal functionals in Eq. ( 61) for the uniform refinement as well as for the adaptive refinement simulation with maximum level 8 or 11. Moreover, Figs. 15 and 16 provide the absolute element-wise errors for the adaptive refinement simulation and for the

x y z W L u = 0 u = 0 u = 0 u = 0 P s y m m sy m m E = 1.0 [MPa] ν = 0.3 [-] P = 4 • 10 -7 [N] L = 1 [mm] W = 1 [mm] t = 10 -3 [mm]
Fig. 12: Geometry and parameters for a square thin plate subject to a point load P in the middle. The plate is fully constrained in every corner. Because the problem is symmetric, only a quarter of the domain is modelled. Hence, symmetry conditions are applied. On the x-aligned symmetry axis, this implies that u y = ∂uz ∂y = ∂ux ∂y = 0 and on the y-aligned symmetry axis this implies that u x = ∂uz ∂x = ∂uy ∂x = 0.

Fig. 13: Deformed surface from the benchmark presented in Fig. 12. The result is the last solution from the adaptive meshing routine with deformation norm goal-functional of which the results are presented in Fig. 14.

uniform refinement series for both considered goal functionals. The contour lines in these error fields represent the vertical deflection of the sheet.

From all provided results, it can be observed that the adaptive mesh provides for both goal functionals an efficient converging mesh, where the accuracy per degree of freedom is higher compared to uniformly refined meshes. However, it can also be seen that the total error ∆L is not strictly decreasing for both goal functionals. The bottom figures of Fig. 14 indicate that this point is closely related to the maximum depth that is reached: the percentage of the total element error e that can still be refined rapidly decreases, meaning that the only elements that are still available for refinement are the ones that have insignificant contribution to the total element error e, deeming refinement of these element meaningless. It can be seen from Figs. 15 and 16 that the maximum element depth is reached in the corner where the sheet is fixed. After the maximum level is reached, the refined elements start distributing over the diagonal of the domain and in the top-left corner where the load is applied. For both goal functionals, Fig. 14 shows that some further decrease in the total error ∆L can be gained after the maximum error is reached, but that most effective is to increase the maximum refinement depth. Comparing the error fields and meshes for 10 -7

#DoFs

Adaptive [START_REF] Hu | Isogeometric analysis of thin Reissner-Mindlin shells: locking phenomena and B-bar method[END_REF] Adaptive ( 11) Uniform ..., blocked ..., blocked Fig. 14: Estimated error convergence (top) and the percentage of the total element error e that is available for refinement (bottom) against the number of degrees of freedom (DoFs) for adaptively and uniformly refined meshes with respect to the goal functionals from Eq. ( 61). The markers labeled with a black border are the markers for which the mesh is plotted in Figs. 15 and16. The filled markers represent points where refinement is not blocked, and empty markers points where refinement is blocked because the maximum refinement depth is reached. Note that the errors are computed before refinement, hence blocked elements in iteration i have effect on the error computation in iteration i + 1.

both goal functionals (see Figs. 15 and16), it can be seen that the second-principal stress-based error field shows narrower error bands in the finest depicted uniform refinement error fields than the displacement-based error estimator. As a consequence, the corresponding adaptive meshes show that the elements indeed tend to be broader distributed along the diagonal of the domain in case of the displacement-driven refinement (Fig. 15).

Concluding, the non-linear shell benchmark shows that an adaptive meshing strategy provides accurate solutions on meshes with a small number of degrees of freedom compared to uniform meshes. Furthermore, the benchmarks shows the importance of refinement levels, meaning that convergence of the adaptive meshing strategy vanishes as soon as the elements contribution to a large extent to the total error are on the lowest allowed level. This stresses the relevance of spline constructions that allow for deep levels of refinement for with computational costs.

Snap-through instability of a cylindrical roof

In order to present results of the presented adaptive isogeometric method for quasistatic problems, hence completing full cycles in Fig. 1, the well-known benchmark of a collapsing cylindrical roof [START_REF] Sze | Popular benchmark problems for geometric nonlinear analysis of shells[END_REF] subject to a point-load is considered. The goal of this benchmark problem is to evaluate whether the presented adaptive isogeometric method provides DoF-wise efficient solutions compared to solutions with uniform refinements.

The geometry for the benchmark problem is presented in Fig. 17. Here, the radius of the roof R = 2.540[m], the angle is θ = 0.1 [rad] and the length and thickness are respectively L = 0.508 [m] and 6.35 • 10 3 [mm]. Moreover, the material properties are E = 3102 [MPa] and ν = 0.3 [-] for a Saint-Venant Kirchhoff material. Only a quarter of the roof is modelled, as depicted in Fig. 17, because of symmetry. The simulation is performed using a Crisfield arc-length method with arc-length ∆L = 25 and a zero force-scaling. The goal functional that is used for error evaluation and adaptivityis based on the norm of the flexural strain tensor over the whole domain L = Ω m(u) dΩ. The jump parameter for admissible meshing is set to m = 2. The mesh will be refined when ∆L > tol r and coarsened when ∆L < tol c . As discussed in Section 5, tol r < tol c such that the mesh is refined and coarsened simultaneously when ∆L ∈ [tol r , tol c ], which is also the condition for termination. The maximum number of mesh adaptivity iterations is set to 5 in this case. The tolerances tol r and tol c are determined based on the results of uniformly refined simulations with 16 × 16 (918 DoFs) and 32 × 32 elements (3366 DoFs) by taking a wide band around the error envelopes in Section 6.5 excluding peaks. The tolerances are (tol r , tol c ) = (10 -10 , 10 -8 ). It should be noted that these tolerances can also be based on requirements in engineering, or they can be determined during the computations; both are beyond the scope of this paper. The refinement parameter is set to ρ r = 0.5 and the coarsening parameter ρ c = 0.05 and the maximum refinement level is 11. The adaptive and uniform meshes are modelled with cubic B-spline p = 3 elements.

In Fig. 19, the results of a simulation of the collapsing roof for the adaptive mesh are plotted in a u , λP , w A -space. Reference solutions for the present benchmark problem are typically given in the λP , w A -space, but since the solution curve is not bijective an alternative coordinate u is used to represent solutions for this benchmark problem. This is motivated by the projection of the solutions λP and w A projected against u in Fig. 19. The results of [START_REF] Sze | Popular benchmark problems for geometric nonlinear analysis of shells[END_REF] are provided as a reference.

In Fig. 18, the error and the number of degrees of freedom are plotted against u . The error envelopes for the uniform meshes show that a large peak occurs around u =, relating to the first limit point of the collapse of the roof, as seen by the markers in Fig. 19. Additionally, it can be seen that the present algorithm providing mesh adaptivity manages to keep the error within specified bounds (see Fig. 18, top), except on the peak just before u where the maximum number of adaptivity iterations is insufficient. Furthermore, it has consistently less degrees of freedom than the uniform mesh with 32 × 32 elements. In Fig. 20 a selection of meshes is provided.

L x y z R θ u z = 0 u z = 0 symm A λP E = 3102 [MPa] ν = 0.3 [-] P = 3 • 10 3 [N] L = 0.508 [m] R = 2.540 [m] θ = 0.1 [rad] t = 6.35 • 10 -3 [m]
Fig. 17: Geometry and parameters for cylindrical roof with a point-load P in the middle. The bottom-right corner of each domain corresponds to the point A: the reference point for which the z-displacements are plotted. The roof is free on the curved edges and simply supported (u z = 0) on the straight edges. As a consequence, the problem is symmetric and a quarter of the domain is modelled. On the x-aligned symmetry axis, this implies that u y = ∂uz ∂y = ∂ux ∂y = 0 and on the y-aligned symmetry axis this implies that u x = ∂uz ∂x = ∂uy ∂x = 0.

The meshes are provided as series of 4 consecutive meshes around the limit points of the solution curve, as indicated in Fig. 19, and are depicted in increasing order from left to right for the first (top) and second (bottom) limit points. The black-bordered markers in Figs. [START_REF] Leonetti | A robust penalty coupling of non-matching isogeometric Kirchhoff-Love shell patches in large deformations[END_REF] and 19 indicate the points of which the meshes are shown. From the first row of meshes in Fig. 20, it can be seen that the first limit point requires relatively fine meshes and that the elements start concentrated around point A and its diagonal opposite and spread out on the bottom symmetry boundary as the snapping takes place. Furthermore, in the bottom row of Fig. 20, it can be seen that the second limit point does not require many elements, hence the number of elements is slowly decreasing throughout this section of the load-displacement curve. For a complete overview of the mesh in each load-step, we refer to Video 1 in the supplementary material of this paper.

Concluding, this example shows that the goal-adaptive meshing procedure is capable of keeping the error in terms of a goal functions within pre-defined bounds for a solution stepping simulation with limit-point instabilities. Throughout the simulation, the procedure keeps a relatively high efficiency per degree of freedom compared to uniform meshes. It should be noted, however, that the adaptive refinement iterations require a higher computational demand. Therefore, the next example provides a procedure where no adaptivity iterations are performed.

Wrinkling analysis

As a last example, the wrinkling analysis of a thin membrane subject to a tensile load is considered. This problem is inspired by [START_REF] Cerda | Wrinkling of an elastic sheet under tension[END_REF][START_REF] Cerda | Geometry and Physics of Wrinkling[END_REF][START_REF] Panaitescu | Birth and decay of tensional wrinkles in hyperelastic sheets[END_REF] and was previously modelled using isogeometric Kirchhoff-Love shells in [START_REF] Verhelst | Stretch-Based Hyperelastic Material Formulations for Isogeometric Kirchhoff-Love Shells with Application to Wrinkling[END_REF], to which we refer for a detailed problem set-up. The goal of this benchmark in the present paper is to demonstrate the use of the goal-adaptive meshing procedure on a bifurcation problem with geometric and material non-linearities.

Given the geometry in Fig. 21, a quarter of the domain is considered with a symmetry boundary condition on Γ 4 and an antisymmetry condition on Γ 3 . Furthermore, the boundary Γ 1 is free and on Γ 2 the x-displacement is constant and a horizontal load is applied on this side. The sheet has dimensions L = 280 [mm], W = 140 [mm] and t = 0.14 [mm] such that L/W = 2 and t/W = 10 3 . Furthermore, the material is modelled using a Mooney-Rivlin material model with strain energy density function

Ψ(C) = c 1 2 (I 1 -3) + c 2 2 (I 2 -3) (62)
and with parameters c 1 = 3.16 • 10 5 [Pa] and c 1 = 1.24 • 10 5 [Pa]. A reference solution is given by [START_REF] Verhelst | Stretch-Based Hyperelastic Material Formulations for Isogeometric Kirchhoff-Love Shells with Application to Wrinkling[END_REF] for isogeometric Kirchhoff-Love shell analysis and using ANSYS and LS-DYNA FEA models. Furthermore, [START_REF] Panaitescu | Birth and decay of tensional wrinkles in hyperelastic sheets[END_REF] provides experimental data of the maximum amplitude with respect to the strain of the sheet. In the present paper, the reference simulations are performed on uniform cubic meshes with 32 × 32, 64 × 64 elements. The top row of meshes corresponds to the first limit point in Fig. 19 and the bottom row corresponds to the second limit points. The meshes are ordered from left to right for increasing solution steps. The elements are coloured according to the squared error e k normalized by the total error ∆L, i.e. e k /∆L 2 .

( For the adaptive simulation, a THB spline mesh with initially 32 × 32 elements is used and mesh adaptivity is activated after wrinkling initiation since the errors in the pre-wrinkling regime are small because of the lack of out-of-plane deformations of the sheet. The goal-functional is a displacement-based functional on the z-component, i.e. L(u) = Ω u • e z dΩ. The tolerances for refinement and coarsening are determined based on the results of the uniform mesh and read tol r = 10 -14 and tol c = 10 -10 . These parameters are chosen based on the error envelope of the uniform refinement. The adaptive meshing parameters are chosen as (ρ r , ρ c ) = (0.5, 0.005). These parameters are chosen based on the behaviour of the global error in the first load-steps afer bifurcation. Contrary to the previous example in Section 6.5, there are no refinement iterations performed within the load step. This means that the refinement and coarsening operations are performed after the load step based on the magnitude of the error compared to the tolerances. Furthermore, when the error is below a tolerance ρ c,min , the initial mesh is used again. This is done to prevent further coarsening after restabilisation of the wrinkles, i.e. the moment when the wrinkles have dissapeared. For the refinement algorithm, a maximum depth of the THB grid is fixed to 11 levels and the jump parameter is set to m = 2. The wrinkling simulation is performed using a Crisfield arc-length method [START_REF] Crisfield | A Fast Incremental/Iterative Solution Procedure That Handles[END_REF] with a quadratic procedure to compute the mode shape at the bifurcation [START_REF] Wriggers | A quadratically convergent procedure for the calculation of stability points in finite element analysis[END_REF]. This procedure further described in [START_REF] Verhelst | Equilibrium Path Analysis Including Bifurcations with an Arc-Length Method Avoiding A Priori Perturbations[END_REF].

In Fig. 22, the results for an adaptive isogeometric wrinkling simulation are provided. The top figure provides the normalized wrinkling amplitude with respect to the strain of the sheet ( ), compared to the IGA and ANSYS SHELL181 element reference results from [START_REF] Verhelst | Stretch-Based Hyperelastic Material Formulations for Isogeometric Kirchhoff-Love Shells with Application to Wrinkling[END_REF] as well as the experimental results from [START_REF] Panaitescu | Birth and decay of tensional wrinkles in hyperelastic sheets[END_REF]. The figure in the middle provides the error estimate in terms of the goal functional L with respect to the strain of the sheet and the bottom figure provides the number of degrees of freedom with respect to the strain of the sheet. Firstly, the results from both uniform meshes show that the error estimate ∆L is close to zero when the wrinkles initiate, since the sheet is perfectly flat. As soon as the wrinkles form, the error estimate becomes non-zero and it peaks at the moment of restabilisation (i.e. the moment when the amplitude vanishes). After restabilisation, the error estimate is low, but slightly higher than before wrinkling, probably because the sheet is not numerically flat.

The adaptive meshing simulations show that even with zero inner-iterations for mesh adaptivity the adaptive mesh provides accurate results for a relatively small number of degrees of freedom. Just after wrinkling initialisation, the adaptive meshing error peaks due to the coarsening of a large number of elements (as can be seen by the drop of degrees of freedom in the bottom figure). However, the mesh adaptively refines until in the gray region in the middle figure, where combined refinement and coarsening imply that the error balances around tol c , i.e. the upper bound of the marked region. Towards the end of the wrinkling phase, the error increases for the uniformly refined mesh, explaining the increase in number of degrees of freedom for the adaptive mesh. Nevertheless, the adaptive mesh provides more accurate results compared to the 64 × 64 uniform mesh with less degrees of freedom.

In Fig. 23 a selection of meshes from both adaptive simulations are provided, specifically for the wrinkling initiation and restabilisation points. For a complete overview of the mesh in each load-step, we refer to Video 2 in the supplementary material of this paper. The evolution of the meshes shows that the mesh elements concentrate around the wrinkles (bottom-right) and in the top-left corner, which represents the corner between the clamped edge (top) and the free edge. The former is expected since the employed goal functional is based on out-of-plane deformations. However, the fact that the mesh concentrates around the top-left corner is non-intuitive, but tells that this area is important to reduce the global error in terms of the out-of-plane deformations. Furthermore, it can be seen that around the restabilisation of the wrinkles the mesh concentrates around the bottom-left corner, indicating that this corner is of importance in accurately modelling the wrinkling amplitudes in the whole domain.

Concluding, the wrinkling benchmark shows the potential of mesh adaptivity for such applications. With fewer degrees of freedom, the THB-spline mesh is able to approximate the solution around a pre-defined error, eventhough the selection of the refinement and coarsening parameters and the tolerances has not been optimized in this study.

Conclusions

This paper presents an adaptive method for isogeometric Kirchhoff-Love shells. The main contributions of this paper are a goal-adaptive error estimator for isogeometric Kirchhoff-Love shells using the Dual-Weighted Residual Method, a slightly modified suitably graded refinement scheme taking into account refined elements in the definition of the coarsening neighborhood.

Using the Dual-Weighted Residual method and given a pre-defined goal functional (e.g the second-principal stress integrated over the domain), an estimator for the error in terms of this goal functional can be defined. The adjoint problem that needs

# DoFs

Fig. 22: The non-dimensional maximal out-of-plane deformation max(u z )/t (top), the goal functional error ∆L (mid) and the number of degrees of freedom of the computational mesh (bottom) with respect to the strain of the sheet . The markers with a black border represent the points of which the meshes are provided in Fig. 23.

The colored lines are the solutions obtained by the present model with a Mooney-Rivlin material model with uniform or adaptive meshes. The solid line in the top figure is a SHELL181 result obtained using ANSYS and the markers represent experimental data obtained by [START_REF] Panaitescu | Birth and decay of tensional wrinkles in hyperelastic sheets[END_REF].

to be solved on the original mesh and on a nested degree-elevated ('enriched') mesh has been defined for the isogeometric Kirchhoff-Love shell. In addition, the operators for modal and linear buckling analysis have been derived, implying an additional generalized eigenvalue problem to be solved on the enriched mesh. For suitable grading, the works of [START_REF] Buffa | Adaptive isogeometric methods with hierarchical splines: Error estimator and convergence[END_REF][START_REF] Buffa | Mathematical Foundations of Adaptive Isogeometric Analysis[END_REF][START_REF] Carraturo | Suitably graded THBspline refinement and coarsening: Towards an adaptive isogeometric analysis of additive manufacturing processes[END_REF][START_REF] Bracco | Refinement algorithms for adaptive isogeometric methods with hierarchical splines[END_REF] have been closely followed. In order to be able to refine and coarsen in the same iteration, refined elements have been added to the original definition of the coarsening neighborhood.

To assess the proposed adaptive isogeometric method for Kirchhoff-Love shells, few numerical benchmark problems have been evaluated. Linear static analysis provided an analytical solution has been used to evaluate the DWR error estimators. The eigenvalue problems for modal and buckling analyses have been evaluated on respectively the problem of circular plate vibration and square plate buckling. Based on the linear, modal and buckling analysis with analytical solutions, it can be concluded that the DWR estimator for Kirchhoff-Love shells can be used with several goal and in several applications, as it provides high accuracy with respect to the exact errors.

Using the problem of a pinched membrane, the error estimator has been used to adaptively refine a mesh in a non-linear example. From this example, it can be concluded that eventhough the challenging non-linear problem, high accuracy per DoF can be obtained compared to uniformly refined meshes based on different goal functionals. Lastly, the adaptive isogeometric method from the present paper has been evaluated in solution-stepping problems for structural instabilities. Firstly, the method was applied to limit-point instability problem of a collapse of a cylindrical roof. Here, inner adaptivity iterations were performed for each load-step until the error was located in a desired interval. Again, the present method provides a high accuracy per degree of freedom. In addition, it was shown that the method is indeed able to provide adaptive meshes with respect to a pre-defined interval for a given goal functional. The last benchmark problem involves a tension-wrinkling bifurcation instability of a thin membrane. In this benchmark problem, adaptive meshing has been applied in the post-buckling regime based on wrinkling amplitudes. Here, no adaptivity iterations within the load steps were performed, showing that the method is still able to provide good adaptivity with respect to the pre-defined tolerances. Also, the results of the wrinkling error estimators show that the error peaks in the restabilisation phase where the highest deviation with experimental results are observed.

Future developments for the present method include the application on multipatch domains, both coupled with penalty methods as well as with globally continuous bases as presented in [START_REF] Farahat | Isogeometric analysis for multipatch structured Kirchhoff-Love shells[END_REF] to handle more complex geometries. Furthermore, structural dynamics have been left out of scope of this paper, since the DWR for dynamic problems requires backwards-in-time evaluation of the adjoint problem, which is ideally combined with parallel-in-time methods like parareal or MGRIT [START_REF] Falgout | Parallel time integration with multigrid[END_REF]. Lastly, future work can be done on the (adaptive) determination of the adaptive meshing parameters. On the one hand, one can apply the present method on real-world engineering applications, taking realistic goal functionals and margins into the adaptivity algorithm. On the other hand, advanced schemes for triggering pure or combined refinement or coarsening together with their parameters can be futher investigated. Supplementary information. Video 1: Mesh evolution of for the collapsing roof. Video 2: mesh evolution for the wrinkling example.

Fig. 2 ,

 2 an illustration is given for a refined Bspline basis (left), a refined HB-spline basis (middle) and a refined THB-spline basis (right). For the (T)HB-spline basis, this picture depicts the refinement of a single basis function, corresponding to the elements in its support. The (T)HB-spline bases show

Fig. 4 :

 4 Fig.4: Given the mesh from Fig.3-(f), the coarsening neighborhoods are evaluated in (a)-(c) of this figure. The cell marks the cell of level that is marked for coarsening to its parent and the cells mark cells that are marked for refinement. The ring around the cell marked for coarsening depicts the region that should be checked for the coarsening neighborhood. That is, it defines the region that should not contain cells of level + 1 (for N c ) or cells of level that are marked for refinement (for N c ). The cells for which N c = ∅ are marked in (d) and the cells with N c = ∅ are marked in (e). The final mesh after refinement and coarsening is depicted in (f). The coarsened elements that satisfy N c ∪ N c = ∅ are marked as coarsened elements.

Fig. 5 :

 5 Fig. 5: A graphical summary of the adaptive meshing flowchart from Fig. 1 used in the present work. The equations which are used in each step are indicated in the blocks.The adaptive meshing iterations are performed within each solution step until the total error ∆L is contained in the interval [tol r , tol c ], tol r < tol c , following the tolerances in Section 4.3. In case of convergence, the solution is advanced, e.g. with an arc-length iteration. Algorithm 1 provides an algorithm corresponding to this flow-chart.

  As a next example, the vibration modes of a circular plate with clamped boundary conditions are computed. The geometry with boundary conditions is illustrated in Fig. 8. The circular plate has unit diameter, a thickness t = 10 -2 [m], Young's modulus E = 10 6 [Pa], Poisson's ratio ν = 0.3 and density ρ = 1 [kg/m 3 ]. The analytical Z x y

10 - 2 Fig. 9 :

 1029 Fig. 9: Modal analysis of a circular plate. The top row provides the mode shapes from mode 1 up to 4. The mid row provides ∆L against the uniform mesh size h and the bottom row provides the efficiency of the error estimators against the mesh size h. The markers represent error estimates computed via the DWR and the lines represent the exact error, i.e. the error of the numerical solution with respect to the analytical solution. The eigenfrequencies for mode 1 up to 4 are, respectively, ω 1 = 9.56•10 -4 [Hz], ω 2 = 4.14 • 10 -3 [Hz], ω 3 = 1.11 • 10 -2 [Hz] and ω 4 = 1.45 • 10 -2 [Hz].

  u z = 0

12 . 10 - 2 Fig. 11 :

 1210211 Fig. 11: Buckling analysis of a square plate with simply supported boundary conditions. The top row provides the mode shapes from mode 1 up to 4. The bottom row provides ∆L against the uniform mesh size h. The markers represent error estimates computed via the DWR and the lines represent the exact error, i.e. the error of the numerical solution with respect to the analytical solution. The critical loads for unique modes 1 up to 4 are, respectively, σ 1,1 c t = σ 1 c t = 1.808 [N/m], σ 2,1 c t = σ 1,2 c t = σ 2 c t = 4.519 [N/m], σ 2,2 c t = σ 3 c t = 7.230 [N/m], σ 3,1 c t = σ 1,3 c t = σ 4 c t = 9.038 [N/m].

5 L

 5 = Ω σ p • e y dΩ 10 1 10 2 10 3 10 4 10 5

  e k /∆L (a) Uniform refinement e k /∆L 2 (b) Adaptive refinement with maximum level 8. e k /∆L 2 (c) Adaptive refinement with maximum level 11.

Fig. 15 :

 15 Fig. 15: Normalized element error values e k /∆L for uniformly (a) and adaptively (b, c) refined meshes using goal-function L(u) = Ω u dΩ. The meshing steps increase from left to right. The contour lines represent the displacement of the membrane, with intervals of 0.1[mm]. The bottom right corner of the pictures indicates the fixed corner and the top left corner is the corner where the load is applied.

  e k /∆L (a) Uniform refinement. e k /∆L 2 (b) Adaptive refinement with maximum level 8. e k /∆L 2 (c) Adaptive refinement with maximum level 11.

Fig. 16 :

 16 Fig. 16: Normalized element error values e k /∆L 2 for uniformly (a) and adaptively (b, c) refined meshes using goal-function L(u) = Ω σ p •e y dΩ. The meshing steps increase from left to right. The contour lines represent the displacement of the membrane, with intervals of 0.1[mm]. The bottom right corner of the pictures indicates the fixed corner and the top left corner is the corner where the load is applied.

Fig. 18 :

 18 Fig.18: Errors (top) and number of degrees of freedom (DoFs, bottom) for the goal functional L = Ω m(u) dΩ with tolerances tol r = 10 -10 and tol c = 10 -8 for the collapsing roof subject to the displacement norm u . The gray region represents the region ∆L ∈ [tol r , tol c ] where refinement and coarsening are performed and where the adaptivity iterations are terminated. Above the gray region only refinement is performed and below the gray region only coarsening is performed.

Fig. 19 : 2 Fig. 20 :

 19220 Fig.19: Projection the adaptively refined result from the commonly used λP, w Aspace[START_REF] Sze | Popular benchmark problems for geometric nonlinear analysis of shells[END_REF] onto the displacement norm u . The solid lines correspond to the results obtained by the adaptively refined mesh, the triangular markers correspond to points of which the mesh is provided in Fig.20and the cross-markers indicate the result from[START_REF] Sze | Popular benchmark problems for geometric nonlinear analysis of shells[END_REF].

  Fig. 20: Meshes corresponding to the points marked in Figs. 18 and 19 by the blackbordered marks. The point A marks the point where the load is applied in Fig. 17.The top row of meshes corresponds to the first limit point in Fig.19and the bottom row corresponds to the second limit points. The meshes are ordered from left to right for increasing solution steps. The elements are coloured according to the squared error e k normalized by the total error ∆L, i.e. e k /∆L 2 .

c 1 = 3 . 13 •Fig. 21 :

 131321 Fig.21: Geometry and parameters for the wrinkling problem. The boundary Γ 1 is free, the boundary Γ 2 is fixed in x and y direction and rotations around the y-axis are fixed, the boundaries Γ 3 and Γ 4 have symmetry conditions applied.

64 0 5 • 10 -

 64510 × 32, (ρ r , ρ c ) = (0.5, 0.005) uniform, 32 × 32 uniform, 64 ×

eFig. 23 :

 23 Fig. 23: Normalized element errors (left) and the normalized out-of-plane deformation (right) for the wrinkling benchmark plotted on the undeformed geometry with the corresponding elements.

Acknowledgments. H.M. Verhelst, M. Möller and J.H. Den Besten are greatful to the faculties of Electrical Engineering, Mathematics and Computer Science (EEMCS) and of Mechanical, Materials and Maritime Engineering (3mE) of Delft University of Technology for the financial support to conduct this research. A. Mantzaflaris @Angelos ADD

Compute the dual solution, ξ h given u h (Eq. 21)

5:

Compute the enriched dual solution, ξh given u h (Eq. [START_REF] D'angella | Efficient extraction of hierarchical B-Splines for local refinement and coarsening of Isogeometric Analysis[END_REF] 6:

Compute the total error estimation ∆L according to Eq. 23 and find the element-wise errors e k and the total element error e.

7:

if ∆L > tol r then 8:

Mark elements for refinement into M r using e, see Eq. 49.

9:

end if

Mark elements for coarsening into M c using e, see Eq. 53.

12:

Refine all Q ∈ M r using THB-splines, see Fig. 2 14:

Coarsen all Q ∈ M c using THB-splines, see Fig. 2 15:

Transfer the solutions required to start the new solution step to the new mesh using Quasi-Interpolation, see Eq. 56.

16:

end while Linear buckling analysis of a square plate (Section 6.3) Analytical critical buckling loads and mode shapes are also known for this case. Therefore, the goal of this benchmark problem is to provide verification for the buckling error estimator from Eq. ( 27) with the buckling operators

Non-linear analysis of a pinched thin plate (Section 6.4) In this example a thin plate with very low bending stiffness subject to a out-of-plane load is analysed. The error estimator is used to provide mesh adaptivity to compare to uniform refinement. The goal of this benchmark problem is to evaluate the performance of the DWR as driver for adaptive meshing in a static load case with geometric non-linearities. Snap-through instability of a cylindrical roof (Section 6.5) The snap-through behaviour of a cylindrical roof are considered in this example. The benchmark problem is a well-known application of arc-length methods and shells. The goal of solving this problem is to test the full adaptive solution stepping procedure from Fig. 5 on a benchmark problem. Wrinkling analysis (Section 6.6) In the last example, the procedure from Fig. 5 is applied to the modeling of membrane wrinkling. This problem contains geometric non-linearities and material non-linearities. The results are compared to uniformly refinements to evaluate the efficiency of adaptive meshing for such applications. The goal of this example is to demonstrate the use of the adaptive meshing routine from Fig. 5 on a complex load-stepping problem with geometric and material non-linearities.