Mixed Integration Scheme for Embedded Discontinuous Interfaces by Extended Finite Element Method

The extended Finite Element Method (XFEM) is derived from the traditional finite element method for discontinuous problems. It can simulate the behavior of cracks, which significantly improves the ability of finite element methods to simulate geotechnical and geological disaster problems. The integration of discontinuous enrichment functions in weak form and the ill-conditioning of the system equations are two major challenges in employing the XFEM in engineering applications. A mixed integration scheme is proposed in this paper to solve these problems. This integration scheme has a simple form and exhibits both the accuracy of the subcell integration method and the well-conditioning of a smeared integration method. The correctness and effectiveness of the proposed scheme were verified through a series of element analyses and two typical examples. For XFEM numerical simulations with unstructured meshes and arbitrary cracks/interfaces, this method guarantees the convergence of nonlinear iterations and yields correct results.


INTRODUCTION
The extended Finite Element Method (XFEM) is an excellent numerical method that can visually simulate discontinuous cracks/interfaces and their evolution without mesh regeneration. By introducing enrichment functions with special local properties into the partition of unit method (PUM) framework, this method can enable construction of discontinuous fields within local elements (Belytschko and Black, 1999;Moës et al., 1999;Belytschko et al., 2001). As it offers several advantages, XFEM has been applied to failure problems in many different fields, including geotechnical and geological engineering (Abdelaziz and Hamounie, 2008;Belytschko et al., 2009;Fries and Belytschko, 2010;Gracie and Graig, 2010;Salimzadeh and Khalili, 2015;Matthew and Caglar, 2016;Li et al., 2018;Wang et al., 2018;Cruz et al., 2019). It can be predicted that XFEM will play a more important role in the modeling, assessing, and preventing of geotechnical and geological disasters. However, the wider application of XFEM is limited by several drawbacks. One of these drawbacks is that with the introduction of discontinuous functions, standard Gaussian integration cannot be applied directly to XFEM-enriched elements, and splitting of elements into subcells is often necessary. XFEM converts the complexity of mesh regeneration into the complexity of the discontinuous function integration within the element (Chin et al., 2017). Another drawback is that due to the same interpolation basis functions used, the enriched degrees of freedom (DOFs) of a node may become almost linearly dependent on the regular DOFs of the node, leading to ill-conditioned system equations. The nearly linear dependency and ill-conditioning are especially apparent for discontinuous enrichment functions when the interface approaches an element node and crack tip singular enrichment functions at a certain distance from the singular point (Sukumar et al., 2015;Ventura and Tessei, 2016;Agathos et al., 2019).
Generally, integrating the weak form within an element with discontinuous interfaces requires element splitting. This approach has strong adaptability and has been widely used in previous studies (Moës et al., 1999;Belytschko et al., 2001;Belytschko et al., 2009;Fries and Belytschko, 2010;Li et al., 2018). However, implementation of the element-splitting algorithm is complicated and not uniform. Additionally, it presents challenges concerning the integration point mapping in nonlinear constitutive models. Studies have been devoted to the integration of non-polynomial enrichment functions (discontinuous functions or singular functions) without element splitting. The representative methods include the equivalent polynomial integration method (Ventura 2006;Ventura and Benvenuti, 2015), the arbitrary polygon integration method (Natarajan et al., 2009;Natarajan et al., 2010), the arbitrary domain integration method (Joulaian et al., 2016), and the Legendre polynomial equivalent integration method (Abedian and Düster, 2019), et al. At present, these methods are usually limited to two-dimensional (2-D) conditions; three-dimensional (3-D) implementations remain a challenge. In addition, these methods cannot address curved or kinked cracks/interfaces within an element. Song et al. (2006) used a partition ratio-weighted single-point reduced integration to integrate the weak form of quadrilateral elements with discontinuous enrichment functions, thereby avoiding element splitting and historical variable mapping. However, single-point reduced integration is equivalent to assuming that the element is constantly strained, which results in poor solution accuracy (Wang et al., 2021).
For unstructured meshes and arbitrary non-planar cracks/ interfaces, cases with cracks/interfaces close to element nodes are inevitable, which lead to ill-conditioned system equations. Previous approaches for solving this problem include directly removing the enrichment of nodes whose enrichment functions have only small supports in the cut element (Daux et al., 2000;Bordas et al., 2007), and adjusting the node coordinates such that the partition ratio is not unreasonably small (Choi et al., 2012). However, the former approach introduces interpolation errors, while the latter is not a generalizable method, especially for 3-D problems. Reusken (2008) proposed a method to impose constraints on the enriched DOFs with a small partition ratio, limiting its value to zero; this is essentially the same as removing the enrichment. Ventura and Tesei (2016) proposed a stabilized method that imposes constraints on the enriched DOFs through a penalty method, reducing the additional error caused by the constraints. Other studies have adopted the approach of preconditioning the ill-conditioned system equations to improve convergence. For example, Sauerland and Fries (2013) studied the Jacobi preconditioner, and Béchet et al. (2005) developed a preconditioner based on Cholesky decomposition. However, this approach becomes cumbersome for nonlinear problems and crack evolution problems because the Jacobian matrix must be re-preconditioned after each update, which significantly increases the amount of calculation.
Based on the research of Song et al. (2006), the partition ratioweighted integration method was improved in this study and is referred to as "smeared" integration method. It was observed that this method could significantly alleviate the ill-conditioning caused by the nearly linear dependence of DOFs enriched by the discontinuous functions. Furthermore, this method (for the integration of the Jacobian matrix) was combined with the subcell (element splitting) integration method (for the integration of residuals) to produce a mixed integration scheme, which not only resulted in well-conditioned system equations but also ensured a high solution accuracy.
The remainder of this paper is structured as follows: Firstly, the XFEM formulation is sketched out, and the weak form of the equilibrium equation is presented, including the interface contact relationship and the integral form of the Jacobian matrix. Secondly, the smeared integration method is presented, and its performance is demonstrated through a series of element analyses. Then, the mixed integration scheme is proposed and verification is made in the next section through two benchmark examples. Finally, the main conclusions of this study are summarized in the last section of this paper.

