Extension of the enhanced assumed strain method based on the structure of polyconvex strain‐energy functions

In this work, two well‐known approaches for mixed finite elements are combined to render three novel classes of elements. First, the widely used enhanced assumed strain (EAS) method is considered. Its key idea is to enhance a compatible kinematic field with an incompatible part. The second concept is a framework for mixed elements inspired by polyconvex strain‐energy functions, in which the deformation gradient, its cofactor and determinant are three principal kinematic fields. The key idea for the novel elements is to treat enhancement of those three fields separately. This approach leads to a plethora of novel enhancement strategies and promising mixed finite elements. Some key properties of the newly proposed mixed approaches are that they are based on a Hu‐Washizu type variational functional, fulfill the patch test, are frame‐invariant, can be constructed completely locking free and show no spurious hourglassing in elasticity. Furthermore, they give additional insight into the mechanisms of standard EAS elements. Extensive numerical investigations are performed to assess the elements' behavior in elastic and elasto‐plastic simulations.

in the early 1990s. Basis for EAS elements is a Hu-Washizu type variational functional and the key idea is to enhance the compatible deformation gradient (computed from the deformations) with an incompatible part including the additional degrees of freedom (DOFs). This concept is historically also a mathematically sound justification for the earlier proposed popular incompatible mode models. 23 EAS elements posses many advantageous properties, for example, they are relatively insensitive to mesh distortion * and fulfill the patch test if simple conditions are met. Another huge benefit of EAS elements is their strain driven format which allows simple implementation of complex material laws. Furthermore, they can be constructed completely locking free, which was first achieved by Simo et al. 4 That approach has recently been improved by Pfefferkorn and Betsch 19 in order to cure the violation of the patch test due to the modifications by Simo et al 4 while maintaining its favorable locking behavior.
The main disadvantage of EAS elements are spurious hourglassing instabilities in geometrically nonlinear simulations, which are already mentioned in the first publication by Simo and Armero. 3 Reese 26 and Wriggers and Reese 27 were the first to thoroughly cover this issue and since then there has been a plethora of approaches trying to overcome this drawback. The probably most successful procedure was proposed by Korelc and Wriggers 6 and corrected for frame-invariance by Glaser and Armero. 9,28 Their approach uses the transpose of the ansatzmatrix for the initially employed Wilson-modes. 29 Through this modification, EAS elements become hourglassing-free if material models without instabilities (such as hyperelastic models with polyconvex strain-energy function) are employed. Unfortunately, there is, with the exemption of simple deformation states such as plane strain uniaxial compression, no mathematically sound proof for unconditional stability of the EAS element proposed by Glaser and Armero. 9 However, extensive numerical simulations by many authors 6,9,10,15,19 strongly suggest that result for stable material laws. A remaining major issue are instabilities that arise if the material model exhibits instabilities. This is, for example, the case in tension for popular elasto-plastic material models as those proposed by Simo 30,31 and Simo. 32 Those instabilities are, for example, covered for EAS elements in the works of Armero 10 and Korelc et al. 15 Other attempts to remove the instability include artificial stabilization parameters 9,17,28,[33][34][35] or employ combination with other mixed finite elements. 10,21,22,36,37 Furthermore, an interesting, mathematically more formal approach for a stable mixed element based on the Hu-Washizu functional is given in Angoshtari et al 38 and Faghih Shojaei and Yavari. 39 A drawback of this method is the high number of additional DOFs. However, there is still no unconditionally stable EAS element without other major drawbacks.
Besides the EAS method, the second fundamental concept for this work is the recently introduced framework for mixed elements based on the structure of polyconvex strain-energy functions. This concept was first proposed by Schröder et al 40 and systematically developed by Bonet et al. 41 It employs the tensor cross product proposed by de Boer, 42 which was recently reintroduced by Bonet et al, 41 in order to redefine the cofactor and determinant and greatly simplify computations. Since its introduction, the tensor cross product and the framework for polyconvexity have been used in many publications, including. [43][44][45][46][47][48][49][50] Solid finite elements based on this framework are usually derived from a seven-field-functional, which uses three strain-like and three stress-like fields in addition to the deformation as primary variables. The concept furthermore employs separate discretization of all seven fields. 45 Mixed elements on this basis are extremely robust 50 and allow simple construction of structure preserving time integrators. 44 Among others, a specific set of three kinematic fields is given by the deformation gradient, its cofactor and determinant which are well-known to govern the maps of infinitesimal line, area and volume elements during the deformation process. Those three fields furnish the foundation for the novel enhancement strategies presented in this work.
The key idea for the newly proposed mixed finite elements presented subsequently is to combine the two concepts: The EAS approach and the framework for polyconvex elasticity (see also the preliminary work 51 ). Thus, instead of enhancing the deformation gradient alone and computing other kinematic quantities from it, we treat enhancement of all three kinematic fields, which arise in the framework for polyconvexity, separately. This gives further insight into underlying mechanisms of EAS-elements and opens up a plethora of new avenues for enhancement. In particular, selective enhancement of the kinematic fields and separate enhancement of the cofactor and determinant are examined in the present work. Those approaches lead to promising novel mixed finite elements with many desirable properties.
For that purpose, the present work is divided into seven sections. Section 2 covers basics and fundamentals, such as the tensor cross product, kinematic fields and the material model used. Subsequent sections cover the novel classes of elements starting with selective enhancement in Section 3. To that end, a brief overview of the classical EAS method is given *Within the limits of MacNeal's theorem. 24,25 in Section 3.1. The variational framework for selective enhancement is covered in Section 3.2 and Section 3.3 introduces FE-approximations for the deformations as well as the enhanced fields. The remainder of Section 3 covers orthogonality of stress and strain and the patch test criterion (Sections 3.4 and 3.5). Further novel classes of elements, which are based on additional cofactor and determinant enhancement, are introduced in Sections 4 and 5, respectively. Their presentation relies on many concepts presented in Section 3 and is structured analogously. Extensive numerical investigations follow in Section 6, where standard benchmarks in elasticity and plasticity are presented to assess the performance of selected well-working novel elements and compare them to well-established elements. Finally, a summary and concluding remarks are given in Section 7. The appendix gives additional information on the topics covered in the main matter. Appendix A elaborates on the tensor cross product presented in Section 2.1.1 and Appendix B covers details of the employed material models. Appendix C includes a novel eigenvalue analysis for 3D eight node elements, which extends the 2D method presented by Glaser and Armero. 9 Results of that analysis are given in Section 6.4. As last part of the work, Appendix D gives some details on the performance of all elements tested, that are not included in the detailed evaluation in Section 6.

CONTINUUM MECHANICS
This section covers some basics and fundamentals concerning continuum mechanics, which are used throughout this publication. In particular, aspects of kinematics employing the recently reintroduced tensor cross product [41][42][43] are presented. A key concept of this approach is considering the kinematic fields for the transformation of infinitesimal line, area and volume elements as independent fields, which was inspired by the structure of polyconvex strain-energy functions. 40,41,43,44 This approach motivates the new EAS methods proposed in the remainder of this work.

Tensor cross product
A key algebraic concept in this work is the tensor cross product, which was proposed in the 1980s by de Boer 42 and has recently been rediscovered by Bonet et al. 41 Since then it has been employed in various publications including, 43-48 since it greatly simplifies notation and implementation of large deformation solid mechanics. The tensor cross product between two arbitrary second-order tensors A and B is defined by where ijk denotes the permutation symbol † and the summation convention applies for pairs of repeated indices. Usually, the cofactor of a second-order tensor A is defined by cofA = det(A)A −T . By employing the tensor cross product, this can be recast in the form Furthermore, the determinant can be rewritten as where the double contraction between two arbitrary tensors A ∶ B = A ij B ij is used. A more comprehensive overview of properties and definitions of the tensor cross product is given in Appendix A. † Also known as Levi-Civita symbol.
F I G U R E 1 Reference and deformed configuration of a deformable body and associated kinematic fields F , H and J

Kinematic fields
A deformable body's motion (see Figure 1) from its reference configuration  0 ∈ R 3 to a deformed configuration  (also known as spatial configuration) is described by a bijective deformation map ∶  0 → R 3 , which maps material points X ∈  0 to their spatial counterparts x ∶= (X) ∈  ‡ . The deformations are prescribed on a portion  0 of the boundary  0 by ∶  0 → R 3 . Overall, the map describing the motion of the body  is given by During the deformation process infinitesimal line, area and volume elements denoted by dX, dA and dV are mapped from the reference configuration to corresponding elements dx, da and dv in the spatial configuration (see Figure 1). These maps are described by and governed by the three kinematic fields F (X), H (X) and J (X) which denote the deformation gradient, its cofactor and determinant. The well known relations to compute those three fields can by employing the tensor cross product (1) be recast in the form In the equations above, index of the fields indicates that the corresponding kinematic quantity is solely computed from the deformation field . § Remark 1. The kinematic relations presented so far can also be used as measures of strain and will be used throughout this publication. However, they are not symmetric which is disadvantageous especially for numerical implementation. Therefore, many formulations rely on symmetric measures of strain, for example, the right Cauchy-Green tensor denoted ‡ Note that for simplicity and readability the arguments of functions are subsequently often omitted. § This will not be the case in subsequent sections.
by C . This second-order tensor its cofactor G and determinant C can be used similarly to the fields (6). According to, for example, Betsch et al 44 they assume the forms in the current framework. Corresponding symmetric spatial fields such as the left Cauchy-Green tensor b , its cofactor g and determinant c = C can be defined analogously.

