Numer. Math.

(2002) 92: 401–419

Digital Object Identifier (DOI) 10.1007/s002110100348 Numerische
Mathematik

Mixed finite elements for elasticity

Douglas N. Arnold1 , Ragnar Winther2
1
Institute for Mathematics and its Applications, University of Minnesota, Minneapolis,
MN 55455, USA; e-mail: [email protected]
2
Department of Informatics, University of Oslo, Oslo, Norway; e-mail: rwinther@ifi.uio.no

Received January 10, 2001 /

Published online November 15, 2001 –

c Springer-Verlag 2001

Summary. There have been many efforts, dating back four decades, to
develop stable mixed finite elements for the stress-displacement formulation
of the plane elasticity system. This requires the development of a compatible
pair of finite element spaces, one to discretize the space of symmetric tensors
in which the stress field is sought, and one to discretize the space of vector
fields in which the displacement is sought. Although there are number of
well-known mixed finite element pairs known for the analogous problem
involving vector fields and scalar fields, the symmetry of the stress field is
a substantial additional difficulty, and the elements presented here are the
first ones using polynomial shape functions which are known to be stable.
We present a family of such pairs of finite element spaces, one for each
polynomial degree, beginning with degree two for the stress and degree one
for the displacement, and show stability and optimal order approximation.
We also analyze some obstructions to the construction of such finite element
spaces, which account for the paucity of elements available.
Mathematics Subject Classification (1991): 65N30, 74S05

1. Introduction

Let σ and u denote the stress and displacement fields engendered by a body
force f acting on a linearly elastic body which occupies a planar region Ω
and which is clamped on ∂Ω. Then σ takes values in the space S = R2×2 sym

The first author was supported by NSF grant DMS-9870399 and the second author by the
Norwegian Research Council (NFR) under grant 135420/431.
Correspondence to: D.N. Arnold
402 D.N. Arnold, R. Winther

of symmetric tensors and u in R2 , and the pair (σ, u) is characterized as the

unique critical point of the Hellinger–Reissner functional


1 
(1.1) J (τ, v) = Aτ : τ + div τ · v − f · v dx.
Ω 2

Here the compliance tensor A = A(x) : S → S is bounded and symmetric

positive definite uniformly for x ∈ Ω, and the critical point is sought among
all τ ∈ H(div, Ω, S), the space of square-integrable symmetric matrix
fields with square-integrable divergence, and all v ∈ L2 (Ω, R2 ), the space
of square-integrable vector fields. Equivalently, (σ, u) ∈ H(div, Ω, S) ×
L2 (Ω, R2 ) is the unique solution to the following weak formulation of the
elasticity system:



(Aσ : τ + div τ · u + div σ · v) dx = f v dx,
Ω Ω
(1.2) (τ, v) ∈ H(div, Ω, S) × L2 (Ω, R2 ).

A mixed finite element method determines an approximate stress field

σh and an approximate displacement field uh as the critical point of J over
Σh × Vh where Σh ⊂ H(div, Ω, S) and Vh ⊂ L2 (Ω, R2 ) are suitable
piecewise polynomial subspaces. Equivalently, the pair (σh , uh ) ∈ Σh × Vh
is determined by the weak formulation (1.2), with the test space restricted
to Σh × Vh . As is well known, the subspaces Σh and Vh cannot be chosen
arbitrarily. To ensure that a unique critical point exist and that it provides
a good approximation of the true solution, they must satisfy the stability
conditions from the theory of mixed methods [6, 7]:
(A1) There exists a positive constant
 c1 such that τ H(div) ≤ c1 τ L2
whenever τ ∈ Σh satisfies Ω div τ · v dx = 0 for all v ∈ Vh .
(A2) There exists a positive constant c2 such that for all nonzero v ∈ Vh
there exists nonzero τ ∈ Σh with Ω div τ ·v dx ≥ c2 τ H(div) vL2 .
Condition (A2) is one of the two stability conditions of [6], while (A1)
implies the other, and essentially equivalent to it in practice. However we
shall establish a second set of conditions, which imply (A1) and (A2) and
also some other useful properties of the mixed method. These conditions
are:
(A1 ) div Σh ⊂ Vh .
(A2 ) There exists a linear operator Πh : H 1 (Ω, S) → Σh , bounded in
L(H 1 , L2 ) uniformly with respect to h, and such that div Πh σ =
Ph div σ for all σ ∈ H 1 (Ω, S), where Ph : L2 (Ω, R2 ) → Vh denotes
the L2 -projection.
Mixed finite elements for elasticity 403

It is clear that condition (A1 ) implies (A1) with c1 = 1. In order to see that
the so-called commuting diagram property (A2 ) implies (A2) we have to
invoke the fact that for each v ∈ L2 (Ω, R2 ) there exists a τ ∈ H 1 (Ω, S)
such that
div τ = v, τ H 1 ≤ c0 vL2 ,
where c0 is independent of v (this is discussed briefly in the following
section). For v ∈ Vh and with this choice of τ we have Πh τ ∈ Σh ,
div Πh τ = Ph div τ = Ph v = v, and


div Πh τ · v dx = v2L2 ≥ c−1
0 τ H 1 vL2 ≥ c2 Πh τ H(div) vL2 ,
Ω

where c2 = c0−1 (1 + Πh 2L(H 1 ,L2 ) )−1/2 .

