The Nonlinear Finite Element Method is a widely used numerical technique to solve engineering problems associated with nonlinear solid continua. Often, real-life engineering problems require immense computational resources. Thus, raising the computational efficiency of the method is a highly desirable goal. A nonlinear global system solver requires constant reevaluation of the element-level internal forces to- gether with the nonlinear stiffness matrix. For this purpose, commercial software employs the standard numerical integration based on the quadrature method (e.g., Gauss points). Consequently, computational complexity increases linearly with the number of integration points. Therefore, the derivation of alternative, highly effi- cient integration approaches is essential, and this is the main goal of our study. Herein, we propose new element-level integration formulae of the internal forces and of the nonlinear tangent stiffness matrix. Our formulae admit the “full” order integration requirement. Moreover, we demonstrate analytically that the computa- tional cost of the proposed schemes is roughly equivalent to One-Point quadrature, irrespective of the element type (e.g., tetrahedral, hexahedral, wedge, etc.,) and ir- respective of the element order (e.g., 8-node brick, 20-node brick, etc.,). Code im- plementation of our formula follows a rather familiar, standard, manner. However, prior to code implementation, it requires sets of coefficients/special weights to be pre-computed or adopted from the literature. Undoubtedly, the proposed integra- tors require significantly more coefficients than the standard numerical integration. Recently, the Semi-Analytical approach has been employed to produce highly effi- cient case-specific integrators for the mass matrices. Here, we generalize and expand this approach to all element-level integrals. To this end, the integrands are decom- posed into two multiplicative parts. The first part includes kinetic and kinematic functions which admit “full” integration criteria for one sampling point (e.g., the centroid). While the second part consists of mesh-independent and displacement- independent polynomial functions. We integrate the above polynomials analytically, to derive coefficient sets. Those coefficients are incorporated in the resulting scheme’s subroutine. In other words, we take advantage of the mathematical structure of the integrand to produce a highly efficient yet case-specific integration formula (e.g., our internal forces rule can’t be used for the mass of stiffness matrices and vice versa). Importantly, code implementation doesn’t involve either meta-programming or a computer algebra system for explicit (closed-form) code generation. Key words: Finite Element Technology; Semi-Analytical Integration; Full Inte- gration; Nonlinear Tangent Stiffness Matrix; Closed-Form; Element-Level Integra- tion; Corresponding author(s): Eli Hanukah, Mechanical Engineering Department, Braude College of Engineering, eli.hanukah@braude.ac.il1 2 1Introduction Nonlinear Finite Element Method (NFEM) is the most widely used numerical tool to solve engineering problems associated with nonlinear continua. Often, practical problems require extensive computational resources. Therefore, raising the compu- tational efficiency of the overall procedure is highly desirable. Fundamentally, NFEM requires a nonlinear solver (e.g., Newton/modified Newton- Raphson, etc.) to achieve a global system solution. Such a solver is required to update the internal and external forces together with the nonlinear stiffness matrix every com- putational step i.e., the solution procedure goes back and forth between element-level computations and global system iterations. All element-level evaluations, such as the mass matrix, internal forces, and the nonlinear tangent stiffness matrix, require in- tegration in the element domain. Most commonly, those are performed using the Standard (ST) numerical integration, based on a specific quadrature scheme. Thus, computational complexity grows linearly with the number of integration points, which increases dramatically the element formation cost. The goal of our study is to radically improve element-level integration efficiency for the solid NFEM. We achieve this goal by proposing new integration formulae, based on the Semi-Analytical approach. The proposed method is sufficiently general i.e., it addresses all isoparametric solid ele- ments of general order, including both, geometrical and material non-linearity. It is known that increasing the efficiency of element-level integration yields signif- icant CPU time savings [ 1]. Noticeable savings can be achieved by ”Closed-Form” code generation using Computer Algebra Systems CAS [ 28,1]. For the linear elastic- ity and regular shape elements (straight-sided), a combination of analytical integration, by means of CAS, together with closed-form code generation, yields dramatic im- provements in CPU time (e.g., [ 29,30]). In addition, special structural formulations based on Cosserat Point Element (CPE) theory, technically do not require integration (e.g., [22,21]). Yet, to the best of our knowledge, for a fully nonlinear FE, there is no alternative systematic approach to dramatically over-perform the standard numerical integration. Namely, there is no alternative method, which can be consistently ap- plied to all solid elements having a general constitutive behavior i.e., whether 1D, 2D, or 3D elements, arbitrary shape & order (e.g., quadrilateral, triangular, bricks, tetra- hedral, wedges, etc.,). This is the gap we attempt to close in the present contribution. Recently, the Semi-Analytical (SA) approach has been employed to derive highly efficient integration formulas for the consistent/lumped mass matrices [ 15,17,18]. We generalize this method and specialize it for the nonlinear tangent stiffness matrix and the internal forces. The resulting integration formulae are case-specific, i.e., a unique formula for the mass matrix, an additional unique formula for the internal forces, and the last one for the nonlinear tangent stiffness matrix. Yet, the result- ing schemes are sufficiently general, in the sense that the same formula addresses all isoparametric solid elements of an arbitrary order (whether a wedge, tetrahedral, hex- ahedral, etc.,). SA integration presented here, admits the “full” order integration requirement. Specifically, it is exact for a regular shape and a constant strain element, while ap- proximation for all the other cases. Importantly, computation-wise, the SA formulae are roughly equivalent to the Standard (ST) One-Point quadrature, regardless of the 3 element order! Thus, the higher the element order, the more pronounced the effi- ciency advantages of our schemes. Throughout the study, we explicitly consider the isoparametric element and a Total Lagrangian formulation. Both of those choices are widely used and accepted. There- fore, they serve as a good framework to demonstrate the main ideas of the contribu- tion. However, the SA method is not restricted by those choices. Incorporation of the Updated Lagrangian formulation follows the same guidelines. Subroutine implementation of the proposed integration rules, do not require spe- cial programming strategies such as meta-programming or the use of CAS for closed- form code generation. The resulting formulae take a familiar form, their code imple- mentation follows conventional tactics. Organization of the paper is as follows: section ( 2) recalls the essential well-established element-level kinematic and kinetic relations and formulae, for one to be able to follow through later derivations. Then, subsection ( 2.1) recalls the Standard (ST) numerical integration approach. Next, subsection ( 2.2) recalls and discusses the ”full” order vs. the ”reduced” order integration concepts. Section ( 3) begins with the general concept of the Semi-Analytical approach; it presents the guidelines that we follow to derive the proposed formulae. In addition, computational equivalence between the full integra- tion SA scheme and the ST One-Point quadrature is briefly discussed. Subsections (3.1),(3.2) and (3.3) contain derivation of the SA integration schemes for the mass ma- trix, internal forces and the nonlinear tangent stiffness matrix respectively. Finally, section ( 4) discusses, highlights, and summarises the important points and takeaways of the study. 2Background In this section we recall the basics of the well-established Total Lagrangian (TL) formulation. In particular, we end up with expressions for the consistent mass matrix, internal forces and the nonlinear tangent stiffness matrix. Those require numerical integration. Subsection ( 2.1) emphasizes the standard quadrature-based numerical integration. Then, subsection ( 2.2) recalls the ”full” and ”reduced” order integra- tion approaches. Even though we confine ourselves to isoparametric TL formulation throughout the section and the study, the proposed SA approach is not limited by those constraints. We keep this section extremely brief since all of the below are well- established relations that can be found in at least several great books e.g. [ 42]. Within the framework of the isoparametric concept, the displacement field uis ap- proximated in the element domain by u≈ue=∑n I=1NI(ξ)uI(t)where ndenotes the number of element nodes, ξ= (ξ, η, ζ )stand for the local coordinates, NI(ξ), (I= 1, .., n )are the shape functions, uIrepresent the nodal displacements, while Ωe refer to the initial element configuration. Position vectors in the initial and current configurations are denoted by Xandx. They are interpolated in the element domain using the same shape functions Xe=n∑ I=1NI(ξ)XI,xe=n∑ I=1NI(ξ)xI(t) (1) 4 Where XIandxIstand for respective nodal positions. Shape functions NIusu- ally take a rather simple form, most frequently a polynomial interpolation. For example, for a standard twenty-node hexahedral element, the first shape function is given by N1(ξ) =−1 8((1−ξ) (1−η) (1−ζ) (2 + ξ+η+ζ)). For later conve- nience, we emphasize that the shape functions are mesh & displacement-independent NI̸=NI(XJ),NI̸=NI(uJ(t))One can easily analytically differentiate those func- tions with respect to local coordinates. Also, those functions and their derivatives can easily be analytically integrated in the parent element domain (parametric space) to result in a number. Further down the line, we’ll use those numbers as coefficients. It is important to distinguish between differentiation of NI(ξ)with respect to local and global coordinates ∇ξNI(ξ) =  NI;ξ(ξ) NI;η(ξ) NI;ζ(ξ)  =  NI;1 NI;2 NI;3  ,(I= 1, .., n ) (2) where semicolon is used to denote partial differentiation with respect to local coordi- nates e.g. NI;3=∂NI ∂ζ. Similarly, gradients with respect to initial global coordinates X1, X2, X3or current global coordinates x1, x2, x3are given by ∇XNI=  NI,1 NI,2 NI,3  ,∇xNI=  NI,1 NI,2 NI,3  ,(I= 1, .., n ) (3) where comma represents partial differentiation with respect to global coordinates. The transformation between the gradients of different configurations is given by ∇XNI=J−T e∇ξNI,∇xNI=j−T e∇ξNI,(I= 1, .., n ) (4) where with the help of ( 1), entries of JACOBI matrices Je=Grad ξXe=∂Xe ∂ξand je=Grad ξxe=∂xe ∂ξtake the next form [Je] = X1;ξX1;ηX1;ζ X2;ξX2;ηX2;ζ X3;ξX3;ηX3;ζ , X m;k=n∑ I=1NI;kXm,I ,(m, k = 1,2,3)(5) [je] = x1;ξx1;ηx1;ζ x2;ξx2;ηx2;ζ x3;ξx3;ηx3;ζ , x m;k=n∑ I=1NI;kxm,I,(m, k = 1,2,3)(6) Relations ( 4) require inverse transform of the above 3×3matrices ( 5),(6). Inverse of an arbitrary real square 3×3matrix Ais given by A−1=1 det(A)adj(A), where determinant and adjoint matrix follow simple rule det(A) =a1,1a2,2a3,3−a1,1a2,3a3,2−a1,2a2,1a3,3+ a1,2a2,3a3,1+a1,3a2,1a3,2−a1,3a2,2a3,1 Adj(A) = a2,2a3,3−a2,3a3,2−a1,2a3,3+a1,3a3,2a1,2a2,3−a1,3a2,2 −a2,1a3,3+a2,3a3,1a1,1a3,3−a1,3a3,1−a1,1a2,3+a1,3a2,1 a2,1a3,2−a2,2a3,1−a1,1a3,2+a1,2a3,1a1,1a2,2−a1,2a2,1 (7) 5 thus we can write [ J−1 e] =1 JeAdj(Je),[ J−T e] =1 JeAdjT(Je),[ j−T e] =1 jeAdjT(je)(8) Je=det(Je), j e=det(je) (9) In terms of ( 5),(6),(8),(9) deformation gradient Feis given by Fe=jeJ−1 e, JF e=det(Fe) =je Je(10) The nonlinear strain-displacement 6×3matrix BI,(I= 1, .., n )is given by BI= F1,1NI,1 F2,1NI,1 F3,1NI,1 F1,2NI,2 F2,2NI,2 F3,2NI,2 F1,3NI,3 F2,3NI,3 F3,3NI,3 F1,1NI,2+F1,2NI,1F2,1NI,2+F2,2NI,1F3,1NI,2+F3,2NI,1 F1,2NI,3+F1,3NI,2F2,2NI,3+F2,3NI,2F3,2NI,3+F3,3NI,2 F1,1NI,3+F1,3NI,1F2,1NI,3+F2,3NI,1F3,1NI,3+F3,3NI,1 (11) Importantly, in the above, shape functions are differentiated with respect to global X1, X2, X3coordinates. Next, GIJ,(I, J = 1, .., n )matrix is recalled GIJ= (∇XNI)TSe(∇XNJ) (12) GIJ= [NI,1, NI,2, NI,3]T S11S12S13 S21S22S23 S31S32S33   NJ,1 NJ,2 NJ,3  , (13) where Seis the Second Piola-Kirchhoff stress tensor. Following ( 12) it becomes ap- parent that differentiation in the above ( 13) is done with respect to global coordinates X1, X2, X3. Finally, we recall the expressions, for efficient integration of which, the present work exists i.e. - the consistent mass matrix, internal forces, and the nonlinear tan- gent stiffness matrix. External forces are omitted since their SA treatment follows identical guidelines proposed below and is quite immediate. The consistent mass matrix is given by MIJ=∫ Ωeρ0NINJdΩe=∫ Ω□ρ0NI(ξ)NJ(ξ)Je(ξ;XK)d□ (I, J, K = 1, .., n )(14) where ρ0denotes the initial density. Since the mass matrix does not depend on the displacement, it needs to be evaluated only once. In that sense, even though its effi- cient evaluation is desirable, it is hardly the ”bottleneck” of the overall procedure. On the contrary, a nonlinear solver (such as Newton/Newton-Raphson) requires the in- ternal forces and the nonlinear stiffness matrix to be evaluated every step of the global system solution. This, in turn, makes their integration efficiency improvements ex- tremely sought after! The internal forces RIare given by RI=∫ ΩeBT ISedΩe,(I= 1, .., n ) (15) 6 whereas here, Seis the Voigt notation of the Second Piola-Kirchhoff stress tensor S that is Se={S11, S22, S33, S12, S23, S13}T. The nonlinear tangent stiffness matrix KT IJis given by KT IJ =∫ Ωe[ GIJI+BT IDB J] dΩe,(I, J = 1, .., n ) (16) where for 3D problem, Iis a3×3identity matrix, strain-displacement matrix BI is given by ( 11),GIJis detailed in ( 13),Dstands for incremental constitutive tensor (which is merely a 6×6Voigt notation of incremental constitutive tensor C=2∂S ∂C, withCbeing the right Cauchy–Green deformation tensor C=FTF). The above (16) contains two additive terms. The first, associated with GIJis sometimes referred to as the ”initial stress” or the ”geometric” stiffness matrix. While the other can be referenced as the ”material” stiffness matrix. Both of them contribute to the non- linearity of KT IJwhich makes analytical integration impossible. 2.1Standard(ST)NumericalIntegration The integrand of the internal forces ( 15) and of the tangent stiffness matrix ( 16) is a nonlinear function, which generally speaking, can not be integrated analytically in the element domain Ωe. Thus, an approximate integration approach must be adopted. Most commonly, commercial FE packages utilize the Gauss quadrature or a similar integration rule based on the use of integration points (e.g. [ 4,37]). Transformation of the volume element admits dv=det(F)dV. Hence, integration in the element domain admits dΩe=Jed□where d□is a parameter space. Consider an arbitrary tensor/vector valued function A(ξ;XI;uJ;c0, .., c k),(I, J = 1, .., n ), where c0, .., c kis a set of material properties. Then, an integral of Ain the element domain is given by ¯A=∫ ΩeAdΩe=∫ Ω□A(ξ;XI;uJ;c0, .., c k)Je(ξ;XK)d□=∫ Ω□AJedξdηdζ (I, J, K = 1, .., n ) In particular, for the hexahedral element it becomes∫+1 −1∫+1 −1∫+1 −1AJedξdηdζ . Fol- lowing the standard numerical integration, the approximation of the integral takes the form ¯AST≈nip∑ q=1WqA(ξq;XI;uJ;c0, .., c k)Je(ξq;XK),(I, J, K = 1, .., n )(17) where nipis the number of integration points, Wq,(q= 1, .., n ip)denotes weights as- sociated with each point – ξq= (ξq, ηq, ζg). It follows immediately from the above (17), that the higher the nipthe more additive terms are included in the summation. Therefore - computational complexity grows linearly with nip. Each quadrature has its order (rank), i.e. scheme of order nwill integrate poly- nomial functions of order nexactly. Typically, the higher the quadrature order the higher nip, the higher the computational cost of element formation. As a general 7 consideration, it is important to lower the element formation cost by lowering nip as much as possible. The most computationally inexpensive integration rule, using which the global system converges (performs sufficiently well), is desirable. 2.2"Full"OrderIntegration Following the previous subsection ( 2.1), an important question arises - What order of quadrature should one use to evaluate ( 14),(15),(16)? The higher the scheme’s or- der the higher the nipand therefore the element formation cost. The lower the order andnip, the less accurate and rank-deficient matrices become. Thus, two important integration strategies have been discovered/formulated and widely adopted — ”Full” order and ”Reduced” order integration. i)Full order integration requires quadrature order to match the integrand order provided regular shape andconstant strain el- ement. ii)Reduced order integration happens when quadrature order is lower than full integration. While full integration could be associated with its own problems such as over-stiffness, the reduced integration results in even worse outcomes such as rank deficient element stiffness matrix, which leads to instability, spurious singular mode, zero-energy, or hourglass mode. As a general guideline, the reduced integration is not recommended for a higher order elements (( 44)pp.62). Broadly speaking, full integration is the default choice, a convergence necessity. There is no need to go higher than full order due to computa- tional costs and still present over-stiffness problems. Also, reduced order integration, if used for linear elements, usually employs a special stabilization technique. Herein, we’ll propose a full integration formula at a significantly lower (than reduced) compu- tational cost. Now it is important to recall and detail the practical mathematical meaning of ”reg- ular shape” & ”constant strain” element. i)Regular shape element has several equiv- alent statements such as - constant metric element or constant Jacobian element i.e. Je(ξ,XI) =const =Je(XI),(I= 1, .., n ), where XIis the initial nodal positions (mesh). Therefore, given a regular shape element, one can sample the Jacobian at the centroid and result in an exact value Je|reg. shape =Je|centroid =J0 e. Visually, for an 8-node brick element, this is true as long as an element is parallelepiped (however skewed). For a 20-node hexahedral element parallelepiped shape, all middle nodes must be in the middle of the edges. In fact, for all higher-order elements, the edges are straight lines while middle nodes are in the middle of the edges. For a wedge element triangular bases are equal and parallel to each other etc. In addition, it im- mediately follows that Jegiven by ( 5) is also coordinate independent Je|reg. shape =J0 e, together with its inverse J−1 e|reg. shape =1 J0eAdj(J0 e) =J0−1 e, (see (8)). ii)Constant strain implies that the strain tensor is coordinate independent, i.e. the deformation gradient ( 10) is also constant Fe=Fe|centroid =F0 e. With the help of (10) and (6) it must hold that je=j0 e,F0 e=j0 eJ0−1 e. In addition, for a constant strain element, the stress tensor must remain constant so that Se|Fe=const =Se|centroid =S0 e. In summary, an integration formula that admits full order integration, is exact pro- 8 vided the below conditions and is an approximation for all the other cases Je|regular shape =const =J0 e,Fe|constant strain =const =Fe|centroid =F0 e Je|Je=const =J0 e,J−1 e|Je=const =J0−1 e =1 J0 eAdj( J0 e) ,je|Fe=const =j0 e Se|Fe=const =Se|centroid =S0 e(18) here and throughout the text, the upper index ”0” denotes evaluation at the centroid. 3Semi-Analytical(SA)Integration In this section, we detail the general concept behind our SA approach to element- level full order integration. Then, we discuss the computational cost of the resulting formula. Next, in subsections ( 3.1),(3.2),(3.3) we specialize it to the element mass matrix, the internal forces and the nonlinear tangent stiffness matrix respectively. Consider an integral ¯I=∫ ΩeI dΩe=∫ Ω□I(ξ;XI;uJ;c0, .., c k)Je(ξ;XK)d□,(I, J, K = 1, .., n )(19) An integrand Ican be either a tensor/vector/scalar-valued function, which depends on the local coordinates ξ= (ξ, η, ζ ), initial nodal positions (mesh) XI, nodal displace- ments uJ(t), a set of material properties c0, .., c k. Most commonly, integrand Iis a nonlinear function with respect to local coordinates ξas well as to other arguments. Our goal is to derive an approximate formula for the above integral ( 19), such that it will be exact for ( 18), i.e., admit full integration (see subsec. ( 2.2), while approximate for all the other cases. To this end, we i)we decompose the integrand of the above ( 19) to two multiplicative parts I=I1(ξ;XI;uJ;c0, .., c k)I2(ξ) ¯I=∫ Ω□Je(ξ;XK)I1(ξ;XI;uJ;c0, .., c k)I2(ξ)d□,(I, J, K = 1, .., n )(20) Given regular shape & constant strain element ( 18), the first part JeI1can be exactly evaluated using ONE sampling point (JeI1)|centroid =J0 eI0 1. The second part I2(ξ) consists of mesh-independent and displacement-independent polynomials with re- spect to local coordinates ξ. Importantly, even though the integrand Iis a nonlinear function with respect to local coordinates ξ, provided ( 18), the above decomposition (20) is always possible, as is revealed by careful examination of the mathematical struc- ture of the element-level computations. ii)Thus, SA approximation of ( 20) takes the form ¯ISA≈¯I|Je=const Fe=const=J0 eI0 1∫ Ω□I2(ξ)d□ Thus, the integration rule becomes ¯ISA=J0 eI0 1I∗ 2, I∗ 2=∫ Ω□I2(ξ)d□ (21) 9 Importantly, I∗ 2is merely a number/coefficient/weight. Prior to code implementa- tion of the above formula ( 21),I∗ 2have to be adopted from the literature (if exists) or to be pre-computed. One needs to pre-compute it only once since it is both mesh- independent and deformation-independent. Practically, I2is either a shape func- tion/shape function multiplication (e.g. NI(ξ),NINJ), or its derivative/derivatives multiplication (e.g. NI;2=∂NI ∂ηorNI;2NJ;3). In some sense, I∗ 2resembles weights Wq, (q= 1, .., n ip), from the standard approach ( 17). However, weights Wqare associated with a specific integration point in a particular quadrature; weights are not related to the integral being evaluated. To the contrary, I∗ 2makes sense only for a case-specific integrator ( 21), since I2(ξ)is a specific polynomial function which appears in the in- tegrand decomposition ( 20). Why do we loosely state that roughly speaking, the SA approach is computationally equivalent to the One-Point quadrature? Originally, we are evaluating ( 19), which can be tackled by numerical integration ( 17) in a straightforward manner ¯I=∫ Ω□JeI d□ST≈W1J1 eI1+W2J2 eI2+..+WnipJnip eInip(22) where Jq e=Je|ξq,Iq=I|ξq,(q= 1, .., n ip). Alternatively, one can apply the ST scheme to the decomposed form ( 20) of the same integral ( 19) i.e. ¯I=∫ Ω□JeI1I2d□ST≈W1J1 eI1 1I1 2+W2J2 eI2 1I2 2+...+WnipJnip eInip 1Inip 2 (23) where Iq 1=I1|ξq,Iq 2=I2|ξq,(q= 1, .., n ip). Clearly, both ( 22) and ( 23) produce the same result, however, strictly speaking, computationally they are not exactly the same. Our full integration SA formula ( 21), i.e., ¯ISA=J0 eI0 1I∗ 2is equivalent to one additive term in the above ( 23). In summary, our formula ( 21) is equivalent to One-Point of (23), while each additive term in ( 23) isroughly equivalent to each term in the original ( 22). What order quadrature should one use to evaluate ( 19) using the full integration ST approach? Following ( 20) and ( 18), integral ( 19) is given by ¯I|Je=const Fe=const=∫ Ω□J0 eI0 1I2(ξ)d□ Apparently, the exact evaluation of the above, requires quadrature order to match the order of I2(ξ). For I2of the 2ndorder vs. the 4thorder, the number of integration points nipcan vary quite significantly, especially in 3D. While for the SA approach, the order of I2is irrelevant to the computational complexity of ( 21), since I∗ 2is still one number. Thus, the higher the order of I2the higher the computational efficiency advantage offered by the proposed method. Regarding code implementation. Often, Semi-Analytical work in the field is asso- ciated with meta-programming or extremely lengthy (yet efficient) codes generated by CAS (closed-form approach). None of this is necessary for the proposed method. To illustrate, following the present case, the integration formula ( 21) is similar to the one-point of the ST quadrature use ( 22). The very same is true regarding subroutine implementation. 10 Why the proposed formula ( 21) to evaluate ( 19) is so efficient with respect to ( 23)? For instance, for I2(ξ)of a 6thorder in 3D, one might need dozens of integration points (depending on the domain), while the proposed ( 21) is still equivalent to only one. Quadrature (cubature), such as Gauss points, is well established, the most effec- tive, and the most general integration scheme, reaching far beyond mechanics and NFEM. They are over-performed due to the highly specialized nature of the pro- posed formulas. Our approach exists in the element-level NFEM domain. It relies strongly on the existence of the decomposition ( 20), which in turn has meaning due to the existence of the full integration conditions ( 18), and the coefficients are pre- computed using analytical integration, which is possible once again due to unique nature of the integrand. It is important to emphasize, that SA formulas are case- specific integrators. We combine an essential convergence criteria ( 18) together with a special mathematical structure of the (element-level) integrals to analytically decom- pose the integrand ( 20). Even though, the resulting formula ( 21) somehow resembles the standard approach ( 22),(23), one can not apply this formula to an integral with a different structure. For instance, our formula for the nonlinear tangent stiffness ma- trix (30) can not be used for the mass matrix or the internal forces and vice versa. In contrast, the standard approach quadrature can tackle a general problem. In addition, while the standard approach is versatile in terms of rank i.e., one can apply a variety of distinct quadrature to evaluate a problem ( 19) with different levels of accuracy, our SA formula admits the full integration criteria, no less & no more. Luckily, it fits the demand of a specific problem being tackled in the study. 3.1ConsistentMassMatrix Semi-Analytical integration for the consistent and lumped (diagonal) mass matrices has been thoroughly discussed in [ 15,17,18]. However, herein we omit the higher- order models and include only the basic One-Point full integration model, which keeps this contribution complete and coherent. The decomposition ( 20) of (14) comes quite naturally (NI(ξ)NJ(ξ)) (ρ0Je(ξ;XK)),(I, J, K = 1, .., n ) where the first part ρ0Jeis coordinate and mesh-dependent; to admit the full integra- tion (18), it is sufficient to evaluate it at the centroid ρ0 0J0 e. The second part NINJis a mesh-independent & displacement-independent polynomial function, the analytical integration of which results in a set of numerical coefficients/weights MSA IJ=M0 IJρ0 0J0 e, M0 IJ=∫ Ω□NINJd□,(I, J = 1, .., n ) (24) where a set of coefficients M0 IJis either adopted from the literature or pre-computed once prior to subroutine implementation. Each element type results in a different set. For example, shape functions NIof a 10-node tetrahedral element, are quadratic with respect to the local coordinates (e.g., N7= 4η(1−ξ−η−ζ)). Coefficients M0 IJ are pre-calculated following the rule M0 7 8=∫+1 0∫1−ζ 0∫1−ζ−η 0N7N8dξdηdζ =2 315, see [11]. 11 3.2InternalForces Internal forces are given by ( 15). With the help of ( 11) & (4), strain-displacement matrix BI,(I= 1, .., n )is decomposed such that BI=3∑ r=1ˆBrNI;r,(I= 1, .., n ) ˆBr=1 Je F1,1ar,1 F2,1ar,1 F3,1ar,1 F1,2ar,2 F2,2ar,2 F3,2ar,2 F1,3ar,3 F2,3ar,3 F3,3ar,3 F1,1ar,2+F1,2ar,1F2,1ar,2+F2,2ar,1F3,1ar,2+F3,2ar,1 F1,2ar,3+F1,3ar,2F2,2ar,3+F2,3ar,2F3,2ar,3+F3,3ar,2 F1,1ar,3+F1,3ar,1F2,1ar,3+F2,3ar,1F3,1ar,3+F3,3ar,1 (25) where arn,(r, n= 1,2,3)are entries of the adjoint matrix arn= [Adj(Je)]rn, see (7). Using the above ( 25), the integrand of ( 15) takes the form∑3 r=1JeBT rSeNI;r(ξ), which admits ( 20) given the ( 18). Then, the SA approximation to ( 15) takes the next form RSA I=3∑ r=1J0 eB0T rS0 eN∗ Ir=J0 e( B0T 1N∗ I1+B0T 2N∗ I2+B0T 3N∗ I3) S0 e N∗ Ir=∫ Ω□NI;r(ξ)d□,(I= 1, .., n , r = 1,2,3)(26) where N∗ Iris a set of coefficients/weights having real numerical values. Those coeffi- cients have to be one-time pre-computed prior to code implementation e.g. for a 20- node brick element N16;3=∂N16 ∂ζ=(ξ−1)(η2−1) 4,N∗ 16 3=∫+1 −1∫+1 −1∫+1 −1N16;3dξdηdζ = 4 3. 3.3NonlinearTangentStiffnessMatrix The nonlinear tangent stiffness matrix is given by ( 16). The integrand consists of two additive terms, the first often referred to as the ”geometric” or ”initial stress” stiffness matrix, while the other as the ”material” stiffness matrix. Following the SA approach described in Sec. ( 3), we decompose both additive terms. GIJis given by (12), with the help of ( 4) & (13),GIJ,(I, J = 1, .., n )is rewritten as GIJ=[1 J2 e[Adj(Je)] [Se][ AdjT(Je)]] •[ (∇ξNI) (∇ξNJ)T] (27) where bold dot between two matrices of the same rank ( 3×3) denotes scalar prod- uct (double contraction) i.e. C= [A]•[B],C=∑3 m,n=1[A]mn[B]mn. Hence, SA approximation takes the next form ∫ ΩeGIJdΩeSA≈GSA IJ=1 J0 e[ Adj(J0 e)] [S0 e][ AdjT(J0 e)] •[ K0 IJ] [K0 IJ]mn=∫ Ω□NI;mNJ;nd□,(I, J = 1, .., n , m, n = 1,2,3)(28) 12 Once again we accentuate, coefficient matrices K0 IJ=∫ Ω□(∇ξNI) (∇ξNJ)Td□are ei- ther one time pre-computed prior to in code implementation or adopted from the lit- erature (if exist). Coefficient matrices are element-specific due to the element-specific nature of the shape functions and integration domain. Decomposition of the ”material” part, uses representation ( 25). Accordingly, the integrand becomes BT IDB J=3∑ r,t=1BT rDB t(NI;rNJ;t),(I, J = 1, .., n ) (29) thus, full-order SA approximation takes the form ∫ ΩeBT IDB JdΩeSA≈3∑ r,t=1J0 eB0T rD0B0 t[ K0 IJ] rt Finally, the overall SA approximation of fully integrated ( 16) is given by KSA T IJ =GSA IJI+3∑ r,t=1J0 eB0T rD0B0 t[ K0 IJ] rt,(I, J = 1, .., n ) (30) where GSA IJis given by ( 28). 4Discussionandconclusions Currently, commercial FEA software packages, mostly use the standard numerical integration for element mass matrix, internal forces, and the nonlinear tangent stiff- ness matrix. For example, ANSYS adopts a 14-point quadrature to evaluate the mass and stiffness matrices of the 20-node hexahedral element, while ABAQUS adopts a 27-point quadrature for full and 8-point scheme for reduced integration. Clearly, the more integration points the higher the element formation cost. In this study, we specialize, for the first time, in the recently suggested Semi-Analytical approach to derive a new, full-order, integration formulae for element mass & stiffness matrices plus the internal forces. Our formulas admit ”full” order integration i.e. they are exact for a regular shape & constant strain element ( 18), regardless of an element shape or an order. Code wise, the suggested schemes ( 24),(26),(30), are implemented in a conventional manner. No meta-programming is required; no explicit code generation (closed- form) by means of the Computer Algebra Systems (CAS) is necessary. The present contribution is neither restricted by a specific element or an element or- der nor by a particular constitutive behavior. All solid elements can be implemented, whether 2D or 3D, whether wedge, triangular, bricks, quadrilateral, tetrahedral, etc., whether linearly elastic, hyper-elastic, or other constitutive behavior. In terms of efficiency, Our formulae are roughly equivalent to One-Point quadra- ture, regardless of an element order! Subsection ( 2.1) elaborates on this loose state- ment. Throughout the study, we adopt the nonlinear Total Lagrangian formulation. How- ever, specialization to the Updated Lagrangian formulation follows exactly the same 13 guidelines developed here. Also, external forces are omitted from the study, yet, given the guidelines, their inclusion follows a straightforward manner. Finally, we’ve explicitly considered the isoparametric elements, due to their widest commercial ac- ceptance. Once again, the SA approach can be specialized for other element concepts if necessary. The higher the element order, that is, the higher the shape functions order, the more pronounced the efficiency advantage of the proposed method, see section ( 3). The consistent mass matrix ( 14) is deformation independent. Consequently, it is sufficient to calculate it only once ( 24). Therefore, even though its efficient evalua- tion is desirable, it is hardly the ”bottleneck” of the global nonlinear system solution. In contrast, the nonlinear stiffness matrix is to be evaluated every step of the global system solution, which makes its computational efficiency improvements extremely in demand. SA formulas ( 24),(26),(30) contain coefficient sets. Numerical values of those coef- ficients can be obtained via literature (e.g.[ 17],[18]) or by one-time pre-computation. Clearly, those coefficients are element-specific, due to the element-specific nature of the shape functions and the integration domain. To illustrate, coefficients for a 10- node tetrahedral element differ from coefficients for a 15-node wedge element. SA integration rules ( 24),(26),(30) are case-specific integrators — mass matrix for- mula ( 24) cannot be applied to the stiffness matrix ( 16) or to the internal forces ( 15) and vice versa.