1 Introduction
It is wellknown that incompressible nonlinear elasticity is useful to describe the mechanical behavior of various soft tissues [1]. Deriving stable numerical methods for modeling such tissues that may undergo large deformations is a challenging task, for example see [2] and references therein. In this paper, we develop fourfield mixed finite element methods for incompressible nonlinear elasticity by extending threefield mixed finite element methods introduced in [3] for compressible nonlinear elasticity. These fourfield methods are based on a mixed formulation with displacement, displacement gradient, the first PiolaKirchhoff stress tensor, and pressure as the independent unknowns. A potential advantage of this mixed formulation is providing accurate approximations of strain and stress, since they are independent unknowns and are explicitly calculated in the solution process, that is, no postprocessing is required to indirectly compute strain and stress from displacement.
Similar to [3], finite elements suitable for the and the operators are respectively employed to discretize displacement gradient and stress in the proposed conformal mixed finite element methods. One can show that these choices follow from the strain compatibility and the balance law [4, 5]. To discretize displacement and pressure,  and conformal finite element spaces are used, respectively.
A stability analysis based on the abstract theory of [6, 7] for Galerkin approximations of nonlinear problems is provided. More specifically, an infsup condition is written which is locally sufficient for the convergence of solutions of nonlinear finite element methods to a regular solution of incompressible nonlinear elasticity. This infsup condition can be interpreted as the necessary and sufficient condition for the wellposedness of the linear problems associated to Newtons’ iterations for solving nonlinear finite element methods. Based on this interpretation, we also write other infsup conditions that are necessary for the stability of Newtons’ iterations. Using these conditions, we derive some relations between the dimension of finite element spaces for different unknowns, which are necessary for the stability of Newtons’ iterations.
By considering choices of simplicial finite elements of degree less than or equal to for obtaining mixed finite element methods, we conclude that choices in D and choices in D do not satisfy all these infsup conditions, in general, and lead to unstable methods. The performance of stable choices in D and D are studied by solving numerical examples.
This paper is organized as follows: The fourfield mixed formulation and the associated conformal finite element methods are discussed in Sections 2 and 3. Newtons’ iterations for solving nonlinear finite element methods are provided in Section 4. In Section 5, different infsup conditions are derived for studying the stability of finite element methods and it is shown that some choices of finite elements lead to unstable methods as they violate the infsup conditions. In Section 6, the infsup conditions are numerically studied. Moreover, by solving two numerical examples, the performance of stable finite element methods are investigated in D and D. Finally, some concluding remarks are given in Section 7.
2 A Mixed Formulation for Incompressible Elasticity
We write a fourfield mixed formulation for incompressible nonlinear elasticity in this section. For simplicity, we consider static problems here, however, this formulation can be readily extended to timedependent problems as well. Let denote the reference configuration of a D or D elastic body and let
be the unit outward normal vector field associated to the boundary
of .Suppose are respectively displacement, displacement gradient, and the first PiolaKirchhoff stress tensor. The constitutive relation of incompressible elastic bodies can be expressed as , where is a given function of , the realvalued function is pressure, and is the identity tensor [8, Section ]. The boundary value problem of incompressible nonlinear elastostatics can be stated as:
Given a body force , a displacement of , and a traction vector field on , find such that
To write a week form based on the above equations, we consider the following spaces: Let be the standard space of squareintegrable functions and let be the Sobolev space of functions with firstorder derivatives of class. The space of vector fields in , , with components of class is denoted by . The space of vector fields that vanish on is written as .
We also consider the spaces and , which are the spaces of secondorder tensor fields with components that their and are respectively of class. Consider the matrix representing the Cartesian components of a secondorder tensor field in . Since and are respectively obtained by applying the standard and of vector fields to each row of , the spaces and can be considered as copies of the corresponding spaces for vector fields.
For obtaining a concise form for the weak formulation, we employ the following notation: The standard inner product of is denoted by “” and the inner products of realvalued functions, vector fields, and tensor fields are denoted by , that is, , , and , where we consider the summation convention on repeated indices. By taking the inner product of the equations (2.1) with suitable test functions and using Green’s formula, one can write the following fourfield mixed formulation for incompressible nonlinear elasticity:
Given a body force , a boundary displacement , and a surface traction vector field on , find such that , on , and
The above weak form is an extension of the formulation of [3] for compressible elasticity to incompressible elasticity. The choices of the spaces and guarantee that physically relevant jump conditions, namely, the Hadamard jump condition for strain and the continuity of traction vector fields at interfaces for stress, will be hold in the weak sense and also on the discretelevel. Following the discussion of [3, Remark 2], it is not hard to show that even for hyperelastic materials, the above weak formulation is not associated to a stationary point of any functional.
3 Conformal Finite Element Methods
Suppose , , , and are finite element spaces such that , , , and . Also let and let
be the interpolation operator associated to
. Then, one can write the following conformal mixed finite element method for (2.3):Given a body force , a boundary displacement , and a surface traction vector field on , find such that , on , and
We employ simplicial finite elements shown in Figure 1 to generate the above finite element spaces in D and D. More specifically, we use the standard simplicial Lagrange elements and discontinuous elements respectively for and . As mentioned earlier, the definitions of and for secondorder tensors imply that in , the tensorial spaces and are equivalent to copies of the associated vectorial spaces in the sense that each row of secondorder tensors of the  and classes can be considered as a vector field beloning to the vectorial spaces and corresponding to the standard and operators of vector fields, respectively. We use the Nédélec elements of the first and the second kinds to generate and the RaviartThomas elements and the BrezziDouglasMarini elements to generate .
4 Newtons’ Iterations
We employ Newton’s method to solve (3.1). Here we mention this approach in some details since, as will be discussed in the following section, it is also useful for stability analysis. Let and . Also let and , with , , and . The nonlinear problem (2.3) can be written as: Find such that
(4.1)  
The th Newton’s iteration for solving the above nonlinear equation reads: Given , find by solving a linear problem and let . This linear problem is obtained by linearizing the above equation and can be stated as: Find such that
(4.2) 
where the bilinear form is given by
(4.3)  
where is the elasticity tensor in terms of the displacement gradient and . By discretizing the linear problem (4.2) using the finite element spaces of the previous section, one obtains the following Newton’s iteration for the finite element method (3.1):
Let , and . Given , find and let , where is obtained by solving the linear finite element method
(4.4)  
Suppose
(4.5) 
and let , denote the total number of degrees of freedom. Then, it is straightforward to see that the stiffness matrix of the above linear problem is of the form
(4.6) 
5 Stability Analysis
The stability of the nonlinear finite element method (3.1) can be studied by using the general theory for the Galerkin approximation of nonlinear problems discussed in [6, 7]. Let be a regular solution of (2.3), or equivalently (4.1), in the sense that the derivative of the nonlinear mapping defined in (4.1) is nonsingular at . This derivative can be expressed in terms of the bilinear form introduced in (4.3). The theory of [6, 7] states that in a neighborhood of and for sufficiently refined meshes with the maximum element diameter , under some additional mild conditions, the nonlinear finite element method (3.1) has a unique solution that converges to as if there exists a meshindependent number such that
(5.1) 
where is the mixed finite element space for the Newton’s iteration (4.4) and is the norm of . Thus, the infsup condition (5.1) is a sufficient stability condition.
The condition (5.1) can be interpreted as a stability condition for the linear finite element method of the Newton’s iteration (4.4) as well. More specifically, suppose that at the th Newton’s iteration, we have . Then, the condition (5.1) implies that the linear problem (4.4) at the the th iteration has a unique solution and the stiffness matrix introduced in (4.6) is nonsingular. One can also approximate the lower bound of (5.1) by using a matrix associated to [9, Proposition 3.4.5]. In particular,
can be approximated by the smallest singular value of the matrix
, where is the symmetric, positive definite matrix associated to the norm .The nonsingularity of the stiffness matrix implies that the submatrices
(5.2)  
(5.3) 
must be full rank. This result can be stated by using some infsup conditions associated to suitable bilinear forms induced by (4.4). The upshot can be stated as follows:

