A C 1 -continuous Trace-Finite-Cell-Method for linear thin shell analysis on implicitly defined surfaces
CComputational Mechanics manuscript No. (will be inserted by the editor) A C -continuous Trace-Finite-Cell-Method for linearthin shell analysis on implicitly defined surfaces Michael H. Gfrerer
Received: date / Accepted: date
Abstract
A Trace-Finite-Cell-Method for the numerical analysis of thin shellsis presented combining concepts of the TraceFEM and the Finite-Cell-Method.As an underlying shell model we use the Koiter model, which we re-derive instrong form based on first principles of continuum mechanics by recasting well-known relations formulated in local coordinates to a formulation independentof a parametrization. The field approximation is constructed by restrictingshape functions defined on a structured background grid on the shell surface.As shape functions we use on a background grid the tensor product of cubicsplines. This yields C -continuous approximation spaces, which are required bythe governing equations of fourth order. The parametrization-free formulationallows a natural implementation of the proposed method and manufacturedsolutions on arbitrary geometries for code verification. Thus, the implementa-tion is verified by a convergence analysis where the error is computed with anexact manufactured solution. Furthermore, benchmark tests are investigated. Keywords finite element method · implicit geometry · Koiter shell · Finite-Cell-Method · TraceFEM
Due to the superior load-carrying capabilities, the mechanical analysis of shellsis of great interest in engineering. A large literature body exists on the for-mulation of shell models. We refer to [32,39,9] and references therein for anoverview. In the present paper, we consider the classical Koiter model [30,
Michael H. GfrererInstitute of Applied Mechanics, Graz University of Technology, Technikerstrasse 4, 8010Graz, AustriaTel.: +43 (316) 873 - 7149Fax: +43 (316) 873 - 7641E-mail: [email protected] a r X i v : . [ c s . C E ] J un Michael H. Gfrerer i.e. represented by a collection of local parametrizations. In contrastto these representations, we consider the case where the mid-surface is repre-sented implicitly as the zero-level set of a scalar function φ ( x ), see Figure 1for the illustration of some examples. We refer to the review article [21] for anoverview of finite element methods for problems on such surfaces. In the classi-cal surface finite element method the discretization of the unknown field relieson the higher order or exact meshing (local parametrization) of the surface[19,26]. In contrast to this, in the proposed method the discretization of thedisplacement field does not rely on parametrizations. Therefore, we provide athroughout derivation of the governing equations based on first principles ofcontinuum mechanics independent of a parametrization. This allows a natu-ral implementation of the method and also the construction of manufacturedsolutions for code verification on arbitrary geometries. Equivalent derivationsrelying on a parametrization can be found in e.g. [4] and [42]. We remark thatmembrane and thin shell formulations without relying on a parametrizationwith a mathematical focus can be found in [27,18]. For a treatment from anengineering perspective we refer to [37,45]. (a) sphere (b) torus (c) cylinder (d) gyroid Fig. 1
Examples of implicitly defined surfaces. The surfaces are defined by the level-setfunctions (a) φ = x + y + z − r , (b) φ = ( x + y + z + R − r ) − R ( x + y ), (c) φ = x + y − r (d) φ = sin x cos y + sin y cos z + sin z cos x One of the main difficulties in developing finite element methods for thinshells is the construction of C -continuous approximation spaces. For generalunstructured meshes it is not possible to ensure C -continuity with only lo-cal polynomial shape functions and the nodal degrees of freedom consist ofdisplacements and slopes only [46]. However, different non-standard triangu-lar for developed thin plate bending are the Argyris element [3,20], the Bellelement [6] or the Clough-Tocher macrotriangle [17]. A further possibility to race-Finite-Cell-Method for thin shell analysis 3 construct C -continuous approximation spaces on general space triangulationrelies on sophisticated techniques from subdivision surfaces [16]. However, ona structured quadrilateral mesh the Bogner-Fox-Schmit element [11] is a sim-ple conforming element. The constraint of a structured quadrilateral meshcan be partially overcome by introducing a smooth mapping of the geom-etry [24]. This idea can be realized in an isoparametric way by the use ofsplines for the geometry mapping and for the discretization of the displacementfield [29]. The general difficulty of constructing C -continuous approximationspaces led to approaches where the C -continuity requirement is circumvented.Among them we mention discrete Kirchhoff elements [5,2] where the Kirch-hoff constraint is enforced only at discrete points, the use of shear-deformable(Reissner-Mindlin) shell theory, were only C -continuity approximation spacesare required, mixed methods [40,33], continuous/discontinuous Galerkin meth-ods [23,28] and others.In the present paper, we combine ideas from unfitted finite element meth-ods and the Bogner-Fox-Schmit element. Following the idea of the TraceFEM[35,34] (see also CutFEM [12,13]) the approximation of the displacement fieldis constructed by restricting shape functions defined on a background mesh onthe shell surface. In particular, we follow the idea of the Finite-Cell-Method[38,44] and use a structured grid where the simple tensor product of threeuni-variant cubic spline shape functions leads to a C -continuous approxima-tion in 3D (like the Bogner-Fox-Schmit element in 2D). Therefore, the shapefunctions for approximation of the displacement field on the shell mid-surfaceare C -continuous. We remark that cut Bogner-Fox-Schmit elements for thinplates were proposed and analyzed in [1]. Therefore, the proposed method canbe seen as an extension of the work [1] from plates to curved shells.One challenge in unfitted finite element methods is the efficient integrationon the problem domain [36]. During the preparation of the present paper itturned out that due to a gowning number of constraints for finer meshes thestrategy developed in [25] is not applicable in the present situation. Therefore,we have implemented the quadrature strategy developed in [43].The implementation of the proposed method is verified by a convergenceanalysis where the error is computed with an exact manufactured solution.Furthermore, the capabilities of the method are shown in two standard andone non-standard benchmark tests. The underlying assumption in shell analysis is that the computational domainhas a small extension with respect to one coordinate. Thus, we assume that itis located around a two-dimensional mid-surface Ω which is embedded in R .In the present paper we assume that the mid-surface is defined implicitly asthe zero-level set of a function φ : R → R inside a cuboid B ⊂ R , Ω = { x ∈ B | φ ( x ) = 0 } . (1) Michael H. Gfrerer
The boundary of Ω is denoted Γ , the surface normal vector is denoted by ννν ,and the normal vector tangential to the surface on a boundary point is denotedby µµµ , see Figure 2. We assume that Ω is regular such that Fig. 2
Illustration and notation of the geometric setting ∇ φ ( x ) (cid:54) = , (2)holds in the neighborhood of Ω , where ∇ denotes the usual gradient of somescalar-valued function f : R → R , ∇ f ( x ) = f ,i e i := ∂f ( x ) ∂x i e i (3)with the Cartesian coordinates x = ( x , x , x ) and the standard Cartesianorthonormal basis { e , e , e } . Here, and in the following, the Einstein sum-mation convention applies. Whenever an index occurs once in an upper po-sition and in a lower position we sum over this index, where Latin indices i, j, . . . take the values 1 , , α, β, . . . take the values1 ,
2. Let T be some tensor space of the form R ⊗ · · · ⊗ R . In the followingwe also use the generalization of the gradient for scalar-valued functions (3)to tensor-valued functions T : R → T , ∇ T ( x ) = T ,i ⊗ e i . (4)2.1 Differential geometry of implicitly defined surfacesGiven a implicit representation (1) of the surface Ω , we can compute the unitnormal vector to the surface by ννν ( x ) = ∇ φ ( x ) ||∇ φ ( x ) || , (5)and are able to define the tangential projector, P = I − ννν ⊗ ννν. (6) race-Finite-Cell-Method for thin shell analysis 5 Furthermore, the extended Weingarten map is given by H = −∇ ννν · P = − P · ∇∇ φ · P ||∇ φ || , (7)and the mean curvature H is defined as H = tr( H ) = H : P . (8)2.2 Differential geometry of parametrized surfacesWe briefly review the differential geometry of parametrized surfaces. For de-tails we refer to e.g. [14]. Although the mid-surface Ω is assumed to be givenimplicitly, at least a local parametric representation is guaranteed to existby the implicit function theorem. This justifies to consider parametrizationsˆ g : U ⊂ R → Ω with the parameter domain U for theoretical considerations.Given the parametrization ˆ g ( θ , θ ), we can define the two covariant base vec-tors ˆ g α := ∂ ˆ g ∂θ α , which span the tangent plane to Ω . With the base vectors wecan define the unit normal vectorˆ ννν ( θ , θ ) = ˆ g ( θ , θ ) × ˆ g ( θ , θ ) || ˆ g ( θ , θ ) × ˆ g ( θ , θ ) || , (9)and the covariant coefficients of the metric ˆ G αβ = ˆ g α · ˆ g β . The contravariantcoefficients of the metric are given by [ ˆ G αβ ] = [ ˆ G αβ ] − , where [ ˆ G αβ ] is thecoefficient matrix. The contravariant base vectors can then be computed byˆ g α = ˆ G αβ ˆ g β . The covariant coefficients of the Weingarten map ˆ H = ˆ h αβ ˆ g α ⊗ ˆ g β are given by ˆ h αβ = − ˆ g α · ˆ ννν ,β , (10)and obey the symmetry relation ˆ h αβ = ˆ h βα . The mean curvature ˆ H is givenby ˆ H = ˆ h αα = ˆ h αβ ˆ G βα . (11)Furthermore, the derivatives of the base vectors are given byˆ g α,β = Γ γαβ ˆ g γ + ˆ h αβ ˆ ννν, ˆ g α,β = − Γ αβγ ˆ g γ + ˆ h αβ ˆ ννν, (12)with the surface Christoffel symbols of the second kind defined byˆ Γ γαβ = ˆ g γ · ˆ g α,β . (13) Remark:
In our notion a hat over a quantity refers to a dependency on theparametric coordinates ( θ , θ ) ∈ U , whereas no hat refers to a dependencyon x ∈ R . Michael H. Gfrerer u ( θ , θ ) defined on the parameter space is related to the field u ( x )defined on the embedding space R byˆ u ( θ , θ ) = u ( x ) ◦ ˆ g ( θ , θ ) = u (ˆ g ( θ , θ )) (14)By applying the chain rule we find that the first and second derivatives arerelated by ˆ u ,θ = ( ∇ u ◦ ˆ g ) · ˆ g θ , (15)and ˆ u ,θτ = ( ∇∇ u ◦ ˆ g ) · ˆ g τ · ˆ g θ + ( ∇ u ◦ ˆ g ) · ˆ g θ,τ = ( ∇∇ u ◦ ˆ g ) · ˆ g τ · ˆ g θ + ( ∇ u ◦ ˆ g ) · ( Γ αθτ ˆ g α + h θτ ννν )= ( ∇∇ u ◦ ˆ g ) : (ˆ g τ ⊗ ˆ g θ ) + ( ∇ u ◦ ˆ g ) · ( Γ αθτ ˆ g α + h θτ ννν ) . (16)Furthermore, we have the following relations summarized in the followinglemma. Lemma 1
The metric tensor ˆ G = ˆ G αβ ˆ g α ⊗ ˆ g β = ˆ g α ⊗ ˆ g α and the projector P are related by ˆ G = P ◦ ˆ g . (17) For the Weingarten map we have the relation ˆ H = H ◦ ˆ g . (18) The proof can be found in Appendix A. f : U → T is given by ∇ Ω ˆ f = ˆ f ,α ⊗ ˆ g α . (19) Lemma 2
For the representation f : R → T the surface gradient is given by ∇ Ω f = ∇ f · P . (20) Proof
Using the relation between the projector and the metric tensor and (15)we have ∇ Ω f ◦ ˆ g = ( ∇ f ◦ ˆ g ) · (ˆ g α ⊗ ˆ g α ) = ˆ f ,α ⊗ ˆ g α . (21) race-Finite-Cell-Method for thin shell analysis 7 T = 1 (cid:112) det ˆ G (cid:16) ˆ T · ˆ g α (cid:112) det ˆ G (cid:17) ,α , (22)where we use the notation det ˆ G = det([ ˆ G αβ ]). The next lemma gives thesimpler representation for the surface divergence in case of a surface embeddedin R . Lemma 3
On a surface Ω ⊂ R parametrized by ˆ g : U → Ω the surfacedivergence of a tensor-valued function represented by ˆ T : U → T is given by div ˆ T = ˆ T ,α · ˆ g α + H ˆ T · ˆ ννν. (23) Furthermore, for the representation T : R → T we have div T = ∇ T : P + H T · ννν. (24) The proof can be found in Appendix A.
In the following lemma we collect product rules for the divergence operator.
Lemma 4
Let v × T be the cross product of a vector v = v i e i and a secondorder tensor T = T lk e l ⊗ e k defined by v × T = v i T lk ( e i × e l ) ⊗ e k , (25) and V ·× T the scalar-cross product of two second order tensors V = V ij e i ⊗ e j and T = T lk e l ⊗ e k defined by V ·× T = V ij T lk ( e j · e l )( e i × e k ) . (26) Then, the following product rules hold div( v × T ) = v × div( T ) + ∇ Ω v ·× T (cid:62) , (27)div( v · T ) = v · div( T ) + ∇ Ω v : T (cid:62) . (28) The proof can be found in Appendix A. T (cid:90) Ω div T d x = (cid:90) Γ T · µµµ d s x . (29)Using (28) and (29), the integration by parts formula for a vector v and asecond order tensor T reads (cid:90) Ω v · div T d x = (cid:90) Γ v · T · µµµ d s x − (cid:90) Ω ∇ Ω v : T (cid:62) d x . (30) Michael H. Gfrerer
In this section we derive the governing equations of linear thin shells from firstprinciples of continuum mechanics. Furthermore, we show the equivalence tothe linear Koiter model formulated as a minimization problem.3.1 Shell kinematicsThe kinematics of the surface Ω is described by the change in metric tensorand the change in curvature tensor. In the present paper we focus on the lineartheory and use the linearized change in metric tensor ˆ γγγ = ˆ γ αβ ˆ g α ⊗ ˆ g β andthe linearized change in curvature tensor ˆ ρρρ = ˆ ρ αβ ˆ g α ⊗ ˆ g β . The respectivecovariant components are given by [14,10]ˆ γ αβ (ˆ u ) = 12 (ˆ u ,β · ˆ g α + ˆ u ,α · ˆ g β ) , (31)and ˆ ρ αβ (ˆ u ) = ˆ ννν · (cid:0) ˆ u ,αβ − Γ σαβ ˆ u ,σ (cid:1) . (32)The next lemma establishes the representations for γγγ = ˆ γγγ ◦ ˆ g − and ρρρ = ˆ ρρρ ◦ ˆ g − . Lemma 5
For the linearized change in metric tensor we have the represen-tation γγγ ( u ) = 12 P · ( ∇ u + ( ∇ u ) (cid:62) ) · P , (33) and for the linearized change in curvature tensor we have the representation ρρρ ( u ) = P · ( ννν · ∇∇ u ) · P − ( ννν · ∇ u · ννν ) H . (34) The proof can be found in Appendix A. t ( µµµ ) on a cut defined by the boundary normal µµµ tangential to the surface. Due to Cauchys theorem we have the representation t ( µµµ ) = σσσ · µµµ, (35)with the stress tensor σσσ . We decompose the stress tensor σσσ in a tangential anda normal part, σσσ = N + ννν ⊗ S , (36)with the tangential stress tensor N = N αβ g α ⊗ g β and the vector S = S α g α related to transverse shear. Analogously, we have a moment vector m ( µµµ ) onthe cut defined by µµµ , which can be expressed as m ( µµµ ) = ννν × ( M · µµµ ) , (37)with the tangential moment tensor M . race-Finite-Cell-Method for thin shell analysis 9 (cid:90) Γ t d s x + (cid:90) Ω b d x = 0 . (38)Applying the surface divergence theorem (29) results in (cid:90) Ω div σσσ + b d x = 0 . (39)Due to the fact that Ω is arbitrary the local force equilibrium readsdiv σσσ + b = 0 . (40)3.4 Equilibrium of momentsThe equilibrium of moments states that the sum of boundary moments, themoments of boundary tractions, and the moments due to surface loads van-ishes, (cid:90) Γ m + x × t d s x + (cid:90) Ω x × b d x = 0 . (41)The following lemma summarizes the consequences of the equilibrium of mo-ments. Lemma 6
For T = T ij e i ⊗ e j let [ T ] × = T ij e i × e j . The equilibrium ofmoments is fulfilled if [ − H · M + N (cid:62) ] × = 0 , (42) and S = P · div( M ) . (43) The proof can be found in Appendix A. M = − t E : ρρρ, (44)and N = t E : γγγ − H · M . (45)The fourth order elasticity tensor E is given by E = λ ( P ⊗ P ) + 2 µ P s , with λ = 4¯ λµ ¯ λ + 2 µ , (46) where ¯ λ and µ are the Lam´e constants of the elastic material constituting theshell and P s the symmetric part of the tangential fourth order identity tensor.The Lam´e constants are related to the Young’s modulus E and Poisson’s ratio ν by ¯ λ = Eν (1 + ν )(1 − ν ) , µ = E ν ) . (47)The constitutive equations can also be write as N = ¯ N − H · M , M = − t
12 ( λ P tr ρρρ + 2 µ ρρρ ) , ¯ N = t ( λ P tr γγγ + 2 µ γγγ ) . (48)We remark that with (45) the condition (42) is fulfilled identically.3.6 Weak form of the governing equationsThe weak form of the governing equations is given by t (cid:90) Ω γγγ ( v ) : E : γγγ ( u ) d x + t (cid:90) Ω ρρρ ( v ) : E : ρρρ ( u ) d x = (cid:90) Ω v · b d x + (cid:90) Γ Ni v i N Ni d s x − (cid:90) Γ Nt ∇ Ω ( ννν · v ) · t M Nt d s x − (cid:90) Γ Nµ ∇ Ω ( ννν · v ) · µµµ M Nµ d s x , (49)where v are appropriate test functions, see Appendix B. The boundary con-ditions which can be prescribed are given by, u · e i = u Di or e i · ( ¯ N + ννν ⊗ S ) · µµµ = N Ni , ∇ Ω ( ννν · u ) · t = ω t or t · M · µµµ = M Nt , ∇ Ω ( ννν · u ) · µµµ = ω µ or µµµ · M · µµµ = M Nµ , (50)where u Di is a given displacement in the direction of e i , ω t and ω µ are givenrotations, N Ni is a given force in the direction of e i , M Nt and M Nµ are givenmoments. If u is prescribed on the boundary, the derivative along the boundary d t ( ννν · u ) = ∇ Ω ( ννν · u ) · t is also prescribed [4]. Thus, in this case, only the normalderivative d µµµ ( ννν · u ) = ∇ Ω ( ννν · u ) · µµµ can be independently prescribed by ω µ . race-Finite-Cell-Method for thin shell analysis 11 E ( v ) = 12 (cid:20)(cid:90) Ω t γγγ ( v ) : E : γγγ ( v ) + t ρρρ ( v ) : E : ρρρ ( v ) d x (cid:21) − (cid:90) Ω b · v d x − (cid:90) Γ Ni v i N Ni d s x + (cid:90) Γ Nt ∇ Ω ( ννν · v ) · t M Nt d s x + (cid:90) Γ Nµ ∇ Ω ( ννν · v ) · µµµ M Nµ d s x . (51)Then, the linear Koiter shell model reads: Find u ∈ V such that E ( u ) = inf v ∈V E ( v ) . (52)Since the variational equations of (51) are the equations given in (49) we haveestablished the equivalence. Therefore, we conclude that the Koiter modelproposed out of purely mechanical and geometrical intuitions can be derivedfrom first principles of continuum mechanics. C -Trace-Finite-Cell-Method For the discretization of the weak form (49) we propose a C -continuity versionof the TraceFEM. Following the TraceFEM concept the ansatz space on thesurface is defined as the restriction (trace) of an outer ansatz space defined ona background mesh. We label the present method also a Finite-Cell methodbecause we use as a background mesh a Cartesian grid. On this structured gridwe are able to construct C -continuity shape functions by the tensor productof univariate cubic Hermite form functions. Locally, they are defined on theon the unit interval by (see Figure 3) ϕ ( ξ ) = 1 + ξ (2 ξ − ,ϕ ( ξ ) = ξ ( ξ ( ξ −
2) + 1) ,ϕ ( ξ ) = − ξ (2 ξ − ,ϕ ( ξ ) = ξ ( ξ − . (53)The local functions (53) are pieced together to global C shape functions N i ( x ) by the standard finite element procedure of relating local degrees offreedom and global degrees of freedom. The degrees of freedom are the valuesand the first derivatives at the vertices of the background grid. Thus, at eachvertex we have 2 = 8 degrees of freedom for the discretization of a scalarfield and 24 degrees of freedom for the vector-valued displacement field. Thus,on a background cell we have 64 local form functions and 192 local degrees offreedom for the displacement field. We denote the vector-valued finite element . . . . . ϕ ϕ ϕ ϕ ξ ϕ i ( ξ ) Fig. 3
Local univariate cubic Hermite form functions space on the background grid by V h . Then, the discrete problem is given by:Find u h ∈ V h such that M h ( u h , w h ) = f D ( w h ) , (54a)holds for all w h ∈ V h and K h ( u h , v h ) = f N ( v h ) , (54b)holds for all v h with M h ( v h , w h ) = 0 for all w h ∈ V h . The linear and bilinearforms are given by M h ( u h , v h ) = (cid:90) Γ Di v h · u h d s x + (cid:90) Γ Dc [ ∇ Ω ( ννν · v h ) · µµµ ][ ∇ Ω ( ννν · u h ) · µµµ ] d s x , (55a) f D ( v h ) = (cid:90) Γ Di v h · u D i d s x + (cid:90) Γ Dc [ ∇ Ω ( ννν · v h ) · µµµ ] u D c d s x , (55b) K h ( u h , v h ) = t (cid:90) Ω γγγ ( v h ) : E : γγγ ( u h ) d x + t (cid:90) Ω ρρρ ( v h ) : E : ρρρ ( u h ) d x , (55c) f N ( v h ) = (cid:90) Ω v h · b d x + (cid:90) Γ Ni v h,i N Ni d s x − (cid:90) Γ Nt ∇ Ω ( ννν · v ) · t M Nt d s x − (cid:90) Γ Nµ ∇ Ω ( ννν · v ) · µµµ M Nµ d s x . (55d)In (54a) the solution u h gets fixed to prescribed values at the Dirichlet bound-ary. However, the matrix M h (and also K h ) is singular by definition because oftwo reason. First, in order to avoid further notation, we take w h ∈ V h resultingin zero rows and columns, which can be easily identified. Secondly, since wedefine the shape functions by restriction of the shape functions defined on thebackground mesh, they are not linear independent. In the standard setting ofa fitted finite element method the respective degrees of freedom in (54a) areeasy to identify and can be determined by interpolation. Here, in the case of an race-Finite-Cell-Method for thin shell analysis 13 unfitted method a linear independent basis is not known explicitly in generaland have to be determine. In (54b) the test functions v h are restricted to thenull-space of M h .In the discrete method the integrals in (55) are evaluated by quadrature,which is described in the next section.4.1 Integral evaluationIn order to evaluate the surface and line integrals in (55) we use the quadratureschema developed in [43]. Here, we outline only the main ingredients and referto [43] for technical details. Following the standard finite element procedurethe integrals are evaluated by summing up background cell (face) contribu-tions where the shape functions are smooth. The individual contributions areevaluated by Gaussian quadrature. The main idea from [43] is to subdividethe background cells (faces) until it is possible to convert the implicitly de-fined geometry into the graph of an implicitly defined height function. Then,a recursive algorithm which requires only one-dimensional root finding andone-dimensional Gaussian quadrature can be set up. In order to chose suitableheight function directions we have to be able to ensure the monotonicity ofthe level-set function in that direction. This can be done by showing that thederivative in that direction is uniform in sign, i.e. place bounds on the valuesattainable by the derivative. In contrast to [43] we use interval arithmetic forthis task.4.2 Solution strategiesIn order to solve (54) we consider the three methods, – Null-space method, – Penalty method, and – Lagrange multiplier method.In the null-space method we first solve M h u Dh = f D , (56)and compute the null-space of M h . We denote the null-space basis by Z h . Ina next step the solution u h of( Z (cid:62) h K h Z h ) u h = Z (cid:62) h ( f N − K h u Dh ) (57)is computed. The overall solution is then given by u h = Z h u h + u Dh . (58)In the penalty method we solve a system of linear equations of the form( K h + αM h ) u h = f N + αf D , (59) with the penalty parameter α >
0. In the Lagrange multiplier method we solvethe system of linear equation (cid:20) K h M h M h (cid:21) (cid:20) u h λ h (cid:21) = (cid:20) f N f D (cid:21) , (60)with the Lagrange multipliers λ h .4.3 ImplementationThe proposed method has been implemented in Matlab. Within the methodthe exact level-set function φ ( x ) is used. For the evaluation of the surfacenormal vector (5) and the Weingarten map (7), the first and second orderderivatives of the level-set function are necessary. In the implementation weuse symbolic differentiation of φ ( x ) to provide these derivatives.In the present work, we have not used any stabilization term which isadded to the weak form. Therefore, in each system of equations (56), (57), (59)and (60) the system matrix is singular by definition. One strategy to wouldbe to add stabilization terms to the bilinear forms (54b) and (54b). We referto [34] for an overview of different possibilities. Although such a stabilizationcan be designed in such a way that the convergence order of the method is notbeen altered, a stabilization decreases the accuracy of the method. Therefore,a strategy is investigated, where no stabilization is necessary. However, wehave observed in numerical experiments that the Matlab backslash operatordoes not give satisfactory results. Due to this reason, we use the direct solversuitable for under-determined linear equation systems from the SuiteSparse project. In this section, numerical results are presented. First, we verify the implemen-tation of the proposed method against exact manufactured solutions. Secondly,we demonstrate that the method works by testing it with two benchmarkproblems (one cylindrical shell and one spherical shell) of the well-known shellobstacle course [7]. Finally, in a fourth example, a complex shell geometry isinvestigated.5.1 Verification exampleThe implementation of the proposed method is verified. We have successfullyrun the method on various geometries, displacement fields and boundary con-dition combinations, but present only the results of two configurations. In allverification examples we used B = [ − . , . × [0 , × [ − . ,
1] and a shell http://faculty.cse.tamu.edu/davis/suitesparse.htmlrace-Finite-Cell-Method for thin shell analysis 15 thickness t = 0 .
25. The material parameters are E = 10 and ν = 0 .
4. Theconsidered geometries are defined by the zero level-set of the functions φ ( x, y, z ) = x + z − . , (61a) φ ( x, y, z ) = 4 x + 0 . y + z − . , (61b)and are illustrated in Figure 4. The level-set function (61a) implies a flat (a) φ ( x, y, z ) = x + z − . φ ( x, y, z ) = 4 x + 0 . y + z − . Fig. 4
Problem geometries of the verification geometry, whereas (61b) implies a surface with varying curvature.As manufactured solutions we consider the two displacement fields u ex ( x, y, z ) = x e x + y x e y + ( xzy + x ( x − y ( y − e z , (62a) u ex ( x, y, z ) = sin(16 y ) cos(16 x z ) e x + cos(16 x y z ) e y + 2 sin(16 x y z ) e z . (62b)Using (40), we compute symbolically the necessary surface force b such that(62) is the respective exact solution of the shell problem. The displacementfield u is chosen as a third order polynomial such that the solution can be rep-resented exactly in the discrete space, whereas u can only be approximated.In the following, we study the behavior of the error e = (cid:115) (cid:82) Ω ( u h − u ex ) d x (cid:82) Ω ( u ex ) d x (63)under uniform mesh refinement for the three different solution methods givenin Section 4.2. Furthermore, for comparison, we also consider the Hermiteinterpolation u inth of the solution on the background grid and the surface L -projection: Find u L h ∈ V h such that (cid:90) Ω ( u ex − u L h ) d x → min . (64) We remark that the Hermite interpolation and the L -projection are onlypossible if the solution is known, and is thus only computed for the verificationexamples in this section. Furthermore, all solutions apart from the Hermiteinterpolation require the solution of a system of linear equations.The numerical results for the flat plate defined in (61a) are visualized inFigures 5 and 6. The refinement level 0 refers to a single background cell. InFigure 5 the errors obtained by the penalty method are given for differentrefinement levels and penalty factors. We remark that for this flat geometry (a) displacement u (b) displacement u Fig. 5
Errors for different penalty factors on the geometry induced by φ the numerical integration is exact up to round of errors. Therefore, for u the sources for errors are the error due to the imposition of the boundaryconditions by the penalty method (depending on the penalty factor) and roundoff errors in the numerical computations. As expected the errors decrease withincreasing penalty parameter up to a point where the ill-conditioning of thelinear system dominates the error. Since u can not be represented exactly inthe discrete space an approximation error limits the overall error. In Figure 6the results for the different solution methods are given. For the penalty methodwe used the lowest errors of the results shown in Figure 5. As no system ofequation has to be solved for the interpolation on the background grid ( volumeinterpolation ) the error is around 10 − for u for all refinement levels. Incontrast to this the other results require the solution of a system of equationsand therefore the errors are between 10 − and 10 − due to round off errors.For u we observe the convergence of all methods with optimal rate. Here,by definition the L -projection gives the lowest error, whereas the Hermiteinterpolation results in the highest error for a fixed refinement level (apartfrom level 5, where the error due to ill-conditioning of the system of equationsdominates).The numerical results for the part of the ellipse defined in (61b) are vi-sualized in Figures 7 and 8. In Figure 7 the results of the penalty method race-Finite-Cell-Method for thin shell analysis 17 Fig. 6
Results for φ and u (left) and u (right)(a) displacement u (b) displacement u Fig. 7
Errors for different penalty factors on the geometry induced by φ for different penalty parameters are given. In contrast to the flat geometry,now the numerical integration is not exact, yielding an additional error, whichdominates for the three coarsest levels.In Figure 8 the errors for the three different solution methods, the Her-mite interpolation and the L -projection are visualized. Again, for the penaltymethod we used the lowest errors of the results shown in Figure 7. Similaras before, for the displacement field u the volume interpolation gives errorsaround 10 − , whereas for the other solutions the errors are between 10 − and10 − due to round off errors. For displacement field u we observe the conver-gence of all methods. However, due to round off errors the accuracy is limited.Nevertheless, we remark that an relative error level of about 10 − is more thansufficient for practical problems. This can also be seen from the visualizationsof the solutions obtained with the null-space method for u in Figure 9 and Fig. 8
Results for φ and u (left) and u (right) for u in Figure 10. For the fine levels no difference in the solutions can beseen by eye. (a) level 0 (b) level 1 (c) level 2(d) level 3 (e) level 4 (f) level 5 Fig. 9
Visualization of the displacement results for φ and u race-Finite-Cell-Method for thin shell analysis 19(a) level 0 (b) level 1 (c) level 2(d) level 3 (e) level 4 (f) level 5 Fig. 10
Visualization of the displacement results for φ and u r = 25) is supported by rigid diaphragms at theends ( x = 0 and x = 50), i.e. u y = u z = 0. The straight edges are free.The geometry and the material parameters are depicted in Figure 11. Thestructure is subjected to gravity loading with b = − e z . We describe theproblem geometry by φ ( x, y, z ) = y + z − r , (65)and B = [0 , × [ − r sin( π ) , r sin( π )] × [10 , . A , which is located in themiddle of one free edge. As a reference solution we use the overkill solution u Az = − . E = 4 . · ν = 0 R = 25 L = 50 t = 0 . Fig. 11
Problem description of the Scordelis-Lo roof problem
Fig. 12
Deformed configuration of the Scordelis-Lo roof (displacements scaled by a factorof 10) φ ( x, y, z ) = x + y + z − R , (66)with R = 10 and B = [ − . , . × [ − . , . × [0 , . race-Finite-Cell-Method for thin shell analysis 21 Table 1
Vertical displacements of the Scordelis-Lo roof at point A . Refernece: u Az = − . E = 6 . · ν = 0 . R = 10 t = 0 . F = 2 Fig. 13
Problem description of the pinched hemisphere problemref. level null-space Lagrange mulitplier penalty method0 0.04490 0.04490 0.044901 0.08989 0.08989 0.089902 0.09222 0.09222 0.092223 0.09239 0.09239 0.092394 0.09241 0.09241 0.092405 0.09241 0.09241 0.09241
Table 2
Displacements of the pinched hemisphere at one loading point. Refernece: u r =0 . of the hemisphere is unconstrained and the four radial forces have alternatingsigns such that the sum of the applied forces is zero. We investigate the radialdisplacement at the loaded points. In [7], the reference displacement of u r =0 . Fig. 14
Deformed configuration of the pinched hemisphere (displacements scaled by a factorof 40) φ ( x, y, z ) = sin( πx ) cos( πy ) + sin( πy ) cos( πz ) + sin( πz ) cos( πx ) . (67)The considered shell lies in the cuboid B = [0 , × [ − . , . × [ − . , . x = 0.We assume a thickness t = 0 .
03. We study the vertical deflection due to avolume load b = − e z at the point [2 , . , − . u z = − . race-Finite-Cell-Method for thin shell analysis 23 E = 7 · ν = 0 . t = 0 . Fig. 15
Geometry of the gyroid problem. The structure is clamped at the gray plane.
Fig. 16
Deformed configuration of the gyroid
Table 3
Displacements u z of the gyroid at the point [2,0.5,-0.25].ref. level null-space Lagrange mulitplier penalty method0 -0.24147 -0.54289 -0.316391 -1.70309 -2.05171 -1.715212 -1.80865 -2.40048 -1.809003 -1.80891 -2.61238 -1.809254 -1.80905 3.77214 -1.80318 We have developed a C -continuous finite element method for thin shells withmid-surface given as the zero level-set of a scalar function. In order to achieve the continuity of the discretization, concepts of the TraceFEM and the Finite-Cell-Method are combined. In particular the shape functions on the shell sur-face are obtained by restriction of tensor-product cubic Hermite splines on astructured background mesh. In order to allow a natural implementation, theunderlying shell model is formulated in a parametrization-free way. Further-more, the strong form of the governing equations are given. This allows toobtain manufactured solutions on arbitrary geometries. Thus, the implemen-tation of the proposed method is verified by a convergence analysis where theerror is computed with an exact manufactured solution.In the present method, the shape functions on the shell surface are lin-early dependent. In order to avoid a singular system matrix, a stabilizationterm can be used. In the presented method such a stabilization is avoided.However, it is necessary to use the direct solver suitable for under-determinedlinear equation systems from the SuiteSparse project. We investigated threestrategies to include the boundary conditions. These are the penalty method,the Lagrange multiplier method, and the null-space method. In the numericalexperiments we have observed that the penalty method and the null-spacemethod give reliable results. However, the Lagrange multiplier method suffersfrom instabilities in some examples, which should be further investigated. Infuture work, it would be also interesting to use iterative solvers in contrast tothe used direct solver.In contrast to thin shells, for Reissner-Mindlin shells only C -continuousshape functions are commonly used. In order to avoid transverse shear locking,in [31,22] an hierarchic concept of shell models is presented. This approach hasthe advantage that transverse shear locking is eliminated on the continuousformulation level, independent of a particular discretization, but requires C -continuous shape functions. An extension of the present work to Mindlin-Reissner shells with implicitly defined mid-surface seems possible and wouldbe worth to investigate. Acknowledgements
The author thanks Thomas-Peter Fries from the Institute of Struc-tural Analysis at TU Graz, and Helmut Gfrerer from the Institute of Computational Math-ematics at JKU Linz for valuable discussions on the topic of the paper.
Conflict of interest
The author declares that he has no conflict of interest. race-Finite-Cell-Method for thin shell analysis 25
A Proofs
Proof (Proof of Lemma 1)
It is sufficient to show that the application of the operators ˆ G and P ◦ ˆ g = I − ˆ ννν ⊗ ˆ ννν to a basis (ˆ g l , ˆ ννν ) gives the same result, G · ˆ g l = ˆ g l , ( P ◦ ˆ g ) · ˆ g l = ˆ g l , G · ννν = , ( P ◦ ˆ g ) · ννν = . Furthermore, direct calculation shows H ◦ ˆ g = − ( ∇ ννν · P ) ◦ ˆ g = − ( ∇ ννν ◦ ˆ g ) · ˆ g α ⊗ ˆ g α = − ˆ ννν ,α ⊗ ˆ g α = ˆ H . Proof (Proof of Lemma 3)
First we establish the relation (cid:112) det ˆ G ,γ = (cid:112) det ˆ G ˆ Γ αγα , (68)Following [14, Theorem 4.4-4], we have (cid:113) det( ˆ G ) = det(ˆ g , ˆ g , ˆ ννν ) (69)and the sought relation follows by (cid:18)(cid:113) det( ˆ G ) (cid:19) ,α = det(ˆ g ,α , ˆ g , ˆ ννν ) + det(ˆ g , ˆ g ,α , ˆ ννν ) + det(ˆ g , ˆ g , ˆ ννν ,α )= det( ˆ Γ ϕ α ˆ g ϕ + ˆ h α ˆ ννν, ˆ g , ˆ ννν ) + det(ˆ g , ˆ Γ ϕ α ˆ g ϕ + ˆ h α ˆ ννν, ˆ ννν )+ det(ˆ g , ˆ g , − ˆ h ϕα ˆ g ϕ )= ( ˆ Γ α + ˆ Γ α ) det(ˆ g , ˆ g , ˆ ννν )= ˆ Γ ββα (cid:113) det( ˆ G ) . (70)With (12), we obtain (cid:18) ˆ g α (cid:113) det( ˆ G ) (cid:19) ,α = ˆ g α,α (cid:113) det( ˆ G ) + ˆ g α (cid:18)(cid:113) det( ˆ G ) (cid:19) ,α = ( − Γ ααγ ˆ g γ + ˆ h αα ˆ ννν ) (cid:113) det( ˆ G ) + ˆ g α ˆ Γ ββα (cid:113) det( ˆ G )= H ˆ ννν (cid:113) det( ˆ G ) . (71)Therefore, the first part of the lemma followsdiv ˆ T = 1 (cid:112) det ˆ G (cid:16) ˆ T · ˆ g α (cid:112) det ˆ G (cid:17) ,α = ˆ T ,α · ˆ g α + (cid:18) ˆ g α (cid:113) det( ˆ G ) (cid:19) ,α (cid:113) det( ˆ G )= ˆ T ,α · ˆ g α + H ˆ T · ˆ ννν. (72)The second part of the lemma can be shown by direct calculation,(div T ) ◦ ˆ g = ( ∇ T ◦ ˆ g ) : (ˆ g α ⊗ ˆ g α ) + H ˆ T · ˆ ννν = (( ∇ T ◦ ˆ g ) · ˆ g α ) · ˆ g α + H ˆ T · ˆ ννν = ˆ T ,α · ˆ g α + H T · ννν = div ˆ T . (73)6 Michael H. Gfrerer Proof (Proof of Lemma 4)
The lemma can be shown by the direct calculationsdiv( v × T ) = ∇ Ω ( v × T ) : P + H v × T · ννν = [( v × T ) ,i ⊗ e i · P ] : P + H v × T · ννν = [( v ,i × T + v × T ,i ) ⊗ e i · P ] : P + H v × T · ννν = ( − T (cid:62) × v ,i ⊗ e i · P + v × ∇ Ω T ) : P + H v × T · ννν = ∇ Ω v ·× T (cid:62) + v × div T , (74)and div( v · T ) = ∇ Ω ( v · T ) : P + H v · T · ννν = [( v · T ) ,i ⊗ e i · P ] : P + H v · T · ννν = [( v ,i · T + v · T ,i ) ⊗ e i · P ] : P + H v · T · ννν = ( T (cid:62) · v ,i ⊗ e i · P + v · ∇ Ω T ) : P + H v · T · ννν = ∇ Ω v : T (cid:62) + v · div T . (75) Proof (Proof of Lemma 5)
For the proof we use the relations (15) and (16). Direct calcula-tion yields for the linearized change in metric tensor γγγ ◦ ˆ g = (cid:20) P · ( ∇ u + ( ∇ u ) (cid:62) ) · P (cid:21) ◦ ˆ g = 12 (ˆ g α ⊗ ˆ g α ) · ( ∇ u ◦ ˆ g + ( ∇ u ) (cid:62) ◦ ˆ g ) · (ˆ g β ⊗ ˆ g β )= 12 ( u ,β · ˆ g α + u ,α · ˆ g β )ˆ g α ⊗ ˆ g β = ˆ γγγ, (76)and for the linearized change in curvature tensor ρρρ ◦ ˆ g = [ P · ( ννν · ∇∇ u ) · P − ( ννν · ∇ u · ννν ) H ] ◦ ˆ g = (ˆ g α ⊗ ˆ g α ) · (ˆ ννν · ∇∇ u ◦ ˆ g ) · (ˆ g β ⊗ ˆ g β ) + ˆ ννν · ( ∇ u ◦ ˆ g ) · ˆ νννh αβ (ˆ g α ⊗ ˆ g β )= ˆ ννν · u ,αβ ˆ g α ⊗ ˆ g β − ˆ ννν · ( ∇ u ◦ ˆ g ) · ( Γ γαβ ˆ g γ + h αβ ˆ ννν )ˆ g α ⊗ ˆ g β + ˆ ννν · ( ∇ u ◦ ˆ g ) · ˆ ννν h αβ (ˆ g α ⊗ ˆ g β )= ˆ ννν · ( u ,αβ − Γ γαβ u ,γ )ˆ g α ⊗ ˆ g β = ˆ ρρρ. (77) Proof (Proof of Lemma 6)
Applying the surface divergence theorem yields (cid:90) Ω div( ννν × M ) + div( x × σσσ ) + x × b d x = 0 . (78)Using the divergence product rule (27) results in (cid:90) Ω ννν × div( M ) + ∇ Ω ννν ·× M (cid:62) + x × div σσσ + ∇ Ω x ·× σσσ (cid:62) + x × b d x = 0 . (79)With (40), and ∇ Ω x = P we have (cid:90) Ω ννν × div( M ) + ∇ Ω ννν ·× M (cid:62) + P ·× σσσ (cid:62) d x = 0 . (80)Due to the definition of the stress tensor (36) we have σσσ T = N T + S ⊗ ννν and it follows P ·× σσσ (cid:62) = ( g α ⊗ g α ) ·× ( N T + S ⊗ ννν )= ( g α ⊗ g α ) ·× ( N βγ g γ ⊗ g β + S γ g γ ⊗ ννν )= N βα ( g α × g β ) + S × ννν = [ N (cid:62) ] × + S × ννν, (81)race-Finite-Cell-Method for thin shell analysis 27and furthermore ∇ Ω ννν ·× M (cid:62) = ( ννν ,α ⊗ g α ) ·× ( M βγ g β ⊗ g γ )= − ( h ϕα g ϕ ⊗ g α ) ·× ( M βγ g β ⊗ g γ )= − h ϕβ M βγ g ϕ × g γ = [ − H · M ] × . (82)Thus, (cid:90) Ω ννν × div( M ) + ∇ Ω ννν ·× M (cid:62) + P ·× σσσ (cid:62) d x = (cid:90) Ω ννν × (div( M ) − S ) + [ − H · M + N (cid:62) ] × d x = 0 . (83)From (83) we deduce the sought conditions[ − H · M + N (cid:62) ] × = 0 , (84)and S = P · div( M ) . (85) B Derivation of the weak form
In this section the weak form of the governing equations is derived. To this end, we multiply(40) with a test function v ∈ V and integrate over the shell surface, − (cid:90) Ω v · div σσσ d x = (cid:90) Ω v · b d x . (86)Here, the function space of the test functions is V = { ηηη : Ω → R | γγγ ( ηηη ) ∈ L ( Ω, R ) , ρρρ ( ηηη ) ∈ L ( Ω, R ) ,ηηη · e i = 0 on Γ D i , ∇ Ω ( ννν · ηηη ) · µµµ = 0 on Γ D c } , (87)where Γ D i , Γ D t , and Γ D µ denote Dirichlet boundaries. On Γ D i the displacement in direc-tion e i is restrained and on Γ D t and Γ D µ the rotation of the shell around the boundarytangent vector and the boundary normal vector is restrained respectively. The correspondingNeumann boundaries are given by Γ N i = Γ \ Γ D i , Γ N t = Γ \ Γ D t , and Γ N µ = Γ \ Γ D µ .Integration by parts of the term on the left side yields (cid:90) Ω ∇ Ω v · σσσ (cid:62) d x = (cid:90) Γ v · σσσ · µµµ d s x + (cid:90) Ω v · b d x . (88)We have σσσ (cid:62) = N (cid:62) + S ⊗ ννν and obtain (cid:90) Ω ∇ Ω v · N (cid:62) d x + (cid:90) Ω ( ννν · ∇ Ω v ) · S d x = (cid:90) Γ v · σσσ · µµµ d s x + (cid:90) Ω v · b d x . (89)Using (43) and integration by parts of the second term on the left yields (cid:90) Ω ( ννν · ∇ Ω v ) · S d x = (cid:90) Γ ( ννν · ∇ Ω v ) · M · µµµ d s x − (cid:90) Ω ∇ Ω ( ννν · ∇ Ω v ) : M d x . (90)Due to (48) we obtain the relation (cid:90) Ω ∇ Ω v : N (cid:62) d x = (cid:90) Ω ∇ Ω v : ¯ N d x − (cid:90) Ω ∇ Ω v : ( M · H ) d x = (cid:90) Ω γγγ ( v ) : E : γγγ ( u ) d x − (cid:90) Ω ( H · ∇ Ω v ) : M d x . (91)8 Michael H. GfrererFurthermore, we have − [ H · ∇ Ω v + ∇ Ω ( ννν · ∇ Ω v )] : M = − ρρρ ( v ) : M , (92)and (cid:90) Γ ( ννν · ∇ Ω v ) · M · µµµ d s x − (cid:90) Γ v · ( H · M ) · µµµ d s x = (cid:90) Γ ∇ Ω ( ννν · v ) · M · µµµ d s x . (93)Therefore, by using (91), (92), and (93) we obtain from (86) the final weak form t (cid:90) Ω γγγ ( v ) : E : γγγ ( u ) d x + t (cid:90) Ω ρρρ ( v ) : E : ρρρ ( u ) d x = (cid:90) Ω v · b d x + (cid:90) Γ Ni v i N Ni d s x − (cid:90) Γ Nt ∇ Ω ( ννν · v ) · t M Nt d s x − (cid:90) Γ Nµ ∇ Ω ( ννν · v ) · µµµ M Nµ d s x . (94) References
1. Cut2. Areias, P.M.A., Song, J.H., Belytschko, T.: A finite-strain quadrilateral shell elementbased on discrete Kirchhoff-Love constraints. International Journal for Numerical Meth-ods in Engineering (9), 1166–1206 (2005)3. Argyris, J.H., Fried, I., Scharpf, D.W.: The TUBA family of plate elements for thematrix displacement method. The Aeronautical Journal (692), 701–709 (1968)4. Basar, Y., Kr¨atzig, W.B.: Mechanik der Fl¨achentragwerke: Theorie, Berechnungsmeth-oden, Anwendungsbeispiele. Vieweg (1985)5. Batoz, J.L., Zheng, C.L., Hammadi, F.: Formulation and evaluation of new triangular,quadrilateral, pentagonal and hexagonal discrete Kirchhoff plate/shell elements. Inter-national Journal for Numerical Methods in Engineering (5-6), 615–630 (2001)6. Bell, K.: A refined triangular plate bending finite element. International Journal forNumerical Methods in Engineering (1), 101–122 (1969)7. Belytschko, T., Stolarski, H., Liu, W.K., Carpenter, N., Ong, J.S.: Stress projection formembrane and shear locking in shell finite elements. Comput Methods Appl Mech Eng (1), 221–258 (1985)8. Bieber, S., Oesterle, B., Ramm, E., Bischoff, M.: A variational method to avoid locking -independent of the discretization scheme. International Journal for Numerical Methodsin Engineering (8), 801–827 (2018)9. Bischoff, M., Bletzinger, K.U., Wall, W., Ramm, E.: Models and finite elements forthin-walled structures. In: E. Stein, R. de Borst, T. Hughes (eds.) Encyclopedia ofComputational Mechanics, vol. 2, chap. 3, pp. 59–137. Wiley Online Library (2004)10. Blouza, A., Le Dret, H.: Existence and uniqueness for the linear Koiter model for shellswith little regularity. Quarterly of Applied Mathematics (2), 317–337 (1999)11. Bogner, F.: The generation of interelement-compatible stiffness and mass matrices bythe use of interpolation formulas. In: Proc. Conf. Matrix Meth. Struct. Mech., pp.397–443. Wright-Patterson AFB (1965)12. Burman, E., C, S., Hansbo, P., Larson, M., Massing, A.: CutFEM: Discretizing geometryand partial differential equations. Int. J. Numer. Methods Eng. (7), 472–501 (2015)13. Burman, E., Hansbo, P., Larson, M.G., Massing, A.: Cut finite element methods forpartial differential equations on embedded manifolds of arbitrary codimensions. ESAIM:Mathematical Modelling and Numerical Analysis (6), 2247–2282 (2018)14. Ciarlet, P.G.: An Introduction to Differential Geometry with Applications to Elasticity,vol. 78. Springer (2006)15. Ciarlet, P.G., Lods, V.: Asymptotic analysis of linearly elastic shells. III. Justification ofKoiter’s shell equations. Archive for Rational Mechanics and Analysis (2), 191–200(1996)race-Finite-Cell-Method for thin shell analysis 2916. Cirak, F., Ortiz, M., Schr¨oder, P.: Subdivision surfaces: a new paradigm for thin-shellfinite-element analysis. International Journal for Numerical Methods in Engineering (12), 2039–2072 (2000)17. Clough, R.W.: Finite element stiffness matricess for analysis of plate bending. In: Proc.of the First Conf. on Matrix Methods in Struct. Mech., pp. 515–546 (1965)18. Delfour, M.C., Zol´esio, J.P.: Differential equations for linear shells: comparison betweenintrinsic and classical models. Advances in mathematical sciences: CRMs , 41–124(1997)19. Demlow, A.: Higher-order finite element methods and pointwise error estimates forelliptic problems on surfaces. SIAM J. Numer. Anal. (2), 805–827 (2009)20. Dominguez, V., Sayas, F.J.: Algorithm 884: A simple matlab implementation of theargyris element. ACM Transactions on Mathematical Software (TOMS) (2), 16:1–16:11 (2008)21. Dziuk, G., Elliott, C.: Finite element methods for surface PDEs. Acta Numer. ,289–396 (2013)22. Echter, R., Oesterle, B., Bischoff, M.: A hierarchic family of isogeometric shell finiteelements. Computer Methods in Applied Mechanics and Engineering , 170–180(2013)23. Engel, G., Garikipati, K., Hughes, T.J., Larson, M.G., Mazzei, L., Taylor, R.L.: Con-tinuous/discontinuous finite element approximations of fourth-order elliptic problemsin structural and continuum mechanics with applications to thin beams and plates, andstrain gradient elasticity. Computer Methods in Applied Mechanics and Engineering (34), 3669–3750 (2002)24. Gfrerer, M., Schanz, M.: Code verification examples based on the method of manufac-tured solutions for Kirchhoff–Love and Reissner–Mindlin shell analysis. Engineeringwith Computers (4), 775–785 (2018)25. Gfrerer, M., Schanz, M.: A high-order fem with exact geometry description for thelaplacian on implicitly defined surfaces. International journal for numerical methods inengineering (11), 1163–1178 (2018)26. Gfrerer, M.H., Schanz, M.: High order exact geometry finite elements for seven-parameter shells with parametric and implicit reference surfaces. Computational me-chanics pp. 1–13 (2018)27. Gurtin, M.E., Murdoch, A.I.: A continuum theory of elastic material surfaces. Archivefor rational mechanics and analysis (4), 291–323 (1975)28. Hansbo, P., Larson, M.G.: Continuous/discontinuous finite element modelling of Kirch-hoff plate structures in R using tangential differential calculus. Computational Me-chanics (4), 693–702 (2017)29. Kiendl, J., Bletzinger, K.U., Linhard, J., W¨uchner, R.: Isogeometric shell analysis withKirchhoff-Love elements. Computer Methods in Applied Mechanics and Engineering (49), 3902–3914 (2009)30. Koiter, W.T.: On the nonlinear theory of thin elastic shells. Proc. Koninkl. Ned. Akad.van Wetenschappen, Series B , 1–54 (1966)31. Long, Q., Burkhard Bornemann, P., Cirak, F.: Shear-flexible subdivision shells. Inter-national Journal for Numerical Methods in Engineering (13), 1549–1577 (2012)32. Naghdi, P.: Finite deformation of elastic rods and shells. In: Proceedings of the IUTAMSymposium on Finite Elasticity, pp. 47–103. Springer (1981)33. Neunteufel, M., Sch¨oberl, J.: The Hellan-Herrmann-Johnson method for nonlinear shells.Computers & Structures , 106109–106120 (2019)34. Olshanskii, M.A., Reusken, A.: Trace finite element methods for pdes on surfaces. In:S.P.A. Bordas, E. Burman, M.G. Larson, M.A. Olshanskii (eds.) Geometrically Un-fitted Finite Element Methods and Applications, pp. 211–258. Springer InternationalPublishing (2017)35. Olshanskii, M.A., Reusken, A., Grande, J.: A finite element method for elliptic equationson surfaces. SIAM Journal on Numerical Analysis (5), 3339–3358 (2009)36. Olshanskii, M.A., Safin, D.: Numerical integration over implicitly defined domains forhigher order unfitted finite element methods. Lobachevskii Journal of Mathematics (5), 582–596 (2016)0 Michael H. Gfrerer37. van Opstal, T., van Brummelen, E., van Zwieten, G.: A finite-element/boundary-element method for three-dimensional, large-displacement fluid–structure-interaction.Computer Methods in Applied Mechanics and Engineering , 637–663 (2015)38. Parvizian, J., Dster, A., Rank, E.: Finite cell method. Computational Mechanics (1),121–133 (2007)39. Pietraszkiewicz, W.: Geometrically nonlinear theories of thin elastic shells. Advancesin Mechanics , 51–130 (1989)40. Rafetseder, K., Zulehner, W.: A new mixed approach to Kirchhoff–Love shells. Com-puter Methods in Applied Mechanics and Engineering , 440–455 (2019)41. Rosenberg, S.: The Laplacian on a Riemannian manifold: An introduction to analysis onmanifolds. No. 31 in London Mathematical Society. Cambridge University Press (1997)42. Sauer, R.A., Duong, T.X.: On the theoretical foundations of thin solid and liquid shells.Mathematics and Mechanics of Solids (3), 343–371 (2017)43. Saye, R.: High-order quadrature methods for implicitly defined surfaces and volumes inhyperrectangles. SIAM J. Sci. Comput. (2), A993–A1019 (2015)44. Schillinger, D., Ruess, M.: The finite cell method: A review in the context of higher-orderstructural analysis of CAD and image-based geometric models. Archives of Computa-tional Methods in Engineering (3), 391–455 (2015)45. Sch¨ollhammer, D., Fries, T.P.: Kirchhoff–Love shell theory based on tangential differ-ential calculus. Computational mechanics64