1 Introduction

# Nonconforming tetrahedral mixed finite elements for elasticity

## Abstract.

This paper presents a nonconforming finite element approximation of the space of symmetric tensors with square integrable divergence, on tetrahedral meshes. Used for stress approximation together with the full space of piecewise linear vector fields for displacement, this gives a stable mixed finite element method which is shown to be linearly convergent for both the stress and displacement, and which is significantly simpler than any stable conforming mixed finite element method. The method may be viewed as the three-dimensional analogue of a previously developed element in two dimensions. As in that case, a variant of the method is proposed as well, in which the displacement approximation is reduced to piecewise rigid motions and the stress space is reduced accordingly, but the linear convergence is retained.

###### Key words and phrases:
mixed method, finite element, linear elasticity, nonconforming.
###### 2000 Mathematics Subject Classification:
Primary: 65N30, Secondary: 74S05
The work of the first author was supported by NSF grant DMS-1115291. The work of the second author was supported by NSF grant DMS-0811052 and the Sloan Foundation. The work of the third author was supported by the Research Council of Norway through a Centre of Excellence grant to the Centre of Mathematics for Applications.

## 1. Introduction

Mixed finite element methods for elasticity simultaneously approximate the displacement vector field and the stress tensor field. Conforming methods based on the classical Hellinger–Reissner variational formulation require a finite element space for the stress tensor that is contained in , the space of symmetric tensor fields which are square integrable with square integrable divergence. For a stable method, this stress space must be compatible with the finite element space used for the displacement, which is a subspace of the vector-valued function space. It has proven difficult to devise such pairs of spaces. While some stable pairs have been successfully constructed in both and dimensions, the resulting elements tend to be quite complicated, especially in dimensions. For this reason, much attention has been paid to constructing elements which fulfill desired stability, consistency, and convergence conditions, but which relax the requirement that the stress space be contained in in one of two ways: either by relaxing the interelement continuity requirements, which leads to nonconforming mixed finite elements, or by relaxing the symmetry requirement, which leads to mixed finite elements with weak symmetry. In this paper we construct a new nonconforming mixed finite element for elasticity in three dimensions based on tetrahedral meshes, analogous to a two-dimensional element defined in [11]. The space of shape functions on a tetrahedral element (which is defined in (3.1) below) is a subspace of the space , the space of symmetric tensors with components which are polynomials of degree at most . It contains and has dimension 42. The degrees of freedom for are the integral of over (this is six degrees of freedom, since has six components), and the integral and linear moments of on each face of (nine degrees of freedom per face). For the displacements we simply take as the shape functions and use only interior degrees of freedom so as not to impose any interelement degrees of freedom. See the element diagrams in Figure 1. We note that, since there are no degrees of freedom associated to vertices or edges, only to faces and the interior, our elements may be implemented through hybridization, which may simplify the implementation. See [5] for the general idea, or [18] for a case close to the present one.

After some preliminaries in section 2, in section 3 we define the shape function space and prove unisolvence of the degrees of freedom. In section 4 we establish the stability, consistency, and convergence of the resulting mixed method. Finally in section 5 we describe a variant of the method which reduces the displacement space to the space of piecewise rigid motions and reduces the stress space accordingly. The results of this paper were announced in [13].

