An improved stress recovery technique for the unfitted finite element analysis of discontinuous gradient fields

Stress analysis is an all-pervasive practice in engineering design. With displacement-based finite element analysis, directly-calculated stress fields are obtained in a post-processing step by computing the gradient of the displacement field—therefore less accurate. In enriched finite element analysis (EFEA), which provides unprecedented versatility by decoupling the finite element mesh from material interfaces, cracks, and structural boundaries, stress recovery is further aggravated when such discontinuities get arbitrarily close to nodes of the mesh; the presence of small area integration elements often yields overestimated stresses, which could have a detrimental impact on nonlinear analyses ( e.g ., damage or plasticity) since stress concentrations are just a nonphysical numerical artifact. In this article, we propose a stress recovery procedure for enhancing the stress field in problems where the field gradient is discontinuous. The formulation is based on a stress improvement procedure (SIP) initially proposed for low-order standard finite elements. Although generally applicable to all EFEA, we investigate the technique with the Interface-enriched Generalized Finite Element Method and compare the procedure to other post-processing smoothing techniques. We demonstrate that SIP for EFEA provides an enhanced stress field that is more accurate than directly-calculated stresses—even when compared with standard FEM with fitted meshes.

Enriched finite element analysis (EFEA) [4][5][6][7][8][9] provides an elegant solution to the use of fitted meshes by decoupling discontinuities from the FE mesh.In these methods, which have become fairly popular in the computational mechanics community because of their versatility and the possibility of using simple (usually structured) FE meshes, the standard FE approximation is extended or enhanced by means of an enriched FE space that properly resolves the discontinuities' kinematics-for example, via a Heaviside function for cracks (strong discontinuities) 5 or distance-based functions for material interfaces (weak discontinuities). 10However, even though complete discontinuity-mesh decoupling is attained, the accuracy of gradient fields greatly depends on the resulting enriched FE space, which may be severely affected by the way elements are cut by discontinuities.
It is well known that finite elements with bad aspect ratios are prone to significant errors in gradient fields. 11But even when meshes are properly constructed, stress fields are usually obtained by taking the gradient of the primal field-therefore less accurate-and are C −1 -continuous at element edges.Consequently, recovery/smoothing techniques that post-process the primal field have been proposed to improve the accuracy of the stress field.Improving gradient fields was first conducted by Brauchli, 12 who obtained a consistent stress field based on the theory of conjugate approximations; 13 because the approximation functions used to interpolate stresses are expressed as a linear combination of linearly independent functions defined over the entire domain, the resulting stress field is continuous across element boundaries.Hinton and Campbell 14 later proposed local and global smoothing techniques based on a least squares procedure, whereby enhanced stresses are obtained by minimizing their difference with respect to directly-calculated stresses.The superconvergent patch recovery (SPR) of Zienkiewicz and Zhu 15,16 was introduced afterwards as a more computationally efficient method that uses an element patch surrounding the nodes where stress is sought; in this technique, the same shape functions used to interpolate the displacement field are used to construct a smooth stress field, whose polynomial form is obtained by solving a least squares problem at superconvergent points (locations where stress accuracy is the highest within the element).Subsequently, Wiberg and Abdulwahab 17 proposed an enhanced version of SPR by taking the equilibrium equation into account (SPRE), where the residual of the governing equilibrium equation obtained with the smooth stress interpolation is considered in the least squares procedure.However, the quality of recovered stresses at boundaries is worse than that obtained in the interior.Wiberg et al. 18 further proposed an improved extension of SPR considering both the equilibrium equation and boundary conditions (SPREB), where the recovered stress field satisfies the appropriate prescribed boundary conditions by adding the corresponding residuals to the least squares fit.A similar idea was also introduced by Blacker and Belytschko, 19 who only focused on natural boundary conditions.Later, Ródenas et al. 20 proposed an improvement of SPR that uses constraint equations (SPR-C); this method shares similarities with SPREB, where stress interpolation polynomials are required to satisfy, in the local patch, equilibrium and compatibility equations, and boundary conditions.With several modifications, this technique was extended to recover the stress field in singular elasticity problems by decomposing the primal field in smooth and singular components; 21,22 while recovering the smooth part follows the standard SPR-C procedure, an expression based on stress intensity factors (SIFs) that describes the asymptotic fields near the crack tip is used to recover the singular part.
Even though SPR-based techniques perform well in recovering gradient fields, they rely on the existence of superconvergent points; these points are not always defined, for instance in quadratic triangular elements. 23Therefore, several post-processing recovery methods that do not rely on superconvergent points were proposed.Tabbara et al. 24 used a moving least square (MLS) method to construct a local interpolation of the displacement field, where the recovered strain field is obtained by taking the gradient of the displacement polynomial interpolants.Later, the recovery by equilibrium in patches (REP) 25,26 was proposed to extract a stress field by satisfying equilibrium equations in a weak sense over patches of elements; a smooth stress field is obtained by using a least squares scheme, where the error between enhanced and directly-calculated stresses is projected onto the finite element strain space over the patch.Instead of minimizing the difference between recovered and directly-calculated stresses over the patch, Ubertini proposed another recovery technique by considering compatibility in patches (RCP), 27 whereby the strain error (difference between the directly-calculated and enhanced strains) is orthogonal to the space of enhanced stresses over each patch.Rather than using the support of a shape function (i.e., union of elements sharing it) as the patch, 27 Benedetti et al. 28 further developed RCP by considering patches composed of a main target element and all its surrounding elements; once the enhanced stress interpolant is obtained over the patch, the stress in the target element is recovered, and without considering results from adjacent patches.Later, Payen and Bathe 29 proposed a stress improvement procedure (SIP) derived from a mixed formulation based on the Hu-Washizu principle, which is effective in obtaining accurate stress predictions.A mixed formulation is also the only option for getting optimal finite element discretizations when solving problems including (almost) incompressible media, thin structures, or multiphysics phenomena (see Reference 29 and references therein).Within SIP, the point-wise relationship between stress and strain is relaxed; two equations are used to obtain the recovered stress field, one that requires the stress error (difference between enhanced and directly-calculated stresses) is orthogonal to the space of self-equilibrated stresses, and another one that enforces equilibrium over the element patch (in weak sense) with the enhanced stress field.This recovery procedure has several advantages over its predecessors: Unlike RCP, this technique does not use a particular stress field defined a priori that satisfies equilibrium within the patch.Moreover, the number of basic equations used is independent of the number and type of elements considered in the stress calculation domain.Finally, as a limitation of using low-order elements to model nearly incompressible materials, a spurious checkerboard pattern of stresses can be avoided by using this recovery technique. 30Recently, Sharma et al. 31 extended this method to recover the stress field in 3D low-order finite element meshes.The preceding procedures were not proposed for problems with discontinuities that exhibit singular stress fields.Xiao and Karihaloo 32 introduced another recovery approach for fracture problems, which they named the statically admissible stress (SAR) recovery scheme; this method, which is based on complete polynomial functions that satisfy the equilibrium equation and traction boundary conditions in each patch, uses MLS to enforce stress continuity between patches; unlike SPR-C, the smooth and singular stress fields are handled together, where the original stress field near the crack tip is obtained by using hybrid crack elements. 33egarding EFEA, many stress recovery procedures have been proposed in the context of the extended/generalized finite element method (X/GFEM).The SAR technique, which was designed to recover the stress field in fracture problems with standard FEM, was later extended to improve the accuracy of the stress approximation in X/GFEM. 34Duflot and Bordas 35 used the global derivative recovery method-which requires the entire model as the calculation domain-to smooth the strain field for linear elastic fracture mechanics problems; the procedure works by minimizing the square of the L 2 -norm of the difference between the directly-calculated and smoothened strain fields, the latter consisting of three parts (i.e., a smooth component interpolated by standard shape functions, a nonsmooth discontinuous part enriched by a Heaviside function, and an enrichment used to capture the singularity near the crack tip).Jin et al. 36 used this technique to quantify the interpolation error, a measure that is then used to drive mesh adaptivity.Later, the technique was used to recover the strain field for problems that contain only weak discontinuities, with an a posteriori error estimate based on the recovered strain that is used for adaptive local mesh refinement. 37,38As an alternative to the global derivative recovery method, a derivative recovery procedure based on extended moving least squares (XMLS) was proposed to obtain a smooth strain field; 39,40 a smoothened displacement field is first constructed with the original nodal displacements and MLS shape functions, where near-tip asymptotic functions are added to capture singularities; the recovered strains are then obtained by applying the gradients to this new displacement field.SPR was also adapted to X/GFEM (SPR-XFEM) 41 with three major modifications: direct calculation of recovered stresses at integration points (not mesh nodes), decomposing the stress field into singular and smooth parts, and using different stress interpolation polynomials at each side of the crack. 19Later, this technique was used to evaluate upper error bounds 42 and to create a recovery based goal-oriented error estimator for XFEM approximations. 43While these works were modifications/extensions of methods proposed for standard FEM, other procedures were also introduced exclusively for X/GFEM.Prange et al. 44 developed a recovery procedure based on a global least squares projection of raw stresses calculated from X/GFEM approximations in fracture problems with arbitrarily distributed cracks and inelastic material behavior.Lins et al. 45 later adapted this approach to improve the stress field within stable GFEM (SGFEM); because SGFEM modifies the X/GFEM enrichment functions to solve the issue of ill-conditioned stiffness matrices, [46][47][48] the enrichment functions used to interpolate the recovered stress field are also modified.Instead of using a global least squares projection, Lins et al. 49 proposed a more computationally efficient recovery technique, whereby a consistent block-diagonal projection operator 50 is used to perform a locally weighted least squares projection of directly-calculated stresses over patches of elements.Sharma and Maute 51 introduced the ghost penalty method 52 into X/GFEM for penalizing the jump in displacement gradients across element boundaries near material interfaces; this method is later combined with a smoothing technique defined over the entire domain to obtain a smooth stress field.
Contrary to X/GFEM, where enrichments are associated to nodes of the original FE discretization, interface-and discontinuity-enriched finite element procedures seek to enhance the approximation by associating enrichments to nodes that are created directly along discontinuities.The Interface-enriched Generalized Finite Element Method (IGFEM), 7 the Hierarchical Interface-enriched Finite Element Method (HIFEM), 53,54 and the Discontinuity-Enriched Finite Element Method (DE-FEM) 8,9,55 thus provide an alternative enriched approach for solving problems with discontinuities.These procedures decouple the background mesh from discontinuities as in X/GFEM, but they do not have many of the latter's drawbacks: by creating enriched degree of freedoms (DOFs) at the intersections between the background mesh and discontinuities, these EFEA techniques are endowed with intrinsic stability with regard to the condition number of stiffness matrices, 54 no loss of accuracy in blending elements (neighboring to cut elements) due to the locality of enrichment functions, 8 straightforward enforcement of nonhomogeneous essential boundary conditions, 56 and smooth recovered traction profiles from Dirichlet boundaries. 56,57Nevertheless, because these methods recover the standard FEM space by means of an enrichment, they also suffer from inaccurate stress fields caused by sliver elements; it has been reported that stresses can be overestimated more than 150%. 58n this article, we propose a stress recovery technique for EFEA based on SIP, which as stated above has several advantages over other recovery techniques.SIP has been shown to enhance stress fields from displacement-based FEM results, with higher convergence rates than those of directly-calculated stresses.While the method is generally applicable and thus could also be used with X/GFEM, we demonstrate the procedure on weak discontinuities resolved by IGFEM.When solving multiphase material problems, only elements with the same material properties are considered when constructing the patch of elements for stress recovery, as initially suggested by Payen and Bathe. 29The performance of the proposed method is studied by means of three numerical examples.With Eshelby's inclusion problem we explore the accuracy of stress fields obtained using standard FEM with fitted meshes and IGFEM with unfitted meshes.The recovery technique coupled to IGFEM provides more accurate stresses, for the same mesh size, than directly-calculated stresses obtained by standard FEM on fitted meshes.The second example studies an elliptical cavity in an infinite plate under remote loading, where we first adjust the ratio between the elliptical axes.The corresponding results show that the proposed recovery technique generally provides more accurate approximations than the direct calculation.Later, we consider a circular cavity under the same boundary condition, where the radius of the circle is adjusted to create different enriched discretizations for the same structured background mesh.In addition to evaluating elemental stresses, we also calculate nodal stresses obtained via the proposed method and other smoothing procedures based on directly-calculated stresses.The proposed stress recovery technique performs best for evaluating the stress concentration factor K t , which is the ratio of the highest stress to the nominal far field stress.Furthermore, a convergence study shows the reliability of the proposed recovery technique with regards to mesh size.Finally, a pressurized sphere example demonstrates that the proposed recovery technique can also avoid overestimated stresses in tiny integration elements in 3D.

