{'fraction_non_alphanumeric': 0.03559084806763975, 'fraction_numerical': 0.014182067468365277, 'mean_word_length': 4.554859182863943, 'pattern_counts': {'":': 0, '<': 4, '<?xml version=': 0, '>': 8, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 3, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Improving efficiency and robustness of enhanced assumed strain elements for nonlinear problems

The enhanced assumed strain (EAS) method is one of the most frequently used methods to avoid locking in solid and structural finite elements. One issue of EAS elements in the context of geometrically nonlinear analyses is their lack of robustness in the Newton–Raphson scheme, which is characterized by the necessity of small load increments and large number of iterations. In the present work we extend the recently proposed mixed integration point (MIP) method to EAS elements in order to overcome this drawback in numerous applications. Furthermore, the MIP method is generalized to generic material models, which makes this simple method easily applicable for a broad class of problems. In the numerical simulations in this work, we compare standard strain‐based EAS elements and their MIP improved versions to elements based on the assumed stress method in order to explain when and why the MIP method allows to improve robustness. A further novelty in the present work is an inverse stress‐strain relation for a Neo‐Hookean material model.


INTRODUCTION
In the early days of the finite element method (FEM), it was soon discovered that low-order purely displacement-based (U) finite elements yield poor results in many cases due to severe locking phenomena. Therefore, a plethora of mixed finite elements has been developed subsequently. Two of the probably most successful classes in linear analyses are enhanced assumed strain (EAS) elements introduced in 1990 by Simo and Rifai 1 as a generalization of the popular incompatible modes elements by Taylor et al. 2 and assumed stress (AS) elements proposed in the 1980s by Pian and Sumihara 3 and Pian and Tong 4 for two-(2D) and three-dimensional (3D) problems, respectively (see also the pioneering work of Fraeijs de Veubeke 5 ). Both classes exhibit excellent behavior in linear simulations. They are completely locking-free if the additional fields are approximated appropriately, they are stable and relatively insensitive to mesh distortion. In fact, there are hardly any drawbacks of using such elements in the linear elastic case. Furthermore, it is shown by Bischoff et al. 6 that for every EAS element an equivalent AS element can be found and vice versa. For nonlinear problems, however, the two approaches are not identical anymore. Here, AS elements are less popular, since they need an inverse stress-strain relation. With a few exceptions for simple material models (see e.g. Wriggers 7 ), this is extremely difficult or even impossible to obtain. In cases, where a solution exists, the recent numerical approach by Viebahn et al. 8 is a simple alternative to the cumbersome analytical inversion. However, there exist cases in which inversion of the stress strain relation is impossible (e.g., Ogden 9 and the present work). EAS elements do not encounter these difficulties, because their strain-based construction naturally fits the typically strain driven format of nonlinear material laws. Thus, they can easily be extended to general material models, also including geometric nonlinearity, see for example References 10-22, among others. Unfortunately, two major open issues of EAS elements remain. First, unphysical instabilities occur, depending on the chosen material model, ansatz functions for the enhanced field and state of strain. This was first discovered by Reese and Wriggers 23 and addressed in many subsequent works (notably References 12,17,18). The effect is also present for AS elements 8 and other locking-free formulations. The second open issue is lack of robustness in the Newton-Raphson (NR) solution algorithm. By the term robustness we herein refer to two properties: maximum size of applicable load steps and number of NR iterations required to find equilibrium. In this sense, robust elements are also efficient, since they require fewer iterations and therefore also fewer time consuming matrix factorizations. In this paper, only the second open issue is addressed, which seems to be a particular feature of EAS elements (AS elements, for instance, are superior in this respect). So far, this kind of algorithmic robustness has received little attention in the development of finite elements and focus has been put onto other topics, such as locking and stability. In the context of EAS elements the only works known to the authors that deal with robustness are References 8,12,[24][25][26][27]. For beam finite elements Garcea et al. 28 discovered that (AS-type) mixed finite elements provide superior robustness compared to displacement-based elements. This result has recently been reproduced by Magisano et al. 29 for solid-shell problems. In the context of continuum elements Viebahn et al. 8 obtained similar results while comparing AS and EAS elements. The authors showed that EAS elements have inferior robustness even if they yield the same equilibrium solution as AS elements. For another class of mixed finite elements, which are based on a Hu-Washizu-type variational functional, superior robustness of the mixed format compared to the displacement form has been observed by Wisniewski et al. 26,27 and Betsch et al. 30 One approach to improve the behavior of strain-driven elements is presented by Mei et al., 31 who propose to transform nonlinear finite element equations to make them more "linear." A simple example for this would be transforming the equation e x = a with the unknown x and a constant a by taking the logarithm on both sides. Unfortunately, this mathematically elegant approach is difficult to implement into practical simulations, since a suitable transformation is not straightforward to determine. Moreover, the main source of nonlinearity has to be known a priori, which becomes difficult if multiple sources of nonlinearity (e.g., geometric and material nonlinearities, constraints, plasticity) are present or if their corresponding impact changes during simulation.
In the context of shell problems, Magisano et al. 32 introduce the mixed integration point (MIP) method (see also References 33,34). The key idea is to introduce independent stresses at the Gauss quadrature points, which ultimately leads to a modification of the stresses used for the geometric stiffness matrix. A special feature of the method is that the residuum is not altered, which means that only the robustness of the method during iteration is (usually) improved without changing the converged result. The principal idea of modifying the stress in the geometric part of the stiffness matrix has already been proposed in the 1980s by Kuo-Mo et al. 35,36 They suggested to use an inconsistent tangent, by using the stress of the previous converged load step for all NR iterations of the next load step in order to compute the geometric stiffness contributions. However, the MIP update algorithm is more sophisticated and improves upon those results. In the context of contact problems, Zavarise et al. 37 proposed a similar philosophy of modifying the geometric stiffness contributions as a way of efficiently dealing with large penetrations and increasing the robustness of contact algorithms.
In the present work we use the idea of the MIP method introduced by Magisano et al. 32 and apply it to solid finite elements instead of structural (solid-)shell and beam elements, which are usually considered in the literature. 32,34 Furthermore, we propose three extensions to the MIP method. First of all, we consider general material models, since being restricted to the simple St. Venant-Kirchhoff material model is one of the current limitations of the MIP method. We discuss under which conditions the MIP method can be applied successfully to more advanced materials and when less beneficial results are to be expected. Secondly, we propose a transformation of the method to use spatial stress and strain measures in order to support computationally more efficient elements and allow simple implementation of many material models. Finally, and most importantly, we apply the method to EAS elements, which significantly increases their robustness and makes this class of elements even more interesting in practical simulations.
Various numerical examples, including elastic and plastic material laws are presented to highlight the properties and performance of the MIP method. With these we show, when and why the method allows the greatest benefit and when only little improvement is to be expected.
A final novelty of this work is an analytic inverse stress-strain relation for a Neo-Hookean material model, which can be used for AS elements. It demonstrates the limits of this class of elements in general nonlinear simulations.
The present work is structured into seven sections. General relations and notation for nonlinear continuum mechanics and finite element formulations of the displacement-based (U) as well as EAS and AS mixed finite elements are given in Section 2. Afterwards, a simple model problem is presented in Section 3 to show basic relations exploited for the MIP method. The MIP method and its extension to general material models are described in Section 4. Additional extensions and more details of the MIP method are given in Section 5. In particular Section 5.1 covers the MIP method for EAS elements while Sections 5.2-5.4 deal with its application to problems based on spatial quantities, the implementation of the method and a simple one-element test. Extensive numerical simulations follow in Section 6 before conclusions are drawn in Section 7. Appendix A covers the St. Venant-Kirchhoff, Neo-Hookean and elasto-plastic material models used throughout this work. Special emphasis is put on the invertibility of the stress-strain relation of the Neo-Hookean model (see Section A.2.2). More details concerning the simple model problem in Section 3 are given in Appendix B.

Continuum mechanics
The deformation of a deformable body from its reference configuration  0 to the current configuration  is described by the bijective deformation map which maps material points X ∈  0 to corresponding spatial points x = (X) ∈ . Therein,  0 denotes the part of the body's boundary  0 on which the deformations are prescribed by ∶  0 → R 3 . Linearization of the deformation map 1 at a point X yields the deformation gradient where the index denotes that F is computed from the deformations alone, which will not necessarily be the case in subsequent sections. The deformation gradient 2 F can be used to define the right Cauchy-Green tensor and Green-Lagrange strain tensor, respectively, as Furthermore, the constitutive second Piola-Kirchhoff stress tensor and its linearization with respect to E are introduced asŜ The notation Δ (•) ( * ) is used to denote linearization of ( * ) with respect to (•). The form of constitutive law forŜ covers a wide range of material models including inelastic behavior via internal variables . Linearization of the constitutive second Piola-Kirchhoff stress with respect to E is governed by the fourth-order material tensorĈ =Ŝ∕ E.
Remark 1. In case of a homogeneous hyperelastic material with strain energy function W(E), relations (4) are given bŷ However, we emphasize that hyperelasticity is no prerequisite for the novel techniques proposed in this work. 1 Subsequently, arguments of functions are often omitted in order to improve readability and notational simplicity. 2 Index is omitted here since the following relations are valid for general deformation gradients.

Displacement-based finite elements
In the sequel, we only consider numerical solutions of nonlinear solid mechanic problems obtained with the FEM. More specifically, a body  is approximated with n el four-node quadrilateral or eight-node brick finite elements Ω e in 2D and 3D, respectively. Furthermore, the isoparametric concept is applied. Approximations of geometry X and deformation within one element Ω e are given by where superscript h denotes an approximation and e marks elementwise quantities. However, these superscripts are usually omitted in the remainder of this work since most of the subsequent presentations are dealing with discretized fields on element level. Thus, if not specified otherwise, every quantity introduced in the sequel is an approximation on element level, for example ≡ h,e . In the equation above, n no is the number of nodes of the element and X e I , e I are the nodal values of the respective quantities. Moreover, N I ( ) denote the Lagrangian shape functions defined on reference elementΩ = [−1, 1] n dim for 2D and 3D approximations, respectively (e.g., Reference 7).
The internal part of the weak form of equilibrium of a purely displacement-based formulation on element level is given by (see e.g. 7,38 ) where is a kinematically admissible test function, approximated in the same way as in (6) (Bubnov-Galerkin approach) and (•) ( * ) denotes the first variation of ( * ) with respect to (•). The displacement-based Green-Lagrange tensor E and constitutive second Piola-Kirchhoff stressŜ are defined through (2), (3), and (4) as E = E(F ) andŜ =Ŝ(E , ), respectively. In these relations, index denotes the purely displacement-based nature of these quantities and index d specifies that the respective variable is part of the displacement-based FEM. 3 The external part G ext ( , ) of the weak form includes all external forces acting on  0 and the boundary t  0 =  0 ⧵  0 . They are not specified in detail here, since only the internal part (7) will be needed for the developments in this paper.
Considering the arbitrariness of and using Gauss quadrature to approximate the integral in (7), the linearized discrete form of (7) on element level reads where denote the integrand of the element residual as well as the integrand of the material and geometric part of the stiffness matrix, evaluated at the integration points (for more details see e.g., Wriggers 7 chapter 4.2.2). Abbreviation (•) g = (•)( g ) denotes evaluation of a quantity at Gauss point g = 1, … , n with weight w g . QuantitiesŜ andĈ in (9a) and (9b) have to be cast in Voigt (vector-matrix) notation, which is not separately marked since it becomes clear from the context when this notation is required. Superscripts in (8) and (9), indicating applied directional derivatives with respect to (w.r.t.) , are not necessary here but are merely introduced to get similar notation as for the mixed formulations presented in the sequel. The nodal operator matrix of the material displacement form is denoted by B d . We emphasize the dependence of G d,g on the constitutive stressesŜ ,g in (8), since this is the key to the MIP method presented in Sections 4 and 5.

EAS method
The first well-known class of multifield finite elements used in this work are EAS elements. The first EAS formulation for nonlinear problems was proposed by Simo and Armero. 10 In the present work, we give only a brief overview and refer to Pfefferkorn and Betsch, 21 among others, for a more thorough summary of many topics concerning EAS elements.
The key idea of the EAS method in the form used herein is to recast the deformation gradient in the form where F andF denote the compatible and incompatible (or enhanced) part of the deformation gradient. The latter includes additional enhanced degrees of freedom (DOFs) , which are used to improve the element's behavior with respect to locking. The variational basis for EAS finite elements is a Hu-Washizu 39 type functional. Usually, the independent stress is eliminated by enforcing the discrete stress and the enhanced strain field to be L 2 -orthogonal. This ultimately yields the approximated weak form 4 formulated using material quantities S and E in the form where the Green-Lagrange strain (3) and constitutive second Piola-Kirchhoff stress (4) are given by E = E(F ) and S =Ŝ(E , ), respectively. On top of that, index e is introduced to mark quantities of the EAS method in analogy to index d used in Section 2.2.1 for the displacement-based method. Note that in (11) only (11b) is equal to zero on element level without global assembly since the enhanced DOFs are introduced elementwise. Approximation of the deformation and compatible deformation gradient is applied as described in Section 2.2.1 and the enhanced deformation gradient can be discretized on element level using Wilson-modes 40 with nine (in 3D) enhanced parameters i byF In the last equation, J 0 and j 0 are the evaluation of the Jacobian of the isoparametric map J and its determinant j = det(J) at the element centroid = 0. Moreover, F 0 = F ( = 0) denotes the compatible deformation gradient at the element centroid and is needed to ensure objectivity of the formulation.
Remark 2. We refer to Pfefferkorn and Betsch, 21 among others, for many other possibilities to approximate both, the compatible and the incompatible part of the deformation gradient and for a thorough explanation of the requirements to be met. Furthermore, an overview of other ansatz functions for the enhanced deformation gradient is given in the aforementioned reference. Note especially the transposed Wilson modes which cure the spurious hourglass instability of the standard EAS element in compression.
The discrete linearized form of (11) in matrix notation is given by [ M e,g + G e,g (Ŝ ,g ) M e,g + G e,g (Ŝ ,g ) M e,g + G e,g (Ŝ ,g ) M e,g + G e,g (Ŝ ,g ) where the integrands of the element residuals and the geometric as well as material parts of the tangent are given by , G e,g = G e (Ŝ ,g ), (14b) , G e,g = G e (Ŝ ,g ).
Here, superscripts and denote the various parts of strain-displacement matrices, tangent and residual. Moreover, C is the material tangent (4) evaluated using E . Like in the previous section, the geometric part of the stiffness matrix depends on the constitutive stressŜ .
Remark 3. Usually, static condensation is performed on element level to eliminate the internal DOFs Δ e from (13). However, in this work the full form of (13) is needed to develop the MIP method, presented in Section 5.1, and static condensation is carried out afterwards.

AS method
The second class of mixed finite elements considered in this work are AS elements, which are based on a Hellinger-Reissner (HR) functional and employ separate stress approximation. However, with a few exceptions for simple material models, a complementary stored energy function, which is required for a HR-functional, does not exist. This is why usually a weak form is used as starting point for AS elements if more complex material models are to be used (cf. References 7,8). In the present work the discretized weak form is given by where the independent (not constitutive) stress tensor S is introduced as function of internal elementwise DOFs and index s denotes the AS method. Common approximations for S can, among others, be found in the work of Viebahn et al. 8 In this work, however, we apply a special stress approximation that is described in Section 4. The expression in (15b), like (11b), is zero on element level due to the elementwise approximations of S . The constitutive strainsÊ introduced in (15) are computed from the inverse of relation (4) such that This requires that (4) is at least locally invertible, which is given for most commonly used material models in the neighborhood of the stress-free reference configuration F = I, as described in chapter 6.2.2 of Ogden. 9 The linearization of this relation can then be computed via the law of differentiation of the inverse by where the tangential compliance matrixD ∶=Ĉ −1 (Ê , ) is the inverse 5 of the material tangent (4). 5 The inverse of a fourth-order tensor can easily be computed in Voigt notation where it reduces to computing the inverse of the corresponding matrix.
Remark 4. While most common material models are invertible in the neighborhood around the stress free configuration, analytic relations of inverse stress-strain relations rarely exist. Appendix A covers analytic inversions. Besides the simply invertible St. Venant-Kirchhoff model, we present a novel inverse stress-strain relation for a Neo-Hookean model (see Appendix A2). More complex models are difficult or even impossible to analytically invert. Recently, a numerical procedure for general hyperelastic models has been proposed by Viebahn et al. 8 Remark 5. Even though inversion is usually possible close to the stress free reference configuration, it is not necessarily so for states of larger strain. This is, for example, the case for the Neo-Hookean model described in Appendix A2. The states of stress, for which unique invertibility is given for that model, are described in detail in Appendix A2.2. These requirements are a severe limitation of the applicability of AS elements since these states can be observed in practical simulations as shown in Section 6.3, where the noninvertibility leads to failure of the NR scheme.

Problem description
In this section, a simple model problem is considered to demonstrate issues that may arise in the iterative solution process and to motivate the MIP method presented in subsequent sections. Following Magisano et al., 29 the simple nonlinear bar-spring system, shown in Figure 1, serves as model problem. The structure consists of a geometrically nonlinear bar with axial stiffness EA and a linear spring with stiffness k. We use a unit length L = 1 for the bar in order to simplify expressions. The upper end of the bar is subjected to a vertical load and a horizontal load c . We summarize here the most important relations for this example. More details are given in Appendix B. For a displacement-based formulation with two DOFs , in analogy to Section 2.1, the kinematic and constitutive relations are where (u, w) is the displacement-based axial Green-Lagrange strain andN(u, w) the constitutive "second Piola-Kirchhoff normal force." With these expressions at hand, the linearized displacement-based weak form of equilibrium (8) is given by For the AS formulation, similar to Section 2.2.3, the linearized version of the weak form of kinematics and equilibrium (15) reads where, for sake of simplicity, the independent "stress" field solely consists of the axial force N of the bar, whereas the linear spring remains unaffected. Furthermore, the constitutive strain̂is given by inverse stress-strain relation = N∕EA. Condensation of N reproduces Equation (19). However, the performance of the incremental solution process can differ significantly if the matrices in (20) are used, which will be shown next.
For both formulations, the results after the first iteration (j = 1) of the standard NR procedure are considered, starting from a stress-free reference configuration u (0) = w (0) = 0, N (0) = 0 and (0) = 0. For both formulations, the first displacement predictor due to an incremental load Δ is the same, While these values are identical for both formulations (which is not the case in general), the axial force may differ significantly. After the first iteration, the updated values for the constitutive and independent normal force arê which reveals the intrinsic problem of the extrapolated constitutive forceN(u, w) compared to the independent field N. At equilibrium points both force values are identical, but during the iteration process this is not necessarily the case. The extra terms Δ 2 2EA + Δ 2 c 2 EA 2k 2 emanate from the fact thatN is computed from the strains via the constitutive law. The strains, however, contain spurious contributions from the poor displacement predictor. For a high stiffness ratio EA/k 2 , the extrapolated forceN may be estimated far more inaccurately compared to the independent stress field N. In fact, N follows a quadratic path while N is extrapolated linearly. Since the extrapolated normal force affects the geometric part of the tangential stiffness matrix in the next iteration step (j + 1), the update of the NR scheme in this step will also be a worse guess in the displacement-based method. This ultimately explains why displacement-based finite elements may require more iterations within a NR scheme in comparison with AS elements, as demonstrated in Section 3.3.
As pointed out in Reference 29, a similar phenomenon can occur in the structural analysis of thin-walled structures, because of the high ratio of membrane stiffness to bending stiffness. This, especially for deformations involving large rotations (which is the case for the model problem herein), may spoil the extrapolated stress quantities tremendously. However, it is not relevant anymore when no rotations are involved, for example, in the case of purely axial stretching. Such parameters, which influence the robustness of the iteration process, will be denoted as critical parameters. They may be the same parameters that are responsible for locking. As shown in Section 5.4, material parameters, such as the Poisson's ratio, can be critical parameters, too. For nearly incompressible problems, small deviations from volume preserving states yield high deviations in stresses due to the activation of nonisochoric deformation modes. Remark 6. The EAS approach of Section 2.2.2 is not applied for this model problem, since it is identical with the displacement based approach (19). Any enhancement applied to a bar problem would vanish because of the constant strains assumed within the bar and the patch test requirements on the enhanced field. The equivalence of EAS and displacement-based approach implies especially, that the EAS method suffers from the same poor extrapolation of the axial force as the displacement-based approach.
Remark 7. In the context of asymptotic stability analyses the authors of Reference 41 denoted the phenomena as "extrapolation locking" (or "bifurcation locking"), since the issues described above can affect the discrete solutions in the sense that critical loads are overestimated. This denomination will not be used here, since for standard boundary value problems the aforementioned issues do, in contrast to "classical" locking phenomena, not affect the converged solution.

Modified displacement-based NR schemes
To improve the convergence properties of the displacement-based formulation (19), without using mixed methods, the principal idea is to modify the update of the stress (normal force, respectively) that is used for the geometric tangent. A simple strategy for a modified NR scheme is presented by Kuo-Mo. 35 The idea is to use the converged stress from the previous load increment for the entire iteration process of the current load increment. This is denoted as previous stress (PS) method in the following. The influence of inaccurately extrapolated stresses is avoided and the method yields a robust iteration process, independent of critical parameters. However, quadratic convergence of the method close to the solution is destroyed. This effect is particularly pronounced for large load increments, since the modified tangent is not close to the consistent one. Thus, in many situations a high number of iterations is necessary which can make the method uneconomical.
This drawback is overcome by another approach recently proposed by Magisano et al. 32 Although for this simple model problem their approach coincides with the AS formulation, it will be used to demonstrate how the method works. Its key idea is to take the AS approach (20) and conduct static condensation of the independent stress increment where ΔN = EA((1 + u)Δu + wΔw) = EAΔ . Then, for iteration (j), the update of the recovered normal force N (j+1) = N (j) + ΔN (j) can be rewritten solely in terms of constitutive stress quantities as which differs from the standard displacement-based stress updateN (j+1) = EA (u (j+1) , w (j+1) ). The constitutive force increment ΔN (j) can be interpreted as a linearized constitutive extrapolated force. It provides a reduced extrapolation error and facilitates the improved robustness. For the given model problem, the reduced linearized form, obtained after static condensation, is equivalent to the displacement-based version (19) with the exception ofN in the tangent being replaced by (24). In analogy to (31) in the next section, we may rewrite (24) in terms of the current displacement DOFs d (j) , which yields where B is the incremental strain-displacement operator matrix of the bar element. The authors of Reference 32 denoted this approach as MIP, since its derivation for shells and solids requires an AS formulation on Gauss point level, as presented in Section 4. For the simple model problem studied in this section, MIP and AS yield identical results, which is why they are used synonymously in (25). However, this is not the case when nonlinear constitutive relations or more advanced structural elements are considered, which will be emphasized in detail in the subsequent sections.
It should also be mentioned that the displacement-based formulation and the AS formulation yield the same equilibrium solution for this model problem. Thus, the MIP approach, which is basically a reformulation of the AS method, converges to this result as well. This means that the tangent close to the equilibrium is the same for both methods and thus quadratic convergence is maintained.

Numerical investigation
Next, a numerical study of the introduced model problem illustrates the differences between the displacement-based (U), previous stress (PS), AS, and MIP methods. For the fixed values c = 0.25 and k = 10 an applied load = 4.5 is still below the critical load factor for EA ≥ 10 2 , see Figure 1. Figure 2 shows the number of total NR iterations for different values of EA, needed for convergence if the load is applied within one load increment (left) and five load increments (right). It can be observed that the robustness 6 of the displacement based formulation depends on the axial stiffness EA (critical parameter). For increasing values of EA, the required number of iterations to fulfill the convergence criterion increases significantly. The results of the PS method are independent of the critical parameter EA, but the required number of iterations is high due to impaired convergence properties. In sharp contrast to that, the results of both the MIP and the AS method are independent of EA and significantly more efficient than the PS method.
These results are confirmed in Table 1, where the normal forces which occur within the iteration process are shown. At the beginning of the iteration process strongly deviating values, w.r.t. the equilibrated force, are obtained for the  displacement-based formulation, especially for large values of EA. An axial stiffness of EA = 10 4 for instance, even yields an artificially high constitutive tensile force in the first iteration step, emanating from the linearized (not yet correct) deformed configuration. In contrast to this, the values provided by the MIP and AS method are close to the converged one from the beginning, independent of the value of EA. Figure 3 visualizes the spaces of the normal forces spanned by u − w, as well as the values for the first 4 iterations (dots) for the case of EA = 10 4 and a single load increment. The black curve corresponds to the equilibrated normal force for 4.5 ≥ ≥ 0. The space of the displacement-based force (in red, left) is defined by (18) and remains the same quadratic surface for all iterations. On the other hand the space of the MIP/AS force (in blue, right) is defined by (25) and consists of different planes for each iteration step (j). Here, a hollow dot indicates the force N (j) MIP/AS , which is extrapolated at the current plane from the displacement-based forceN (j−1) (full dot) as shown in (25). Note especially, that the MIP/AS has visually converged to the black dot (equilibrium), while the displacement-based formulation is still far off that point. This shows the high robustness of the MIP/AS formulation and is in line with the results in Table 1. The essence of Figure 3 is, that the linear extrapolated force in every step of the MIP/AS approach is a better guess than the nonlinear extrapolation of the displacement-based formulation.

MIP METHOD FOR DISPLACEMENT-BASED FINITE ELEMENTS AND GENERAL MATERIAL MODELS
In this section, we first show the extension of the MIP concept to displacement-based solid finite elements, following Magisano et al. 32 Next, it is applied to general nonlinear material models, thus overcoming one of the major limitations of the MIP concept. The discrete weak form for AS elements (15) serves as basis for the MIP method. Usually, the stresses in (15) are approximated with elementwise continuous ansatz functions. In contrast to this, the key idea of the MIP method for shell and solid problems is to choose pointwise stress approximations. Therefore, independent stresses S g are  32 ). Using Gauss quadrature with weights w g , the discrete weak form on element level can be written as where n is the number of Gauss points and m denotes the weak form of the MIP method. The nodal operator matrix B d,g is the standard strain-displacement matrix for displacement-based elements, see Section 2.2.1. The constitutive strainŝ E g are given by (16), while E ,g is computed via (3) from the deformations alone. Due to the pointwise definition of the independent stresses, (26b) yields n independent equations (one at each Gauss point) of the form since S g is arbitrary. Inserting this relation into (26a) would lead to the pure displacement-based formulation (8). However, herein system (26) is used to increase robustness of displacement-based elements in the NR scheme, as demonstrated by Magisano et al. 32 To that end, the linearized form of (26) is needed. After considering that the variations e and S g are arbitrary this is given by with the tangential complianceD g , defined in (17). The geometric part of the tangent G d,g has already been introduced in (9c), but here it is evaluated from the independent stresses S g , which in general differ fromŜ ,g . This is the key to the MIP method derived below.
In the next step, static condensation of the independent stresses S g is performed in analogy to the derivations in Section 3. The last n rows of (28) yield withĈ g =Ĉ(Ê g ) =D −1 g . Inserting this result into the first row of (28) leads to which is rather similar to (8). It includes, however, a major disadvantage, which is the necessity of computingÊ g from the inverse stress-strain relation (e.g., Appendix A2.2). For linear stress-strain relations, i.e. St. Venant-Kirchhoff material, (29) and (30) can be greatly simplified. In this case,Ĉ g =Ĉ ,g andĈ g (E ,g −Ê g ) =Ŝ ,g − S g hold, which yields for the update procedure of the independent stresses at every Gauss point within iteration (j) of the global NR scheme. Moreover, (30) reduces accordingly to which is identical to (8), with the exception of S g being used for computing the geometric tangent instead ofŜ ,g . The right-hand side remains unchanged, which means that the equilibrium solutions of (32) and (8) coincide. Only the tangent is modified, which ultimately leads to improved robustness in the NR scheme, as demonstrated in the numerical simulations in Section 6. In contrast to the simple example from Section 3, the MIP method for solid finite elements is not equivalent to the standard HR formulation due to the point-wise, instead of elementwise, approximations of the stresses.
Equations (31) and (32) are the core of the MIP method proposed by Magisano et al., 32 who are taking linear material behavior into account. In case of general material models, we assume that (31) and (32) still hold, at least approximately. This implies that assumptionsĈ g ≈Ĉ ,g andĈ g (E ,g −Ê g ) ≈Ŝ ,g − S g (33) are made. As the right-hand side in (32) does not require inversion of the stress-strain relation, it is easily applied to general material models. There are three cases in which Assumptions (33) are fulfilled. First assumptions (33) are fulfilled in the first NR iteration of any load step since in equilibrium (27) holds and thus also S g =Ŝ ,g . Second, the assumptions are true for NR iterations close to equilibrium (i.e., in the range of quadratic convergence) since then ||Δ (j) || ≪ 1 and we get from (31) S g ≈Ŝ ,g respectively E ,g ≈Ê g . A final case, for which (33) (33) are studied in the numerical simulations presented in Section 6. There it is shown that the proposed MIP method (32) works best as long as the problem is limited to small strains. As soon as large strain occur, the method becomes less beneficial but still provides advantageous behavior in many cases.
Remark 8. The point-wise stress ansatz chosen in (28) and its condensed version (30) increase robustness, but do not allow the design of locking-free finite elements due to the limitation principle by Fraeijs de Veubeke 5 (cf. chapter 10.2 of Zienkiewicz et al. 42 ).

EAS finite elements
In this section, the MIP method is extended to EAS elements in order to overcome one of the major drawbacks of this highly popular element class, which is the lack of robustness in the NR scheme (cf. References 8,25).
In a first step, the weak forms of the EAS method (11) and the AS method (15) are combined by replacing the constitutive stressesŜ in (11) with the independent stresses S from (15). Additionally, (15b) is utilized to enforce the constitutive law. The resulting weak form of the three-field formulation is given by and marked with indexm. In a second step, we approximate the independent stresses S , as in (26), in a point-wise manner at the Gauss points and employ Gauss quadrature to get which is the linearized form of (34) on element level. The static condensation procedure, as presented in Section 4, leads in this case to for the update of the independent stresses at the Gauss points in step (j) of the NR scheme. This finally yields the linearized weak form of the MIP method for EAS elements given by n ∑ g=1 w g [ M e,g + G e,g (S g ) M e,g + G e,g (S g ) M e,g + G e,g (S g ) M e,g + G e,g (S g ) where the only difference compared to (13) is the use of the independent stresses S g instead of constitutive stresseŝ S ,g to compute the tangent. In analogy to Section 4, the right-hand side (i.e., the residual) is not altered. Thus, the converged solutions of (37) and (13) are identical. However, (37) improves robustness of the EAS method in NR iterations as shown in Section 6. Static condensation of the enhanced DOFs e can be performed on element level as usual. The extension of the MIP method for EAS elements to general material models is straightforward and analogous to Section 4. This implies especially that assumptions (33) are made in order to avoid the need for an inverse material law.

Spatial formulation of the MIP method
For many material laws, that are formulated on the basis of spatial measures, finite element formulations in the current configuration are advantageous due to higher numerical efficiency and simpler implementation (cf. e.g., Reference 7). Because of these benefits, in this section, the MIP method is transformed from the reference configuration (Sections 4 and 5.1) to its spatial form.
A push forward of the AS weak form (15), which is the basis for the displacement-based MIP method (see Section 4), yields where ∶= F S F T is the (independent) Kirchhoff stress tensor. Furthermore, e = 1 2 (I − F −T F −1 ) andê = F −TÊ F −1 denote the displacement-based and constitutive Almansi strain tensor, respectively. However, in the approximated regime two ways to introduce the pointwise stresses at the Gauss points emerge. One could, as in previous sections, approximate the material second Piola-Kirchhoff stress S g or instead define new independent stresses g ∶= F S g F T and use these as unknowns. Our numerical experiments have shown that the latter approach is inferior to the pointwise approximation of S, which is why only this approach is considered here.
Proceeding similarly with the linearized form of (38) and performing the steps presented in Section 4 ultimately yields the MIP weak form and the stress update formula in the current configuration given by where b d,g are the spatial nodal operator matrices andĉ denotes the spatial material tangent, which is given in index notation byĉ abcd = F aA F bB F cC F dDĈABCD . Moreover, performs the transformation g = F S g F T in vector-matrix form. All in all, the only necessary changes to get from material to spatial MIP form is exchanging all material variables in (31) and (32) with their spatial counterparts and perform the transformation using (41). We emphasize that (39) and (40) are equivalent to (32) and (31), respectively. A spatial version is also easily obtained for the MIP method for EAS elements. The same steps lead to a similar result, where all quantities in the reference configuration in (37) and (36) have to be replaced by corresponding spatial variables. For the sake of brevity, we omit corresponding derivations.

Implementation
The implementation of the MIP method described in the previous two sections is relatively simple. Only a few changes have to be made in existing codes. A possible global code structure in pseudo-code format is given in Algorithm 1 for the EAS element with MIP extension. The corresponding element routine is shown in Algorithm 2. In those algorithms, all changes and additional operations necessary for the MIP method are marked in red. The Algorithms 1 and 2 are generally applicable regardless of material model or solution routine. For instance, they can be applied for elasto-plastic simulations even if a line-search modified NR scheme is used (see Section 6.6).
end if compute residual and tangent: R e,g (Ŝ ,g ); R e,g (Ŝ ,g ); M e,g + G e,g (S g ); M e,g + G e,g (S g ); M e,g + G e,g (S g ); M e,g + G e,g (S g ); save info for update of S g :Ŝ ,g ; (Ĉ ,g B e,g ); (Ĉ ,g B e,g ); end for static condensation of e end for

One-element test
A one-element test is presented in this section to demonstrate the effects of the MIP method on the stresses and convergence properties of solid elements. In analogy to the simple example in Section 3, the aim is to highlight basic effects of the MIP method for solid problems and reveal differences compared to standard methods. More complex numerical investigations follow in Section 6. The problem setup of the one-element test investigated here is illustrated in Figure 4. The free end of the cantilever beam structure is subjected to a load F = (Et 3 )/(L 3 ) and a St. Venant-Kirchhoff material model under plane strain conditions is used (see Appendix A1). The evaluations presented in the following are carried out on Gauss point level. In particular, the spatial Kirchhoff -stresses are examined at the lower left Gauss point, as indicated in Figure 4.
In a first step, the AS element Q1/S5 3 is compared to the standard EAS element Q1/E4 1,10 (see Section 6 for more detailed descriptions of the elements). These two elements are free of locking for the present undistorted case and are equivalent in linear problems (e.g., Bischoff et al. 6 ). In the geometrically nonlinear case, the elements are no longer identical, but still can be expected to provide similar results. Figure 5 shows the development of the axial stress component in x 1 -direction, 11 , for different values of Poisson's ratio . The stress 11 is normalized with the converged stress after all five load steps denoted by end 11 . The stresses plotted are determined at the beginning of every iteration, which explains why the stress in the first NR iteration of a load step is the same as the previously converged result. On the horizontal axis, load steps with converged results are marked at integers with vertical lines. Every NR iteration in between is marked with an "x". For example, the third Newton iteration out of a total of five necessary iterations in the second load step is marked at 2 + 3/5 − 1 = 1.6. For = 0, Figure 5 (left) shows that the stresses converge much faster for Q1/S5, which are almost converged after the first NR iteration within each load increment. This is in line with the observations of Magisano et al. 29,32 and the 2-DOF example presented in Section 3. The converged results of Q1/E4 are visually identical. However, after the first NR iteration the stresses are far worse, which results in an additionally required iteration in load steps 2-5. This difference becomes worse for the nearly incompressible case with = 0.499 shown in Figure 5 (right). In this case, the stresses partly have the wrong sign and are off by a factor of ≈36 after the first iteration in the first load step. This results in one to two additionally required iterations per load increment, while Q1/S5 is insensitive to the critical parameter . The MIP method allows to overcome this drawback which is illustrated in Figure 6. Therein, standard EAS element Q1/E4 is compared to Q1/E4-MIP, which uses the MIP method as described in Section 5.1. The MIP stresses g,11 at the Gauss-point (36) show almost identical behavior as the stresses of Q1/S5, while the constitutive stresseŝ1 1 (dashed line) used for the residuum are still far off in the first iterations. Using the altered g,11 for the geometric tangent enables the improved convergence behavior of elements. Furthermore note that the two kinds of stresses converge quickly to each other and are equivalent in case of convergence. The Q1/E4-MIP element requires the same number of NR iterations as AS element Q1/S5.

Overview
This section covers various numerical investigations examining the properties of the MIP method presented in Sections 4 and 5. The main concern of all investigations is robustness of the simulations, where in the context of the present work we characterize robustness by the size of applicable load steps and number of NR iterations needed to find a solution as well as the sensitivity w.r.t. to critical parameters, for example, the elements' aspect ratio or Poisson's ratio. In this sense, robustness implies also efficiency, since fewer load steps with fewer NR iterations yield the desired result. It will be shown that the MIP method improves robustness in many situations but also that there exist scenarios in which little or no positive effects can be observed. Other classical topics of element technology such as locking and convergence with mesh refinement are not covered here. All simulations are performed using 2D plane strain quadrilaterals or 3D brick elements, that is, hexahedra. The considered element formulations are: • Qp: Isoparametric Lagrangian displacement-based, quadrilateral with polynomial degree p. Corresponding hexahedra are denoted as Hp.
• Q1/E4: EAS element proposed by Simo and Armero, 10 which employs the four Wilson modes. In 3D the element is labeled H1/E9 and has nine enhanced modes.
• HA1/E12T: The 3D version of QA1/E4T. In addition to the modified quadrature rule it uses three additional enhanced modes (compared to H1/E9) and employs a special evaluation of the compatible deformation gradient.
• Q1/S5: AS element with a five-parameter stress interpolation, as proposed for linear kinematics by Pian and Sumihara. 3 The extension to nonlinear kinematics can, for example, be found in Viebahn et al. 8 In 3D, the element is denoted H1/S18 and has 18 stress modes (cf. for example, References 4,8,43). • Q1/E4-MIP*: Element formulation as defined above but without considering the simplifying assumptions (33).
The material models used throughout this Section are a St. Venant-Kirchoff, Neo-Hooke and logarithmic strain-based elasto-plastic model, which are described more detailed in Appendix A.
As mentioned above, robustness is the major concern of the present work which is highly influenced by settings chosen for the NR procedure. For all simulations we use a convergence criterion based on the residual norm in the form ||R|| < NR . Note that no scaling of this norm is applied as, for example, suggested by Belytschko et al. 44 If not mentioned otherwise, tolerance NR is set to 10 −8 for all simulations. Failure of the NR procedure is determined if either ||R|| > 10 14 or more than 20 iterations are necessary to find a solution within one load step.

Clamped beam
The setup of the first numerical example is shown in Figure 7 and consists of a rectangular block under plane strain conditions, meshed with 1 × 10 elements. Two elastic material models are considered for this example. The left edge is fixed horizontally and the vertical displacement is fixed only at one single point, to avoid artificial constraint stresses in y-direction. The Singularity at this point has no influence for the coarse mesh considered here. On the right-hand side the structure is subjected to a uniformly distributed vertical force resultant which is scaled such that the vertical tip displacement of a corresponding geometrically linear infinitely wide plate is always u lin = 4.0. In case of = 0 it is furthermore identical to a thin beam solution (Bernoulli). Value u lin clearly deviates from the solutions for a thick beam or when a geometrically nonlinear setting is considered. However, it is a way to keep the deformation (degree of nonlinearity) in a similar range for varying problem parameters, which is useful for the subsequent investigations. The displacement-based elements that will be considered in the following are the bi-quadratic Q2 element and its MIP version. They are chosen instead of the Q1 elements, since they are less susceptible to locking and thus better comparable to the other elements tested. These are the standard Q1/E4 and its MIP version and finally Q1/S5. The latter always performs most robust for this problem setup (this is not the case in general, see Section 6.3) and thus its number of necessary NR iterations is considered as "target" for all other elements.  First, the influence of critical geometric parameters, such as the beam's slenderness L/t are investigated. In order to separate different effects, Poisson's ratio is set to zero and only the thickness is varied (t = {1, 0.5, 0.2, 0.1, 0.05, 0.02, 0.01}). The load F is applied within one single load increment. The diagrams in Figure 8 show the cumulative number of iterations n NR that are required to fulfill the convergence criterion, plotted versus the slenderness. For both materials, it can be observed that for Q2 and Q1/E4 the number of required iterations grows or the NR scheme even diverges with an increasing slenderness. This indicates that the parameter dependency of the robustness is not a unique phenomenon of EAS elements but rather an intrinsic behavior of constitutive-based stress updates. On the other hand, the MIP versions of both the Q2 and Q1/E4 elements exhibit robust behavior, independent of the slenderness, with only minor differences compared to Q1/S5. This test confirms the results of Reference 32, where an improved robustness was observed for the numerical analysis of thin-walled structures using solid-shells. For the case t = 0.05, a detailed summary of the residual norm as well as the converged tip displacement û is shown in Table 2. Whether MIP is applied or not does not affect the converged displacements, since the residual is not modified (cf. Section 4).  Figure 9 shows the required cumulative number of NR iterations n NR over the bulk modulus when F is applied within one (top) and five (bottom) load increment(s), respectively. For the St. Venant-Kirchhoff material a similar behavior as before is observed. All MIP elements and the Q1/S5 perform robustly, whereas Q2 and Q1/E4 show, for high bulk moduli, an increasingly less robust convergence behavior. However, this outstanding performance of the MIP elements cannot be reproduced completely when a Neo-Hookean material is used. Although the average response of the MIP formulations is still superior to Q2 and Q1/E4, convergence is not achieved in all cases. The reason for this bad performance can be explained by the stronger violation of assumptions (33). First, because larger strains occur compared to the thin beam case and second, the nonlinearity of the NH material directly depends on degree of incompressibility (cf. term 2 ln 2 J in (A3)). However, for Q1/E4-MIP*, which is the "correct" MIP formulation based on (35) (i.e., the simplifying assumptions (33) are not made), the same results as for Q1/S5 are obtained. Unfortunately though, the inverse stress-strain relation is needed again, which destroys the major advantage of the MIP method.

Elastic strip
The next example is the elastic strip test proposed by Korelc and Wriggers 12 (see also Reference 45), which is one of the few examples where robustness of EAS elements has been examined. Its initial square geometry (a = 10) and the deformed state are shown in Figure 10 together with the boundary conditions. The strip is loaded by prescribed displacements u p = 10 applied on the right edge. A regular FE-mesh with 10 × 10 elements is used for all simulations. As material model we consider the Neo-Hookean model described in Appendix A and chose the elasticity constants to = 24 and = 6. The St. Venant-Kirchhoff material is not considered in the following for this example, since it leads to unphysical results with artificial boundary layers along the free edges as shown in Figure 10.
Results of the elastic strip example in 2D are shown in Table 3. In contrast to the clamped beam example in Section 6.2 there are only minor improvements in the total number of NR iterations n NR for Q1/E4 if the MIP method is used. For other elements, namely Q1, Q2, and Q1/E4T, even more iterations are required if the MIP method is applied. This behavior follows from the very high strains occurring in this example. For these states the Neo-Hookean and St. Venant-Kirchhoff model differ noticeably which implies that assumptions (33) are strongly violated. Thus, the MIP approach is, for this kind of problems, a poorer approximation of an AS method and therefore less efficient. Furthermore, the MIP technique is more favorable in bending dominated problems instead of the uniaxial problem considered here. In general, it can be observed that the MIP method looses efficiency in case of large strains and general material models. However, in less extreme cases than the example considered here, there are usually still improvements due to the MIP strategy (see subsequent examples).

TA B L E 3 Results of the elastic strip example
Very interesting are the results obtained with Q1/S5. This element employs the inversion of the stress-strain relation of the Neo-Hookean material model presented in Appendix A2.2. While it works well and converges fast for smaller displacements, failure of the simulations can be observed at u p ≈ 5.36. At that level of deformation the inverse stress-strain relation looses its uniqueness (see Section A.2.2 Remark 9) starting with elements close to the necking zone. Thus, the NR routine aborts at that point as a direct consequence of failure of the material routine. This behavior is almost independent of the number of load steps and mesh refinement. Furthermore, it is not an artifact of intermediate states during the NR procedure, since using the very high number of n steps = 6400 load steps, which leads to almost immediate convergence with n NR ≤ 3 in every step, fails as well. Even the numerical procedure to invert the stress-strain law proposed by Viebahn et al. 8 finds no solution in those states. All in all, it can be concluded this is an actual mathematical problem of the inverse stress-strain relation of the Neo-Hookean material and not related to the numerical procedures.
All results presented for the 2D case can also be qualitatively observed in 3D simulations.

Thin circular ring
The first 3D example in this work is the thin circular ring shown in Figure 11, which was introduced for the analysis of shells by Reference 46 (see also References 32,47). In the present work, we consider the example for solid elements as described by Korelc et al. 13 The ring shown in Figure 11   apply. To complete the setup, a surface dead load q = 6.67 ⋅ 10 −3 is applied in z-direction on face F 2 and the elasticity constants are chosen to = 10,500 and = 0. Results of this numerical example are summarized in Table 4 for a selected set of elements. Note, that H1 and H1-MIP are excluded from Table 4 due to severe locking, which yields underestimated displacements and thus also few necessary NR iterations (cf. Section 6.2). For all other elements, the MIP method greatly improves convergence of the NR procedure. It enables to apply the complete load for all EAS elements within one step instead of 3 to 5 without MIP. Furthermore, the number of total necessary NR-iterations n NR is reduced by a factor of at least 5.9 for the St. Venant-Kirchoff material and 4.7 for the Neo-Hookean model, with the greatest improvement observed for element HA1/E12T. As demonstrated before, AS element H1/S18 shows again superior robustness for both material models and is on top of that closest to the converged result of w P = 10.26.

Spherical shell with opening
The second 3D example is the spherical shell problem with opening shown in Figure 12  The required numbers of load steps n steps and NR iteration n NR as well as the reaction force in z-direction R z are shown in Table 5 for both the St. Venant-Kirchhoff and Neo-Hookean material model.
In a first step we examine the results for t = 0.5. For this setup, all elements exhibit reaction forces within a 2% margin of the converged result for both materials. More interesting are the required number of load steps and NR iterations. It can be observed for the St. Venant-Kirchhoff material, that, even though the number of load steps differs, almost the same amount of total NR iterations is needed regardless of which element is used. However, for the Neo-Hookean material larger differences can be observed. For this material model the AS element H1/S18 excels by still requiring only a total of n NR = 27 NR iterations. All other elements need at least twice the amount of iterations. Better results are obtained with the MIP-modified versions, which require approximately 15% less NR steps compared to the standard EAS elements. For the thin thickness t = 0.05 more load steps and iterations are necessary for all elements and the reaction forces differ substantially. However, improvements due to the MIP method are more pronounced, which is in line with the results in Section 6.3. Since smaller strains occur for the thinner shell, (33) is less severely violated which leads to the improved behavior of the MIP method. All in all, this example confirms that the MIP does not always yield such high improvements as shown in Sections 6.2 and 6.4, depending on the magnitude of the strains. However, it usually improves robustness of strain driven elements if strains are not too high.

6.6
Elasto-plastic circular bar The final example of the present work is an elasto-plastic simulation based on the material model proposed by Simo,51  The test covers necking of a circular bar with radius R = 6.413 and a total length of 2L = 53.334 (e.g., References 10,11,13,16,22,48,51,[53][54][55]. Due to symmetry, only one-eighth of the bar has to be considered, which is shown in Figure 13. To initiate necking, the radius is linearly reduced from R to R = R − 0.07 along the length of the bar. Two quarter cylinders with 480 elements each are used to mesh the specimen such that the lower fifth of the bar is refined (see Figure 13). Symmetry boundary conditions u i (X i = 0, X j , X k ), i, j, k ∈ {1, 2, 3} apply and the structure is loaded by prescribed displacements u(X = L, Y , Z) =ū. During the simulationsū is gradually increased toū = 7.0, where the first half of load steps is used untilū = 5.6 and the second half covers the remaining Δū = 1.4. This is a standard procedure (cf. References 10,22) since the final steps of this simulation are especially demanding due to the softening of the material. Usually, a line-search (LS) algorithm is used to stabilize the NR procedure in elasto-plastic simulations (e.g., References 10,22,48,51,53). In the present work we use the method described in Bonet and Wood. 38 The total number of NR iterations is shown in Table 6 for element H1/E9 and various numbers of load steps. Without MIP method and line search (LS) the simulation requires at least 28 load steps and is volatile in a sense that it is not guaranteed that a higher number of load steps leads to a solution. The robustness can be improved with either the MIP method or a line search algorithm. Both reduce the number of necessary load steps and ensure (with one exemption) that the NR method converges for higher number of load steps. The best results are obtained by combination of both methods. This allows the lowest number of load steps and also needs slightly fewer iterations than the line search approach without MIP method.

CONCLUSION
The present work covers extensions to the MIP method, which has recently been proposed by Magisano et al. 32 This method has so far been used to improve numerical robustness (meaning number of required NR iterations and size of applicable load steps) of displacement-based and mainly structural finite elements using the St. Venant-Kirchhoff material. Herein, we considered solid finite elements and proposed three extensions.
First of all, we suggested a simple extension to general material models. As for the standard MIP method, we do not interfere with the residuum of the nonlinear FE problem and only modify the tangent, which means that the converged solution of standard FE simulations and the MIP-modified versions are identical. Only the robustness of the MIP-modified elements is increased as we showed in many numerical simulations. Second, a MIP version using a spatial formulation of FE equations was introduced in order to simplify implementation of more complex material models such as the elasto-plastic model considered in this work. Finally, we extended the method to EAS finite elements and were able to overcome one of the major drawbacks of these elements in many applications. In particular, we showed that the proposed MIP approach cures their lack of robustness when compared to AS elements. With the simple modifications necessary for the MIP method, EAS elements become much more efficient and robust, which makes them even more interesting for practical simulations.
On top of that, we proposed a novel inverse stress-strain relation for a Neo-Hookean material, which had to the best knowledge of the authors not been proposed before. With this relation we were able to show that AS elements exhibit superior robustness compared to classical strain driven element formulations. Unfortunately, in accordance with literature, we also showed that this approach is not generally applicable in all states of strain. In fact, this observation provides another argument in favor of using much simpler strain-driven finite elements together with the MIP modification. This approach allows the benefit of both, favorable robustness and simple implementation of complex material models.
However, there are still a few open issues. While the MIP method greatly increases robustness of finite elements for general material models and small strains, less improvement or in some special cases even slightly disadvantageous behavior can be observed if large strains occur. A thorough explanation and cure of this issue should be the goal of further investigations. Another line of research could follow creating a MIP approach based on a Hu-Washizu functional instead of the Hellinger-Reissner approach followed thus far. form given by where and are the Lamé constants. Note furthermore that (A2) is straightforward to invert by simply computing the inverse of constant C SVK in vector-matrix form. Thus, the St. Venant-Kirchhoff material can easily be used for AS-elements introduced in Section 2.2.3.

A.2. Neo-Hooke A.2.1 Standard form
The second material law considered in this work is a Neo-Hooke law with strain-energy function where > 0 and > −2∕3 are the Lamé constants and J = det F. Derivation of this definition with respect to C yields the second Piola-Kirchhoff stresses

A.2.2 Inverse Neo-Hookean law
The inverse relation of (A4) presented here has to the best knowledge of the authors never been proposed before. Only a similar inverse relation for a different Neo-Hookean model is given in Wriggers. 7 Similar to the derivations there several cases have to be considered which are laid out subsequently.

Case = 0
In this case, which corresponds to = 0, it is straightforward to obtain as inverse stress-strain relation from (A4). Thus, I − −1 S must be invertible which can be examined with relations det(I + A) = 1 + I 1 (A) + I 2 (A) and det(I + A) = 1 + I 1 (A) + I 2 (A) + I 3 (A) holding for an arbitrary tensor A with Invariants I i (A) in 2D and 3D, respectively (e.g., Reference 56). From these relations, condition det(I − −1 S) ≠ 0 for invertibility and requirement J 2 = det(C) > 0 necessary for physically meaningful results, we get that the eigenvalues S i of S must fulfill either S 1 , S 2 > or S 1 , S 2 < in the 2D case. For 3D problems either restrictions S i , S j > , S k < or S i , S j , S k < where i, j, k are permutations of {1, 2, 3} apply. See also chapter 6.2.2 of Ogden 9 for similar results.
Case ≠ 0 In this case simple rearranging of (A4) yields where auxiliary variables have been introduced. Tensor A has to be invertible which implies that the eigenvalues of S have to fulfill the same requirements as in the case = 0. The inverse stress-strain relation of (A4) is given by (A6) 2 where the only unknown is (J).
In case of a = 0 (A6) 1 yields = 0 since C has to be invertible for physically meaningful results. Thus only A = 0 ⇔ S = I would be allowed in order to fulfill (A6) 1 . In that case, however, C is not uniquely defined and thus no inversion of the stress-strain relation is possible.
If a ≠ 0, the next step is to take the determinant of (A6) 1 . This yields for both the 2D plane strain and 3D case an equation for J given by where d ∈ {2, 3} is the spatial dimension. Lengthy computations using the Lambert-W function (x) finally yield where requirements have to be met. Note that always branch  0 of the Lambert-W function is needed, which follows from the solution for a > 0 and continuity requirements. Furthermore, (A8) automatically ensures J > 0 which allows using the logarithm in (A9).
Remark 9. Note that regardless of the case = 0 or ≠ 0 there are restrictions on which stress tensors S allow a unique inversion of the stress-strain relation. These states actually occur in practical simulations with large strains as shown in Section 6.3. In that case these states lead to failure of the computation. This is not a problem of the numerical procedure but results form the physical equations, that is, (A4) as shown in Section 6.3.

A.2.3 Logarithmic strain-based von Mises elasto-plasticity
The final material model considered in the present work is the elasto-plastic model proposed by Simo. 51 This eigenvalue-based formulation is widely used in the context of finite element development (e.g., References 11,18,22,25) and based on the standard multiplicative split F = F e F p into elastic and plastic parts. 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 nonlinear isotropic hardening with saturation and the associative flow rule. More information on the material model and algorithms for standard elements are given in the work of Simo. 51 For the numerical implementation of the model an eigenvalue perturbation technique according to Miehe 57 is applied to avoid treatment of duplicate stretches.

APPENDIX B. TWO DOF EXAMPLE
This appendix gives additional details on the simple 2 DOF example presented in Section 3. Focus is put on the assumptions and derivations needed to obtain (19) and (20). First of all, the kinematic relation in (18) follows from the assumption of linear displacements u and v along the axis of the bar. This implies that the Green-Lagrange strain is constant and given by where L = 1 and l are the original and deformed length of the bar, respectively. From the assumption of a linear relation between Green-Lagrange strains and the normal force N (i.e., St. Venant-Kirchhoff material) we get the constitutive relationN where A is the constant reference cross section of the bar and E denotes its Young's modulus. The variational functional associated with the 2 DOF example can be cast in the form Therein, k is the stiffness of the linear spring and is the external force. The internal potential of the bar Π bar is given by in the displacement-based and AS case, respectively. For the AS formulation an independent stress field N is introduced. Furthermore, the Legendre-transformation of the internal energy W( ) =N given by U( ) = N̂− W(̂) has been used. Therein,̂are strains obtained from the inverse stress strain relation. Variation of the internal part Π bar of functional Π 2DOF yields for the displacement-based and AS formulation. Variations (B5) are closely related to the continuum formulations (7) and (15). In fact in case of the bar we have a uniaxial stress state, where the only nonzero stress component is S 11 , and furthermore constant stress and strain along the bar. Imposing these restrictions on (7) and (15) directly yields (B5). Note that this close relation allows the transfer of the results concerning robustness from the simple bar problem to the continuum formulation. Imposing the stationary condition Π 2DOF = 0 and subsequent linearization yields the residuum and tangent of both the displacement-based and AS form of the 2 DOF example given in (19) and (20), respectively.