As mentioned, conforming mixed finite elements for elasticity tend to be quite complicated. The earliest elements, which worked only in two dimensions, used composite elements for stress [22, 7]. Much more recently, elements using polynomial shape functions were developed for simplicial meshes in two [10] and three dimensions [1, 4] and for rectangular meshes [3, 14]. Heuristics given in [10] and [4] indicate that it is not possible to construct significantly simpler elements with polynomial shape functions and which preserve both the conformity and symmetry of the stress. Many authors have developed mixed elements with weak symmetry [2, 6, 26, 28, 27, 24, 8, 9, 17, 20, 15, 19], which we will not pursue here. For nonconforming methods with strong symmetry, which is the subject of this paper, there have been several elements proposed for rectangular meshes [29, 30, 21, 12, 23], but very little work on simplicial meshes. A two-dimensional nonconforming element of low degree was developed by two of the present authors in [11]. As shape functions for stress it uses a dimensional subspace of the space of all quadratic symmetric tensors, while for the displacement it uses piecewise linear vector fields. A second element was also introduced in [11], for which the stress shape function space was reduced to dimension and the displacement functions reduced to the piecewise rigid motions. In [18] Gopalakrishnan and Guzmán developed a family of simplicial elements, in both two and three dimensions. As shape functions they used the space of all symmetric tensors of polynomial degree at most , paired with piecewise polynomial vector fields of dimension , for . Thus, in two dimensions and in the lowest degree case, they use an dimensional space of shape functions for stress, while in three dimensions, the space has dimension . Gopalakrishnan and Guzmán also proposed a reduced variant of their space, in which the displacement space remains the full space of piecewise polynomials of degree , but the dimension of the stress space is reduced to in two dimensions and to in three dimensions. However, their reduced spaces have a drawback, in that they are not uniquely defined, but for each edge of the triangulation require a choice of a favored endpoint of the edge. In particular, in two dimensions, the reduced space of [18] uses the same displacement space as the non-reduced space of [11], uses a stress space of the same dimension, and uses identical degrees of freedom, but the two spaces do not coincide (since the space of [11] does not require a choice of favored edge endpoints).

The elements introduced here may be regarded as the three-dimensional analogue of the element in [11]. Again, they have the same displacement space and the same degrees of freedom as the reduced three-dimensional elements of [18], but the stress spaces do not coincide. Also, as in the two-dimensional case, our reduced space is of lower dimension than any that has been heretofore proposed.

## 2. Preliminaries

Let be a bounded domain. We denote by the space of symmetric matrices and by and the space of square-integrable vector fields and symmetric matrix fields on , respectively. The space consists of matrix fields with row-wise divergence, , in . The Hellinger–Reissner variational formulation seeks such that

 (2.1) ∫Ω(Aσ:τ+divτ⋅u)dx =0, τ∈H(div,Ω;S) ∫Ωdivσ⋅vdx =∫Ωf⋅vdx, v∈L2(Ω;Rn).

Here denotes the Frobenius inner products of matrices and , and denotes the compliance tensor, a linear operator which is bounded and symmetric positive definite uniformly for . The solution solves the Dirichlet problem for the Lamé equations and so belongs to . If the domain is smooth and the compliance tensor is smooth, then and

 (2.2) ∥σ∥1+∥u∥2≤c∥f∥0.

with a constant depending on and . The same regularity holds if the domain is a convex polyhedron, at least in the isotropic homogeneous case. See [25].

We shall also use spaces of the form where is a finite-dimensional vector space and is a nonnegative integer, the Sobolev space of functions for which all derivatives of order at most are square integrable. The norm is denoted by or .

To discretize (2.1), we choose finite-dimensional subspaces and . Assuming that consists of matrix fields which are piecewise polynomial with respect to some mesh of , we define by applying the (row-wise) divergence operator piecewise. A mixed finite element approximation of (2.1) is then obtained by seeking such that:

 (2.3) ∫Ω(Aσh:τ+divhτ⋅uh)dx =0, τ∈Σh ∫Ωdivhσh⋅vdx =∫Ωf⋅vhdx, v∈Vh.

If this is a conforming method, otherwise, as for the elements developed below, it is nonconforming. We recall that a piecewise smooth matrix field belongs to if and only if whenever two tetrahedra in meet in a common face, the jump of the normal components across the face vanish.

## 3. Definition of the new elements

We define the finite element spaces and in the usual way, by specifying spaces of shape functions and degrees of freedom. The space is simply the space of all piecewise linear vector fields with respect to the given tetrahedral mesh of (which we therefore assume is polyhedral). Thus the shape function space on an element is simply , the space of polynomial vector fields on of degree at most . For degrees of freedom we choose the moments with weights . Since no degrees of freedom are associated with the proper subsimplices of , no interelement continuity is imposed on . The associated projection is the projection.

To define the space we introduce some notation. If is a unit vector, let be the orthogonal projection onto the plane orthogonal to . Then is given by the symmetric matrix . For a tetrahedron , let denote the subsimplices of dimension (vertices, edges, faces, and tetrahedra) of . For an edge let be a unit vector parallel to and, for a face , let be its outward unit normal. We can then define the shape function space

 (3.1) ΣK ={σ∈P2(K;S) | QseσQse|e∈P1(e;S)∀e∈Δ1(K)}.