IGFEM-based analysis
Consider an open domain Ω ⊂ R d with closure Ω that is composed of two disjoint regions Ω 1 and Ω 2 , which are occupied by different phases.For simplicity, the 2D case is shown in Figure 1.We denote by u the displacement field and by u i its restriction to the ith domain Ω i , that is, u i = u| Ω i .A similar notation is used for other subscripted quantities.The domain's boundary Γ ≡ Ω = Ω ⧵ Ω is smooth and consists of two disjoint regions Γ N and Γ D , where surface tractions t and displacements u are prescribed, respectively.We denote by n the outward unit normal to the boundary Ω.The material interface between domains Ω 1 and Ω 2 is denoted as Γ 12 and has unit outward normal n 12 .The strong form of the elastostatics boundary value problem is: find the displacement field u such that

F I G U R E 1
with interface conditions where ∇⋅ denotes the divergence operator,  is Cauchy's stress tensor, and b the body force.We assume a linear elastic material behavior and thus  = C ∶ (u), where C is the elastic modulus tensor, and  = 1 2 (∇u + ∇u ⊺ ) the infinitesimal strain tensor.The weak form of the boundary value problem is: where  and  are vector-valued function sets defined as with  2 denoting the Hilbert space of square-integrable functions and  1 the first-order Sobolev space, respectively.The bilinear a(u, v) and linear (v) forms are given, respectively, by and In order to solve the problem above, Ω is discretized by a finite element mesh-not necessarily fitting to material interfaces-such that the discretized domain Ω h = ⋃ E i=1 e i , with e i denoting the ith finite element and E the total number of elements.Then, the finite-dimensional form of Equation ( 2) is expressed as where u h ∈  h ⊂  and v h ∈  h ⊂  are trial and test functions, respectively.In IGFEM, the displacement field u h is expressed as where the first term is the standard FEM approximation, with  h denoting the index set of all standard nodes of the background mesh (dark circles in Figure 1), associated with standard shape functions N i and degree of freedoms (DOFs) U i ; the second term represents the enrichment, where  w is the index set of enriched nodes that are associated with enrichment functions  i and enriched DOFs  i .Enriched nodes (red circles in Figure 1) are created at the intersections between element edges and material interfaces.Cut background elements are split into integration elements, which as the name suggests, are used for the numerical quadrature of local stiffness and force vector arrays; in addition to be used for constructing enrichment functions  i (by means of Lagrange shape functions) and to ensure that the lowest number of integration points is used for their numerical quadrature (since  i ∈ C 0 ), integration elements are also helpful in the post-processing stage to properly visualize the primal field. 9Enrichment functions attain their maximum value at the location of enriched nodes and decrease linearly to zero at other nodes in the cut elements (see Figure 2).Background elements that are not cut by discontinuities follow the standard FEM procedure for constructing their corresponding stiffness matrix k e and force vector f e .For integration elements, these arrays are computed as B ⊺ CB de, and where is the strain-displacement matrix, N and  are the element's standard shape functions and enrichment functions, respectively, and the differential operator  defined in 2D or 3D is given, respectively, by or The global stiffness matrix K and force vector F are then obtained by where A is the standard finite element assembly operator.For more details on IGFEM's formulation see References 7,54.
IGFEM not only retains the most salient feature of X/GFEM-decoupling the background mesh from discontinuities-but also keeps the attractive properties of standard FEM: As enrichment functions are constructed with Lagrange shape functions of integration elements, they are exactly zero at original mesh nodes; DOFs associated with these therefore keep their physical meaning.Standard shape functions retain the Kronecker- property, that is, N i (x j ) =  ij for any standard node x j .Furthermore, prescribing nonzero essential boundary conditions along discontinuities is as straightforward as in the standard FEM, 56 with smooth recovered tractions (reactions). 56,57,59The computer F I G U R E 2 Enrichment function  5 created with the aid of Lagrange shape functions in integration elements, which are conforming to the material interface (marked with a red curve); the function is nonzero only in cut elements e 1 (with nodes x 1 , x 2 , and x 4 ) and e 2 (with nodes x 2 , x 3 , and x 4 ).The enrichment function, which is associated with enriched DOFs  5 , attains a maximum unit value at node x 5 and decreases linearly to zero at all other nodes implementation is therefore considerably simpler than that of X/GFEM since it only requires a few modification to a displacement-based FE code. 9Most importantly, intrinsic stability is attained in IGFEM by scaling enrichment functions with a proper scaling factor as interfaces approach mesh nodes. 54With such scaling the rate of growth of the condition number is the same as that of standard FEM on fitted meshes, that is, (h −2 ), where h is the mesh size.Furthermore, even without such scaling factor-which causes the condition number of the global stiffness matrix to h grow unbounded as interfaces get arbitrarily close to mesh nodes-a simple Jacobi-like preconditioner recovers stability. 54

