A $C^1$-continuous Trace-Finite-Cell-Method for linear thin shell analysis on implicitly defined surfaces

A Trace-Finite-Cell-Method for the numerical analysis of thin shells is 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 in strong form based on first principles of continuum mechanics by recasting well-known relations formulated in local coordinates to a formulation independent of a parametrization. The field approximation is constructed by restricting shape functions defined on a structured background grid on the shell surface. As shape functions we use on a background grid the tensor product of cubic splines. This yields $C^1$-continuous approximation spaces, which are required by the governing equations of fourth order. The parametrization-free formulation allows a natural implementation of the proposed method and manufactured solutions on arbitrary geometries for code verification. Thus, the implementation is verified by a convergence analysis where the error is computed with an exact manufactured solution. Furthermore, benchmark tests are investigated.

, which is one of the most popular models for thin shells. The Koiter shell model has been justified by asymptotic analysis in [15], in the sense of being a reasonable approximation to the full 3D problem of a thin shell-like body. Existence and uniqueness results can be found in [10,14].
Classically, for theoretical treatment it is assumed that the shell midsurface is represented by a global parametrization. However, for the numerical treatment typically the mid-surface is approximated by, possibly curved, finite elements, i.e. represented by a collection of local parametrizations. In contrast to these representations, we consider the case where the mid-surface is represented implicitly as the zero-level set of a scalar function φ(x), see Figure 1 for the illustration of some examples. We refer to the review article [21] for an overview of finite element methods for problems on such surfaces. In the classical surface finite element method the discretization of the unknown field relies on the higher order or exact meshing (local parametrization) of the surface [19,26]. In contrast to this, in the proposed method the discretization of the displacement field does not rely on parametrizations. Therefore, we provide a throughout derivation of the governing equations based on first principles of continuum mechanics independent of a parametrization. This allows a natural implementation of the method and also the construction of manufactured solutions for code verification on arbitrary geometries. Equivalent derivations relying on a parametrization can be found in e.g. [4] and [42]. We remark that membrane and thin shell formulations without relying on a parametrization with a mathematical focus can be found in [27,18]. For a treatment from an engineering perspective we refer to [37,45]. Examples of implicitly defined surfaces. The surfaces are defined by the level-set functions (a) φ = x 2 + y 2 + z 2 − r 2 , (b) φ = (x 2 + y 2 + z 2 + R 2 − r 2 ) 2 − 4R 2 (x 2 + y 2 ), (c) φ = x 2 + y 2 − r 2 (d) φ = sin x cos y + sin y cos z + sin z cos x One of the main difficulties in developing finite element methods for thin shells is the construction of C 1 -continuous approximation spaces. For general unstructured meshes it is not possible to ensure C 1 -continuity with only local polynomial shape functions and the nodal degrees of freedom consist of displacements and slopes only [46]. However, different non-standard triangular for developed thin plate bending are the Argyris element [3,20], the Bell element [6] or the Clough-Tocher macrotriangle [17]. A further possibility to construct C 1 -continuous approximation spaces on general space triangulation relies on sophisticated techniques from subdivision surfaces [16]. However, on a structured quadrilateral mesh the Bogner-Fox-Schmit element [11] is a simple conforming element. The constraint of a structured quadrilateral mesh can be partially overcome by introducing a smooth mapping of the geometry [24]. This idea can be realized in an isoparametric way by the use of splines for the geometry mapping and for the discretization of the displacement field [29]. The general difficulty of constructing C 1 -continuous approximation spaces led to approaches where the C 1 -continuity requirement is circumvented. Among them we mention discrete Kirchhoff elements [5,2] where the Kirchhoff constraint is enforced only at discrete points, the use of shear-deformable (Reissner-Mindlin) shell theory, were only C 0 -continuity approximation spaces are required, mixed methods [40,33], continuous/discontinuous Galerkin methods [23,28] and others.
In the present paper, we combine ideas from unfitted finite element methods and the Bogner-Fox-Schmit element. Following the idea of the TraceFEM [35,34] (see also CutFEM [12,13]) the approximation of the displacement field is constructed by restricting shape functions defined on a background mesh on the 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 three uni-variant cubic spline shape functions leads to a C 1 -continuous approximation in 3D (like the Bogner-Fox-Schmit element in 2D). Therefore, the shape functions for approximation of the displacement field on the shell mid-surface are C 1 -continuous. We remark that cut Bogner-Fox-Schmit elements for thin plates were proposed and analyzed in [1]. Therefore, the proposed method can be seen as an extension of the work [1] from plates to curved shells.
One challenge in unfitted finite element methods is the efficient integration on the problem domain [36]. During the preparation of the present paper it turned out that due to a gowning number of constraints for finer meshes the strategy 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 convergence analysis where the error is computed with an exact manufactured solution. Furthermore, the capabilities of the method are shown in two standard and one non-standard benchmark tests.

