A Consistent Higher-Order Isogeometric Shell Formulation
AA Consistent Higher-Order Isogeometric ShellFormulation
Daniel Schöllhammer , Benjamin Marussig , Thomas-Peter Fries December 23, 2020 , Institute of Structural AnalysisGraz University of TechnologyLessingstraße 25/II, 8010 Graz, Austria Institute of Applied MechanicsGraz University of TechnologyTechnikerstraße 4/II, 8010 Graz, Austria [email protected] , [email protected] , [email protected] a r X i v : . [ c s . C E ] D ec Abstract
Shell analysis is a well-established field, but achieving optimal higher-order con-vergence rates for such simulations is a difficult challenge. We present an isogeomet-ric Kirchhoff-Love shell framework that treats every numerical aspect in a consistenthigher-order accurate way. In particular, a single trimmed B-spline surface providesa sufficiently smooth geometry, and the non-symmetric Nitsche method enforces theboundary conditions. A higher-order accurate reparametrization of cut knot spansin the parameter space provides a robust, higher-order accurate quadrature for (mul-tiple) trimming curves, and the extended B-spline concept controls the conditioningof the resulting system of equations. Besides these components ensuring all require-ments for higher-order accuracy, the presented shell formulation is based on tangentialdifferential calculus, and level-set functions define the trimming curves. Numericalexperiments confirm that the approach yields higher-order convergence rates, giventhat the solution is sufficiently smooth.
Keywords: Kirchhoff-Love shells; higher-order; Trimming; IsogemetricAnalysis;
ONTENTS Contents R . . . . . . . . . . . . . . . . . . . . . . . 314.3 Clamped circular shell . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 34 Introduction
Shells are components of utmost importance in numerous engineering applications, andhence, the development of more accurate mechanical models and numerical analysis schemesis a long-established [8, 84], yet everlasting [13, 112], challenge. The finite element method(FEM) is the most common technique for the numerical treatment of shell formulations.In this context, isogeometric analysis (IGA) [60] has proven itself as an efficient and robustparadigm that outperforms traditional structural mechanics simulations [28, 68, 75]. Theunderlying idea of IGA is quite simple: using the functions of spline-based geometry rep-resentations of computer-aided design (CAD) directly for performing analysis. Althoughthis concept was – and still is – aimed at mitigating the profound inefficiencies in the con-ventional interaction of CAD and simulation tools, it turned out that IGA enables severalcomputational benefits. The most important advantages for shells are that it provides (i)the most accurate representation of the geometry and (ii) high continuity between ele-ments. Indeed, the latter rejuvenated research on Kirchhoff-Love (KL) shell formulations,because establishing the required C continuity is straightforward using splines.Since the pioneering work for KL shells from Kiendl et al. [68], IGA shell analysis gaineda notable amount of research activities. The curvilinear coordinates, which are implied bythe map from the parameter space to the physical space of the NURBS patch, and thehigh continuity are a perfect fit for the discretization of the well-known shell equationswith IGA. It is worth noting that IGA has also been considered for Reissner-Mindlin (RM)shells, which take transverse shear deformations into account. A successful discretization ofRM shells with isogeometric analysis is given in, e.g., [10]. Furthermore, an approach withexactly calculated directors in the frame of IGA is presented in [35, 36]. Hierarchical shellsthat eliminate transverse shear locking are shown in [38, 86, 87]. A mixed displacementapproach that also eliminates membrane locking is elaborated in [12, 114]. Isogeometriccollocation methods for KL shells and RM shells are considered in [69, 79].In this paper, we employ the rotation-free Kirchhoff-Love shell formulated in the frameof the tangential differential calculus (TDC) [97, 98]. The resulting shell formulation isindependent of the concrete geometry definition and enables shell analysis on implicitlydefined surfaces where a parametrization is not available nor needed, see [99]. However,herein, the midsurface is explicitly specified by a single NURBS patch. The benefits ofthe TDC-based shell formulation in our context are that the obtained shell equations areformulated in the global Cartesian coordinate system and may be presented in a morecompact and intuitive fashion employing a symbolic notation, while the classical formula-tion is generally expressed in curvilinear coordinates using index notation. As a result, thetypically rather cumbersome effective boundary forces, see, e.g., [11], needed later on forNitsche’s method, are more compact, which simplifies the implementation.The desire to perform shell analysis on practical engineering CAD models leads to the chal-lenge of trimmed geometries. Trimming refers to a procedure that is used to specify visibleareas of an object independent from the object’s parameter space, and it is paramount forfundamental geometrical operations denoted as Boolean operations (union, intersection,and differencing). This concept provides great flexibility in defining arbitrary shapes overspline surfaces but introduces several computational issues as detailed in [77]. In essence,trimming results in cut parameter spaces, and thus, the analysis of trimmed geometriesfaces the same challenges as fictitious domain methods (FDMs) [20, 59, 95]. That is, thepresence of cut elements complicates1. the application of boundary conditions,2. the conditioning of system matrices, and3. accurate integration.In the context of IGA with trimmed geometries, the RM shell theory was adopted in [65–67, 111]. The first KL shell formulation was presented in [96]. The work of Breitenbergeret al. [14, 15] develops the concepts further, enabling direct use of CAD data for thesimulation. While these works rely on a penalty approach for the application of boundaryconditions and a boundary-fitted quadrature scheme, the contributions of Guo et al. [50,51, 53, 54] employ the symmetric or non-symmetric Nitsche’s method and a non-boundary-fitted sub-cell quadrature; except for [51] where this quadrature is used only for cut elementsthat cannot be handled by the blending function method presented in [70]. The sub-cellquadrature originates from a FDM and includes the application of a fictitious stiffness in theexterior domain to mitigate the ill-conditioning issues due to cut elements. The literatureon isogeometric formulations for trimmed shells barely covers this conditioning aspect.Besides the fictitious stiffness, Coradello et al. [27] use diagonal scaling for preconditioning.As mentioned above, trimming and FDMs are closely related from a computational pointof view. For instance, membranes have been introduced for Cut FEM [21] and for TraceFEM [45]. Furthermore, linear RM shells with Trace FEM have been presented in [99]. A Introduction recent pre-print where a C -continuous Trace FEM approach for KL shells with boundaryconforming background meshes is presented in [46].The motivation of the present work is the realization of a higher-order KL shell analysis,which is an “open” topic because the requirements (also on the geometry) are high. Whilesetting up higher-order discretizations is straightforward using IGA, it is important to notethat achieving higher-order accuracy involves more than the utilization of basis functionswith a high polynomial degree. The analysis has to consider all numerical aspects ina consistent higher-order fashion: starting from the discretization over the applicationof boundary conditions to the construction of the system matrices. To the best of theauthors’ knowledge, there exists no finite element shell model in the literature that canclaim complete success in this regard. In fact, there is no higher-order accurate shellformulation in the context of trimming. It is, however, noted that the already mentionedworks from Fries et al. [45] on membranes and Schöllhammer et al. [99] on linear RM shellsare presented in the framework of higher-order Trace FEM.An aggravating factor for the development of higher-order shell formulations is that thetraditional benchmarks tests [9], which mark the current gold standard for validating shellmodels, do not allow for higher-order convergence rates. First of all, the solutions of thesebenchmarks are not smooth enough (e.g., due to point forces or geometric corners), andsecond, the reference solutions are often reported with only a few numbers of significantdigits (see e.g., [68]). A recent work [11] demonstrates that these benchmarks are “un-able to assess the order of accuracy” and maybe passed by even variationally inconsistentformulations.In this paper, we present a higher-order accurate KL shell formulation based on the iso-geometric paradigm with trimmed geometries. Trimming allows the definition of arbitraryshapes on single-patch representations. We do not consider multi-patch models becausethey are inherently low-order due to the non-smoothness of the solution (across edges be-tween patches, which is well-documented in KL-shell theory). In general, any corner inthe shell’s boundary – even for single-patch geometries – may typically result in singulari-ties, immediately hindering higher-order convergence rates in the analysis. Using trimmingalso provides a means to prevent such discretizations. As listed above, simulations withtrimmed geometries involve three key challenges, which we address in a consistent higher-order accurate form. In particular, the following methods are utilized:1. The non-symmetric Nitsche method enforces the boundary conditions.2. The extended B-spline concept controls the conditioning.3. Higher-order Lagrange elements provide a reparameterization of cut elements whichallows for proper distribution of integration points.It is worth noting that different approaches for dealing with each of these problems exist.They may be used instead of the proposed ones, as long as they can establish higher-order accuracy; if not, the overall convergence will be sub-optimal. We are using level-setfunctions for the definition of the trimmed domain and formulate the KL shell in theframe of tangential differential calculus (TDC) [97, 98]. These two particular choices are,however, no necessary ingredients for the realization of higher-order accuracy and may besubstituted by alternative concepts preferred by the reader.Three numerical examples validate the performance of the proposed KL shell model: (i) atraditional shell benchmark, (ii) the flat shell embedded in R from [97], and (iii) a curved,clamped circular shell. The latter addresses the three trimming challenges, as well as therepresentation of the geometry and the boundary conditions. Hence, it allows an overallassessment of the order of accuracy. In particular, we measure the error of this problemin the relative L -norm of the equilibrium in strong form, which requires the evaluation offourth-order surface derivatives of the midsurface displacement. Consequently, the theoret-ical optimal order of convergence is O ( p − Geometry representation and treatment of trimmed domains B B B C C
111 2 3 444 (a)
Polynomial segments B s of a B-spline B , B , B , B , B ,
111 2 3 444 (b)
B-spline basis
Fig. 1:
Quadratic B-splines defined by Ξ = { , , , , , , , } : (a) polynomial segments of thefunction B , , which is part of (b) the entire B-spline basis. In (a), the resulting B-splineis illustrated by solid lines, whereas dotted lines indicate extensions of its segments. We are aiming for a surface discretization that allows for (i) at least C continuity, (ii)arbitrary degree, and (iii) arbitrary boundaries. The continuity requirement stems fromthe Kirchhoff-Love shell theory (see Section 3.1), and high-degree basis functions are aprecondition for higher-order approximation power. Splines are a natural choice to fulfillthese two properties, but regular spline surfaces are defined by a tensor product structure,which results in a four-sided topology. Besides the related geometric limitation, it is alsonoted that any corner in the shell boundary typically results in singularities, which mayhinder the optimal convergence. Hence, the concept of trimming is used to specify arbitraryshapes over spline surfaces. To be precise, the discretization employs trimmed B-splinetensor product splines. In the following, we recall the basic properties of these functions.The interested reader is referred to [30, 77, 89] for a detailed discussion.A B-spline B i,p is described by piecewise polynomial segments B si of degree p , as shownin Figure 1(a). These B si and the continuity between them are specified by the knotvector Ξ , which is a non-decreasing sequence of parametric coordinates ξ j (cid:54) ξ j +1 . Thevalues of these knots ξ j define the location where adjacent B si join. The continuity at these breakpoints is C p − m , with m denoting the multiplicity of the corresponding knot value, i.e., ξ j = ξ j +1 = · · · = ξ j + m − . The knot span s refers to the half-open interval [ ξ s , ξ s +1 ), and ifits size is non-zero, it marks the valid region of B si . These regions also impose the elements .1 Geometry representation and discretization Ξ defines not only a single function, but an entire set oflinearly independent B-splines { B i,p } ni =0 , as illustrated in Figure 1(b). Each B i,p has localsupport, supp { B i,p } , specified by the knots { ξ i , . . . , ξ i + p +1 } , and each knot span s contains p + 1 non-zero B-splines represented by B si with i = s − p, . . . , s .Basis functions of B-spline surfaces are obtained by computing the tensor product of uni-variate B-splines B i,p and B j,p of degrees p and p , which are defined by separate knotvectors Ξ and Ξ for the parametric directions ξ and ξ , respectively. A correspondingsurface is given by Γ( ξ ) = I − X i =0 J − X j =0 B i,p ( ξ ) B j,p ( ξ ) c i,j (2.1)where ξ is a vector containing all parametric coordinates ξ k with k = { , } , and c i,j arecontrol points in physical space R . Following the isogeometric paradigm, we employ thesame basis also for the test and trial functions. Later, the surface’s unit normal vector n Γ plays an essential role. It is determined by the cross-product n Γ = t × t k t × t k (2.2)of the tangent vectors t ( ξ ) = I − X i =0 J − X j =0 B (1) i,p ( ξ ) B j,p ( ξ ) c i,j (2.3) t ( ξ ) = I − X i =0 J − X j =0 B i,p ( ξ ) B (1) j,p ( ξ ) c i,j (2.4)where B (1) denotes the basis function’s first derivative. Furthermore, we introduce theco-normal vector n ∂ Γ along the surface’s boundary ∂ Γ, which is given by n ∂ Γ = n Γ × t ∂ Γ (2.5)with t ∂ Γ being the tangent vector along the boundary curve. Note that Equation (2.5)guarantees that n ∂ Γ points out of the surface and lies in its tangent plane. The resultinglocal triad ( t ∂ Γ , n ∂ Γ , n Γ ) will be used to properly represent the boundary conditions ofthe shell.Representation (2.1) allows full control over the continuity and degree of the discretiza-0 Geometry representation and treatment of trimmed domains Γ( ξ ) (a) Regular surface C t A v ξ ξ (b) Trimmed para-metric space Γ( ξ ) (cid:12)(cid:12) A v (c) Trimmed surface
Fig. 2:
Trimmed surfaces: (a) regular surface as defined by Equation (2.1), (b) the relatedtrimmed parametric space where trimming curves (thick lines) delineates the visiblepart A v of (c) the resulting trimmed surface considered for the analysis. tion, but it yields surfaces with an intrinsic four-sided topology due to its tensor productstructure, as illustrated in Figure 2(a). Trimming provides a remedy to this restrictionby defining the visible area A v of arbitrary shape over a surface Γ( ξ ). As indicated inFigure 2(b), curves in the parametric space specify A v . Usually, these trimming curves C t are also represented by B-splines [77], where the curve’s orientation determines whichregions are visible. In this paper, we follow a different approach and define a C t as thezero-isoline of a level-set function φ ( ξ ). This choice is mainly motivated by the simplifieddetection of A v during the analysis, i.e., A v := { ( ξ ) | φ ( ξ ) ≥ } . (2.6)We want to emphasize that the chosen representation of C t has no impact on the proposedmethodology for achieving higher-order accuracy.It remains to define the tangent vector along the boundary t ∂ Γ when ∂ Γ is given by atrimming curve ˜ C t . In this case, t ∂ Γ is no longer aligned with a parametric direction ofthe surface. Hence, we need a curve representation of C t in physical space, ˜ C t ∈ R , tocompute t ∂ Γ . It is common practice to represent ˜ C t by an approximation because theexact description by the expansion of Γ ◦ C t is of degree p t ( p + p ) where p t refers to thedegree of C t , which is usually un-practically high. As shall be detailed in Section 2.3, knotspans cut by a trimming curve are decomposed into higher-order Lagrange elements. Theseelements are used to define quadrature rules in the trimmed domain and, furthermore, theiredges provide higher-order accurate approximations of ˜ C t , which we use to evaluate thetangent t ∂ Γ . The normal vector n Γ and the co-normal vector n ∂ Γ , on the other hand, are .2 Conditioning considerations Analysis schemes that employ cut elements may suffer from sever conditioning problems.In particular, ill-conditioning occurs when only a small portion of a basis function’s supportcontributes to the overall system, which leads to small eigenvalues of the system matrix.For a detailed analysis of the conditioning issue, the interested reader is referred to [34].There are several remedies to this problem proposed in the literature. The perhaps mostsimplistic approaches are diagonal scaling, e.g., [2, 27], and the use of a fictitious domainstiffness, e.g., [29, 95]. They are appealing due to their ease of implementation, but theydo not address the actual source of the conditioning problems. More adequate strategiesare tailored preconditioning, e.g., [33, 34], the ghost penalty, e.g., [16, 18, 20], and themodification of cut basis functions, e.g., [59, 62, 76, 78, 80]. Each of these different methodshas its advantages, but a detailed comparison is beyond the scope of this paper. In thefollowing, we focus on the latter approach as detailed in [78], as it is particularly well-suited for spline bases. This so-called extended B-spline concept removes basis functionswith small support and replaces them by extensions of neighboring basis functions. This isobviously a change of the approximation space, however, optimal approximation power ismaintained, see [59]. As will be shown shortly, the extension procedure is done algebraicallyand is hence independent of the problem formulation and the structure of the resultingsystem matrix. The subsequent sections provide an introduction to extended B-splinesand their application to a system of equations, while appendix A gives details on theimplementation.
First, we classify the basis functions of a trimmed basis as either stable , exterior , or de-generate , where the latter is referring to those that may cause conditioning problems. Therelated identification utilizes the Greville abscissa (see [76, 78]) or the intersection of thesupport with the domain, i.e., S v i := supp { B i,p } ∩ A v . Here, we employ the latter strategyand label a B-spline B i,p as• Stable if S v i ≥ α supp { B i,p } ,2 Geometry representation and treatment of trimmed domains • Degenerate if S v i < α supp { B i,p } ,• Exterior if S v i = ∅ ,where α ∈ (0 ,
1] is a user-defined threshold for the minimal relative support size consideredto be stable. Usually, we choose α between 0 . . extensions of stable ones. These extensionsare provided by the polynomial segments B i of the closest element that contains only stableB-splines. The accumulation with all other stable functions yields the final extended B-spline basis. B B B B B A v (a) B-splines of the initial basis S v0 < α supp { B } ⇒ J = { }A v t (b) Detection of degenerate B-splines B B B A v t (c) Extensions of polynomial segments B e B e B e B A v t (d) Extended B-splines
Fig. 3:
Basic procedure to get from (a) conventional to (d) extended B-splines: (b) determi-nation of degenerate B-splines and substitution of trimmed polynomial segments by (c)extensions of non-trimmed ones. In (b), J is the index set of all degenerate B-splines. In this construction, the extended polynomial segments B si of the non-trimmed knot span s are defined by the initial B-splines of the knot span t that contains degenerate functionsby B si ( ξ ) = t X j = t − p B j,p ( ξ ) e i,j ξ ∈ [ ξ t , ξ t +1 ) . (2.7) .2 Conditioning considerations extrap-olation weights e i,j . If the index j corresponds to a stable B-spline, the weights are trivial,i.e., e i,i = 1 ⇔ B si ( ξ ) ≡ B i,p ( ξ ) , ξ ∈ [ ξ s , ξ s +1 ) , ∀ i ∈ { s − p, . . . , s } , (2.8) e i,j = 0 ⇔ B si ( ξ ) = B j,p ( ξ ) , ξ ∈ [ ξ s , ξ s +1 ) , ∀ j ∈ { s − p, . . . , s }\ i. (2.9)The remaining extrapolation weights related to degenerate B-splines can be determined by e i,j = 1 p ! p X k =0 ( − k ( p − k )! β p − k k ! α k (2.10)with the coefficients α and β denoting the constants of the polynomials B si ( ξ ) = p X k =0 α k ξ k and ψ j,p ( ξ ) = p Y m =1 ( ξ − ξ j + m ) = p X k =0 β k ξ k . (2.11)The former is the power basis form of the extended polynomial segment, and the Newtonpolynomials ψ j,p result from a quasi interpolation procedure, i.e., the de Boor–Fix func-tional [30, 31]. For more details, we refer to [59, 78] and appendix A which provides furtherinformation on the determination of the extrapolation weights (2.10).By taking all extrapolations weights, see Eqs. (2.8)–(2.10), into account, an extended B-spline is defined as B ei,p = B i,p + X j ∈ J i e i,j B j,p (2.12)where B i,p is the stable B-spline from which the extension originates, and J i is the index setof all degenerate B-splines related to B ei,p . Definition (2.12) applies to multivariate basisfunctions as well. The corresponding extrapolation weights are obtained by the tensorproduct of the univariate components. Remark . The degree of the basis affects the extended B-splines regarding the balance ofimproved conditioning and quality of the solution along the trimming curve. The parameter α provides a certain level of control, but the combination with local refinement, e.g., assuggested in [76], is a more comprehensive solution which is, however, not utilized in thispaper.4 Geometry representation and treatment of trimmed domains
A very convenient feature of extended B-splines is that they do not need to be explicitlyevaluated during an analysis. Suppose we set up the following linear system of equationsusing the initial – potentially ill-conditioned – B-splines basis K u = f where u ∈ R m , f ∈ R n and K ∈ R n × m with m > n. (2.13)Here, n refers to the number of stable B-splines, while m includes the degenerate onesas well. In order to get a well-conditioned square system matrix, K is multiplied by an extension matrix E ∈ R m × n [59]. The entries of E are the extrapolation weights e i,j of allextended B-splines. The trivial weights e i,i = 1 are stored as well, even if a stable B-splinehas no related degenerate ones. The matrix entries of E are assembled such that columnsof the i th row of K are distributed according to the definition (2.12) of the associatedextended B-spline B ei,p . The final system of equations due to the extended B-spline basisis determined by K st u st = f with K st = KE , K st ∈ R n × n . (2.14)The solution vector u st ∈ R n corresponds to the extended B-spline basis. Its relation tothe original basis can also be expressed by the extension matrix as u = E u st . Remark . The example given above considers a scalar problem. The generalization tovector-valued problems merely involves an appropriate repetition of the entries of E . In isogeometric analysis, the integration of the weak form is classically performed throughan outer loop over all knot spans with an inner loop over quadrature points. Reducedquadrature concepts following e.g. [5, 56, 61, 63, 94] use a tailored quadrature for thetensor-product B-splines or NURBS bases and make a profit of the increased continuity ofthese bases. Thus, they have shown a significantly higher efficiency in standard GalerkinIGA analyses. However, their application in the context of trimming is delicate as trimmedbasis functions lack the required continuity over the whole support being a fundamentalrequirement of the reduced quadrature concepts.Hence, we follow the classical path and use standard Gauss quadrature in knot spans .3 Numerical integration in cut knot spans C t requires special attention. This is closely relatedto the context of fictitious domain and extended finite element methods where domainboundaries or interfaces freely cut the elements and various alternatives for the integrationin such elements exist. They differ in their ability to enable higher-order accuracy andcapture complex geometries including corners.The first approach to the quadrature in cut elements relies on the decomposition intointegration cells. Those may be polygonal cells often obtained in recursive refinements ofthe cut element until the desired accuracy is obtained [1, 37, 81]. Nevertheless, the numberof integration cells is often immense, even more so in a higher-order context [37, 72, 74,107, 113] which has direct implications to IGA. It was already noted in [25, 42, 73, 92]that a decomposition into sub-elements with curved, higher-order edges is a strategy thatconsistently enables the generation of higher-order accurate quadrature rules with only amodest number of integration points. This approach has been investigated thoroughly in[41–43, 88] and is followed herein. Other higher-order alternatives based on decompositionsare reported in [71, 106]. Methods that do not decompose the cut elements are, e.g., foundin [64, 82, 83, 109, 110]. There, the number of integration points is also small, however,the generation of proper integration weights is often involved and requires the solution ofsmall systems of equations, e.g., in the context of moment fitting and Lasserre’s technique.For the standard context where the trimming curve is defined by means of a B-splineor NURBS in the parameter space [77, 102], this explicit definition is converted to animplicit level-set function φ ( ξ ) by computing signed distances which is called implicitation [24, 57, 58, 101, 103, 105]. The trimming curve may also be directly given in the implicitform. In both cases, the level-set function is positive in the domain of interest, φ ( ξ ) > φ ( ξ ) <
0. Most importantly, the trimming curve coincides with thezero-isoline of φ ( ξ ) in the parameter space, i.e., C t = ∂ Γ = { ξ : φ ( ξ ) = 0 } .Here, we adopt the procedure outlined in [41–43, 88]. The basic idea is to first reconstruct the zero-isoline by curved, higher-order, interpolatory, and one-dimensional finite elementsand then decompose the cut knot spans into higher-order, two-dimensional finite elementsserving as integration cells in the cut knot spans.6 Geometry representation and treatment of trimmed domains - +0
Fig. 4:
Lagrange elements are placed in cut knot spans of the parameter space. The level-setdata is evaluated at the nodes in order to reconstruct integration cells.
Starting point is a level-set function φ ( ξ ) in the parameter space whose zero-level set ∂ Γ isthe trimming curve. Then, the following procedure is applied to each knot span: The nodes ξ FE k of a two-dimensional, quadrilateral Lagrange finite element of order p are defined in thecurrent knot span, see Figure 4. These are n FE = ( p + 1) nodes, hence, k = 1 , . . . , n FE and p matches the order of the underlying B-spline (or NURBS) basis in the parameter space.One may use equally spaced nodes [42, 43] for simplicity or special node distributions suchas Gauss-Lobatto or Chen-Babuška nodes [22, 23, 41]. The level-set function is evaluated atthe nodes of the finite element resulting in φ k = φ (cid:16) ξ FE k (cid:17) . For explicit trimming curves, suchfunction evaluations are basically signed distance computations at each node. Throughthe basis functions of the element, N FE k ( ξ ), these nodal values define an approximation φ h ( ξ ) = P n FE k =1 N FE k ( ξ ) · φ k of the exact function φ ( ξ ), see Figs. 5(a) and 6(a). Also thecorresponding zero-level set ∂ Γ h of φ h ( ξ ) is a higher-order accurate approximation of ∂ Γ .Next, sample grid points ξ grid m are introduced in the knot span, typically related to theorder p such that there are n grid , = p · ( q + 1) + 1 nodes in each direction of the knotspan resulting in a total of n grid = (cid:16) n grid , (cid:17) grid points, hence, m = 1 , . . . n grid . q ∈ N +0 is a user-defined number related to the complexity (e.g., curvature) of φ ( ξ ) and is oftenchosen as q = 3. The level-set function φ h ( ξ ) in the element is evaluated at the grid pointsgiving φ grid m = φ (cid:16) ξ grid m (cid:17) = P n FE k =1 N FE k (cid:16) ξ grid m (cid:17) · φ (cid:16) ξ FE k (cid:17) , see Figs. 5(b) and 6(b).Then, the cut situation in the element is determined based on the level-set data at thesample grid points φ grid m . The situation is called valid if (i) each element edge is only cutonce, (ii) the overall number of cut edges is two, and (iii) if no edge is cut then the element .3 Numerical integration in cut knot spans (a) level-set function φ h (b) sample grid (c) reconstruction (d) decomposition Fig. 5:
Two neighboring edges of a higher-order finite element ( p = 3) located in a knot spanare cut (topology type 1): (a) the level-set function φ h implied by the nodal level-setdata, (b) the level-set data sampled at the grid points ( q = 3), (c) the reconstructedone-dimensional finite element representing ∂ Γ h , (d) the quadrature cells with resultingintegration points. (a) level-set function φ h (b) sample grid (c) reconstruction (d) decomposition Fig. 6:
Two opposite edges of a higher-order finite element ( p = 3) located in a knot span arecut (topology type 2): (a) to (d) as in Figure 5. Geometry representation and treatment of trimmed domains is completely uncut. Otherwise, the situation is invalid , however, a valid situation maythen always be recovered by recursive (quad-tree like) refinements of the element. Forfurther details, we refer to [42, 43]. For valid situations, each (possibly refined) elementfalls into two topologically different cases: When two neighboring edges are cut (case 1),there results a triangle and a pentagon, and for two opposite edges being cut (case 2),there result two quadrilaterals on each side, see Figs. 5(d) and 6(d).
Remark . It is noted that finite element operations such as the evaluation of shape func-tions are evaluated in the (quadrilateral) reference element. However, the map to eachknot span is affine, even linear in each direction ξ i . Therefore, we avoid to complicate thesituation notationally by distinguishing finite element operations in the reference elementversus the knot span. Remark . The described procedure may be generalized in various ways. Instead of ap-proximating a given level-set function φ using interpolatory finite elements, one could alsouse the B-spline (or NURBS) basis in the parameter space to determine φ h based on an L -projection. The boundary of the trimmed shell geometry would then be implied bylevel-set values at the control points. Let us assume that the nodal level-set data in a cut knot span is valid in the above sense.First, the two intersections ξ ? and ξ ? of the level-set function φ h ( ξ ) with the two relevantelement edges are determined for which a Newton-Raphson procedure is employed, seeFigs. 5(c) and 6(c). Next, p + 1 nodes are distributed along the straight line going through ξ ? and ξ ? yielding starting points for another Newton-Raphson loop which, for each innernode on the line, finds the intersection with the zero-level set of φ h ( ξ ) using the gradient ∇ φ h at the starting points as the search direction, see Figs. 5(c) and 6(c). Hence, the p + 1 nodes are placed on ∂ Γ h up to some user-define tolerance, usually in the range of10 − . . . − . These nodes imply a curved, one-dimensional element being a higher-orderapproximation of ∂ Γ h (and ∂ Γ ). Concerning the reconstruction, many other choices ofthe start values and search directions are possible, see [42, 43] for details.The reconstructed element representing the zero-isoline and, hence, the trimming curve, isalso used to determine normal and tangent vectors of the boundary in the parameter space.After the isogeometric map to the real shell geometry in R , this supports the definitionof the full triad ( t ∂ Γ , n ∂ Γ , n Γ ) on the shell boundary, see Figure 7. .3 Numerical integration in cut knot spans (a) situation in cut knotspan (b) situation in R Fig. 7:
Normal and tangential vectors at the trimming curves in (a) parameter space and (b)the shell geometry in R . The final task is to decompose the element in the knot span based on the reconstructed one-dimensional element representing the zero-isoline. The resulting two-dimensional elementsdepend on the topology of the cut. For case 1, where two neighboring edges are cut, thetwo sides are decomposed into four triangles featuring straight edges except for the oneedge coinciding with the reconstructed zero-isoline, see Figure 5(d). For case 2, wheretwo opposite edges are cut, the two sides fall into two quadrilateral elements, each withone curved edge being the reconstructed zero-isoline, see Figure 6(d). Standard Gaussrules in triangles and quadrilaterals are mapped to the decomposed elements that areon the positive side of φ h . Those on the negative side are neglected. Concerning themapping of integration points to the decomposed elements (with one curved, higher-orderside representing ∂ Γ h ), several approaches may be used [47–49, 100, 104, 108]. Herein,we use the transfinite maps suggested by Šolín in [104] as they apply to all finite elementshapes. See Figs. 5(d) and 6(d) for the resulting integration points, only the red ones havepositive level-set data and are considered; the blue ones are neglected as they are outsidethe domain of interest.It is noted that the reconstruction and decomposition may also fail for valid level-set dataat the sample grid points. For example, the Newton-Raphson loops for placing nodes onthe zero-isoline may not converge or result in points outside the knot span or the Jacobiansin the decomposed elements may be negative. Such cases are then also treated by recursiverefinements until the reconstruction and decomposition are successful, see [42, 43]. Thistypically only applies to complex trimming curves and only in very few knot spans.0 Governing equations
In the following, the fundamentals of the employed shell model are presented. In particular,the rotation-free KL shell [97, 98] formulated in the frame of the tangential differentialcalculus (TDC) is used herein. It has been shown in [99] that the TDC-based formulationis more general than the classical outline based on curvilinear coordinates as it also appliesto implicitly defined shell geometries. Anyhow, herein, the midsurface is explicitly definedby a single NURBS patch, see Section 2.1, and the classical shell formulations as presentedin, e.g., [8, 13, 68], are also applicable and lead to equivalent results. Nevertheless, also inthe present context, one may find the shell formulation based on the TDC also beneficialbecause the obtained shell equations are formulated in the global Cartesian coordinatesystem and may be presented in a more compact and intuitive fashion employing symbolicnotation whereas the classical formulation is typically formulated in curvilinear coordinatesusing index notation.The whole displacement field u Ω ( x ) of a material point within the shell of thickness t isdecomposed into the displacement of the midsurface u and the rotation of the normal vectorwhich is modelled with a difference vector approach w , see Figure 8. With the absence of ζ n Γ P ( x Γ , ζ )Undeformedmidsurface Γ ζζ ¯ n Γ ¯ P ( x Γ , ζ )Deformedmidsurface ¯Γ yzx x Γ x Γ u ζ wu Ω Fig. 8:
Displacement field u Ω of the KL shell. transverse shear deformations, the rotation of the normal vector is only a function of themidsurface displacement w = − h ∇ dirΓ u + ( ∇ dirΓ u ) T i · n Γ = H · u − ∇ Γ ( u · n Γ ) = − u dir ,x · n Γ u dir ,y · n Γ u dir ,z · n Γ , (3.1) .1 Kirchhoff-Love shell based on the TDC ∇ dirΓ u is the directional surface gradient, u dir ,i are partial directional surface deriva-tives w.r.t. i with i = { x, y, z } and H = ∇ covΓ n Γ is the Weingarten map. Based on that,the displacement field u Ω ( x ) takes the form u Ω = u + ζ w = u − ζ u dir ,x · n Γ u dir ,y · n Γ u dir ,z · n Γ , (3.2)and is only a function of the midsurface displacement. The parameter ζ is the coordinate inthickness direction | ζ | ≤ t / . The linear strain tensor ε Γ is defined by the symmetric part ofthe surface gradient of the displacement field, i.e., ε Γ = sym h ∇ dirΓ u Ω i . The in-plane strains ε PΓ are obtained with a double projection of ε Γ onto the tangent space of the midsurface Γ ε PΓ = P · ε Γ · P = ε PΓ , M + ζ ε PΓ , B , (3.3)with ε PΓ , M = sym [ ∇ covΓ u ] , (3.4) ε PΓ , B = − h u cov ,ij · n Γ i , (3.5)where ε PΓ , M are membrane and ε PΓ , B are bending strains, respectively. The superscript “cov”refers to covariant derivatives. The covariant gradient of u is obtained with an additionalprojection of the directional gradient onto the tangent plane, i.e., ∇ covΓ u = P · ∇ dirΓ u . Thetransverse shear strains ε SΓ ( x ) are vanishing due to the kinematic assumption from above ε SΓ = Q · ε Γ + ε Γ · Q = . (3.6)A linear elastic material governed by Hooke’s law with the modified Lamé constants µ = E ν ) , λ = Eν − ν in order to eliminate the normal stress in thickness direction shall beemployed herein. Based on that, the linear in-plane stress tensor yields σ Γ ( x ) = 2 µ ε PΓ ( x ) + λ tr[ ε PΓ ( x )] I . (3.7)With the assumption of a constant shifter in thickness direction, an analytical pre-integrationw.r.t. the thickness is possible and the stress resultants, such as moment tensor m Γ and2 Governing equations effective normal force tensor ˜ n Γ are identified as m Γ = Z t / − t / ζ P · σ Γ · P d ζ = t σ Γ ( ε PΓ , B ) = P · m dirΓ · P , (3.8)˜ n Γ = Z t / − t / P · σ Γ · P d ζ = t σ Γ ( ε PΓ , M ) = P · ˜ n dirΓ · P , (3.9)where the superscript “dir” indicates that only directional derivatives of u are used and,therefore, the physical stress resultants ( m Γ , ˜ n Γ ) can be computed without the need ofcovariant derivatives which yields a significant simplification in the implementation [97].The moment and effective normal force tensors are symmetric, in-plane tensors and theirtwo non-zero eigenvalues are the principal moments or forces, respectively. Note that inthe case of curved shells, the physical normal force tensor is n realΓ = ˜ n Γ + H · m Γ and is, ingeneral, not symmetric, but features one zero eigenvalue just as ˜ n Γ .Based on the stress resultants, the force equilibrium for the KL shell in strong form becomesdiv Γ n realΓ + n Γ div Γ q Γ + H · div Γ m Γ = − f , (3.10)where f is the load vector per area and q Γ = P · div Γ m Γ is the transverse shear force vectorobtained by equilibrium considerations. The equilibrium in strong form is a fourth-ordersurface PDE for the unknown field u . With suitable boundary conditions, the completeboundary value problem (BVP) for KL shells in the frame of the TDC is defined.For the unknown displacement field u and the rotation ω t ∂ Γ along the boundary, there existtwo non-overlapping parts at the boundary of the shell ∂ Γ. In particular, the Dirichletboundary ∂ Γ D ,i and the Neumann boundary ∂ Γ N ,i with i = { u , ω } . The correspondingboundary conditions are u = ˆ g u on ∂ Γ D , u , e p = ˆ p on ∂ Γ N , u ,ω t ∂ Γ = − [ ∇ dirΓ u T · n Γ ] · n ∂ Γ = ˆ g ω on ∂ Γ D ,ω ,m t ∂ Γ = n ∂ Γ · m Γ · n ∂ Γ = ˆ m on ∂ Γ N ,ω , (3.11)where ˆ g u are the displacements, ˆ p are the conjugated, effective forces to u , ˆ g ω is the ro-tation along the boundaries and ˆ m is the conjugated bending moment to ω t ∂ Γ at theircorresponding part of the boundary. In Figure 9, the possible boundary conditions ex-pressed in terms of the local triad ( t ∂ Γ , n ∂ Γ , n Γ ) and their conjugated, effective forces .1 Kirchhoff-Love shell based on the TDC e p t ∂ Γ , e p n ∂ Γ , e p n Γ ) and bending moment m t ∂ Γ are visualized. Furthermore, in Figure 9(b),the Kirchhoff forces at corners are shown in green. u t ∂ Γ u n ∂ Γ u n Γ u n Γ ω t ∂ Γ Γ ∂ Γ D (a) displacements and rotation e p t ∂ Γ e p n ∂ Γ e p n Γ e p n Γ m t ∂ Γ F c F c F c F c Γ ∂ Γ N (b) effective forces and bending mo-ment Fig. 9:
Boundary conditions for the KL shell in terms of the local triad ( t ∂ Γ , n ∂ Γ , n Γ ): (a) de-composition of the midsurface displacement u and the rotation ω t ∂ Γ along the boundary,(b) conjugated, effective forces, bending moment and corner forces at the boundary. The effective boundary forces e p are e p = e p t ∂ Γ t ∂ Γ + e p n ∂ Γ n ∂ Γ + e p n Γ n Γ (3.12)with e p t ∂ Γ = ( n realΓ · n ∂ Γ ) · t ∂ Γ | {z } p t ∂ Γ + ( H · t ∂ Γ ) · t ∂ Γ m n ∂ Γ , (3.13) e p n ∂ Γ = ( n realΓ · n ∂ Γ ) · n ∂ Γ | {z } p n ∂ Γ + ( H · t ∂ Γ ) · n ∂ Γ m n ∂ Γ , (3.14) e p n Γ = q Γ · n ∂ Γ | {z } p n Γ + ∇ Γ m n ∂ Γ · t ∂ Γ , (3.15)where m n ∂ Γ = ( m Γ · n ∂ Γ ) · t ∂ Γ . The Kirchhoff forces are defined as F c = n Γ m n ∂ Γ (cid:12)(cid:12)(cid:12) − C + C , (3.16)where + C and − C are points close at a corner C . Further details regarding the derivationof the effective boundary forces and the Kirchhoff forces are presented in [97].4 Governing equations
The enforcement of boundary conditions on trimmed surfaces is a non-trivial task, in par-ticular if higher-order accuracy is desired. The treatment of Neumann boundary conditionis, compared to Dirichlet (essential) boundary conditions, less challenging because theinvolved terms are already incorporated in the weak form and “only” a higher-order accu-rate integration, see Section 2.3, along the corresponding Neumann boundary remains. Thesituation for Dirichlet boundary conditions is significantly different and there are variousapproaches for their enforcement.In principal, essential boundary conditions may be enforced in a strong sense by directlyprescribing nodal values or in a weak sense by modifying the weak form with additionalconstraints. In the frame of IGA, a strong enforcement of essential boundary conditionsfor the rotation-free KL shell is only possible in special cases and would limit the generalityof the approach significantly. Therefore, we follow the same rationale as in [54, 91] andemploy a weak enforcement of the essential boundary conditions.A popular strategy for weakly enforcing essential boundary conditions is the penaltymethod [7]. The main advantages of the penalty method are the built-in linear indepen-dence of the constraints which maintains the positive definiteness of the stiffness matrix andthe straight forward implementation. On the other hand, the overall approach is variation-ally inconsistent and suffers from the interplay between the accuracy and violation of theconstraint conditions, in particular, if optimal higher-order convergence rates are desired.Furthermore, the penalty parameter has direct influence on the conditioning of the stiffnessmatrix. Another approach is the Lagrange multiplier method [6, 19, 97]. Although themethod is variationally consistent, additional degrees of freedom are introduced and thepositive definiteness of the augmented system of equations is not guaranteed. Furthermore,in the context of fictitious domain methods (FDMs) the discretization of the Lagrange mul-tiplier fields along the inner-element boundaries is not an easy task and depends on thecut scenarios. Due to these shortcomings, Nitsche’s method [85] has been developed to bea standard choice in FDMs because Nitsche’s method is variationally consistent, suitablefor higher-order and does not require the discretization of auxiliary fields. The originalapproach from [85] has been adopted to various applications for enforcing essential bound-ary conditions, see, e.g., [39, 40, 45, 55, 91, 99] and coupling, see, e.g., [3, 52, 54, 93]. Inprincipal, there are two different versions of Nitsche’s method. The symmetric version ofNitsche’s method requires an additional stabilization to ensure positive definiteness [40]. .2 Enforcement of boundary conditions and weak form L -norm of the primal variables arepossibly suboptimal [4, 17]. Nevertheless, we prefer the non-symmetric version of Nitsche’smethod because our numerical studies showed optimal behaviour even in the L -norm (inagreement to the results in [17, 93]), and the absence of an additional stabilization termcompared to the symmetric Nitsche method is very beneficial. In addition, a direct solveris employed so that the non-symmetric system of equations is not a problem.The non-symmetric version of Nitsche’s method for the KL shell formulated in curvilinearcoordinates is presented for the linear shell in [54] and for large deformations in [50]. Herein,we propose the non-symmetric version of Nitsche’s method for KL shell in the frame of theTDC in symbolic notation which results in a more compact formulation and the effectiveboundary forces are not needed explicitly.The following function spaces are introduced in order to convert the equilibrium in strongform, see Equation (3.10), to the weak form using the non-symmetric version of Nitsche’smethod for the enforcement of essential boundary condtions S := { u : Γ → R | u ∈ [ H (Γ)] : u · n Γ ∈ H (Γ) } , (3.17) V := S . (3.18)A detailed discussion regarding the function space for the KL shell with Nitsche’s methodis given in [11]. Following the usual procedure as presented in [97], the weak form of theequilibrium reads as follows: The task is to find u ∈ S such that a ( u , v ) − Z ∂ Γ D , u v · e p ( u ) d s + v · n Γ m n ∂ Γ ( u ) (cid:12)(cid:12)(cid:12) − D + D (3.19) − Z ∂ Γ D ,ω ω t ∂ Γ ( v ) m t ∂ Γ ( u ) d s = h F , v i ∀ v ∈ V , (3.20)6 Governing equations with a ( u , v ) = Z Γ ∇ dirΓ v : ˜ n Γ − ε dirΓ , B ( v ) : m Γ d A, h F , v i = Z Γ v · f d A + Z ∂ Γ N , u v · ˆ p d s + Z ∂ Γ N ,ω ω t ∂ Γ ( v ) ˆ m d s − v · ˆ F c , where D is a corner within the Dirichlet boundary for the displacements ∂ Γ D , u and ˆ F c arecorner forces at the corresponding Neumann boundary ∂ Γ N , u . The additional terms alongthe Dirichlet boundaries are caused by the fact that the test functions do not vanish there.The Nitsche terms which have to be added to the weak form in order to enforce thedisplacements ˆ g u and rotations ˆ g ω at the Dirichlet boundaries are the energy conjugatedterms for the displacements and rotation at their Dirichlet boundaries Z ∂ Γ D , u (ˆ g u − u ) · e p ( v ) d s + (ˆ g u − u ) · n Γ m n ∂ Γ ( v ) (cid:12)(cid:12)(cid:12) − D + D , (3.21) Z ∂ Γ D ,ω [ˆ g ω − ω t ∂ Γ ( u )] m t ∂ Γ ( v ) d s . (3.22)Adding these terms to the RHS of Equation (3.20) and reordering terms yields a ( u , v ) − Z ∂ Γ D , u v · e p ( u ) d s + v · n Γ m n ∂ Γ ( u ) (cid:12)(cid:12)(cid:12) − D + D + Z ∂ Γ D , u u · e p ( v ) d s − u · n Γ m n ∂ Γ ( v ) (cid:12)(cid:12)(cid:12) − D + D − Z ∂ Γ D ,ω ω t ∂ Γ ( v ) m t ∂ Γ ( u ) d s + Z ∂ Γ D ,ω ω t ∂ Γ ( u ) m t ∂ Γ ( v ) d s = h F , v i + Z ∂ Γ D , u ˆ g u · e p ( v ) d s + ˆ g u · n Γ m n ∂ Γ ( v ) (cid:12)(cid:12)(cid:12) − D + D Z ∂ Γ D ,ω ˆ g ω m t ∂ Γ ( v ) d s ∀ v ∈ V . (3.23)Applying integration by parts on the effective normal force e p n Γ we can get rid of thecorner forces. Furthermore, a shift of the rather cumbersome derivative on the drillingmoment m n ∂ Γ in e p n Γ , see Equation (3.15), simplifies the obtained terms in Equation (3.23)7significantly to a ( u , v ) − Z ∂ Γ D , u ω n ∂ Γ ( v ) m n ∂ Γ ( u ) + v · p ( u ) d s | {z } bound. terms due to v ∈ V N + Z ∂ Γ D , u ω n ∂ Γ ( u ) m n ∂ Γ ( v ) + u · p ( v ) d s | {z } Nitsche terms for displ. on LHS − Z ∂ Γ D ,ω ω t ∂ Γ ( v ) m t ∂ Γ ( u ) d s | {z } bound. terms due to v ∈ V N + Z ∂ Γ D ,ω ω t ∂ Γ ( u ) m t ∂ Γ ( v ) d s | {z } Nitsche term for rot. on LHS = h F , v i + Z ∂ Γ D , u ˆ g u · p ( v ) + ω n ∂ Γ (ˆ g u ) m n ∂ Γ ( v ) d s | {z } Nitsche terms for displ. on RHS + Z ∂ Γ D ,ω ˆ g ω m t ∂ Γ ( v ) d s | {z } Nitsche term for rot. on RHS ∀ v ∈ V , (3.24)with ω n ∂ Γ ( u ) = − [ ∇ dirΓ u T · n Γ ] · t ∂ Γ . Circumventing the corner forces and shifting thetangent derivative on the drilling moment has already been presented for the Kirchhoffplate in [93] and is, herein, straightforwardly extended to the KL shell. Compared to [54]this may lead to a more compact implementation of the Nitsche terms because the effectiveboundary forces are not needed explicitly in the implementation. Note that in case of slipsupports, e.g., membrane support, or displacement constraints in a selected, arbitrary unitdirection d with the magnitude ˆ G d , the involved terms for the displacement along ∂ Γ D , u in Equation (3.24) have to be converted accordingly. In this section, the proposed approach is tested on a set of benchmark examples consistingof the well-known Scordelis-Lo roof [9], a flat shell embedded in R [97] and a curved,clamped circular shell loaded by gravity and inhomogeneous Dirichlet boundary conditions.All examples are computed on trimmed patches. The trimming curves are defined in theparameter space by means of level-set functions, see Section 2.1. The convergence criteriaor error measurement becomes more strict from the first to last test case. This is usefulin order to confirm the higher-order accuracy of the overall approach. In particular, inthe Scordelis-Lo roof problem, the convergence criterion is a comparison between a given8 Numerical results reference displacement at a certain point. In the second example, relative L -norms in thedisplacements and stress resultants ( n realΓ , m Γ ) are investigated. In the last test case, wherethe geometry and boundary conditions are more challenging, the residual and energy errorsare computed due to the fact that an analytical solution is not available. Note that thedefinitions of the residual and energy errors are given below. In addition to the error plots,the estimated condition numbers for all test cases are shown.The employed basis functions are NURBS as defined in Section 2.1. Following the isopara-metric concept, the discrete midsurface displacement is then defined as u h = u h,i E i , with E i being Cartesian base vectors, with i = 1 , , u h,i = B T · ˆ u i . The degrees of free-dom (DOFs) ˆ u i are stored at the control points c . In the presented examples, the orderof the basis functions is varied between 3 ≤ p ≤
6. Therefore, the discrete midsurfacedisplacement u h is in the Sobolev space [ H p (Γ)] ⊂ S with p ≥
3. Consequently, theadditional continuity requirements, i.e., third-derivatives in the transverse shear forces atthe Dirchlet boundary ∂ Γ D , u , resulting from the additional Nitsche terms are easily metwhen cubic NURBS are employed. The element scale factor n , which is the number ofknot spans (elements) in one direction, is varied as 4 ≤ n ≤
64. The element scale factor n is proportional to the element size h ∼ / n . The presented numerical results are computedwith a threshold for the minimal relative support size α ∈ [0 . , .
6] in the context of theextended B-splines, see Section 2.2.1.
The first example is the Scordelis-Lo roof, which is part of the well-known shell obstaclecourse, and is taken from [9]. The whole problem including the definition of the param-eter space and its corresponding trimming curves is defined in Figure 10. As a referencedisplacement for the KL shell, we choose the converged displacement presented in [12] asreference displacement.Employing symmetry boundary conditions only one fourth of the whole geometry is con-sidered in the numerical simulation. The symmetry boundaries are modelled as trimmingcurves in the parameter space and depending on the number of knot spans n , the knotspans may be cut or aligned at the knot spans as visualized in Figure 11(a). In Fig-ure 11(b), the whole, undeformed patch is plotted in the physical space. Furthermore, thedeformed midsurface with magnified displacements by a factor 10 is added to the figure.The colors on the surface are the Euclidean norm of the displacement field u h . In this .1 Scordelis-Lo roof Geometry: Cylindrical shell L = 50 R = 25 φ = 80° t = 0 . r, s ) ∈ [0 , Trimming curves: φ ( r ) = − ( r + 1) , φ ( r ) = − ( s + 1)Material parameters: E = 4 . × ν = 0 . f = [0 , , − T Support: Rigid diaphragms at it endsRef. displacement: | u z,i, Ref | = 0 . x i = ( − R · sin( φ / ) , , R · cos( φ / )) T Fig. 10:
Definition of Scordelis-Lo roof problem. (a) trimmed parameter space (b) displacements
Fig. 11: (a) The grey domain is the trimmed parameter space A v defined by φ and φ , (b)deformed domain with scaled displacement u by a factor 10. particular example, the element scale factor n in the convergence study is varied between { , , , } and the threshold α = 0 . / n . It can be concluded that the accuracy at eachlevel increase for higher orders and the overall behaviour of the converge is in agreementwith the results presented in , e.g., [12, 68, 97]. An interesting phenomenon occurs at thecoarsest level for the orders p ≥
4, where the displacements are overestimated. This maybe due to large errors in the boundary conditions caused by a combination of rather coarsemeshes together with the elimination of unstable DoFs at the boundaries and Nitsche’s0
Numerical results method. However, this is only relevant on very coarse meshes and decreases drasticallyfor finer levels when the resolution allows a reasonable boundary discretization. Due to (a) convergence (b) condition number
Fig. 12: (a) Normalized convergence to the reference displacement u z,i, Ref = − . x i , (b) normalized, estimated condition numbers, the reference value is 8 . × ,which is the estimated condition number at n = 6 , p = 3. the lack of significant digits in the given reference solution, the error plot is not in doublelogarithmic scale. However, the style of the presentation follows those in many otherreferences, e.g., [9, 12, 26, 68, 97]. In Figure 12(a), the normalized estimated conditionnumbers of the stiffness matrix are plotted. The condition numbers are computed usingthe “condest” function of MATLAB. The condition numbers are bounded for all ordersand levels of refinement. The rather large jump between the orders 4 and 5 indicates thatthere might be potential for improvement in the choice of the threshold α for the higherorders. However, a further improvement is not mandatory since the resulting conditionnumbers are moderate. As a result of this, it can be concluded that the particular choiceof the threshold is rather flexibel.Unfortunately, this test case does not feature a sufficiently smooth solution and, therefore,optimal higher-order convergence rates cannot be achieved. Consequently, the higher-order accuracy of the proposed approach cannot be numerically confirmed with this testcase. Nevertheless, the Scordelis-Lo roof problem is a well-known and accepted benchmarkexample in the field of shell analysis. .2 Simply supported flat shell in R R The second example is taken from [97] and is a flat shell embedded in the physical space R . The problem and parameter space including trimming curves is defined in Figure 13.The simply supported shell is loaded in tangential and normal direction and in this par- Geometry: Quadrilateral flat shell L = 1 t = 0 . n Γ = [ − / , − √ / , √ / ] T Parameter space: r ∈ [0 , T − / + 0 . s ∈ [0 , T − / + 0 . φ ( r ) = s, φ ( r ) = − ( r + 1), φ ( r ) = r, φ ( r ) = − ( s + 1)Material parameters: E = 10 000 ν = 0 . f t and f n from [97]Support: Simple support on all edges Fig. 13:
Definition of flat shell problem. ticular case, an analytical solution is available. The exact solution is given in [97]. InFigure 14(a), the trimmed parameter space A v and the corresponding trimming curves arevisualized. The deformed domain is illustrated in Figure 14(b) and Figure 14(c) with scaled (a) trimmed parameter space (b) front view (c) rotated view Fig. 14: (a) The grey domain is the trimmed parameter space A v defined by φ i with i = ∈ [1 , u by two orders of magnitude, (c) rotatedview of deformed shell displacements by a factor 100. Analogously as above, the color on the surface indicatesthe Euclidean norm of the displacement field u h .2 Numerical results
The results of the convergence studies are presented in Figure 15. The element scale factor n is varied between { , , , , } and α = 0 .
4. In Figure 15(a), the relative L -norm ofthe primal variable, i.e., the midsurface displacement u h is plotted. Optimal higher-orderconvergence rates O ( p + 1) w.r.t. to the element size h are achieved which is an agreementto results shown in [17, 93]. Note that, theoretically, sub-optimal convergence rates couldhave been expected [17]. Optimal convergence rates are also achieved in the relative L -norms of the normal forces O ( p ), see Figure 15(b), and the bending moments O ( p − L -errors clearly emphasises thehigher-order accuracy of the proposed approach for the flat shell problem. In Figure 15(d),the estimated condition numbers are plotted. In this test case, the condition numbersincrease with O (4) which is expected for fourth-order PDEs. The jumps in the conditionnumbers between each order is also well-known in the context of higher-order finite elementapproaches, see, e.g., [44, 99]. The overall behaviour of the condition numbers comparedto the first test case is different which indicates that the choice of the threshold α seemsto be almost optimal (with the exception for p = 3 where it might be slightly too high).Nevertheless, a significant influence on the results in the relative L -norms has not beenobserved. .2 Simply supported flat shell in R (a) rel. L -norm displacements (b) rel. L -norm normal forces (c) rel. L -norm bending moments (d) condition number Fig. 15:
Results of convergence studies for simply supported flat shell: (a) Relative L -norm ofdisplacements u , (b) relative L -norm of internal normal forces n realΓ , (c) relative L -norm of bending moments m Γ , (d) normalized condition numbers, the reference valueis 1 . × , which is the condition number at n = 4 , p = 3. Numerical results
In the last example, a more general geometry and inhomogeneous boundary conditionsare considered. In particular, the shape of the shell is defined through a non-linear mapfrom the parameter space to the physical domain. The boundary is defined by a circulartrimming curve in the parameter space. The details of the problem definition are givenin Figure 16. For this example, an analytical solution is not available and the concept of
Geometry: Mapped circular shell x ( r ) = 2 + r − s − . rs + 0 .
75 sin(2 r + 0 . y ( r ) = 1 + r + s + 0 . rs + 0 . r + 1 . s ) z ( r ) = − . . r + 0 . s + sin( rs ) + 0 . x − Parameter space: ( r, s ) ∈ [ − . , . T Trimming curves: φ ( r ) = 0 . − k r k Material parameters: E = 10 000 ν = 0 . f = [0 , , − T Inhomogeneous essential boundary condition: ˆ g u = − . n Γ Support: Clamped edge
Fig. 16:
Definition of the clamped, circular shell problem. manufactured solutions is not practical for the defined problem. As an alternative errormeasure, we employ the concept of residual errors [97]. The residual error is the relative L -error of the equilibrium in strong form, see Equation (3.10), and is defined as ε = R A v h div Γ n realΓ + n Γ div Γ q Γ + H · div Γ m Γ + f i d A R A v f d A . (4.1)The evaluation of the strong form requires fourth-order surface derivatives of the midsurfacedisplacement u h . Consequently, the theoretical optimal order of convergence is O ( p − n is varied as { , , , } and thethreshold α = 0 . p = 3) does not convergence in the residual .3 Clamped circular shell (a) trimmed parameter space (b) displacements Fig. 17: (a) The grey domain is the trimmed parameter space A v defined by φ , (b) deformeddomain with scaled displacement u by a factor 1 . errors due to O ( p −
3) = O (0). For the order p ≥ O ( p − (a) force equilibrium (b) energy error (c) condition number Fig. 18: (a) Residual error ε rel,residual , (b) relative energy error, the reference value is5 .
857 19 × , (c) normalized, estimated condition numbers, the reference value is4 . × , which is the condition number at n = 4 , p = 3. the estimated condition numbers of the stiffness matrices in the convergence analyses areplotted and it can bee seen that the behaviour is similar to the second test case which isas expected for a fourth-order PDE. Analogously as above, the threshold α = 0 . p = 36 Conclusion may be slightly reduced in order to decrease the condition number.
The use of higher-order basis functions for shell models is no novelty, especially since theintroduction of isogeometric analysis (IGA). So it may come as a surprise that achievinghigher-order convergence rates is still the exception. This circumstance is rooted in thefact that the related requirements on the geometry and the analysis are high. Alreadyminor discrepancies in the discretization, application of boundary conditions, or the set-upof the system of equations diminish the accuracy of the results significantly.We present a Kirchhoff-Love (KL) shell formulation that treats all these numerical aspectsin a consistent higher-order fashion. From a geometric point of view, a B-spline tensor-product surface defines the geometry. It is at least C to fulfil the requirements of theKL theory. Furthermore, the trimming concept allows the specification of arbitrary shapesover this single-patch representation. Multi-patch models are not considered because thesolution across different patches is non-smooth in general, and thus, these models are in-herently low-order. It is worth noting that this argument also holds for corners in theboundary of a shell for most cases. Using a trimmed single-patch representation provides aremedy to this problem. Regarding the analysis itself, the non-symmetric Nitsche methodenforces the boundary conditions. The proposed definition of the Nitsche terms circum-vents corner forces and shifts the tangent derivative on the drilling moment, which enablesa more compact and less cumbersome implementation. The system of equations is set upby a conventional Galerkin approach. However, the treatment of cut elements, which arisedue to trimming, requires special attention with respect to conditioning and integration.We employ the extended B-spline concept and a reparameterization based on higher-orderLagrange elements to address these points. Combining these ingredients listed so far en-ables the consistent high-order accurate numerical treatment of KL shells. Finally, ourframework is completed by formulating the KL shell in the frame of tangential differentialcalculus (TDC) and using level-set functions as trimming curves. These two particularchoices are, however, not necessary for the achievement of higher-order accuracy.Numerical experiments demonstrate that the approach proposed yields optimal conver-gence rates. The Scordelis-Lo roof problem is investigated as a classic shell benchmark.Although the expected convergence behavior is shown, this example – just like other tra-7ditional shell benchmark tests – does not provide insight if higher-order accuracy has beenachieved. On the other hand, the results of a simply supported flat shell subjected to a man-ufactured solution show optimal higher-order convergence rates in the relative L -norm.This optimal behavior is not restricted to flat geometries, as confirmed by the problemof a general clamped shell with inhomogeneous boundary conditions. Note that for sucha complex problem, the concept of manufactured solutions is not practical. Hence, weemploy the residual error, i.e., the relative L -error of the equilibrium in strong form. Thismeasure requires the evaluation of fourth-order surface derivatives, which is admittedly achallenge implementation-wise. Nevertheless, it allows a reliable assessment of higher-orderaccuracy.The numerical experiments also confirm that the extended B-spline concept controls thecondition number of the system matrices. For coarse high-order discretizations, however,the involved extension operation may affect the Nitsche method. A possible solution maybe the application of extended truncated hierarchical B-splines or another strategy for theconditioning. In general, the substitution of the techniques employed for each numericalaspect by an alternative scheme that allows a consistent higher-order treatment is a futuretopic worth exploring. A Computation of extrapolation weights
The main task in the computation of the non-trivial extrapolation weights (2.10) is thedetermination of the polynomial coefficients α and β . The latter corresponds to the Newtonpolynomial ψ j,p = P pk =0 β k ξ k and can be explicitly expressed as β k = ( − p − k L X l =1 Y m ∈ T p − k,l ξ m with L = pp − k ! (A.1)where (cid:16) pm (cid:17) is the binomial coefficient defined as pm ! := p !( p − m )! m ! , (A.2)and the sum over T m,l denotes all m -combinations with repetition of the knots appearingin the definition (2.11) of ψ j,p , i.e., ξ j +1 , . . . , ξ j + p .8 Computation of extrapolation weights
The coefficients α k of the segments B si of B i,p can be obtained by Taylor expansion B si ( ξ ) = p X k =0 B ( k ) i,p (cid:16) ˜ ξ (cid:17) k ! (cid:16) ξ − ˜ ξ (cid:17) k = p X k =0 α ∗ k (cid:16) ξ − ˜ ξ (cid:17) k , ˜ ξ ∈ [ ξ s , ξ s +1 ) (A.3)and subsequent conversion to power basis form B si ( ξ ) = p X k =0 α k ξ k with α k = p X m = k mk ! α ∗ m (cid:16) − ˜ ξ (cid:17) m − k . (A.4)Alternatively, Algorithm (1) provides a derivation of the α k , which avoids the derivativesin Equation (A.3). It employs a MATLAB based syntax, and is inspired by the conversionof a B-spline to a matrix representation presented in [90].Algorithm 1: Computation of the polynomial coefficient matrix A ( k, i ) = α ik . function A = M a t r i x R e p r e s e n t a t i o n O f B s p l i n e S e g m e n t ( s, p, Ξ ) % Input : % s ... knot span % p ... polynomial degree of the spline space % Ξ ... knot vector correspond to ξ % Output : % A ... all polynomial coeffients α ik related to r A = ones (1 ,1) ; for i = 1 : p k = 1;
13 L = zeros (i , i +1) ;
14 R = zeros (i , i +1) ; for j = s -i +1 : s den = Ξ ( j +1) - Ξ ( j ) ;
17 L (k , k ) = 1. - ( Ξ ( s ) - Ξ ( j ) ) / den ;
18 L (k , k +1) = ( Ξ ( s ) - Ξ ( j ) ) / den ;
19 R (k , k ) = -( Ξ ( s +1) - Ξ ( s ) ) / den ;
20 R (k , k +1) = ( Ξ ( s +1) - Ξ ( s ) ) / den ; k = k + 1; end i = A ; A = zeros ( i +1 , i +1) ; A (1:end -1 ,:) = A (1:end -1 ,:) + A i * L ; A (2:end ,:) = A (2:end ,:) + A i * R ; end The algorithm computes the coefficients α ik for all non-zero B-splines B i,p of the knot spanand stores them in a matrix A ( k, i ) = α ik . Using A , Equation (A.4) reads B si ( r ) = (1 , r, . . . , r p ) A (: , i ) with r = ξ − ξ s ξ s +1 − ξ s . (A.5)Note that the intrinsic coordinate is now r ∈ [0 ,
1) for each knot span. Hence, the entries of A provide the α ik sought only if the knot span s is defined from 0 to 1. This condition canbe easily achieved by a linear transformation of the knots such that r s = 0 and r s +1 = 1.When these transformed knots r i are also employed for the Newton polynomials ψ j,p , thereis no alteration of the resulting extrapolation weights.Algorithm (2) computes the non-trivial extrapolation weights (2.10) for the stable B-splinesof a particular knot span s . We apply Algorithm (2) for each knot span that providesextensions. Note that r j in line 18 refer to the vector of all transformed knots r i of ψ j,p .The resulting matrix E s contains the non-trivial extrapolation weights, where the stableB-spline of the knot span s and the associated degenerate ones given by J s correspond tothe matrix’s rows and columns, respectively. Finally, assembling all E s according to theglobal degrees of freedom, togehter with the trivial extrapolations weights (2.8) and (2.9),yields the overall extension matrix E .Algorithm 2: Extrapolation weights for a particular knot span function E s = C o m p u t e E x t r a p o l a t i o n W e i g h t s ( s, J s , p, Ξ ) % Input : % s ... knot span providing the extensions % J s ... index set of all degenerate B - splines related to s % p ... polynomial degree of the spline space % Ξ ... knot vector corresponding to ξ % Output : % E s ... extrapolation weights for the stable B - splines of s % polynomial coefficients of the extensions A = M a t r i x R e p r e s e n t a t i o n O f B s p l i n e S e g m e n t ( s, p, Ξ ) % knot values mapped to intrinsic coordinate r Ξ r = ( Ξ - Ξ ( s ) ) ./ ( Ξ ( s +1) - Ξ ( s ) ) ; % coefficients of the Newton polynomial based on Ξ r REFERENCES
16 B = zeros ( p +1 , length ( J s ) ) ; for n = 1 : length ( J s ) r j = Ξ r ( J s ( n ) +1 : J s ( n ) + p ) ; for k = 0 : p -1 L = nchoosek ( p , p -k ) ; if L == 1 T = r j ; else T = nchoosek ( r j , p -k ) ; end tmp = 0; for l = 1 : L tmp = tmp + prod ( T (l ,:) ) ; end
32 B ( k +1 , n ) = ( -1) ^( p -k ) * tmp ; end end % extrapolation weights E s = zeros ( size ( A ,2) , length ( J s ) ) ; kFactorial = [1 cumprod (1: size ( A ,1) -1) ]; % pre - computation of k ! for j = 1 : size ( E s ,2) for i = 1 : size ( E s ,1) for k = 0 : p E s (i , j ) = E s (i , j ) + ( -1) ^ k * ... kFactorial (end-k ) * B (end-k , j ) * ... kFactorial ( k +1) * A ( k +1 , i ) ; end end end E s = E s ./ kFactorial (end) ; References [1] Abedian, A.; Parvizian, J.; Düster, A.; Khademyzadeh, H.; Rank, E.: Performanceof different integration schemes in facing discontinuities in the Finite Cell Method.
EFERENCES Int. J. Comput. Methods , , 1–24, 2013.[2] Antolin, P.; Buffa, A.; Martinelli, M.: Isogeometric Analysis on V-reps: First results. Comp. Methods Appl. Mech. Engrg. , , 976–1002, 2019.[3] Apostolatos, A.; Schmidt, R.; Wüchner, R.; Bletzinger, K. U.: A Nitsche-type for-mulation and comparison of the most common domain decomposition methods inisogeometric analysis. Internat. J. Numer. Methods Engrg. , , 473–504, 2013.[4] Arnold, D. N.; Brezzi, F.; Cockburn, B.; Marini, L. D.: Unified analysis of discontin-uous Galerkin methods for elliptic problems. SIAM J. Numer. Anal. , , 1749–1779,2002.[5] Auricchio, F.; Calabrò, F.; Hughes, T. J. R.; Reali, A.; Sangalli, G.: A simple algo-rithm for obtaining nearly optimal quadrature rules for NURBS-based isogeometricanalysis. Comp. Methods Appl. Mech. Engrg. , , 15–27, 2012. Higher OrderFinite Element and Isogeometric Methods.[6] Babuška, I.: The Finite Element Method with Lagrangian Multipliers. Numer.Math. , , 179–192, 1973.[7] Babuška, I.: The Finite Element Method with penalty. Math. Comput. , , 221–228,1973.[8] Başar, Y.; Krätzig, W. B.: Mechanik der Flächentragwerke . Vieweg+Teubner Verlag,Braunschweig, 1985.[9] Belytschko, T.; Stolarski, H.; Liu, W. K.; Carpenter, N.; Ong, J. S. J.: Stressprojection for membrane and shear locking in shell finite elements.
Comp. MethodsAppl. Mech. Engrg. , , 221–258, 1985.[10] Benson, D. J.; Bazilevs, Y.; Hsu, M. C.; Hughes, T. J. R.: Isogeometric shell analysis:The Reissner-Mindlin shell. Comp. Methods Appl. Mech. Engrg. , , 276–289, 2010.[11] Benzaken, J.; Evans, J. A.; McCormick, S.; Tamstorf, R.: Nitsche’s method for linearKirchhoff–Love shells: Formulation, error analysis, and verification. Comp. MethodsAppl. Mech. Engrg. , , 113544, 2021.2 REFERENCES [12] Bieber, S.; Oesterle, B.; Ramm, E.; Bischoff, M.: A variational method to avoidlocking–independent of the discretization scheme.
Internat. J. Numer. Methods En-grg. , , 801–827, 2018.[13] Bischoff, M.; Ramm, E.; Irslinger, J.: Models and Finite Elements for Thin-WalledStructures. In Encyclopedia of Computational Mechanics (Second Edition) . (Stein,E.; Borst, R.; Hughes, T. J.; Hughes, T. J., Eds.), John Wiley & Sons, 2017.[14] Breitenberger, M.:
CAD-integrated Design and Analysis of Shell Structures . PhDthesis, Technische Universität München, 2016.[15] Breitenberger, M.; Apostolatos, A.; Philipp, B.; Wüchner, R.; Bletzinger, K. U.:Analysis in computer aided design: Nonlinear isogeometric B-Rep analysis of shellstructures.
Comp. Methods Appl. Mech. Engrg. , , 401–457, 2015. IsogeometricAnalysis Special Issue.[16] Burman, E.: Ghost penalty. Comptes Rendus Mathematique , , 1217–1220, 2010.[17] Burman, E.: A penalty-free nonsymmetric Nitsche-type method for the weak impo-sition of boundary conditions. SIAM J. Numer. Anal. , , 1959–1981, 2012.[18] Burman, E.; Claus, S.; Hansbo, P.; Larson, M. G.; Massing, A.: CutFEM: Discretiz-ing geometry and partial differential equations. Internat. J. Numer. Methods Engrg. , , 472–501, 2015.[19] Burman, E.; Hansbo, P.: Fictitious domain finite element methods using cut ele-ments: I. A stabilized Lagrange multiplier method. Comp. Methods Appl. Mech.Engrg. , , 2680–2686, 2010.[20] Burman, E.; Hansbo, P.: Fictitious domain finite element methods using cut ele-ments: II. A stabilized Nitsche method. Appl. Numer. Math. , , 328–341, 2012.[21] Cenanovic, M.; Hansbo, P.; Larson, M. G.: Cut finite element modeling of linearmembranes. Comp. Methods Appl. Mech. Engrg. , , 98–111, 2016.[22] Chen, Q.; Babuška, I.: Approximate optimal points for polynomial interpolation ofreal functions in an interval and in a triangle. Comp. Methods Appl. Mech. Engrg. , , 405–417, 1995. EFERENCES
Comp. Methods Appl. Mech. Engrg. , , 89–94,1996.[24] Chen, X. D.; Yong, J. H.; Wang, G.; Paul, J. C.; Xu, G.: Computing the minimumdistance between a point and a NURBS curve. Computer-Aided Des. , , 1051–1054,2008.[25] Cheng, K. W.; Fries, T. P.: Higher-order XFEM for curved strong and weak discon-tinuities. Internat. J. Numer. Methods Engrg. , , 564–590, 2010.[26] Cirak, F.; Ortiz, M.; Schröder, P.: Subdivision surfaces: A new paradigm for thin-shell finite-element analysis. Internat. J. Numer. Methods Engrg. , , 2039–2072,2000.[27] Coradello, L.; Antolin, P.; Vázquez, R.; Buffa, A.: Adaptive isogeometric analysis ontwo-dimensional trimmed domains based on a hierarchical approach. Comp. MethodsAppl. Mech. Engrg. , , 2019.[28] Cottrell, J. A.; Reali, A.; Bazilevs, Y.; Hughes, T. J. R.: Isogeometric analysis ofstructural vibrations. Comp. Methods Appl. Mech. Engrg. , , 5257–5296, 2006.[29] Dauge, M.; Düster, A.; Rank, E.: Theoretical and Numerical Investigation of theFinite Cell Method. J. Sci. Comput. , , 1039–1064, 2015.[30] de Boor, C.: A practical guide to splines , Vol. 27,
Applied Mathematical Sciences .Springer, New York, 2001.[31] de Boor, C.; Fix, G. J.: Spline approximation by quasiinterpolants.
J. Approx.Theory , , 19–45, 1973.[32] de Prenter, F.; Lehrenfeld, C.; Massing, A.: A note on the stability parameter inNitsche’s method for unfitted boundary value problems. Comput Math Appl , ,4322–4336, 2018.[33] de Prenter, F.; Verhoosel, C. V.; van Brummelen, E. H.; Evans, J. A.; Messe, C.;Benzaken, J.; Maute, K.: Multigrid solvers for immersed finite element methods andimmersed isogeometric analysis. Comput. Mech. , , 807–838, 2020.4 REFERENCES [34] de Prenter, F.; Verhoosel, C. V.; van Zwieten, G. J.; van Brummelen, E. H.: Condi-tion number analysis and preconditioning of the finite cell method.
Comp. MethodsAppl. Mech. Engrg. , , 297–327, 2017. Special Issue on Isogeometric Analysis:Progress and Challenges.[35] Dornisch, W.; Klinkel, S.: Treatment of Reissner-Mindlin shells with kinks withoutthe need for drilling rotation stabilization in an isogeometric framework. Comp.Methods Appl. Mech. Engrg. , , 35–66, 2014.[36] Dornisch, W.; Klinkel, S.; Simeon, B.: Isogeometric Reissner-Mindlin shell analysiswith exactly calculated director vectors. Comp. Methods Appl. Mech. Engrg. , ,491–504, 2013.[37] Dréau, K.; Chevaugeon, N.; Moës, N.: Studied X-FEM enrichment to handle materialinterfaces with higher order finite element. Comp. Methods Appl. Mech. Engrg. , ,1922–1936, 2010.[38] Echter, R.; Oesterle, B.; Bischoff, M.: A hierarchic family of isogeometric shell finiteelements. Comp. Methods Appl. Mech. Engrg. , , 170–180, 2013.[39] Embar, A.; Dolbow, J.; Harari, I.: Imposing Dirichlet boundary conditions withNitsche’s method and spline-based finite elements. Internat. J. Numer. MethodsEngrg. , , 877–898, 2010.[40] Fernández-Méndez, S.; Huerta, A.: Imposing essential boundary conditions in mesh-free methods. Comp. Methods Appl. Mech. Engrg. , , 1257–1275, 2004.[41] Fries, T. P.: Higher-order accurate integration for cut elements with Chen-Babuškanodes. In Advances in Discretization Methods: Discontinuities, virtual elements, fic-titious domain methods . (Ventura, G.; Benvenuti, E., Eds.), Vol. 12,
SEMA SIMAISpringer Series , Springer, Berlin, 245–269, 2016.[42] Fries, T. P.; Omerović, S.: Higher-order accurate integration of implicit geometries.
Internat. J. Numer. Methods Engrg. , , 323–371, 2016.[43] Fries, T. P.; Omerović, S.; Schöllhammer, D.; Steidl, J.: Higher-order meshing ofimplicit geometries—part I: Integration and interpolation in cut elements. Comp.Methods Appl. Mech. Engrg. , , 759–784, 2017. EFERENCES
Comp. Methods Appl. Mech. Engrg. , , 270–297,2017.[45] Fries, T. P.; Schöllhammer, D.: A unified finite strain theory for membranes andropes. Comp. Methods Appl. Mech. Engrg. , , 113031, 2020.[46] Gfrerer, M. H.: A C -continuous Trace-Finite-Cell-Method for linear thin shell anal-ysis on implicitly defined surfaces. arXiv e-prints , 2020. arXiv: 2007.00075.[47] Gockenbach, M. S.: Understanding and Implementing the Finite Element Method .SIAM, Philadelphia, 2006.[48] Gordon, W. J.; Hall, C. A.: Construction of curvi-linear co-ordinate systems andapplications to mesh generation.
Internat. J. Numer. Methods Engrg. , , 461–477,1973.[49] Gordon, W. J.; Hall, C. A.: Transfinite element methods: blending function inter-polation over arbitrary curved element domains. Numer. Math. , , 109–129, 1973.[50] Guo, Y.; Do, H.; Ruess, M.: Isogeometric stability analysis of thin shells: Fromsimple geometries to engineering models. Internat. J. Numer. Methods Engrg. , ,433–458, 2019.[51] Guo, Y.; Heller, J.; Hughes, T. J. R.; Ruess, M.; Schillinger, D.: Variationallyconsistent isogeometric analysis of trimmed thin shells at finite deformations, basedon the STEP exchange format. Comp. Methods Appl. Mech. Engrg. , , 39–79,2018.[52] Guo, Y.; Ruess, M.: Nitsche’s method for a coupling of isogeometric thin shells andblended shell structures. Comp. Methods Appl. Mech. Engrg. , , 881–905, 2015.[53] Guo, Y.; Ruess, M.: Weak Dirichlet boundary conditions for trimmed thin isogeo-metric shells. Comput Math Appl , , 1425–1440, 2015.[54] Guo, Y.; Ruess, M.; Schillinger, D.: A parameter-free variational coupling approachfor trimmed isogeometric thin shells. Comput. Mech. , , 693–715, 2017.6 REFERENCES [55] Hansbo, A.; Hansbo, P.: An unfitted finite element method, based on Nitsche’smethod, for elliptic interface problems.
Comp. Methods Appl. Mech. Engrg. , ,5537 – 5552, 2002.[56] Hiemstra, R. R.; Calabrò, F.; Schillinger, D.; Hughes, T. J. R.: Optimal and reducedquadrature rules for tensor product and hierarchically refined splines in isogeometricanalysis. Comp. Methods Appl. Mech. Engrg. , , 966–1004, 2017. Special Issue onIsogeometric Analysis: Progress and Challenges.[57] Hoffmann, C. M.: Geometric and Solid Modelling: An Introduction . Morgan Kauf-mann, San Francisco, CA, 1989.[58] Hoffmann, C. M.: Implicit curves and surfaces in CAGD.
IEEE Computer Graphicsand Applications , , 79–88, 1993.[59] Höllig, K.: Finite Element Methods with B-Splines , Vol. 26,
Frontiers in AppliedMathematics . SIAM, 2003.[60] Hughes, T. J. R.; Cottrell, J. A.; Bazilevs, Y.: Isogeometric analysis: CAD, finiteelements, NURBS, exact geometry and mesh refinement.
Comp. Methods Appl. Mech.Engrg. , , 4135–4195, 2005.[61] Hughes, T. J. R.; Reali, A.; Sangalli, G.: Efficient quadrature for NURBS-basedisogeometric analysis. Comp. Methods Appl. Mech. Engrg. , , 301–313, 2010.[62] Höllig, K.; Reif, U.: Nonuniform web-splines. Comput. Aided Geom. Design , ,277–294, 2003.[63] Johannessen, K. A.: Optimal quadrature for univariate and tensor product splines. Comp. Methods Appl. Mech. Engrg. , , 84–99, 2017. Special Issue on IsogeometricAnalysis: Progress and Challenges.[64] Joulaian, M.; Hubrich, S.; Düster, A.: Numerical integration of discontinuities onarbitrary domains based on moment fitting. Comput. Mech. , , 979–999, 2016.[65] Kang, P.; Youn, S. K.: Isogeometric analysis of topologically complex shell structures. Finite Elem. Anal. Des. , , 68–81, 2015.[66] Kang, P.; Youn, S. K.: Isogeometric shape optimization of trimmed shell structures. Struct. Multidiscip. Optim. , , 825–845, 2016. EFERENCES
Finite Elem. Anal. Des. , , 18–40, 2016.[68] Kiendl, J.; Bletzinger, K. U.; Linhard, J.; Wüchner, R.: Isogeometric shell analysiswith Kirchhoff-Love elements. Comp. Methods Appl. Mech. Engrg. , , 3902–3914,2009.[69] Kiendl, J.; Marino, E.; De Lorenzis, L.: Isogeometric collocation for the Reissner-Mindlin shell problem. Comp. Methods Appl. Mech. Engrg. , , 645–665, 2017.[70] Kudela, L.; Zander, N.; Bog, T.; Kollmannsberger, S.; Rank, E.: Efficient and accu-rate numerical quadrature for immersed boundary methods. Advanced Modeling andSimulation in Engineering Sciences , , 1–22, 2015.[71] Kudela, L.; Zander, N.; Kollmannsberger, S.; Rank, E.: Smart octrees: Accuratelyintegrating discontinuous functions in 3D. Comp. Methods Appl. Mech. Engrg. , ,406–426, 2016.[72] Laborde, P.; Pommier, J.; Renard, Y.; Salaün, M.: High-order extended finite ele-ment method for cracked domains. Internat. J. Numer. Methods Engrg. , , 354–381,2005.[73] Legay, A.; Wang, H. W.; Belytschko, T.: Strong and weak arbitrary discontinuitiesin spectral finite elements. Internat. J. Numer. Methods Engrg. , , 991–1008, 2005.[74] Legrain, G.; Chevaugeon, N.; Dréau, K.: High order X-FEM and levelsets for complexmicrostructures: Uncoupling geometry and approximation. Comp. Methods Appl.Mech. Engrg. , , 172–189, 2012.[75] Lipton, S.; Evans, J. A.; Bazilevs, Y.; Elguedj, T.; Hughes, T. J. R.: Robustness ofisogeometric structural discretizations under severe mesh distortion. Comp. MethodsAppl. Mech. Engrg. , , 357 – 373, 2010. Computational Geometry and Analysis.[76] Marussig, B.; Hiemstra, R.; Hughes, T. J. R.: Improved conditioning of isogeometricanalysis matrices for trimmed geometries. Comp. Methods Appl. Mech. Engrg. , ,79–110, 2018.[77] Marussig, B.; Hughes, T. J. R.: A Review of Trimming in Isogeometric Analysis:Challenges, Data Exchange and Simulation Aspects. Archive Comp. Mech. Engrg. , , 1059–1127, 2018.8 REFERENCES [78] Marussig, B.; Zechner, J.; Beer, G.; Fries, T.-P.: Stable Isogeometric Analysis ofTrimmed Geometries.
Comp. Methods Appl. Mech. Engrg. , , 497–521, 2016.[79] Maurin, F.; Greco, F.; Coox, L.; Vandepitte, D.; Desmet, W.: Isogeometric colloca-tion for Kirchhoff–Love plates and shells. Comp. Methods Appl. Mech. Engrg. , ,396–420, 2018.[80] Mößner, B.; Reif, U.: Stability of tensor product B-splines on domains. J. Approx.Theory , , 1–19, 2008.[81] Moumnassi, M.; Belouettar, S.; Béchet, É.; Bordas, S. P. A.; Quoirin, D.; Potier-Ferry, M.: Finite element analysis on implicitly defined domains: An accurate rep-resentation based on arbitrary parametric surfaces. Comp. Methods Appl. Mech.Engrg. , , 774–796, 2011.[82] Mousavi, S. E.; Sukumar, N.: Numerical integration of polynomials and discontin-uous functions on irregular convex polygons and polyhedrons. Comput. Mech. , ,535–554, 2011.[83] Müller, B.; Kummer, F.; Oberlack, M.: Highly accurate surface and volume integra-tion on implicit domains by means of moment-fitting. Internat. J. Numer. MethodsEngrg. , , 512–528, 2013.[84] Naghdi, P.: The Theory of Shells and Plates. Truesdell C. (eds) Linear Theories ofElasticity and Thermoelasticity , Springer, Berlin, Heidelberg, 425–640, 1973.[85] Nitsche, J.: Über ein Variationsprinzip zur Lösung von Dirichlet-Problemen bei Ver-wendung von Teilräumen, die keinen Randbedingungen unterworfen sind.
Abhand-lungen aus dem Mathematischen Seminar der Universität Hamburg , , 9–15, 1971.[86] Oesterle, B.; Ramm, E.; Bischoff, M.: A shear deformable, rotation-free isogeometricshell formulation. Comp. Methods Appl. Mech. Engrg. , , 235–255, 2016.[87] Oesterle, B.; Sachse, R.; Ramm, E.; Bischoff, M.: Hierarchic isogeometric largerotation shell elements including linearized transverse shear parametrization. Comp.Methods Appl. Mech. Engrg. , , 383–405, 2017.[88] Omerović, S.; Fries, T. P.: Conformal higher-order remeshing schemes for implicitlydefined interface problems. Internat. J. Numer. Methods Engrg. , , 763–789, 2017. EFERENCES
The NURBS book . Springer-Verlag New York, Inc., New York,NY, USA, 2 edition, 1997.[90] Qin, K.: General matrix representations for B-splines.
The Visual Computer , ,177–186, 2000.[91] Ruess, M.; Schillinger, D.; Bazilevs, Y.; Varduhn, V.; Rank, E.: Weakly enforcedessential boundary conditions for NURBS-embedded and trimmed NURBS geome-tries on the basis of the finite cell method. Internat. J. Numer. Methods Engrg. , ,811–846, 2013.[92] Sala-Lardies, E.; Fernández-Méndez, S.; Huerta, A.: Optimally convergent high-order X-FEM for problems with voids and inclusion. Proceedings of the ECCOMAS2012 , Vienna, Austria, 2012.[93] Schillinger, D.; Harari, I.; Hsu, M. C.; Kamensky, D.; Stoter, S. K. F.; Yu, Y.; Zhao,Y.: The non-symmetric Nitsche method for the parameter-free imposition of weakboundary and coupling conditions in immersed finite elements.
Comp. Methods Appl.Mech. Engrg. , , 625–652, 2016.[94] Schillinger, D.; Hossain, S. J.; Hughes, T. J. R.: Reduced Bézier element quadraturerules for quadratic and cubic splines in isogeometric analysis. Comp. Methods Appl.Mech. Engrg. , , 1–45, 2014.[95] Schillinger, D.; Ruess, M.: The Finite Cell Method: A Review in the Contextof Higher-Order Structural Analysis of CAD and Image-Based Geometric Models. Archive Comp. Mech. Engrg. , , 391–455, 2015.[96] Schmidt, R.; Wüchner, R.; Bletzinger, K. U.: Isogeometric analysis of trimmedNURBS geometries. Comp. Methods Appl. Mech. Engrg. , , 93–111, 2012.[97] Schöllhammer, D.; Fries, T. P.: Kirchhoff-Love shell theory based on tangentialdifferential calculus. Comput. Mech. , , 113–131, 2019.[98] Schöllhammer, D.; Fries, T. P.: A unified approach for shell analysis on explicitlyand implicitly defined surfaces. Form and Force: Proceedings of the IASS Symposium2019 - Structural Membranes 2019, C. Lázaro, K.U. Bletzinger, E. Oñate (eds.) ,Barcelona, 750–757, 2019.0
REFERENCES [99] Schöllhammer, D.; Fries, T. P.: A higher-order Trace finite element method for shells.
Internat. J. Numer. Methods Engrg. , 2020.[100] Scott, R. L.:
Finite Element Techniques for Curved Boundaries . PhD-Dissertation,Massachusetts Institute of Technology, Cambridge, MA, USA, 1973.[101] Sederberg, T. W.; Anderson, D. C.; Goldman, R. N.: Implicit representation ofparametric curves and surfaces.
Computer Vision, Graphics and Image Processing , , 72–84, 1984.[102] Sederberg, T. W.; Finnigan, G. T.; Li, X.; Lin, H.; Ipson, H.: Watertight trimmedNURBS. ACM Transactions on Graphics , , 1–8, 2008.[103] Selimovic, I.: Improved algorithms for the projection of points on NURBS curvesand surfaces. Comput. Aided Geom. Design , , 439–445, 2006.[104] Šolín, P.; Segeth, K.; Doležel, I.: Higher-order finite element methods . CRC Press,Boca Raton, FL, 2003.[105] Stanford, J. W.; Fries, T. P.: A higher-order conformal decomposition finite elementmethod for plane B-rep geometries.
Computers & Structures , , 15–27, 2019.[106] Stavrev, A.; Nguyen, L. H.; Shen, R.; Varduhn, V.; Behr, M.; Elgeti, S.; Schillinger,D.: Geometrically accurate, efficient, and flexible quadrature techniques for the tetra-hedral finite cell method. Comp. Methods Appl. Mech. Engrg. , , 646–673, 2016.[107] Stazi, F. L.; Budyn, E.; Chessa, J.; Belytschko, T.: An extended finite elementmethod with higher-order elements for curved cracks. Comput. Mech. , , 38–48,2003.[108] Szabó, B.; Düster, A.; Rank, E.: The p-Version of the Finite Element Method ,Chapter 5, 119–139. John Wiley & Sons, Chichester, 2004.[109] Ventura, G.: On the elimination of quadrature subcells for discontinuous functionsin the eXtended finite element method.
Internat. J. Numer. Methods Engrg. , ,761–795, 2006.[110] Ventura, G.; Gracie, R.; Belytschko, T.: Fast integration and weight function blend-ing in the extended finite element method. Internat. J. Numer. Methods Engrg. , ,1–29, 2009. EFERENCES
Comp. Methods Appl. Mech. Engrg. , , 1–15,2013.[112] Wu, B.; Pagani, A.; Chen, W. Q.; Carrera, E.: Geometrically nonlinear refinedshell theories by Carrera Unified Formulation. Mechanics of Advanced Materials andStructures , 1–21, 2019.[113] Zi, G.; Belytschko, T.: New crack-tip elements for XFEM and applications to cohe-sive cracks.
Internat. J. Numer. Methods Engrg. , , 2221–2240, 2003.[114] Zou, Z.; Scott, M. A.; Miao, D.; Bischoff, M.; Oesterle, B.; Dornisch, W.: Anisogeometric Reissner–Mindlin shell element based on Bézier dual basis functions:Overcoming locking and improved coarse mesh accuracy. Comp. Methods Appl. Mech.Engrg. ,370