Material model
The material model used throughout the theoretical parts of the present work is a general homogeneous, isotropic, hyperelastic model based on a polyconvex strain-energy function of the form where all three kinematic fields ¶ F, H and J are separate arguments of the polyconvex function W. 41 Note that hyperelasticity is no prerequisite for the mixed finite element methods presented in Sections 3-5, but only employed to simplify notation. The hyperelastic and elasto-plastic material models actually used for the numerical simulations in Section 6 are described in Appendix B.

EAS METHOD WITH SELECTIVE ENHANCEMENT OF KINEMATIC FIELDS
This section covers the first of three extensions of the classical enhanced assumed strain (EAS) method we propose in this work. To that end we briefly summarize important principles of the classical EAS method in Section 3.1 and deal with the novel approach termed selective enhancement in subsequent sections. After the variational framework for that method is introduced in Section 3.2, other vital properties and relations concerning approximation of the enhanced deformation gradient, orthogonality of the discrete stresses and the patch test requirement are covered in Sections 3.3-3.5.

Classical EAS method
The classical EAS method was first introduced for nonlinear kinematics by Simo and Armero. 3 Its key idea is to recast the deformation gradient in the form where the F andF denote the compatible and enhanced part of the deformation gradient, respectively. Therein, the compatible part is defined in Equation (6a) and is computed from the deformations alone. A general form of the enhanced part of the deformation gradient isF which was proposed by Simo et al. 4,9,19 Thus,F is assumed to be a function of the deformations ∈  and n enhF enhanced The basis for EAS finite elements is a Hu-Washizu 52 type variational functional of the form ¶ Index is omitted here since the material model is evaluated with the actual kinematic fields as introduced in the remainder of this work.
Therein, P is a Lagrange multiplier representing the stresses and ensures that conditionF = 0 is fulfilled in a continuous setting. However, in a discrete settingF = 0 is only approximately fulfilled which can ultimately lead to improved numerical solutions.

Variational framework
A variational framework for the novel class of selectively enhanced EAS-elements is presented in this section. For that matter, we start by introducing enhanced kinematic fields in addition to the fields (6), which are based on displacements alone. Thus, all strain-type fields used in this section are given by where fields based on displacements are denoted by index and enhanced fields are identified by index E. The enhanced version of the deformation gradient F E is computed analogously to (9) with the enhanced part of the deformation gradient F given in (10). The choice, whether the standard or enhanced versions of the fields presented in (12) are employed, is described by respectively. They assume value = 0 if the respective field is in its displacement based form and = 1 if the enhanced version is selected.
Remark 2. Note that setting F = H = J = 0 recovers the displacement formulation and F = H = J = 1 relates # to the EAS-element presented by Simo and Armero. 3

Variational functional
The Hu-Washizu type variational functional, on which the modified elements in this section rely, is given by where  0 refers to the body defined in Section 2.1.2 and the deformations is given in (4). The modified strain-energy function W mix is based on the strain-type fields (12) as well as strain-energy function (8) and reads Furthermore, conditionF = 0 is enforced to ensure that the displacement formulation can be recovered in a continuous, non-discrete setting. This is done in (14) similar to (11) with a stress-like Lagrange multiplier Finally, external forces acting on the body are assumed to be dead loads for simplicity. The associated potential reads # Depending on the shape functions used forF.
with prescribed body forces b acting on the volume of  0 and prescribed stresses t on the Neumann boundary  0 ⊂  0 . Conditions  0 ∪  0 =  0 and  0 ∩  0 = ∅ for the boundaries apply as usual. 1

Stationary conditions
Stationary conditions of functional (14) with respect to the variables ∈  ,F ∈ andF ∈  F are derived in this section. The corresponding admissible variations are ∈ ,F ∈ andF ∈  F with the space of admissible variations of the deformation given by Application of the chain rule yields together with definition (10) sinceF was introduced as function of and . With this information at hand, the stationary conditions of (14) are given by The two constitutive stress tensorsP E andP introduced above summarize derivatives of the strain energy function and readP Remark 3. Note that (20c) impliesF = 0. Inserting this result into (20a) and (20b) reduces (20) to a pure displacement formulation, showing consistency with usual continuum mechanics. However, this result is only valid in a continuous setting. Employing enhancement withF in a discrete regime leads to a different, possibly improved numerical solution.

FE-approximation
An approximate solution of (20) can be found with the finite element method (FEM). To that end, body  0 is split into n el finite elements Ω e (see Figure 2), where in the present work, only 4-node quadrilateral (2D) and 8-node brick (3D) elements are considered. The corresponding reference elementsΩ ∶= [−1, 1] n dim are a bi-unit square and cube, respectively. 1 In subsequent sections, shape functions, etc. are only given for 3D elements. However, they can easily be reduced to their 2D equivalents.

Geometry and compatible deformation gradient
We approximate geometry and compatible deformations using the same shape functions, which is referred to as isoparametric concept. The approximation ‖ of geometry within one finite element Ω e is given by where X e I are the nodal reference coordinates of the element and N I denote the standard tri-linear, lagrangian shape functions. They are defined on the reference element in coordinate system = [ , , ] T and read Therein, ( I , I , I ) are the coordinates of the eight vertices of the bi-unit cube. Defining discrete deformations h,e analogously to (22) yields where e I are the nodal deformations. As usual, a Bubnov-Galerkin method is applied, meaning that the variations are approximated in the same way as the deformations . This yields h,e and e I (X e I ) = 0, X e I ∈  h,e for the discrete ansatz spaces of the deformations and variations thereof (for more details see, eg, Wriggers 1 ). With these approximations at hand, the transformation between the reference element and material configuration of a finite element Ω e (see Figure 3) is described by which denote the Jacobian matrix J h,e of the transformation and its determinant j h,e , respectively. 1,2 Finally, the compatible deformation gradient (6a), mapping between material and spatial configuration (see Figure 3), assumes the simple form where the Jacobian (27) is used to compute derivatives of the shape functions with respect to X from the derivatives with respect to via well-known relation

Enhanced deformation gradient
In this work, we approximate the enhanced deformation gradient (10) in accordance with Glaser and Armero. 9 Their general structure for the enhanced deformation gradient is given bỹ which is an element-wise ansatz in order to enable static condensation and therefore a more efficient element because of Moreover, the enhanced deformation gradientF is split multiplicatively into a part F dependent on the enhanced parameters and F 0 dependent on the deformations. The latter, namely F 0 , is the compatible deformation gradient evaluated at the element centroid = 0 and needed for objectivity (see, eg, Pfefferkorn and Betsch 19 or Glaser and Armero 9 for more details). By introduction of ∇ 0 N I ∶= ∇ X N I | =0 tensor F 0 can be cast in the form which has the same structure as (29). The second tensor F e in (31) includes the enhancement parameters . A general form for arbitrary shape functions is where I are the n enhF enhancement parameters per element and M I are shape functions already transformed from the reference element to the physical space. They are computed from local shape functionsM I ( ) defined onΩ via transformation This relation was proposed by Glaser and Armero 9 and is superior to the original, slightly different transformation proposed by Simo et al. 4 In the last equation the expressions J 0 and j 0 are defined as J 0 ∶= J h,e ( = 0) and j 0 ∶= j h,e ( = 0), (36) which are the jacobian matrix J h,e (27) and its determinant j h,e (28) evaluated at the element centroid. This is necessary in order to fulfill the patch test as shown in Section 3.5.
Remark 5. The transformations presented in (33) and (35) are not the only possible ones for EAS elements. In fact, there is a plethora of possibilities as shown by Pfefferkorn and Betsch. 19 However, the transformations used in this work exhibit good performance, which was proven in the aforementioned reference.
An important requirement on the space of the discrete enhanced gradient  h is that conditioñ is fulfilled. This means that there is no intersection of the enhanced deformation gradient with the space of Grad[ h ], which contains the gradient of h,e and therefore the compatible deformation gradient F h,e . Condition (37) ensures stability and avoids rank deficiency of the resulting finite element in the linear case. 2,53 Unfortunately, condition (37) is not sufficient for unconditionally stable elements in the geometrically nonlinear case. Non-physical hourglassing instabilities of the original EAS element proposed by Simo and Armero 3 have for instance already been mentioned in their seminal work 3 and were first thoroughly described by Wriggers and Reese. 27 In the present work, we consider the Wilson modes first proposed for the EAS method by Simo and Rifai 2 and the transpose thereof proposed by Korelc and Wriggers 6 as shape functions for the enhanced deformation gradient. The latter is chosen since it completely removes hourglassing in compression for hyperelastic materials with polyconvex strain-energy function. 6,10 The classical Wilson modes incorporate 9 enhanced modes ** and are given in the present framework by ] . (38) With this information at hand, especially (29) and (31), the approximation of both the deformation based F and enhanced version of the deformation gradient F E presented in (12a) are fully defined. Additionally introducing as modified gradient of the shape functions, enables to cast the discrete enhanced deformation gradient F h,e E into the form which has again the same structure as (29) and (33) and allows simple implementation. Straightforward computations using (29) and (40) yield the remaining discrete versions of the kinematic fields defined in (12).