Notation and geometric preliminaries
The underlying assumption in shell analysis is that the computational domain has a small extension with respect to one coordinate. Thus, we assume that it is located around a two-dimensional mid-surface Ω which is embedded in R 3 . In the present paper we assume that the mid-surface is defined implicitly as the zero-level set of a function φ : (1) 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 denoted by µ µ µ, see Figure 2. We assume that Ω is regular such that holds in the neighborhood of Ω, where ∇ denotes the usual gradient of some scalar-valued function f : with the Cartesian coordinates x = (x 1 , x 2 , x 3 ) and the standard Cartesian orthonormal basis {e 1 , e 2 , e 3 }. Here, and in the following, the Einstein summation convention applies. Whenever an index occurs once in an upper position and in a lower position we sum over this index, where Latin indices i, j, . . . take the values 1, 2, 3 whereas Greek indices α, β, . . . take the values 1, 2. Let T be some tensor space of the form R 3 ⊗ · · · ⊗ R 3 . In the following we also use the generalization of the gradient for scalar-valued functions (3) to tensor-valued functions T : R 3 → T , 2.1 Differential geometry of implicitly defined surfaces Given a implicit representation (1) of the surface Ω, we can compute the unit normal vector to the surface by and are able to define the tangential projector, Furthermore, the extended Weingarten map is given by and the mean curvature H is defined as

Differential geometry of parametrized surfaces
We briefly review the differential geometry of parametrized surfaces. For details we refer to e.g. [14]. Although the mid-surface Ω is assumed to be given implicitly, at least a local parametric representation is guaranteed to exist by the implicit function theorem. This justifies to consider parametrizationŝ g : U ⊂ R 2 → Ω with the parameter domain U for theoretical considerations. Given the parametrizationĝ(θ 1 , θ 2 ), we can define the two covariant base vectorsĝ α := ∂ĝ ∂θ α , which span the tangent plane to Ω. With the base vectors we can define the unit normal vector and the covariant coefficients of the metricĜ αβ =ĝ α ·ĝ β . The contravariant coefficients of the metric are given by is the coefficient matrix. The contravariant base vectors can then be computed bŷ g α =Ĝ αβĝ β . The covariant coefficients of the Weingarten mapĤ =ĥ αβĝ α ⊗ g β are given byĥ and obey the symmetry relationĥ αβ =ĥ βα . The mean curvatureĤ is given byĤ =ĥ α α =ĥ αβĜ βα .
Furthermore, the derivatives of the base vectors are given bŷ with the surface Christoffel symbols of the second kind defined bŷ Remark: In our notion a hat over a quantity refers to a dependency on the parametric coordinates (θ 1 , θ 2 ) ∈ U , whereas no hat refers to a dependency on x ∈ R 3 .

Relations between parameter space and embedding space
The fieldû(θ 1 , θ 2 ) defined on the parameter space is related to the field u(x) defined on the embedding space R 3 bŷ By applying the chain rule we find that the first and second derivatives are related byû Furthermore, we have the following relations summarized in the following lemma.
For the Weingarten map we have the relation The proof can be found in Appendix A.