Stress improvement procedure (SIP)
SIP is derived from a mixed formulation based on the Hu-Washizu principle. 60It relaxes the stress-strain relationship point-wise, but enhances the fulfillment of equilibrium over the patch.Following the procedure proposed by Payen and Bathe, 29 the recovered stress  e of the eth element should satisfy the equilibrium in weak form over the calculation domain  (i.e., a patch of elements), and the projection of the difference between enhanced and directly-calculated stresses onto the space of self-equilibrated stresses should be zero.To wit, and where denoting by P k is the space of complete polynomials of degree k in the patch,  e belongs to the vector-valued space [P 1 ] 2 , σe is any element in the tensor-valued space of the self-equilibrated stress defined as and  h e is the directly-calculated stress in the eth element, that is,  h e = CBU e , where U e is the local element DOF vector.As shown in Figure 3A, for evaluating the recovered stress of the hatched element, a patch of elements (in darker shade) is considered as the calculation domain.The quantities above can be interpolated as where σ, ζ, and σ are coefficient vectors, and E  , E  , and E  are interpolation matrices as defined in Appendix A. Noteworthy, a complete polynomial of degree 2 is used to construct the interpolation matrices E  and E  , thereby having an expected order of convergence (h 2 ), where h is the mesh size. 29Equations ( 13) and ( 14) can then be expressed as which is solved for coefficient vector σ of the target element, after which its corresponding recovered enhanced stress is simply  e = E  (x) σ.Gauss quadrature is used to evaluate Equation (16).Considering the interpolation matrices and the differential operator, the highest-order term is found in E ⊺  E  .Since the maximum polynomial order is quartic, six integration points are then used to exactly integrate the integrands.
A nodal stress field can be constructed by recovering stresses in all elements and averaging their values at nodes.In order to obtain the recovered stress  j at the jth node, a nodal patch  j containing N j elements sharing the node is first assembled (see Figure 3B).Then, the nodal stress  j is simply computed as the average, that is, where  e (x j ) = E  (x j ) σe is the stress of the eth element evaluated at the jth node.
In IGFEM-based analysis, elements at either side of material interfaces are assigned different material properties.Therefore, the calculation domain considers only a patch of elements with the same material properties.Figure 4 shows three cases of element patches used for recovering the stress near a material interface (marked in red) in both 2D and 3D, where the union of elements connected to a target element is set as the stress calculation domain.A similar strategy is used to obtain the enhanced nodal stress (see Figure 5).There is a special scenario for recovering the stress of a node at a corner, where a larger patch of elements is used for making the nodal stress distribution smoother when constructing the calculation domain (see Figures 5A and 5D for 2D and 3D, respectively).Figures 5C and 5F show a situation where the patch of a node, which is not located along the material interface, is composed of both background and integration elements.It is worth noting that the stress calculation domain (the patch) is tied to the finite element discretization, and therefore refining the original background mesh creates a discretization where discontinuities cut more background elements; however, the configuration of the element patch used to recover the stress distribution does not change for structured meshes regardless of the mesh size.