Strong and Weak Form of Boundary Value Problem
A quasi-static boundary value problem with internal interfaces was assumed, as shown in Figure 1. In the solution domain Ω, there exists a strong discontinuous interface Γ d , the two sides of which are denoted as Γ + d and Γ − d ; the corresponding normal vectors are n + and n − , respectively. b is the body force in the domain; Γ t is the force boundary, and Γ u is the displacement boundary.
The force equilibrium conditions must be satisfied in the continuum, and the corresponding Dirichlet or Neumann boundary conditions must be satisfied at the boundaries. Thus, the governing equations of the boundary value problem are Contact forces t + and t − are defined on the internal interfaces Γ + d and Γ − d , respectively. They satisfy the equilibrium condition and can be expressed as a unified variable: correspond to the displacement vectors on the two sides of the internal interface; the relative displacement on the interface is defined as The kinematics admissible (trial) space and the weighting (test) space of the displacement are defined as For u(x) ∈ U and any δu(x) ∈ U 0 , the equivalent integral weak form of the governing equations can be expressed as where δW int , δW ext , and δW c are the work of the internal force, external force, and interface contact force, respectively, and are expressed as where  s is the symmetric gradient operator. According to the definitions of t (Eq. 2) and w (Eq. 3), the expression of δW c in Eq. 7 can be simplified as The constitutive relationship of the continuum is expressed as where D is the constitutive tensor. The traction-displacement relationship of the interface can also be expressed in a constitutive form as where D Γ is the interface constitutive tensor, implemented by imposing constraints on the interface in practice; Q is the transfer matrix from local coordinates to global coordinates.

