Mesh distortion insensitive and locking-free Petrov–Galerkin low-order EAS elements for linear elasticity

One of the most successful mixed finite element methods in solid mechanics is the enhanced assumed strain (EAS) method developed by Simo and Rifai in 1990. However, one major drawback of EAS elements is the highly mesh dependent accuracy. In fact, it can be shown that not only EAS elements, but every finite element with a symmetric stiffness matrix must either fail the patch test or be sensitive to mesh distortion in bending problems (higher order displacement modes) if the shape of the element is arbitrary. This theorem was established by MacNeal in 1992. In the present work we propose a novel Petrov–Galerkin approach for the EAS method, which is equivalent to the standard EAS method in case of regular meshes. However, in case of distorted meshes, it allows to overcome the mesh-distortion sensitivity without loosing other advantages of the EAS method. Three design conditions established in this work facilitate the construction of the element which does not only fulfill the patch test but is also exact in many bending problems regardless of mesh distortion and has an exception-ally high coarse mesh accuracy. Consequently, high quality demands on mesh topology might be relaxed.

isogeometric analysis, see Reference 3). An explanation for this can be found in a landmark paper published by MacNeal 11 in 1992 (see also the preliminary work 12 ), which has in the opinion of the authors not gotten the attention it deserves. In this work MacNeal proves the theorem that a finite element with arbitrary shape cannot simultaneously satisfy the constant strain patch test 2 and be exact for higher order modes under the premise of a symmetric stiffness matrix. This fundamental limit to an element's perfectibility applies regardless of internal methodologies, that is, which of the remedies to cure locking is used. Thus, the final conclusion drawn by MacNeal 11 is that it is impossible to substantially improve elements beyond what has already been achieved with the standard (Bubnov-Galerkin) approach. Indeed, to the best knowledge of the authors, no element could so far break the limits of MacNeal's theorem * .
In particular, most elements have symmetric stiffness matrices, which come of course with many advantages † . Moreover, they are usually designed to fulfill the patch test, see for example, References 8,9 since this comes with a lot of benefits 2 as well. Consequently, in accordance with MacNeal's theorem, they perform poorly in (higher order) bending problems if meshes are distorted. Vice versa, the less frequent choice of (deliberately or not) sacrificing the ability to pass the patch test enables construction of elements with increased bending accuracy. An example for this approach is the incompatible mode model by Wilson et al. 16 which famously fails the patch test, and the recently proposed method of reverse adjustment to the patch test by Hu et al. 17,18 (which ironically does not pass the strong patch test).
All in all, improvement of element performance beyond what has already been achieved almost inevitably leads to an unsymmetric stiffness matrix. MacNeal even briefly mentions this possibility but then deems this option "abhorrent for many reasons". 11 However, the works of for example, Rajendran et al., 19,20 Xie et al., 21 and the present contribution show that there is much to be gained with unsymmetric stiffness matrices.
The first element with unsymmetric stiffness matrix has been proposed by Rajendran and Liew, 19 who chose a Petrov-Galerkin approach for higher order displacement elements instead of the usual Bubnov-Galerkin ansatz, and named their method unsymmetric finite element method. The key idea is to use metric 1,22 shape functions, which are constructed in the physical space, as ansatz for the trial function of the displacement while the usual isoparametric functions are employed for the test function. Unfortunately, merely using physical coordinates leads to frame-dependent ansatz functions as noted by Ooi et al. 23 A cure for this issue, which does not induce other problems such as anisotropies, has been proposed by Xie et al. 21 It involves skew coordinates, which are affine-equivalent 24 to physical coordinates and objective. The skew frame has first been proposed for assumed stress finite elements (cf. Yuan et al. 25 and Wisniewski and Turska 26 ) and has subsequently also been employed for other mixed elements [26][27][28][29][30] and unsymmetric finite elements 21,31,32 .
While higher order unsymmetric finite elements work well and pass higher order patch tests regardless of element shape (cf. the works [19][20][21][33][34][35], there is a fundamental issue for low-order elements: Consider for example, a plane four-node quadrilateral element with eight displacement degrees of freedom. It is then impossible to represent complete quadratic polynomials since these would require at least 12 degrees of freedom. Thus, it is crucial to choose the higher order modes carefully in order to get the best performance. Consequently, the unsymmetric low-order elements presented so far have either extremely complex ansatz functions which are material dependent 13,21,32,36,37 or require higher order integration and many internal degrees of freedom. 31 Furthermore, the elements are often not straightforward to extend to the 3D case. However, the advantages of low-order elements with regard to mesh generation, (stress) singularities and bandwidth (sparsity) of the stiffness matrix make low-order unsymmetric finite elements desirable. Hence, we propose a novel low-order unsymmetric mixed element based on a Petrov-Galerkin approach for the enhanced assumed strain (EAS) method introduced by Simo and Rifai. 9 We will show that this approach allows to construct low-order unsymmetric elements without the drawbacks described above.
We choose to start our developments from the EAS method since the commonly used symmetric (Bubnov-Galerkin) version of it features several desirable properties: it fulfills the patch test, is locking-free in case of undistorted meshes, and, most importantly, it is straightforward to extended to large deformations and general material * Cen et al. 13 claim to break through MacNeal's theorem. However, their element has an unsymmetric stiffness matrix and does therefore not violate the limits proven by MacNeal. 11 Another interesting candidate is the element by Wu and Chueng. 14 It has a symmetric stiffness matrix and is exact in a bending problems at some nodes. Unfortunately, the displacement is not exact at all nodes and it does thus also not break through MacNeal's theorem (see the work of Sze 15 ). † For example, the reduced computational cost to solve the linear equation system and decreased memory consumption since only half of the sparse matrix has to be stored. models due to its strain-driven format (cf. Simo et al. 38,39 and Glaser & Armero 40 ). Consequently, it is one of the most frequently used mixed methods in research and application (cf. e.g., the works 9,27,28,[38][39][40][41][42][43][44][45][46][47][48][49]. With the newly proposed Petrov-Galerkin approach we show that it is possible to overcome the sensitivity to mesh distortion of existing EAS elements. To that end, we first establish three design conditions required in order for the element to be exact for a specific displacement mode. If these conditions are met, the element is exact for that mode in the sense that the nodal displacements coincide with the analytic solution (nodally exact response). In particular, we choose to fit the novel element to a modified version of the assumed stress modes proposed by Pian and Sumihara 8 and Pian and Tong 50 in 2D and 3D, respectively. These stress modes include besides the important patch test states also bending modes which are of the utmost importance in many engineering applications. The resulting element is therefore exact for constant stress and bending problems regardless of element shape. It is also locking-free, frame-indifferent, isotropic and its stiffness matrix is integrated exactly by the standard Gauss-quadrature rule. Moreover, the newly proposed unsymmetric EAS element exhibits an substantially increased coarse mesh accuracy. Finally, in case of regular meshes, the element coincides with the original well working EAS element by Simo and Rifai. 9 The design of the novel element with all these desirable properties is made possible by combining ingredients from a multitude of previously developed element formulations. Besides the obvious EAS framework 9,39 we also employ ideas from assumed stress approaches, 8,26,50 the skew coordinate frame, 21,25,26 incompatible mode elements, 7,16,31,51 and others. 5,29,52 Naturally, the unsymmetric incompatible mode 2D element proposed by Huang et al. 31 is closely related to the present Petrov-Galerkin EAS element. However, the novel EAS approach is more general, requires less internal degrees of freedom, and allows to examine the underlying mechanisms in a deeper way. Interestingly, a violation of one of the three design conditions for exact solutions mentioned above is the reason why the 2D element by Huang et al. 31 is not straightforward to extend to 3D.
The present work is structured into six sections. In Section 2 we revisit MacNeal's theorem since it is key to the methods developed in the remainder of this work. We aim at giving a simpler approach to the proof in Reference 11 and present some extensions to the original theorem of MacNeal. To that end we compare the finite element approximation to the continuum description and examine both formulations in detail in Sections 2.1 and 2.2, respectively. After that, conclusions of MacNeal's theorem are drawn in Section 2.3. Some generalizations of the proof in Section 2 can be found in Appendix A and the key findings of MacNeal's theorem are summarized in Section 2.4. Section 3 covers the weak form for EAS elements and the three design conditions required for nodally exact solutions. Afterwards, we determine ansatz spaces which fulfill the design conditions in Section 4. To that end we first introduce the skew coordinate frame in Section 4.1 and describe the analytic modes for which the element is optimized in Section 4.2. The actual ansatz spaces for the Petrov-Galerkin EAS elements are covered in Section 4.3 for the 2D case and in Section 4.4 for the 3D case, respectively. Numerical simulations comparing the novel element to established ones are presented in Section 5 before conclusions are drawn in Section 6.

MACNEAL'S THEOREM
MacNeal's theorem 11 states that a finite element cannot simultaneously pass the patch test (see e.g., Reference 2) and be exact for higher order displacement modes in case of arbitrarily distorted meshes if its stiffness matrix is symmetric.
To prove this proposition we consider, following MacNeal, 11 a linear elastic continuum  and a single finite element Ω e ⊂  of arbitrary shape and arbitrary number of nodes N which is embedded into . The two domains interact at the element's boundary Ω e . This connection is described by displacement u(x) and traction t(x) ‡ for x ∈ Ω e , where x are Cartesian (or physical) coordinates. Comparing the continuum solutions u * and t * with the finite element aliases u h,e and t h,e on Ω e allows to evaluate how well the element replaces the respective continuum domain. § To that end we examine both formulations in depth in the sequel.
Example Q1. Figure 1 shows exemplarily the case of a plane quadrilateral element with straight edges and four nodes (Q1). This simple case is used throughout the remainder of this work to clarify relations. ‡ In order to improve readability, arguments of functions are frequently omitted in the sequel. § Throughout the remainder of the work we denote continuum solutions with an asterisk (•) * and finite element approximations with superscript (•) h . Given a prescribed displacement field u * (x) for x ∈ , the usual relations for a linear elastic continuum * ∶= ∇ s determine the corresponding strain * , stress * , and traction t * . Therein, C denotes the symmetric positive definite fourth order elasticity tensor, ∇ s x (•) symbolizes the symmetric part of the gradient and n is the unit outward normal on a surface. Furthermore, the strong form of equilibrium is given by where b * is the field of body force which is readily determined by the prescribed displacement field u * via (1) and (2). By virtue of the given displacement u * and (1), it is also straightforward to compute appropriate Dirichlet boundary conditions u * = u * on  u ⊂  along with appropriate Neumann boundary conditions t * = t * on  t ⊂  on the body's boundary . Here, the standard conditions  u ∩  t = ∅ and  u ∪  t =  hold. The corresponding weak form (or principle of virtual work) is given by where is an arbitrary test function. Moreover, in the sequel the left and right hand side of (3) are abbreviated by G * int and G * ext and identify the internal and external part of the weak form, respectively.

Displacement modes on Ω e
where the components K * mn constitute the modal stiffness matrix. The major symmetry of C implies the symmetry of K * mn . Since the weak form (3) holds for every subdomain of , relation needs to be satisfied. Now the arbitrariness of m yields K * mn = F * mn . In particular, this establishes the symmetry of F * mn . Finally, by use of Gram-Schmidt orthogonalization any set of displacement modes u * m can be uncoupled such that matrix F * mn can be made diagonal. In the following we assume that the modes u * m are uncoupled and thus F * mn is diagonal. A final conclusion drawn from (8) is that f * m = 0 for rigid body modes u * r m since * r m = ∇ s x u * r m = 0 holds.

Connection between nodes and modes
In the next step we consider the finite element approximation. In contrast to the continuum description, the finite element approximation is based on nodes rather than modes. It is therefore crucial to determine the relation between the node based finite element displacement u h,e and the displacement modes introduced in (4). For standard (interpolatory) finite elements the displacement on element level is given by where N is the number of nodes and u e i is the displacement at node i. The corresponding scalar ansatz functions ¶ M e i , which may (for now) differ from element to element, have the Kronecker-delta property M e i (x j ) = ij , where x j is the position of node j. With these definitions at hand, the relation between modes and nodes can be described by where is the value of mode u * m at node x i . This ensures that the finite element interpolation u h,e (11) coincides with the prescribed displacement field u * at the nodes, that is, u h,e (x j ) = u * (x j ). However, at other points x ∈ Ω e the values do generally not coincide.
Remark 2. Assuming a displacement field of the form (12) is rather restrictive. In general, M e i are not necessarily scalar and the weights u e i need not correspond to the nodal displacements. Such non-interpolatory approaches include the widely used hierarchical higher order elements 2 and the isogeometric analysis. 3 These generalizations are covered in Appendix A.3.

2.2.2
General finite element framework Regardless of internal element methodologies, the relation between nodal displacements u e i and nodal forces P h,e i of any linear finite element can be written as where K h,e ij denotes a partition of the element's nodal stiffness matrix K h,e . This format is very general and includes displacement-based elements, reduced integrated elements as well as mixed methods with internal degrees of freedom # . The left-hand side of (14) refers to the internal part associated with elastic deformations while the right-hand side contains the external nodal forces. Defining the element's modal stiffness matrix and modal external forces bỹ allows to recast (14) by virtue of (12), pre-multiplication with (u * ,e im ) T and summation over i in the form The last equation allows to determine the maximum number of modes M for which an element can exactly reproduce the continuum response in the sense that the nodal displacements are exact. In order for (16) to be solvableK h,e mn must have proper rank. To that end we assume that our finite element fulfills the design imperative of stress-free response to ¶ In anticipation of the Petrov-Galerkin approach proposed in the sequel we denote the nodal shape functions in (11) by M e i instead of the more common N e i . # Internal degrees of freedom can be statically condensed on element level to achieve form (14). rigid body motions. For the corresponding M r rigid body modes f * m = 0 as shown in Section 2.1.2 while for all other modes f * m ≠ 0.
For the finite element to correctly represent rigid body motions we thus require accordingly f h,e m = 0 and rank[K h,e mn ] = M − M r . Naturally, the nodal stiffness matrix has to correctly account for rigid body motion as well. Its rank is then, provided that there are no spurious instable modes, N DOF − M r where N DOF = N d is the number of degrees of freedom and d is the spatial dimension. Rewriting the element's modal stiffness matrix (15a) in vector-matrix notation yields K h,e = (u * ,e ) T K h,e u * ,e . Thus, the element's modal stiffness matrix has at best the same rank as the nodal stiffness matrix || .
Therefore, a finite element can in general exactly represent no more than M = N DOF modes including M r rigid body modes.
Example Q1. In case of the quadrilateral element there are M r = 3 rigid body modes, N = 4 nodes, and N DOF = 8 degrees of freedom which determines M = 8. Thus, considering the three rigid body modes, the element can only be adapted to a maximum of five displacement modes. Thereof, three are usually chosen to be the constant stress modes (patch test, represented by linear displacement modes). The two remaining higher order modes should therefore be chosen carefully to get the best performance (cf. Cen et al. 13 ).
Remark 4. The total of eight modes for Q1 is especially not sufficient to represent fully quadratic displacement modes, which require a total of 12 modes (1, x, y, xy, x 2 , y 2 in both displacement components). It is therefore futile to try and improve quadrilateral elements to exactly represent fully quadratic displacement fields regardless of the applied technique. This observation is the essence of the principle of limitation by Fraeijs de Veubeke 54 .
So far, we focused on the internal part of the weak form leading to the stiffness matrix. We now turn to the external part of the weak form which yields the nodal forces P h,e i (14). In particular, the nodal forces are related to the traction t h,e acting on the element's boundary Ω e . We assume that the analytic traction acts on the finite element, that is, t h,e = t * . For the element's response to be nodally exact we ultimately have to show that assumption t h,e = t * is consistent with u e i = u * (x j ) (see Section 2.3). For the finite element discretization of the external part of the weak form (9) we get where the approximation of the test functions assumes the form (11) with nodal shape functions now denoted by N e i . Inserting this result into definition (15b) of the element's external modal forces f h,e m yields The modal forces f h,e m can be viewed as discrete counterpart of the modal forces f * m introduced in (7) and establish the connection between the continuum and the discrete finite element formulation. Comparing (18) and (7) motivates the introduction of In this way, the components F h,e mn have the same structure as the components F * mn defined in (7). In particular, u h,e m introduced in (19) can be regarded as the finite element alias of the displacement mode u * m . We emphasize again that in general both quantities coincide only at the nodes as pointed out in Section 2.2.1.

2.2.3
Nodal equilibrium All investigations so far focused on a single finite element Ω e . However, the correct reproduction of displacement modes u * m on element level is not sufficient to ensure exact solutions for patches of elements. We show subsequently, that it is crucial to additionally fulfill equilibrium of nodal forces. || Theorem rank(AB) ≤ min(rankA, rankB) holds for arbitrary matrices A and B (see e.g., the work 53 ). F I G U R E 2 Two finite elements Ω a and Ω b with shared boundary Ω ab separated by Δx (left) and a patch of finite elements embedded into the linear elastic continuum (right) To that end we consider two finite elements Ω a and Ω b embedded into the linear elastic continuum  as shown in Figure 2 on the left. Assume that their adjacent faces Ω ab and Ω ab are congruent and separated by Δx. In case of an infinitesimal Δx clearly Ω ab → Ω ab and t * → −t * .
Furthermore, the total contribution from boundary Ω ab to the nodal force at nodeζ → i is given by where use has been made of (17). For nodal equilibrium we require P h i,ab = 0, which ensures that, after assembly, P h i = 0 at every node i inside a finite element patch (see Figure 2 on the right). Moreover, contributions from internal faces to nodes at the boundary of an element patch are likewise zero (see e.g., node j in Figure 2 on the right). Thus, if P h i,ab = 0, all interior contributions vanish and it is possible to consider a complete patch of elements as a single element (with many nodes). Conversely, if nodal equilibrium is fulfilled, it is justified to investigate only a single finite element since results can then be transferred to arbitrary meshes.
To fulfill condition P h i,ab = 0 for an arbitrary traction t * , it is obvious in view of (20) that a valid choice is given by With this choice, nodal equilibrium is not only fulfilled for the considered M displacement modes u * m but for arbitrary traction t * . Condition (21) is the classical conformability requirement or, likewise, the C 0 inter-element continuity requirement and severely restricts the choice of ansatz functions N i . Suitable functions can be created by mapping parent elements from the parametric space to the physical space ** . This approach was indeed introduced to fulfill nodal equilibrium and mapped elements are so successful that they have been used almost exclusively since their first occurrence in the 1960s. All in all, (21) is a usual assumption and not as serious of a restriction in practical applications.
Even less restrictive choices to fulfill nodal equilibrium than (21) still put a restriction on the design of N e i . This ultimately establishes MacNeal's theorem as shown in the next section.
Example Q1. To demonstrate the restrictions imposed by nodal equilibrium we consider again a quadrilateral element: The only sensible choice of parent element for a quadrilateral is a square. All other forms would induce anisotropies into the formulation. Furthermore, the four linearly independent monomials in the parametric space have to be 1, , , since any other choice would either be incomplete or anisotropic. These restrictions lead directly to the well-known bilinear Lagrangian ansatz functions.
Thus, there seems to be only one reasonable choice for the ansatz functions of the displacement test function if conformability is to be fulfilled. This illustrates the severe restriction imposed by (21).

Design problem
Having summarized properties of the continuum and the finite element formulation allows to finally determine the optimal stiffness matrix for which the finite element is capable to exactly represent M = N DOF displacement modes.  (19). Now, the right-hand side of (22) can be recast in the form where only F h,e,LL mn is diagonal and equal to the corresponding part of F * mn . Equation 23 clearly shows that F h,e mn is unsymmetric which is induced by the coupling between lower and higher order modes.
Accordingly, satisfaction of design problem (22) would require the modal stiffness matrixK h,e mn to be unsymmetric. This in turn would mean an unsymmetric nodal stiffness matrix † † as can be concluded from (15a). Conversely, a symmetric stiffness matrix makes it impossible to simultaneously fulfill the patch test and correctly represent higher-order displacement modes, which is the central statement of MacNeal's theorem. It is worth mentioning that MacNeal's theorem does not rule out special cases such as in the examples below.
Example Q1. A quadrilateral element can simultaneously satisfy the patch test and exactly represent quadratic displacement modes under certain conditions. An example for this is the EAS element proposed by Simo and Rifai 9 . This element is capable to exactly capture pure bending of a cantilever beam (see sec. 6.1.1 in Reference 9 and the 3d version of this example in Section 5.4). However, the exact response is limited to the special case of rectangular elements. In that case the integrals in (7) and (19) have the same value and thus F h,e mn = F * mn holds even though u h,e m ≠ u * m . For other element shapes the results deteriorate rapidly.
Example Q2. Another exemption can be observed for quadratic displacement-based nine-node quadrilaterals (Q2) (cf. MacNeal, 1 Ch. 3.5 and Zienkiewicz et al., 2 Ch. 5.7). These elements are capable of representing element shapes with quadratic borders. However, if we restrict the shape to bilinear forms (i.e., the geometric capabilities of quadrilateral elements), the elements are capable of correctly representing full quadratic polynomials in the physical space. Thus, for bilinear shapes u h,e m = u * m holds for all quadratic modes and the element is exact for bending problems.
Higher order modes (cubic and quartic) incorporated in the base of Q2 are unfortunately not represented correctly. Furthermore, it is not advisable to generally restrict the geometry to bilinear shapes in case of curved domains since this degrades the order of convergence with mesh refinement (cf. Braess,55 Ch. II. §1). † † The option of an unsymmetric finite element formulation was already recognized by MacNeal 11 . At that time, however, it seemed "abhorrent for many reasons" (MacNeal 11 ) to consider such elements. However, previous works 13,19,31 and the present contribution show, that there is much to be gained with unsymmetric stiffness matrices.

Summary of MacNeal's theorem
Before we start to develop Petrov-Galerkin EAS elements we summarize the key conclusions of the previous pages. First, a finite element can in general only be "exact" for a maximum of M = N DOF modes, where N DOF is the number of displacement degrees of freedom. By exact we denote a finite element solution which coincides with the continuum solution at the nodes (nodally exact response). Second, it is essential to fulfill nodal equilibrium to extend the solution capabilities of a single element to a whole patch of elements. This is ensured if the ansatz functions for the test function of the displacement are conforming. Finally, to reach exact representation of a maximum of M = N DOF displacement modes for general element geometries, the element's stiffness matrix K e needs to be unsymmetric. In case of a symmetric K e there is always a trade-off between lower and higher order modes (cf. MacNeal 11 ). All these conclusions hold for a broad class of finite elements. In particular, they apply to all well-known finite element discretization schemes regardless of their internal design. That is, regardless of whether they are based on mixed formulations, higher-order methods or reduced integration.

PETROV-GALERKIN EAS FRAMEWORK
Our investigations based on MacNeal's theorem presented in the last section indicate that the best performance for a given set of nodes can only be achieved with an unsymmetric stiffness matrix. However, the question on how to construct such a stiffness matrix has not been addressed.
In order to answer this question, we propose an Petrov-Galerkin extension of the enhanced assumed strain (EAS) method originally proposed by Simo and Rifai. 9 Our subsequent examinations of this basis yields three design conditions that enable the construction of high performance EAS elements.

Weak form
The key idea of the EAS method proposed by Simo and Rifai 9 is to recast the strain in the form = ∇ s x u +̃. Therein ∇ s x u and̃denote the compatible and incompatible part of the strain, respectively. Inserting this ansatz into the three-field Hu-Washizu 56 functional yields the variational framework for the EAS method (see Simo and Rifai 9 for details). In the present work we start with the discrete weak form and assume that the independent stress has already been eliminated via a suitable L 2 -orthogonality condition. The weak form is then given by 9 where v h ∈  h , u h ∈  h as well as̃h ∈  h ,̃h ∈  h denote the finite element approximations of test and trial functions of the displacement and incompatible strain, respectively. Furthermore, define the constitutive stresŝh as well as the contribution of the external forces to the weak form G h ext in analogy to G * ext given in (3). The external forces b h and t h are the finite element approximations of the continuum counterparts b * and t * .
Since the ansatz for the enhanced assumed strain,̃h ∈  h ,̃h ∈  h , does not need to be inter-element continuous, it is as usual introduced element-wise, which later allows static condensation of the internal degrees of freedom on element level. The weak form (23) is so far exactly the same as the one presented in Simo and Rifai. 9 The only difference in the present work is that we subsequently apply a Petrov-Galerkin approach instead of the usual Bubnov-Galerkin method. Thus, we use different ansatz functions for the test and trial spaces. The corresponding ansatz spaces  h ,  h ,  h,e ,  h,e will be specified in Section 4.3.

Design conditions for exact nodal solutions
The goal of this section is to find suitable design conditions for EAS elements which ensure that the solution of (23) for an analytic displacement state u * is nodally exact. That is, the nodal displacements coincide with the analytic displacement state according to (12). The first design condition stems directly from Section 2.2.3 where it has been shown that it is crucial to fulfill nodal equilibrium by proper choice of the test function for the displacement v h,e . This allows to consider only a single finite element in analogy to Section 2, since results may be generalized to larger patches of elements as shown in Section 2.2.3.
To determine further design conditions we start with the analytical weak form for a single finite element given in (9) and add the body force contribution, which has been neglected in (9).
Similar to the investigations in Section 2.2.2 we set the external loads in (26) . Finally substituting this result into (24a) and enforcing the weak form for a single finite element gives which determines the next condition for an exact finite element solution. Since we have seen that nodal equilibrium puts severe restrictions on the choice of the test functions v h,e (cf. the first design condition), it is in general only possible to fulfill (28) if̂h ,e = * holds pointwise, where, according to (1b), * = C ∶ ∇ s x u * . On the other hand, the constitutive stress (25) is given bŷh ,e = C ∶ (∇ s x u h,e +̃h ,e ). Thus, condition̂h ,e = * can be fulfilled by choosing the ansatz spaces for u h,e and̃h ,e appropriately. Furthermore,̂h ,e = * must be possible under the premise that the nodal displacements are exact. Otherwise it would be possible to fulfill (28) without having exact nodal displacements.
The third and final condition emerges from the second equation of the weak form. Substitutinĝh ,e = * (i.e., the second condition) into (24b) yields which determines that the test functions of the incompatible strain has to be L 2 -orthogonal to the stress * . Interestingly, this is an extended version of the patch-test condition of the classical EAS element by Simo and Rifai. 9 However, instead of requiring (29) only for constant stress fields we demand L 2 -orthogonality for all modes for which the element should be exact. Summarizing the above we established three conditions, which have to be fulfilled in order to construct a Petrov-Galerkin EAS element on the basis of (23) that is exact for a chosen displacement field u * regardless of element shape: C1 The test functions for v h,e have to ensure that nodal equilibrium is fulfilled regardless of element geometry and neighboring elements (see Section 2.2.3), C2 the discrete constitutive stress tensor̂h ,e , which is computed from u h,e and̃h ,e via (25), must include * under the premise that the nodal displacements are exact (see (28)) and C3 the test functions for the enhanced straiñh ,e must be L 2 -orthogonal to * (see (29)).

DESIGN OF MESH-DISTORTION INSENSITIVE EAS ELEMENTS
The final theoretical part of this contribution concerns the actual finite element design. We have to specify the analytic modes, for which the element should be exact, and then construct the ansatz spaces for all four fields introduced in (23) such that conditions C1 to C3 are fulfilled.
F I G U R E 3 Geometry maps and coordinate systems for the four node quadrilateral element. Isoparametric map x h,e ∶Ω  → Ω e (see (30)), skew map ∶ Ω e  → Ω (see (32)) and additional map ∶Ω  → Ω (see (35)) We start with the introduction of the skew coordinate frame in Section 4.1 since it is the key for the analytic stress modes in Section 4.2 and the finite element discretization in Section 4.3. For simplicity, we focus on the 2D case in Sections 4.1 to 4.3 and cover the extension to 3D in Section 4.4.

Skew coordinate frame
An element Ω e is described as usual by mapping a reference elementΩ transformation (see Figure 3) Therein, J h,e is the Jacobian of the transformation x h,e ∶Ω  → Ω e , x e i are the nodal coordinates and denote the usual bi-linear Lagrangian shape functions. They are defined with help of the nodal positions [ ] T in the parametric space.
In the context of assumed stress elements Yuan et al. 25 and Wisniewski and Turska 26 proposed to use a skew coordinate frame for an alternative description of geometry. The skew coordinates are introduced via the elementwise affine map ∶ Ω e  → Ω defined by is the Jacobian at the element center x 0 = x h,e| | =0 . The inverse of transformation (32) is given by Since the covariant base vectors g 0 = x h,e ∕ | | =0 and g 0 = x h,e ∕ | | =0 constitute the columns of J 0 , (33) can also be written as x h,e = g 0 + g 0 + x 0 . Accordingly, the coordinates and refer to the base vectors g 0 and g 0 , which in general span a skew basis (see Figure 3). It can further be observed from Figure 3 that the skew element Ω plays the role of an intermediate configuration between the reference elementΩ and the physical element Ω e .
Two important properties of the skew coordinates follow directly from (32). First, since the map (32) from Ω to Ω e is affine, the two elements Ω and Ω e are affine-equivalent (see, e.g., Ch. 12 in Reddy 24 ). This implies in particular that complete polynomials defined on Ω are mapped to complete polynomials in Ω e . This property will be exploited below for the construction of complete ansatz functions regardless of the geometry of element Ω e .
Second, the skew coordinates are frame-indifferent. This can be shown by considering superposed rigid motions corresponding to nodal position changes of the form x e# i = Rx e i + c, where R ∈ SO(2) is a rotation tensor and c ∈ R 2 is an arbitrary vector. Substituting x e# i for x e i in (30) yields x h,e# = Rx h,e + c and J h,e# = RJ h,e . With regard to (32), we then obtain which establishes the frame-indifference of the skew coordinates. This facilitates subsequently the construction of a frame-indifferent finite element formulation.
Remark 8. Although (32) links the skew coordinates to the parametric coordinates via x h,e and (30), it is convenient to have a direct relation between and . In particular, the map ∶Ω  → Ω is given by where T are the hourglass function and vector, respectively (see Appendix B for further details). It is straightforward to verify that c 1 = 0 for a parallelogram shaped element and therefore = in case of a constant Jacobian.

Prescribed analytic modes
The investigations of Section 2.2.2 show that the four node quadrilateral element with eight degrees of freedom can be fitted to a maximum of M = N DOF = 8 prescribed modes. Three of which are necessarily the rigid body modes, since elements that are not exact for those would lead to completely unphysical results. For the remaining five modes we propose to use the assumed stress modes developed by Pian and Sumihara. 8 These modes are suited perfectly since they include the important constant stress modes as well as the bending modes. The former are necessary to fulfill the patch test 2 and the latter are of the utmost importance in many engineering applications. Furthermore, the almost unchanged form of the assumed stress modes for more than 35 years suggests their optimality for assumed stress elements.
One of the few successful modifications of the Pian-Sumihara stress modes is formulating them in the skew coordinate frame. This has the advantage that equilibrium, that is, div h,e = 0, is fulfilled a priori regardless of element shape (see References 25,26). We therefore choose the analytic stress modes in the skew system, which are in vector matrix form given by as prescribed analytic modes. These stress modes have to be transformed to the physical frame to use them in (28). The usual assumption of a contravariant stress field yields the transformation * = J 0 * J T 0 .
The last two equations could be used to compute the corresponding displacement modes u * m , m = 1, … , 5 via the system of coupled PDEs (1). Fortunately, it is not necessary to perform this tedious task. To fulfill conditions (28) and (29) we only need the prescribed stress * . By proper choice of the ansatz space of the displacement and the incompatible strain we can ensure that the nodal displacements are exact without knowing the actual value of the analytic displacement at the node.

Ansatz spaces
The final task is to find ansatz spaces that fulfill conditions C1 to C3 for the prescribed modes (37) such that we obtain an element that is nodally exact for these deformation states.

Test function for the displacement
Condition C1, which concerns nodal equilibrium, is straightforward to fulfill. The only requirement is a conforming ansatz (21) for the test functions of the displacement v h,e . Since this is exactly what the isoparametric concept was designed for, we choose which employs the same bi-linear shape functions (31) as the ansatz for the geometry (30).

Trial function for the displacement
Condition C2 determines that the discrete stresŝh ,e has to include the prescribed stress modes * given in (37). Since the chosen stress modes (36) (25), the discrete constitutive stresŝh ,e is also a complete linear polynomial and contains in particular * . This approach is more general than explicitly needed for condition C2, since we can represent not only (36) but complete linear stress modes. However, its independence of the material model (i.e., C) is a crucial advantage.
Although the isoparametric concept ensures that complete linear displacement modes are correctly represented in the physical space, this is in general not the case for higher-order displacement modes. This deficiency is caused by the nonlinearity of the map (30).
Example Q1. The components of the map (30) can be written as x = a 0 + a 1 + a 2 + a 3 and y = b 0 + b 1 + b 2 + b 3 with constants a i , b i . Unless a 3 = b 3 = 0 (a parallelogram) the expression xy contains monomials up to 2 2 not included in the base of the isoparametric interpolation. The expression xy can therefore not be represented by the bilinear Lagrangian shape functions. Similar relations can be found for all isoparametric elements.
Complete polynomials in the physical space can be obtained by using shape functions constructed in the physical space, which are termed metric shape functions (see the book by MacNeal 1 ‡ ‡ ). The construction of metric shape functions starts by choosing an appropriate set of N monomials. For the four node quadrilateral (N = 4) the only sensible choice is 1, x, y, xy. Furthermore, it is necessary to construct the shape functions using the skew coordinates instead of the physical coordinates x to maintain frame-indifference (cf. Xie et al. 21 ). Despite using the skew frame, the shape functions can still be regarded as metric shape functions since the two elements Ω and Ω e are affine-equivalent (see Section 4.1). Following these considerations, we choose  The inverse matrix required to compute M e i exists as long as there are no coincident nodes. It is also well-conditioned even for elements with high aspect ratios because of the skew coordinates which scale the element appropriately. Naturally, both displacement components of u h,e have to be approximated using the same M e i to maintain isotropy. The final approximation of the trial function for the displacement reads in analogy to (38). Starting from (40) it is now straightforward to compute the derivative of the displacement component u h,e with respect to which is given by The last equation yields expressions for the derivatives of the metric shape functions ∇ M e i , i = 1, … , 4. Derivatives with respect to the physical coordinates follow from application of the chain rule together with (32) such that The derivatives ∇ M e i are obviously not complete linear polynomials which prohibits fulfilling C2 using only the displacement ansatz (41). It is further worth noting that the inverse matrices in (40) and (42) required to calculate the metric shape functions M e i and their derivatives ∇ M e i are the same and need to be computed only once per element. We summarize a few properties of the metric ansatz described above which are straightforward to verify using (40) before we conclude this section. Functions M e i are by construction frame-indifferent and isotropic. 21 The latter implies that M e i are invariant with respect to the node-numbering, that is, the direction of the skew basis vectors g 0 i . Furthermore, they are a partition of unity

4.3.3
Trial function for the enhanced strain Relation (42) reveals that the derivatives of the metric shape functions are not complete linear polynomials and have to be complemented in order to fulfill condition C2. Inspection of (40) shows that the missing expressions for a complete quadratic displacement field arem e 1 = 2 andm e 2 = 2 . A simple way to add these terms is considering two elementwise incompatible mode ansatz functions with internal weightsũ e j . In this way we obtain a displacement field complete in P 2 (Ω) given by However, directly choosingM e j =m e j does not fulfill the second part of condition C2 which states that identitŷ h,e = * must hold under the premise that the nodal displacements u e i associated with M e i are exact. In other words, the incompatible modes may not contribute to the displacement at the nodes. Choosing as ansatz in (44) ensuresM e j = 0 at the element's nodes while maintaining the quadratic completeness of the displacement field (44).
It is now straightforward to obtain a suitable ansatz for the trial function of the incompatible strain from (44) using the strain-displacement relationship of the form (1a). Accordingly, taking the derivative of the second term on the right-hand side of (44) results iñh where the incompatible displacement vectorsũ e j contain the four internal degrees of freedom.
Remark 4. The shape functionsM e j presented above are a more compact form of the ansatz proposed by Huang and Li 57 in a Bubnov-Galerkin frame. Furthermore, a similar idea of using metric incompatible modes has recently been suggested by Huang et al. 31 Their idea is to employ a higher order parent element which automatically ensuresM e j = 0 at the nodes. The drawbacks of that method are the higher number of internal degrees of freedom (at least twice as much), the higher-order of Gauss-quadrature required and the difficulty to extend it to 3D problems due to a violation of condition C3. All of these drawbacks are overcome with the present approach.

Test function for the enhanced strain
The remaining task is to find appropriate test functions for the incompatible straiñh ,e such that condition C3 holds. To that end, we start with the incompatible strain field of the classical EAS approach 9 where j h,e = det J h,e and e v contains the classical Wilson-modes in vector-matrix notation. The transformations chosen in (47) and (37) allow to recast condition C3, that is, (29), in the form This condition is fulfilled in case of the special choice (47). Thus condition C3 is immediately fulfilled by choosing h,e =ẽ h,e . However, in case of general ansatz functions, this is not necessarily the case. See for example, the 3D case described in Section 4.4.3.
Remark 5. Interestingly, (48) offers an explanation why transformation (47), which works well in linear simulations, 9 is the best choice of transformation of the enhanced modes in nonlinear simulations even though it is not the obvious choice if the deformation gradient is enhanced (cf. Pfefferkorn and Betsch 47 ): Since any nonlinear element should be close to the linear case for moderate deformations we may conclude that (48) is required for good bending accuracy even in nonlinear simulations. In case of rectangular elements (48) is fulfilled for the standard (transposed) Wilson ansatz functions if transformation (47) is used. This holds only approximately for distorted meshes but still improves performance. Other transformations lead to more severe violations of (48) which implies worse overall performance.
This effect is described in numerous works 40,47,58 with numerical investigations but could, to the best knowledge of the authors, not be explained so far.

Exact numerical integration
Similar to other Petrov-Galerkin finite elements 20 standard Gauss-quadrature allows to evaluate the integrals in the weak form (23) exactly.
To show this we first note that the constitutive stresŝh ,e contains only polynomials of x, y since we choose metric interpolations of u h,e and̃h ,e . In view of (35), these expressions are also simple polynomials of , . Next, transformation of the integral in (24b) to the reference element is possible using dV = j h,e dΩ. Due to the specific transformation in (47) j h,e immediately cancels out which establishes the polynomial form of the integrand in (24b). The same holds for (24a). Here we can recast the gradient of the isoparametric shape functions in the form where cof (•) denotes the cofactor. 59 This shows that j h,e cancels analogously after transformation of the integral in (24a) and the remaining expression again contains only polynomials of , . By inspection of the integrands we find that the highest order monomials are 2 , 2 2 , and 2 , all of which are integrated correctly by standard 2 × 2 Gauss-quadrature. Thus, standard integration is exact for the present element regardless of element shape.

Connection to the standard EAS element
If an element is parallelogram shaped, the present approach reduces to the standard EAS method presented in the work of Simo and Rifai. 9 In that case the Jacobian of the isoparmetric map J h,e = J 0 is constant and (35) implies = . We then obtain from (40) that M e i = N i . Finally, (45) become the classical Wilson-bubble modes which establishes the equivalence.

Eliminated orthogonal stress field
The weak form introduced in (23) is a reduced form of the full Hu-Washizu framework as mentioned in Section 3.1. The complete weak form on element level is given by 9 where in contrast to (23) the stress field h with its test functions h has not yet been eliminated via L 2 -orthogonality condition. It is straightforward to get from (49) to (23). Simply choosing h,e = * for the independent stress field on element level, where * is given by (37) and (36), eliminates h from (50b) due to (29). Moreover, it is simple to construct an ansatz for the test functions for the independent stress h such that (50c) is fulfilled by construction. To that end we simply start by taking (37) and replace the skew coordinates in (36) with the parametric coordinates . This set of ansatz functions can then be orthogonalized tõh ,e with a Gram-Schmidt orthogonalization similar to Section 4.4.2. Thus, it is possible to use (23) instead of (49).

Three dimensional problems
So far, we described the 2D case in a detailed way. However, most of the ideas and methods presented are straightforward to transfer to the 3D setting. We list differences and specifics for the 3D case in this section.

Prescribed stress modes
Similar to Section 4.2 we use the stress modes from assumed stress elements as prescribed analytic modes. They were first proposed in the work of Pian and Tong 50 in the parametric space. In skew coordinates (32), the assumed stress modes are defined by using (37) Again, the analytic modes contain the patch test ( 1 to 6 ) and bending modes ( 7 to 12 ). Additionally, the three remaining linear modes ( 13 to 15 ) concern the shear components and cover states of pure torsion. § § The final three modes are bilinear and were first proposed by Pian and Chen. 60 They are required to suppress zero energy modes in assumed stress elements. These modes fulfill equilibrium (2) with b * = 0 but do not lead to a compatible strain field in case of isotropy. ¶ ¶ This is why they are labeled incompatible stress modes.
It would of course be better to find suitable compatible stress modes. Unfortunately, these do, to the best knowledge of the authors, not exist. Compatible modes would need either full quadratic fields or a priori knowledge of the material model. Both of which are not satisfactory which is why (51) is employed in most works (e.g., References 50,61,62). Other approaches 63,64 perform either worse or at best similar compared to (51).

Trial functions for displacement and incompatible strain
The 3D metric shape functions M e i , i = 1, … , 8 for a eight node hexahedral element are defined analogously to (40) and contain the monomials 1, , , , , , , and . Again we obtain complete quadratic fields in skew coordinates by complementing the metric shape functions with quadratic incompatible modesM As consequence of the complete quadratic displacement ansatz all linear stress modes (51) up to 15 can be exactly represented regardless of the material model.
Approaching the remaining bilinear stress modes 16 to 18 in the same way would require nine additional cubic incompatible displacement modes which are by far too many for an efficient element. Furthermore, the bilinear stress modes are incompatible and of subordinate importance due to their higher order nature. We therefore propose to follow the approach of Simo et al. 39 instead and focus on volumetric locking.
Andelfinger et al. 65 first observed that using only the three quadratic incompatible modes described above is not enough to eliminate volumetric locking. Following Simo et al. 39 this can be explained by looking at the trace of the strain tensor, which assumes the form tr( h,e ) = tr(∇ s x u h,e +̃h ,e ) where C i andC i are constants depending on the displacement and enhanced degrees of freedom, respectively. The last equation immediately reveals that the incompressibility condition tr( h,e ) = 0 can only be met if C 4 = C 5 = C 6 = 0. This restriction leads to a mild form of locking which can be alleviated using three additional enhanced modes h,e = J −T 0 ( 10 I + 11 I + 12 I)J −1 0 (53) § § See for example, the work of Reference 52 for analytic solutions associated with the linear modes. ¶ ¶ It is not possible to find a displacement field that leads to * = 16 yz(e x ⊗ e x ). The compatibility condition * 11,23 = * 13,21 + * 12,31 − * 23,11 is violated for the simple isotropic material which ultimately prohibits finding a suitable u * . in the spirit of Simo et al. 39 It is straightforward to verify this by observing that (52) now includes the additional terms C 4 +C 5 +C 6 which relaxes the incompressibility conditions from C m = 0 to C m = −C m , m = 1, 2, 3. In contrast to Simo et al. 39 we use skew coordinates in (53) and do not apply the scaling by determinants. With these three additional modes there are now in total 12 enhanced modes. Remark 6. Other approaches for volumetric locking-free EAS elements require either more internal degrees of freedom 18,41,66,67 or are only slight modifications of (53) 46 .
All in all, the three quadratic incompatible Wilson-modes formulated in terms of skew coordinates together with (53) enable exact representation of arbitrary linear stress fields and volumetric locking-free elements.

Test function for the incompatible strain
The test functions for the incompatible strain are chosen in exactly the same way as in the work of Simo et al., 39 which are given by along with the transformation formula in (47). Unfortunately however, (48) is not a priori fulfilled for the bilinear modes in (51) and (54). It is thus necessary to find a proper L 2 -orthogonal set of test functions for the strain modes. This can be achieved by starting with an arbitrary e v and a two step Gram-Schmidt orthogonalization procedure. The stress modes are orthogonalized to each other in a first step before these modified stress modes are used to create an orthogonal enhanced field. It is crucial to conduct these operations using (48) instead of (29). This is because the latter leads to coupling between the different stress components which in turn induces frame dependence to the orthogonalized stress. The final orthogonalization in vector-matrix notation is described by where v,i and e v,j denote columns of the ansatz matrices of * v and e v , respectively. Furthermore, n and n are the number of stress and incompatible strain modes and k = 1, … , n and j = 1, … , n . The final ansatz for the test functions of the incompatible strain is given bỹh

Additional features
It is straightforward to show that the 3D element presented reduces to the linear version of the EAS element proposed by Simo et al. 39 in case of regular meshes. Furthermore, as in Section 4.3.5, all integrands in (23) and (55) are of pure polynomial nature. However, standard 2 × 2 × 2 Gauss-quadrature is not sufficient for exact integration of the higher order expressions in (51), (53), and (54) in case of distorted elements. Fortunately though, all lower order modes are integrated exactly by standard 2 × 2 × 2 Gauss-quadrature and the under-integration of the higher order modes does not seem to affect the results as our numerical experiments show. It is furthermore interesting that increasing the order of quadrature does not lead to over-integration and too stiff response which is in contrast to standard finite elements.
Remark 7. In case of generally distorted elements with non-constant Jacobians 4 × 4 × 4 and 3 × 3 × 3 Gauss-quadrature is required to exactly integrate the expressions of the orthogonalization procedure (55) and the weak form (23), respectively. Using different order of integration for the weak form and the orthogonalization procedure is not advisable as it leads to numeric violations of (29). The required order of Gauss-quadrature can be reduced to 3 × 3 × 3 in the orthogonalization procedure (55) if one replaces the skew coordinates in the bilinear modes in (51) with parametric ones, which does not seem to influence the numeric results.
However, we stress that despite not giving exact results, standard 2 × 2 × 2 Gauss points is sufficient for convergence and we thus employ it if not stated otherwise.

NUMERICAL INVESTIGATIONS
This final section covers several numerical benchmarks comparing the novel finite element formulation to established ones. Special emphasize is put on mesh distortion sensitivity. To keep the descriptions of setup and results as concise as possible, we drop the 2D case in favor of the more complex 3D simulations. However, if not mentioned otherwise, results can qualitatively be transferred to 2D.
In the sequel we denote the present element Q1U/E4 and H1U/E12 in 2D and 3D, respectively. Therein, U denotes the unsymmetry of the stiffness matrix and Ex accounts for the number of enhanced modes. The standard elements we use for comparison are • H1 standard isoparametric displacement-based element with Lagrangian shape functions.
• H1/E9 EAS element based on the nine Wilson-modes. The 2D version of this element, Q1/E4, was proposed in the seminal work of Simo and Rifai. 9 A 3D extension can for example, be found in the work of Andelfinger and Ramm. 41 • HA1/E12 linearized form of the EAS element proposed by Pfefferkorn and Betsch 47 which is closely related to the improved EAS element by Simo et al. 39 It uses three additional volumetric enhanced modes, a special nine point quadrature rule and a modified evaluation of the compatible displacement gradient.
• H1/S18 assumed stress element with 18 stress modes as proposed by Pian and Tong 50 .
• H1/P0 the three field mixed pressure element with constant pressure and dilatation approximation as proposed by Simo et al. 68 • H1U/IM-S the 3D extension of the unsymmetric incompatible mode element by Huang et al. 31 We use a serendipity approach instead of the proposed Lagrange ansatz and perform the orthogonality procedure (55) to fulfill condition C3, which has not been considered in Huang et al. 31 Other low-order unsymmetric finite element formulations as the ones proposed in Cen et al., 13 Xie et al., 21 and Huang et al. 31 (only 2D) have also been tested. Their accuracy is very similar to the present element in most tests. However, they suffer from other severe drawbacks as outlined in the introduction and are therefore not included in this section. The only exemption is H1U/IM-S (see above) in the Cook's-membrane example in Section 5.5.

Patch test
The first benchmark is the classical patch test which is performed here as described in Pfefferkorn and Betsch 47,48 .
A block  0 = [0, 1] 3 with Young's modulus E = 10 4 and Poisson's ratio = 0.3 is subjected to boundary conditions u i (x i = 0, x j , x k ) = 0, i, j, k ∈ {1, 2, 3}, i ≠ j ≠ k and loaded via displacement u 3 (x 3 = 1, x 1 , x 2 ) = u. A regular mesh with 3 × 3 × 3 elements as well as a distorted mesh with seven elements according to MacNeal and Harder 69 is used to discretize the block. Figure 5 shows the deformed configuration and stress distribution for both types of meshes. By considering the stress distribution we can establish if the elements pass the patch test, that is, if the linear displacement states lead to the correct constant stress state. The nodal stresses of the novel H1U/E12 at a displacement state of 6.000000 6.000000 6.000000 6.000000 6.000000

F I G U R E 6 Setup of the isotropy and frame-indifference test in 3D
u = 0.6 average to = 6000 for both meshes, which is the analytic solution. Moreover, the standard deviation of below 10 −10 establishes the constant stress state. Thus, the novel element passes the patch test, which is not surprising since H1U/E12 was designed to pass the patch test by choice of the analytic stress modes (see Section 4.2).
Remark 8. Other quadrature rules than the standard 2 × 2 × 2 rule might induce failure of the patch test. Using for example, 2 × 2 × 2 and 3 × 3 × 3 Gauss-points for integration of the residual (23) and orthogonalization procedure (55), respectively, leads to a violation of the patch test. An equal number of Gauss-points for both the residual and the orthogonalization (e.g., the 2 × 2 × 2 for H1U/E12) ensures satisfaction of the patch test. Alternatively, an exact integration, which requires at least 3 × 3 × 3 and 4 × 4 × 4 Gauss-points, respectively, also ensures that the patch test is passed.

Isotropy and frame-indifference
With this simple one element test we check whether the novel element formulation fulfills two basic properties. First, we test for frame-indifference with respect to rotation of the coordinate system and second for isotropy with respect to choice of the local element axes. This test was proposed in the works of References 21,64 and is slightly modified here ## . We consider a single finite element of general shape whose geometry in the local x y z coordinate frame is described in Figure 6. Dirichlet boundary conditions u = 0 apply on nodes A, D, E, H and two nodal forces act on nodes C and G as shown in Figure 6. The material properties are E = 1500 and = 0.25.
To test for frame-indifference we consider the case of no rotation, that is, the local x y z-frame coincides with the global xyz-frame, and a case where the element is rotated by 65 • , 15 • , and 25 • around the x, y, and z axis, respectively. Furthermore, we use different element node-numberings (see Figure 6) to change the local element axes with respect to the coordinate frame x y z. If this does not influence the result the element is isotropic.
## In particular, we use a different geometry to ensure generality of the test. The nodes chosen in references 21,64 are symmetric to the xz-plane which should be avoided.
All standard elements and the novel element H1U/E12 pass this test and are indifferent to rotation as well as node-numbering. Specifically, the novel element exhibits a displacement of u = 1.3696 for all six combinations of rotation and node-numbering. The standard deviation of the six displacements is 5.4 ⋅ 10 −16 .
Remark 9. Isotropy and frame-indifference of H1U/E12 are a direct consequence of using skew coordinates as described in Section 4.1 which are extensively used for both trial fields u h,e and̃h ,e .
If we used the centroid coordinates x = x h,e − x 0 instead of the skew coordinates (32), the element would loose its frame-indifference. In our numerical test this is established by the increased standard deviation of 1.7 ⋅ 10 −3 . This frame-dependence of the centroid coordinates has first been cured by a rotation. 23 Unfortunately, that approach induces anisotropy (cf. Xie et al. 21 ). So far, the only proper choice is using the skew coordinates which simultaneously scale and rotate the centroid coordinates, (see (32)).

Eigenvalue analysis
Next we perform the eigenvalue analysis of the stiffness matrix as proposed by Simo et al. 39 A single finite element with regular or distorted shape as shown in Figure 7 is considered for this test. The shear and bulk modulus are = 1 and = 10 9 which corresponds to nearly incompressible material behavior and causes volumetric locking for standard elements.
After static condensation, the stiffness matrix of the single finite element K e is always of the size 24 × 24. Thus, K e has 24 eigenvalues which can be grouped into four sets. First, zero eigenvalues i = 0 which correspond to rigid body modes. Second, locking modes with i → +∞. In the numeric examples infinity is not reached and therefore modes with "high" eigenvalues are considered as locking modes. The remaining eigenpairs are normal modes if their eigenvalues i > 0 are positive and instable if i < 0.
The desired result is that the element exhibits the correct number of six rigid body modes, has no instable modes i < 0 and exactly one locking mode which corresponds to pure dilatational deformation and accounts for the incompressibility. The other modes should be normal modes with i > 0.
The number of modes in each category for the elements tested are shown in Table 1. Apart from elements H1 and H1/E9 all elements are free from locking in this test and exhibit the desired behavior. Failure to reproduce the correct number of modes immediately establishes that the element suffers from locking. For example, the four and five locking modes of H1/E9 in the regular and distorted case result in a mild form of locking which influences for example, limit loads in elasto-plastic simulations as first reported by Andelfinger et al. 65 Unfortunately, passing this test is not sufficient to prove that an element is locking free. It is for instance well-known that H1/P0 suffers from shear locking yet it passes the test. It is much more difficult to detect shear locking using this simple test. Even if the block is scaled in one direction to mimic a shell-like geometry there are no additional locking modes for H1/P0. The only difference compared to shear locking-free elements is that a few eigenvalues are slightly higher but not clearly distinguishable.
The novel element H1U/E12 passes this test by exhibiting the correct number of modes in each category. It is also shear locking-free which we show with subsequent examples. However, due to the unsymmetric stiffness matrix it is not guaranteed that all eigenvalues are real-valued. Indeed there is one pair in the distorted case with i = 2.177 ± 0.02374i where the imaginary part is not negligibly small. So far, there seems to be no consequence of this in all the simulations we conducted. Interestingly, no imaginary eigenvalues occur in the 2D case even though the stiffness matrix is unsymmetric as well.

Mesh distortion
To examine the influence of mesh distortion we consider a 3D version 47 of the widely used mesh distortion test (e.g., the works 9,13,31,41 ). The cantilevered beam-like structure is subjected to pure bending as shown in Figure 8. The specimen of length l = 10, height h = 2, and width b = 1 is meshed with only two elements and distortion is applied via parameter s. At the clamped end of the beam boundary conditions u(x = 0, y, z) = 0, v(x = 0, y = 0, z) = 0, and w(x = 0, y, z = 0) = 0 apply. The free end is subjected to a bending moment M = 20 which is applied via traction boundary condition (z) = 30 ⋅ (1 − z). We use material parameters E = 1500 and = 0.25. An analytic solution exists for this problem and can for example, be found in the work of Jabareen and Rubin. 52 The exact value for displacement in z-direction at point (10, 0.5, 2) is = 1.
Results of the mesh distortion test are shown in Figure 9 where we compare the numerically computed displacement for increasing skew s. We limit the degree of skew to s = 4.9 instead of the usual s = 5 since the matrix needed to compute the metric shape functions in (40) becomes singular for coincident nodes.
Apart from element H1/P0, which suffers from shear locking, all elements are exact for a regular mesh with s = 0. However, for distorted meshes the computed displacement of the standard mixed finite elements drops to less than 60% of the required = 1 for distortions s ≥ 2. Of the standard finite elements the best results can be obtained with the assumed stress element H1/S18. There exist modified symmetric assumed stress elements that exhibit even better accuracy and show exact displacement for some nodes. We mention the elements proposed by Wu and Chueng 14 and Sze. 15 However, we stress that since the displacement is not exact at all nodes (see Table 6 in the work 15 ), these elements do not breakthrough MacNeal's theorem.
In contrast to the symmetric elements the novel H1U/E12 is exact regardless of the mesh distortion parameter s. The second plot in Figure 9 shows that there is only a small numerical round-off error compared to the analytic solution. This error slightly increases for higher distortions since the matrix needed to compute the metric shape functions is less well conditioned for highly distorted meshes. Moreover, the newly proposed element does not only reproduce the correct displacement at every node but is even pointwise exact. This holds for the displacement, the stress and the strain.  As pointed out in Section 2.2.2, it is not possible to obtain exact solutions to quadratic problems with low-order elements for all geometries. We can observe this by scaling the distorted mesh from the patch test in Section 5.1 for this test (see Figure 10, mesh 1). In contrast to the mesh in Figure 8 the skew axes of the elements do not align for this mesh. As a consequence, the displacement = 1.004 obtained with H1U/E12 differs marginally from the analytic value = 1. However, it is still a much better approximation than any of the standard elements can provide. Of those H1/S18 shows the best accuracy with = 0.146, which is more than 85% off in contrast to the 0.4% of H1U/E12 (see Figure 11). Furthermore, H1U/E12 exhibits the analytic results in case of the slightly regularized mesh 2 which is shown in Figure 10 on the right. For that mesh the skew axes in longitudinal direction of the beam align which implies that pure bending around the y-axis is again part of the ansatz space.
All in all, the novel Petrov-Galerkin EAS element H1U/E12 was designed to be accurate for bending problems and therefore it performs extraordinarily well in this test.

5.5
Cook's membrane  Figure 12 is clamped (u(x = 0, y, z) = 0) on the left and subjected to equally distributed shear stress = 100 in y-direction on the right face. We assume near incompressibility and set the material parameters to E = 2261 and = 0.4955. There are always two elements in direction of the thickness while we gradually refine the mesh in the other two directions using n el = {2, 4, 8, 16} elements per side. Figure 13 shows the displacement u of the top right corner and von Mises stress at the mid point of the lower surface in dependence of the number of elements. The stress is extrapolated from the Gauss-points to the nodes using an L 2 -smoothing procedure given by Therein, N i are the standard parametric shape functions, j the stress at the nodes and̂the constitutive stress |||| within the element. We use the same Gauss-points for the stress extrapolation as for the integration of element stiffness matrices and residuals.
The plots in Figure 13 show that the novel element exhibits the most accurate displacement and stress for all mesh sizes. It is even slightly better than the other unsymmetric element H1U/IM-S. Especially the accuracy for very coarse meshes is extremely high compared to the other elements. For finer meshes the difference becomes smaller because the deformation state within single elements converges towards the patch test state for which all elements are exact. Furthermore, the principle of limitation discussed in Remark 2 prohibits higher order convergence beyond (h 2 ) in the limit of fine meshes.
Nevertheless, the present example shows that the novel element can also improve results if the sought displacement state is not part of the ansatz space. However, the improvements are in general not as impressive as in bending dominated examples. After all, it is a low-order numerical method with its limitations.
Another aspect we examine with the Cook's membrane benchmark is the computational effort needed for the newly proposed H1U/E12 in comparison with the symmetric element HA1/E12. To that end we focus on the time taken to solve the global linear equation system since it usually requires most of the computing time of a finite element program. 71 We |||| For H1/S18 we use the independent stress field. employ Matlab's 72 "∖"-operator, which automatically selects a suitable efficient and highly optimized algorithm. In the present case it uses a sparse Cholesky and LU-factorization for the symmetric and unsymmetric stiffness matrix, respectively. Other parts of the finite element program are of subordinate importance and are either equivalent for both methods (assembly, pre-& post-processing) or hardly comparable due to non-optimized code (element-routine, computation of shape functions). These parts are therefore not considered here. Figure 14 shows the median relative computation time of H1U/E12 compared to HA1/E12 of 50 runs for different mesh sizes with up to 114,075 degrees of freedom. The run times are highly volatile in particular for fine meshes. There are many outliers which is why we considered the median instead of the mean value. Nevertheless, Figure 14 gives an idea of the additional computational effort which seems to be relatively independent of the mesh size and is about four times the effort of a symmetric method for the same mesh size. However, as Figure 13 and the other examples of the present work suggest, less elements might be required with the novel approach.

CONCLUSION
We designed a novel high performance Petrov-Galerkin mixed finite element based on the EAS framework. The rather uncommon Petrov-Galerkin approach is motivated by MacNeal's theorem, 11 which is revisited in the present work and in particular states that the optimal element performance can only be obtained with an unsymmetric stiffness matrix. Within the newly-proposed Petrov-Galerkin EAS framework, we devised three design conditions which ensure that the finite element is capable of exactly reproducing nodal displacements corresponding to specific deformation modes. The proposed element is designed to fulfill these conditions for rigid body, patch test, and bending modes, which are of the utmost importance in many engineering applications. For the construction of the ansatz spaces skew coordinates play a crucial role to obtain the novel frame-indifferent element with exceptional coarse mesh accuracy and low sensitivity to mesh distortion. In case of the patch test and bending problems, for which the element is specifically crafted, the element's response is even point wise exact not only for the displacement but also for the stress and strain. Furthermore, in case of regular meshes it is equivalent to the standard EAS approach.
Previously developed low-order unsymmetric elements suffer from major drawbacks that could be overcome with the present Petrov-Galerkin EAS-framework. In particular, the new approach does neither involve ansatz functions dependent on material parameters 13,21 nor does it require many internal degrees of freedom and/or higher order numerical integration. 31 Its extension to 3D is also straightforward in contrast to many existing elements.
In practical simulations the key feature of the novel EAS-framework is probably its reduced mesh-distortion sensitivity. This might help to reduce the effort to generate finite element meshes since the high quality demands on mesh topology can be reduced.
A drawback of unsymmetric elements in general is the increased numerical effort compared to standard methods. This concerns in particular the computation of the ansatz functions and the factorization of the unsymmetric stiffness matrix. However, this additional effort is in the opinion of the authors more than compensated by the increased accuracy.
Future work should first focus on extending the method presented to nonlinear problems. Even though it is not possible to find general solutions for the large variety of material models, our preliminary results show significant improvements, especially regarding sensitivity to mesh distortion. Other interesting lines of research are the extensions to shell problems as well as dynamics.

ACKNOWLEDGMENTS
This work was funded by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation)-project number 466086399 (BE 2285/16-1). This support is gratefully acknowledged.

DATA AVAILABILITY STATEMENT
The data that support the findings of this study are available from the corresponding author upon reasonable request. The first simplification made in Section 2 was considering only modes without body force (5). However, higher order modes can in general not be represented without appropriate b * computed using (2). These non-zero body forces b * ≠ 0 have to be considered in (6). Fortunately, it is still possible to write the external part of the weak form as where F * mn now includes contributions form t * and b * . Since there is no change of the internal part (8), all conclusions on the structure of F * mn remain the same. Thus it is symmetric and by proper choice of u * m diagonal. The external part of the weak form from the finite element solution (17) changes accordingly with the only necessary change being the addition of the body force integral. Furthermore, there is no change in the conclusions drawn for nodal equilibrium. This follows from the fact that the body force concerns only the internal of the element while nodal equilibrium is determined by balancing the forces on the surface. All in all, it is straightforward to consider modes with b * ≠ 0 by simply adding the missing integral in all external parts of the weak form.

A.2 PETROV-GALERKIN APPROACH
Choosing the Bubnov-Galerkin approach in Section 2.1.2 does not limit the generality of the present proof. Any proper set of M linear independent modes v * m for the test function yields a generally unsymmetric F * mn as (6) and (8) show. However, by use of the singular value decomposition and transformation of the weights n , m using the unitary matrices of said decomposition, the virtual work of the external forces can always be expressed in the form (6) with a diagonal matrix F * mn . Thus, it is possible to get the same structure with the Petrov-Galerkin approach as with the Bubnov-Galerkin approach and all further conclusions apply.

A.3 GENERALIZED DISPLACEMENT APPROXIMATIONS
As mentioned in Section 2.2.1 the approximation (11) is not very general and especially excludes cases of high practical importance where M e i do not have the Kronecker-delta property and thus u e i are not the nodal displacements. Widely used examples of such non-interpolatory approaches are hierarchical higher order elements 2 and the isogeometric analysis. 3 Fortunately, it is still possible to find unique u * ,e im such that u e i can be computed according to (12). Considering an even more general displacement ansatz with matrix-valued ansatz functions such that still allows to compute u e i from equation system u h,e (x j ) = u * (x j ), j = 1, … , N if the fixed locations x j are chosen appropriately (usually the nodes). Examples for matrix valued ansatz functions for the displacement can be found for example, for Petrov-Galerkin finite elements in References 13,21,32,37 or as special cases of incompatible mode elements in Reference 51.
The general approximation (A2) can also be applied to the test function. If this is done, the forces P h,e i defined in (17) should not be termed nodal force anymore since the do not directly correspond to a single node anymore due to the lost Kronecker-delta property. Accordingly, the expression nodal equilibrium should not be used. We use the expression generalized nodal force and equilibrium instead. Apart from this terminology, however, there are no changes.
In view of the restrictions due to equilibrium of generalized nodal forces, which apply accordingly, it is questionable if matrix valued functions are sensible for the test functions of the displacement. As shown in Section 2.2.3 nodal equilibrium severely restricts the choice of ansatz functions and one does not gain much flexibility with the matrix valued functions if standard assumption (21) is made. Furthermore, the matrix valued functions are restricted due to isotropy requirements and must be chosen such that changes of coordinate system do not alter the results. Ultimately, we conclude that using scalar functions for the test function of the displacement is most likely the only reasonable choice.

APPENDIX B. SKEW COORDINATE FRAME
This Appendix describes how the direct link between the skew coordinates and the isoparametric coordinates given in (35) can be established. We present the 3D case which can easily be reduced to the corresponding 2D formulation. By virtue of (32) and (30) we get in a first step .

( B 1 )
Now we use an alternative representation of shape functions frequently used in the context of hourglass-stabilization (see e.g., the works 47,71