Cut Bogner-Fox-Schmit Elements for Plates

We present and analyze a method for thin plates based on cut Bogner-Fox-Schmit elements, which are $C^1$ elements obtained by taking tensor products of Hermite splines. The formulation is based on Nitsche's method for weak enforcement of essential boundary conditions together with addition of certain stabilization terms that enable us to establish coercivity and stability of the resulting system of linear equations. We also take geometric approximation of the boundary into account and we focus our presentation on the simply supported boundary conditions which is the most sensitive case for geometric approximation of the boundary.


Introduction
The Bogner-Fox-Schmit (BFS) element [6] is a classical C 1 thin plate element obtained by taking tensor products of cubic Hermite splines and removing the interior degrees of freedom that are zero on the boundary. In this paper, we consider a variant where we retain these degrees of freedom to obtain a C 1 version of the Q 3 approximation [23]. This element is only C 1 on tensor product (rectangular) elements, which is a serious drawback since it severely limits the applicability of the resulting finite element method. However, on geometries allowing for tensor product discretization it is generally considered to be one of the most efficient elements for plate analysis, cf. [24, p. 153]. It is also a reasonably low order element for plates which is very simple to implement, in contrast with triangular elements which either use higher order polynomials, such as the Argyris element [1], or macro element techniques, such as the Clough-Tocher element [10]. The construction of curved versions of these elements for boundary fitting can also be cumbersome, see, e.g., [5]. It should be noted that the use of straight line segments for discretizing the boundary is not to be recommended, not only because of accuracy issues but also due to Babuška's paradox for simply supported plates, see [3].
An alternative to C 1 approximations for Kirchhoff plates is to either use discontinuous Galerkin methods [8,15,17,22], or to use mixed finite elements for the Reissner-Mindlin model with small plate thickness [2,4,12]. These C 0 methods alleviate the problem of boundary approximation. In this paper we present an alternative idea where C 1 continuity is retained: we develop a cut finite element version, allowing for discretizing a smooth boundary which may cut through the tensor product mesh in an arbitrary manner. Adding stabilization terms on the faces associated with elements that intersect the boundary, we obtain a stable method with optimal order convergence. We prove a priori error estimates which also take approximation of the boundary into account. The focus of the analysis is on simply supported boundary conditions, the computationally most challenging case.
The paper is organized as follows. In "The Kirchhoff plate" section, we recall the thin plate Kirchhoff model; in "The finite element method" section we formulate the cut finite element method; in "Error estimates" section we present the analysis of the method starting with a sequence of technical results leading up to a Strang Lemma and an estimate of the consistency error and finally a priori error estimates in the energy and L 2 norms. In "Numerics" section, we present some numerical illustrations, and in "Concluding remarks" section some concluding remarks are included.

The Kirchhoff plate
Consider a simply supported thin plate in a domain ⊂ R 2 with smooth boundary ∂ .
where σ (∇v) is the stress tensor with (∇v) the strain tensor and κ the parameter with E the Young's modulus, ν the Poisson's ratio, and t the plate thickness. Since 0 ≤ ν ≤ 0.5 both κ and ν(1 − ν) −1 are uniformly bounded.
We shall focus our presentation on simply supported boundary conditions where the moment tensor M is defined by and M ab = a · M · b for a, b ∈ R 2 . Other conditions, such as clamped boundaries, can be handled using the same techniques as in the following, cf. [15]. The weak form of (1) and (5) takes the form: where and l(v) = (f, v) . The form a is symmetric, continuous, and coercive on V equipped with the H 2 ( ) norm and it follows from the Lax-Milgram theorem that there exists a unique solution in V to (7). Furthermore, for smooth boundary and f ∈ L 2 we have the elliptic regularity The finite element method