Stress-like quantities
The remaining field, which needs to be discretized for a complete discrete form of (20), is the stress-likeF defined in (16). Its approximation is done similarly to (31) and is given by where h,ẽ F are frame-invariant interpolation functions, which have to include at least constant stresses in order to fulfill the patch test (see Section 3.5 and Simo and Rifai 2 ). However, the exact form of h,ẽ F is not needed, because the stress-like field is eliminated via orthogonality condition in Section 3.4.
The discrete constitutive stress tensorsP h andP h E are simply computed by using the analytical definitions (21) and evaluate them with the discrete kinematic fields.

Orthogonality condition
Stress-like field h F is eliminated from the discrete version of (20) by setting  h F L 2 -orthogonal to h . Since both h F and F h are discontinuous across element boundaries, L 2 -orthogonality implies for every finite element Ω e . With this relation at hand, and by defining as usual that the variationsF h,e ∈ h and h,ẽ F ∈  h F employ the same shape functions as the actual fields, equation system (20) can be reduced to in a discrete setting. Those last two equations are the basis for a finite element implementation. The relations covering orthogonality and elimination of the discrete stress-like field presented here are based on the works of Simo et al. 2,3

Patch test requirement
A crucial demand on any finite element is that it fulfills the patch test requirement, meaning that the response of a solid meshed with an arbitrary patch of elements and subject to a constant state of strain is exactly reproduced. 54 Ultimately, this means for the classical EAS method that the shape functions of the enhanced deformation gradient have to fulfill simple conditions. These requirements have first been defined by Simo and Rifai 2 and extended to non-linear kinematics by Simo and Armero. 3 In this work, the presentation of the patch test condition for the newly proposed elements is mainly based on the work of Armero. 10,19 Below we present the necessary and sufficient conditions for selective enhancement elements. We start by noting two important conclusions that can be drawn from the imposed constant strains F h,e = F 0 . First, they imply constant (constitutive) stressesP for homogeneous materials, which are assumed in the present work (see Section 2.2). Secondly, follows † † from condition (37) since the space of compatible deformations  h (25) includes constant strains ‡ ‡ . Inserting (45) into (43) and comparing it with a pure displacement formulation yields relation † † Note that the variationsF h,e ≠ 0 since they are arbitrary. ‡ ‡ It is well known that the isoparametric displacement formulation fulfills the patch test.
for the overall stressesP. Together with (44) this implies that both constitutive stress tensorsP h E = P E,0 andP h = P ,0 are constant under the assumption of constant strains.
Ultimately, (43) has to be fulfilled exactly to satisfy the patch test. Imposing conditions (44), (45) and (46) on (43a) yields for all admissible variations h . This equation is satisfied as it is equivalent to a pure displacement formulation, which fulfills the patch test. Inserting the same relations into (43b) and splitting the integral into integrals across single finite elements yields In a next step, (34) and the transformation (35) are inserted into the last result. Furthermore, the integral is transformed to the reference elementΩ. This leads to for arbitrary I . As F 0 , J 0 and j 0 are constant within an element and j h,e can be cancelled out, the final condition for satisfaction of the patch test emerges as This has to hold for all shape functionsM I of the enhanced part of the deformation gradient and is equivalent to the condition for standard EAS elements. 9 It is simple to show that this condition is fulfilled by (38) and its transpose.

EAS METHOD WITH COFACTOR ENHANCEMENT
We present a second novel class of finite elements based on the EAS method in this section. It includes a separate enhancement of the cofactor (6b) additionally to enhancement of the deformation gradient (see Section 3.1). The newly proposed method relies on many concepts presented in Section 3, which is why only theory different for or unique to cofactor enhancement will be presented in this section.

Variational framework
For the matter of introducing a variational framework for EAS elements based on cofactor enhancement, we introduce the kinematic fields and their enhancement in a first step. They are given by § § § § Note that indices E and denoting enhanced and deformation based fields in (12) are omitted here, since the fields are always enhanced.
where ∈  (4),F ∈ (10) and the enhanced cofactorH is introduced as This additional enhanced field is assumed to be a function of the deformations ∈  and enhanced parameters ∈ { ∶  0 → R n enhH | i ∈ L 2 } in analogy to (10). Note that the determinant (51c) is computed according to (6c) using the enhanced versions of deformation gradient and cofactor. This allows the benefit from cofactor enhancement even for material models not based on the cofactor. Furthermore, there are many possible ways of enhancing the cofactor, three of which are examined in this work as introduced in (51b). The three different enhancement types of the cofactor are denoted by indices. However, these subscripts are omitted if no specific version but the cofactor in general is addressed (see, eg, Equation (51c)). Furthermore, the auxiliary variableH ∈ {1, 2, 3} is used to denote the enhancement options of H. That is, H = H k for k ∈H.

Variational functional
Similarly to (14) a Hu-Washizu type variational functional 52 is introduced in this section. It is a five field functional of the form where the kinematic fields F, H and J are given above and ∈  (4). The Lagrange multiplierF ∈  F and the potential of external forces Π ext are given in (4), (16), and (17), respectively. The second Lagrange multiplier ensures that conditionH = 0 is fulfilled. Together with (16) it allows to recover the displacement formulation in a continuous setting (see Section 4.1.2 below).

Stationary conditions
The stationary conditions of functional (53) are derived by computing the first variations with respect to ∈  ,F ∈ , (19), and the admissible variations ∈ ,F ∈ ,F ∈  F ,H ∈,H ∈  H the stationary conditions are given by Therein, the constitutive stress tensorŝF and̂H have been introduced for notational simplicity. They summarize derivatives of the strain energy function W (8) with respect to the kinematic fields and assume the formŝF Remark 6. Note that (56c) and (56e) implyF =H = 0. Inserting these results into the remaining equations of (56) yields a pure displacement formulation, proving accordance with continuum mechanics. However, this is not possible in a discrete setting (see subsequent sections), where the numerical solution can be improved by (cofactor) enhancement.

Finite-element approximation
This section covers the FE-Approximation of stationary conditions (56). The discretization of geometry X, deformations , enhanced deformation gradientF and stressesF as well as some general concepts are presented in Section 3.3. The relations described there, especially (22), (25), (26), (31) and (41), are applied without changes for the approximation of the current formulation. This leaves only the enhanced cofactorH and the corresponding stressesH to be discretized.

Enhanced cofactor
The discrete enhanced cofactorH h is defined similarly toF h in (31) yielding where e is given in (32) and F 0 in (33). The approximation is again element-wise due to e and F 0 ensures objectivity of the elements. This can be shown by examining the transformation of the cofactor under superposed rigid body motions * = Q + c with rotation matrix Q ∈ (3) and constant vector c ∈ R 3 . From this we get with (A20) for the deformation based version of the cofactor since Q −1 = Q T and det(Q) = 1. A similar relation has to hold for each of the three alternative versions of the enhanced cofactor introduced in (51b) for the overall method to be objective. In particular, the cofactor enhancement introduced in (51b) has to satisfỹ for all three version, which is the case in (58) due to F 0 incorporated there. The remaining field in (58) to be introduced is H e which is assumed to have the form in accordance with (34) and (35). The shape functionsL I for the enhanced cofactorH are defined locally onΩ and again, the Wilson modes (38) and their transpose are used. However, it is also possible to enhance, for example, only the diagonal or non-diagonal entries of H resulting in Additionally to the (Wilson-)modes introduced above, the purely volumetric shape functions are considered as shape functionsL I for H e in the present work. An important condition to ensure stability and avoid rank deficiency of the resulting element is in analogy to (37), which is fulfilled by the present shape functions. We provide this condition here without sound mathematical proof, but numerical investigations suggest the importance of (64).

Stress-like quantities
The last field to be discretised are the Lagrange multipliersH defined in (54). Following the descriptions presented in Section 3.3.3 yields where h,ẽ H is a frame-invariant interpolation containing at least constant stresses in order to fulfill the patch test (see Section 4.4). Again, the exact form of h,ẽ H is not needed since the discrete stress-like fields are eliminated via orthogonality condition in Section 4.3.

Orthogonality condition
The stress-like quantities h F and hH are eliminated from the discrete version of (56) by setting h and h L 2 -orthogonal to  F and  H , respectively. This implies for every element as the enhanced fields and stress-like quantities are defined element-wise. Setting as usualF h ∈ h ,