Discretization and Jacobian Matrix
For simplicity, only a case with elements completely separated by the interface was considered, and only the Heaviside enrichment function was used in this study. Nevertheless, the conclusions of this study also hold for a case that the end of the interface stays inside the element, for which treatments have been provided by Zi and Belytschko (2003), Kumar and Bhardwaj (2018). The displacement approximation in a shifted form is expressed as where the first summation term represents the standard FE approximation, and the second summation term represents the XFEM-enriched discontinuous approximation; I is the set of all nodes in the discrete domain; I*∈I is the set of enriched nodes; N i is the regular finite element basis function of node i; u i and q i are the regular DOF vector and the enriched DOF vector of node i, respectively; and H is the Heaviside enrichment function, which is expressed as where f(x) is a signed distance function depending on the definition of the positive normal of the discontinuous interface, as shown in Figure 1.
Considering the variation in the discrete approximation of the displacement field given by Eq. 11, the testing field can be obtained: According to the definition of the interface relative displacement w given by Eq. 3, the continuous term of the displacement field has no contribution to w; instead, w is determined by the discontinuous term, as follows: Its variational form is expressed as By substituting the testing fields expressed in Eqs 13, 15 into the equivalent integral weak form of Eq. 6, and considering the arbitrary nature of these fields, the following discrete form of the equilibrium equations can be obtained: where f int u and f ext u are the nodal internal force vector and nodal external force vector, respectively, corresponding to the regular DOFs; f int q and f ext q are the nodal internal force vector and nodal external force vector, respectively, corresponding to the enriched DOFs; and f c is the nodal interface contact force vector. All these vectors are integrated from the nodal force vectors of related elements, and their specific expressions are presented in Eq. 17: where B is the discrete strain-displacement operator, and H and H I are the matrix forms of the Heaviside function and its nodal values, respectively. Taking the derivative of Eq. 16 with respect to the displacement fields and considering Eq. 17, the Jacobian (stiffness) matrix of the finite element discretization equations can be obtained as When the same basis function is used, the Jacobian matrix expressed as Eq. 18 may have a large condition number, especially when the interface is close to the enriched nodes. Thus, the enriched DOFs and the standard DOFs are almost linearly dependent, and the matrix tends to be severely illconditioned.

Smeared Integration Method
Based on Eqs 12, 19 and considering H 2 H, the numerical integration of the discontinuous part (multiplied with the Heaviside function) for each term in the stiffness matrix can be generalized in the following form: where p(x) is an arbitrary polynomial function. Evidently, the discontinuity in the integrand arises from the discontinuity in H(x). As shown in Figure 2, a directional interface cut the element into two subdomains Ω A and Ω B , denoting the positive and negative parts of the signed distance function, respectively. The integration of the discontinuous function over the entire element is converted to a summation of continuous function integrations over the two subdomains [H(x) is constant on each side]: (21) Noted that H(x) equals to unit on Ω A and vanishes on Ω B , Eq. 20 can be simplified as As the subdomains can be arbitrary polygons, the calculation of Eq. 22 often requires further division of the subdomains into triangular subcells, as illustrated in Figure 3, which brings much burden for calculation.
To avoid element decomposition, Song et al. (2006) integrated Eq. 20 over the entire element domain via single-point reduced integration with hourglass control. Their starting point was the assumption that all elements are constant-strain elements [p(x) C]; thus, the integration over the subdomain could be replaced with integration over the entire element domain with a coefficient representing the subdomain contribution: where R is the partition ratio of the element cut by the interface and is defined as Clearly, R is related to the location and configuration of the interface; thus, it can reflect the distribution characteristics of the discontinuous function defined on the element to some extent. The rightmost integral in Eq. 23 is performed over the entire element domain; thus, standard Gaussian integration can be performed directly. However, the assumption of constantstrain elements reduces the integration accuracy and prevents an accurate reflection of the constitutive behaviors of nonlinear materials.
Based on the method of Song et al. (2006), single-point reduced integration was replaced with quadratic Gaussian integration in this study, as shown in Figure 2B. This alteration obviated the constant-strain assumption and the Frontiers in Earth Science | www.frontiersin.org January 2022 | Volume 9 | Article 829203 need for hourglass control. The integral domain remained unchanged, and the influence of H(x) was reflected approximately by R: This integration method is simple and can be easily implemented. The integration points remain consistent before and after the discontinuity interface is introduced; thus, no historical variable mapping is required. However, this method ignores the shape information of the real integration domain and introduces information from the integrand function outside the domain. This approach is similar to a smearing/homogenizing operation; it can be referred to as "smeared" integration and is expressed as Although this method suffers from certain approximation errors, it offers several advantages: (i) Compared with subcell integration, the smeared integration method avoids element splitting and significantly reduces the number of required integration points. (ii) In each element, the Gaussian integration of B T DB needs to be performed only once. Only K uu in Eq. 19 must be integrated; while the integration of K uq and K qq can be achieved by multiplying K uu with corresponding coefficients related to the partition ratio R, which further reduces the calculation complexity. The detailed relations between these matrixes are derived in the Supplementary Appendix. (iii) The stiffness matrix conditions are improved, and the nonlinear iteration stability is enhanced. These features are demonstrated in detail in the following sections.

