Weakly symmetric stress equilibration for hyperelastic materialmodels

A stress equilibration procedure for hyperelastic material models is proposed andanalyzed in this paper. Based on the displacement-pressure approximation computed with a stable finite element pair, it constructs, in a vertex-patch-wise manner, an $H(div)$-conforming approximation to the first Piola-Kirchhoff stress. This is done in such a way that its associated Cauchy stress is weakly symmetric in the sense that its anti-symmetric part is zero tested against continuous piecewise linear functions. Our main result is the identification of the subspace of test functions perpendicular to the range of the local equilibration system on each patch which turn out to be rigid body modes associated with the current configuration. Momentum balance properties are investigated analytically and numerically and the resulting stress reconstruction is shown to provide improved results for surface traction forces by computational experiments.


INTRODUCTION
This paper is concerned with a stress equilibration procedure for hyperelastic material models in nonlinear solid mechanics. It extends the approach proposed and studied in our earlier work [1] to the case of geometrically and materially nonlinear elasticity in the form of a hyperelastic material law. Due to the fact that the symmetry condition does not hold for the first Piola-Kirchhoff stress (which is the result from the reconstruction process) but for the Cauchy stress, the use of symmetric stress elements is not feasible anymore in the hyperelastic case. The weak symmetry condition from linear elasticity can, however, be generalized to a suitable constraint for the Piola-Kirchhoff stress as is done in this contribution. To the best of our knowledge, our contribution is the first attempt to develop a stress equilibration procedure for the hyperelastic situation. Our hope is that this will be of use for the development of an a posteriori error estimator for hyperelastic problems in the future. The issue of a posteriori error estimation and adaptive refinement is, however, beyond the scope of this contribution.
Expressing the internal forces of a material, the components of the stress-tensor are crucial for the prediction of the weakening of a material, including plastic behavior or damage. A specific application area where this is an issue is associated with implant shape design which constitutes an optimal control problem. [2] Therefore, the accurate approximation of the stress-tensor is of strong importance in numerous applications and in particular in the hyperelastic material model this paper is concerned with. The mathematical foundations of hyperelastic material models in solid mechanics are covered, for example, in the books by Marsden and Hughes [3] and Ciarlet. [4] The numerical treatment of the associated variational problems is investigated in detail by Le Tallec. [5] Specifically for incompressible hyperelasticity, issues connected to the use of displacement-pressure formulations are discussed. [6] A priori analysis of numerical methods is available under restrictive assumptions, see Carstensen and Dolzmann [7] and, for a least-squares finite element approach, Müller et al. [8] Common displacement-based approaches or, in the incompressible regime, mixed displacement-pressure formulations for this model lead to approximations of the stresses that are not H(div)-conforming, that is, have discontinuities of the normal components on the interface between two elements. In particular, this means on the one hand that they do not control momentum conservation and on the other hand that the normal component of the boundary traces is not well defined implying that the approximation of the surface traction forces can also not be guaranteed. In contrast to variational principles involving a direct approximation of the stress in an H(div)-conforming space (see chapter 9 of the monograph [9] for an overview), this paper proposes an algorithm to obtain an H(div)-conforming approximation of the stress-tensor by postprocessing the displacement-based approximation.
The idea of reconstructing the matrix-valued stress and vector-valued flux goes back to the hypercircle theorem by Prager and Synge [10] (see also section III.9 in Braess' book [11] for a presentation in modern mathematical language). Besides the accurate approximation in an H(div)-conforming space, the stress or flux reconstruction builds the basis of an a posteriori error estimator, which was actually already one of the motivations of Prager and Synge. [10] Over the years, a posteriori error estimators based on flux reconstruction were explored in detail in many contributions. [12][13][14][15][16] An important algorithmic innovation was given by Braess and Schöberl [17] by the equilibration procedure which is completely local and provides the link to residual error estimation. An important aspect of the use of reconstruction-based error estimation of the above type is that it provides guaranteed upper bounds for the error with accessible constants. Another important aspect is that these a posteriori error estimators are valid for any approximation that is inserted into the procedure. In particular, it does not assume that the underlying finite-dimensional variational problems are solved to high precision. The extension of reconstruction strategies to linear elasticity was the subject of a number of contributions in the last two decades, [18][19][20][21][22] stress reconstruction in the context of Stokes flow was also studied recently. [23] More recently, a posteriori error estimation based on the reconstruction of weakly symmetric stresses was investigated in our earlier work. [1,24] In particular, the stress equilibration procedure considered in our recent contribution [1] serves as a point of departure for our treatment of hyperelastic material models in the present paper. The recent paper by Botti and Riedlbeck [25] should also be mentioned here. It treats nonlinear elasticity restricted to a geometrically linear situation. In that case, the (Piola-Kirchhoff) stress is still symmetric which allows the use of symmetric stress elements as it is done in the approach by Botti and Riedlbeck. [25] Besides the fact that the symmetry condition for the stresses becomes more complicated in the geometrically and materially nonlinear situation associated with hyperelastic models which was already mentioned above, other challenging issues arise if one wants to extend the stress equilibration procedure from our recent work [1] to that case. The stresses computed directly from displacement and, possibly, pressure approximations are no longer piecewise polynomial due to the nonlinearity of the model. Therefore, in order to get a stress reconstruction in an appropriate H(div)-conforming finite element space, a suitable projection to piecewise polynomial stresses needs to be carried out first. Another problem is concerned with the subspace of test functions which are perpendicular to the range of the local equilibration systems for vertex patches not connected to the Dirichlet boundary. The main result of this contribution is the identification of these subspaces as associated with rigid body modes in the current configuration, that is, involving the displacement approximations. The right-hand sides arising from straightforward piecewise polynomial projections of the stresses are shown to have components outside of these ranges which means that the local equilibration systems possess an additional compatibility error. We propose a remedy involving a more complicated test space to overcome this problem. This leads to compatible local problems and thus to a truly equilibrated stress reconstruction. The development of an a posteriori error estimator based on the stress equilibration for hyperelastic material models and, in particular, its analysis are expected to be rather involved and to require rather restrictive assumptions. After all, it is well known that the solution of the variational problem may not be unique (see the examples in chapter 5 by Ciarlet [4] ). Other approaches to the direct finite element approximation of stresses in geometrically nonlinear elasticity can be found. [8,26] The outline of this paper is as follows. We start with the variational formulation of elastic deformations governed by hyperelastic material models and the weakly symmetric stress reconstruction in Section 2. Section 3 presents the local equilibration algorithm. The solvability of the local problems on vertex patches is analyzed in Section 4. In particular, the subspace of test functions orthogonal to the range of the local operators associated with equilibration is identified and this result is used for the investigation of the compatibility of the right-hand side. Section 5 proposes our remedy to deal with this problem and derives a more complicated test space which leads to compatible local equilibration systems for which an inf-sup condition holds. The improved accuracy of the surface forces associated with the equilibrated stresses will be the topic of Section 6. Finally, computational results illustrating the properties of the equilibrated stresses are collected in Section 7.

HYPERELASTICITY AND WEAKLY SYMMETRIC STRESS RECONSTRUCTION
The hyperelastic problems under our consideration are based on an open, bounded, and connected domain Ω ⊂ R d (d = 2, 3) with Lipschitz-continuous boundary which constitutes the reference configuration of the undeformed state. The boundary is divided into two disjoint nonempty subsets Γ D and Γ N . On Γ D , homogeneous displacement boundary conditions u = 0 are imposed, while surface traction forces P ⋅ n = g are prescribed on Γ N . For an appropriate subspace V ⊂ 1 Γ (Ω) , the boundary value problem of hyperelasticity then consists in the variational problem of finding u ∈ V such that holds for all v ∈ V. Here, P(u) = F (B) denotes the first Piola-Kirchhoff stress tensor with respect to the stored energy function ∶ R × sym → R, where the deformation gradient is given by F(u) = I + ∇u and the left Cauchy-Green strain tensor is defined as B(u) = F(u)F(u) T . Simple brackets (⋅, ⋅) as in (1) will from now on always abbreviate the inner product in L 2 (Ω) with respect to the reference configuration; f and g stand for volume and surface loads, transformed back to the reference configuration. An example of a stored energy function which we will also use later in our computations in Section 7 is associated with the Neo-Hookean model ) . ( In this case, the Piola-Kirchhoff stress tensor is given by and V = 1,4 Γ (Ω) would be sufficient for the variational problem (1) to be properly defined. In order to deal with materials in the incompressible parameter regime ( ≫ ), a pressure-like variable may be introduced, for example, by setting p = (det(F(u)) − 1). Note that with this choice, p does not really stand for the physical pressure but that it is possible to obtain the pressure from p even in the incompressible limit. The above choice is motivated from the fact that it turns into the constraint p = div u, familiar from linear elasticity, in the small strain limit. Other options for the definition of p are possible and may have advantages. The Piola-Kirchhoff stress is now given in terms of u and p which, in the Neo-Hookean example, reads due to the fact that det(B(u)) − 1 = det(F(u)) 2 − 1 = (det(F(u)) − 1)(det(F(u)) + 1) holds. With a pressure space Q (Q = L 4/3 (Ω) would be appropriate in the Neo-Hooke case), the variational problem turns into one of saddle point type which consists in finding u ∈ V and p ∈ Q such that with Q ′ denoting the dual space of Q (Q ′ = L 4 (Ω) with the above choices for the Neo-Hooke case). For k ≥ 1, let V h ⊂ V be the subspace of continuous piecewise polynomials of degree k + 1 with respect to a triangulation  ℎ for each component of V h . Our finite-dimensional variational problem for hyperelasticity consists in finding u h ∈ V h such that holds for all v h ∈ V h . In the incompressible regime, a discrete pressure space Q h consisting of continuous piecewise polynomials of degree k may be used to define a corresponding discrete saddle point problem. It consists in finding is satisfied. The direct use of P(u h ) or, in the incompressible regime, P(u h , p h ) as an approximation for the Piola-Kirchhoff stress, has, however, certain deficiencies which are already known from the linear elasticity situation. Most importantly, P(u h ) ⋅ n is not continuous at interfaces between elements of the underlying triangulation implying that traction forces are not well-defined. It also means that P(u h ) is not H(div)-conforming and that the conservation of momentum is not controlled. This motivates the need to construct an H(div)-conforming stress reconstruction P ℎ with all these desired properties.
The idea of equilibration is to compute the reconstructed stress P ℎ in the H(div)-conforming Raviart-Thomas space of degree k as an additive correction to P(u h ). This is done using the broken Raviart-Thomas space of degree k for each row leading to where P k (T) denotes the space of polynomials of degree k on the triangle (d = 2) or tetrahedron (d = 3) T. In other words, each row of the stress tensor P ℎ ∈ Δ ℎ is element-wise given by a function in the Raviart-Thomas space. Unfortunately, in contrast to the linear elasticity situation, P(u h ) ∈ Δ does not hold, in general, due to the nonlinearity of the stress-strain relation. Obviously, for the Neo-Hookean model in (3), P(u h ) is not even piecewise polynomial. Therefore, P(u h ) needs to be projected first to an elementP ℎ (u ℎ ) ∈ Δ ℎ . An obvious candidate would be to setP ℎ (u ℎ ) =  ℎ P(u ℎ ), where  ℎ denotes the component-wise and element-wise L 2 -orthogonal projection onto P k (T). We will stick with this choice ofP ℎ (u ℎ ) for the moment until we present an alternative one in Section 5 as a remedy for certain deficiencies associated with it.
Following the weakly symmetric equilibration procedure from Bertrand et al., [1] we perform the construction for the difference P Δ ℎ ≔ P ℎ −P ℎ (u ℎ ) between the reconstructed and the projected original stress. Recall that the extension of the hypercircle theorem to linear elasticity requires a symmetric reconstruction satisfying the equilibration condition div P Δ ℎ = −f − divP ℎ (u ℎ ) in each triangle and the jump condition allowing P ℎ to be H(div)-conforming. In order to write this jump condition in a precise way, let  ℎ denote the set of all sides (edges in 2D and faces in 3D) of the triangulation  ℎ and  * ℎ the set of sides not contained in Further, for all sides ∈  ℎ , let n be the normal direction associated with S (depending on its orientation), T + and T − the elements adjacent to S (such that n points into T + ) and the jump of P h over S defined by For sides S ⊂ Γ N located on the Neumann boundary we assume that n points outside of Ω and define the jump by In order to use the same formulas also for patches adjacent to the Neumann boundary Γ N we define the auxiliary jump by With this, the jump condition for the correction reads Similarly as in Bertrand et al., [1] the symmetry condition will be imposed weakly in order to obtain a reconstructed stress with reasonable symmetry properties. In the hyperelastic setting, symmetry does not hold for P(u) but instead for the related Cauchy stress tensor (u) = P(u)F(u) T /det(F(u)) which adequately describes stresses in the deformed configuration. Rewritting the equilibration and jump conditions in a weak form and applying the weak symmetry condition to P ℎ F(u) leads to the following conditions for P Δ ℎ : Z h may be chosen to be the space of discontinuous d-dimensional vector functions which are piecewise polynomial of degree k, and X h may stand for the continuous d(d − 1)/2-dimensional vector functions which are piecewise polynomial of degree k with J( ) being defined by for every d(d − 1)/2-dimensional vector . This choice is motivated by the inf-sup stability of the corresponding combination with the use of Raviart-Thomas element of degree k ≥ 1 as stress approximation space in the Hellinger-Reissner formulation. [27]

LOCAL STRESS EQUILIBRATION ALGORITHM
For the sake of the efficient computation of the stress reconstruction, we localize the problem using a partition of unity. The commonly used partition of unity with respect to the set  ℎ of all vertices of  ℎ , consists of continuous piecewise linear functions̃. In this case, the support of̃is restricted tõ In analogy to the stress equilibration procedure described by Bertrand et al. [1] for the linear elasticity case, we modify this classical partition of unity in order to exclude patches formed by vertices z ∈ Γ N , where the local problems may possess to few degrees of freedom to be solvable. To this end, let  ′ ℎ = { ∈  ℎ ∶ ∉ Γ } denote the subset of vertices which are not located on a side (edge/face) of Γ N . The modified partition of unity is defined by For ∈  ′ ℎ not connected by an edge to Γ N the function z is equal tõ. Otherwise, the function z has to be modified in order to account for unity at the connected vertices on Γ N . For each z N ∈ Γ N one vertex z I ∉ Γ N connected by an edge with z N is chosen and̃is extended by the value 1 along the edge from z I to z N to obtain the modified function . The support of z is denoted by ≔ ∪{ ∈  ℎ ∶ = 1 for at least one vertex of }.
For the partition of unity (14) to hold, we require the triangulation  ℎ to be such that each vertex on Γ N is connected to an interior edge. For the localized equilibration algorithm, we will also need the local subspaces for all ∈  ′ ℎ . Moreover, we need to work with the local sets of sides  ℎ, ≔ { ∈  ℎ ∶ ⊂ } and the restrictions Z h,z and X h,z to z of the test spaces Z h and X h , respectively. The conditions in (10) can be restated for a sum of patch-wise contributions leading, for each ∈  ′ ℎ , to the following minimization problem: subject to the constraints , ⟩ for all ∈ ( ) , ∈  ℎ, , Similarly to the linear elasticity case treated in our earlier work, [1] solutions to (18) are not expected to be unique. This is due to the fact that the number of linearly independent constraints in (18) is less than the dimension of the space Δ ℎ, , in general. At this point, we may introduce the local orthogonal projections  ℎ, ∶ 2 ( ) → Z ℎ, and  ℎ, ∶ 2 ( ) → ( ) which means that the first two conditions in (18) can be written shortly as ).
For each ∈  ′ ℎ , (18) constitutes a low-dimensional quadratic minimization problem with linear constraints for which standard methods are available for the efficient solution. Note that it is not guaranteed at this point that (18) has a solution at all. In fact, it does not, in general, as will become clear from the results of the next section. This is the reason why we will modify the test space in Section 5 in order to have well-posed local patch problems.
To get an idea about the structure of the system (18) and as a motivation for the result in the next section, we consider its underlying continuous problem. On the continuous level, the system (18) constitutes the stress-based dual formulation of the variational problem (1) restricted to z . With a suitable subspace V ⊂ 1 Γ ∩ ( ) this means that z ∈ V z is sought such that holds. On vertex patches with Obviously, all constants are contained in V • . Moreover, since (P(z), ∇v) = (P(z)F(z) , ∇vF(z) −1 ) holds and since P(z)F(z) T is a symmetric matrix, also all v with ∇vF(z) −1 being skew-symmetric will be contained in V • . In two dimensions, we arrive at ) , being contained in V • , and for d = 3 this is true for These are exactly the rigid body modes associated with the current configuration deformed by (x) = x + z which we would like to denote by RM(z) from now on. From the above derivation, it should not be surprising that the corresponding rigid body mode spaces RM(u h ) will appear in the investigation of the well-posedness of the discrete local problems (18) in the following section.

