Symmetric Matrix Fields in the Finite Element Method

The theory of elasticity is used to predict the response of a material body subject to applied forces. In the linear theory, where the displacement is small, the stress tensor which measures the internal forces is the variable of primal importance. However the symmetry of the stress tensor which expresses the conservation of angular momentum had been a challenge for finite element computations. We review in this paper approaches based on mixed finite element methods.


Introduction
The unknowns in the elasticity equations are the stress, which encodes the internal forces, and the displacement, which measures how points in the body move.In the linear theory of elasticity, displacements are assumed to be small and hence the stress becomes the primary unknown.The outline of the paper is as follows.The next section has four parts in which we introduce the linear theory of elasticity following [1], the concept of finite element discretizations and some notation.Then, we consider several variational formulations of the elasticity equations and make the case of the use of mixed finite elements.The main section is devoted to a survey of mixed finite elements for elasticity.We explain the Arnold-Winther construction and its generalizations to three-dimensional elements and nonconforming elements.We conclude with some remarks.

Linear Theory of Elasticity
When forces are applied to a solid body, occupying a domain Ω in R n , it will undergo deformations: its shape and volume will change.A point initially at a position x with coordinates x i , i = 1, . . .n will move to a new position x ′ i , i = 1 . . .n.The vector u with components x i − x ′ i is called the displacement vector.The deformation of the body is completely determined if u, as a function of x, is specified.Let dl and dl ′ denote the distance between two points very close to each other before and after the deformation respectively.We have and for small deformations [1] dl where is the strain tensor.We also have ϵ(u) = 1/2(∇u + (∇u) T ), where ∇u is the matrix field with rows the gradient of each component of u.
When a deformation occurs, the arrangement of molecules is changed and internal forces arise to return the body to its original state of equilibrium.They are encoded in the stress tensor σ.The ith component of the force acting on the element of surface ds with normal n is given by ∑ k σ ik n k .The conservation of angular momentum implies that the stress tensor is a symmetric matrix field [2], Page 53.
We now give the field equations of elasticity for an isotropic material, a material which has the same mechanical properties in all directions.From thermodynamic considerations, we have the following linear relationship, called Hooke's law, between the strain tensor and the stress tensor for small deformations Aσ = ϵ(u) (1) where Aσ = 1 2µ σ − λ 2µ(2µ+nλ) (tr σ)I and I is the n × n identity matrix.The constants λ and µ are called Lame coefficients.They satisfy 0 < µ 1 < µ < µ 2 and 0 < λ < ∞ for some constants µ 1 , µ 2 .Finally, the equilibrium condition is encoded in the equation For simplicity, we will assume homogeneous boundary conditions, u = 0 on the boundary ∂Ω of Ω.

Finite Elements
The finite element method is a widely used technique for approximating the solutions of boundary value problems arising in science and engineering.In its most general form, the starting point is a variational formulation of the equation.By analogy, a vector in R n is characterized by its coordinates or its inner product with vectors of the canonical basis.Here the solution u of the problem is sought in an infinite dimensional (closed) subspace U of a Hilbert space H and characterized as where F is a linear continuous functional on U and a is linear in the second argument.Under appropriate assumptions a solution exists and is unique.The finite element proceeds in constructing finite dimensional subspaces U h of U , where h is a parameter accumulating to zero associated to a decomposition of the domain of computation Ω into elements, known as triangulation.The discrete problem consists in solving for r h ∈ U h , the problem Again, under appropriate assumptions which will be stated below, the error between r and r h will be quantified.We will assume for simplicity that the computational domain for n = 2 and n = 3 has piecewise planar boundary so that it can be subdivided into triangles, for plane problems, and tetrahedra for three-dimensional problems.Furthermore, we will assume that the intersection of any two triangles is either empty or is a common edge or is a common vertex.Similarly, we require the intersection of two tetrahedra to be either empty or a common vertex, edge or face.We will also discuss rectangular finite elements, where the domain is assumed to be subdivided into rectangles or cubes, satisfying the analogous conditions.In all the above cases, h can be taken as the maximum of the diameter of the elements.Following [3], a finite element is given by a triple (K, U K , Θ K ) where 1.K is a closed bounded subset of R n with a nonempty interior and a Lipschitz continuous boundary, 2. U K is a finite dimensional space of vector-valued or matrix-valued functions defined over the set K, and 3. Θ K is a finite set of linearly independent linear functionals, θ i , i = 1, . . ., N referred to as degrees of freedom of the finite element, defined over the set U K .
It is assumed that the set Θ K is U K -unisolvent in the sense that The finite element space U h is defined as the space of piecewise polynomials which restrict to an element of U K on each element K.The degrees of freedom give the global continuity property of the finite element space.For example, if K is a triangle and we take as U K , the space of piecewise linear functions, a set of degrees of freedom can be defined as the values at the vertices.A piecewise polynomial on a triangulation with these degrees of freedom is necessarily continuous, since a linear polynomial on an edge is uniquely determined by the values at the vertices.Let {ϕ i , i = 1 . . ., N } ⊂ U K be a basis of U K , dual to Θ K , i.e. θ i (ϕ j ) = δ j i .The natural interpolation operator Π U h associated to U h is defined by