Error Analysis of Smeared Integration
If the continuous integrand p(x) is a constant, the smeared integration is accurate, as expressed in Eq. 23. If p(x) is a linear or higher-order polynomial, the smeared integration has a certain approximation error that is, directly related to the order of p(x) and the location of the discontinuous interface relative to the element (represented by the partition ratio R).
Defining the approximation error of the smeared integration (Eq. 25) by comparing with the exact discontinuous function integration (Eq. 21): (27) For low-order continuous integrand p(x), there exists a point x i in the domain Ω A , satisfying:  Similarly, there exists a point x j in the domain Ω B , satisfying: Where S A and S B are the areas of Ω A and Ω B , respectively. Noted that, S A RpS, S B (1-R)*S, where S is the area of the intact element Ω, Eq. 27 can be rewritten as: Several conclusions can be drawn from this expression: (i) For a fixed partition ratio R, a higher order of p(x) produces a greater difference between p(x i ) and p(x j ), thus greater approximation error in the smeared integration. Order of p(x) decides the upper bound of the error. For the bilinear/ trilinear elements used in this study, the errors were limited and were influenced more strongly by the location of the interface. (ii) R behaves as a scale factor of the error. When R → 0 or R → 1, the error approaches zero with the same order with R or 1-R. Actually, these two conditions are equivalent as they are interconvertible by changing the definition of the interface direction in Figure 2. (iii) When R → 0.5, the smeared integration has a relatively large error because the shape information of the real integration domain is ignored. For the three cases illustrated in Figure 4, the same stiffness matrix is calculated from the smeared integration; thus, smeared integration is not recommended in this condition.

Influence of Smeared Integration on Element Stiffness Condition
To understand the effects of partition ratio R and the configuration of the discontinuous interface on the efficacy of smeared integration, element tests for four typical cases, as shown in Figure 5 1 , were performed. The element stiffness matrices were calculated via subcell integration and smeared integration. The square element has a side length of a 1 and an elastic modulus of E 1.0. Only four integration points are required for each cracked element by the smeared integration while the subcell integration involved the third-order Hammer integration with four integration points in each triangle subcell.
Interface contact was fixed in a bonded state, and the penalty method was employed to impose a spring-like constraint on the interface. The penalty parameter k was assigned a value ten times of the elastic modulus. The integration of tractions on the interface was unified as two-point Gaussian integration. On each element node, there were two regular DOFs and two enriched DOFs. The element stiffness matrix was formed into a 16 × 16 square matrix. Three rigid-body displacement DOFs need to be removed when calculating the condition number.
The contents of Table 1 2 and Figure 6 correspond to the partition patterns shown in Figure 5 and present the condition numbers of the element stiffness matrix obtained via smeared integration and subcell integration, respectively. With the same partition pattern, the condition numbers of the stiffness matrix obtained through smeared integration were smaller than those obtained through subcell integration. For partition pattern (A), smeared integration produced the most significant reduction in the condition number; when R was relatively small, the difference reached several orders of magnitude. In double logarithmic coordinates, the condition number of the stiffness matrix exhibited a nearly linear relationship with R. The smeared integration exhibited a slope of approximately −1.1, compared with −1.9 for the subcell integration. This result is consistent with the analysis in the last section; when R is relatively small, smeared integration replaces the small quantities with larger ones, thus alleviates the ill-conditioning of the element stiffness matrix. For partition patterns (B), (C), and (D), the effect of smeared integration was successively weakened. The stiffness matrix of this bilinearly interpolated square element was integrated in the space of span (x, y, x 2 , y 2 , xy); partition patterns (B), (C), and (D) reduced the order of the integrand function. In partition pattern (D), the variable x was fixed and exhibited no dependency on R; partition patterns (B) and (C) represent the intermediate transition states.