SOLVABILITY OF THE LOCAL PROBLEMS ON VERTEX PATCHES
We turn our attention to the solvability of the local minimization problem ‖P Δ ℎ, ‖ 2 → min! subject to the constraints (18). To this end, we need to guarantee that for every right hand side, a function P Δ ℎ, ∈ Δ ℎ, exists such that the constraints (18) are satisfied. The left-hand side in (18) defines a linear operator  ℎ, ∶ Π Δ ℎ, → Z ′ ℎ, × S ′ ℎ, × X ′ ℎ, , where S ℎ, = { ∈ ( ) ∶ ∈  ℎ, } denotes the trace space on the interior sides and ( ⋅ ) ′ stands for the dual space. The subspace R ⊥ ℎ, ⊆ Z h, z × S h,z × X h, z orthogonal to the range of  h,z , that is, the null space of its adjoint  * ℎ, , is obviously of interest for the solvability since the linear functionals on the right-hand side in (18) need to vanish on R ⊥ ℎ, . To anticipate the results of this section we will show and justify the nonsolvability of (18). This is done by further investigation of the subspace R ⊥ ℎ, which can be characterized as follows. Proposition 1. The subspace that is, the null space of the adjoint operator  * ℎ, associated with the constraints (18) can be characterized as follows: Here, | ⋅ | denotes the d − 1-dimensional measure of boundary curves or surfaces, respectively.