For , is a quadratic polynomial on taking values in the -dimensional subspace of . As illustration, for and , we have

 QseσQse=⎛⎜⎝σ11σ120σ12σ220000⎞⎟⎠.

Thus the requirement that belong to represents linear constraints on , and so . We shall now specify degrees of freedom (linear functionals) and show unisolvence, i.e., that if all the degrees of freedom vanish for some , then vanishes. This will imply that , and so the dimension is exactly .

The degrees of freedom we take are:

 (3.2) ∫fσnf⋅vds,v∈P1(f;R3), f∈Δ2(K), (36 degrees of freedom), (3.3) ∫Kσdx, (6 degrees of freedom).

The following lemma will be used in the proof of unisolvence.

###### Lemma 3.1.

Let and be the faces of opposite two distinct vertices and and let be their common edge, with endpoints and . Given , there exists a unique satisfying the following four conditions (see Figure 2):

1. ,

2. ,  ,

3. , ,

4. .

Moreover .

###### Proof.

For uniqueness we must show that if satisfies (1)–(4) with , then vanishes. Certainly, from (1) and (2), vanishes on , and then, using (3), vanishes on and . Therefore where is the barycentric coordinate function equal to on and at , similarly for , and is a constant. Integrating this equation over and invoking (4) we conclude that does indeed vanish.

To show the existence of , we simply exhibit its formula in terms of barycentric coordinates:

 p=βλ2k+(β+γ)λkλl+γλ2l+32(β+γ)(λ2i+λ2j)+(−5β−γ)(λi+λj)λk+(−β−5γ)(λi+λj)λl+3(β+γ)λiλj.

That this function satisfies (1)–(4) follows from the elementary formula

 ∫Tλα=α1!⋯αd+1!d!(|α|+d)!|T|,α∈Nd+10,

for the integral of a barycentric monomial over a simplex of dimension , which can be established by induction (see, e.g., [16]). ∎

We are now ready to prove the claimed unisolvence result.

###### Theorem 3.2.

The degrees of freedom given by (3.2) and (3.3) are unisolvent for the shape function space defined by (3.1): if the degrees of freedom all vanish for some , then .

###### Proof.

Let be the gradient of the th barycentric coordinate function. Thus is an inward normal vector to the face with length where is the distance from the th vertex to . Note that any three of the form a basis for and that .

For , define . We shall show that if and all the degrees of freedom vanish, then on for all . This is sufficient, since, fixing and varying , we conclude that , and, then, since this holds for each , that .

If is an edge of the faces and of , which may or may not coincide, then . Thus, from the definition (3.1) of the space , is linear on . In particular, is linear on each edge of . Thus is a quadratic polynomial on whose restriction to each edge of is linear. Therefore, on the boundary of , coincides with its linear interpolant, and, since a quadratic function on a triangle is determined by its boundary values, is linear. Thus is actually a linear polynomial on , and, in view of the degrees of freedom (3.2), we conclude that vanishes on .

For any pair of distinct indices (that is, and ), define

 (3.4) βlk=σij(vk),βkl=σij(vl),

where are the two indices unequal to and . Now is linear on the common edge of and , and, because of the vanishing degrees of freedom of , is orthogonal to on and on and has integral on . Therefore, by Lemma 3.1 applied with , it is sufficient to show that and both vanish in order to conclude that vanishes. In fact, we shall show that the 12 quantities , corresponding to the 12 pairs of distinct indices, satisfy a nonsingular homogeneous system of 12 equations, and so vanish.

The lemma also tells us that . Interchanging and gives

 σik(vk)=32(βlj+βjl).

Also, by definition,

 (3.5) βjk=σil(vk).

Combining (3.4)–(3.5) gives

 σij(vk)+σik(vk)+σil(vk)=32(βlj+βjl)+(βlk+βjk).

But , which vanishes on and so, in particular, at the vertex . Thus we have established the equation

 (3.6) a(βlj+βjl)+b(βlk+βjk)=0,

where , .

For each of the 12 pairs of distinct indices, we let and be the remaining indices and consider the equation (3.6). In this way we obtain a system of 12 linear equations in 12 unknowns. If we number the pairs of distinct indices lexographically, the matrix of the system is:

 ⎛⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜⎝0000000ba0ba0000ba0000ab0000ab0ab000000000b0ab0a0ba000000a0b0ab000a0b000000b0a000ba0b0a000000ab0a0ba0b000000000ba0ba0000ba0000ab0000ab0ab0000000⎞⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟⎠.