Despite four decades of effort, very few choices of finite element spaces
have been constructed which satisfy conditions (A1) and (A2). In fact, the
only ones known use composite elements, in which Vh consists of piece-
wise polynomials with respect to one triangulation of the domain, while Σh
consists of piecewise polynomials with respect to a different, more refined,
triangulation [4, 12, 13, 19]. Because of the lack of suitable mixed elasticity
elements, several authors have resorted to the use of Lagrangian functionals
which are modifications of the Hellinger–Reissner functional given above
[1,3,5,15–18], in which the symmetry of the stress tensor is enforced only
weakly or abandoned altogether. Another modification is analyzed in [14],
where the solution space is altered such that Σh is only required to be a
subspace of L2 (Ω, S).
In this paper we present a family of mixed finite element spaces for the
unmodified Hellinger–Reissner formulation. The spaces consist of piece-
wise polynomials with respect to a single arbitrary triangular subdivision of
Ω and satisfy conditions (A1 ) and (A2 ) with constants c1 and c2 indepen-
dent of the triangulation (assuming, as usual, uniform shape regularity of
the triangulations). The space Vh we use to approximate the displacement
simply consists of all piecewise polynomials of degree k with no interele-
ment continuity constraints. The degree k may be any positive integer (with
different values of k determining different finite element spaces on the same
triangulation). The space Σh is, of course, more complicated. Restricted to a
single simplex the elements of Σh consist of all piecewise polynomial matrix
fields of degree at most k + 1 together with the divergence-free matrix fields
of degree k + 2. Degrees of freedom are specified on each element which
determine the interelement continuity and ensure that Σh ⊂ H(div, Ω, S).
In order to present the main ideas as clearly as possible, we shall first
analyze the lowest order family of elements. After some preliminaries in
Sect. 2, in Sect. 3 we construct the spaces Σh and Vh in the case k = 1, by
404 D.N. Arnold, R. Winther

describing their restrictions to each triangle and giving a unisolvent set of

local degrees of freedom. We establish a relation between these space and
the Hermite quintic C 1 finite element space, and explain in general the rela-
tionship between stable finite elements for the Hellinger–Reissner principle
and C 1 finite elements. This provides a major obstruction to the construc-
tion of the former. Our element involves vertex degrees of freedom, a fact
which, as we shall show, is unavoidable. Because of this, the interpolation
operator into Σh associated with the degrees of freedom is not bounded on
H 1 (Ω, S), and so, even though it satisfies the commutativity property of
condition (A2 ), it does not establish that condition. In the following section
we modify the interpolation operator, maintaining the commutativity, and
establish its boundedness and approximation properties. Then in Sect. 5 we
complete the error analysis for the k = 1 case. In Sect. 6, we briefly describe
the case k ≥ 1, which is altogether analogous to the case k = 1. In a few
brief final paragraphs we describe a variant of the lowest order element with
fewer degrees of freedom, and mention some simple extensions.

2. Notation and preliminaries

We denote by H k (T, X) the Sobolev space consisting of functions with

domain T ⊂ R2 , taking values in the finite-dimensional vector space X, and
with all derivatives of order at most k square-integrable. For our purposes,
the range space X will be either S, R2 , or R. In the latter case we may
write simply H k (T ). We will generally write  · k or  · H k instead of
 · H k (T,X) . We similarly denote by Pk (T, X) the space of polynomials on
T with degree at most k.
If τ is a symmetric matrix field then its divergence, div τ , is the vector
field obtained by applying the ordinary divergence operator to each row. The
symmetric part of the gradient of a vector field v, denoted  v, is given by
 v = [grad v + (grad v)T ]/2.
Throughout the paper we assume that the elastic domain Ω is a simply
connected polygonal domain in R2 . Any smooth vector field on Ω may be
realized as the divergence of a smooth symmetric matrix field. E.g., we may
extend the vector field smoothly to a larger smoothly bounded domain, and
then solve the equations of elasticity there with the extended vector field as
body forces. The same argument shows that any vector field in L2 (Ω, R2 )
may be realized as the divergence of a matrix field in H 1 (Ω, S), a fact we
already invoked in the introduction and will use as well in the sequel. Further,
a symmetric matrix field τ on a simply connected domain is divergence-free
if and only if τ admits a potential, or Airy stress function,
 2 
∂ q/∂y 2 −∂ 2 q/∂x∂y
τ = Jq := .
−∂ 2 q/∂x∂y ∂ 2 q/∂x2
Mixed finite elements for elasticity 405

The potential q is determined by τ up to addition of a linear polynomial.

These considerations are summarized by the statement that the following se-
quence is exact (i.e., that the range of each map is the kernel of the following
one):
⊂
0 −−−→ P1 (Ω) −−−→ C ∞ (Ω)
J div
−−−→ C ∞ (Ω, S) −−−→ C ∞ (Ω, R2 ) −−−→ 0.
This exact sequence is related, although rather indirectly, to the de Rham
sequence for the domain the Ω [11]. We have stated it in terms of infinitely
differentiable functions, but analogous results hold with less smoothness.
E.g., the sequence
⊂ J
0 −−−→ P1 (Ω) −−−→ H 2 (Ω) −−−→ H(div, Ω, S)
(2.1)
div
−−−→ L2 (Ω, R2 ) −−−→ 0
is also exact. There is a polynomial analogue of the sequence as well: for
any integer k ≥ 0 the sequence
⊂ J
0 −−−→ P1 (Ω) −−−→ Pk+3 (Ω) −−−→ Pk+1 (Ω, S)
div
−−−→ Pk (Ω, R2 ) −−−→ 0
is exact. (To verify the surjectivity of the final divergence, it suffices to count
dimensions and use the exactness of the sequence at the other points.) As we
shall see a key ingredient in the development and analysis of mixed methods
for elasticity is a discrete analogue of the exact sequence (2.1).

3. The finite element method in the lowest order case