Proof.
The proof is carried out for d = 3; the two-dimensional case is much easier and can be derived from the three-dimensional one in the usual way by setting u 3 ≡ 0 and all other functions to be independent of x 3 (with appropriate modifications of operators such as div, , curl, etc.). 1st Step. We start by showing that the component h,z of R ⊥ ℎ, in (22) needs to satisfy Let us restrict ourselves to the H(div)-conforming subspace of Δ ℎ, , that is, with the property that [[P Δ ℎ, ⋅ n]] = 0 for all ∈  ℎ, . Then, the condition in (22) for the definition of R ⊥ ℎ, turns into By definition, we can write We may restrict ourselves further to divergence-free P Δ ℎ, with P h,z ⋅n = 0 on the entire boundary z . These stress approximations can be written as P Δ ℎ, = curl ℎ, with h,z in the Nédélec space ( ℎ ) (cf. Boffi et al., [9] Corollary 2.3.2) with boundary conditions n × h,z = 0 on z . Inserting this into (25) and integrating by parts leads to 0 = (P Δ ℎ, , J( ℎ, )F(u ℎ )) = (P Δ ℎ, , 1 where we used the fact that curl 1 = curl 2 = curl 3 = 0. It can be shown that (27) can only hold for all h,z if ∇ 1 = ∇ 2 = ∇ 3 = 0 in the following way: In the lowest-order case k = 1, one may insert as test functions h,z with tangential component h, z ⋅ t E ≡ e i for i = 1, 2, 3 on an interior edge E ⊂ z ⧵ z and ℎ, ⋅ t ′ ≡ 0 on all the other interior edges E ′ . If (#E) z denotes the number of interior edges in z , this gives 3(#E) z linearly independent conditions for the 3(#E) z constant values (∇ i ) ⋅ t E for i = 1, 2, 3 on all interior edges E. Therefore, (27) implies that the tangential derivatives of 1 , 2 , and 3 vanish along all interior edges E which implies that 1 , 2 , and 3 are themselves constant. For the higher-order case, each increase of the polynomial degree from k − 1 to k gives additional degrees of freedom to be controlled: for each of the three components, one per edge (including edges on z ), additionally k − 2 per face (including faces on z ) and additionally (k − 2)(k − 3)/2 per tetraeder. These degrees of freedom are forced to vanish by the additional test functions available in the Nédélec space ( ℎ ), see (Boffi et al., [9] Proposition 2.3.5) so that (27) still implies that 1 , 2 , and 3 must remain constant. Finally, the fact that 1 , 2 , and 3 need to be constant implies that (26) can be written as (24). 2nd Step. Inserting (24) into (25) and, restricting ourselves to P Δ ℎ, ∈ Δ ℎ, with, in addition to [[P h, z ⋅ n]] S = 0 for all ∈  ℎ, , P h,z ⋅ n = 0 on all of z (which is automatically satisfied if z ∩ Γ D = ∅), integration by parts leads to with an arbitrary constant a ∈ R 3 . The range of the divergence operator satisfies {div P Δ ℎ, ∶ P ℎ, ∈ Δ ℎ, with [[P Δ ℎ, ⋅ n]] = 0 for all ∈  ℎ and P ℎ, ⋅ n = 0 on } = {z ℎ, ∈ Z ℎ, ∶ (z ℎ, , e ) = 0 for = 1, … , } ≕ Z 0 ℎ, .

Remark 1.
In the linear elasticity case, Proposition 1 turns into the corresponding result from our earlier work, [1] where Basic linear algebra tells us that the right-hand side of the linear system (18) is in the range of the operator  h, z if it is orthogonal to R ⊥ ℎ, , the null space of  * ℎ, . Using Proposition 1 this is obviously the case for patches z with | z ∩ Γ D | > 0 since R ⊥ ℎ, only contains zero in that case. In the case of interior patches z in the sense that | z ∩ Γ D | = 0, we may insert the representation of R ⊥ ℎ, into the right-hand side of (18). This leads to for all (z h,z , s h,z , h,z ) ∈ R ⊥ ℎ, . The first two terms vanish since (P ℎ (u ℎ )), ∇(  ,0 ℎ, )) = (P ℎ (u ℎ )), ∇(  ℎ, )) = (f,  ℎ, ) holds due to the definition ofP ℎ (u ℎ ) as projection onto piecewise polynomials of degree k and the Galerkin condition (6) which holds for piecewise polynomials of degree k + 1 (of which  ℎ, is a fine specimen). Therefore, for each ( we end up with the expression for the inconsistency of the right-hand side in (18). This motivates the choice of a modified test space such that the term in (30) actually vanishes which will be done in the next section.

A MODIFICATION LEADING TO EQUILIBRATED STRESSES
Our construction so far is based on using the simple component-wise L 2 (Ω)-projectionP ℎ (u ℎ ) =  ℎ P(u ℎ ) onto the space of piecewise polynomials of degree k. Due to the incompatibility of the right-hand sides in the local equilibration systems (18) on interior vertex patches, these problems do not possess a solution, in general, as we have shown in Proposition 1. It is certainly possible to solve these systems in a least-squares sense but that would mean that we do not get equilibrated stresses from this procedure. In particular, this means that momentum conservation would not be satisfied locally on each element. We will therefore take up our findings from Section 4 and derive a modification ofP ℎ (u ℎ ) such that the right-hand side in (18) becomes compatible. In view of (29), it is reasonable to choose the test space Z h,z as well as the test functions on the sides ∈  ,ℎ in such a way that they contain the rigid body modes ∈ RM(u h ) of the deformed configuration. With this choice,  ℎ, =  ℎ, = and the sum over the sides in (29) vanishes. A straightforward way to do this consists in building the test spaces on the basis of piecewise polynomials in the deformed variables (x) = x + u h (x) instead of x. This choice also makes sense in view of the fact that the quantities div P ℎ and [[P ℎ ⋅ n]] which are actually tested are mappings from the reference configuration to (forces in) the current configuration. In fact, this modification of the test spaces is only needed for the subspace of polynomials of degree 1 and one can use a hierarchical construction where the enrichment to polynomials of higher degree is again based on the reference coordinates from x.
Let us assume, for the moment, that also the test space in the Galerkin formulation (6) would contain the rigid body modes of the deformed configuration. Then, the compatibility condition in (29) would turn into ifP ℎ (u ℎ ) is defined as the L 2 ( z )-orthogonal projection with respect to piecewise polynomials in the deformed coordinates. The compatibility term in (31) does indeed miraculously cancel out, if z is assumed to be in the test space of the Galerkin formulation (6). Using such a test space is not as far-fetched as one might think. It would ensure invariance with respect to the rigid body modes in the deformed configuration which is not fulfilled for the use of standard polynomial-based finite elements. However, one may argue that such an approach is too complicated for practical use. This is the motivation for the derivation of a suitable choice forP ℎ (u ℎ ) leading to a compatible right-hand side in the absence of this ideal situation. We consider the following slightly more general formulation of the minimization problem (18): , ⟩ for all ∈̃1( ) , ∈  ℎ, , , J( ℎ, )) for all ℎ, ∈ X ℎ, .
An appropriate choice would bef =  ℎ f, the construction ofP ℎ (u ℎ ) will be explained below. As test space in the first equation of (32), could be chosen, where again denotes the mapping from the reference to the (approximated) deformed configuration given by (x) = x + u h (x). The test space for the second equation in (32) would then be given component-wise by transformed polynomials of the form̃( ) = {q• ∶ q ∈ ( ( ))}.
However, in order to make sure that the rigid body modes associated with the deformed configuration RM(u h ) are contained in the test space, it is sufficient to replace the original undeformed rigid body modes RM(0) in the piecewise polynomial test space by RM(u h ) and leave the remaining part unchanged. The test space X h,z for the third equation in (32), the weak symmetry condition, may remain unchanged since only constant rotations appear in the compatibility conditions resulting from Proposition 1. For these spaces, the compatibility condition for all (z h,z , s h,z , h,z ) ∈ R ⊥ ℎ, is therefore, due to (29), equivalent to for all ∈ RM(u h ). A sufficient condition for (36) to hold for all ∈  ′ ℎ is that is satisfied for all ∈ RM(u h ), for all z ∈ T, and for all ∈  ℎ . On each element ∈  ℎ , (37) constitutes d(d + 1) 2 /2 conditions (d(d + 1)/2 rigid body modes in RM(u h ), d + 1 vertices). These conditions can be fulfilled byP ℎ (u ℎ ) ∈ ( ) × of dimension at least d 2 (d + 1) corresponding to k = 1 (which exceeds the number of conditions). A reasonable elimination of the spare degrees of freedom consists in minimizing ‖P ℎ − P(u ℎ )‖ among allP ℎ ∈ 1 ( ) × satisfying the constraints (37). Note that these are separate low-dimensional minimization problems on all elements ∈  ℎ . Thus, our final reconstruction procedure on a vertex patch consists in solving the minimization problem (32) with the test spaces defined in (33) and (34) and withP(u ℎ ) constructed as above.
We end this section with a remark on the inf-sup stability of the system (32) which follows along the same lines as by Boffi et al. [27] for the linear elasticity formulation. It is easy to see that the null space associated with the first and second equation in (32) remains unchanged by the modification of the test spaces, that is, where h,z is the subspace of Nédélec elements (of the first kind) on z with vanishing tangential trace on z . All that is left to show for the inf-sup stability of (32) is therefore that holds with a constant > 0. If we define ℎ, ∶ ( ) → R × by h,z • = h,z F(u h ) −1 , then, according to the transformation rule of the curl operator (cf. Boffi et al., [9] section 2.1.3), h,z ∈ H(curl , ( z )) and where curl denotes the curl with respect to the mapped coordinates. The inf-sup condition (40) is therefore equivalent to the existence of a constant > 0 such that with the mapped Nédélec space h,z holds. This is exactly the inf-sup condition for the original spaces from Boffi et al. [9] in mapped coordinates using parametric Raviart-Thomas elements [28] for the stress approximation. The combination of the inf-sup stability of the system (32) with the fact that our right-hand side is guaranteed to be in its range ensures that there is a correction P Δ ℎ, in the broken Raviart-Thomas space leading to an equilibrated stress P ℎ in the end.