Alternative smoothing formulations
We compare the proposed stress recovery technique for both elemental and nodal stresses to other two stress smoothing procedures.The first one averages directly-calculated stresses on a target element or node considering a patch; the averaged stress is computed as where N e is the number of the element in the patch, which is taken to be the same as that used for the proposed recovery technique.We use  h e and  h i to denote the averaged stresses of the eth element and ith node, respectively.We also consider an area-weighted average stress, computed as where A e is the area of the eth element.Similarly,  h⊳ e and  h⊳ i denote the area-weighted averaged stresses of the eth element and ith node, respectively.

NUMERICAL EXAMPLES
In this section three examples are investigated to demonstrate the effectiveness of the proposed approach.The finite element analysis is performed under plane strain conditions, and no units are used so results can be interpreted in any consistent unit system.We perform a convergence study for both examples to investigate the accuracy of the proposed approach with increasingly finer meshes, and we use the H 0 -norm of the stress to quantify the error: where  ex is the exact analytical stress.This global error measure is evaluated by summing up the contribution of all standard and integration elements.For the error in directly-calculated stresses,  e is replaced by  h e .

Eshelby's inclusion problem
As shown in Figure 6, a circular inclusion with radius r i is embedded into a matrix with radius r o = 2.The mismatch in material properties at the interface between matrix and inclusion is responsible for the discontinuous gradient field.Young's moduli and Poisson ratios for the inclusion and matrix are, respectively, E 1 = 1,  1 = 0.25 and E 2 = 10,  2 = 0.3.Dirichlet boundary conditions u r = r o and u  = 0 are prescribed on the matrix's outer boundary (along r = r o ).The exact solution for the displacement field in polar coordinates is given by 10,56,59 where the function f (r) = In this equation, C is a function of material properties: where  i ,  i , i = 1, 2 are the Lamé constants, which can be obtained as a function of E i and  i as The stress field for this problem is given by where  ij is the Kronecker delta, i, j ∈ {r, } and  kk is the trace of the strain tensor ( rr +   ), and the exact strain field for this problem is given by A 2 × 2 square computational domain is chosen, as shown in Figure 6, and the domain is discretized by a finite element mesh composed of 60 × 60 × 2 constant strain triangular elements.The exact displacement given in Equation ( 21) is prescribed along the square's boundary.For studying the performance of the proposed recovery technique, the internal radius r i is set to increase from r i = 0.35 to r i = 0.42 with step size Δr i = 0.035, which ensures the creation of integration elements with bad aspect ratios and/or tiny areas along the interface.Furthermore, we compare results with standard FEM using conforming meshes, where cut elements are replaced by standard FEM elements with the same geometry as integration elements in IGFEM.][63] Figure 7A shows the maximum von Mises stress as a function of the internal radius r i .Directly-calculated and recovered stresses for standard FEM and IGFEM are compared to the analytical solution.It can be seen that IGFEM recovered stresses are more accurate than directly-calculated ones for both standard FEM on matching meshes and IGFEM with a fixed unfitted mesh.As IGFEM can fully recover the approximation space of standard FEM, curves of directly-calculated and recovered stresses obtained under standard FEM on the same conforming discretizations overlap with those of IGFEM.The stresses closest to the analytical solution are obtained by the recovery technique applied to standard FEM on fitted meshes at the expense of losing the versatility of using a mesh that is fully decoupled from the interface.Figure 7B further compares the proposed procedure to the other two smoothing techniques discussed in Section 2.3.It is shown that the maximum von Mises stress computed with averaged stress fields given by Equations ( 18) and ( 19) is worse than that obtained by the proposed approach.The von Mises stress distributions for r i = 0.3605 are given in Figure 8.The fields for directly-calculated and recovered stresses using standard FEM on good-quality fitted meshes are given in Figures 8A and 8B.Figures 8C-F show the results for IGFEM-based analysis.Figure 8C displays an element along the material interface with the peak stress, which is caused by the tiny elemental area.However, the recovery technique proposed removes the stress overestimation and makes the stress distribution smoother around that region (see Figure 8D).The von Mises stress fields obtained with the averaged and area-weighted elemental stresses are displayed on Figures 8E and 8F, which show the averaging equations yield an underestimated stress.Figures 8G and 8H show von Mises distributions evaluated with directly-calculated and recovered stresses with CDFEM, which are the same to those of IGFEM.F I G U R E 9 H 0 -norm of the error in stress as a function of mesh size h with r i = 0.371.The curves for recovered stress  e and directly-calculated stress  h e show that the former is not only more accurate for any given mesh size, but also that it converges at a faster rate Finally, four background meshes with 30 × 30 × 2, 60 × 60 × 2, 120 × 120 × 2, and 240 × 240 × 2 linear triangular elements are used for the convergence analysis with r i = 0.371.Figure 9 shows the global error defined in Equation ( 20) as a function of mesh size h.It can be seen that the recovered stress  e convergences faster (with a rate of 1.20) than the directly-calculated stress  h e .The recovered stress with a coarse mesh can reach the same level of accuracy of the directly-calculated stress with a refined mesh.