The mesh and finite element space
We begin by introducing the following notation.
• Let T h , h ∈ (0, h 0 ], be a family of partitions of R 2 into squares with side h. Let V h be the Bogner-Fox-Schmit space consisting of tensor products of cubic Hermite splines on T h . Note that V h | T = Q 3 (T ), with Q 3 (T ) the tensor product P 3 (I 1 ) ⊗ P 3 (I 2 ) of cubic polynomials where T = I 1 × I 2 ⊂ R 2 . • Let ρ be the signed distance function, positive on the outside and negative on the inside, associated with ∂ and let U δ (∂ ) = {x ∈ R 2 : |ρ(x)| < δ} be the tubular neighborhood of ∂ of thickness 2δ. Then there is δ 0 > 0 such that the closest point mapping p : Furthermore, we assume that for each element T such that ∂ h intersects the interior of T , i.e. int(T ) ∩ ∂ h = ∅, the curve segment ∂ h ∩ T is smooth and intersect the boundary ∂T of T in precisely two different points. Let X h be the set of all points where ∂ h is not smooth and note that the number of elements We can always satisfy this assumption by enlarging the active mesh T h if necessary.
We illustrate some of these quantities in Fig. 1 consists of all element sides on the grey elements excluding those without neighbouring elements.

The finite element method
The method reads: find u h ∈ V h such that The forms are defined by where T = (M · ∇) n + ∇ t M nt (18) with sub-indices n and t indicating scalar product with the normal and tangent of ∂ h , and β, γ are positive parameters which are proportional to κ. Here s h is a stabilization form involving jumps of the second and third normal derivatives in the direction n F , the normal to the element face, with which provides necessary stability at the cut elements, see (23). The bilinear form, apart from the stabilization terms, stems from Nitsche's method [20], first analyzed for plates in a discontinuous Galerkin setting in [15].

The energy norm
Define the following energy norm on where and we employ the standard notation v 2 where ∇ j v is the tensor of all j:th order derivatives.

Stabilization
The stabilization term provides us with the following bound which follows from the standard estimate where T 1 and T 2 are elements that share the face F , and v| T i ∈ P p (T i ), the space of polynomials of order p. See for instance [16,19] for further details.

Continuity and coercivity
The form A h is continuous which follows directly from the Cauchy-Schwarz inequality, and for γ large enough coercive

Verification of (25)
We first recall the cut trace inequality which combined with (26) give the cut inverse trace inequality Now, using the inverse trace inequality (28), the inverse inequality (27), followed by the stabilization estimate (22) we obtain and thus there is a constant C * such that As a consequence, |||v||| h and |||v||| h, * , where are equivalent norms on V h . We then have and we find that taking δ small enough to guarantee that

Poincaré inequality
We have the following Poincaré inequality

Verification of (36)
Using the stabilization estimate (22) and the fact that T h,I is covered by h we have Next to estimate v 2 O h we again use the stabilization estimate (22) and the fact that T h,I is covered by , Let P 1, : L 2 ( ) → P 1 ( ) be the L 2 projection onto the space of linear functions on . Then for v ∈ H 2 ( ), and in particular and using the trace inequality Note that the constants are independent of the mesh parameter since is fixed. We then have v 2 which together with (37) proves (36). Here we estimated term i using the stabilization (22), Finally, to estimate ii we recall the following technical estimate, which we prove in the "Appendix" of this paper see also the appendix of [9], where δ ∼ h 4 , bounds the distance between ∂ and ∂ h , see (10). Note that the finite element functions are defined on O h , that contains both and h , and (40), the stabilization estimate (22), we obtain where we finally used (44) to estimate ∇ 2 v and the fact that This completes the verification.

Interpolation
Let To construct an interpolation operator for cut elements we recall that given v ∈ H s ( ) there is an extension operator E : for all s > 0, cf. [21]. For simplicity we will often use the notation w = Ew for the extension of w ∈ H s ( ) to R 2 . We define the interpolation operator Combining (52) with (53) we obtain the interpolation error estimate For the energy norm we have the estimate

Verification of (56)
Let η = v − π h v and recall that The first term is directly estimated using the interpolation error estimate (55), For the second term we employ the trace inequality to conclude that For the third term we use the cut trace inequality (26) and the interpolation estimate (55), where T h (F h,B ) ⊂ T h is the set of elements with a face in F h,B . Finally, the fourth term is estimated in the same way as the third, which completes the verification of (56).

Consistency error estimate
Lemma 1 Let u be the exact solution to (1) with boundary conditions (5), and u h the finite element approximation defined by (13), then Proof Adding and subtracting an interpolant we obtain Using coercivity we can estimate the second term on the right hand side as follows where we added and subtracted u in the numerator and for the first term used the estimate A h (π h u−u, v) |||π h u−u||| h |||v||| h and for the second used (13) to eliminate u h . Combining the estimates the desired result follows directly.