Surface gradient
The surface gradient of a tensor-valued function represented with respect to parametric coordinates by the mapf : U → T is given by Lemma 2 For the representation f : R 3 → T the surface gradient is given by Proof Using the relation between the projector and the metric tensor and (15) we have

Surface divergence
We define the surface divergence as the adjoint operator to the surface gradient [41]. Therefore, on an Riemannian manifold we have in local coordinates where we use the notation detĜ = det([Ĝ αβ ]). The next lemma gives the simpler representation for the surface divergence in case of a surface embedded in R 3 .

Lemma 3
On a surface Ω ⊂ R 3 parametrized byĝ : U → Ω the surface divergence of a tensor-valued function represented byT : U → T is given by Furthermore, for the representation T : 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 second order tensor 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 Then, the following product rules hold The proof can be found in Appendix A.

Integral identities
For further use we introduce the surface divergence theorem for a tensor-valued function T Using (28) and (29), the integration by parts formula for a vector v and a second order tensor T reads

The linear thin shell problem
In this section we derive the governing equations of linear thin shells from first principles of continuum mechanics. Furthermore, we show the equivalence to the linear Koiter model formulated as a minimization problem.

Shell kinematics
The kinematics of the surface Ω is described by the change in metric tensor and the change in curvature tensor. In the present paper we focus on the linear theory and use the linearized change in metric tensorγ γ γ =γ αβĝ α ⊗ĝ β and the linearized change in curvature tensorρ ρ ρ =ρ αβĝ α ⊗ĝ β . The respective covariant components are given by [14,10] The next lemma establishes the representations for γ γ γ =γ γ γ •ĝ −1 and ρ ρ ρ =ρ ρ ρ•ĝ −1 .
Lemma 5 For the linearized change in metric tensor we have the representation and for the linearized change in curvature tensor we have the representation The proof can be found in Appendix A.

Stress and moment tensors
We define the traction vector t(µ µ µ) on a cut defined by the boundary normal µ µ µ tangential to the surface. Due to Cauchys theorem we have the representation with the stress tensor σ σ σ. We decompose the stress tensor σ σ σ in a tangential and a normal part, 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(µ µ µ) on the cut defined by µ µ µ, which can be expressed as with the tangential moment tensor M.

Equilibrium of forces
The equilibrium of forces states that the sum of the resulting force of boundary traction and the resultant force from the surface loading vanishes, Applying the surface divergence theorem (29) results in Due to the fact that Ω is arbitrary the local force equilibrium reads

Equilibrium of moments
The equilibrium of moments states that the sum of boundary moments, the moments of boundary tractions, and the moments due to surface loads vanishes, The following lemma summarizes the consequences of the equilibrium of moments.
The proof can be found in Appendix A.

Constitutive equations
In the present paper we assume linear constitutive equations of the form and The fourth order elasticity tensor E is given by whereλ and µ are the Lamé constants of the elastic material constituting the shell and P s the symmetric part of the tangential fourth order identity tensor. The Lamé constants are related to the Young's modulus E and Poisson's ratio ν byλ The constitutive equations can also be write as We remark that with (45) the condition (42) is fulfilled identically.

Weak form of the governing equations
The weak form of the governing equations is given by where v are appropriate test functions, see Appendix B. The boundary conditions which can be prescribed are given by, where u D i is a given displacement in the direction of e i , ω t and ω µ are given rotations, N N i is a given force in the direction of e i , M N t and M N µ are given moments. 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 normal derivative d µ µ µ (ν ν ν · u) = ∇ Ω (ν ν ν · u) · µ µ µ can be independently prescribed by ω µ .

Equivalence to the Koiter model
In this section the equivalence of the classical linear Koiter model [30,14] and (49) is outlined. To this end, we define the energy functional (51) Then, the linear Koiter shell model reads: Find u ∈ V such that Since the variational equations of (51) are the equations given in (49) we have established the equivalence. Therefore, we conclude that the Koiter model proposed out of purely mechanical and geometrical intuitions can be derived from first principles of continuum mechanics.

