Goal-Adaptive Meshing of Isogeometric Kirchhoff-Love Shells

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) [1] 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) 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 [2].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 [3][4][5][6], moderately thick Reissner-Mindlin shells [7][8][9][10][11] or thicker solid-like shells [12,13] 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 [14][15][16] or patch coupling methods [17][18][19][20][21] 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 the Partial Differential Equation (PDE) at hand, estimate element-wise error contributions, mark regions for refinement, refine (coarsen) marked regions [22], see Fig. 1.Here, localised regions can be defined element-wise or based on knot spans.The AIGM 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 [23].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 [24][25][26][27][28].
Solve The solve block contains the partial differential equation (PDE) at hand.It can be a physics-based problem, e.g., to solve shell [14,29,30], linear elasticity [31] or free-surface flow [32] problems.Alternatively, the solve step can involve a non-physics PDE, e.g. for mesh generation [33].Estimate Determination of localised errors is done in the estimate block.In the works [14,29,30], 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, e.g., by the Dual-Weighted Residual (DWR) method.This method has been applied in the FEA context in various works [34][35][36][37][38][39][40][41] and was used in the works [32,42] for Poisson and free-boundary problems, in [43] for a geometrically non-linear rod, in [33] PDE-based 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 regions for refinement based on a specific rule.The Refine block transforms the current mesh to a new mesh, where regions 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.domain parametrisations and in [31] for micromechanical modelling of trabecular bone.Goal-oriented refinement in general provides localised error estimates by solving a linear adjoint problem on the current space and an enriched space.Mark As soon as localised error contributions are known, regions can be marked for refinement.This marking is mostly done using the Dörfler marking strategy [44], as in [22,25,35], which involves marking the regions with the largest error contributions until their sum exceeds a certain percentage of the total error.An alternative is to mark the regions with an error higher than a threshold (an absolute threshold based on the maximum error) [29,45] or based on a relative threshold taking a fixed percentage of the total number of cells for refinement.In [45] 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) [46], Truncated Hierarchical B-splines (THB-splines) [45,47] or T-splines [48] amongst other spline constructions, which are reviewed in [49].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 [47].For (T)HB-splines suitable grading to generate admissible meshes should be taken into account in order to guarantee a bounded error [22], for which algorithms have been presented in [50].In the work of [51], 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 [22].Besides for adaptive meshing for solving PDEs [24][25][26][27][28]33], THB-splines have also been succesfully applied in the context of fitting [52,53].Transfer The transfer from previous time/load-steps onto a new mesh can be done using different methodologies.In the work of [27] different least-squares approaches are provided.Furthermore, quasi-interpolation [54,55] 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 [34,56] 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 [54,55].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 [3] 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 [59], 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 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: Given the second fundamental form b αβ = a 3 • a α,β = −a 3,β • a α and the metric coefficients defined as 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 ⊗g i , the deformation tensor C can be obtained: 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 one 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 [3,4]: (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: 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βδ [60].For non-linear 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 where T = [−t/2, t/2] is the through-thickness domain.For more details on hyperelastic material models, the reader is referred to [4,57] and specifically for stretch-based ones to [5].

Variational Formulation
The shell internal and external equilibrium equations in variational form are derived by the principle of virtual work [3,4].The weak formulation follows from the principle of virtual work with virtual displacements φ: 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 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 non-linear weak formulation Eq. ( 8) in a discrete setting [58]: where ε (u, φ, ψ) and κ(u, φ, ψ) are the second variations of the membrane and bending strains and 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 [3,58,60].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: 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 [61][62][63] 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 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 [64,65].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 [66] can be used.In our previous work, we provide more details on arc-length continuation for post-buckling analysis without providing a priori perturbations [67].
In the case of modal analysis and buckling analysis, a generalised eigenvalue problem needs to be solved.These eigenvalue problems have the general form 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: Where ρ is the density function over the surface.

Dual-Weighted Residual method
This section elaborates on the Dual-Weighted Residual (DWR) method [68,69], 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 [37,68,69].For the sake of completeness, we provide a brief overview of the DWR here.Consider the following non-linear problem to solve where W(u) is a semi-linear operator, u is the solution, φ 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) where u h and φ h are the discrete counterparts of u and φ, respectively, 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 while updating the discrete solution.From Eqs. ( 15) and ( 16), and the fact that φ h ∈ S p h ⊂ S the Galerkin orthogonality is found: where the last equality follows from the fact that the operator W(•, •) is linear in its second argument.Let us now define a non-linear and differentiable goal functional L(u) or quantity of interest, such that and let the residual be defined as If ξ ∈ S p h is the solution to the following dual problem with L (u, ξ) being the linearisation of the goal functional L(u) with discrete counterpart then the following error representation formula can be obtained (see Theorem 4.1 of [37]): 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.
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 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 [35,41,56,[68][69][70] to give a brief overview of the DWR for eigenvalue problems.
Let us consider the following eigenvalue problem Here, A and B are bi-linear operators.Typically, discretizing the system gives the following: where the eigenpairs vh = (µ h , v h ) are the solutions of the eigenvalue problem.For uniqueness of the problem, the discrete eigenvectors v h are normalised by the condition To derive the DWR method for the eigenvalue problem in Eq. ( 26), the functional V(•, •) is defined, such that the following problem should be solved where the normalisation condition from Eq. ( 28) is enforced weakly.The discrete counterpart of this equation reads: Furthermore, a goal-function for the eigenvalues is defined as follows: giving Using the non-linear 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: 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 such that the adjoint eigenvalue problem becomes (35) which is equivalent to solving [56] Find Assuming the operator B is self-adjoint, it follows that the adjoint eigenvalue problem is solved by η = µ and when normalizing the adjoint eigenvectors ψ using It can be seen that the adjoint solution (η, ψ) is obtained by solving the eigenvalue problem using Eqs.( 15), ( 23) and (32) 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 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 [38] different choices for constructing Sp h are given, including an h-refinement and a p-refinement.As in the work of [33], 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 .

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: 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.
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 [25,50]

(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 [46,47].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 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 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 [71].that for THB-splines a truncation is performed to ensure partition of unity, which is discussed in more detail in [47].

Admissible Meshing
The concept of admissible meshing was discussed in [22,23,50].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 [22,23,50].In Fig. 3, we illustrate a simple mesh together with the refinement neighborhood of some selected elements.The T- 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. [25].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 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 in [25], 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.
with Q ∈ S(Q , − 1) .(45) 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 coarse labeling can be performed conforming with the 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 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 [44] is defined as the elements Q i ∈ Q, i = 0, ..., k such that the sum of their respective errors is smaller than a fraction ρ r of the total element error e = i e i : 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 whose span contains elements Q k ∈ Q and their refinement neighborhoods: 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: such that the Dörfler marking including refinement neighborhoods is 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 : 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 .
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: such that the Dörfler marking strategy taking into account coarsening admissibility is defined as 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.: 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 [54,55].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 where the coefficients λ i, are suitable linear functionals on C(Ω 0 ).Across all levels = 0, . . ., N − 1, the interpolant for the function becomes: 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 summarised workflow for adaptive isogeometric shell analysis is depicted in Fig. 5 and Algorithm 1.
The Solve block involves solving the non-linear 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. (54).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. 21  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.

5:
Compute the enriched dual solution, ξh given u h (Eq.24) 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.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 A(v, φ) = W (u L , v, φ) and B(v, φ) = M(v, φ).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.
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 [71].

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 is applied, based on the manufactured solution given by 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 spline 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
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 where γ n is the n th root of the equation ( following from a separation of variables solution [67], 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 [72] and reads: 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: 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 normalised with the analytical value of the critical buckling load.As can be seen in this figure, 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 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 is similar to the behaviour observed in [73] 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 generalised eigenvalue problem for buckling.For the results presented in Fig. 11, the load σt = 10 −4 [N/m] was used to compute the buckling linearisation.

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.
case are based on displacements as well as on principal stresses: 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 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 plots 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 it is most effective to increase the maximum refinement depth.Comparing the error fields and meshes for 10 −7

#DoFs
Adaptive (8) 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 and 16.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 and 16), 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 show 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 with moderate computational costs.

Snap-through instability of a cylindrical roof
In order to present results of the adaptive isogeometric method developed in this paper for quasi-static problems, hence completing full cycles in Fig. 1, the well-known benchmark of a collapsing cylindrical roof [74] 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 t = 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 [65] with arc-length ∆L = 25 and a zero force-scaling.The goal functional that is used for error evaluation and adaptivity is 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, the coarsening parameter to ρ c = 0.05 and the maximum refinement level is 11.The adaptive and uniform meshes are modelled with bi-cubic B-spline basis functions (i.e.p = 3).
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 [74] 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 = 0.2, 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 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 uniform mesh with 32 × 32 elements.In Fig. 20 a selection of meshes is provided.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.18 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 and final 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 [75][76][77] and was previously modelled using isogeometric Kirchhoff-Love shells in [5], 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 and with parameters c 1 = 3.16 • 10 5 [Pa] and c 1 = 1.24 • 10 5 [Pa].Reference solutions are given by [5] for isogeometric Kirchhoff-Love shell analysis and using ANSYS and LS-DYNA FEA models, respectively.Furthermore, [77] 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  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 normalised 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 due to 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 tol r = 10 −14 and tol c = 10 −10 , respectively, and they 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 re-stabilisation 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 [65] with a quadratic procedure to compute the mode shape at the bifurcation [66].This procedure further described in [67].
In Fig. 22, the results for an adaptive isogeometric wrinkling simulation are provided.The top figure provides the normalised wrinkling amplitude with respect to the strain of the sheet ( ), compared to the IGA and ANSYS SHELL181 element reference results from [5] as well as the experimental results from [77].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 re-stabilisation (i.e. the moment when the amplitude vanishes).After re-stabilisation, 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 initiates, 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 re-stabilisation 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 re-stabilisation 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 optimised in this study.

# 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, the dotted line in the top figure is a result obtained using the fully integrated shell in LS-DYNA and the dashed line with markers represents experimental data obtained by [77].

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 and 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 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 generalised eigenvalue problem to be solved on the enriched mesh.For suitable grading, the works of [22,23,25,50] 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 functionals 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

Fig. 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 :
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 :
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.
if ∆L < tol c then 11:

17 :
Advance the solution to the next solution step 18: 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 :
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].

12 . 10 − 2 Fig. 11 :
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 sise 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].

Fig. 15 :
Fig. 15: Normalised 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.

Fig. 16 :
Fig. 16: Normalised 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 :
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 :
Fig.19: Projection the adaptively refined result from the commonly used λP, w Aspace[74] 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[74].
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 normalised by the total error ∆L, i.e. e k /∆L 2 .

c 1 = 3 . 13 •Fig. 21 :
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.