Lemma 2
Let ϕ ∈ H 4 (R 2 ) and v ∈ V + V h , then with n ± h and t ± h the left and right limits to tangent and normal to the discrete boundary Proof Using the simplified notation M = M(ϕ) and T = T (ϕ) for brevity we obtain by integrating by parts Splitting ∇v in tangent and normal contributions on ∂ h , we have the identity where we integrated by parts along the curve segments ∂ h ∩ T , and ν is the exterior unit tangent vector to ∂ h ∩ T . Summing over all elements that intersect ∂ h , we obtain the identity Combining (75) and (78), we obtain and setting T = (M · ∇) n + ∇ t M nt we obtain the desired result.

Lemma 3 Let u be the exact solution to (1) with boundary conditions (5), then there is a constant such that for all v ∈ V h ,
where |||v||| h, is the norm Proof Using the definition (14), the fact that s h (u, v) = 0 for u ∈ H 4 ( ), and the partial integration identity (71) we obtain Before turning to the estimates of I −IV we first note that for w ∈ H 1 0 ( ), with extension to R 2 also denoted by w, we may apply (45) and the stability (53) of the extension that with δ ∼ h 4 . Note that here we do not need to restrict the norms to O h as in (45) since the extended function is defined on R 2 . See [9, Appendix] for detailed derivations. Furthermore, for more regular functions such that w ∈ H 2 ( ) with w = 0 on ∂ , we may strengthen the estimate as follows where ∂ t = {x ∈ R 2 |ρ(x) = t} for t ∈ (−δ 0 , δ 0 ), are level sets of ρ, and again δ ∼ h 4 . In the last step we used a version of (45) to conclude that where we used a trace inequality on and the stability (53) of the extension operator.
I. Using Cauchy-Schwarz followed by (89) with w = M nn (u), we get Here we used the estimate which we derive by applying (45) with w = ∇v, where we used the trace inequality ∇v ∂ ∇v H 1 ( ) v H 2 ( ) , and at last the Poincaré inequality (36).
II. Using the assumption on the accuracy of the discrete normal (11) we have for each x ∈ X h , where the first term on the right hand side can be estimated as follows We then have where we used the fact that the number of elements, denoted by |X h |, in X h satisfies |X h | ∼ h −1 , and the Sobolev inequality [13] followed by the stability (53) of the extension operator to obtain and the estimate To verify (100) we shall employ an inverse inequality locally using a linear approximation of the boundary. To that end consider x ∈ X h and let B r (x) be a ball of radius r ∼ h centred at x. Let T x ∈ T h be one of the elements such that x ∈ ∂T x ∩ ∂ h and given v ∈ V h let v x ∈ Q 3 (R 2 ) be the extension to Such a line exists since ∂ is smooth and linear approximation is of second order locally. We then have where we used the fact that | x ∩ B r (x)| ∼ h. In order to employ a local version of (45) we define the cylindrical tubular neighborhood then we added and subtracted v in the first term, used the triangle inequality, the inverse estimate ∇v x 2 U δ 1 ,x ∇v x T x which holds since v x is a polynomial, and U δ 1 ,x ⊂ B r (x) for some ball of radius r ∼ h, and finally we noted that v x = v on T x . Inequality (103) is an application of (154). Inserting (105) in (101) we get which establishes (100). Here we used the fact that the number of cylinders U δ 1 ,y , y ∈ X h , that intersect U δ 1 ,x is uniformly bounded independent of x ∈ X h and h ∈ (0, h 0 ]. We also used certain estimates of terms i and ii, which we verify next. which is a local version of (22) on the patch T h,x , we get . ii. We first note that where N h (T h (∂ h )) is the set of elements that share a node with an element in T h (∂ h ), and we used the stabilization. We can then choose δ 2 ∼ h such that We now proceed in a similar way as in (89), where we used the estimate (90) in the last step. Inserting (114) in (112) and using that δ 1 ∼ h 2 and δ 2 ∼ h, we obtain III. Using (89) with w = u and recalling δ ∼ h 4 , IV. Proceeding in the same way as for Term III, Combining the estimates we find that which completes the proof.

Theorem 1 The finite element solution defined by (13) satisfies
Proof Using the second bound of (81) in (66) followed by the interpolation error bound (56) we directly get the desired estimate.