C 1 -Trace-Finite-Cell-Method
For the discretization of the weak form (49) we propose a C 1 -continuity version of the TraceFEM. Following the TraceFEM concept the ansatz space on the surface is defined as the restriction (trace) of an outer ansatz space defined on a background mesh. We label the present method also a Finite-Cell method because we use as a background mesh a Cartesian grid. On this structured grid we are able to construct C 1 -continuity shape functions by the tensor product of univariate cubic Hermite form functions. Locally, they are defined on the on the unit interval by (see Figure 3) The local functions (53)  space on the background grid by V h . Then, the discrete problem is given by: holds for all w h ∈ V h and The linear and bilinear forms are given by In (54a) the solution u h gets fixed to prescribed values at the Dirichlet boundary. However, the matrix M h (and also K h ) is singular by definition because of two reason. First, in order to avoid further notation, we take w h ∈ V h resulting in zero rows and columns, which can be easily identified. Secondly, since we define the shape functions by restriction of the shape functions defined on the background mesh, they are not linear independent. In the standard setting of a fitted finite element method the respective degrees of freedom in (54a) are easy to identify and can be determined by interpolation. Here, in the case of an unfitted method a linear independent basis is not known explicitly in general and have to be determine. In (54b) the test functions v h are restricted to the null-space of M h . In the discrete method the integrals in (55) are evaluated by quadrature, which is described in the next section.

Integral evaluation
In order to evaluate the surface and line integrals in (55) we use the quadrature schema developed in [43]. Here, we outline only the main ingredients and refer to [43] for technical details. Following the standard finite element procedure the integrals are evaluated by summing up background cell (face) contributions where the shape functions are smooth. The individual contributions are evaluated by Gaussian quadrature. The main idea from [43] is to subdivide the background cells (faces) until it is possible to convert the implicitly defined geometry into the graph of an implicitly defined height function. Then, a recursive algorithm which requires only one-dimensional root finding and one-dimensional Gaussian quadrature can be set up. In order to chose suitable height function directions we have to be able to ensure the monotonicity of the level-set function in that direction. This can be done by showing that the derivative in that direction is uniform in sign, i.e. place bounds on the values attainable by the derivative. In contrast to [43] we use interval arithmetic for this task.

Solution strategies
In 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 and compute the null-space of M h . We denote the null-space basis by Z h . In a next step the solution u 0 h of is computed. The overall solution is then given by In the penalty method we solve a system of linear equations of the form with the penalty parameter α > 0. In the Lagrange multiplier method we solve the system of linear equation with the Lagrange multipliers λ h .

Implementation
The proposed method has been implemented in Matlab. Within the method the exact level-set function φ(x) is used. For the evaluation of the surface normal vector (5) and the Weingarten map (7), the first and second order derivatives of the level-set function are necessary. In the implementation we use symbolic differentiation of φ(x) to provide these derivatives.
In the present work, we have not used any stabilization term which is added 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 would be to add stabilization terms to the bilinear forms (54b) and (54b). We refer to [34] for an overview of different possibilities. Although such a stabilization can be designed in such a way that the convergence order of the method is not been altered, a stabilization decreases the accuracy of the method. Therefore, a strategy is investigated, where no stabilization is necessary. However, we have observed in numerical experiments that the Matlab backslash operator does not give satisfactory results. Due to this reason, we use the direct solver suitable for under-determined linear equation systems from the SuiteSparse 1 project.

Numerical results
In this section, numerical results are presented. First, we verify the implementation of the proposed method against exact manufactured solutions. Secondly, we demonstrate that the method works by testing it with two benchmark problems (one cylindrical shell and one spherical shell) of the well-known shell obstacle course [7]. Finally, in a fourth example, a complex shell geometry is investigated.