Patch test requirement
The important patch test requirement is thoroughly described in Section 3.5 and has to be fulfilled by any finite element. 54,55 For cofactor enhanced elements we start by noting similar to (44) and (45) and with the same reasoning that conditions follow from the imposed constant strains D ∶= F 0 . Condition (68a) furthermore implies that the constitutive stresseŝ h F and̂hH are constant. This means that the constitutive stresses can be cast into the form where P 0 ,F ,0 , andH ,0 , are constant. Ultimately, (67) has to be fulfilled exactly in order to fulfill the patch test. Thus, we impose the conditions above on (67a) which yields The last equation is equivalent to the displacement formulation and therefore passes the patch test. Straightforward algebraic manipulations of (67b) and (67c) using conditions (68) and the approximations described in Section 3.3.2 as well as in (58) ultimately lead to requirements on the shape functions used for the enhanced fieldsF andH (cf. Equation (50)). Thus, the same requirements as for the standard EAS method have to hold for both enhanced fields in the novel class of cofactor enhancement elements.

EAS METHOD WITH DETERMINANT ENHANCEMENT
A third and final class of novel finite elements is presented in this section. It is motivated by the relatively good performance of the element employing pure determinant enhancement presented in Section 3. Thus, it is based on separate enhancement of the determinant (6c) in addition to standard enhancement of the deformation gradient (see Section 3). Many concepts already presented in Sections 3 and 4 can therefore be directly transferred to this class of elements.

Variational framework
The variational framework covered in this section is the basis for the novel class of elements called determinant enhancement. Similarly to previously introduced elements we start with the kinematic fields and their enhancement which are given by where ∈  andF ∈ are defined in (4) and (10), respectively. The determinant is enhanced with scalar field which is assumed to be a function of the enhanced parameters ∈ { ∶  0 → R n enhJ |( ) i ∈ L 2 }. Again, the deformation gradient and its enhancement are employed in (72) as given in (9) and (10). This enhanced form of F is then used via relations (6) for cofactor and determinant as well. This means that enhancement of the deformation gradient is done in the same way as for the originally proposed geometrically non-linear EAS element. 3 On top of that, the determinant is enhanced with an additional scalar fieldJ.
Remark 7. Note that the enhanced field for the determinant is not dependent on the deformations. This is because of its scalar form which means that in contrast to enhancement of the deformation gradient (31) and cofactor (58) no objectivity requirements apply.

Variational functional
The five field Hu-Washizu type variational functional, 52 on which the element class at hand is based, has the form Therein, kinematic fields are defined above and Lagrange multiplierF ∈  F is given in (16). The potential of external forces Π ext is defined in (17). Moreover, the additional scalar Lagrange multiplier ensures that conditionJ = 0 is fulfilled. This allows together withF = 0, which is enforced byF, that the displacement formulation can be recovered.

Stationary conditions
The stationary conditions of (74) with respect to the fields ∈  ,F ∈ ,F ∈  F ,J ∈ , Λ̃J ∈  J are given below. With relation (19) andJ =J the stationary conditions of functional (74) read where ∈ ,F ∈ ,F ∈  F ,J ∈ and Λ̃J ∈  J are the admissible variations. The constitutive stress-like quantitiesP andΛ̃J contain derivatives of strain energy function (8) and are given bŷ Remark 8. Note similar to previous sections, that (76c) and (76e) yieldF = 0 andJ = 0 in a continuous setting. Inserting these results into the remaining equations of (76) leads again to a pure displacement formulation ensuring accordance with continuum mechanics.

FE-approximation
The FE-Approximation of (76) is covered subsequently. Many of the approximations presented in Section 3.3 can be used for the element class at hand as well. Note especially the discrete forms of X h (22), h (25), h (26),F h (31) and h F (41) presented there. The remaining fields left for approximation are the enhanced determinantJ and the corresponding Lagrange multiplier Λ̃J.

Enhanced determinant
The enhanced determinantJ is approximated element-wise similar to the enhancement of deformation gradient (31) and cofactor (58). This allows static condensation in order to obtain an efficient element. However, the deformation gradient at the element center F 0 is not needed in contrast to (31) since the determinant is scalar valued and objectivity is thus easily maintained. With that information at hand and e defined in (32) the discrete enhanced determinant can be cast in the formJ where the element-wise enhanced determinant J e is given by It uses the same scaling with the determinants j 0 (36) and j h,e (28) as the transformation for the enhanced deformation gradient given in (35) in order to fulfill the patch test (see Section 5.4).
Two different sets of local shape functionsK I defined on the reference elementΩ are considered and examined within the present work. They are given by which are linear and bi-linear, respectively. Note that the second set of shape functions (81) is inspired by the work of Simo et al. 4 An important requirement for stability of the elements is in analogy to (37) and (64) which is fulfilled by the shape functions presented above.

Stress-like quantities
where Λ h,ẽ J is a frame-invariant function containing at least constant stresses. This last condition is necessary to ensure that the patch test requirement is fulfilled.

Orthogonality condition
The set of Equation (76) is reduced in the discrete case by imposing an orthogonality condition on the stress-like quantities h F and Λ h J . More specifically, the discrete ansatz spaces  h F and  h J are set L 2 -orthogonal to the spaces h and h , respectively. This leads, due to the element-wise definition of the discrete stress-like fields and enhanced strains, to Together with the variations approximated as usual with the same shape functions as the actual fields, H , this ultimately yields the reduced discrete set of equations given by

Patch test requirement
The important patch test requirement described in Section 3.5 is covered in this section for elements with determinant enhancement. The imposed constant strains yield together with (82) conditionŝ in analogy to (44) and (45). Furthermore, condition Λ̃J = const. follows from (86a). Eventually, the discrete set of Equation (85) has to be fulfilled analytically to satisfy the patch test requirement. Imposing conditions (86) on (85a) yields in the same manner as in (47), which is equal to a pure displacement formulation. Thus, (87) is fulfilled since the isoparametric displacement element passes the patch test. Conditions on the shape functions of the enhanced fields arise from imposing (86) on the remaining equations of (85). They are of the form where straightforward algebraic manipulations using the approximations presented in Section 3.3.2 as well as 5.2.1 have been applied. A more detailed description of the steps used to reach the conditions above is given in Section 3.5.
Remark 9. Note that the patch test requirement is the reason for incorporating the fraction with Jacobian-determinants in (79) as j h,e cancels out after transformation of the integral from domain Ω e toΩ.

NUMERICAL INVESTIGATIONS
We present standard numerical tests in elasticity and plasticity to evaluate the performance and properties of the novel elements proposed in Sections 3, 4 and 5. Objectivity, sensitivity to mesh distortion, convergence behavior, stability of the elements (hourglassing) and the patch test requirement are some of the properties of interest. Furthermore, locking in bending dominated situations as well as in the incompressible limit are examined.
To that end, the considered elements are listed here before we present the numerical simulations in Sections 6.1-6.9. If not stated otherwise all elements employ the standard eight point (four point in 2D) quadrature rule. The newly proposed elements are in all tests compared to a set of standard elements denoted by: • H1: isoparametric displacement formulation 8-node element. The 2D ¶ ¶ version is denoted Q1.
• H1/E9: EAS element proposed by Simo and Armero 3 using the Wilson modes (38) for the enhanced field. In 2D the element has only 4 enhanced modes and is labeled Q1/E4.
• HA1/E12T: Element proposed by Pfefferkorn and Betsch 19 using the transposed Wilson modes and 3 additional enhanced modes together with a special nine point quadrature rule which was proposed by Simo et al. 4 Furthermore, this element uses a special form of the compatible deformation gradient, which improves the locking behavior in distorted meshes (see Pfefferkorn and Betsch 19 for more information). The only difference of the corresponding 2D element QA1/E4T to Q1/E4T is the five point quadrature rule.
• HA1/E12T-F − T 0 : Same as HA1/E12T but using F −T 0 instead of F 0 in (31). This approach was recently proposed by Pfefferkorn and Betsch 19 and improves the bending performance of the element.
Note that H1/E9T was chosen because it completely eliminates spurious hourglassing from the original H1/E9 for polyconvex materials. This instability occurs under pressure and was first thoroughly investigated by Reese 26  where the first type of elements employs Wilson modes (38) for the enhanced deformation gradient and "E9T" indicates use of the transposed Wilson modes. Which kind of selective enhancement is chosen, is denoted by suffix "Sfhj" where "S" denotes selective enhancement and "fhj" represents enhancement of F, H and J, respectively. A zero indicates usage of the displacement form and a one stands for the enhanced version listed in (12). As example to illustrate the naming pattern, a 2D element based on Wilson-modes and solely employing enhancement of the deformation gradient is identified as Q1/E4-S100. Note that the standard elements H1, H1/E9 and H1/E9T are recovered by the elements H1/E9-S000, H1/E9-S111 and H1/E9T-S111, respectively. The second class of new EAS-elements is presented in Section 4 and is based on separate enhancement of the cofactor H in addition to enhancement of the deformation gradient. Identification of these elements follows structure

H1/FH FfHh
where subscript indicates which cofactor enhancement version is chosen from (51b) and suffix "FfHh" identifies the shape functions used for the enhanced field of the deformation gradientF (31) and cofactorH (58), respectively. More precisely, the upper case letter indicates the type of shape functions and the lower case character identifies the structure of the enhancement as defined in Table 1. The 2D version of the elements are labeled with Q1 instead of H1. With this notation at hand, H1/FH -XxXx, H1/FH 3 -WaXx and H1/FH 3 -TaXx are aliases for the standard elements H1, H1/E9 and H1/E9T. ¶ ¶ All 2D elements presented in this work are subject to the plane strain condition. ## Note that an earlier proposed element using the transposed Wilson modes as well but lacks frame invariance was proposed by Korelc and Wriggers. 6 TA B L E 1 Denomination of shape functions for enhanced deformation gradientF and cofactorH