Theorem 2
The finite element solution defined by (13) satisfies Proof Adding and subtracting an interpolant and using the interpolation error estimate (55) we have the estimate In order to estimate π h u − u h h we let φ h ∈ V h be the solution to the discrete dual problem Setting v = π h u − u h we obtain the error representation Here φ is the solution to the continuous dual problem extended to R 2 using the stable extension operator, see (53).
I. Since φ h is the finite element approximation of φ we have the error estimate where we used elliptic regularity (9), which directly gives II. Using the fact that s h (π h u − u, φ) = 0 since φ ∈ H 4 ( ), the partial integration formula (71), the Cauchy-Schwarz inequality, and the interpolation error estimates we obtain Here we used (88) with δ ∼ h 4 followed by the elliptic regularity (9) to conclude that and using (89) we obtain III. Using (81) we obtain We used the estimate where we used (89) to conclude that T (φ) 2 Collecting the estimates of Terms I-III completes the proof.

Implementation
We consider two higher order approximations of the boundary: a piecewise cubic C 0 approximation or a piecewise cubic C 1 approximation. The steps to create the approximate boundary are as follows. Note that the approximation of the boundary may partly land outside the element. In such cases, the basis functions of the element containing the straight segment is used also outside of the element.

Example
We consider a circular simply supported plate under uniform load p. The plate is of radius R = 1/2 and has its center at x = 1/2, y = 1/2. Defining r as the distance from the midpoint we then have the exact solution see, e.g., [18]. The constitutive parameters were chosen as E = 10 2 , ν = 0.3, t = 10 −1 , and the stabilization parameters as β = 10 −2 , γ = 10 2 (2κ + 2κν(1 − ν) −1 ). We compare the convergence in normalized (||u − u h ||/||u||) L 2 , H 1 and piecewise H 2 norms in Fig. 6. These norms are computed on the discrete geometry, for simplicity the straight segment geometry. The solid lines indicate second, third, and fourth order convergence, respectively from top to bottom, and we note that we observe a slightly higher than optimal rate of convercence of about O(h 1/2 ) in all norms. We note that the continuity of the approximation of the boundary seems not to be crucial as the convergence curves are very close. In Fig. 7 we show an elevation of the solution on one of the meshes in the sequence. We also show the influence of the parameter β on a fixed mesh (coarse, 386 active nodes). In Fig. 8 we show the condition number as β increases from its critical number, the number for which the system matrix is singular. For lower values of β there are negative eigenvalues in the system matrix. We note that as β increases, the condition number will eventually increase again after an initial drop. In Fig. 9 we show the H 1 error which is more sensitive to the increase in β. We remark that this effect, however, does not affect the convergence rate.

Concluding remarks
We have proposed and analyzed a cut finite element method for a rectangular plate element, allowing for curved boundaries. The analysis shows that the method is optimally order convergent and stable. Two different approximations of the boundary have been tested, a standard cubic interpolation of the exact boundary and a cubic spline approximation leading to a continuously differentiable approximation of the boundary. Numerical results and theory indicate that the continuity of the boundary approximation is not

Appendix: Some inequalities
Let ω ⊂ ∂ and define the cylindrical tubular neighborhood U δ (ω) = {x ∈ U δ (∂ ) : p(x) ∈ ω} δ ∈ (0, δ 0 ] (147) Let 1 and 2 be two surface segments in U δ (ω) with unit normals n i , such that the closest point mapping p : i → ω is a bijection with inverse p −1 i : ω → i and there is constant such that 1 min where n(x) = n • p(x). We can then define a bijection q : 1 → 2 by q = p −1 2 p. For each x ∈ 1 let I x be the line segment with endpoints x ∈ 1 and q(x) ∈ 2 . We then have where we integrate along I from q(x) to x. Using the Cauchy Schwarz inequality we get v 2 (x) v 2 (q(x)) + Here S = ∪ x∈ 1 I x is the domain between the surfaces 1 and 2 , and we changed coordinates to integration over 2 and S equipped with Euclidian measure. We conclude that v 2 In applications, it is of convenient to simply replace S by the larger domain U δ (ω). Typical applications include taking 2 = ω = ∂ and 1 = ∂ h , which gives v 2 For For v ∈ H 1 ( ), with extension to R 2 also denoted by v, we have S ⊂ U δ (∂ ) and we get v 2 In this paper we have recall that ∂ h ⊂ U δ (∂ ), with δ ∼ h 4 , see (10).