In this section we present a mixed finite element method based on piecewise

linear displacements and piecewise quadratic stresses, the latter augmented
by some cubic shape functions. First we describe the finite elements on a
single triangle T ⊂ Ω. Define
ΣT = P2 (T, S) + { τ ∈ P3 (T, S) | div τ = 0 }
(3.1) = { τ ∈ P3 (T, S) | div τ ∈ P1 (T, R2 )},
VT = P1 (T, R2 ).
The space VT has dimension 6 and a complete set of degrees of freedom are
given by the value of the two components at the three nodes interior to T . The
space ΣT clearly has dimension at least 24, since the dim P3 (T, S) = 30
and the condition that div τ ∈ P1 (T, R2 ) represents six linear constraints.
406 D.N. Arnold, R. Winther

We now exhibit 24 degrees of freedom ΣT → R and show that they vanish

simultaneously only when τ = 0. This implies that the dimension of ΣT
is precisely 24 (which could also be established directly using the fact that
div P3 (Ω, S) = P2 (Ω, R2 )), and that the degrees of freedom are unisolvent.
The degrees of freedom are
– the values of three components of τ (x) at each vertex x of T (9 degrees
of freedom)
– the values of the moments of degree 0 and 1 of the two normal compo-
nents of τ on each edge e of T (12 degrees of freedom)
– the value of the three components of the moment of degree 0 of τ on T
(3 degrees of freedom)
Otherwise stated, we  determine τ ∈ ΣT by giving its values at the ver-
tices,
 the values of e (τ n) ds and e (τ n)s ds for all edges, and the value of
T τ dx. (Here s is a parameter giving the distance to one of the end points
of e and n is one of the unit vectors normal to e.) Note that the degrees of
freedom associated to an edge and its end points determine τ n on that edge.
This is just the condition required to obtain a conforming approximation of
H(div, Ω, S). The element diagrams in Fig. 1 are mnemonic of the degrees
of freedom.

Fig. 1. Element diagrams for the lowest order stress and displacement elements

Lemma 3.1. If the 24 degrees of freedom just given all vanish for some
τ ∈ ΣT , then τ = 0.
Proof. We immediately have that τ n vanishes on each edge. Letting v =
div τ , a linear vector field on T , we get




2
(3.2) v dx = − τ :  v dx + τ n · v ds = 0
T T ∂T

since the integral of τ vanishes as well as τ n. Thus τ is divergence-free

and hence τ = Jq for some q ∈ P5 (T ). Adjusting by a linear function we
may take q to vanish at the vertices. Now ∂ 2 q/∂s2 = τ n · n = 0 on each
edge, whence q is identically zero on ∂T . This implies that the gradient of
Mixed finite elements for elasticity 407

q vanishes at the vertices. Since ∂ 2 q/∂s∂n = −τ n · t = 0 on each edge

(with t a unit vector tangent to the edge), we conclude that ∂q/∂n vanishes
identically on ∂T as well. Since q has degree at most 5, it must vanish
identically. 

Having given a unisolvent set of degrees of freedom for VT and ΣT , our
finite element space is assembled in the usual way. Let Th be some triangu-
lation of Ω, i.e., a set of closed triangles with union Ω̄ and such that any two
distinct non-disjoint elements of Th meet in a common edge or vertex. The
associated finite element element space Vh is then the space of all piece-
wise linear vector fields with respect to this triangulation, not subject to any
interelement continuity conditions. The space Σh is the space of all matrix
fields which belong piecewise to ΣT , subject to the continuity conditions
that the normal components are continuous across mesh edges and all com-
ponents are continuous at mesh vertices. The first condition is necessary
to ensure that Σh ⊂ H(div, Ω, S), but the second condition is a further
restriction not implied by the inclusion in H(div, Ω, S). (This is analogous
to the condition of continuity of second derivatives at mesh vertices for the
Hermite quintic C 1 finite element. We shall see below that these two phe-
nomena are in fact related, and we shall argue that the restriction to vertex
continuity is unavoidable.) Note that the stability condition (A1 ) holds by
construction.
The global degrees of freedom for the assembled finite element space
Σh are the values of all three components at all the mesh vertices, the values
of the moments of degrees 0 and 1 of the normal components on all the
mesh edges, and the values of the moments of degree 0 for all components
on all the mesh triangles. These functionals extend naturally to C(Ω, S) and
so determine a canonical interpolation operator C(Ω, S) → ΣT . However,
because of the vertex degrees of freedom, the canonical interpolation op-
erator is not bounded with respect to the norm in H 1 (Ω, S) and so cannot
be used to establish the stability condition (A2 ). Therefore we define an
alternate interpolation operator which is bounded on H 1 (Ω, S). To do so,
we require bounded linear operators Ehx : H 1 (Ω, S) → S for each vertex x
of the triangulation. Given such operators, we define Πh : H 1 (Ω, S) → Σh
by
(3.3) Πh τ (x) = Ehx τ for all vertices x,


(3.4) (τ − Πh τ )n · v ds = 0 for all edges e and all v ∈ P1 (e, R2 ),
e


(3.5) (τ − Πh τ ) dx = 0 for all triangles T .
T
IfEhxwere simply the evaluation operator τ → τ (x), we would obtain the
canonical interpolation operator, but this choice of Ehx does not fulfill the
408 D.N. Arnold, R. Winther

required boundedness on H 1 (Ω, S). A choice that is certainly bounded is