Its determinant is , as may be verified with a computer algebra package. In particular, when , , the system is nonsingular. Thus all the vanish as claimed, and the proof is complete. ∎

Having established unisolvency, the assembled finite element space is defined as the set of all matrix fields such that for all and for which the degrees of freedom (3.2) have a common value when a face is shared by two tetrahedra in . If , then the jump of across such an interior face need not vanish, but it is orthogonal to . The normal component is, by the definition of the shape function space, linear on each edge of so belongs to , and thus

 (3.7) \llbracketn′fτnf\rrbracket=0on f,

for any interior face of the triangulation.

## 4. Error analysis

In this section, we show that the pair of spaces , give a convergent finite element method. The argument follows the one given in [11] for the two-dimensional case. As usual, we suppose that we are given a sequence of tetrahedral meshes indexed by a parameter which decreases to zero and represents the maximum tetrahedron diameter. We assume that the sequence is shape regular (the ratio of the diameter of a tetrahedron to the diameter of its inscribed ball is bounded), and the constants which appear in the estimates below may depend on this bound, but are otherwise independent of .

We start by observing that, by construction,

 (4.1) divhΣh⊂Vh.

The degrees of freedom determine an interpolation operator by

 ∫f(Πhτ−τ)n⋅ vds =0,v∈P1(f), f∈Δ1(Th), ∫K(Πhτ−τ)dx =0,K∈Th,

where . Since

 ∫K(divΠhτ−divτ)⋅vdx=−∫K(Πhτ−τ):ϵ(v)dx+∫∂K(Πhτ−τ)n⋅vds=0,

for , , , we have the commutativity property

 (4.2) divhΠhτ=Phdivτ,τ∈H1(Ω;S).

Since maps onto , (4.2) implies that maps onto . An immediate consequence is that the finite element method system (2.3) is nonsingular. Indeed, if , then the choice of test functions and implies that and then, choosing with , we get .

For the error analysis we also need the approximation and boundedness properties of the projections and . Obviously, for the projection, we have

 (4.3) ∥v−Phv∥0≤chm∥v∥m,0≤m≤2.

Since is defined element-by-element and preserves piecewise linear matrix fields, we may scale to a reference element of unit diameter using translation, rotation, and dilation, and use a compactness argument, to obtain

 (4.4) ∥τ−Πhτ∥0≤chm∥τ∥m,m=1,2,

where the constant depends only on the shape regularity of the elements. See, e.g., [10] for details. Taking and using the triangle inequality establishes boundedness of :

 (4.5) ∥Πhτ∥0≤c∥τ∥1.

The final ingredient we need for the convergence analysis is a bound on the consistency error arising from the nonconformity of the elements. Define

 (4.6) Eh(u,τ)=∫Ω[ϵ(u):τ+divhτ⋅u]dx,u∈\ringH1(Ω;R3), τ∈Σh+H(div,Ω;S).

If , then , by integration by parts. In general,

 Eh(u,τ)=∑K∈Th∫∂KτnK⋅uds=∑f∈Δ2(Th)∫f\llbracketτnf\rrbracket⋅uds,

where, again, denotes the jump of across the face . Only the interior faces enter the sum, since vanishes on . Now , so

 Eh(u,τ) =∑f∈Δ2(Th){∫f\llbracketQnf(τnf)\rrbracket⋅uds+∫f\llbracketn′fτnf\rrbracket(n′fu)ds} =∑f∈Δ2(Th)∫f\llbracketQnf(τnf)\rrbracket⋅uds,

where the last equality follows from (3.7).

We let be the subspace of the displacement space consisting of continuous functions which are zero on the boundary of . In other words, is the standard piecewise linear subspace of . For any the jumps, , are orthogonal to , so for any .

###### Lemma 4.1.

We may bound the consistency error

 (4.7) |Eh(u,τ)|≤ch(∥τ∥0+h∥divhτ∥0)∥u∥2,τ∈Σh, u∈\ringH1(Ω;R3)∩H2(Ω;R3).