An elliptical cavity in an infinite plate under remote loading
Figure 10 shows a traction-free elliptical hole in an infinite domain that is subjected to a distant tensile stress  0 in the x direction, where the elliptical axes 2a and 2b are aligned with the x and y coordinate axes, respectively.The analytical stress field for this problem is given by: 64 5

F I G U R E 10
Schematic of an elliptical cavity with semiprincipal axes a and b in an infinite plate under remote loading  0 in the x direction.For the analysis, only the shaded area is considered due to symmetry (with appropriate symmetry boundary conditions) and exact tractions obtained from Equation ( 26) are applied on the right and top sides where ] , ] , The outward unit normal vector n is expressed as where ) .
Due to symmetry, only the 5 × 5 region (darker shade) is considered as the analysis domain, and is discretized with a mesh of 60 × 60 × 2 constant strain triangles.Symmetry boundary conditions are then imposed along left and bottom edges, and tractions obtained from the analytical equation with  0 = 1 are prescribed on top and right sides.In order to avoid an ill-conditioned stiffness matrix, a material with Young's modulus 10 −9 is assigned to the void part; the solid part has Young's modulus E = 1 and Poisson's ratio  = 0.3.For this problem we first adjust the ratio between the semimajor axis a and semiminor axis b to investigate the performance of the proposed recovery technique, where a and b are chosen from the set of values {2.49, 1.245, 0.498, 0.249, 0.1245}, resulting in a∕b ratios ranging from 0.05 to 20.In addition, we solve the problem on three different discretizations, that is, 60 × 60 × 2, 120 × 120 × 2, and 240 × 240 × 2 triangular elements with corresponding mesh sizes h = 0.1348, 0.0674, and 0.0337, respectively.For simplicity, the maximum von Mises stress  vM obtained via direct calculation and the recovery technique is normalized by the analytical result for the corresponding ratio.
Figure 11 shows the normalized von Mises stress values as a function of a∕b, where it can be seen that the proposed recovery technique generally provides more accurate approximations than the direct calculation.However, as the ratio a∕b is decreased, the cavity progressively morphs into a crack perpendicular to the remote loading.As a result, the stress field becomes increasingly singular, and both directly-calculated and recovered stresses do not behave well.The corresponding errors in the H 0 norm are also given in the figure, where the difference between the recovery technique and the direct calculation becomes small for small and large values of a∕b.
We now set a = b for a circular cavity in an infinite plate with the same boundary conditions.The analyses are performed with different values of radius a, increasing from a = 1.895 to a = 2.095 with step Δa = 0.01. Figure 12 shows the stress concentration factor K t -that is, the ratio of the highest stress to the nominal far field stress 65 -evaluated by means of Equations ( 17), (18), and (19) for different radii a.For this problem K t = 3 (Equation ( 26) for (x, y) = (0, a)) and thus it can be seen that all numerical values fluctuate around this value, since different values of radius a also imply different finite element discretizations.Moreover, the interface approaches a standard node under certain values of the radius, leading to tiny integration elements with peak von Mises directly-calculated stresses.The curve of the recovered stresses, however, does not show peak values and is smoother with respect to variations of the radius.Therefore, we conclude that the proposed recovery method is more robust with respect to changes in the location of interfaces.Averaged and area-weighted averaged stresses show large differences between the maximum and minimum values, and thus they are greatly effected by the discretization.For instance, there is a huge drop in the stress concentration factor as a increases from 1.915 to 1.925.With a = 1.915, Figure 13A shows the element patch (in darker shade) used to evaluate K t at the location of a node (red circle).Since patch elements are close to the expected coordinate, they provide an accurate approximation of the stress concentration factor.As the radius a increases to 1.925, a larger patch is considered as the calculation domain (see Figure 13B).However, considering more distant elements in the patch, which have lower stress, yields an underestimated stress concentration factor.The proposed recovery technique predicts the best approximation of the analytical value as its curve is the most stable one.Notice that the same patches are used in the proposed recovery technique, which is therefore less sensitive to changes in the discretization; Figures 13C and 13D show the same configuration of element patches within a refined mesh used to recover the stress of nodes shown in Figures 13A and 13B.Furthermore, we perform a convergence study by using background meshes with increasingly smaller mesh size h to investigate the behavior of each formulation.In addition to the mesh with 60 × 60 × 2 elements used earlier, we consider extra meshes with 30 × 30 × 2, 120 × 120 × 2, and 240 × 240 × 2 linear triangular elements.For a given mesh size, minimum, maximum, and average values were obtained for the same radii range, that is, a = [1.895,2.095].Results are reported in Figure 14, where the curves correspond to average values and the shades span from minimum to maximum values of the stress concentration factor K t .While all formulations approach the analytical value as the mesh is refined, the stress recovery technique shows the most accurate prediction, as its shaded area is the smallest one among all formulations.With the same background meshes, we also calculate the global H 0 -norm of the error for a = 1.895.Note that this study looks again at elemental stresses for computing Equation (20). Figure 15 shows that the convergence rate of the recovered stress is 1.34, which is in agreement with the value obtained when applying the SIP recovery technique on standard FEM tetrahedral elements.F I G U R E 15 H 0 -norm of the error in stress as a function of mesh size h for a = 1.895.As in Figure 9, the recovered stress  e converges faster than the directly-calculated stress  h e , this time with a higher convergence rate of 1.34