simply Ehx τ = 0 for all τ and all x, and in fact this choice is sufficient
for verifying the stability condition (A2 ). We shall denote the resulting
interpolation operator by Πh0 . However, in the next section we shall choose
Ehx τ to be an approximate evaluation operator (a weighted average of τ
near x), and in that way obtain an interpolation operator Πh with better
approximation properties.
We have for any τ ∈ H 1 (Ω, S), and T ∈ Th , and any v ∈ VT that


div(τ − Πh τ ) · v dx
T



= − (τ − Πh τ ) :  v dx + (τ − Πh τ )n · v ds.
T ∂T

The right hand side vanishes in view of (3.5) and (3.4). This verifies the
commutativity property

(3.6) div Πh τ = Ph div τ,

where Ph : L2 (Ω, R2 ) → Vh is the orthogonal projection. We note that

this property holds no matter how the choice of the Ehx are made. A useful
consequence of (3.6) is that div Σh = Vh . Indeed, given any v ∈ Vh we
may find τ ∈ H 1 (Ω, S) such that div τ = v, and then Πh τ ∈ Σh satisfies
div Πh τ = Ph v = v.
Let Qh = { q ∈ H 2 (Ω) | Jq ∈ Σh }. It is easy to identify Qh . Its el-
ements are piecewise quintic polynomials and belong globally to C 1 (Ω).
Moreover, they are C 2 at the vertices of the mesh, since Jq ∈ Σh . Further,
any piecewise quintic with these continuity properties is mapped by J into
Σh , so belongs to Qh . We conclude that Qh is precisely the space of C 1
piecewise quintics which are C 2 at the vertices, that is, the well-known Her-
mite quintic or Argyris finite element [2]; cf., also [8, § 9]. The relationship
between Qh and Σh is even more intimate. Define a projection operator
Ih : C ∞ (Ω) → Qh by requiring that the vertex values of Ih q, the vertex
values of grad Ih q, and the edge moments of degree 0 of ∂(Ih q)/∂n all
be equal to the corresponding values for q, and that the Hessian of Ih q at
each vertex h be given by JIh q = Ehx Jq. Then the commutativity property
JIh q = Πh Jq can be verified easily. All these considerations are summa-
rized in the following commutative diagram with exact rows:

⊂ J div
0 −−−−−→ P1 (Ω) −−−−−→ C ∞ (Ω) −−−−−→ C ∞ (Ω, S) −−−−−→ C ∞ (Ω, R2 ) −−−−−→ 0
   
 I Π P
id  h  h  h
⊂ J div
0 −−−−−→ P1 (Ω) −−−−−→ Qh −−−−−→ Σh −−−−−→ Vh −−−−−→ 0
Mixed finite elements for elasticity 409

It is instructive to examine which properties of our H(div, Ω, S) finite el-

ement space Σh led us to an H 2 finite element space. Suppose we have
any finite element space Σh ⊂ H(div, Ω, S) and an interpolation operator
Πh : C ∞ (Ω, S) → Σh (e.g., the interpolation operation determined by a
set of degrees of freedom for Σh ) satisfing:
(a) Πh τ is divergence-free whenever τ is divergence-free,
(b) (Π
 h τ )n on any edge  by the restriction of τ n to that edge,
 is determined
(c) e (Πh τ )n ds = e τ n ds and e (Πh τ )n · n s ds = e τ n · n s ds.
The first property is satisfied if a commutativity property like (3.6) holds,
and the second is necessary if the interpolation operator maps into the space
H(div, Ω, S) and is determined by local degrees of freedom. The third prop-
erty is usually required to verify the commutativity property (3.6). Under
these assumptions we can define a corresponding subspace Qh of H 2 (Ω)
as the inverse image of Σh under J (so, in particular, the elements of Qh are
piecewise polynomials of degree two greater than those of Σh ), and define
Ih : C ∞ (Ω) → Qh as follows. Given q ∈ C ∞ (Ω), Πh Jq is divergence-free
by the second property, and so is equal to Jr for some function r determined
up to addition of a linear polynomial. Thus for any specified vertex x we de-
termine Ih q uniquely by the conditions that JIh q = Πh Jq, Ih q(x) = q(x),
∇Ih q(x) = ∇q(x). Now (Jq)n·n = ∂ 2 q/∂s2 and (Jq)n·t = −∂ 2 q/∂s∂n.
Since Ih q preserves the moments of Jq indicated in property (c), and it also
preserves the value and the gradient of q at one vertex, we can integrate on
edges to conclude that Ih q = q and ∇Ih q = ∇q at all vertices. We can
then use property (b) in a similar way to show that Ih q and ∂Ih q/∂n is
determined on any edge by the the restriction of q and ∂q/∂n on that edge.
This gives us the essential ingredients of a conforming H 2 finite element
space determined by local degrees of freedom. Roughly speaking, we have
shown that whenever one can construct an H(div, Ω, S) finite element space
Σh with expected properties, one can construct an H 2 finite element space
related by Qh = J −1 (Σh ). Our element is so related to the Hermite quintic
element, the element of [13] is so related to the Clough–Tocher composite
H 2 element, and the element family of [4] is so related to the higher order
composite H 2 elements of [10].
The relationship between H(div, Ω, S) finite elements and H 2 finite
elements just outlined presents a major obstruction toward the construction
of the former, since the difficulty in constructing conforming H 2 elements
is well-known. There is an analogous obstruction to the development of
H(div, Ω, R2 ) finite elements, but it is, in contrast, very minor. It requires
only that the inverse image of such a space under the curl operator must
be a conforming finite element discretization of H 1 (Ω), as, for example,
the Lagrange finite elements of degree k + 1 are the inverse image of the
Raviart–Thomas elements of order k.
410 D.N. Arnold, R. Winther