Analysis of Coupling With Interface Contact
For frictional contact problems with interface constraints imposed using the penalty method, the accuracy and stability depend on the penalty parameter k. In this study, a penalty coefficient α was used to define the penalty parameter: For the case shown in Figure 5A, some interesting phenomena are observed when studying the condition number with different penalty parameters. Here, three rigid-body DOFs must be removed from the enriched DOFs to ensure that the analyzed stiffness matrix is non-singular in the absence of interface constraints.  FIGURE 6 | Relationship curves of condition number with regard to partition ratio corresponding to four partition patterns (see Figure 5) obtained using two integration methods.
Frontiers in Earth Science | www.frontiersin.org January 2022 | Volume 9 | Article 829203 The influence of interface constraints imposed by the penalty method on element stiffness matrix includes two aspects: 1) The interface constraints restrict relative movements between the two sides of the interface, hence improving the stiffness matrix's condition.
2) The condition number of the stiffness matrix may increase significantly with a large penalty parameter 2 | Condition numbers obtained using two integration methods with different penalty coefficients and partition ratios.   Frontiers in Earth Science | www.frontiersin.org January 2022 | Volume 9 | Article 829203 8 involved (usually much larger than other parameters in the stiffness matrix).

Subcell integration Smeared integration
Table 2; Figure 7 present comparisons of the stiffness matrix condition numbers obtained using smeared and subcell integration with different penalty coefficients and partition ratios. The penalty coefficient α 0 indicates that no penalty constraint has been imposed. It is observed that as R approached 0, the condition number increased in all cases in a nearly linear manner; however, the condition number corresponding to smeared integration was much smaller than that corresponding to subcell integration, and the increasing slope was much smaller. For subcell integration, the condition number showed almost no change after an interface constraint with a penalty coefficient α 10 was imposed; the condition number increased significantly with further increase in α, indicating that aspect 2) in the pre-analysis played a dominant role. However, for smeared integration, after an interface constraint with a penalty coefficient α 10 was imposed, the condition number exhibited a slight decrease; only marginal increase was observed with further increase in α. This indicates that aspect 1) in the pre-analysis played a dominant role. With the coupling of interface contact and internal stiffness, the smeared integration method exhibits a coupling effect that can maintain the stability of the overall condition number, thereby improving the solution stability and the convergence of nonlinear iterations.

MIXED INTEGRATION SCHEME FOR XFEM
When using Newton's method to solve nonlinear problems, the Jacobian matrix must be updated in each iteration step; the condition number of the Jacobian matrix can significantly affect the stability and convergence of the nonlinear iterations. The iterations can hardly converge with an ill-conditioned Jacobian matrix (Fries and Belytschko, 2010;Belytschko et al., 2014). In practice, an approximate but well-conditioned Jacobian   matrix could be better than an accurate but ill-conditioned Jacobian matrix for iteration stability. An approximate Jacobian matrix is commonly employed for problems with an asymmetric Jacobian matrix, such as non-associated elastoplastic constitutive and frictional contact (Chen and Cheng, 2011). Based on the element analysis of smeared integration presented in the previous section, a mixed integration scheme for XFEM was proposed. A critical partition ratio R* was set, and the stiffness of the element cut by the interface was integrated with respect to the element partition ratio R, as shown in Figure 8. When R was less than R* or greater than 1-R* (the situations in which illconditioning occurs), smeared integration was used to improve the stiffness condition; when R was in the middle range, subcell integration was used to obtain a higher integration accuracy. For all residual integrations, subcell integration was used to obtain accurate residuals and ensure accurate results. Based on previous element analyses and the example analyses in next section, the value of R* should be set in the range of 0.05-0.2 to achieve a suitable balance between convergence rate and convergence stability. The partition ratio R, along with the areas or volumes of the subcells and other crack information for a cracked element, are computed in the crack processing procedure once the element reaches the cracking criteria, and are stored for reutilization thereafter. For fixed interface problems mainly discussed in this paper, this information is computed and stored as initial information at the very beginning of the solution.