Pressurized sphere
A hollow sphere with internal radius a = 1 and external radius b = 2 is considered in this example, where a uniform pressure p = 1 is applied to the internal surface (see Figure 16A).The analytical displacements and stresses for this problem are given by: 66 respectively, where r = √ x 2 + y 2 + z 2 is the distance from the sphere center,  = arctan(y∕x), and  = arccos(z∕r).Due to symmetry, we only consider one eighth of the sphere for the analysis (see Figure 16B); consequently, symmetric boundary conditions are prescribed on the planar surfaces in the x, y, and z directions.The sphere portion is then immersed into a cubic domain with dimensions 2.25 × 2.25 × 2.25, which is discretized by a structured mesh composed of 10 × 10 × 10 × 6 linear tetrahedral elements.
According to the analytical solution, the maximum value of von Mises stress is  vM = 1.71.However, the maximum von Mises stress obtained with directly calculated stresses is  vM = 2.11 (see Figure 17A), where peak stresses are the result of several tiny integration elements (marked in gray).Figure 17B shows that the recovered stresses are much closer to the exact solutions.The corresponding error distributions are shown in Figures 17C and 17D, where the maximum value of the latter is one order of magnitude smaller than that of the former.In addition, background meshes with 5 × 5 × 5 × 6, 10 × 10 × 10 × 6, 20 × 20 × 20 × 6, and 40 × 40 × 40 × 6 tetrahedral elements are used to perform a convergence study, and we compute the H 0 -norm of the error for both directly-calculated and recovered stresses.Figure 18 shows a much faster convergence rate for recovered stresses (1.46 vs. 0.89).The convergence rate is on par with the value 1.49 obtained by the recovery technique applied on standard FEM using fitted meshes by Sharma et al.