Notation
Let Ω be a simply connected domain in R n , n = 2, 3 with a shape regular family of meshes with mesh size h decreasing to zero.We will assume that the triangulation is regular in the sense that the aspect ratio [3], of the elements is bounded by a fixed constant.Thus no element can become too thin or too long.Let S denote the space of n × n symmetric matrix fields.We let L 2 (K, X) be the usual space of square integrable functions with domain K ⊂ R n and with values in the finite dimensional space X.For k ≥ 1 we denote by H k (K, X) the space of functions on K, taking values in X, and with all derivatives of order at most k in L 2 (K, X).For our purposes, X will be either S, R n , or R, and in the latter case, we will simply write Here the divergence of a matrix field is the vector field obtained after taking the divergence of each row.We let (, ) denote the L 2 (Ω, X) inner product.The space of infinitely differentiable functions on Ω with values in X is denoted by C ∞ (Ω, X).
We use the usual notations of P k (K, X) for the space of polynomials of degree at most k and P k 1 ,k 2 (K, X) for the space of polynomials on K of degree at most k 1 in x and of degree at most k 2 in y.We write P k and P k 1 ,k 2 respectively when X = R and it is clear that they are considered on K.We recall that the dimension of . For two elements σ and τ ∈ S, σ : τ = ∑ n i,j=1 σ ij τ ij is the Kronecker product of σ and τ .

Variational Formulations of the Elasticity Equations
There are several variational formulations on which one can base a finite element discretization of the elasticity equations.We describe briefly three of them following [4].In the primal variational formulation, the displacement is the unique critical point over The main drawback of this approach is that computing the stress from the displacement using (1) involves approximate derivatives and hence less accuracy.In the dual variational principle, the stress field can be computed directly in principle as the unique critical point of the complementary energy functional However enforcing constraints in the finite element method is a challenging task.On the other hand, ( 1) is an overdetermined system if one wants to compute the displacement from the stress.For example, for n = 2, two unknowns determine the displacement vector while (1), which is an equality between two 2 × 2 symmetric matrices, gives three equations.We now describe the mixed variational principle in which the stress field and the displacement form the unique critical point of the Hellinger-Reissner functional The mixed weak formulation is then given by: Find Alternatively, the above equations can be obtained directly from ( 1) and ( 2) by integration by parts ) and using the boundary condition for u.It is also possible to write them in the abstract form (3) for an unknown (σ, u).Clearly, the stress is approximated directly in the mixed formulation.But there are other reasons behind the popularity of mixed methods.In certain situations, e.g. for viscoelastic materials with an incompressibility constraint [5], it is not possible to eliminate the stress from the equations, hence the development of stable mixed elements for linear elasticity is a promising important step towards the resolution of related problems.On the other hand, mixed methods are in general robust in certain limiting situations [6][7][8], in the incompressibility limit for the elasticity equations.
It can be shown that the following conditions hold: where and These conditions in turn imply existence and uniqueness of a solution of the variational problem (4) [9].The condition ( 8) is known as inf-sup condition.The main disadvantages of mixed finite elements is that they lead to indefinite system of equations and their construction, especially for symmetric matrix fields, is challenging as Brezzi's stability conditions for mixed methods, discrete analogues of ( 6)-( 8), must be satisfied.

