A new family of efficient conforming mixed finite elements on both rectangular and cuboid meshes for linear elasticity in the symmetric formulation
aa r X i v : . [ m a t h . NA ] J a n A NEW FAMILY OF EFFICIENT CONFORMING MIXED FINITE ELEMENTS ONBOTH RECTANGULAR AND CUBOID MESHES FOR LINEAR ELASTICITY INTHE SYMMETRIC FORMULATION
JUN HU
Abstract.
A new family of mixed finite elements is proposed for solving the classical Hellinger–Reissnermixed problem of the elasticity equations. For two dimensions, the normal stress of the matrix-valued stressfield is approximated by an enriched Brezzi–Douglas–Fortin–Marini element of order k , and the shear stressby the serendipity element of order k , the displacement field by an enriched discontinuous vector-valued P k − element. The degrees of freedom on each element of the lowest order element, which is of first order,is 10 plus 4. For three dimensions, the normal stress is approximated by an enriched Raviart–Thomaselement of order k , and each component of the shear stress by a product space of the serendipity elementspace of two variables and the space of polynomials of degree ≤ k − Q k − element. The degrees of freedom oneach element of the lowest order element, which is of first order, is 21 plus 6. A family of reduced elementsis also proposed by dropping some interior bubble functions of the stress and employing the discontinuousvector-valued P k − (resp. Q k − ) element for the displacement field on each element. As a result the lowestorder elements have 8 plus 2 and 18 plus 3 degrees of freedom on each element for two and three dimensions,respectively.The well-posedness condition and the optimal a priori error estimate are proved for this family of finiteelements. Numerical tests are presented to confirm the theoretical results. Introduction
The first order system of equations, for the symmetric stress field σ ∈ Σ := H (div , Ω , S ) and the dis-placement field u ∈ V := L (Ω , R n ), reads: Given f ∈ L (Ω , R n ) find ( σ, u ) ∈ Σ × V such that( Aσ, τ ) L (Ω) + (div τ, u ) L (Ω) = 0 , (div σ, v ) L (Ω) = ( f, v ) L (Ω) , (1.1)for any ( τ, v ) ∈ Σ × V . Here and throughout this paper, the compliance tensor A ( x ) : S → S is bounded andsymmetric positive definite uniformly for x ∈ Ω with S := R n × n sym the set of symmetric tensors. The space H (div , Ω , S ) is defined by H (div , Ω , S ) := { τ ∈ L (Ω , S ) and , div τ ∈ L (Ω , R n ) } equipped with the norm k τ k H (div , Ω) := k τ k L (Ω) + k div τ k L (Ω) . The stress-displacement formulation within the Hellinger-Reissner principle for the linear elasticity is onecelebrated example of (1.1).
Date : November 6, 2018.
Key words and phrases.
Mixed method; enriched Brezzi–Douglas–Fortin–Marini element; enriched Raviart–Thomas ele-ment; serendipity elementAMS Subject Classification: 65N30, 65N15, 35J25.The author was supported by the NSFC Projects 11271035, 91430213 and 11421101.
Compared with the mixed formulation of the Poisson equation, see for instance, [16], there is an ad-ditional symmetric requirement on the stress tensor. Such a constraint makes the stable discretization ofthe piecewise polynomials extremely difficult. Then one idea that may be come up with is to enforce thesymmetry condition weakly, which in fact leads to Lagrange multiplier methods [2, 6, 13, 33, 34, 35, 36]. Asan alternative method, composite elements were proposed by Johnson and Mercier [31], and Arnold, Dou-glas Jr., Gupta, [7]. That idea might be motivated by the Hsieh-Clough-Tocher element for the biharmonicproblem [21]. Indeed, there is an observation in [31] that the discrete divergence free space therein is therange of the Airy stress function of the Hsieh-Clough- Tocher plate element space, see a similar observationin [7]. Given a scalar field q , the Airy stress function reads Jq := ∂ q∂y − ∂ q∂x∂y − ∂ q∂x∂y ∂ q∂x ! . Unfortunately, this observation was not further explored until more than twenty years later its importancewas realized by Arnold and Winther [9]. In that landmark paper, it was found that to design a stablediscrete scheme is to look for a discrete differential complex with the commuting diagram which reads, fortwo dimensions,0 −−−−→ P (Ω) ⊂ −−−−→ C ∞ (Ω) J −−−−→ C ∞ (Ω , S ) div −−−−→ C ∞ (Ω , R ) −−−−→ y id y I h y Π h y P h −−−−→ P (Ω) ⊂ −−−−→ Q h J h −−−−→ Σ h div h −−−−→ V h −−−−→ Q h is some conforming or nonconforming finite element space for the biharmonic equation; J h anddiv h are the discrete counterparts of the Airy operator J and the divergence operator div, respectively, withrespect to some regular triangulation T h of Ω; Σ h and V h are some finite element approximations of Σ and V , respectively; I h and Π h are canonical interpolation operators for the spaces Q h and Σ h , respectively; P h is the L projection operator from V onto V h . In particular, this commuting diagram implies the FortinLemma [16]. See, Arnold, Awanou, and Winther [5] for the corresponding theory in three dimensions. Basedon those fundamental theories, conforming mixed finite elements of piecewise polynomials on both simplicialand product meshes can then be developed for both 2D and 3D [1, 3, 5, 9]; see [17, 18] for the implementationof the lowest order method of [9]. To avoid complexity of conforming mixed elements, several remedies areproposed, see, [8, 22, 25, 26] for new weak-symmetry finite elements, [10, 24, 30, 32, 38] for non-conformingfinite elements. See also [11, 19] for the enrichment of nonconforming elements of [30, 32] to conformingelements. In a recent paper [29], a family of first order nonconforming mixed finite elements on productmeshes is proposed for the first order system of equations in any dimension, which was extended to a familyof conforming mixed elements in [28].This paper presents a family of conforming mixed elements for both two and three dimensions ( n = 2 , H (div)-conformity of the normal stress and the H -conformity of two corresponding variables for eachcomponent of the shear stress; see also [12] and [28, 29] for a similar observation in two dimensions. Fortwo dimensions, in these elements, an enriched Brezzi-Douglas-Fortin-Marini (hereafter BDFM) element oforder k is proposed to approximate the normal stress, the serendipity element of order k [4, 14, 21] is used toapproximate the shear stress. This discrete space for the stress and an enriched discontinuous P k − elementfor the displacement space are able to form a stable discretization of the two dimensional problem underconsideration. In the first order method which is the two dimensional element of [28], of this family, the ONFORMING MIXED FINITE ELEMENTS 3 total degrees of freedom is, | E | + 6 | K | + | P | , with | E | the number of edges, and | K | the number of elements,and | P | the number of vertices of the partition T h . Note that the total degrees of freedom of the first orderconforming mixed element method on rectangular meshes in [19] is, 3 | E | + 6 | K | + | P | . For three dimensions,an enriched Raviart–Thomas element of order k is constructed to approximate the normal stress, and eachcomponent of the shear stress is approximated by a product space of the serendipity element of order k withrespect to two associated variables and the P k − element with respect to the rest variable. An enriched Q k − element space is taken as the space for the displacement. In the first order method which is the threedimensional element of [28], of this family, the total degrees of freedom is, | E | + | F | + 9 | K | , with | E | thenumber of edges, | F | the number of faces, and | K | the number of elements, of the partition T h . Note thatthe total degrees of freedom of the first order conforming mixed element method on cuboid meshes in [11]is, 2 | E | + 8 | F | + 18 | K | . A family of reduced elements is also proposed by dropping interior bubble functionson each element. As a result the lowest order elements have 8 plus 2 and 18 plus 3 degrees of freedom oneach element for two and three dimensions, respectively, which were announced independently in [20] afterthe first version of this paper was submitted.These spaces of this paper are perfectly and tightly matched on each element. However, the analysisof the discrete inf-sup conditions for these elements has to overcome the difficulty of not using directlythe Fortin Lemma, the key ingredient for the stability analysis of the mixed finite element method for theelasticity problem, see, for instance, [1, 3, 5, 9]. For pure displacement boundary problem, the remedy is anexplicitly constructive proof of the discrete inf-sup condition, which can be regarded as a generalization tothe more general case of the idea due to [29]; see also [12] and [28]. For the more general case, in particularthe pure traction boundary problem, we prove that the divergence space of the H (div) bubble functionspace is identical to the orthogonal complement space of the rigid motion space with respect to the discretedisplacement space on each macro-element. As we shall see in Section 4, the proof for such a result is verydifficult and complicated. One important technique is to use two classes of orthogonal polynomials, namely,the Jacobi polynomials and the Legendre polynomials. As a second step, we construct a quasi–interpolationoperator to control macroelementwise rigid motion for k >
1. Then the discrete inf–sup condition follows.For the first order methods with k = 1, we succeed in proposing a new macroelement technique to finallyestablish the discrete inf–sup condition, which can be regarded as an extension to the more general case ofthat from [34, 37].This paper is organized as follows. In the following two sections, we present the new mixed elements fortwo dimensions and analyze their properties including the well-posedness. In section 4, we consider the puretraction boundary problem and prove the well-posedness of the discrete problem. In section 5 we define thenew mixed elements for three dimensions. In section 6, we present a family of reduced elements by droppingsome interior bubble functions on each element. In section 7 we briefly summarize the error estimates of thediscrete solutions and present two numerical examples, one for the pure displacement boundary problem,and the other for the pure traction boundary problem.2. Mixed finite element approximation in two dimensions
For approximating Problem (1.1) by the finite element method, we introduce a rectangular triangulation T h of the rectangular domain Ω ⊂ R such that S K ∈T h K = ¯Ω, two distinct elements K and K ′ in T h areeither disjoint, or share the common edge e , or a common vertex. Let E denote the set of all edges in T h with E K,V the two vertical edges of K ∈ T h and E K,H the two horizontal edges. Let E IV and E IH denote thesets of all the interior vertical and horizontal edges of T h , respectively, and V I be the set of all the internalvertices of T h . Given vertex A ∈ V I , let E ( A ) be the set of edges that take A as one of their endpoints. JUN HU
Given any edge e ∈ E we assign one fixed unit normal ν with ( ν , ν ) its components, also let t = ( − ν , ν )denote the tangential vector.For each K ∈ T h , we introduce the following affine invertible transformation F K : ˆ K → K, x = h x,K ξ + x ,K , y = h y,K η + y ,K with the center ( x ,K , y ,K ), the horizontal and vertical edge lengthes h x,K and h y,K , respectively, and thereference element ˆ K = [ − , . Given any integer k , let P k ( ω ) denote the space of polynomials over ω oftotal degrees not greater than k , let Q k ( ω ) denote the space of polynomials of degree not greater than k in each variable. Let P k ( X ) be the space of polynomials of degree not greater than k with respect to thevariable X , and P k ( X, Y ) be the space of polynomials of degree not greater than k with respect to thevariables X and Y .For the symmetric fields σ = σ σ σ σ ! ∈ S , we refer to σ n := ( σ , σ ) T as the normal stress and σ as the shear stress.Before defining the space for the stress, we introduce new mixed finite elements for the second orderPoisson equation and the serendipity element of [4, 14, 21]. Given K ∈ T h and an integer k ≥
1, the newmixed finite element space of order k for the second order Poisson equation reads : H k ( K ) := ( P k ( K )) \ span { (0 , x k ) T , ( y k , T } ⊕ E k ( K ) , where E k ( K ) := span { ( x k +1 , T , (0 , y k +1 ) T , ( x y k − , T , (0 , y x k − ) T } . To define the degrees of freedom of the space H k ( K ), we introduce the well–known Jacobi polynomials:(2.1) J ℓ ( ξ ) := (( ℓ + 1)!) ℓ X s =0 s !( ℓ + 1 − s )!( s + 1)!( ℓ − s )! (cid:18) ξ − (cid:19) ℓ − s (cid:18) ξ + 12 (cid:19) s , for any ξ ∈ [ − , Z − (1 − ξ ) J l ( ξ ) J m ( ξ ) dξ = 82 l + 3 (( l + 2)!) ( l + 3)! l ! δ lm . We also need the Legendre polynomials L ℓ ( ξ ) := 12 ℓ ℓ ! d ℓ ( ξ − ℓ dξ ℓ for any ξ ∈ [ − , . The Legendre polynomials satisfy the orthogonality condition:(2.3) Z − L l ( ξ ) L m ( ξ ) dξ = 22 l + 1 δ lm . Lemma 2.1.
The vector-valued function (ˆ q , ˆ q ) T =: ˆ q ∈ H k ( ˆ K ) can be uniquely determined by the followingconditions: (1) R ˆ e ˆ q · ˆ ν ˆ pd ˆ s for any ˆ p ∈ P k − (ˆ e ) and any ˆ e ⊂ ∂ ˆ K , (2) R ˆ K ˆ q J k − ( ξ ) dξdη , and R ˆ K ˆ q J k − ( η ) dξdη , (3) R ˆ K ˆ q L k − ( η ) dξdη , and R ˆ K ˆ q L k − ( ξ ) dξdη , (4) R ˆ K ˆ q · ˆ pdξdη for any ˆ p ∈ ( P k − ( ˆ K )) . ONFORMING MIXED FINITE ELEMENTS 5
Proof.
Since the dimension of the space H k ( ˆ K ) is equal to the number of these conditions, it suffices toprove that ˆ q ≡ q · ˆ ν ∈ P k − (ˆ e ), the first condition (1) implies thatˆ q = (1 − ξ )(ˆ g + c J k − ( ξ ) + b L k − ( η )) , and ˆ q = (1 − η )(ˆ g + c J k − ( η ) + b L k − ( ξ )) , where ˆ g , ˆ g ∈ P k − ( ˆ K ), and c , c , b , b are four interpolation parameters, and J k − and L k − are theJacobi and Legendre polynomials of degree k −
1, respectively. We first consider the case k ≥
2. It followsfrom (2.2) that Z ˆ K (1 − ξ )(ˆ g + b L k − ( η )) J k − ( ξ ) dξdη = 0and Z ˆ K (1 − η )(ˆ g + b L k − ( ξ )) J k − ( η ) dξdη = 0 . Therefore, by the condition (2), c = c = 0 . The condition (2.3) implies Z ˆ K (1 − ξ )ˆ g L k − ( η ) dξdη = 0and Z ˆ K (1 − η )ˆ g L k − ( ξ )) dξdη = 0 . This and the condition (3) yield b = b = 0 . Hence the final result follows from the condition (4). For the case k = 1, the condition (2) is identical tothe condition (3). A similar argument above completes the proof. (cid:3) Remark 2.2.
The space H k ( K ) is an enrichment of the BDFM element space from [15] . Hence we callthis new mixed element as the enriched BDFM element. The global space of the enriched BDFM element reads H k ( T h ) := { q ∈ H (div , Ω , R ) , q | K ∈ H k ( K ) for any K ∈ T h } . Note that, for any q ∈ H k ( T h ), the first component of q is continuous across the interior vertical edges of T h while the second component of q is continuous across the interior horizontal edges of T h .To get a stable pair of spaces, we propose to use the serendipity element of order k from [4, 14, 21] toapproximate the shear stress, which reads S k ( x, y ) := P k ( x, y ) + span { x k y, xy k } . Given any τ ∈ S k ( x, y ), it can be uniquely determined by the following conditions [4]:(1) the values of τ at four vertices of K ,(2) the values of τ at k − K ,(3) the moments R K τ pdxdy for any p ∈ P k − ( K ). JUN HU
The global space of the serendipity element of order k is defined as S k ( T h ) := { τ ∈ H (Ω) , τ | K ∈ S k ( x, y ) for any K ∈ T h } . Note that the space S ( T h ) is the usual H -conforming bilinear element space.The discrete space of the element is combined from the enriched BDFM element space and the serendipityelement space: Σ k ( K ) := { τ ∈ S , τ n ∈ H k ( K ) , τ ∈ S k ( x, y ) } . The degrees of freedom are inherited from the enriched BDFM element and the serendipity element:(1) the moments of degree not greater than k − K for σ n · ν ,(2) the moments of degree not greater than k − K for σ n ,(3) the values R K ( σ n ) J k − (2( x − x ,K ) /h x,K ) dxdy , and R K ( σ n ) J k − (2( y − y ,K ) /h y,K ) dxdy where( σ n ) is the first component of σ n , and ( σ n ) is the second component of σ n ,(4) the values R K ( σ n ) L k − (2( y − y ,K ) /h y,K ) dxdy , and R K ( σ n ) L k − (2( x − x ,K ) /h x,K ) dxdy ,(5) the values of σ at four vertices of K ,(6) the values of σ at k − K ,(7) the moments of degree not greater than k − K for σ .The definitions of the enriched BDFM element and the serendipity element imply that these conditions areunisolvent for the space Σ k ( K ). The degrees of freedom for the lowest order element is illustrated in Figure1. ✲✛ ✻❄ ①①①① ①① ①① ①① e e e e Figure 1.
Element diagram for the lowest order stress and displacementThe global space of order k is defined as(2.4) Σ k ( T h ) := { τ ∈ Σ , τ | K ∈ Σ k ( K ) for any K ∈ T h } . On each element K , the space for the displacement is taken as V k ( K ) := ( P k − ( K )) ⊕ span { ( x k , T , (0 , y k ) T , ( xy k − , T , (0 , x k − y ) T } . Then the global space for the displacement reads(2.5) V k ( T h ) := { v ∈ V, v | K ∈ V k ( K ) for any K ∈ T h } . ONFORMING MIXED FINITE ELEMENTS 7
Remark 2.3.
The lowest order element (k=1) of this family has 10 stress and 4 displacement degrees offreedom per element, which is the two dimensional element of [28] , see degrees of freedom in Figure 1.
It follows from the definitions of the spaces Σ k ( T h ) and V k ( T h ) that div Σ k ( T h ) ⊂ V k ( T h ); in the followingsection, we shall prove the converse V k ( T h ) ⊂ div Σ k ( T h ). This indicates the well-posedness of this familyof elements.The mixed element methods can be stated as: Find ( σ k,h , u k,h ) ∈ Σ k ( T h ) × V k ( T h ) such that( Aσ k,h , τ ) L (Ω) + (div τ, u k,h ) L (Ω) = 0 , (div σ k,h , v ) L (Ω) = ( f, v ) L (Ω) , (2.6)for any ( τ, v ) ∈ Σ k ( T h ) × V k ( T h ).3. Well-posedness of discrete problem for pure displacement boundary problem in twodimensions
In this section, we analyze the well-posedness of the discrete problem (2.6). From the mixed theory of[16], we need the following two assumptions(1) K-ellipticity. There exists a constant
C >
Aτ, τ ) L (Ω) ≥ C k τ k H (div , Ω) for any τ ∈ Z k ( T h ) := { τ ∈ Σ k ( T h ) , (div τ, v ) L (Ω) = 0 for all v ∈ V k ( T h ) } . (2) Discrete B-B condition. There exists a positive constant C independent of the meshsize withsup = τ ∈ Σ k ( T h ) (div τ, v ) L (Ω) k τ k H (div , Ω) ≥ C k v k L (Ω) for any v ∈ V k ( T h ) . Herein and throughout, C denotes a generic positive constant, which may be different at the differentoccurrence but independent of the meshsize h . It follows from div Σ k ( K ) ⊂ V k ( K ) for any K ∈ T h thatdiv τ = 0 for any τ ∈ Z k ( T h ). This implies the K-ellipticity condition.To prove the discrete B-B condition, the usual idea in the literature is to use the Fortin Lemma [16].More precisely, a bounded interpolation operator Π K : H ( K, S ) → Σ k ( K ) is constructed such that thefollowing commuting diagram property holds(3.1) div Π K σ = P K div σ for any σ ∈ H ( K, S ) , where P K is the projection operator from L ( K, R ) onto V k ( K ). So far, most of stable mixed finite elementmethods for the linear elasticity problem within the Hellinger-Reissner principle are designed with such aproperty, see, for instance, [1, 3, 5, 9]. However, such a technique can not be used directly herein since thereare not enough local degrees of freedom for this family of elements under consideration. The idea is to makea construction proof. More precisely, given v ∈ V k ( T h ), we find explicitly τ ∈ Σ k ( T h ) such that(3.2) div τ = v and k τ k H (div , Ω) ≤ C k v k L (Ω) . Such an idea is motivated by the stability analysis of the Raviart–Thomas element for the Poisson equation inone dimension, which is first explored to analyze the stability of a family of first order nonconforming mixedfinite element methods on the product mesh for the linear elasticity problem with the stress-displacementformulation in any dimension in a recent paper [29]. Therein, the discrete displacement is a piecewiseconstant vector, which implies that the τ of (3.2) can be directly given so that div τ = v for any v . In thispaper we use the form from [28] to construct τ . JUN HU
For convenience, suppose that the domain Ω is a unit square [0 , which is triangulated evenly into N elements, { K ij } . This implies that h x,K = h y,K = h := 1 /N for any K ∈ T h . For any v ∈ V k ( T h ), it can bedecomposed as a sum, v := ( v , v ) T = N X i =1 N X j =1 v ij ϕ ij ( x ) , where ϕ ij ( x ) is the characteristic function on the element K ij and v ij = ( v ij , v ij ) T = v | K ij . Before theconstruction of τ ∈ Σ k ( T h ) with properties of (3.2), we need a decomposition of v . We define the space V y,ij := span { , y, · · · , y k − } , which introduces the following decomposition: P k − ( K ij ) ⊕ span { x k , xy k − } = V y,ij ⊕ ( x − ih ) (cid:18) P k − ( K ij ) ⊕ span { x k − , y k − } (cid:19) . This implies that there exist unique v x,ij ∈ P k − ( K ij ) ⊕ span { x k − , y k − } and v y,ij ∈ V y,ij such that(3.3) v ij = ( x − ih ) v x,ij + v y,ij . Theorem 3.1.
It holds that sup = τ ∈ Σ k ( T h ) (div τ, v ) L (Ω) k τ k H (div , Ω) ≥ r k v k L (Ω) for any v ∈ V k ( T h ) . Proof.
Given v = ( v , v ) T ∈ V k ( T h ), we define T as the integration of v along the rectangles along the x direction:(3.4) τ ( x, y ) = Z x v ( t, y ) dt. On the element K ij := [( i − h, ih ] × [( j − h, jh ], 1 ≤ i, j ≤ N , by (3.3) and (3.4), it is straightforwardto see that τ ∈ P k ( K ) \ span { y k } ⊕ span { x k +1 , x y k − } , for ( x, y ) ∈ K ij . Similarly we can define τ as(3.5) τ ( x, y ) = Z y v ( x, t ) dt. Then by (3.4), τ is continuous in the x direction, by (3.5), τ is continuous in the y direction. Hence weget an H (div) field τ = τ τ ! ∈ Σ k ( T h ) . By the definition of τ , it follows that(3.6) div τ = v. ONFORMING MIXED FINITE ELEMENTS 9
It remains to bound the L norm of τ . We first consider the L norm of the first component τ : k τ k L (Ω) = N X i =1 N X j =1 Z K ij (cid:18) Z x v dt (cid:19) dxdy ≤ N X i =1 N X j =1 Z jh ( j − h Z ih ( i − h (cid:18) x Z x v dt (cid:19) dxdy ≤ N X i =1 N X j =1 Z jh ( j − h (cid:18) Z ih v dt (cid:19)(cid:18) Z ih ( i − h xdx (cid:19) dy = N X i =1 N X j =1 h (2 i − Z jh ( j − h Z ih v dtdy = h N X j =1 N X i =1 N X ℓ = i ( ℓ −
12 ) k v k L ( K ij ) = h N ( N − N X j =1 N X i =1 k v k L ( K ij ) ≤ k v k L (Ω) . A similar argument proves that k τ k L (Ω) ≤ k v k L (Ω) . Hence k τ k H (div , Ω) = k div τ k L (Ω) + k τ k L (Ω) ≤ k v k L (Ω) . This completes the proof. (cid:3) The pure traction boundary problem
This section considers the pure traction boundary problem, i.e., the stress space is subject to zeroNeumann boundary condition while no boundary condition on the displacement. In practice, part of theelasticity body should be located, i.e, the displacement has a Dirichlet boundary condition on some non-zero measure boundary. But the pure traction boundary problem is the most difficult one in mathematicalanalysis. A similar proof for Theorems 4.5 and 4.7 can prove them for partial displacement boundaryproblems.Let RM be the rigid motion space in two dimensions, which readsRM := span (cid:26) ! , ! , y − x ! (cid:27) . Consider a pure traction boundary problem:div σ = f in Ω := (0 , ,σν = 0 on ∂ Ω , ( u, v ) = 0 for any v ∈ RM , (4.1)where σ := A − ǫ ( u ) for u ∈ H (Ω , R ). By the same discretization of the uniform square grid T h with h = 1 /N as in the previous section, the finite element equations remain the same except the spaces are changed with boundary and rigid-motion free conditions:( Aσ h , τ ) L (Ω) + (div τ, u h ) L (Ω) = 0 for all τ ∈ Σ k, ( T h ) , (div σ h , v ) L (Ω) = ( f, v ) L (Ω) for all v ∈ V k, ( T h ) , (4.2)where Σ k, ( T h ) = { τ = τ τ τ τ ! ∈ Σ k ( T h ) , τ ν = 0 on ∂ Ω } ,V k, ( T h ) = { v = v v ! ∈ V k ( T h ) , ( v, w ) L (Ω) = 0 for all w ∈ RM } . (4.3)The earlier analysis remains the same except the discrete B-B condition as the stress space Σ k, ( T h ) issmaller than Σ k ( T h ). To prove the discrete B-B condition for the pair (Σ k, ( T h ) , V k, ( T h )), we introducethe concept of a macro-element, i.e., a union of four rectangles, see Figure 2. K K K K e e e e e e e e e e e e Figure 2.
MacroelementGiven a macroelement M , we define finite element spacesΣ k, ( M ) := { τ ∈ H (div , Ω M , S ) , τ | K ∈ Σ k ( K ) for any K ⊂ M, τ ν = 0 on ∂ Ω M } and V k ( M ) := { v ∈ L (Ω M , R ) , v | K ∈ V k ( K ) for any K ⊂ M } , where Ω M := S K ⊂ M int(K). Define the orthogonal complement space of the rigid motion space RM withrespect to V k ( M ) by RM ⊥ ( M ) := { v ∈ V k ( M ) , ( v, w ) L (Ω M ) = 0 for any w ∈ RM } . Discrete inf–sup conditions for higher order elements with k ≥ . In this subsection, we shallprove that, for k ≥
2, div Σ k, ( M ) = RM ⊥ ( M ) , which helps to establish the discrete inf–sup conditions. To this end, we define the discrete kernel space ofthe divergence operator on the macroelement M by(4.4) N M := { v ∈ V k ( M ) , (div τ, v ) L (Ω M ) = 0 for any τ ∈ Σ k, ( M ) } . We shall show that N M = RM for k ≥ J i ( ξ ) = Z ξ − J i ( s ) ds, i ≥ , ξ ∈ [ − , . ONFORMING MIXED FINITE ELEMENTS 11
For the four elements K i , i = 1 , · · · ,
4, in the macroelement M (see Figure 2), we recall the following affinemapping: ξ i = 2 x − x ,K i h x,K i , η i = 2 y − y ,K i h y,K i , ( x, y ) ∈ K i , where ( x ,K i , y ,K i ) is the center of K i , h x,K i and h y,K i are the horizontal and vertical edge lengthes of K i ,respectively. We also need the following spaces:Σ ,e : = (cid:26) τ ∈ H (Ω M ) , τ = q (1 − η ) × ( (1 + ξ ) on K (1 − ξ ) on K , q ∈ P k − ( η ) , τ = 0 on K , K (cid:27) , Σ ,e : = (cid:26) τ ∈ H (Ω M ) , τ = q (1 − ξ ) × ( (1 + η ) on K (1 − η ) on K , q ∈ P k − ( ξ ) , τ = 0 on K , K (cid:27) , Σ ,e : = (cid:26) τ ∈ H (Ω M ) , τ = q (1 − η ) × ( (1 − ξ ) on K (1 + ξ ) on K , q ∈ P k − ( η ) , τ = 0 on K , K (cid:27) , Σ ,e : = (cid:26) τ ∈ H (Ω M ) , τ = q (1 − ξ ) × ( (1 − η ) on K (1 + η ) on K , q ∈ P k − ( ξ ) , τ = 0 on K , K (cid:27) . The restriction space on M of the space S k ( T h ) of the serendipity element reads S k, ( M ) := { τ ∈ H ( M ) , τ | K i ∈ S k ( x, y ) , i = 1 , · · · , } . Note that Σ ,e i ⊂ S k, ( M ), i = 1 , , , Lemma 4.1.
For k ≥ , suppose that ( v , v ) T ∈ V k ( M ) is of the form (4.6) v | K i = a − ,i + k − X ℓ =0 a ℓ J ℓ ( η i ) and v | K i = b − ,i + k − X ℓ =0 b ℓ J ℓ ( ξ i ) , with a − , = a − , , a − , = a − , , b − , = b − , , b − , = b − , , and that Z Ω M ∂τ ∂y v + ∂τ ∂x v dxdy = 0 for any τ ∈ S k, ( M ) , then a − , = a − , = a − , = a − , , b − , = b − , = b − , = b − , , a h y,k i = − b h x,k i , a ℓ = b ℓ = 0 , ℓ = 1 , · · · , k − . Proof.
An integration by parts yields0 = Z Ω M ∂τ ∂y v + ∂τ ∂x v dxdy = − X i =1 Z K i τ (cid:0) ∂v ∂y + ∂v ∂x (cid:1) dxdy + Z e ∪ e τ ( a − , − a − , ) dx + Z e ∪ e τ ( b − , − b − , ) dy. (4.7)We take τ in (4.7) such that τ | K i ∈ (1 − ξ i )(1 − η i ) span { J ( η i ) , · · · , J k − ( η i ) , J ( ξ i ) , · · · , J k − ( ξ i ) } . This leads to 2 a h y,k i = − b h x,k i , a ℓ = b ℓ = 0 , ℓ = 1 , · · · , k − . To show these four parameters a k − , a k − , b k − and b k − to be zero, we turn to the case where k = 4. Since J ( ξ i ) = J ( η i ) = 1, J ( ξ i ) = 2 ξ i and J ( η i ) = 2 η i , we take τ ∈ Σ ,e with q = J ( η ) in (4.7). Since R e (1 − η )(1 + ξ ) J ( η ) dy = 0, this yields a = 0. Similarly, the choice of τ ∈ Σ ,e with q = J ( ξ ) shows b = 0. Then the choice of τ ∈ Σ ,e with q = J ( η ) in (4.7), yields a − , = a − , ; while the choiceof τ ∈ Σ ,e with q = J ( ξ ) in (4.7), leads to b − , = b − , . Hence we choose τ ∈ Σ ,e with q = J ( η ) ,τ ∈ Σ ,e with q = J ( ξ ) , in (4.7), respectively, to show a = b = 0. Next we consider the case where k > τ ∈ Σ ,e with q = 1 in (4.7). This leads to a − , = a − , . A similar argument with τ ∈ Σ ,e and q = 1 gets b − , = b − , . Therefore, the choices of τ ∈ Σ ,e with q = J k − ( η ) and q = J k − ( η ) ,τ ∈ Σ ,e with q = J k − ( ξ ) and q = J k − ( ξ ) , in (4.7), respectively, to prove a k − = b k − = a k − = b k − = 0 . This completes the proof. (cid:3)
Lemma 4.2.
For k = 2 , , suppose that ( v , v ) T ∈ V k ( M ) is of the form (4.8) v | K i = a − ,i + k − X ℓ =0 a ℓ J ℓ ( η i ) and v | K i = b − ,i + k − X ℓ =0 b ℓ J ℓ ( ξ i ) , with a − , = a − , , a − , = a − , , b − , = b − , , b − , = b − , , and that Z Ω M ∂τ ∂y v + ∂τ ∂x v dxdy = 0 for any τ ∈ S k, ( M ) , then a − , = a − , = a − , = a − , , b − , = b − , = b − , = b − , , a h y,k i = − b h x,k i , a ℓ = b ℓ = 0 , ℓ = 1 , · · · , k − . Proof.
We only present the details for the case where k = 3 since the proof for the case k = 2 is similar andsimple. An integration by parts yields0 = Z Ω M ∂τ ∂y v + ∂τ ∂x v dxdy = − X i =1 Z K i τ (cid:0) ∂v ∂y + ∂v ∂x (cid:1) dxdy + Z e ∪ e τ ( a − , − a − , ) dx + Z e ∪ e τ ( b − , − b − , ) dy. (4.9)For such a case, we have ∂v ∂y + ∂v ∂x | K i = 2 a h y,K i + 2 b h x,K i + 2 a J ( η i ) h y,K i + 2 b J ( ξ i ) h x,K i . Let τ ∈ Σ ,e with q = J ( η ) and τ ∈ Σ ,e with q = J ( ξ ) in (4.9), respectively. This yields a = 0and b = 0, respectively. The choices of τ ∈ Σ ,e with q = J ( η ) and τ ∈ Σ ,e with q = J ( ξ ) yield,respectively, 2 a h y,K i + 2 b h x,K i + a − , − a − , = 0 and 2 a h y,K i + 2 b h x,K i + b − , − b − , = 0 . Now we let τ | K = (1 + ξ )(1 + η ) (with appropriate definitions in K , K , and K ) in (4.9) to obtain2 a h y,K i + 2 b h x,K i + a − , − a − , + b − , − b − , = 0Finally we solve these three equations to show the desired result. (cid:3) ONFORMING MIXED FINITE ELEMENTS 13
Lemma 4.3.
It holds, for k ≥ , that (4.10) div Σ k, ( M ) = RM ⊥ ( M ) . Proof.
Since it is straightforward to see that div Σ k, ( M ) ⊂ RM ⊥ ( M ), we only need to prove that N M = span (cid:26) ! , ! , y − x ! (cid:27) . Any v = ( v , v ) ∈ V k ( M ) can be expressed as, for ℓ = 1 , · · · , v | K ℓ = X i + j ≤ k − a ( ℓ ) i,j x i y j + a ( ℓ ) k, x k + a ( ℓ )1 ,k − xy k − and v | K ℓ = X i + j ≤ k − b ( ℓ ) i,j x i y j + b ( ℓ )0 ,k y k + b ( ℓ ) k − , yx k − . We choose τ such that τ = τ = 0 and τ | K ℓ ◦ F − K ℓ ∈ (1 − ξ ) × ˆΣ ,K ℓ := span { J k − ( ξ ) , L k − ( η ) } ⊕ P k − ( ˆ K ) , ℓ = 1 , · · · , . The condition for N M implies that0 = (cid:18) ∂τ ∂x , v (cid:19) L (Ω M ) = − X ℓ =1 (cid:18) τ | K ℓ , ∂ ( v | K ℓ ) ∂x (cid:19) L (Ω M ) . Since ∂ ( v | Kℓ ) ∂x ◦ F − K ℓ ∈ ˆΣ ,K ℓ , this yields a ( ℓ ) k, = a ( ℓ )1 ,k − = a ( ℓ ) i,j = 0 for all i ≥ i + j ≤ k − . Hence v is of the form(4.11) v | K ℓ = X j ≤ k − a ( ℓ )0 ,j y j The continuity and degrees of τ across the edges e and e show (using the moments of degree not greaterthan k − e and e for τ )(4.12) a (1)0 ,j = a (2)0 ,j and a (3)0 ,j = a (4)0 ,j , j = 0 , · · · , k − . A similar argument for v shows that v is of the form(4.13) v | K ℓ = X i ≤ k − b ( ℓ ) i, x i and(4.14) b (1) i, = b (4) i, and b (2) i, = b (3) i, , i = 0 , · · · , k − . To decide these parameters a ( ℓ )0 ,j and b ( ℓ ) i, , we propose to use the degrees of freedom for the shear stresscomponent τ . We choose τ such that τ = τ = 0 and τ = τ ( ℓ )12 ∈ Σ ,e ℓ , ℓ = 1 , · · · ,
4. The conditionfor N M , and (4.11)–(4.14), produce0 = Z e ℓ τ ( ℓ )12 [ v ( ℓ +1 mod 2) ] ds − Z K ℓ τ ( ℓ )12 (cid:18) ∂v | K ℓ ∂y + ∂v | K ℓ ∂x (cid:19) dxdy − Z K ( ℓ +1 mod 4) τ ( ℓ )12 (cid:18) ∂v | K ( ℓ +1 mod 4) ∂y + ∂v | K ( ℓ +1 mod 4) ∂x (cid:19) dxdy, (4.15) where [ · ] denotes the jump of piecewise functions across edge e ℓ . Since v (resp. v ) is a piecewise polynomialwith respect to variable y (resp. x ), the symmetries of τ ( ℓ )12 with the edges e ℓ , ℓ = 1 , · · · ,
4, lead to(4.16) Z K τ (1)12 ∂v | K ∂y dxdy = Z K τ (1)12 ∂v | K ∂y dxdy, Z K τ (3)12 ∂v | K ∂y dxdy = Z K τ (3)12 ∂v | K ∂y dxdy and(4.17) Z K τ (2)12 ∂v | K ∂x dxdy = Z K τ (2)12 ∂v | K ∂x dxdy, Z K τ (4)12 ∂v | K ∂x dxdy = Z K τ (4)12 ∂v | K ∂x dxdy. Next let all τ ( ℓ )12 ∈ Σ ,e ℓ , ℓ = 1 , · · · ,
4, be defined by, up to the variable x or y , and some transformation(s),the same polynomials q of one variable of degree ≤ k −
2. Since [ v ] | e = [ v ] | e and [ v ] | e = [ v ] | e , thesymmetries of τ ( ℓ )12 imply additionally that(4.18) Z e τ (1)12 [ v ] ds = Z e τ (3)12 [ v ] ds and Z e τ (2)12 [ v ] ds = Z e τ (4)12 [ v ] ds. A substitution of equations (4.16) through (4.18) into (4.15) shows that(4.19) Z K τ (1)12 ∂v | K ∂y dxdy = Z K τ (1)12 ∂v | K ∂y dxdy = Z K τ (3)12 ∂v | K ∂y dxdy = Z K τ (3)12 ∂v | K ∂y dxdy and(4.20) Z K τ (2)12 ∂v | K ∂x dxdy = Z K τ (2)12 ∂v | K ∂x dxdy = Z K τ (4)12 ∂v | K ∂x dxdy = Z K τ (4)12 ∂v | K ∂x dxdy. Since both ∂v | K ∪ K ∂y and ∂v | K ∪ K ∂y are polynomials of degree ≤ k − y , the conditions of(4.19) for all τ (1)12 ∈ Σ ,e and τ (3)12 ∈ Σ ,e , the conditions of (4.20) for all τ (2)12 ∈ Σ ,e and τ (4)12 ∈ Σ ,e show that ( v , v ) T ∈ V k ( M ) is of the form(4.21) v | K i = a − ,i + k − X ℓ =0 a ℓ J ℓ ( η i ) and v | K i = b − ,i + k − X ℓ =0 b ℓ J ℓ ( ξ i ) , i = 1 , · · · , , with a − , = a − , , a − , = a − , , b − , = b − , , b − , = b − , . Hence it follows from Lemmas 4.1 and 4.2that a − , = a − , = a − , = a − , , b − , = b − , = b − , = b − , , a h y,k i = − b h x,k i , a ℓ = b ℓ = 0 , ℓ = 1 , · · · , k − . This completes the proof. (cid:3)
Lemma 4.4.
For any v h ∈ V k, ( T h ) , there exists a τ h ∈ Σ k, ( T h ) such that (4.22) Z M (div τ h − v h ) · wdx = 0 for any w ∈ RM and any macro–element M and (4.23) k τ h k H (div , Ω) ≤ C k v h k L (Ω) . Proof.
It is standard that there exists a τ := τ τ τ τ ! ∈ H (Ω , S ) such that(4.24) div τ = v h and k τ k H (Ω) ≤ C k v h k L (Ω) . It follows from the degrees of freedom for the enriched BDFM element in Lemma 2.1 and for the serendipityelement that there exist ( τ ,h , τ ,h ) T ∈ H k ( T h ) and τ ,h ∈ S k ( T h ) such that, for edges of macro-element M (see Figure 2 for notation), Z e ∪ e ( τ − τ ,h ) pdy = Z e ∪ e ( τ − τ ,h ) qdy = 0 for any p ∈ P ( e ∪ e ) , q ∈ P ( e ∪ e ) , ONFORMING MIXED FINITE ELEMENTS 15 Z e ∪ e ( τ − τ ,h ) pdx = Z e ∪ e ( τ − τ ,h ) qdx = 0 for any p ∈ P ( e ∪ e ) , q ∈ P ( e ∪ e ) , Z e ∪ e ( τ − τ ,h ) dy = Z e ∪ e ( τ − τ ,h ) dy = Z e ∪ e ( τ − τ ,h ) dx = Z e ∪ e ( τ − τ ,h ) dx = 0 . Let τ h = τ ,h τ ,h τ ,h τ ,h ! . We additionally have k τ h k H (div , Ω) ≤ C k τ k H (Ω) . This completes the proof. (cid:3)
We are now ready to establish the following inf–sup condition.
Theorem 4.5.
For k ≥ , there exists a positive constant C independent of the meshsize with sup = τ ∈ Σ k, ( T h ) (div τ, v ) L (Ω) k τ k H (div , Ω) ≥ C k v k L (Ω) for any v ∈ V k, ( T h ) . Proof.
Given v ∈ V k, ( T h ), it follows from Lemma 4.4 that there exists a τ ∈ Σ k, ( T h ) such that(4.25) Z M (div τ − v ) · wdx = 0 for any w ∈ RM and any macro–element M and(4.26) k τ k H (div , Ω) ≤ C k v k L (Ω) . By Lemma 4.3, there exists a τ ∈ Σ k, ( T h ) such that(4.27) div τ = div τ − v and k τ k H (div , Ω) ≤ C k div τ − v k L (Ω) . Then we have div( τ + τ ) = v and k τ k H (div , Ω) ≤ C k v k L (Ω) . (cid:3) Discrete inf–sup condition for the first order element with k = 1 . Since the analysis in theprevious subsection can not be applied to the current case, it needs a separate analysis. The ingredient isa modified macroelement technique. We also note that the macroelement technique from [34] can not beused directly here since the semi-norm | · | ,h,M there is not equivalent to the semi–norm | · | M there for thepresent case, see [34, Theorem 4.1]. To overcome this difficulty, for v ∈ V ( M ), we propose the followingmesh dependent semi–norm, see Figure 2 for notation, | v | ,h,M = X i =1 k ǫ ( v ) k ,K i + h − e k [ v ] k L ( e ) + h − e k [ v ] k L ( e ) + h − e k [ v ] k L ( e ) + h − e k [ v ] k L ( e ) + (( v | K − v | K )( M ( e )) + ( v | K − v | K )( M ( e ))+ ( v | K − v | K )( M ( e )) + ( v | K − v | K )( M ( e ))) , (4.28)where M ( e i ), i = 1 , · · · ,
4, denote the midpoints of edges e i , and [ · ] denote the jump of piecewise functionsover edge. Define a global seminorm(4.29) | v | ,h = X M | v | ,h,M for all macro-elements consisting of four elements like that in Figure 2 . It is straightforward to see that | · | ,h defines a norm over V , ( T h ). For τ ∈ Σ , ( T h ), we define the followingmesh dependent norm:(4.30) k τ k ,h = k τ k L (Ω) + X e ∈E IV h e k τ k L ( e ) + X e ∈E IH h e k τ k L ( e ) + X A ∈V I X e ∈E ( A ) h e τ ( A ) . Lemma 4.6.
For any macroelement M illustrated in Figure 2, it holds that N M = span (cid:26) ! , ! , ǫ y, − x (cid:27) , where ǫ y, − x := − ! on K , − − ! on K , − ! on K , ! on K . Proof.
Any v = ( v , v ) ∈ V ( M ) can be expressed as, for ℓ = 1 , · · · , v | K ℓ = a ( ℓ )0 + a ( ℓ )1 x and v | K ℓ = b ( ℓ )0 + b ( ℓ )1 y. We choose τ such that τ = τ = 0 and τ | K ℓ ◦ F − K ℓ = 1 − ξ , ℓ = 1 , · · · , . The condition for N M implies that a ( ℓ )1 = 0 , ℓ = 1 , · · · , . Similarly, b ( ℓ )1 = 0 , ℓ = 1 , · · · , . Hence we can use the degrees on the edges e and e for τ (using the moments of degree zero on e and e for τ ) to shows that a (1)0 = a (2)0 and a (3)0 = a (4)0 . A similar argument proves b (1)0 = b (4)0 and b (2)0 = b (3)0 . At the end we use the degree of τ at the interior vertex of M to complete the proof. (cid:3) We need another seminorm for the space V ( M ):(4.31) | v | M = sup = τ ∈ Σ , ( M ) (div τ, v ) L (Ω M ) k τ k ,h,M It follows from Lemma 4.6 that the seminorm | · | ,h,M is equivalent to the seminorm | · | M . This allows forfollowing a similar argument of [34] and the references therein to prove the discrete inf–sup condition. Theorem 4.7.
There exists a positive constant C independent of the meshsize with sup = τ ∈ Σ k, ( T h ) (div τ, v ) L (Ω) k τ k ,h ≥ C | v | ,h for any v ∈ V k, ( T h ) . ONFORMING MIXED FINITE ELEMENTS 17 Mixed finite element for three dimensions
We define a family of conforming mixed finite element methods in three dimensions in this section. Tothis end, let T h be a cuboid triangulation of the cuboid domain Ω ⊂ R such that S K ∈T h K = ¯Ω. Onelement K ∈ T h , for k ≥
1, we define an enriched Raviart–Thomas element space by H k ( K ) = ( P k,k − ,k − ( K ) ⊕ E k,x ) × ( P k − ,k,k − ( K ) ⊕ E k,y ) × ( P k − ,k − ,k ( K ) ⊕ E k,z ) , where P k,k − ,k − ( K ) = P k ( x ) × P k − ( y ) × P k − ( z ) ,P k − ,k,k − ( K ) = P k − ( x ) × P k ( y ) × P k − ( z ) ,P k − ,k − ,k ( K ) = P k − ( x ) × P k − ( y ) × P k ( z )and E k,x = x k +1 ( P k − ( y ) + P k − ( z )) ,E k,y = y k +1 ( P k − ( z ) + P k − ( x )) ,E k,z = z k +1 ( P k − ( x ) + P k − ( y )) . To construct the degrees of freedom of the space H k ( ˆ K ), we defineΨ k − ( ˆ K ) : = P k − ,k − ,k − ( ˆ K ) × P k − ,k − ,k − ( ˆ K ) × P k − ,k − ,k − ( ˆ K ) , J k − ( ξ ) : = J k − ( ξ )( P k − ( η ) + P k − ( ζ )) , J k − ( η ) : = J k − ( η )( P k − ( ξ ) + P k − ( ζ )) , J k − ( ζ ) : = J k − ( ζ )( P k − ( ξ ) + P k − ( η )) , for any ( ξ, η, ζ ) ∈ ˆ K := [ − , . We recall that J k − ( ξ ) is the Jacobi polynomial of degree k − ξ , and L k − ( ξ ) is the Legendre polynomial of degree k − ξ . Lemma 5.1.
The vector-valued function (ˆ q , ˆ q , ˆ q ) T =: ˆ q ∈ H k ( ˆ K ) can be uniquely determined by thefollowing conditions: (1) R ˆ e ˆ q · ˆ ν ˆ pd ˆ s for any ˆ p ∈ Q k − (ˆ e ) and any ˆ e ⊂ ∂ ˆ K , (2) R ˆ K ˆ q ˆ pdξdηdζ for any ˆ p ∈ J k − ( ξ ) , (3) R ˆ K ˆ q ˆ pdξdηdζ for any ˆ p ∈ J k − ( η ) , (4) R ˆ K ˆ q ˆ pdξdηdζ for any ˆ p ∈ J k − ( ζ ) , (5) R ˆ K ˆ q · ˆ pdξdηdζ for any ˆ p ∈ Ψ k − ( ˆ K ) .Proof. Since the dimensions of the space H k ( ˆ K ) is equal to the number of these conditions, it suffices toprove that ˆ q ≡ q · ˆ ν ∈ Q k − (ˆ e ), the first condition (1) implies thatˆ q = (1 − ξ )(ˆ g + ˆ f ) , ˆ q = (1 − η )(ˆ g + ˆ f ) , ˆ q = (1 − ζ )(ˆ g + ˆ f ) , where (ˆ g , ˆ g , ˆ g ) T ∈ Ψ k − ( ˆ K ), ˆ f ∈ J k − ( ξ ), ˆ f ∈ J k − ( η ), and ˆ f ∈ J k − ( ζ ). Note that Z ˆ K (1 − ξ )ˆ g ˆ pdξdηdζ = 0 for any ˆ p ∈ J k − ( ξ ) . By the condition (2), this shows ˆ f = 0. A similar argument (using the conditions (3)–(4)) yieldsˆ f = ˆ f = 0 . Finally the condition (5) proves ˆ g = ˆ g = ˆ g = 0, which completes the proof. (cid:3) To design finite element spaces for the components σ , σ , and σ , of the shear stress, we need thefollowing space S k ( X, Y ) × P k − ( Z ) for any ( X, Y, Z ) ∈ ˆ K := [ − , . Lemma 5.2.
Given any τ ∈ S k ( X, Y ) × P k − ( Z ) , it can be uniquely determined by the following conditions: (1) the values of τ at k distinct points on each edge of ˆ K that is perpendicular to the ( X, Y ) -plane, (2) the values of τ at k ( k − distinct points in the interior of each face of ˆ K that parallels to the Z -axis, (3) the moments R ˆ K τ p k − dXdY dZ for any p k − ∈ P k − ( X, Y ) × P k − ( Z ) .Here the points in the second term are chosen in this way so that they lie in the sam plane as the points inthe first term.Proof. Since S k ( X, Y ) is the space of the serendipity element of order k with respect to the variables X and Y , on each rectangle that parallels to the ( X, Y )-plane, τ can be uniquely determined bythe values of τ at four vertices of the rectangle , the values of τ at k − , the moments of order k − τ on the rectangle . Then the desired result follows from the fact that S k ( X, Y ) × P k − ( Z ) is a product space. (cid:3) Then, on element K , the space for the stress can be defined asΣ k ( K ) := { σ ∈ H (div , K, S ) | σ n ∈ H k ( K ) , σ ∈ S k ( x, y ) × P k − ( z ) ,σ ∈ S k ( x, z ) × P k − ( y ) , σ ∈ S k ( y, z ) × P k − ( x ) } . (5.1)The global space is defined asΣ k ( T h ) := { τ ∈ Σ , τ | K ∈ Σ k ( K ) for any K ∈ T h } . On each element K , the space for the displacement is taken as V k ( K ) := ( Q k − ( K )) ⊕ { ( Q k,x , , T } ⊕ { (0 , Q k,y , T } ⊕ { (0 , , Q k,z ) T } , where Q k,x = x k ( P k − ( y ) + P k − ( z )) ,Q k,y = y k ( P k − ( z ) + P k − ( x )) ,Q k,z = z k ( P k − ( x ) + P k − ( y )) . Then the global space for the displacement reads V k ( T h ) := { v ∈ V, v | K ∈ V k ( K ) for any K ∈ T h } . ONFORMING MIXED FINITE ELEMENTS 19
Remark 5.3.
The lowest order element (k=1) of this family has 21 stress and 6 displacement degrees offreedom per element, which is the three dimensional element of [28] , see some degrees of freedom in Figure3. ✚✚✚✚✚✚✚✚✚ s ❝ ❝ σ : 1 , x, y ✚✚✚✚✚✚✚✚✚ s s σ : 1 , x,y, xy s❝ ✚✚✚✚✚✚✚✚✚ ss σ : 1 , y,z, yz s❝ ✚✚✚✚✚✚✚✚✚ ❝ ❝ u : 1 , x Figure 3.
Some nodal degrees of freedom.To discretize the pure traction boundary problem, we introduce the rigid motion spaceRM := span (cid:26) , , , − yx , − z x , − zy (cid:27) , which defines Σ k, ( T h ) = { τ ∈ Σ k ( T h ) | τ ν = 0 on ∂ Ω } ,V k, ( T h ) = { v ∈ V k ( T h ) | ( v, w ) L (Ω) = 0 for all w ∈ RM } . (5.2)It follows from the definitions of the spaces Σ k ( T h ) (resp. Σ k, ( T h )) and V k ( T h ) (resp. V k, ( T h )) thatdiv Σ k ( T h ) ⊂ V k ( T h ) (resp. div Σ k, ( T h ) ⊂ V k, ( T h )). Similar arguments of Theorems 3.1, 4.5 and 4.7 canprove the converses V k ( T h ) ⊂ div Σ k ( T h ) and V k, ( T h ) ⊂ div Σ k, ( T h ), respectively. In fact, to extend theresult of Theorem 3.1 to the present case, we only need essentially three one-dimension-arguments usedin Theorem 3.1; while to get a generalization of Theorems 4.5 and 4.7, we only need essentially threetwo-dimension-arguments used in Theorems 4.5 and 4.7. In particular, in this way, we can get three two-dimension-rigid motion spaces, which proves, on a macroelement consisting of eight elements, the kernelspace N M (see (4.4) for the definition in two dimensions) is the rigid motion space in three dimension. Thisin turn implies a similar result of Lemma 4.3. Finally, this indicates the well-posedness of this family ofelements. 6. Reduced elements in both two and three dimensions
In this section we present a family of reduced elements for these in Sections 2 and 5. To this end, weintroduce Airy’s stress function for a scalar field q ( X, Y ) as follows J X,Y ( q ( X, Y )) := ∂ q∂Y − ∂ q∂X∂Y − ∂ q∂X∂Y ∂ q∂X ! . Throughout this section we let (
X, Y, Z ) denote permutations of ( x, y, z ).6.1.
The reduced elements in two dimensions.
We define the shape function space for the BDFMelement [15] as
BDF M k ( K ) := ( P k ( K )) \ span { (0 , x k ) , ( y k , } . The stress space of the reduced element of order k is defined asΣ Rk ( K ) := { τ ∈ S , τ n ∈ BDF M k ( K ) , τ ∈ P k ( x, y ) } ⊕ E k ( K ) , where E k ( K ) := span { J x,y ( x k +1 y ) , J x,y ( x y k +1 ) } . The degrees of freedom for the stress are inherited from the BDFM element and the serendipity element:(1) the moments of degree not greater than k − K for σ n · ν ,(2) the moments of degree not greater than k − K for σ n ,(3) the values of σ at four vertices of K ,(4) the values of σ at k − K ,(5) the moments of degree not greater than k − K for σ .The global space for the stress of order k is defined as(6.1) Σ Rk ( T h ) := { τ ∈ Σ , τ | K ∈ Σ Rk ( K ) for any K ∈ T h } . On each element K , the space for the displacement is taken as V Rk ( K ) := ( P k − ( K )) . Then the global space for the displacement reads(6.2) V Rk ( T h ) := { v ∈ V, v | K ∈ V Rk ( K ) for any K ∈ T h } . Remark 6.1.
Let RT k ( K ) = RT k ( K ) = P k,k − ( K ) × P k − ,k ( K ) . One can also define the space for thestress as ˆΣ k ( K ) := { τ ∈ S , τ n ∈ RT k ( K ) , τ ∈ P k ( x, y ) } ⊕ E k ( K ) . The space for the displacement in this case is ˆ V k ( K ) := ( Q k − ( K )) . The reduced elements in three dimensions.
On element K ∈ T h , for k ≥
1, we define theRaviart–Thomas element space by RT k ( K ) = P k,k − ,k − ( K ) × P k − ,k,k − ( K ) × P k − ,k − ,k ( K ) , where P k,k − ,k − ( K ) = P k ( x ) × P k − ( y ) × P k − ( z ) ,P k − ,k,k − ( K ) = P k − ( x ) × P k ( y ) × P k − ( z ) ,P k − ,k − ,k ( K ) = P k − ( x ) × P k − ( y ) × P k ( z ) . Given a scalar field q ( X, Y ) and the corresponding Airy’s function J X,Y ( q ( X, Y )), we define τ ( J X,Y q ( X, Y )) ∈ H (div , K, S ) such that τ X,X = ( J X,Y q ( X, Y )) X,X , τ
Y,X = τ X,Y = ( J X,Y q ( X, Y )) X,Y , τ
Y,Y = ( J X,Y q ( X, Y )) Y,Y
ONFORMING MIXED FINITE ELEMENTS 21 and the rest entries are zero. This notation allows to define E k ( K ) : = span (cid:26) τ ( J x,y ( x k +1 y )) , τ ( J x,y ( x y k +1 )) (cid:27) P k − ( z ) ⊕ span (cid:26) τ ( J x,z ( x k +1 z )) , τ ( J x,z ( x z k +1 )) (cid:27) P k − ( y ) ⊕ span (cid:26) τ ( J y,z ( y k +1 z )) , τ ( J y,z ( y z k +1 )) (cid:27) P k − ( x ) . Then, on element K , the space for the stress can be defined asΣ Rk ( K ) := { σ ∈ H (div , K, S ) | σ n ∈ RT k ( K ) , σ ∈ P k ( x, y ) × P k − ( z ) ,σ ∈ P k ( x, z ) × P k − ( y ) , σ ∈ P k ( y, z ) × P k − ( x ) } ⊕ E k ( K ) . (6.3)The stress τ ∈ Σ Rk ( K ) can be uniquely determined by the following conditions:(1) R e τ n · νpds for any p ∈ Q k − ( e ) and any e ⊂ ∂K ,(2) R K τ n · pdxdydz for any p ∈ Ψ k − ( K ),(3) the values of τ XY at k distinct points on each edge of K that is perpendicular to the ( X, Y )-plane,(4) the values of τ XY at k ( k −
1) distinct points in the interior of each face of K that parallels to the Z -axis,(5) the moments R K τ XY p k − dXdY dZ for any p k − ∈ P k − ( X, Y ) × P k − ( Z ).The proof for unisolvence of these degrees of freedom follows directly from those the RT element and theserendipity element, which is omitten herein; c.f. similar proofs in Lemmas 5.1 and 5.2.The global space is defined asΣ Rk ( T h ) := { τ ∈ Σ , τ | K ∈ Σ Rk ( K ) for any K ∈ T h } . On each element K , the space for the displacement is taken as V Rk ( K ) := ( Q k − ( K )) . Then the global space for the displacement reads V Rk ( T h ) := { v ∈ V, v | K ∈ V Rk ( K ) for any K ∈ T h } . To discretize the pure traction boundary problem, we defineΣ
Rk, ( T h ) = { τ ∈ Σ Rk ( T h ) | τ ν = 0 on ∂ Ω } ,V Rk, ( T h ) = { v ∈ V Rk ( T h ) | ( v, w ) = 0 for all w ∈ RM } . (6.4)It follows from the definitions of the spaces Σ Rk ( T h ) (resp. Σ Rk, ( T h )) and V Rk ( T h ) (resp. V Rk, ( T h )) thatdiv Σ Rk ( T h ) ⊂ V Rk ( T h ) (resp. div Σ Rk, ( T h ) ⊂ V Rk, ( T h )). Similar arguments of Theorems 3.1, 4.5 and 4.7can prove the converses V Rk ( T h ) ⊂ div Σ Rk ( T h ) and V Rk, ( T h ) ⊂ div Σ Rk, ( T h ), respectively. This indicates thewell-posedness of this family of elements. Remark 6.2.
The lowest order element (k=1) of this family has 8 stress and 2 displacement, and 18 stressand 3 displacement degrees of freedom per element for two and three dimensions, respectively, which wereannounced independently by Chen and his collaborators [20] after the first version of the paper was submitted. The error estimate and numerical results
The error estimate.
The section is devoted to the error analysis of the approximation defined by(2.6). It follows from (1.1) and (2.6) that( A ( σ − σ k,h ) , τ h ) L (Ω) + (div τ h , ( u − u k,h )) L (Ω) = 0 for any τ h ∈ Σ k ( T h ) , (div( σ − σ k,h ) , v h ) L (Ω) = 0 for any v h ∈ V k ( T h ) . (7.1)Let P h be the L projection operator from L (Ω , R n ) onto V k ( T h ). Since div σ k,h ∈ V k ( T h ), the secondequation of (7.1) yields k div( σ − σ k,h ) k L (Ω) = k div σ − P h div σ k L (Ω) ≤ Ch m | div σ | H m (Ω) for any 0 ≤ m ≤ k. (7.2)It follows from the K-ellipticity, Theorem 3.1 and the approximation properties of Σ k,h and V k,h that k σ − σ k,h k L (Ω) + k u − u k,h k L (Ω) ≤ C (cid:0) inf τ h ∈ Σ k,h k σ − τ h k H (div , Ω) + inf v h ∈ V k,h k u − v h k L (Ω) (cid:1) ≤ Ch m ( | σ | H m +1 (Ω) + | u | H m (Ω) ) for any 0 ≤ m ≤ k. (7.3)A similar error estimate holds for the pure traction boundary problem studied in section 4 and the reducedelements in Section 6. Remark 7.1.
By using the mesh dependent norm in Subsection 4.2, cf. (4.28) and (4.30) , we can get animproved error estimate: k σ − σ k,h k L (Ω) ≤ Ch m | σ | H m (Ω) for any ≤ m ≤ k. The numerical result.
The first example is presented to demonstrate the second order method (with k = 2) for the pure displacement boundary problem with a homogeneous boundary condition that u ≡ ∂ Ω; see [28] for numerical examples for k = 1. Assume the material is isotropic in the sense that Aσ = 12 µ (cid:18) σ − λ µ + 2 λ tr( σ ) δ (cid:19) , where δ is the identity matrix, and µ and λ are the Lam´e constants such that 0 < µ ≤ µ ≤ µ and0 < λ < ∞ . In the numerical example, these parameters are chosen as λ = 1 µ = 12 . Let the exact solution on the unit square [0 , be(7.4) u = (sin πx sin πy, sin πx sin πy ) T . In the computation, the level one grid is the given domain, a unit square or a unit cube. Each grid isrefined into a half-size grid uniformly, to get a higher level grid, see the first column in Table 1.As the second example, we compute the pure traction boundary problem with the exact solution(7.5) u = (cid:20) x (1 − x ) y (1 − y ) − (cid:21) − ! . The matrix A is same as that in the first example. Our new finite element has no problem in solving thepure traction boundary problems. The convergence results are listed in Table 2. ONFORMING MIXED FINITE ELEMENTS 23
Table 1.
The error and the order of convergence. k u − u ,h k rate k σ − σ ,h k rate k div( σ − σ ,h ) k rate1 0.3156 0.0 2.0116 0.0 7.8083 0.02 0.0693 2.2 0.4465 2.2 1.9752 2.03 0.0166 2.1 0.1134 2.0 0.4760 2.14 0.0041 2.0 0.0285 2.0 0.1175 2.05 0.0010 2.0 0.0071 2.0 0.0293 2.06 2.5408e-004 2.0 0.0018 1.9 0.0073 2.07 6.3503e-005 2.0 4.4605e-004 2.0 0.0018 2.0 Table 2.
The errors and the order of convergence for the pure traction boundary problem k u − u ,h k rate k σ − σ ,h k rate k div( σ − σ ,h ) k rate2 0.0264 0.0 0.2516 0.0 2.4645 0.03 0.0107 1.3 0.0804 1.6 0.7090 1.84 0.0029 1.9 0.0211 2.0 0.1807 2.05 7.2940e-004 2.0 0.0054 2.0 0.0453 2.06 1.8315e-004 2.0 0.0013 2.0 0.0113 2.07 4.5836e-005 2.0 3.3684e-004 2.0 0.0028 2.0 References [1] S. Adams and B. Cockburn, A mixed finite element method for elasticity in three dimensions, J. Sci. Comput. 25 (2005),no. 3, 515–521.[2] M. Amara and J. M. Thomas, Equilibrium finite elements for the linear elastic problem, Numer. Math. 33 (1979),367–383.[3] D. N. Arnold and G. Awanou, Rectangular mixed finite elements for elasticity, Math. Models Methods Appl. Sci. 15(2005), 1417–1429.[4] D. N. Arnold and G. Awanou, The serendipity family of finite elements, Found. Comput. Math. 11(2011), 337–344.[5] D. Arnold, G. Awanou and R. Winther, Finite elements for symmetric tensors in three dimensions, Math. Comp. 77(2008), no. 263, 1229–1251.[6] D. N. Arnold, F. Brezzi and J. Douglas, Jr., PEERS: A new mixed finite element for plane elasticity, Jpn. J. Appl. Math.1 (1984), 347–367.[7] D. N. Arnold, J. Douglas Jr., and C. P. Gupta, A family of higher order mixed finite element methods for plane elasticity,Numer. Math. 45 (1984), 1–22.[8] D.N. Arnold, R. Falk and R. Winther, Mixed finite element methods for linear elasticity with weakly imposed symmetry,Math. Comp. 76 (2007), no. 260, 1699–1723.[9] D. N. Arnold and R. Winther, Mixed finite element for elasticity, Numer. Math. 92 (2002), 401–419.[10] D. N. Arnold and R. Winther, Nonconforming mixed elements for elasticity, Math. Models. Methods Appl. Sci. 13 (2003),295–307.[11] G. Awanou, Two remarks on rectangular mixed finite elements for elasticity, J. Sci. Comput. 50 (2012), 91–102.[12] E. B´ecache, P. Joly and C. Tsogka, A new family of mixed finite elements for the linear elastodynamic problem , SIAMJ. Numer. Anal., 39(2002), pp. 2109–2132.[13] D. Boffi, F. Brezzi and M. Fortin, Reduced symmetry elements in linear elasticity, Commun. Pure Appl. Anal. 8 (2009),no. 1, 95–121.[14] S. C. Brenner and L. R. Scott, The mathematical theorey of finite element methods, Springer-Verlag, 1996.[15] F. Brezzi, J. Douglas, Jr., M. Fortin, L. D. Marini, Efficient rectangular mixed finite elements in two and three spacevariables, A. I. R.O., Mode. Math. Anal. Numer., 21(1987): 581–604.[16] F. Brezzi and M. Fortin, Mixed and hybrid finite element methods, Springer, 1991. [17] C. Carstensen, M. Eigel and J. Gedicke, Computational competition of symmetric mixed FEM in linear elasticity,Comput. Methods Appl. Mech. Engrg. 200 (2011), no. 41-44, 2903–2915.[18] C. Carstensen, D. G¨unther, J. Reininghaus, J. Thiele, The Arnold–Winther mixed FEM in linear elasticity. Part I:Implementation and numerical verification, Comput. Methods Appl. Mech. Engrg. 197 (2008) 3014–3023.[19] S.-C. Chen and Y.-N. Wang, Conforming rectangular mixed finite elements for elasticity, J. Sci. Comput. 47 (2011), no.1, 93–108.[20] S. C. Chen, Presentation in the workshop on “Finite element methods and its applications”, Beijing, China, December7 2013.[21] P. G. Ciarlet, The finite element method for elliptic problems, North–Holland, 1978; reprinted as SIAM Classics inApplied Mathematics, 2002.[22] B. Cockburn, J. Gopalakrishnan and J. Guzm´an, A new elasticity element made for enforcing weak stress symmetry,Math. Comp. 79 (2010), no. 271, 1331–1349.[23] V. Girault and P. A. Raviart, Finite element methods for Navier-Stokes equations: theory and algorithms, Springer-Verlag, Berlin, Heidelberg 1986.[24] J. Gopalakrishnan and J. Guzm´an, Symmetric nonconforming mixed finite elements for linear elasticity, SIAM J. Numer.Anal. 49 (2011), no. 4, 1504–1520.[25] J. Gopalakrishnan and J. Guzm´an, A second elasticity element using the matrix bubble, IMA J. Numer. Anal. 32 (2012),no. 1, 352–372.[26] J. Guzm´an, A unified analysis of several mixed methods for elasticity with weak stress symmetry, J. Sci. Comput. 44(2010), no. 2, 156–169.[27] Jason S. Howell, Noel J. Walkington, Inf-sup conditions for twofold saddle point problems. Numer. Math. 118(2011), No.4, 663–693.[28] J. Hu, H. Y. Man, and S. Y. Zhang, A simple conforming mixed finite element for linear elasticity on rectangular gridsin any space dimension, J. Sci. Comput. DOI 10.1007/s10915-013-9736-6.[29] J. Hu, H. Y. Man, and S. Y. Zhang, A minimal mixed finite element method for linear elasticity in the symmetricformulation on n -rectangular grids, arXiv:1304.5428[math.NA] (2013).[30] J. Hu and Z. C. Shi, Lower order rectangular nonconforming mixed elements for plane elasticity, SIAM J. Numer. Anal.46 (2007), 88–102.[31] C. Johnson and B. Mercier, Some equilibrium finite element methods for two-dimensional elasticity problems, Nu-mer.Math. 30 (1978), 103–116.[32] H.-Y. Man, J. Hu and Z.-C. Shi, Lower order rectangular nonconforming mixed finite element for the three-dimensionalelasticity problem, Math. Models Methods Appl. Sci. 19 (2009), no. 1, 51–65.[33] M. Morley, A family of mixed finite elements for linear elasticity Numer. Math. 55 (1989), no. 6, 633-666.[34] R. Stenberg, On the construction of optimal mixed finite element methods for the linear elasticity problem, Numer.Math. 48 (1986), 447–462.[35] R. Stenberg, Two low-order mixed methods for the elasticity problem, In: J. R. Whiteman (ed.): The Mathematics ofFinite Elements and Applications, VI. London: Academic Press, 1988, 271–280.[36] R. Stenberg, A family of mixed finite elements for the elasticity problem, Numer. Math. 53 (1988), no. 5, 513–538.[37] R. Stenberg, A technique for analysing finite element methods for viscous incompressible flow, Internat. J. Numer.Methods Fluids, 11 (1990), pp. 935–948.[38] S. Y. Yi, A New nonconforming mixed finite element method for linear elasticity, Math. Models Methods Appl. Sci. 16(2006), 979–999.-rectangular grids, arXiv:1304.5428[math.NA] (2013).[30] J. Hu and Z. C. Shi, Lower order rectangular nonconforming mixed elements for plane elasticity, SIAM J. Numer. Anal.46 (2007), 88–102.[31] C. Johnson and B. Mercier, Some equilibrium finite element methods for two-dimensional elasticity problems, Nu-mer.Math. 30 (1978), 103–116.[32] H.-Y. Man, J. Hu and Z.-C. Shi, Lower order rectangular nonconforming mixed finite element for the three-dimensionalelasticity problem, Math. Models Methods Appl. Sci. 19 (2009), no. 1, 51–65.[33] M. Morley, A family of mixed finite elements for linear elasticity Numer. Math. 55 (1989), no. 6, 633-666.[34] R. Stenberg, On the construction of optimal mixed finite element methods for the linear elasticity problem, Numer.Math. 48 (1986), 447–462.[35] R. Stenberg, Two low-order mixed methods for the elasticity problem, In: J. R. Whiteman (ed.): The Mathematics ofFinite Elements and Applications, VI. London: Academic Press, 1988, 271–280.[36] R. Stenberg, A family of mixed finite elements for the elasticity problem, Numer. Math. 53 (1988), no. 5, 513–538.[37] R. Stenberg, A technique for analysing finite element methods for viscous incompressible flow, Internat. J. Numer.Methods Fluids, 11 (1990), pp. 935–948.[38] S. Y. Yi, A New nonconforming mixed finite element method for linear elasticity, Math. Models Methods Appl. Sci. 16(2006), 979–999.