EXAMPLE VERIFICATION
Two typical benchmark examples were used to verify the robustness and practicability of the mixed integration scheme.   interface is set at θ tan −1 (0.2). The square has an elastic modulus of E 1,000 MPa, a Poisson's ratio of ] 0.3, and an interface friction coefficient of μ 0.19. The penalty coefficient is set as α 10,000. This example considers a regular mesh and an irregular mesh, as shown in Figures 9B,C, respectively. The regular mesh consists of 40 × 40 square elements. To prevent the interface from passing exactly through the element nodes, the interface is shifted upward by 0.01 m. The irregular mesh consists of 1857 arbitrary quadrilateral elements.
The partition ratio distributions of all cut elements were plotted for both sets of meshes, arranged from small to large, as shown in Figure 10. In both meshes, there were several cut elements with an R value tending to extreme conditions (R < 0.1 and R > 0.9), indicating poor conditions of the Jacobian matrix.
The condition numbers of the overall Jacobian matrix for different penalty parameters when using the mixed integration scheme and the total subcell integration scheme are presented in Table 3. For the mixed integration scheme, the critical partition ratio R* is set as 0.1. The mixed integration scheme significantly reduced the condition number in both meshes; the reductions in the regular and irregular meshes were greater than one order and two orders of magnitude, respectively.
A negligible difference was present between the calculation results of the mixed integration scheme and the total subcell integration scheme. A comparison of the calculation results for the regular and irregular meshes also showed a negligible difference. Thus, Figure 11 shows only the displacement contours of the irregular mesh calculated using the mixed integration scheme.
The convergence profiles of the two sets of meshes when using the mixed and total subcell integration schemes are shown in Figure 12. For the regular mesh, the total subcell scheme required 14 iteration steps, whereas the mixed scheme required 36 iteration steps to converge to the minimum error limit. This indicates that owing to its approximate Jacobian matrix, the mixed scheme could not reach the optimal Newton convergence rate but can still guarantee convergence. For the irregular mesh, the total subcell scheme could not converge within the Newton iterations, due to the ill-conditioned Jacobian matrix. However, the mixed scheme converged to the minimum error limit after 148 iterations; thus, this scheme can reduce the Jacobian matrix condition number and improve convergence. Martin et al. (2015) analyzed a layered 3-D beam under uniaxial tensile loading using XFEM. The beam dimensions and boundary conditions are shown in Figure 13. An interface parallel to the xoy plane divided the beam into upper and lower layers; the tensile stresses applied to the right end of the upper and lower layers were T up 300 MPa and T low 100 MPa, respectively. The material was an elastoplastic body with linear strain hardening behavior, an elastic modulus of E 0 195 GPa, a Poisson's ratio of ] 0, a plastic yield strength of σ s 180 MPa, and a strainhardening modulus (the tangent modulus of the plastic zone) of E p 1.95 GPa. This setting resulted in obvious deformation and stress discontinuity in the element cut by the interface. The material above the interface entered the plastic deformation stage, while that below the interface remained elastic. No friction acted on the interface (i.e., the friction coefficient was μ 0). The normal constraints (non-penetration condition) were imposed using the penalty method, and the penalty coefficient was set as α 10,000. Firstly, a structured hexahedral mesh consisting of 20 × 5 × 3 300 elements was used. The splitting of the hexahedral element into tetrahedrons for subcell integration was performed based on the method used by Martin et al. (2015).