Shape function type
Enhanced entries X No enhancement (n enhF = 0 or n enhH = 0, respectively) The final class of novel elements is presented in Section 5 and is based on enhancement of the determinant in addition to the standard enhancement of the deformation gradient. Naming of these elements follows convention

H1/FJ-FfJ
where "Ff" identifies the shape functions used to approximate F in the same way as for the cofactor enhanced elements (see above and Table 1). The type of shape functions used to enhance the determinant is denoted by the final "J". Abbreviations for the ansatz functions used in this work are listed in Table 2. Again, it is possible to name the standard elements H1, H1/E9 and H1/E9T with the presented structure. They are given by H1/FJ-XxX, H1/FJ-WaX and H1/FJ-TaX, respectively.
Having introduced the structure of names for the novel elements, the selected combinations are presented in a next step. Note that only a few well working elements are included for the detailed evaluation in this section. Refer to Appendix D for a complete list of the 69 tested combinations and a short explanation in which test(s) they fail or perform poorly in. The chosen, well working elements are • H1/E9-S001 • H1/FH 3 -TnVv and their respective 2D versions. Elements presented in Section 3 do in general not perform well, which is why only one thereof is included in this detailed evaluation. Furthermore, hardly any of the well-working elements is based on Wilson-modes. This is because these elements often exhibit instabilities in the stability test (see Section 6.4), which is deemed to be a disqualifying criterion in the present work.
Remark 10. Some of the novel elements are not suitable for material models that do not directly depend on the cofactor H. This is for instance the case for the widely used Neo-Hookean model. A good example for an improper element is H1/E9-S010 which uses only enhancement of the cofactor. If this element is combined with a material model not incorporating H the tangential stiffness matrix for the enhanced degrees of freedom (in the undeformed state) will become singular which ultimately causes failure of the static condensation process. In deformed states, when nonlinear effects have to be taken into account, this issue may possibly not arise but it is still a serious limitation for an element. A list of elements encountering this problem is given in Appendix D. Note that these elements may still exhibit great performance if material models such as the Mooney-Rivlin model (B1) are used. Two material models, which are presented in Appendix B are used for the numerical simulations. First, a Mooney-Rivlin material model with strain energy function (B1) and second an elasto-plastic model presented in-depth in Appendix B.2. The material parameters for the Mooney-Rivlin model are defined for each test separately. However, in all elasto-plastic simulations in the present work we use the standard material parameters given in Table 3. [2][3][4]9,10,15,31,32 Note that we employ the line-search algorithm described in the textbook of Bonet and Wood 56 in elasto-plastic simulations to stabilize the Newton-Raphson algorithm.

Patch test
The first important property of the elements we check in a numerical test, is the patch test requirement. 24,54 In this test a state of deformation is applied such that constant strains occur and the required outcome is constant stress for homogeneous materials (see Section 3.5). It is an important condition, that any finite element has to fulfill. 54 are applied, yielding restraint free bearings. The state of constant strain is introduced by imposing an additional boundary condition on the upper surface X 3 = 1 (analogously in 2D with u 2 and X 2 = 1). Displacement u is increased in steps of Δu = 0.05 until either the Newton-Raphson scheme fails or a displacement of u = 0.95 is reached. The material parameters for the Mooney-Rivlin material (B1) are chosen to a = 1.538 ⋅ 10 5 , b = 7.692 ⋅ 10 4 and c = 2.692 ⋅ 10 5 which correspond to E = 10 6 and = 0.3 in linear theory. The block is discretized both with a regular 3D mesh with 3 × 3 × 3 elements and an initially distorted ‖‖ mesh (see Figure 4). Corresponding 2D meshes are chosen accordingly. Every element tested passes this numerical patch test, which is the expected result since fulfillment of the patch test is proven analytically for all elements: • Analytic investigations for the novel elements proposed in this work are given in Sections 3.5, 4.4 and 5.4. 19 give proofs for all standard EAS elements elements *** .

• Pfefferkorn and Betsch
• H1 is well known to pass the patch test. 54 The numeric fulfillment of the patch test can be verified by looking at the mean value and standard deviation of the von Mises stresses in the domain. At a given deformation state u = 0.95 in the elastic case, all elements are able to reproduce ‖‖ The geometry of the distorted mesh is taken from MacNeal and Harder. 57 ***Note that Simo and Armero 3 and Glaser and Armero 9 proposed the first proofs for element H1/E9 and H1/E9T, respectively. (left) and final configurations (right) for the elasto-plastic simulations with the Q1/E4 element. Reference configuration plotted with dotted lines in plot of final configurations the correct mean value 8.570 ⋅ 10 6 and show standard deviations below 1 ⋅ 10 −6 † † † . Thus, the standard deviation is at least 12 orders of magnitude lower than the mean value which establishes the result. In the elasto-plastic case only about 6 orders of magnitude difference can be observed, which is due to the additional inaccuracy induced by the eigenvalue perturbation technique. 58

Objectivity test
The next test covers another crucial property of finite elements which is frame invariance. It was originally proposed by Glaser and Armero 9 (see also other works 15,19 ). The version shown here is identical to the one investigated in Pfefferkorn and Betsch 19 with the exception of the additional elasto-plastic simulations in the present work. Figure 5 shows the test's geometric setup.
The considered bar has a length of L = 1.0 and a square cross section with h = b = 0.1. Dirichlet boundaries are applied on left and right end of the undeformed beam-like structure  0 , respectively. The final deformations state is reached in n load steps during which both the scalar displacement u i = 2h ⋅ i∕n and the rigid body motions governed by rotation matrix are gradually increased. Note that (91) together with (92) ensures that the rotation path is followed exactly, which is important in elasto-plastic simulations to avoid that spurious plastic strains destroy objectivity. 15 The final angles of rotation are chosen to = 0, 15, 30, 45, 60, 75, 90 • . 9 In elastic simulations the number of load steps can be adapted in a way that fewer steps n = ∕3 + 1 are necessary for smaller . However, the same varying step size cannot be applied for elasto-plastic simulations because of the path-dependency of the material model. In that case it is necessary to maintain a constant n = 60 for all angles to avoid influence of errors from integration of the flow rule ‡ ‡ ‡ . † † † Due to round-off errors. ‡ ‡ ‡ This dependency on number of load steps is, for example, illustrated in Simo 32 F I G U R E 6 Regular (left) and distorted (right) element in 3D 19 The finite element discretization with 6 elements and all final configurations are shown in Figure 5. To complete the setup, the material parameters for the Mooney-Rivlin material (B1) are chosen to a = 40, b = c = 10, which correspond to = 50 and = 100 or equivalently E = 233.3, = 0.1667 in linear theory. In plasticity Table 3 applies.
Frame invariance of the elements is finally verified by examining the axial force N, shear force V and bending moment M in a beam like manner at the bearings. In a coordinate system rotated by , there may of course be no change for the element to be objective. We quantify this behavior by computing the standard deviation of the reaction forces for all final angles . All elements tested pass this test with standard deviations below the Newton tolerance 1 ⋅ 10 −8 for elastic as well as elasto-plastic § § § material. This is inline with the results presented in previous parts of the present work.

Linearized eigenvalue analysis
This purely elastic test is designed to investigate the locking behavior of finite elements in the undeformed state and has first been proposed for EAS elements by Simo et al. 4,8,13,19,59 The simulation is performed in the incompressible limit to include volumetric locking. Since incompressibility cannot be exactly enforced with the elements at hand, it is approximated by choosing for the Mooney-Rivlin material (B1), which corresponds to a ratio of the bulk modulus K to the shear modulus of 1 ⋅ 10 9 in linear theory. Figure 6 shows the distorted and regular one element mesh on which the eigenvalue analysis is performed in the stress free reference configuration. The outcome of the spectral analysis of the stiffness matrix are 24 eigenvalues in 3D and 8 in 2D, which resemble the number of respective displacement DOFs of a single element. No additional modes occur due to the enhanced DOFs since those are eliminated on element level by static condensation. Three groups of eigenvalues i can be formed for every element: First, rigid body modes with i = 0, of which there have to be exactly six in 3D and three in 2D representing the number of independent motions of an unconstrained rigid body. More rigid body modes (or modes with i < 0), which, for example, occur for element-material combinations mentioned in Remark 10, indicate an instability. The second set of modes are locking modes with i → ∞ ¶ ¶ ¶ . An ideal element would have only one locking mode representing pure volumetric deformations. 4,13 Finally, all remaining modes with finite eigenvalues denote soft modes and govern the deformation behavior of the elements.
For some of the tested 3D elements, the results are shown in Table 4. All elements are able to recover the right amount of rigid body modes but, unfortunately, the ideal number of modes in the other categories is in general not obtained. The only elements showing exactly the desired behavior with only one locking mode in both meshes are HA1/E12T and HA1/E12T-F −T 0 . All other elements show more or less severe locking revealed by additional modes with i → ∞. Especially many locking modes occur for the displacement element H1 which shows seven and eight too stiff modes for the regular and distorted mesh, respectively. The standard EAS element H1/E9T has three locking modes less but is, as many of the newly proposed elements, not completely locking free. This is nevertheless a great improvement in practical simulations and is therefore sometimes referred to as mild locking. 4 The best of the newly proposed elements is H1/FJ-TaB which shows superior behavior by fulfilling the regular test and exhibiting only one extra locking mode in the distorted § § § Without considering a constant number of load steps for all elasto-plastic simulations, deviations would go up to approximately ≈ 1 ⋅ 10 −2 . ¶ ¶ ¶ Since incompressibility is only approximately enforced ∞ is not exactly reproduced. Thus, we consider, in this benchmark, modes with "high" eigenvalues as locking modes. test. This additional mode in the irregular mesh could be removed by modifying the compatible part of the deformation gradient as shown by Pfefferkorn and Betsch. 19 In 2D, the only element subject to locking is Q1. All other elements exhibit only one mode with i → ∞ and are thus completely locking free.