Furthermore, for any

 (4.8) |Eh(u,Πhρ)|≤ch2∥ρ∥1∥u∥2,u∈\ringH1(Ω;R3)∩H2(Ω;R3).
###### Proof.

For any we have , where is the piecewise linear interpolant of . Referring to the definition (4.6), we obtain

 |Eh(u,τ)|≤c(∥divhτ∥0∥u−uIh∥0+∥τ∥0||ϵ(u−uIh)||0≤ch(∥τ∥0+h∥divhτ∥0)∥u∥2,

which is (4.7). For the second estimate we use that , which implies that

 Eh(u,Πhρ)=∑K∈Th∫Kdivh(Πhρ−ρ)⋅(u−uIh)dx+∫K(Πhρ−ρ):ϵ(u−uIh)dx.

Utilizing the estimate (4.4), the bound

 |Eh(u,Πhρ)|≤c(∥divρ∥0∥u−uIh∥0+∥Πhρ−ρ∥0||ϵ(u−uIh)||0≤ch2∥ρ∥1∥u∥2

is an immediate consequence. ∎

###### Remark 4.2.

The consistency error estimate (4.7) holds for any satisfying for each , provided one replaces with the broken norm

With these ingredients assembled, error bounds for the finite element method now follow in a straightforward fashion.

###### Theorem 4.3.

Let be the solution of (2.1) and the solution of (2.3). Then

 ∥σ−σh∥0 ≤ch∥u∥2, (4.9) ∥divσ−divhσh∥0 ≤chm∥divσ∥m,0≤m≤2, ∥u−uh∥0 ≤ch∥u∥2.

Furthemore, if problem (2.1) admits full elliptic regularity, such that the estimate (2.2) holds, then

 ∥u−uh∥0≤ch2∥u∥2.
###### Proof.

Subtracting the first equations of (2.1) and (2.3) and invoking the definition (4.6) of the consistency error, we get the error equation

 (4.10) ∫Ω[A(σ−σh):τ+(u−uh)⋅divhτ]dx=Eh(u,τ),τ∈Σh.

Comparing the second equations in (2.1) and (2.3), we obtain , which immediately gives the claimed error estimate on . Using the commutativity (4.2), we find that . Choosing in (4.10), we get

 ∫ΩA(σ−σh):(Πhσ−σh)dx=Eh(u,Πhσ−σh),

which implies that

 ∥σ−σh∥2A≤∥σ−Πhσ∥2A+2Eh(u,Πhσ−σh),

where . Combining with (4.4) and (4.7) we conclude that

 ∥σ−σh∥≤ch(∥σ∥1+∥u∥2)≤ch∥u∥2,

which is the desired error estimate for .

To get the error estimate for , we choose such that and . Then, in light of the commutativity property (4.2) and the bound (4.5), satisfies and . Hence, using (4.1), (4.10), and (4.7), we get

 ∥Phu −uh∥20=∫Ωdivhτ⋅(Phu−uh)dx=∫Ωdivhτ⋅(u−uh)dx (4.11) =−∫ΩA(σ−σh):τdx+Eh(u,τ)≤c(∥σ−σh∥0+h∥u∥2)∥Phu−uh∥0.

This gives , and then, by the triangle inequality and (4.3), the error estimate for .

To establish the final quadratic estimate for in the case of full regularity, we use a duality argument. Let , where solves the problem . It follows from (2.2) that

 (4.12) ∥ρ∥1+∥w∥2≤c∥Phu−uh∥0.

By introducing as the piecewise linear interpolant of , we now obtain from (4) that

 ∥Phu−uh∥20 =−∫ΩA(σ−σh):Πhρdx+Eh(u,Πhρ) =−∫ΩA(σ−σh):(Πhρ−ρ)dx+Eh(u,Πhρ)−∫Ω(σ−σh):ϵ(w−wIh)dx,

where the final equality follows since

 ∫Ω(σ−σh):ϵ(wIh)dx=−∑K∈Th∫Kdivh(σ−σh)⋅wIhdx+Eh(wIh,σ−σh)=0.

However, by utilizing (4.4), (4.8), the estimate for given in (4.3), combined with the approximation property of the interpolant , we obtain from the representation of above that

 ∥Phu−uh∥20 ≤c(h2∥ρ∥1∥u∥2+∥σ−σh∥∥ϵ(w−wIh)∥0) ≤ch2∥u∥2(∥ρ∥1+∥w∥2)≤ch2∥u∥2∥Phu−uh∥0,

where we have used (4.12) to obtain the final inequality. This gives . As above, the desired estimate for now follows from (4.3) and the triangle inequality. ∎

###### Remark 4.4.

Although , we have only shown first order convergence of the finite element solution: . The lower rate of convergence is due to the consistency error estimated in (4.7).

## 5. The reduced element

As for the two-dimensional element in [11], there is a variant of the element using smaller spaces. Let

 T(K)={v∈P1(K;R3)|v(x)=a+b×x, a,b∈R3},

be the space of rigid motions on . In the reduced method we take instead of as the space of shape functions for displacement, so the dimension is reduced from 12 to 6. As shape functions for stress we take

 ~ΣK={τ∈ΣK|divhτ∈T},

so . As degrees of freedom for we take the face moments (3.2) but dispense with the interior degrees of freedom (3.3).

Let us see how the unisolvence argument adapts to these elements. If with vanishing degrees of freedom, then , and for all ,

 ∫K(divτ)vdx=−∫Kτ:ϵ(v)dx+∫∂Kτn vds=0,

using the degrees of freedom and the fact that . Thus on and for all ,

 ∫Kτ:ϵ(v)dx=−∫K(divτ)vdx+∫∂Kτn vds=0.

This shows that , so all degrees of freedom (3.3) vanish as well. Therefore the previous unisolvence result applies, and gives .

A similar argument establishes the commutativity of the projection into (the analogue of (4.2)), and the analogue of the inclusion (4.1) obviously holds. The space still contains so the approximability (4.4) still holds, but the approximability of is of one order lower, i.e., in (4.3) can be at most . As a result, the error estimates given by (4.3) in Theorem 4.3 carry over, except that is limited to in the error estimate for .

### References

1. Scot Adams and Bernardo Cockburn, A mixed finite element method for elasticity in three dimensions, J. Sci. Comput. 25 (2005), no. 3, 515–521. MR 2221175 (2006m:65251)
2. Mohamed Amara and Jean-Marie Thomas, Equilibrium finite elements for the linear elastic problem, Numer. Math. 33 (1979), no. 4, 367–383. MR 553347 (81b:65096)
3. Douglas N. Arnold and Gerard Awanou, Rectangular mixed finite elements for elasticity, Math. Models Methods Appl. Sci. 15 (2005), no. 9, 1417–1429. MR 2166210 (2006f:65112)
4. Douglas N. Arnold, Gerard Awanou, and Ragnar Winther, Finite elements for symmetric tensors in three dimensions, Math. Comp. 77 (2008), no. 263, 1229–1251. MR 2398766 (2009b:65291)
5. Douglas N Arnold and Franco Brezzi, Mixed and nonconforming finite element methods: implementation, postprocessing and error estimates, RAIRO-M2AN Modelisation Math et Analyse 19 (1985), no. 1, 7–32. MR 813687 (87g:65126)
6. Douglas N. Arnold, Franco Brezzi, and Jim Douglas, Jr., PEERS: a new mixed finite element for plane elasticity, Japan J. Appl. Math. 1 (1984), no. 2, 347–367. MR 840802 (87h:65189)
7. Douglas N. Arnold, Jim Douglas, Jr., and Chaitan P. Gupta, A family of higher order mixed finite element methods for plane elasticity, Numer. Math. 45 (1984), no. 1, 1–22. MR 761879 (86a:65112)
8. Douglas N. Arnold, Richard S. Falk, and Ragnar Winther, Finite element exterior calculus, homological techniques, and applications, Acta Numer. 15 (2006), 1–155. MR 2269741 (2007j:58002)
9. by same author, Mixed finite element methods for linear elasticity with weakly imposed symmetry, Math. Comp. 76 (2007), no. 260, 1699–1723 (electronic). MR 2336264 (2008k:74057)
10. Douglas N. Arnold and Ragnar Winther, Mixed finite elements for elasticity, Numer. Math. 92 (2002), no. 3, 401–419. MR 1930384 (2003i:65103)
11. by same author, Nonconforming mixed elements for elasticity, Math. Models Methods Appl. Sci. 13 (2003), no. 3, 295–307, Dedicated to Jim Douglas, Jr. on the occasion of his 75th birthday. MR 1977627 (2004f:65176)
12. Gerard Awanou, A rotated nonconforming rectangular mixed element for elasticity, Calcolo 46 (2009), no. 1, 49–60. MR 2495247 (2010c:74057)
13. by same author, Symmetric matrix fields in the finite element method, Symmetry 2 (2010), no. 3, 1375–1389. MR 2804836 (2012e:74013)
14. Shao-Chun Chen and Ya-Na Yang, Conforming rectangular mixed finite elements for elasticity, Journal of Scientific Computing 47 (2010), no. 1, 93–108. MR 2804836 (2012e:74013)
15. Bernardo Cockburn, Jayadeep Gopalakrishnan, and Johnny Guzmán, A new elasticity element made for enforcing weak stress symmetry, Math. Comp. 79 (2010), no. 271, 1331–1349. MR 2629995 (2011m:65276)
16. Martin A. Eisenberg and Lawrence E. Malvern, On finite element integration in natural co-ordinates, International Journal for Numerical Methods in Engineering 7 (1973), no. 4, 574–575.
17. Richard S. Falk, Finite element methods for linear elasticity, Mixed finite elements, compatibility conditions, and applications (Daniele Boffi and Lucia Gastaldi, eds.), Lecture Notes in Mathematics, vol. 1939, Springer-Verlag, Berlin, 2008, Lectures given at the C.I.M.E. Summer School held in Cetraro, June 26–July 1, 2006. MR 2459075 (2010h:65219)
18. Jayadeep Gopalakrishnan and Johnny Guzmán, Symmetric nonconforming mixed finite elements for linear elasticity, SIAM J. Numer. Anal. 49 (2011), no. 4, 1504–1520. MR 2831058
19. by same author, A second elasticity element using the matrix bubble, IMA J. Numer. Anal. 32 (2012), no. 1, 352–372. MR 2875255
20. Johnny Guzmán, A unified analysis of several mixed methods for elasticity with weak stress symmetry, J. Sci. Comput. 44 (2010), no. 2, 156–169. MR 2659794 (2011h:74021)
21. Jun Hu and Zhong-Ci Shi, Lower order rectangular nonconforming mixed finite elements for plane elasticity, SIAM J. Numer. Anal. 46 (2007/08), no. 1, 88–102. MR 2377256 (2008m:65321)
22. Claes Johnson and Bertrand Mercier, Some equilibrium finite element methods for two-dimensional elasticity problems, Numer. Math. 30 (1978), no. 1, 103–116. MR 0483904 (58 #3856)
23. Hong-Ying Man, Jun Hu, and Zhong-Ci Shi, Lower order rectangular nonconforming mixed finite element for the three-dimensional elasticity problem, Math. Models Methods Appl. Sci. 19 (2009), no. 1, 51–65. MR 2484491 (2009m:74012)
24. Mary E. Morley, A family of mixed finite elements for linear elasticity, Numer. Math. 55 (1989), no. 6, 633–666. MR 1005064 (90f:73006)
25. Serge Nicaise, Regularity of the solutions of elliptic systems in polyhedral domains, Bull. Belg. Math. Soc. Simon Stevin 4 (1997), no. 3, 411–429. MR 1457079 (98k:35044)
26. Rolf Stenberg, On the construction of optimal mixed finite element methods for the linear elasticity problem, Numer. Math. 48 (1986), no. 4, 447–462. MR 834332 (87i:73062)
27. by same author, A family of mixed finite elements for the elasticity problem, Numer. Math. 53 (1988), no. 5, 513–538. MR 954768 (89h:65192)
28. by same author, Two low-order mixed methods for the elasticity problem, The mathematics of finite elements and applications, VI (Uxbridge, 1987), Academic Press, London, 1988, pp. 271–280. MR 956898 (89j:73074)
29. Son-Young Yi, Nonconforming mixed finite element methods for linear elasticity using rectangular elements in two and three dimensions, Calcolo 42 (2005), no. 2, 115–133. MR 2158594 (2006h:74071)
30. by same author, A new nonconforming mixed finite element method for linear elasticity, Math. Models Methods Appl. Sci. 16 (2006), no. 7, 979–999. MR 2250030 (2007e:65127)