The finite element method for elasticity based on our spaces Σh and Vh
has many features in common with popular mixed finite element methods
for scalar elliptic problems, and we shall use these in Sect. 5 to give an
error analysis. However there are some notable differences as well. One of
these is the presence of vertex nodes. This leads to technical complications
in the analysis, since it adds an additional regularity requirement for func-
tions to belong to the domain of Πh . However, as we now explain, vertex
nodes are unavoidable whenever continuous shape functions are used to con-
struct an H(div, Ω, S) finite element space. To see why, imagine building an
H(div, Ω, S) finite element space from spaces ΣT of continuous symmetric
matrix fields, imposing interelement continuity only by means of quantities
defined on the edges (and so shared by only two neighboring elements).
This means that the degrees of freedom associated with each edge must
determine the normal component on the edge. Now consider two edges of a
triangle meeting at a common vertex x. If n1 is the normal to the first edge
and n2 the normal to the second edge then the degrees of freedom on the
first edge must determine τ n1 there and similarly the degrees of freedom
on the second edge must determine τ n2 there. Since τ is continuous, we
have in particular that the degrees of freedom on the first edge determine
τ (x)n1 · n2 , and those on the second edge determine τ (x)n2 · n1 . But these
quantities are equal since τ is symmetric. This is a contradiction, since the
degrees of freedom on the two edges are necessarily independent.
This argument indicates that it is at least necessary to take the quantity
τ (x)n1 · n2 as a degree of freedom associated to the vertex node x. But
the node x is shared by other triangles which will have other values for
the edge normal vectors ni . For this reason, except if we restrict to very
special triangulations, we are forced to take all three components of τ (x)
as degrees of freedom associated to the node x1 . Note that the composite
H(div, Ω, S) finite elements of [13] and [4] avoid the necessity of vertex
degrees of freedom because they use discontinuous shape functions.

4. Approximation properties

In this section we analyze the approximation properties of the interpolation

operator Πh : H 1 (Ω, S) → Σh . Recall that this operator is defined by (3.3)–
(3.5), where Ehx is a suitable approximate evaluation operator. We start by
specifying our choice of Ehx precisely. Let Σh = Σh ∩ H 1 (Ω, S) be the
space of continuous piecewise quadratic symmetric matrix fields, and let
Rh : L2 (Ω, S) → Σh ⊂ Σh be a Clement interpolant [9] satisfying
1
A variant of this argument can be used to show that an H 2 finite element which uses C 1
shape functions must include second derivatives at the vertices among its degrees of freedom
Mixed finite elements for elasticity 411

(4.1) Rh τ − τ j ≤ chm−j τ m , 0 ≤ j ≤ 1, j ≤ m ≤ 3,

where h denotes the mesh size (i.e., the maximum diameter of a triangle in
Th ), and the constant c depends only on the shape regularity of the triangu-
lation. We then specify the interpolation operator Πh via (3.3)–(3.5) with
Ehx τ = Rh τ (x).
Recall also that Πh0 : H 1 (Ω, S) → Σh denotes the simpler interpolation
operator obtained with the choice Ehx τ = 0. It is straightforward to verify
the identity

Πh = Πh0 (I − Rh ) + Rh .

Hence the error in Πh , is given by

(4.2) I − Πh = (I − Πh0 )(I − Rh ).

The operator Πh0 is completely local with respect to the triangulation Th

and we may analyze it triangle-by-triangle, mapping each triangle to a fixed
reference element and applying standard scaling arguments. We can then
establish error estimates for the interpolation operator Πh by combining
estimates for Πh0 with the estimates (4.1) for the Clement interpolant.
Let F : T̂ → T be an affine isomorphism of the form F x̂ = B x̂ + b.
Given τ̂ : T̂ → S, define τ : T → S by the matrix Piola transform

(4.3) τ (x) := B τ̂ (x̂)B T ,

where x = F x̂. Clearly this sets up a one-to-one correspondence between

L2 (T, S) and L2 (T̂ , S). A direct computation also shows that

(4.4) div τ (x) = B div τ̂ (x̂).

It follows that τ ∈ H(div, T, S) if and only if τ̂ ∈ H(div, T̂ , S) and that

τ ∈ ΣT if and only if τ̂ ∈ ΣT̂ . Indeed, ΣT is the space of τ polynomial of
degree 3 such that div τ is polynomial of degree 1, and, by (4.3) and (4.4),
these properties clearly transform.
Let ΠT0 : H 1 (T, S) → ΣT be the restriction of Πh0 to a single triangle,
i.e.,
– ΠT0 τ (x) = 0 for all vertices x of T ,
– e (τ − Πh0 τ )n · v ds = 0 for all edges e of T , and v ∈ P1 (e, R2 ),
– T (τ − Πh0 τ ) dx = 0.
We claim that

(4.5) ΠT0 τ (x) = BΠT̂0 τ̂ (x̂)B T ,

412 D.N. Arnold, R. Winther

where τ and τ̂ are related by (4.3) and x = F x̂. To verify (4.5) we first note
that
BΠT̂0 τ̂ (x̂)B T = 0
for each vertex of T . Furthermore,



0 T
BΠT̂ τ̂ (x̂)B dx = (det B)B ΠT̂0 τ̂ (x̂) dx̂B T
T T̂


= (det B)B τ̂ (x̂) dx̂B T
T̂


= τ (x) dx.
T

Hence, it only remains to verify that BΠT̂0 τ̂ (x̂)B T has the edge moments
required of ΠT0 τ (x). Let e be an edge of T , ê the corresponding edge of T̂ ,
and let v ∈ P1 (e, R2 ). Since B T ne is normal to ê it follows that