SUMMARY AND CONCLUSIONS
In this article we extended the SIP proposed by Payen and Bathe 29 to recover the stress field in problems with discontinuous gradient fields, whose solution is obtained by means of an enriched finite element analysis.We investigated the proposed recovery technique in conjunction with IGFEM, although we foresee no further developments when coupling the technique to other EFEA such as X/GFEM.In IGFEM material interfaces subdivide mesh elements into subdomains to which different material properties are assigned.In order to smooth the stress field close to the interface, the calculation domain used is carefully constructed by selecting a patch of (integration) elements with the same material properties.We also used the recovery technique to evaluate nodal stresses along discontinuities, where only integration elements at the same side of discontinuities are used as the element patch.For corner nodes along the discontinuity, a larger patch of elements that involves standard uncut elements is used as the calculation domain.
The technique was demonstrated by means of 2D and 3D numerical examples.With Eshelby's inclusion problem, we showed that recovered IGFEM stresses are more accurate than directly-calculated ones for both standard FEM and IGFEM.Therefore, for IGFEM the stress recovery technique provides an elegant means to approximate the stress field compared to standard FEM on fitted meshes.In addition, an elliptical cavity in an infinite plate under remote loading was used to investigate the performance of the recovery method while changing the ratio between semimajor and semiminor axes.As the ratio is decreased the cavity morphs into a crack, leading to an increasingly singular field and consequently to a less accurate recovered stress field.Later, the circular cavity example was used to investigate the recovery of nodal stress fields.Averaging and area-weighted averaging smoothing techniques were also used to evaluate the nodal stresses.In comparison, the proposed recovery technique does not only provide a more accurate result for the value of the stress concentration factor, but is also less sensitive to the choice of discretization (see Figure 12).For the 3D pressurized sphere example, it was shown that the proposed recovery technique can avoid overestimated stresses in tiny integration elements.More importantly, the convergence rate of the recovered stresses associated with mesh size h is close to the value 1.49 obtained from standard FEM using the recovery approach on fitted meshes by Sharma et al. 31 Note that although in all examples we used analytical solutions, the proposed methodology could also be used for a posteriori error estimation when solving problems without closed-form solutions (similarly to other recovery approaches such as SPR).
Even though this work focused on static material interfaces, the proposed procedure does not suffer any modification for problems with evolving interfaces, for example, in stress-based topology optimization, for which an accurate approximation of the stress field is paramount. 67In fact, the proposed method has been used to recover the stress field in an incoming article that designs structures with tailored brittle fracture resistance, 68 which hinges on an accurate evaluation of energy release rates-which in turn greatly depend on the approximation of the stress field.
As the proposed recovery technique aims to improve the stress approximation for a linear field using a quadratic interpolant in a target element, future work may develop the proposed methodology to recover stress fields of higher order.This could be done via increasing the order of the approximation-for both standard and enrichment functions-and using an order higher interpolant to recover the stress field.The proposed technique can also be extended straightforwardly to recover stress fields in problems with strong discontinuities such as fracture, for which DE-FEM has been proposed recently; 8,9 DE-FEM is a generalization of IGFEM for the seamless treatment of both weak and strong discontinuities with a unique formulation, and inherits all of IGFEM's advantages.Nevertheless, the DE-FEM formulation does not use singular enrichments so a modification would be required to the enriched formulation if accurate singular stress fields are sought.For instance, following the work of Duflot and Bordas, 35 the stress field could be decomposed into smooth, nonsmooth nonsingular, and nonsmooth singular components.While the first two could be handled using the same strategy mentioned in this paper, the last one could use an analytical formulation based on stress intensity factors that describes the asymptotic fields near the crack tip.This could greatly improve the recovered stress at crack fronts.
As a final note, although the proposed recovery approach can greatly improve the accuracy of the stress field obtained with IGFEM, it is just a post-processing technique that does not solve the fundamental issue at its core: Because IGFEM recovers-through an enriched procedure-the same standard finite element space, the recovery of gradient fields is still susceptible to the way elements are cut by interfaces.Therefore, addressing this issue at the root would require for alternative means to create the enriched finite element space.and A domain Ω consists of two parts Ω 1 and Ω 2 with a smooth boundary Ω = Γ D ∪ Γ N .Dirichlet boundary conditions are prescribed on Γ D and surface tractions on Γ N .For the discretized model, enriched nodes (marked in red) are created at the intersections between material interfaces (red segments) and element edges; integration elements (red triangles) are also created near the interface