is full rank if and only if there exists such that
(5.4) 
is full rank if and only if there exists such that
(5.5) 
is full rank if and only if there exists such that
(5.6) where
(5.7) 
is full rank if and only if there exists such that
(5.8) where
(5.9) (5.10) 
is full rank if and only if there exists such that
(5.11) where
(5.12) (5.13) 
is full rank if and only if there exists such that
(5.14) where
(5.15) (5.16) (5.17)
A simple inspection of the number of rows and columns of the above matrices together with the ranknullity theorem yield the following necessary condition for the validity of (5.1) and the wellposedness of (4.4), where we use the notation of (4.5).
Theorem 1.
By using this result, one can specify some unstable combinations of finite elements. For example, consider the D elements of Figure 1 and suppose , , and are respectively the number of vertices, edges, and elements of a D simplicial mesh. One can write the relations
(5.18) 
where is the number of holes and is the number boundary edges [10]. By using the notation of Figure 1, we obtain the following relations on a D simplicial mesh:
(5.19)  
(5.20) 
A simple varification of the inequalities of Theorem 1 by using the relations (5.18) and (5.19) yields the following result.
6 Numerical Results
Some numerical examples are discussed in this section to study the stability and the convergence of the mixed finite element methods (3.1). We consider the simplicial finite elements of Figure 1 and use FEniCS [11] to implement these numerical examples. The underlying domains are assumed to be a unit square in D and a unit cube in D with simplicial meshes shown in Figure 2. We consider incompressible NeoHookean materials with the stored energy function , , and the constitutive equation
(6.1) 
The term of the bilinear form (4.3) then simply reads .
Notation.
To refer to any choice of the elements of Figure 1 for the discretization of (2.3), we use the abbreviated names shown in Figure 1. For example, indicates the choice of the Lagrange element of degree , the degree Nédélec element of the first kind, the RaviartThomas element of degree , and the discontinuous element of degree respectively for the discretization of .
6.1 Stability Study
The infsup conditions introduced earlier hold only if the corresponding matrices are full rank. Therefore, a simple approach for computationally studying these infsup conditions is to compute the rank deficiency of the matrices , , , , , . Some choices of finite elements that violate these infsup conditions are shown in Figure 3. In this figure, the fullrankness of the above matrices are plotted for several choices of finite elements by using D and D meshes of Figure 2, where the fullrankness of a matrix is the rank of divided by its maximum possible rank. Thus, , if is full rank and , otherwise.
The results of Figure 3 are computed using the incompressible NeoHookean constitutive equation (6.1) with near the reference configuration, that is, the matrices are associated to the first Newton’s iteration (4.4) starting at the reference configuration. The rank deficiency of some finite elements such as for , for , for , and for , is a simple consequence of the number of columns of the associated submatrix being smaller than the number of its rows as mentioned in Theorem 1. For some finite elements such as for , the associated submatrix is nearly square and as shown in Figure 3, it is also nearly full rank. In this situation, the socalled locking phenomena may occur in the sense that the associated submatrix may represent an injective operator, and thus, it only admits a unique solution.
Our computations suggest that out of possible choices of simplicial elements of Figure 1 for displacement, displacement gradient, stress, and pressure, choices in D and choices in D satisfy all infsup conditions. More specifically, the choices that violate each infsup condition are as follow:

In D and D, the infsup condition (5.4) does not hold for the choices of elements and for displacement and stress.

The infsup condition (5.5) does not hold with the following choices of elements for displacement gradient and pressure:

2D: and ;

3D: and , .


The infsup condition (5.6) does not hold with the following choices of elements for displacement, displacement gradient, and stress:

2D: , , , and , ;

3D: and , .


In D and D, the infsup condition (5.8) does not hold with the following choices of elements for displacement, displacement gradient, stress, and pressure: , , , and , , .

The infsup condition (5.11) does not hold with the following choices of elements for displacement, displacement gradient, and pressure:

2D: and , ;

3D: and , , .


The infsup condition (5.14) does not hold with the following choices of elements for displacement gradient, stress, and pressure:

2D: and , ;

3D: , , , and , .

The infsup conditions of Section 5 are not independent. For example, the infsup conditions (5.4) and (5.6) are closely related. Such dependencies can be read off from the above numerical results. Common strategies may be employed to stabilize choices that violate the above infsup conditions such as using different underlying meshes for different unknowns, enriching trial and test spaces using bubble functions, and using perturbed formulations [9].