3-D Beam Under Uniaxial Tensile Loading
The partition ratio of the elements cut by the interface are related with the interface location. For the structured hexahedral mesh, the interface was parallel to one face of the element, and the partition patterns for all elements cut by the interface were the same. The cases with R 0.5, 0.1, and 0.01, i.e., the interface was located in the middle, upper 90%, and upper 99% of the element, respectively, were analyzed. The critical partition ratio R* of the mixed integration scheme was set as 0.5. Therefore, the element stiffness of the cut elements in all cases was calculated using the smeared integration method. Figure 14 shows the calculation results for displacement and stress using the mixed integration scheme, which were consistent with the theoretical solution; the residual errors were less than E-10. Table 4 presents the condition numbers of the Jacobian matrix with different partition ratios, as obtained using the two integration schemes, when the upper structure entered the plastic stage. At this point, the stiffness difference between the structures above and below the interface was relatively large, leading to relatively large Jacobian matrix condition numbers. The condition numbers increased with a decrease in partition ratio R. For all cases, the mixed integration scheme produced smaller condition numbers than the total subcell scheme.
The convergence profiles of the two integration schemes with different partition ratios are shown in Figure 15. Owing to accurate integration of the Jacobian matrix, for all three partition ratios, only three iterations were required to converge to a residual error less than E-10 when using the total subcell scheme. However, when using the mixed scheme, 17, 15, and 11 iterations were required to converge to the same error limit for R 0.5, 0.1, and 0.01, respectively. This result indicates that although the mixed scheme could reduce the Jacobian matrix condition number, it could not achieve the optimal convergence rates, due to approximation of the Jacobian matrix; however, it could guarantee convergence.
To determine the adaptability of the mixed integration scheme to a more generalized mesh, an unstructured tetrahedral mesh was analyzed. However, it was difficult to avoid the case with the interface very close to some element nodes, which resulted in larger Jacobian matrix condition numbers.  Frontiers in Earth Science | www.frontiersin.org January 2022 | Volume 9 | Article 829203 The mesh consisted of 3608 quadratic tetrahedral elements, which yielded a non-constant Jacobian matrix. The vertex nodes were enriched, whereas the edge nodes were not. The interface was described using linearly interpolated level sets defined on the vertex nodes.
A total of 405 elements were cut by the interface; the partition ratio distribution of these elements, arranged from small to large, is shown in Figure 16. For different critical partition ratios R*, the Jacobian matrix condition numbers under the two integration schemes are presented in Table 5; Figure 17. A smaller R* indicates a smaller portion of cut elements integrated by the smeared integration method in the mixed scheme. In other words, the smeared integration method plays a less significant role in the mixed scheme; thus, the mixed scheme is more similar to the total subcell scheme. With an increase in R*, the Jacobian condition number decreases sharply at first, as shown in Table 5; Figure 17. However, after R* reaches approximately 0.1, the decrease becomes gentle; this indicates that almost all the ill-conditioned elements are included in the smeared integration portion. Thus, the mixed integration scheme can be considered to have reached its limit in terms of the ability to improve the condition of the matrix. Thus, R* 0.1 was used in the final calculation.
The calculation results for stress and displacement using the mixed integration scheme are shown in Figure 18; these results are consistent with the theoretical solution. The mixed scheme converged to the given minimum error limit after 69 iterations, as shown in Figure 19; however, the subcell scheme did not converge within the Newton iterations due to the illconditioning. These results prove that the mixed integration scheme exhibits good convergence stability for strongly nonlinear problems with unstructured meshes.

CONCLUSION
This study introduced a smeared integration method that can integrate the weak form of elements with discontinuous enrichment without element splitting. The method also avoids variable mapping in the nonlinear constitutive equations. Compared with conventional methods, this method is simple in form and easy to implement, and achieves excellent performance in terms of improving the condition of system equations. The performance of the method was verified through a series of element analyses. A mixed integration scheme was developed by leveraging the advantages of the smeared integration method and the subcell integration method. To improve the condition and convergence, the smeared integration method was used for the stiffness integration of elements with discontinuous enrichment that may be illconditioned; moreover, the subcell method was used for residual integration to ensure the accuracy of the final solution. In the example calculations, the mixed integration scheme significantly improved the system Jacobian condition and provided a good balance between convergence rates and convergence stability, especially with an unstructured mesh.
In general, the proposed integration scheme ensures the convergence of Newton iterations with Jacobian matrix embedded with discontinuous functions, hence promoting the practical application of the extended finite element method in geotechnical and geological disaster assessment and prevention analysis.

DATA AVAILABILITY STATEMENT
The raw data supporting the conclusion of this article will be made available by the authors, without undue reservation.

AUTHOR CONTRIBUTIONS
PY developed the method and wrote the manuscript. QH performed the engineering application. XW designed the research and participated in writing the manuscript. YY revised and edited the final version. ZZ helped perform the analysis with constructive discussions.

SUPPLEMENTARY MATERIAL
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/feart.2021.829203/ full#supplementary-material FIGURE 19 | Convergence profile of the unstructured mesh for two integration schemes.