Mixed Finite Elements for Elasticity
A conforming mixed finite element approximation of (4) consists in choosing finite dimensional spaces We must therefore choose finite dimensional spaces for each of the n(n + 1)/2 components σ ij , i ≥ j of the stress field and each of the n components for the displacement u.To satisfy the condition V h ⊂ L 2 (Ω) n , we simply take spaces of piecewise discontinuous polynomials.However, the condition Σ h ⊂ H(div, Ω, S) is equivalent to the condition that σn is continuous across interior edges (for n = 2) or interior faces (for n = 3) of the triangulation [10].It is not easy to satisfy both requirements of the Brezzi conditions.Heuristically, the first one would be easier to satisfy if K h were small.That translates into Σ h small and V h big.For the second, one wants V h small since the condition must be satisfied by all v ∈ V h .And the bigger Σ h is, the bigger the supremum will be.So Σ h small and V h big for the first one and Σ h big but V h small for the second one.Another set of sufficient conditions which imply the Brezzi conditions [11], and are easier to satisfy are given by: • There exists a linear operator to h, and such that with commutes.We recall that a commutative diagram is a diagram with spaces and maps such that when picking two spaces, one can follow any path through the diagram and obtain the same result by composition.We then have the optimal error estimate with γ independent of h.The estimate (13) says that the errors in approximating the stress and the displacement depend only on the approximation property of the spaces Σ h and V h .The following approximation properties of finite element spaces [12], under the assumption of shape regularity, allows to have a rough estimate of the order of approximation under higher regularity assumptions.We have whenever Moreover the L 2 (Ω, R n ) projection operator onto a space of piecewise discontinuous polynomial of degree k satisfies Among the reasons to incorporate the symmetry of the stress field in finite element computations, it has been argued that it is necessary to obtain good free surface conditions [13].Perhaps also motivated by the need to preserve the structure of the continuous problem and to reduce the number of unknowns, researchers have unsuccessfully tried for several decades to construct stable conforming mixed finite elements for elasticity using polynomial shape functions.While for the corresponding problem for the mixed formulation of the Poisson equation several compatible pairs were known, e.g. the Raviart-Thomas and BDM elements [10], the symmetry of the stress field is a substantial difficulty in the case of the elasticity equations.Earlier researchers circumvent the difficulty of the symmetry condition by using composite elements, that is the stress and displacement are approximated on different grids [14][15][16][17].
The first conforming elements with symmetric stress fields were introduced in [18] on triangular meshes.In this section, we first review the Arnold-Winther construction, then discuss the various extensions which have appeared in the literature.

The Arnold-winther Construction
We recall that an exact sequence is a sequence of spaces and maps between them such that the image of one map is the kernel of the next.A key ingredient in the design is the use of differential complexes which are exact sequences where the maps are differential operators [18].In two dimensions, the elasticity complex is where J is the Airy stress operator defined by An analogous sequence with less smoothness is The above sequences are exact if the domain Ω is contractible [18].Intuitively, a domain is contractible if it can be continuously shrunk to a point.One way to look at the Arnold-Winther construction is that the diagram (12) can be extended to a commutative diagram where I h is a projection operator.The choice of a space Q h which makes the bottom sequence exact provides a discrete version of the elasticity sequence.The difficulty of choosing compatible pairs Σ h and V h is then shifted to the choice of a conforming finite element subspace of H 2 (Ω), a much less difficult task.
Recall that the objective is to choose finite dimensional spaces of polynomials on each element for each of the three components of the stress field and each of the two components of the displacement in such a way that σn is continuous across interior edges, div Σ h ⊂ V h and P h div τ = div Π h τ, τ ∈ H(div, Ω, S).We choose Q h as the Argyris finite element space [19] with Q K = P 5 (K) and the degrees of freedom given by the values of the derivatives up to order two at the vertices and the average of the normal component on each edge.Let V, E and R respectively denote the number of vertices , edges and triangles in the triangulation.We have dim Q h = 6V + E. We will write analogously dim Σ h = xV + yE + zR and dim V h = wR since V h is a space of piecewise discontinuous polynomials.We are then led to consider a discrete sequence which suggest possible choices Σ K = P 3 (K, S) and V K = P 2 (K, R 2 ).For the above sequence to be exact, the alternating sum of the dimensions must equal 0. On the other hand, by Euler's formula, We have A solution of the above equation is given by x = 3, y = 4 and z = w − 3. It now requires some guess to choose the degrees of freedom.Since there are 3 components for the stress field, it is natural to choose as degrees of freedom the value of each of them at the vertices.It is not difficult to see that a polynomial of degree k on an edge is uniquely determined by k + 1 independent conditions.Next, the requirement that σn is continuous across interior edges leads to the choice of degrees of freedom (1) the values of each component of τ (x) at the vertices of K (9 degrees of freedom), (2) the first two moments of each component of τ n on each edge (12 degrees of freedom).
But recall that for stability, we need div Σ h ⊂ V h and the commutativity property Sufficient conditions are v ∈ span{1, x} on each edge.We therefore choose V K = P 1 (K, R 2 ) which give w = 6 and z = 3.The lowest order Arnold-Winther construction is given by with the last set of degrees of freedom given by (3) the values of The space Q h also enters in the proof of unisolvency for the degrees of freedom of Σ h .If one assumes that all degrees of freedom of τ ∈ Σ K vanish, an important step in the proof of unisolvency is to show that div τ = 0. Next, the exactness of the sequence (17) give τ = Jq for q ∈ P 5 (K).One then has to be able to show that q = 0 using the degrees of freedom of τ .
We have described the lowest order Arnold-Winther element.A family of elements indexed by k ≥ 1 has been constructed with convergence rate O(h k+2 ) for the stress and O(h k+1 ) for the displacement.The interpolation operator Π h was chosen as a suitable modification of the natural interpolation operator of Σ h .These elements have the unusual feature of involving vertex degrees of freedom, compared to the BDM elements or Raviart-Thomas elements, and involve a fairly large number of degrees of freedom.There was a reduced element with 21 degress of freedom for the stress and the three-dimensional space of rigid body motions, {(a−cy, b+cx), a, b, c ∈ R} , to approximate the displacement.The convergence rate in that case is O(h 2 ) for the stress and O(h) for the displacement.With the Argyris interpolation operator I h , it can be shown that the diagram ( 16) commutes.We will illustrate the construction of higher order elements and reduced elements for the tetrahedral elements.