Stability analysis
The first stability analysis for EAS elements was presented by Reese 26 to examine the spurious hourglassing behavior of the original EAS element proposed in. 3 This phenomenon was already mentioned in the very first publication on geometrically nonlinear EAS elements by Simo and Armero 3 and has first been thoroughly described by Reese 26 and Wriggers and Reese. 27 It is a major drawback of EAS elements based on Wilson modes, which becomes even more apparent when the recent work by Sussman and Bathe 16 is considered. They show that this instability arises even at sates of small strain if the element's aspect ratio is high. Note that this benchmark is often performed analytically 1,9,10,15,16,26,27 in contrast to the numerical investigations presented here and in de Souza Neto et al. 60 However, both analytic and numeric simulations lead to the same results. The test presented here is mainly taken from the work of Glaser and Armero 9 and Armero, 10 where the stability analysis is conducted on a single unconstrained rectangular element. This is in contrast to the first test proposed by Reese,26 who suggested a test with an element subject to Dirichlet boundary conditions. The drawback of the constrained benchmark is that instabilities of the elasto-plastic simulation cannot be investigated. Furthermore, Armero 9 showed that it is not necessary to consider initially rectangular shapes since it does not change if instabilities occur but only at which level of strain. 16 Thus, we only consider the square element shown in Figure 7 which coincides with the reference element to get especially simple transformations. The deformation state shown in Figure 7 implies a constant diagonal deformation gradient of the form This means, that principal directions of the system coincide with the coordinate axes shown in Figure 7.
A deformed configuration of the given system is fully described by prescribing the stretches 1 in the 2D plane strain case ( 3 = 1), since it allows to compute the missing stretch 2 from condition 2 = 0 and the given material law. Despite the nonlinearity of the problem it is still possible to find an analytic solution for the stretch 2 in case of the Mooney-Rivlin material (B1), whose material parameters are set to a = 9, b = 1, c = 99996 in this test (corresponding to = 20 and = 10 5 in linear theory). For the elasto-plastic material described in Appendix B.2 only a semi-analytic solution is possible since the return-mapping procedure requires a Newton-Raphson scheme. Analogous computations can also be carried out in the uniaxial stress 3D case where again 1 is prescribed and condition 2 = 3 = 0 is used to compute 2 = 3 .
Thus, the complete deformation state of the system is defined by setting 1 and computing the other principal stretches. This (semi-)analytic deformation state is then imposed on the system in order to compute the corresponding tangential stiffness matrix K e of the element.
In 2D, eigenvalues hour i , i = 1, 2 associated with the two hourglass eigenvectors arranged in (see Figure 7) can simply be computed via relation [ hour 1 0 0 Interestingly, (95) are the hourglass eigenvectors of all 2D elements presented in this work regardless of specific formulations or ansatz functions as long as the element's shape remains rectangular, which can be verified numerically and analytically. 9,10 The procedure described above, can only be used for a four node element in a 2D plane strain problem. However, there is a similar, slightly more complex method for the 3D case, which has to the best knowledge of the authors not been proposed so far and is given in Appendix C.
Ultimately, the goal of the present test is to show that no negative hourglass eigenvalues hour i occur, which ensures that no spurious hourglass instabilities arise in the present deformation state. Note that avoiding negative hourglass stiffness in this test does not ensure a completely stable element in all deformation states. However, extensive numerical investigation show at least greatly improved stability throughout many tests. 15 The results of the 2D stability analysis are shown in Figure 8. The hyperelastic simulations reveal that the only element prone to hourglassing is Q1/E4 due to its vanishing eigenvalue hour 2 at 1 ≈ 0.61. All other elements use the transpose Wilson modes which are known to cure this spurious hourglassing behavior of the standard EAS element in compression. 9 Note that displacement element Q1 and mixed element Q1/FJ-TnB are not shown in Figure 8A since their eigenvalues are too high for the range depicted. Such high eigenvalues occur due to locking.
In the elsto-plastic case similar results can be observed in compression. Again only element Q1/E4 shows an instability for 1 < 1. However, in tension, an instability occurs for Q1/E4, Q1/E4T and Q1/FJ-TaB which ultimately leads to unphysical hourglassing patterns in elasto-plastic simulations observed in Section 6.8. This instability occurs due to an instability of the considered elasto-plastic material model in tension, which carries over to the hourglass modes. 9 Results of the 3D stability analysis are shown in Figures 9 and 10 and are computed with the procedure described in Appendix C. First, all 12 modes with nonconstant strains are depicted in Figure 9. In contrast to the 2D simulation, not all of these modes are classical hourglass modes since more complex deformation patterns can also be observed. The usual hourglass modes are in this context mode 1 and 5 to 10. Moreover, three pairs of the twelve modes have the same eigenvalue due to the symmetry of the considered uniaxial stress problem. This concerns pairs 2-3, 7-9, and 8-10, which can be seen from the simple rotation of the corresponding eigenmodes by 90 • around the x-axis.
The actual eigenvalues associated with the 12 nonconstant modes are shown in Figure 10 for the Mooney-Rivlin material. Note that duplicate modes as described above and modes with high eigenvalues (mode 1 and 2 for every element tested here) are not included in the diagrams. Mode 12 exhibits negative stiffness for every element, even for the displacement based H1. Fortunately, this is no problem, since mode 12 cannot lead to global unphysical hourglass patterns due  to its incompatibility with neighboring elements. Again, we observe that the element based on Wilson-modes exhibits instabilities in the compression range while H1/E9T avoids those. Thus, we get similar results as in the 2D analysis. The same can be observed for all other standard elements and novel elements tested in the present work, with the exception of H1/FH 3 -TnVv which exhibits negative eigenvalues in compression for mode 7. However, this instability does not seem to affect the results in any of the other simulations. All in all, the 3D analysis confirms the results of the 2D analysis.

Large mesh stability test
On top of the simple one-element stability analysis presented in Section 6.4, we perform another stability test on larger FE-meshes in this section. The test shown here is similar to the one introduced in Pfefferkorn and Betsch 19 and was inspired by the work of Glaser and Armero. 9 The benefit from this benchmark is that it reveals hourglassing of elements in combination with their neighbors, which naturally cannot be seen in the previously introduced one-element analysis. Moreover, it also shows the ability of the element to depict physically correct instabilities. For this test we consider an elastic cube (square in 2D) with an edge length of a = 50 modeled with the same hyperelastic material as in Section 6.4. It is meshed with 12 elements per side in a regular and distorted manner as shown in Figure 11. The distorted mesh is created by shifting all coordinates of each node in the regular mesh by a random value Δ i ∈ [−1.25, 1.25], such that the outer shape of the elements is not changed ### . Note that the same distorted mesh is used for all elements tested. Dirichlet boundary conditions u X = u Y = u Z = 0 are applied on lower surface, edge and corner as shown in Figure 11. Furthermore, prescribed displacements u Z = u i (compression) are applied on the upper surface. This setup allows again to compute an analytic homogeneous solution as described in Section 6.4. Subsequently, we impose this deformation state on the system and perform an eigenvalue analysis on the tangential stiffness matrix K red in which the constrained DOFs (Dirichlet boundaries) have been eliminated. By gradually increasing the displacement u i in steps of Δ i = 0.01 singular points with eigenvalues i = 0 are approximated until either four such instabilities are found or a displacement of u i = 0.8a is reached. Plotting the mode shapes X i associated with i allows to determine, whether the corresponding instability is physical or unphysical (mainly hourglassing).
### This means that if we assume the square/cube's center to coincide with the coordinate origin, any coordinate with |x i | = a∕2 = 25 remains unchanged.

F I G U R E 12
First two eigenmodes for regular (left) and distorted mesh (right). Results for element H1/E9 (top) and H1/E9T (bottom)

F I G U R E 13
Second eigenmode of element QA1/E4T-F −T 0 for regular and distorted mesh (Spurious behavior especially pronounced in second mode shape) Some of the considered elements are not able to exhibit any instability. This concerns displacement element H1 in both distorted and regular meshes and determinant enhanced element H1/FJ-TnB in regular meshes. All other elements show four instabilities and can be grouped into two major categories: Elements that show physical eigenmodes and ones that exhibit hourglassing. Typical behavior of elements in those categories is shown in Figure 12 for the standard EAS elements H1/E9 and H1/E9T. As suggested by the investigations shown in Section 6.4, the EAS element based on Wilson-modes exhibits severe hourglassing, while the element based on transposed Wilson-modes only reveals physically correct eigenmodes. Note that in regular meshes, spurious behavior of H1/E9 appears at a state of u i = 0.39 (corresponding to 1 = 0.61) as predicted in Section 6.4. In distorted meshes spurious behavior of H1/E9 is restricted to smaller parts of the domain and appears at higher displacement levels u i = 0.51. This is because distortion slightly "stabilizes" H1/E9.
The first physical mode exhibited by H1/E9T appears at u i = 0.62 in the distorted and u i = 0.46 in the regular case, respectively. These results converge towards u i = 0.46 in finer meshes, which indicates the deficient performance of H1/E9T in 3D distorted meshes. In contrast to that, the present discretization level is sufficient in 2D to get almost identical results with both types of meshes.
All other elements (with the exception of HA1/E12T-F −T 0 ) give qualitatively the same physically correct results as H1/E9T. However, they exhibit slightly different levels at which the instabilities occur. Superior behavior is exhibited by H1/FJ-TaB which yields almost identical results u i = 0.45 and u i = 0.44 in the distorted and regular case, showing its insensitivity to mesh distortion.
A special case are elements QA1/E4T-F −T 0 and HA1/E12T-F −T 0 , which use a special transformation to ensure objectivity of the enhanced part of the deformation gradient as proposed by Pfefferkorn and Betsch. 19 While these elements show excellent behavior with regular meshes as shown in, 19 they exhibit a spurious instability in the distorted case identified by very locally confined spurious eigenmodes. This behavior is depicted in Figure 13. Because of this, we do not recommend to use these elements in generally distorted meshes.

Mesh distortion
This widely used 2,8,15,19,35,40,59,[61][62][63] benchmark tests the element's behavior in distorted meshes. It is performed on a cantilever like structure depicted in Figure 14 with length L = 10 and cross-sectional dimension of h = 2 and t = 1. The  Figure 15 where the normalized top edge displacement is plotted against the degree of skew s. Due to the material and geometric nonlinearity, no analytic solution of the problem at hand can be found. Thus, normalization is conducted with a converged result obtained with 40 × 8 Q1/E4 elements.
The first result to be drawn from Figure 15 is the severe locking of the displacement element Q1, which only yields about 30% of the required deformation in the undistorted (s = 0) case and exhibits even worse results in distorted meshes. This outcome is only slightly improved by the novel elements Q1/E4-S001, Q1/FH 1 -TaXx and Q1/FH 2 -TnTd. The best results are obtained with elements Q1/E4T and Q1/FJ-TaB. Those elements are able to give almost the exact displacement for undistorted meshes with s = 0 and are generally relatively insensitive to mesh distortion ‖‖‖ . Element QA1/E4T−F −T 0 exhibits slightly more distortion sensitivity due to the additional Gauß-point. 19 Finally, element Q1/FH 3 -TnVv performs well in the mesh distortion test and even improves upon some of the results shown by the standard H1/E9T element. However, it exhibits too high displacements for skew s = 0 when compared with reference solution. This means that it does not give the correct result for only a few elements, but will still converge to the analytic solution for finer meshes (see Section 6.7). The same problem occurs with elements proposed, for example, in publications. 10,15

Cook's membrane
The final test for elastic problems is the well-known Cook's membrane test. 2,3,9,10,13,14,19,23,40,44,45 It is used to examine shear and volumetric locking in distorted meshes, coarse mesh accuracy, convergence behavior with mesh refinement (h-convergence) and robustness of the elements. Both, the elastic and elasto-plastic material described in Appendix B are considered, with material parameters a = 126, b = 252 and c = 81661 (corresponding to E = 2261 and = 0.4955 in ‖‖‖ Note that there is a limit to improving the elements behavior in distorted meshes as shown by MacNeal. 24 linear theory) chosen for the Mooney-Rivlin material (B1) and the standard parameters (Table 3) used in the elasto-plastic simulations.
The trapezoidal geometry and the boundary conditions of the test are shown in Figure 16. The structure is clamped on the left side and load elast = 100 or plast = 0.26 is applied as constant shear stress on the right side in the elastic and plastic case, respectively. This load is applied in n steps until the maximum load is reached. While a constant number n steps = 10 is chosen in elasto-plastic simulations, we determine the smallest number n steps necessary for the Newton-Raphson scheme to converge in elastic simulation in order to examine the robustness of the elements.
To complete the setup of this test, a mesh with two elements in direction of the thickness and various numbers n el = {2, 4, 8, 16} of elements per side is considered to investigate convergence of the displacement u (see Figure 16) with mesh refinement. The mesh and deformed configurations are shown in Figure 16.
In the first evaluation step of this test we examine the displacements u = u y of the top right corner for various numbers of elements. The results are shown in Figure 17 for both the elastic and elasto-plastic simulation. In the elastic case severe locking of the displacement element H1 and relatively good performance of the standard EAS elements H1/E9 and H1/E9T can be observed. In the group of novel elements, selective enhancement element H1/E9-S001 performs already a lot better than the displacement element, even though only volumetric locking is affected by solely enhancing J. However, it does not improve upon the standard EAS elements. The same is valid for the cofactor enhanced elements H1/FH 1 -TaXx and H1/FH 2 -TnTd. Of the novel elements, H1/FJ-TaB, H1/FH 3 -TnVv and H1/FJ-TnB show the best performance and even outperform the standard EAS elements, which is especially interesting for the latter two, since they both use only nine additional DOFs.
In the elasto-plastic simulation, similar results to those of the elastic case considered above are obtained for most elements. An interesting result is the performance of H1/FH 3 -TnVv and H1/FJ-TnB, which is now worse than the standard EAS element's. This happens due to locking induced by the volume preserving nature of the von Mises yield condition. Mainly volumetric parts of J (and not F) are enhanced by the enhancement structure of the affected elements. This type of enhancement has unfortunately no influence on the constraint and ultimately leads to the observed behavior.
After having examined h-convergence, we investigate the robustness of the proposed elements in a second step. It is determined by identifying the required number of load steps n steps for the Newton-Raphson scheme to converge. For H1/E9 and H1/E9T it amounts to 6 and 7 steps, respectively. All other elements, which exhibit equal or better performance F I G U R E 18 Geometry, mesh and boundary conditions for the plane strain necking test (rotated 90 • )

F I G U R E 19
Results of the plane strain necking test. Necking zone depicted for various elements and u = 5.6. Colors show the accumulated plastic strains where blue denotes a low and yellow a high value in the convergence analysis, need the same amount of steps. Only elements showing poor performance in the analysis above need fewer load steps since their final displacements are smaller. Note, however, that the EAS elements are not known to be very robust. 12

Necking plane strain
The final two examples concern necking problems in elasto-plasticity. First, we consider a rectangular bar with a length of L = 53.334 and width of B = 12.826 subject to plane strain. 3,9,10,15,28,31,60,62 A small geometric imperfection, which is needed to initiate necking, is introduced by linearly reducing the width from the ends to the middle of the specimen by ΔB = 0.14. Due to symmetry only one quarter of the specimen has to be considered, which is shown in Figure 18 together with the boundary conditions and mesh. Symmetry conditions apply on the left and lower boundary. The structure is loaded by prescribing displacements u = 7, which are applied in 200 load steps **** , on the ends of the specimen. Furthermore, we employ the line-search algorithm given in Bonet and Wood 56 to stabilize the Newton-Raphson scheme. A mesh with 200 elements, where the lower fifth of the bar is refined, completes the test's setup.
The main issue to be addressed with this test is hourglassing in tension already mentioned in the one element stability test in Section 6.4. 9,10 As described there, hourglassing occurs due to a material instability of the elasto-plastic model, which unfortunately affects the hourglass modes. Figure 19 shows the necking zone of some elements for a displacement state with u = 5.6. The elements can be grouped in three sets. Group one contains Q1 and Q1/FJ-TnB which are not able to display the necking behavior as consequence of locking. Determinant enhanced element Q1/FJ-TaB, and the three standard EAS elements form the second group characterized by hourglassing patterns. Note that QA1/E4-F −T 0 exhibits ****Note that this many steps are necessary in order to capture the hourglassing behavior described below. Elements not prone to hourglassing would need far less steps.

Necking of a circular bar
The final test in this work is a necking example of a circular bar 3,4,15,28,32,34,35,62 which is similar to the plane strain necking test presented in the previous section. It is conducted on a circular bar of length L = 53.334 and Radius R = 6.413 of which only one eighth is considered due to symmetry. The mesh with 960 eight node elements is shown in Figure 20. Boundary conditions and the imperfection to initialize necking are applied as described in Section 6.8. Since no hourglassing occurs in this benchmark, less load steps are sufficient. We use 15 load steps until u = 5.6 is reached and another 15 up to u = 7.0 as that range is more demanding. Figure 21 shows the necking zone of the circular bar for various elements at u = 7.0. Again, the elements can be put into two groups: Elements that show locking, which is revealed by large remaining cross-sections, and elements that exhibit only little or no locking at all. This is easily determined by the radii R P and R Q of points P and Q (see Figure 20) in the deformed state at u = 7.0. 15 In the first group, the mean ratio of radii is always above 0.55, which indicates necking of little account. Within this group, slightly better results are exhibited by H1/E9-S001, H1/FH 2 -TnTd and H1/FH 3 -TnVv, whereas elements H1, H1/FH 1 -TaXx and H1/FJ-TnB show almost no necking. The other elements are not prone to locking, which can be established from q R ≤ 0.40. Of the newly proposed elements, H1/FJ-TaB performs best in this simulation. Furthermore, it is interesting to investigate the ratio R Q ∕R P . This reveals a mesh dependency of H1/E9 and H1/E9T since these elements have R P ∕R Q ≈ 1.06. The other two well-working elements exhibit a much smaller difference of radii showing that they are capable of dealing with the distorted mesh of the present example. Thus, it is possible to reduce mesh dependency with the determinant enhancement proposed for H1/FJ-TaB which is inline with the results from Section 6.5.

CONCLUSION
Three novel classes of enhanced assumed strain elements have been introduced in the present work. All of them rely on a framework for polyconvex strain energy functions presented by Bonet et al 41,43 and the recently reintroduced tensor cross product. A key concept of the framework for polyconvexity is that three kinematic measures, namely the deformation gradient its cofactor and determinant, are principal variables. This approach facilitates a plethora of possibilities for novel enhancement strategies in the context of enhanced assumed strain elements. The three versions • selective enhancement of the kinematic fields (Section 3), • cofactor enhancement (Section 4) and • determinant enhancement (Section 5) are proposed in the present work. All newly proposed element types are based on standard principles and meet important requirements for finite elements. First, Hu-Washizu type variational functionals are presented in Sections 3-5 for all three classes, which ensures accordance with (displacement based) continuum mechanics. Second, it is shown both numerically and analytically that all of the elements satisfy the crucial patch test condition, provided that simple requirements on the shape functions of the enhanced fields are fulfilled. Furthermore, all formulations are constructed such that they are frame-invariant regardless of chosen ansatz functions. Another advantage of the elements presented herein is the strain driven format inherited from EAS elements, which allows straightforward implementation of complex material laws † † † † . The elements are also efficient due to the L 2 -orthogonality imposed on the discrete enhanced strain and stress fields as well as the possibility to statically condense the internal DOFs. Finally, all novel elements are free of spurious hourglassing in elastic simulations. Only in elasto-plastic examples some of the newly proposed elements show hourglassing while others are prone to mild locking. While problems arise in plasticity as described above, the elements exhibit very promising results in the extensive elastic numerical studies conducted in this work. They perform equally or even better than standard EAS elements in all simulations and are very well suited for the type of problems considered. The overall best performance is exhibited by elements H1/FJ-TaB and H1/FH 3 -TnVv. The first one is completely locking free which is especially advantageous for the portrayal of physical instabilities (see Section 6.5) and computation of limit loads in elasto-plasticity. 64 The latter element is able to overcome the spurious hourglassing of EAS elements in elasto-plastic simulations (see Section 6.8). Unfortunately, mild locking occurs for that element in those simulations. Besides the presentation of new elements, we propose a novel eigenvalue analysis for 3D eight node elements to examine spurious hourglassing. This method extends the 2D procedure introduced by Glaser and Armero 9 and gives further insight into the elements' behavior.
In future works other shape functions and further approaches of enhancement such as simultaneously enhancing all three fields could be considered. Another line of research could follow combination of the newly proposed procedure with other mixed elements. Especially the combination with Q1/P0, 65 which has been done for EAS elements by Armero, 10 seems promising. Furthermore, an extension of the methods to triangular and tetrahedral meshes would be desirable to increase applicability of the elements in engineering problems. This extension has been made for mixed EAS elements in Caylak and Mahnken 20,22 and could be adapted for the present elements as well. However, we believe that the most important issue to be addressed is the behavior in plasticity, where the novel elements show either hourglassing or locking. This is not only an issue of the elements proposed herein but also the key issue of the highly popular EAS elements. Resolving this issue remains a major challenge in element technology.
Therein a, b and c denote the three independent material parameters and d = 2a + 4b to ensure a stress free reference configuration. Strain-energy function (B1) is polyconvex if a, b, c > 0 holds. 41 In the current framework (B1) is especially interesting, since there is no coupling between the kinematic fields. This ultimately leads to simple expressions for the second derivatives of (B1) which are needed for the Newton-Raphson solution scheme. 41 Remark 11. Note that (B1) can easily be transformed to use the symmetric kinematic fields (7), which yields the equivalent form W(C, G, C) = a (tr(C) − 3) + b (tr(G) − 3) + c 2 . (B2)

B.2 Elasto-plastic material model
The second material model considered in the present work is the elasto-plastic model proposed by Simo 32 which is widely used in the context of finite element development 4,10,12 and is an eigenvalue based formulation. Its elastic response is governed by a Hencky strain-energy function which employs the logarithmic principal stretches. The plastic part of the model is governed by the von Mises yield condition with isotropic nonlinear hardening and the associative flow rule. More information on the material model and algorithms for standard elements are given in the work of Simo. 32 In the present appendix the model is presented in-depth for the polyconvexity inspired framework. We start by assuming the multiplicative splits of the kinematic fields into plastic and elastic parts. Therein, (B3a) is the by now standard 66 multiplicative split of the deformation gradient and the other two relations follow from the first relation if compatible kinematic fields (ie, Equation (6)) are assumed. In our general framework, however, (B3c) and (B3b)  and C e = J 2 det(C p−1 ) is the elastic determinant ‡ ‡ ‡ ‡ . Moreover, history variable C p−1 , which is the inverse plastic right Cauchy-Green tensor, is introduced. Note, that the elasto-plastic model only uses symmetric fields in order to enable the compact eigenvalue based formulation. The plastic part of the model is governed by the von Mises yield condition with isotropic nonlinear hardening given by where H, Y 0 , Y∞ and denote material constants and is the history variable for hardening representing the accumulated plastic strains. Furthermore, , a = 1, 2, 3 ( B 8 ) denote the principal stresses of the Kirchhoff stress tensor = ∑ 3 a=1 a n b,e a ⊗ n b,e a . This stress tensor is in general not compatible with the constitutive stress-like quantities presented throughout this work, which may lead to unsymmetric Kirchhoff stresses as pointed out by Bonet et al. 41 However, symmetry of the stress tensor used for the yield condition is vital in order to maintain a simple implementation of the elasto-plastic model. For this reason, the formula above is employed for the yield condition, while the residuals are computed with the constitutive stresses (21), (57) and (77). The ‡ ‡ ‡ ‡ Note that C e must not be computed from det(b e ) since it is in general not equal to J 2 det(C p−1 ) due to separate treatment of the enhancement of the kinematic fields. This is a direct consequence of (B3). plastic response is completed by assumption of a associative flow rule where is the plastic multiplier.
For the computer implementation of the model we use the radial return map presented by Simo, 32 which employs the exponential tensor map and exactly preserves the volume of the plastic deformation. For the given material model, it is especially simple since it leads to a single scalar equation to determine the incremental plastic multiplier Δ .
The stresses and algorithmic tangent we obtain from the descriptions above are based on symmetric measures and b e . Thus, they have to be transformed to account for the unsymmetric fields P and F we use for the element formulations throughout this work. This can be done in a straightforward manner.

APPENDIX C. EIGENVALUE ANALYSIS OF EIGHT NODE 3D ELEMENT
This appendix contains an analytic eigenvalue analysis for the 3D uniaxial stress case necessary for the stability analysis of 3D elements, which is presented in Section 6.4. It extends the 2D analysis first proposed by Glaser and Armero 9,10 which is briefly summarized in Section 6.4. The 3D version has, however, to the best knowledge of the authors never been proposed before.
Note that the analysis presented here has in analogy to the 2D analysis introduced by Glaser and Armero 9 only a very limited validity: It only covers the simple uniaxial stress case for an undistorted single element. Nevertheless, many conclusions can be drawn from the eigenvalues computed in this way.
The eigenmodes of interest here are only those which yield non-constant strains since the ones with constant states of strain are always exactly reproduced as the patch test is fulfilled by all elements presented throughout this work. The modes with non-constant strains can be identified by looking at an alternative representation of the trilinear Lagrangian shape functions N I defined in (23), which is frequently used in the context of hourglass stabilization. This vectorial form is given by Note that all eight vectors a I and h A are orthogonal and therefore a basis of R 8 . Since strains are computed from the derivatives of the shape functions (see (29)) and the first two terms of (C1) are linear, we conclude that the subspace spanned by the four hourglass vectors h A includes all modes with non-constant strains of the eight node element in its reference state. In a next step, the hourglass vectors are arranged in matrix