Verification example
The implementation of the proposed method is verified. We have successfully run the method on various geometries, displacement fields and boundary con- and are illustrated in Figure 4. The level-set function (61a) implies a flat . 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 1 (x, y, z) = x 3 e x + y x 3 e y + (xzy 2 + x(x − 1)y(y − 1))e z , (62a) u ex 2 (x, y, z) = sin(16y) cos(16x z)e x + cos(16x y z)e y + 2 sin(16x 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 displacement field u 1 is chosen as a third order polynomial such that the solution can be represented exactly in the discrete space, whereas u 2 can only be approximated.
In the following, we study the behavior of the error under uniform mesh refinement for the three different solution methods given in Section 4.2. Furthermore, for comparison, we also consider the Hermite interpolation u int h of the solution on the background grid and the surface L 2projection: We remark that the Hermite interpolation and the L 2 -projection are only possible if the solution is known, and is thus only computed for the verification examples in this section. Furthermore, all solutions apart from the Hermite interpolation require the solution of a system of linear equations. The numerical results for the flat plate defined in (61a) are visualized in Figures 5 and 6. The refinement level 0 refers to a single background cell. In Figure 5 the errors obtained by the penalty method are given for different refinement levels and penalty factors. We remark that for this flat geometry Errors for different penalty factors on the geometry induced by φ 1 the numerical integration is exact up to round of errors. Therefore, for u 1 the sources for errors are the error due to the imposition of the boundary conditions by the penalty method (depending on the penalty factor) and round off errors in the numerical computations. As expected the errors decrease with increasing penalty parameter up to a point where the ill-conditioning of the linear system dominates the error. Since u 2 can not be represented exactly in the discrete space an approximation error limits the overall error. In Figure 6 the results for the different solution methods are given. For the penalty method we used the lowest errors of the results shown in Figure 5. As no system of equation has to be solved for the interpolation on the background grid (volume interpolation) the error is around 10 −16 for u 1 for all refinement levels. In contrast to this the other results require the solution of a system of equations and therefore the errors are between 10 −8 and 10 −4 due to round off errors. For u 2 we observe the convergence of all methods with optimal rate. Here, by definition the L 2 -projection gives the lowest error, whereas the Hermite interpolation results in the highest error for a fixed refinement level (apart from level 5, where the error due to ill-conditioning of the system of equations dominates). The numerical results for the part of the ellipse defined in (61b) are visualized in Figures 7 and 8. In Figure 7 the results of the penalty method for different penalty parameters are given. In contrast to the flat geometry, now the numerical integration is not exact, yielding an additional error, which dominates for the three coarsest levels.
In Figure 8 the errors for the three different solution methods, the Hermite interpolation and the L 2 -projection are visualized. Again, for the penalty method we used the lowest errors of the results shown in Figure 7. Similar as before, for the displacement field u 1 the volume interpolation gives errors around 10 −16 , whereas for the other solutions the errors are between 10 −7 and 10 −2 due to round off errors. For displacement field u 2 we observe the convergence of all methods. However, due to round off errors the accuracy is limited. Nevertheless, we remark that an relative error level of about 10 −5 is more than sufficient for practical problems. This can also be seen from the visualizations of the solutions obtained with the null-space method for u 1 in Figure 9 and for u 2 in Figure 10. For the fine levels no difference in the solutions can be seen by eye.

Scordelis-Lo roof
We consider the classical Scordelis-Lo roof problem, which is one example from the shell obstacle course [7]. It is a popular benchmark test to assess the performance of finite elements regarding complex membrane strain states. The cylindrical roof (radius r = 25) is supported by rigid diaphragms at the ends (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. The structure is subjected to gravity loading with b = −90 e z . We describe the problem geometry by φ(x, y, z) = y 2 + z 2 − r 2 , and B = [0, 50] × [−r sin( 40 180 π), r sin( 40 180 π)] × [10, 31.25]. We study the vertical displacement of point A, which is located in the middle of one free edge. As a reference solution we use the overkill solution u A z = −0.3006 from [8] obtained by an isogeometric formulation using fifthorder NURBS and a mesh of 48 control points in each direction. The results for different meshes obtained with the presented methods are given in Table 1and the deformed geometry is depicteed in Figure 12. We observe that the nullspace method and the penalty method are able to reproduce the reference displacement found in literature accurately. However, the results obtained by the Lagrange multiplier method show some instability of the method. The investigation of the origin of these instabilities is topic of further research.

Pinched hemisphere
In this example, we consider the pinched hemisphere problem from the shell obstacle course in [7]. We describe the spherical mid-surface by  of the hemisphere is unconstrained and the four radial forces have alternating signs such that the sum of the applied forces is zero. We investigate the radial displacement at the loaded points. In [7], the reference displacement of u r = 0.0924 is given. The results obtained by the presented methods are given in Table 2 and the deformed configuration is depicted in Figure 14. We observe that here all three methods give nearly the same results. The values obtained for the finest levels are in very good agreement with the reference value found in literature.

Gyroid
In this example, we consider the deformation of a shell structure with a complex geometry. The mid-surface is part of a gyroid which is given by the levelset function φ(x, y, z) = sin(πx) cos(πy) + sin(πy) cos(πz) + sin(πz) cos(πx). The geometry and the material parameters are depicted in Figure 15. The shell structure is clamped at the boundary curve which is in the plane x = 0. We assume a thickness t = 0.03. We study the vertical deflection due to a volume load b = −10 7 e z at the point [2, 0.5, −0.25]. The deformed geometry is depicted in Figure 16. The results of the proposed methods are summarized in Table 3. We observe that the results obtained by the null-space method and the penalty method are nearly the same and that they are in good agreement with the reference displacement u z = −1.8812 given in [26]. We remark that the reference solution was obtained for a seven-parameter shell model including more physical effects and thus leading to a slightly larger displacement. Therefore, the deviation in the deflection is acceptable. However, the results obtained by the Lagrange multiplier method are incorrect. This issue should be further investigated in future work.

Conclusions
We have developed a C 1 -continuous finite element method for thin shells with mid-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 surface are obtained by restriction of tensor-product cubic Hermite splines on a structured background mesh. In order to allow a natural implementation, the underlying shell model is formulated in a parametrization-free way. Furthermore, the strong form of the governing equations are given. This allows to obtain manufactured solutions on arbitrary geometries. Thus, the implementation of the proposed method is verified by a convergence analysis where the error is computed with an exact manufactured solution.
In the present method, the shape functions on the shell surface are linearly dependent. In order to avoid a singular system matrix, a stabilization term can be used. In the presented method such a stabilization is avoided. However, it is necessary to use the direct solver suitable for under-determined linear equation systems from the SuiteSparse project. We investigated three strategies to include the boundary conditions. These are the penalty method, the Lagrange multiplier method, and the null-space method. In the numerical experiments we have observed that the penalty method and the null-space method give reliable results. However, the Lagrange multiplier method suffers from instabilities in some examples, which should be further investigated. In future work, it would be also interesting to use iterative solvers in contrast to the used direct solver.
In contrast to thin shells, for Reissner-Mindlin shells only C 0 -continuous shape functions are commonly used. In order to avoid transverse shear locking, in [31,22] an hierarchic concept of shell models is presented. This approach has the advantage that transverse shear locking is eliminated on the continuous formulation level, independent of a particular discretization, but requires C 1continuous shape functions. An extension of the present work to Mindlin-Reissner shells with implicitly defined mid-surface seems possible and would be worth to investigate.

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 0 and integrate over the shell surface, Here, the function space of the test functions is V 0 = {η η η : Ω → R 3 | γ γ γ(η η η) ∈ L 2 (Ω, R 3 ), ρ ρ ρ(η η η) ∈ L 2 (Ω, R 3 ), η η η · e i = 0 on Γ D i , ∇ Ω (ν ν ν · η η η) · µ µ µ = 0 on Γ Dc }, where Γ D i , Γ D t , and Γ Dµ denote Dirichlet boundaries. On Γ D i the displacement in direction e i is restrained and on Γ D t and Γ Dµ the rotation of the shell around the boundary tangent vector and the boundary normal vector is restrained respectively. The corresponding Neumann 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 We have σ σ σ = N + S ⊗ ν ν ν and obtain Using (43) and integration by parts of the second term on the left yields Due to (48) we obtain the relation