The Three-dimensional Tetrahedral Elements
In three dimensions, the elasticity sequence is given by Here T is the six dimensional space of rigid body motions, i.e., the space of linear polynomial functions of the form x → a+b×x for some a, b ∈ R 3 .One can require less smoothness in the following analogue of the above sequence, where curl curl * is a second order differential operator defined as follows: Given a symmetric matrix field, we first replace each column by its curl, and then replace each row of the resulting matrix with its curl.We recall that for a three-dimensional vector field, Analogous to the definition of H(div, Ω, S), we define A polynomial analogue of the above sequence is The above sequences are exact if the domain Ω is contractible [20,21].In an attempt to mimic the polynomial sequence at the discrete level, we face the difficulty that both Q h and Σ h are spaces of symmetric matrix fields.Instead we will draw in a first step analogies with the two-dimensional triangular elements.For integers s, k ≥ 0, we postulate that V K is a space of discontinuous piecewise polynomials of degree at most s, i.e.V K = P s (K, R 3 ), and Analogous to the triangular elements, the following degrees of freedom for the stress field will be chosen • vertex degrees of freedom, • degrees of freedom to ensure continuity of τ n, (18).
It can be shown that [21], when all the above degrees of freedom vanish for a given T ∈ Σ K , then T ∈ M k (K) := {τ ∈ Σ K , div τ = 0, τ n = 0 on ∂K}.The last set of degrees of freedom can then be chosen as ∫ K τ : ϕ, dxϕ ∈ M k (K).Since there are 4 vertices for a tetrahedron and 6 components for the stress field, we need 6 × 4 = 24 vertex degrees of freedom.Next, we specify the edge degrees of freedom to ensure continuity of τ n across edges.For each edge with unit normals n − and n + , and unit tangential vector t, , can be specified independently.This requires 5(k − 1) degrees of freedom.For a face T of the triangulation, we recall that dim P k (T, R) = (k+2)(k+1) 2 .For a given face, we have specified all 3 vertex and 3(k − 1) edge degrees of freedom.Hence it remains to specify (k+2)(k+1) 2 R) degrees of freedom on each face to ensure continuity of τ n.The face degrees of freedom are then chosen as We now describe the lowest order element of the family of tetrahedral elements with k = 4 and s = 1.We take We have dim Σ K = 162 and dim V K = 12, with the stress degrees of freedom (d.o.f.)That M 4 (K) has the required dimension 6 was proved in [21].An analogous element was obtained in [22], where the last set of degrees of freedom was replaced by ∫ e T t • t dx for each of the six edges.The other elements of the family can be described as The degrees of freedom are defined in a fashion similar to the lowest order case [21].We obtained O(h k+2 ) convergence rate for the stress and O(h k+1 ) convergence rate for the displacement.It should be mentioned that the main difficulty was to prove that M 4 (K) has dimension 6 and construct bases for M k (K) for k ≥ 4, [21].Moreover, analogous to the two-dimensional elements, it can be shown that the following diagram, with suitably defined projection operators, commutes [23], It is possible to somewhat reduce the number of degrees of freedom for the lowest order element by taking as V K the six dimensional space of rigid body motions [21].The stress space then has local dimension 156.