F
I G U R E 3 (A) 2D Patch of elements (darker shade) used to recover the stress of the target element (hatched); (B) Nodal patch (darker shade) used to compute the enhanced average stress at a node (red circle)

4 5
Computational domain used to compute the enhanced element stress considers only elements on the same side of the interface (shown in red).(A-C) 2D element patches shown in darkest shade; (D-F) 3D element patches shown as opaque Computational domain used to compute the enhanced node stress considers only elements on the same side of the interface (shown in red).(A-C) 2D element patches shown in darkest shade; (D-F) 3D element patches shown as opaque

ν 2 2 2F I G U R E 6 F I G U R E 7
Schematic of Eshelby's inclusion problem with outside radius r o = 2 and inside radius r i .For the numerical analysis, a square domain of size 2 × 2 (darker shade) is considered, with the analytical displacement u prescribed on the boundary.For the material properties, E 2 ∕E 1 = 10 Maximum von Mises stress associated with the internal radius r i .(A) Stress evaluated by standard FEM with fitted meshes, IGFEM with a structured background mesh, and CDFEM with the conforming discretizations created in IGFEM; (B) Stress obtained with the averaged and area-weighted smoothing formulations using IGFEM

F
I G U R E 8 von Mises stress fields for an internal radius r i = 0.3605: (A,B) Standard FEM on a fitted mesh; (C-F) IGFEM on an unfitted mesh; (G,H) Standard FEM on a conforming discretization created using IGFEM's integration elements (CDFEM); The fields correspond to directly-calculated (A,C,G), recovered (B,D,H), averaged (E), and area-weighted (F) stresses

F I G U R E 11
Maximum normalized von Mises stress  vM (top) and H 0 -norm of the error in stress (middle) obtained with directly-calculated and recovered stresses as a function of a∕b.Three different discretizations are used with 60 × 60 × 2, 120 × 120 × 2, and 240 × 240 × 2 elements.The corresponding mesh sizes are h = 0.1348, 0.0674, and 0.0337, respectively.The cavity configurations for different values of a∕b are shown at the bottom a K t F I G U R E 12 Stress concentration factor K t associated with radius a that is evaluated by averaged and area-weighted average stresses, recovery technique, and analytical solution, respectively

F
I G U R E 13 Element patches (darker shade) used to evaluate the stress concentration factor K t at the location of node (red circle).a = 1.915 for (A,C) and a = 1.925 for (B,D).The patches' shape remains the same as meshes are refined h K t F I G U R E 14 Stress concentration factor K t as a function of mesh size h.In this convergence study, meshes with 30 × 30 × 2, 60 × 60 × 2, 120 × 120 × 2, and 240 × 240 × 2 linear triangular elements were used.For each mesh size, the figure reports minimum, average, and maximum values of K t for radii in the range U R E 16 (A) Schematic of an internally pressurized sphere with internal radius a and external radius b where the uniform pressure p is imposed on the internal surface; (B) Due to symmetry, an 1/8 sphere is considered and immersed into a background mesh with 10 × 10 × 10 × 6 linear tetrahedral elements and

31
U R E 17 (A,B) von Mises stress fields obtained with (A) directly-calculated and (B) recovered stresses, where the former shows several tiny integration elements (marked with red) with peak stresses and the latter displays a more smooth distribution.(C,D) Element-wise error distributions in the H 0 -norm under (A) directly-calculated and (B) recovered stresses, where the maximum value of the latter is one magnitude smaller than that of the former F I G U R E 18 H 0 -norm of the error in stress as a function of mesh size h, where the recovered stress  e converges faster than the directly-calculated stress  h e , this time with a higher convergence rate of 1.46 How to cite this article: Zhang J, Aragón AM.An improved stress recovery technique for the unfitted finite element analysis of discontinuous gradient fields.Int J Numer Methods Eng.2022;123(3):639-663. https://doi.org/10.1002/nme.6825APPENDIXA. EXPRESSION OF E  , E  , AND E Interpolation matrices E  , E  , and E  used in Equation (16) for recovering stress distributions are given by 3D, respectively, where e  = [1 x y z xy yz zx x 2 y 2 z 2 ].