IMPROVED APPROXIMATION OF SURFACE TRACTION FORCES
One of the motivations for the construction of equilibrated stresses is that this leads to approximations of the surface traction forces with an ensured convergence rate. The divergence theorem implies that holds for all v ∈ H 1 (Ω) d . If we assume that (P − P ℎ ) ⋅ n = 0 on Γ N and div(P − P ℎ ) = 0 in Ω holds (eg, since f and g are piecewise constant), then (43) turns into This implies that is satisfied which means that the approximation of the surface traction forces, measured in the H −1/2 (Γ) norm, converges at least as fast as the stress approximation in the L 2 (Ω) norm. Since, by construction, ||P ℎ − P ℎ (u ℎ )|| is expected to be locally an O(h 2 )-approximation, the term on the right-hand side in (45) will converge at the same order as ‖P − P(u h )‖, in general. If we insert v ∈ RM(u h ), the rigid body modes in the deformed configuration, into the numerators in the middle of (45), then ⟨(P − P ℎ ) ⋅ n, v⟩ Γ = (P − P ℎ , ∇v) = ((P − P ℎ )F(u ℎ ) , (∇v)F(u ℎ ) −1 ) = 0 (46) since vF(u h ) −1 = J( ), which constitutes a global version of (24), and (P − P ℎ )F(u ℎ ) is weakly symmetric in the sense of (32).