Rectangular Elements
With Arnold, we extended the triangular elements to rectangular meshes [24].The polynomial differential complexes are based on the P k,l spaces.Our lowest order element has 36 degrees of freedom for the stress and 3 degrees of freedom for the displacement with O(h) convergence rate for both the stress and displacement.It was recently observed in [25], that on rectangular meshes it is not necessary to impose vertex degrees of freedom for all components of the stress field.This lead to a different family with the lowest order element having 17 degrees of freedom for the stress field and 4 degrees of freedom for the displacement.A low order conforming rectangular elements in three dimensions has recently been constructed in [26] with 72 degrees of freedom for the stress and 12 degrees of freedom for the displacement.

Nonconforming Elements
The high number of degrees of freedom of the conforming elements has prompted the search for more practical elements.One possibility to achieve this is to relax the conformity condition.In this section, we consider nonconforming mixed elements for the elasticity problem in the sense that the stress space Σ h is not contained in H(div, Ω, S) because the normal component of the discrete stress field is not continuous across element edges.The displacement space V h satisfies V h ⊂ L 2 (Ω, R n ).In addition to their relative simplicity, the nonconforming elements do not involve any vertex degrees of freedom.
The framework described above extends to the nonconforming cases except that a so-called consistency error [27] has to be shown to be at least O(h).In the two-dimensional case, it takes the form where E h is the set of edges of the triangulation and [q] denotes the jump of the quantity q across the given edge.Arnold and Winther, [27] obtained a 12 d.o.f.element for the stress with 3 d.o.f. for the displacement on triangular meshes.The analogous nonconforming rectangular elements [28], uses 16 d.o.f. for the stress and 3 d.o.f. for the displacement.Other nonconforming rectangular elements have also been constructed.In [29], a 13 d.o.f. for the stress and 4 d.o.f. for the displacement was constructed.A similar element with the same number of degrees of freedom was given in [30].These elements exhibit O(h) convergence rate for both the stress and displacement and can be written as part of a discrete version of the elasticity sequence.In [31], an element with 19 d.o.f. for the stress and 6 d.o.f. for the displacement was introduced.The construction does not use the commutative diagram, but instead the Brezzi conditions are proved directly.These elements also extend to three-dimensional rectangular elements with 60 d.o.f. for the stress and 12 d.o.f. for the displacement.They have O(h 2 ) convergence rate for the displacement and O(h) convergence rate for the stress.Another three-dimensional element on rectangular meshes was introduced in [32].It uses 54 d.o.f. for the stress and 12 d.o.f. for the displacement with O(h) convergence rate for both the stress and displacement.We recently constructed three-dimensional nonconforming tetrahedral elements with 36 d.o.f for the stress and 6 d.o.f. for the displacement with linear convergence rates in both variables [33].None of the three-dimensional nonconforming rectangular elements have been shown to be part of a discrete version of the elasticity sequence, that is an analogue of (16).

Conclusions
Other researchers have abandoned the symmetry condition, used a different variational formulation, or weakened the regularity conditions on the stress fields in various ways [34][35][36][37][38][39].Stable mixed finite elements with weakly imposed symmetry have been introduced in [17,[40][41][42][43][44], but the simplest ones have only been recently discovered [13,[45][46][47].The construction of discrete versions of the elasticity complex has been made further systematic in [20,45,46] where a procedure to construct discrete versions of the elasticity complex from discrete versions of the de Rham complex, an exact sequence connecting spaces of differential forms, was made explicit.This process has been further improved in the finite element exterior calculus [20,48].Finally, we note that despite their relative complexity, conforming mixed finite elements with symmetric stress field are useful in certain situations [49].