0 T |e|
[BΠT̂ τ̂ (x̂)B ]ne · v(x) ds = B ΠT̂0 τ̂ (x̂)B T ne · v̂(x̂) dŝ
e |ê|

ê
|e|
= B τ̂ (x̂)B T ne · v̂(x̂) dŝ
|ê| ê


= τ (x)ne · v(x) ds,
e

where v̂(x̂) = v(x). The equation (4.5) is therefore verified.

The operator ΠT̂0 is bounded from H 1 (T̂ , S) to L2 (T̂ , S). Therefore, a
standard scaling argument using (4.5) gives
(4.6) Πh0 τ 0 ≤ c(τ 0 + hτ 1 )
where the constant c depends only on the shape regularity of the triangula-
tion. From (4.1) and (4.6) we obtain
Πh0 (I − Rh )τ 0 ≤ c((I − Rh )τ 0 + h(I − Rh )τ 1 ≤ chm τ m
for 1 ≤ m ≤ 3. Hence, it follows from (4.1) and (4.2) that
(4.7) Πh τ − τ 0 ≤ chm τ m , 1 ≤ m ≤ 3.
It follows in particular that
(4.8) Πh τ 0 ≤ cτ 1 ,
and so the stability property (A2 ) holds. We also note that for the orthogonal
projection operator Ph , we have the obvious error estimates
(4.9) v − Ph v0 ≤ chm vm , 0 ≤ m ≤ 2.
Mixed finite elements for elasticity 413

5. Error analysis

Having established the stability properties (A1 ) and (A2 ) for the spaces
Σh and Vh , we conclude that the Hellinger–Reissner functional has a unique
critical point over Σh × Vh and obtain the quasioptimal estimate [6, 7]

σ − σh H(div) + u − uh L2 ≤ c inf (σ − τ H(div) + u − vL2 ).

τ ∈Σh
v∈Vh

The infimum on the right hand side is easily seen to be O(h2 ) for smooth
solutions. However we shall now state and prove more precise estimates,
with an analysis very analogous to that for second order elliptic problems.

Theorem 5.1. Let (σ, u) denote the unique critical point of the Hellinger–
Reissner functional over H(div, Ω, S) × L2 (Ω, R2 ) and let (σh , uh ) denote
unique critical point over Σh × Vh , where Σh and Vh are the spaces defined
in Sect. 3. Then

σ − σh 0 ≤ chm σm , 1 ≤ m ≤ 3,
div σ − div σh 0 ≤ chm div σm , 0 ≤ m ≤ 2,
u − uh 0 ≤ chm um+1 , 1 ≤ m ≤ 2.

Proof. Recall that the exact solution (σ, u) ∈ H(div, Ω, S) × L2 (Ω, R2 ) is

the unique solution of (1.2), i.e.,


(Aσ : τ + div τ · u) dx = 0, τ ∈ H(div, Ω, S),
Ω



div σ · v dx = f · v dx, v ∈ L2 (Ω, R2 ).
Ω Ω

The mixed method solution (σh , uh ) ∈ Σh × Vh satisfies the corresponding

discrete system


(Aσh : τ + div τ · uh ) dx = 0, τ ∈ Σh ,
Ω



div σh · v dx = f · v dx, v ∈ Vh .
Ω Ω

Subtracting we obtain error equations, which, since div Σh ⊂ Vh , we may

write as


(5.1) [A(σ − σh ) : τ + div τ · (Ph u − uh )] dx = 0, τ ∈ Σh ,
Ω


(5.2) div(σ − σh ) · v dx = 0, v ∈ Vh .
Ω
414 D.N. Arnold, R. Winther

The second equation immediately implies that

(5.3) div σh = Ph div σ = div Πh σ,
where the last equality comes from (3.6). Hence, from (4.9), we immediately
obtain the desired error bound for div σ:
(5.4)
div σ − div σh 0 = (I − Ph ) div σ0 ≤ chm div σm , 0 ≤ m ≤ 2.
Taking τ = Πh σ − σh in (5.1) and invoking (5.3) we obtain


A(σ − σh ) : (Πh σ − σh ) dx
Ω


= div(Πh σ − σh ) · (Ph u − uh ) dx = 0,
Ω

from which it easily follows that

σ − σh A ≤ σ − Πh σA ,

where τ 2A := Aτ : τ dx. Since the norm  · A is equivalent to the
L2 –norm it follows from (4.7) that
(5.5) σ − σh 0 ≤ chm σm , 1 ≤ m ≤ 3.
In order to establish the error estimate for the displacement, we recall
that
u − Ph u0 ≤ ch2 u2 .
Therefore, the desired estimate will follow from the bound
(5.6) Ph u − uh 0 ≤ chm um+1 1 ≤ m ≤ 2.
Let τ ∈ H 1 (Ω, S) satisfy div τ = Ph u − uh with τ 1 ≤ cPh u − uh 0 .
Then, in light of the commutativity property, (3.6), and the bound (4.8),
div Πh τ = Ph u − uh and Πh τ 0 ≤ cPh u − uh 0 . Hence,


Ph u − uh 20 = div Πh τ · (Ph u − uh ) dx
Ω


=− A(σ − σh ) : Πh τ dx
Ω
≤ cσ − σh 0 Ph u − uh 0 .
Thus
Ph u − uh 0 ≤ cσ − σh 0 ≤ chm σm ≤ chm um+1 , 1 ≤ m ≤ 3,
which contains (5.6). 

Mixed finite elements for elasticity 415

6. Higher order elements

The foregoing considerations generalize in a straightforward manner to el-

ements of higher order. Let k be a positive integer, the case k = 1 being
that considered heretofore. Generalizing (3.1), we define the finite element
spaces on a single triangle:

ΣT = Pk+1 (T, S) + { τ ∈ Pk+2 (T, S) | div τ = 0 }

= { τ ∈ Pk+2 (T, S) | div τ ∈ Pk (T, R2 )},
VT = Pk (T, R2 ).

Clearly dim VT = (k + 1)(k + 2), and

dim ΣT ≤ dk : = dim Pk+2 (T, S) − [dim Pk+1 (T, R2 )−dim Pk (T, R2 )]

= (3k 2 + 17k + 28)/2.

We shall show that in fact dim ΣT = dk and exhibit a unisolvent set of

degrees of freedom. To this end, we define

Mk (T ) = { τ ∈ Pk+2 (T, S) | div τ = 0 and τ n = 0 on ∂T }.

In the proof of Lemma 3.1, we showed that for k = 1, Mk (T ) = 0. The

same argument shows that for k ≥ 2, Mk (T ) = J(b2T Pk−2 (T )), where bT
is the cubic bubble function on T (the unique cubic polynomial achieving a
maximum value of unity on T and vanishing on ∂T ). Thus dim Mk (T ) =
(k 2 − k)/2 for k ≥ 1. The space [Pk (T, R2 )] has dimension (k + 1)(k +
2) − 3 and is clearly L2 -orthogonal to Mk (T ). Thus the space Nk (T ) :=
[Pk (T, R2 )] + Mk (T ) has dimension (3k 2 + 5k − 2)/2. The degrees of
freedom we take for ΣT are
– the values of three components of τ (x) at each vertex x of T (9 degrees
of freedom)
– the values of the moments of degree at most k of the two normal com-
ponents of τ on each edge e of T (6k + 6 degrees of freedom)
– the value of the moments T τ : φ dx, φ ∈ Nk (T ) ((3k 2 + 5k − 2)/2
degrees of freedom)
We have specified dk degrees of freedom. Thus to show unisolvence it suf-
fices to show that if all the degrees of freedom vanish for some τ ∈ Σk (T ),
then τ vanishes. The proof is just as for Lemma 3.1. From the vanishing
of the first two sets of degrees of freedom we  conclude that τ n vanishes
on ∂T . This, together with the vanishing of T τ : φ dx, φ ∈ [Pk (T, R2 )]
implies that div τ vanishes as well, so τ ∈ Mk (T ), and, using the remain-
ing degrees of freedom, we conclude that τ vanishes identically. Element
416 D.N. Arnold, R. Winther

Fig. 2. Element diagrams for the stress and displacement elements in the case k = 2

diagrams for the case k = 2 are shown in Fig. 2. In this case N2 (T ) is the
span of P1 (T , S) and the additional matrix field J(b2T ).
Having specified finite element spaces on individual triangles and the
degrees of freedom, we have determined the finite element spaces Σh and
Vh for any triangulation. The latter space is simply the space of piecewise
polynomial vector fields of degree at most k, while the former space consists
of H(div, Ω, S) matrix fields which belong to ΣT on each T ∈ Th and
in addition are continuous at the mesh vertices. In addition to the natural
interpolation operator determined by the degrees of freedom, we define
Πh : H 1 (Ω, S) → Σh by using as degrees of freedom at the vertices the
values of a suitable Clement interpolant. In this way we obtain an operator
which satisfies the commutativity property (3.6) (with Ph the L2 -projection
into Vh ), and the error estimates (4.7), now for 1 ≤ m ≤ k + 2. The analysis
generalizes directly from the case k = 1, resulting in the following analogue
of Theorem 5.1.

Theorem 6.1. Let (σ, u) denote the unique critical point of the Hellinger–
Reissner functional over H(div, Ω, S) × L2 (Ω, R2 ) and let (σh , uh ) denote
the unique critical point over Σh × Vh , where Σh and Vh are the spaces
defined above for some integer k ≥ 1. Then

σ − σh 0 ≤ chm σm , 1 ≤ m ≤ k + 2,
div σ − div σh 0 ≤ chm div σm , 0 ≤ m ≤ k + 1,
u − uh 0 ≤ chm um+1 , 1 ≤ m ≤ k + 1.

7. A simplified element of low order

There does not exist an element pair for k = 0 in our family. Indeed,
there could not, since there does not exist an H 2 finite element of degree 4.
However, there is a variant of the lowest order (k = 1 element), with fewer
degrees of freedom, which we now describe briefly. To do so we denote
by RM (T ) the space of infinitesmal rigid motions on T , i.e., the span of
P0 (T, R2 ) and the linear vector field (−x2 , x1 ). This is precisely the kernel
Mixed finite elements for elasticity 417

Fig. 3. Element diagrams for the elements of reduced order

of the symmetric gradient operator . For our simplified element we take VT

to be RM (T ) and
ΣT = { τ ∈ P3 (T, S) | div τ ∈ RM (T ) }.
Then dim ΣT = 21 and a unisolvent set of degrees of freedom are
– the values of three components of τ (x) at each vertex x of T (9 degrees
of freedom)
– the values of the moments of degree 0 and 1 of the two normal compo-
nents of τ on each edge e of T (12 degrees of freedom)
See Fig. 3. In this case P0 (T, R2 )  VT  P1 (T, R2 ) and P1 (T, S) 
ΣT  P2 (T, S). Therefore the final error estimates are
σ − σh 0 ≤ chm σm , 1 ≤ m ≤ 2,
div σ − div σh 0 ≤ chm div σm , 0 ≤ m ≤ 1,
u − uh 0 ≤ chu2 .
The analysis follows closely that of the previous element family. One dif-
ference is that the space ΣT for this element is not invariant under the Piola
transform. Consequently a different argument is required to prove the ap-
proximation properties of Σh . This can be done, for example, by scaling to
a similar element of unit diameter using translation, rotation, and dilation,
and using a compactness argument.

8. Final remarks

In the interest of clarity we have considered an elastic body clamped all

around its boundary and described by a uniformly positive definite com-
pliance tensor. Both restrictions can be loosened. The extension to traction
boundary conditions on part of the boundary is straightforward. For the
Hellinger–Reissner variational form, traction boundary conditions are es-
sential, and thus must be imposed in the stress space Σh . When traction
boundary conditions are imposed on both edges meeting at a corner, then
418 D.N. Arnold, R. Winther

the entire stress tensor must vanish there, and so all three degrees of freedom
at the corner set equal to zero. At other boundary points traction boundary
conditions imply two linear relations among three degrees of freedom. (Such
linear relation boundary conditions can be implemented by modifying the
relevant columns and rows of the unconstrained stiffness matrix, maintain-
ing symmetry.)
Another generalization that can easily be handled is the extension to
nearly incompressible or incompressible elastic materials. In the homoge-
neous isotropic case the compliance tensor is given by Aτ = [τ − λ/(2µ +
2λ) tr τ I]/2µ, where µ > 0, λ ≥ 0 are the Lamé constants. For our mixed
method, as for most methods based on the Hellinger–Reissner principle,
one can prove that the error estimates hold uniformly in λ (the incom-
pressible
 limit being λ → ∞). In the analysis above we used the fact that
Aτ : τ dx ≥ c0 τ 20 for some positive c0 . This estimate degenerates
(c0 → 0) as λ → ∞. However the estimate remains true with c0 > 0 de-
pending only on Ω and µ if we restrict τ to functions for which div τ ≡ 0,
Ω tr τ dx = 0, and this is enough to carry through the analysis. See [4]
where the details are presented for a composite mixed element.

References

1. Mohamed Amara, Jean-Marie Thomas: Equilibrium finite elements for the linear elastic
problem. Numer. Math. 33(4), 367–383 (1979)
2. John H. Argyris, Isaac Fried, Dieter W. Scharpf: The TUBA family of plate elements
for the matrix displacement method. Aero. J. Roy. Aero. Soc. 72, 514–517 (1968)
3. Douglas N. Arnold, Franco Brezzi, Jim Douglas, Jr.: PEERS: a new mixed finite element
for plane elasticity. Japan J. Appl. Math. 1(2), 347–367 (1984)
4. Douglas N. Arnold, Jim Douglas, Jr., Chaitan P. Gupta: A family of higher order mixed
finite element methods for plane elasticity. Numer. Math. 45(1), 1–22 (1984)
5. Douglas N. Arnold, Richard S. Falk: A new mixed formulation for elasticity. Numer.
Math. 53(1-2), 13–30 (1988)
6. Franco Brezzi: On the existence, uniqueness and approximation of saddle-point prob-
lems arising from Lagrangian multipliers. Rev. Française Automat. Informat. Recherche
Opérationnelle Sér. Rouge 8(R-2), 129–151 (1974)
7. Franco Brezzi, Michel Fortin: Mixed and Hybrid Finite Element Methods. New York:
Springer-Verlag, 1991
8. Philippe G. Ciarlet: The finite element method for elliptic problems. Amsterdam: North-
Holland, 1978
9. Philippe Clément: Approximation by finite element functions using local regularization.
Rev. Française Automat. Informat. Recherche Opérationnelle Sér. Rouge 9(R-2), 77–84
(1975)
10. Jim Douglas, Jr., Todd Dupont, Peter Percell, Ridgway Scott: A family of C 1 finite
elements with optimal approximation properties for various Galerkin methods for 2nd
and 4th order problems. RAIRO Anal. Numér. 13(3), 227–255 (1979)
11. Michael Eastwood: A complex from linear elasticity. Rend. Circ. Mat. Palermo (2)
Suppl. (2000), no. 63, 23–29
Mixed finite elements for elasticity 419

12. Baudoiun M. Fraejis de Veubeke: Displacement and equilibrium models in the finite
element method. Stress analysis, pp. 145–197 (O.C Zienkiewics, G.S. Holister, eds.),
New York, Wiley, 1965
13. Claes Johnson, Bertrand Mercier: Some equilibrium finite element methods for two-
dimensional elasticity problems. Numer. Math. 30(1), 103–116 (1978)
14. Alain L. Mignot, Claude Surry: A mixed finite element family in plane elasticity. Appl.
Math. Modelling 5, 259–262 (1981)
15. Erwin Stein, Raimund Rolfes: Mechanical conditions for stability and optimal conver-
gence of mixed finite elements for linear plane elasticity. Comput. Methods Appl. Mech.
Engrg. 84(1), 77–95 (1990)
16. Rolf Stenberg: On the construction of optimal mixed finite element methods for the
linear elasticity problem. Numer. Math. 48(4), 447–462 (1986)
17. Rolf Stenberg: A family of mixed finite elements for the elasticity problem. Numer.
Math. 53(5), 513–538 (1988)
18. Rolf Stenberg: Two low-order mixed methods for the elasticity problem. The mathemat-
ics of finite elements and applications, VI (Uxbridge, 1987), Academic Press, London,
1988, pp. 271–280
19. Vernon B. Watwood Jr., B. J. Hartz: An equilibrium stress field model for finite element
solution of two–dimensional elastostatic problems. Internat. Jour. Solids and Structures
4, 857–873 (1968)