COMPUTATIONAL RESULTS
We tested our stress equilibration procedure based on (32) for the well-known Cook's membrane example with a quadrilateral geometry. The corners of the domain are located at (0, 0), (0.48, 0.44), (0.48, 0.6), and (0, 0.44) and the boundary is divided into the left line segment Γ D and the lower, right, and upper segments which together form Γ N . Figure 1 shows this geometry and the triangulation  3 which is the result of three levels of uniform refinement. The surface traction force on the right boundary segment is g = (0, ) T with different values > 0, while the upper and lower boundary parts are traction-free; the volume forces f are set to zero. In order to test the robustness of our approach with respect to the incompressibility, we set = 1 and = ∞ in the Neo-Hookean law (4) and use the displacement-pressure approximation from (7) as starting point for our stress equilibration procedure. This means that our stress equilibration procedure really starts from P(u h , p h ) andP ℎ (u ℎ , ℎ ) but the dependence on p h does not change anything in the algorithm. All our computations are for the lowest-order case k = 1 using the Taylor-Hood combination of finite element spaces.
Of particular interest is the distribution of the traction forces on the left boundary Γ D including the singularity with infinite stress components at the upper left corner. The approximation of this distribution of the normal traction force along the left F I G U R E 1 Cook's membrane and triangulation  3 after three uniform refinement steps F I G U R E 2 Normal traction forces n ⋅ (P(u ℎ , ℎ ) ⋅ n) (left) and n ⋅ (P ℎ ⋅ n) on Γ D for  5 ( = 0. 2) boundary is shown in Figure 2. The vertical axis corresponds to the location on Γ D while the horizontal axis represents the value of the (dimensionless) normal traction force. The graphs in this figure are for a load value of = 0.2 on the triangulation  5 which results from two further uniform refinements of  3 . The left graph shows the values for n ⋅ (P ℎ (u ℎ , ℎ ) ⋅ n), corresponding to the projected Piola-Kirchhoff stress from the Galerkin approximation. The right graph shows n ⋅ (P ℎ ⋅ n) for the reconstructed stress. Both pictures represent piecewise affine traction force distributions along the vertical axis. At a first glance, one may get the impression that the left distribution "looks better than" the right one. However, at closer inspection it becomes obvious that the reconstructed stress in the right graph is better able to represent the singular behavior at the upper end. More importantly, the surface forces obtained from the reconstructed Piola-Kirchhoff stress P ℎ recover the correct resultant force obtained from integrating along Γ D . This is a consequence of the divergence theorem which implies .
(48) The approximations  , (P ℎ (u ℎ , ℎ )) and  , (P ℎ ) are shown for the two triangulations  3 and  5 and several values of in Tables 1 and 2. Apparently, the values produced byP ℎ (u ℎ , ℎ ) are not exact but decrease as the mesh is refined. On the other hand, those coming from the stress reconstruction differ from the correct value zero only in the range of machine precision. Tables 1 and 2 showed the approximation quality of very specific quantities, namely the integrated normal force along Γ D , of which the exact value is known. However, one would rather be interested in the convergence of the approximations to the normal force distribution along Γ D with respect to a suitable norm. Table 3 compares the convergence of ‖(P − P ℎ ) ⋅ n‖ −1∕2,Γ versus ‖(P −P(u ℎ , ℎ )) ⋅ n‖ −1∕2,Γ on a sequence of meshes. Since we do not know the exact values of P ⋅ n on Γ D , we access the convergence behavior by the computation of ‖(P ℎ − P 2ℎ ) ⋅ n‖ −1∕2,Γ and ‖(P(u ℎ , ℎ ) −P(u 2ℎ )) ⋅ n‖ −1∕2,Γ , respectively. The norm is evaluated approximately by where V h denotes the space of continuous piecewise linear functions on  ℎ . The values in Table 3 indicate that the convergence for the equilibrated stresses is quite a bit faster than the O(h )-behavior with ≈ 0.544 expected from the regularity of the problem. It can also be seen that the convergence rate is much higher than the one obtained for the original stresses. The reference and the deformed configuration are shown in Figure 3 for = 0.2. The picture clearly indicates that this example is well inside the geometrically nonlinear regime.

CONCLUSIONS
In this paper, a stress equilibration procedure for hyperelastic material models was proposed and investigated. It is necessarily based on a weakly symmetric stress formulation and treats geometrically and materially nonlinear elasticity problems. Our main contribution is the identification of the subspace of test functions perpendicular to the range of the equilibration system on vertex patches not attached to the Dirichlet boundary. This result is used to propose an appropriate test space for the reconstruction of the Piola-Kirchhoff stress and a suitable projection in order to get compatible patch problems. For the moment, this stress F I G U R E 3 Reference and deformed configuration for = 0.2 equilibration procedure is used for its own sake, for example, in order to obtain better approximations of traction forces. Our future goal will be to develop an a posteriori error estimator on the basis of stress equilibration for hyperelastic material models.
Clearly, this will only be possible under restrictive assumptions excluding all the known situations where uniqueness of the solution does not hold.