Linear smoothed extended finite element method
Abstract
The extended finite element method (XFEM) was introduced in 1999 to treat problems involving discontinuities with no or minimal remeshing through appropriate enrichment functions. This enables elements to be split by a discontinuity, strong or weak and hence requires the integration of discontinuous functions or functions with discontinuous derivatives over elementary volumes. Moreover, in the case of open surfaces and singularities, special, usually nonpolynomial functions must also be integrated. A variety of approaches have been proposed to facilitate these special types of numerical integration, which have been shown to have a large impact on the accuracy and the convergence of the numerical solution. The smoothed extended finite element method (SmXFEM) bordas2010strain (), for example, makes numerical integration elegant and simple by transforming volume integrals into surface integrals. However, it was reported in bordas2010strain (); bordas2011performance () that the strain smoothing is inaccurate when nonpolynomial functions are in the basis. This is due to the constant smoothing function used over the smoothing domains which destroys the effect of the singularity. In this paper, we investigate the benefits of a recently developed Linear smoothing procedure francis2015linear () which provides better approximation to higher order polynomial fields in the basis. Some benchmark problems in the context of linear elastic fracture mechanics (LEFM) are solved to compare the standard XFEM, the constantsmoothed XFEM (SmXFEM) and the linearsmoothed XFEM (LSmXFEM). We observe that the convergence rates of all three methods are the same. The stress intensity factors (SIFs) computed through the proposed LSmXFEM are however more accurate than that obtained through SmXFEM. To conclude, compared to the conventional XFEM, the same order of accuracy is achieved at a relatively low computational effort.
keywords:
Smoothed finite element method, linear smoothing, extended finite element method, numerical integration, fracture mechanics.1 Introduction
The major difficulty associated with solving problems involving evolving discontinuities is the meshing and remeshing required as the discontinuities evolve in time. This difficulty is exacerbated when singularities are also present, as is the case in crack growth simulations. These difficulties are somewhat alleviated by the introduction of enrichment functions to represent the discontinuities and the singularities at the patch level, in finite elements or meshfree methods. A first approach to treat discontinuities without an explicit meshing was proposed as early as 1995 in oliyer1995continuum (). A much more versatile approach was presented a few years later in the form of the extended finite element method (XFEM) belytschko1999elastic () dolbow1999finite () by exploiting the partition of unity property identified by Melenk and Babuška melenk1996partition (). Partition of unity enrichment for problems with discontinuous solutions is now widely used both in academia and in industrial practice bordasmoran2006 (); bordasconley2007 (); wyartcoulon2007 () and is known under various names, including the generalized finite element method (GFEM) duartehamzeh2001 (); babuskabanerjee2012 () and the extended finite element method (XFEM). The approach has also been widely used in the form of enriched meshfree methods rabczukbordas2007 ().
Another problem associated with partition of unity methods involving nonpolynomial basis functions is to integrate the resulting fields accurately. These enriched methods, also carry along the element mapping involved in building the system matrices. The regularity and positive definiteness of the isoparametric mapping poses a number of restrictions on the allowable shapes of the finite elements: for example, the element should be convex. Meshfree methods also have to face such problems associated with the regularity, distortion and clustering in the point cloud. Under large distortions, meshfree methods face numerical instabilities and low accuracy belytschkoguo2000 (). Nodal integration also leads to instabilities in cases where higher order shape functions are used. This is due to the fact that in the meshfree methods each node would be associated with a support domain. And the shape functions derivatives would be integrated in this support domain. This implies that each integration domain would be associated with only one integration point (i.e the node). Hence when only one integration is point is considered for higher order functions (other than constant strain) the integration would be similar to the inadequate reduced integration which in turn causes instabilities.
To alleviate these instabilities, the strain smoothing concept was introduced for meshfree methods chen2001stabilized (). The basic idea of strain smoothing is to transform numerical integration over volumes to integration over surfaces, thereby removing instabilities due to the evaluation of the shape function derivatives at the nodes. This approach was later extended to finite element methods by Liu et al liu2007smoothed (). The resulting method was coined the Smoothed finite element method (SFEM), was discussed in a number of papers bordas2010approximation (); bordas2011performance (); liu2009edge (); liu2009node (); nguyenrabczuk2008 (); nguyen2008smooth () and applied to a wide variety of problems. The SFEM reduces the mesh sensitivity to some extent by avoiding the necessity of evaluating the Jacobian. Since the derivatives are not needed, the isoparametric mapping is also avoided.
The SFEM involves computation of a smoothed strain from the standard compatible strain field. The smoothed strain is evaluated as a spatial average of the standard strain field over smoothing subcells which cover the domain andthat can be chosen independently from the mesh structure. These smoothing cells are typically constructed from the mesh following different approaches. This gives rise to a number of methods including cellbased SFEM (CSFEM) liu2007smoothed (); dai2007n (); nguyen2008smooth (); bordas2010approximation (), nodebased SFEM (NSFEM) liu2009node (), edgebased SFEM (ESFEM) liu2009edge (), facebased SFEM (FSFEM).
The method was later extended to solve problems with field discontinuities, both strong and weak, by Bordas et al bordas2010strain (). This was achieved by extending strain smoothing to the partition of unity framework babuvska1997partition (); melenk1996partition (). Though the smoothed FEM did make the integration procedure more elegant, it was also reported in bordas2011performance () that the error levels are higher due to the inaccurate approximation of the near tip singular fields. Similar errors were also encountered in higher order elements and polygons natarajan2014convergence (). It is noteworthy that similar difficulties are also present in meshfree methods, as addressed in duan2012second () by employing the Discrete Divergence Consistency (DDC) requirement by establishing a compatibility relation between the shape function and its derivatives. This was recently extended for the FEM in francis2015linear () and named: Linear smoothing (LS) scheme. It was also reported that the linear smoothing scheme provides an improved accuracy compared to the standard constantbased smoothing of the SFEM.
The present paper aims at investigating how the use of higherorder smoothing, in particular linear smoothing, resolves the lack of accuracy observed when constant smoothing is employed with nonpolynomial bases functions. The paper is organized as follows. After presenting the governing equations for elastostatics, a brief overview of the extended finite element method is given in Section 2. Section 3 presents the linear smoothing technique. A few standard benchmark problems in linear elastic fracture mechanics, solved by using the developed method are presented and the accuracy, convergence and the efficiency of the proposed method are discussed in Section 4, followed by concluding remarks in the last section.
2 Theoretical formulation
2.1 Governing equations for elastostatics
Consider a linear elastic body as shown in Figure 1, with a discontinuity. Let the domain be bounded by . In this case the boundary has three parts namely , where the displacement boundary conditions are applied, , where the tractive boundary conditions are applied and , which is the traction free surface representing the discontinuity, such that and .
The governing equation to be solved is
(1) 
The boundary conditions are as follows
(2a)  
(2b)  
(2c) 
where, is the gradient operator, is the Cauchy stress tensor, is the body force per unit volume, is the unit outward normal and is the applied tractive stress. For a body undergoing small displacements and strains, the strain displacement equation reads as
(3) 
where, is the symmetric part of the gradient operator. By substituting the constitutive relations and the straindisplacement relations the weak form of the above Equation (1) becomes Equation(4) in the absence of the body forces: find such that
(4) 
where, and are the trial and the test function spaces, respectively. For the aforementioned problem, the function spaces includes functions that are discontinuous across .
where the space includes linear displacement fields. The domain is partitioned into elements , and on using shape functions that span at least the linear space, we substitute vectorvalued trial and test functions and , respectively, into Equation (4) and apply a standard Galerkin procedure to obtain the discrete weak form: find such that
(6) 
which leads to the following system of linear equations:
(7a)  
(7b)  
(7c) 
where is the assembled stiffness matrix, the assembled nodal force vector, the assembled vector of nodal displacements, is the matrix of shape functions, is the constitutive matrix for an isotropic linear elastic material, and is the straindisplacement matrix that is computed using the derivatives of the shape functions.
2.2 eXtended Finite Element Method
With the regular FEM, the mesh topology has to conform to the discontinuous surface. The introduction of the XFEM has alleviated these difficulties by allowing the discontinuities to be independent of the underlying discretization. Within the framework of the eXtended Finite Element Method (XFEM), the trial functions take the following form:
(8) 
where is the set of all the nodes in the system, is the set of nodes which are completely cut by the crack, is the set of nodes which contain the crack tips as shown in Figure 2. are the standard shape functions associated with the standard DOF , is the Heaviside function associated with the enriched DOF, and are the tip enrichment functions associated with the DOF, that span the near tip asymptotic fields:
(9) 
Following the Galerkin procedure, this modification to the trial function space leads to an enlarged problem to solve:
(10) 
where
(11) 
where, the superscript refers to standard FEM components, refers to the Heaviside enrichment terms, refers to the asymptotic enrichment terms and other terms can be defined similarly. The augmentation of nonpolynomial functions to the trial function space, makes the numerical integration nontrivial. This has been of particular interest among research community, for example, equivalent polynomial approach by Ventura ventura2006elimination () and Ventura et al., ventura2009fast (), conformal mapping natarajan2010integrating (), Duffy transformation mousavisukumar2010 (), generalized Gaussian quadrature mousavisukumar2010a (), strain smoothing technique bordas2011performance (), exponentially convergent mapping minnebo2012 (),polar mapping parkpereira2009 () and very recently by using Euler’s homogeneous function theorem and Stoke’s theorem chinlasserre2016 (). In bordas2011performance (), the strain smoothing technique was combined with the XFEM, coined as smoothed XFEM (SmXFEM) to integrate over elements intersected with discontinuous surface. The main advantages of the SmXFEM are that no subdivision of the split elements is required and that the derivatives of the shape functions (including the enrichment functions) are not required. However, it was observed that the error level was greater when compared to the conventional XFEM, whilst yielding similar convergence rates.
3 Linear smoothing in the XFEM
The strain smoothing was introduced in chen2001stabilized () for the meshfree methods, which was later extended to the FEM by Liu and coworkers liu2007smoothed (). The basic idea is to compute a strain field, referred to as ‘smoothed’ strain field by evaluating the weighted average of the standard (or compatible) strain field. The support domain of the associated material point can be chosen based on surrounding cells, nodes or edges. In this paper, we restrict our discussion only to the cell based strain smoothed FEM. Within this framework, at a point in element the smoothed strain is given below
(12) 
In terms of the standard element shape function derivatives, , the smoothed derivatives are given by
(13) 
where, is the smoothing function and . By invoking the Divergence theorem Eq. 2 can be written as
(14) 
This form of the strain has the following advantages

Domain integration is reduced to a boundary integration along the smoothing cells

Does not require the derivatives of the shape functions and hence does not need the Jacobian

Does not need isoparametric mapping there by giving a leverage on the distortion level of the mesh
The choice of the smoothing function and the integration order used, decide the accuracy of the smoothed strain field. If a constant smoothing function is used, the method is termed the SFEM. It was shown in bordas2011performance (); natarajan2014convergence () that the gradient becomes inaccurate for nonpolynomials and higher order polynomial functions. Same issue was also faced with in the context of meshfree approximations and in duan2012second () this inaccuracy was addressed by introducing an additional domain integral term which ensures consistency between the shape functions and their derivatives. This modified equation was termed the Divergence Consistency (DC). It was also shown that such consistency requirement is implicitly satisfied, if linear field is used. It can be seen that Equation(15) would reduce to Equation(14), if is a constant.
(15) 
where, is the contour of the smoothing cell. Here the domain integral term vanishes as the smoothing function is constant over the domain. Since we assumed linear displacement functions, the strain would be a constant and a unique value can be computed using a single equation. Hence requiring just one interior Gauß points. This can be written as
(16a)  
(16b)  
(16c)  
(16d) 
where, is the number of subcells in an element, is the volume of the subcell, is the number of boundary surfaces of the subcell, and are the area and Gauß point of the boundary surface . The smoothing technique has been very efficient for polyhedral elements since the polyhedrons can be divided into number of subcells (usually tetrahedrons) and the stiffness matrix is summed up over each subcell. It can be seen in Equation(16) that the derivatives of the shape functions are not needed in order to evaluate the strains. Hence, the computation of Jacobian is avoided. This also avoids the associated isoparametric mapping. The stiffness matrix is evaluated as in the regular finite element method by replacing the terms in the strain gradient matrix with the terms evaluated above and summing it up over the subcells. The constant smoothing technique when applied to elements other than Constant Strain elements (3noded triangles and 4noded tetrahedrons) yields results which are bounded by the results of reduced integration procedure (smoothing over one subcell) and full integration procedure (smoothing over infinite number of subcells). The method is hence not âvariationally consistentâ for any number of subcells other than 1 and bordas2010approximation (), whereas the linear smoothing procedure is variationally consistent. The constant smoothing and linear smoothing schemes differ in the choice of the smoothing function. In the linear smoothing scheme the basis function used is in case of hexahedral subcells and if tetrahedral subcells are used. Figure 3 shows one possible division of hexahedral elements into tetrahedral elements (also referred as subcells in the literature) for the purpose of numerical integration. The number of terms in the basis function should be consistent with the number of Gauß points to obtain a unique solution. Since a linear basis function is being used the domain integral term which results as a consequence of the DC does not vanish and hence it has to be computed by using the appropriate order of Gaussian integration. In the case of tetrahedral subcells the system of equations for a linear basis would be
(17a)  
(17b)  
(17c)  
(17d) 
for
(18a)  
(18b)  
(18c)  
(18d) 
for .
(19a)  
(19b)  
(19c)  
(19d) 
for . Here represents the shape function associated with the node of the parent element. It is independent of the subcell. The location of the gauss points for the boundary integration and domain integration in a tetrahedral subcell are shown in Figure 4. Let the natural coordinates of the interior gauss point of a subcell be and its associated gauss weight be ; coordinates of the boundary of the subcell be and the associated weights be . The unit outward normal associated with the gauss point of the boundary of the subcell is denoted by . The smoothed derivatives are computed numerically as follows
(20) 
(21) 
(22) 
(23) 
(24) 
The smoothed derivative of the shape function evaluated at the four interior Gauß points of a tetrahedral subcell is given by
(25)  
The same procedure is to be repeated for all the shape functions of the parent element. For the interior gauss point of a subcell of a sided polygon the smoothed strain displacement matrix is given by
where,  (26) 
(27) 
For the displacement approximation given by Equation (8), the compatible strain field is given by:
(28) 
where is the vector of degrees of freedom, , and contains the strain displacement terms corresponding to the regular finite element part, Heaviside enriched part and the tip enriched part. The components of the compatible strain field are:
(29) 
where,
The smoothed counterpart of the above compatible strain field can be obtained by following the procedure outlined earlier. The elements that are intersected by the discontinuous surface is divided into tetrahedra and a linear smoothing basis, is chosen to evaluate the smoothed strain.
Remark 1
In case of two dimensions, the subcell is a triangle and the smoothing procedure can be derived from the linear basis
(30) 
with derivative
(31) 
4 Numerical Examples
In this section, the accuracy and the convergence properties of the proposed formulation is numerically studied within the framework of linear elastic fracture mechanics (LEFM) in both two and three dimensions. The domain is discretized with four noded quadrilateral and eight noded hexahedral elements in two and three dimensions, respectively. The numerical results from the present formulation is compared with the conventional XFEM and the SmXFEM bordas2011performance (). The following convention is adopted to compute the stiffness matrix within the framework of the smoothing technique:
Type of element  SmXFEM  LSmXFEM  

two dimensions  Standard elements  4  1 
Tip enriched elements  5  5  
Split enriched elements  8  8  
three dimensions  Standard elements  6  1 
Tip enriched elements  24  24  
Split enriched elements  12  12 
For the conventional XFEM, the elements that are intersected by the discontinuous surface is triangulated and a higher order triangular quadrature scheme is adopted. And for the standard elements, 22 Gauß quadrature rule is adopted. To estimate the error and to study the convergence properties, the norm and the seminorm is used.
4.1 Infinite plate with center crack under farfield tension
Consider an infinite plate with a centre crack subjected to far field tension under plane strain condition has been considered. Consider an infinite plate as shown in Figure 5 A small section ABCD has been solved. The effect of the farfield stress has been accounted by prescribing equivalent displacements as given by following closed form solution Equation (32) in polar coordinates centered at the crack tip.
(32) 
where , the mode I stress intensity factor, is the Poisson’s ratio, is the Youngâs modulus. The dimension has been chosen as 10 x 10 mm. âaâ is chosen as 100 mm. The convergence of the relative error in the displacement ( norm) and the stress intensity factor is shown in Figure 6. It can be seen that in general all the three methods yields a rate of convergence of 1 in the norm and 0.5 in the seminorm. For a given dof, the conventional XFEM yields slightly accurate results than the SmXFEM or the LSmXFEM but the errors are within the same order. Moreover it is noted that in the XFEM, 13 integration points per triangle (for the tip element) is used when compared to three integration points in case of LSmXFEM and one integration point in case of SmXFEM. The suboptimal rate of convergence is due to the fact that we are employing topological enrichment scheme as opposed to geometric enrichment.
4.2 Finite Plate with an edge crack subject to tensile and shear stresses
Next, consider a finite element with an edge crack subjected to tensile and shear stresses as shown in Figure 7. A consistent system of units is used for the analysis.
Plate subjected to tensile stress
In this case, the dimension of the plate is 1 x 2 units. The Youngs’ modulus, E and Poissonâs ratio, are taken as 1000 and 0.3 respectively. A state of plane strain condition is assumed. The crack width is taken as 0.5 units. The obtained SIF are compared with the reference empirical solutionewalds1984fracture ():
(33) 
where, , is the crack width ratio, is the halfcrack width and is the plate width. The convergence of the relative error in the stress intensity factor is shown in Figure 8. It can be seen that the all the three methods converge at almost identical rates ( 0.5). The results of LS scheme are better than the CS scheme and are almost equal with the conventional XFEM.
Plate subjected to shear stresses
In this case, the dimensions of the plate are taken as 7 units and 16 units. The plate is subjected to shear stress on the top edges, while the displacements are constrained at the bottom edge. The crack width is taken as 3.5 units. The Youngs’ modulus, and Poissonâs ratio, are taken as 3 10 and 0.25 respectively. Plane strain condition is assumed. The reference SIF is taken from ewalds1984fracture (), which is = 34 units, = 4.55 units. The convergence of the and the are presented in Figure 9. It is again seen that all the three methods have similar convergence rates. The LS scheme is also more accurate than the CS scheme with a very minor additional computational expense. It is again recalled that the additional integration points still require only the shape function values which can be obtained by linear interpolation along the boundary. The error can be attributed to the inadequate approximation space in the local crack tip region, i.e, the asymptotic fields are approximated by a linear field.
4.3 Plate with an inclined center crack subject to tensile stresses
Next, to illustrate the efficacy of the formulation SIFs in case of mixedmode loading conditions, consider an inclined center crack subjected to far field tension (see Figure 10). The dimensions of the plate are taken as 20 20. The crack width, is chosen as 2 units. A far field uniform tensile stress, 110 units is applied with Young’s modulus, 110 and Poisson’s ratio, 0.3. The accuracy of the numerically computed SIFs are compared with analytical SIFs given by:
(34) 
where is the inclination of the crack measured anticlockwise from the axis. Based on a progressive refinement, it was observed that a structured mesh of 100 100 quadrilateral mesh is adequate. The influence of the crack angle and different modelling approaches, viz., XFEM, SmXFEM, LSmXFEM on the SIFs are shown in Figure 11. It can be seen that the results from the proposed approach are accurate and comparable with the conventional XFEM and slightly more accurate than the SmXFEM.
4.4 Finite plate with a throughthickness edgecrack subject to tensile stresses
As a last example, the linear smoothing technique is extended to three dimensional domain with a throughthethickness edge crack subjected to uniform tensile stress as shown in Figure 12 with dimensions 1 and 3. The displacement at the bottom face is constrained in all directions and a uniform tensile stress 110 is applied on the top face. The material properties are: Young’s modulus 110 and Poisson’s ratio 0.3. The domain is discretized with structured eight noded hexahedral elements and the normalised SIF from saputra2015computation () is taken as the reference solution.
The smoothed strain field over a standard element is computed without any further subdivisions and with as a smoothing function. For the elements that are intersected by the discontinuous surface, the element is subdivided into tetrahedra and is chosen as the smoothing function. In case of the LSmXFEM a total of 96 Gauß points are used in case of tip enriched elements where as 300 Gauß points are used in the conventional XFEM. In the case of Heaviside enriched elements 48 Gauß points are used in case of LSmXFEM and 60 Gauß points are used in case of the conventional XFEM. The convergence of the relative error in the normalised stress intensity factor is shown in Figure 13. It can be seen that the LSmXFEM is more accurate than the SmXFEM and is in good agreement with the conventional XFEM.
5 Conclusions
In this paper the Linear smoothing (Second order smoothing) was discussed and a method to couple it with the extended finite element method was presented. The developed method was used to solve problems with discontinuities and singularities in both two and three dimensions. The method also involves a rational integration procedure employing the Greenâs theorem. The performance of the linear smoothing scheme for enriched approximation space was studied. Through numerical examples it was shown that the Linear smoothing scheme is more accurate than its constant counterpart. The linear smoothing scheme leads to almost identical results to the standard extended finite element method, but it requires fewer quadrature points, viz., approximately onefourth to what is required with the conventional XFEM.
The constant smoothing and the linear smoothing technique is extended to three dimensions for the first time. Although the presented example in three dimensions is for straight crack, it can be easily extended to other crack profiles. The superior accuracy of the linear smoothing technique is also obtained in the three dimensional case. These results are attributed to the superior approximation properties of the linear smoothing compared to the constant strain smoothing, which is immediately apparent for problems involving complex, nonpolynomial, enrichment functions. The remaining, incompressible, error level is attributed to the inadequate approximation space in the smoothed strain, i.e. to the inability of a linear smoothed strain to approximate the singular strains provided by the enriched approximations. Future, ongoing work includes the enrichment of the smoothing space with suitable enrichment functions in order to investigate any additional accuracy improvements as well as the introduction of the approach in recently developed stable extended finite element schemes agathoschatzi2016 (); agathoschatzi2016a ().
References
 [1] Bordas SP, Rabczuk T, Hung NX, Nguyen VP, Natarajan S, Bog T, Hiep NV, et al.. Strain smoothing in fem and xfem. Computers & structures 2010; 88(23):1419–1443.
 [2] Bordas S, Natarajan S, Kerfriden P, Augarde CE, Mahapatra DR, Rabczuk T, Pont SD. On the performance of strain smoothing for quadratic and enriched finite element approximations (xfem/gfem/pufem). International Journal for Numerical Methods in Engineering 2011; 86(45):637–666.
 [3] Francis A, OrtizBernardin A, Bordas S, Natarajan S. Linear smoothed polygonal and polyhedral finite elements 2015; .
 [4] Oliyer J. Continuum modelling of strong discontinuities in solid mechanics using damage models. Computational Mechanics 1995; 17(12):49–61.
 [5] Belytschko T, Black T. Elastic crack growth in finite elements with minimal remeshing. International journal for numerical methods in engineering 1999; 45(5):601–620.
 [6] Dolbow J, Belytschko T. A finite element method for crack growth without remeshing. International journal for numerical methods in engineering 1999; 46(1):131–150.
 [7] Melenk JM, Babuška I. The partition of unity finite element method: basic theory and applications. Computer methods in applied mechanics and engineering 1996; 139(1):289–314.
 [8] Bordas S, Moran B. Enriched finite elements and level sets for damage tolerance assessment of complex structures. Engineering Fracture Mechanics 2006; 73:1176–1201.
 [9] Bordas S, Conley JG, Moran B. A simulationbased design paradigm for complex cast components. Engineering with Computers 2007; 23(1):25–37.
 [10] Wyart E, Coulon D, Duflot M, Pardoen T, Remacle JF, Lani F. A substructured FE shell/XFE 3D method for crack analysis. International Journal for Numerical Methods in Engineering 2007; 72(7):757–779, doi:10.1002/nme.2029.
 [11] Duarte CA, Hamzeh ON, Liszka TJ, Tworzydlo WW. A generalized finite element method for the simulation of threedimensional dynamic crack propagation. Computer Methods in Applied Mechanical and Engineering 2001; 190:2227–2262.
 [12] Babuška I, Banerjee U. Stable generalized finite element method (SGFEM). Computer Methods in Applied Mechanical and Engineering 2012; 201–204:91–11.
 [13] Rabczuk T, Bordas S, Zi G. A threedimensional meshfree method for continuous crack initiation, nucleation and propagation in statics and dynamics. Computational Mechanics August 2007; 40(3):473–495, doi:10.1007/s0046600601221.
 [14] Belytschko T, Guo Y, Liu WK, Xiao SP. A unified stability analysis of meshless particle methods. International Journal for Numerical Methods in Engineering 2000; 48:1359–1400.
 [15] Chen JS, Wu CT, Yoon S, You Y. A stabilized conforming nodal integration for galerkin meshfree methods. International journal for numerical methods in engineering 2001; 50(2):435–466.
 [16] Liu G, Dai K, Nguyen T. A smoothed finite element method for mechanics problems. Computational Mechanics 2007; 39(6):859–877.
 [17] Bordas S, Natarajan S. On the approximation in the smoothed finite element method (sfem). International Journal for Numerical Methods in Engineering 2010; 81(5):660–670.
 [18] Liu G, NguyenThoi T, Lam K. An edgebased smoothed finite element method (esfem) for static, free and forced vibration analyses of solids. Journal of Sound and Vibration 2009; 320(4):1100–1130.
 [19] Liu G, NguyenThoi T, NguyenXuan H, Lam K. A nodebased smoothed finite element method (nsfem) for upper bound solutions to solid mechanics problems. Computers & structures 2009; 87(1):14–26.
 [20] Nguyen NT, Rabczuk T, NguyenXuan H, Bordas S. A smoothed finite element method for shell analysis. Computer Methods in Applied Mechanics and Engineering 2008; 198:165–177.
 [21] NguyenXuan H, Bordas S, NguyenDang H. Smooth finite element methods: convergence, accuracy and properties. International Journal for Numerical Methods in Engineering 2008; 74(2):175–208.
 [22] Dai K, Liu G, Nguyen T. An nsided polygonal smoothed finite element method (nsfem) for solid mechanics. Finite elements in analysis and design 2007; 43(11):847–860.
 [23] Babuška I, Melenk JM. The partition of unity method. International journal for numerical methods in engineering 1997; 40(4):727–758.
 [24] Natarajan S, Ooi ET, Chiong I, Song C. Convergence and accuracy of displacement based finite element formulations over arbitrary polygons: Laplace interpolants, strain smoothing and scaled boundary polygon formulation. Finite Elements in Analysis and Design 2014; 85:101–122.
 [25] Duan Q, Li X, Zhang H, Belytschko T. Secondorder accurate derivatives and integration schemes for meshfree methods. International Journal for Numerical Methods in Engineering 2012; 92(4):399–424.
 [26] Ventura G. On the elimination of quadrature subcells for discontinuous functions in the extended finiteelement method. International Journal for Numerical Methods in Engineering 2006; 66(5):761–795.
 [27] Ventura G, Gracie R, Belytschko T. Fast integration and weight function blending in the extended finite element method. International journal for numerical methods in engineering 2009; 77(1):1–29.
 [28] Natarajan S, Mahapatra DR, Bordas S. Integrating strong and weak discontinuities without integration subcells and example applications in an xfem/gfem framework. International Journal for Numerical Methods in Engineering 2010; 83(3):269–294.
 [29] Mousavi S, Sukumar N. Generalized duffy transformation for integrating vertex singularities. Computational Mechanics 2010; 45:127–140.
 [30] Mousavi S, Sukumar N. Generalized gaussian quadrature rules for discontinuities and crack singularities in the extended finite element method. Computer Methods in Applied Mechanical and Engineering 2010; 199:3237–3249.
 [31] Minnebo H. Threedimensional integration strategies of singular functions introduced by the XFEM in the LEFM. International Journal for Numerical Methods in Engineering 2012; 92:1117–1138.
 [32] Park K, Pereira J, Duarte C, Paulino G. Integration of singular enrichment functions in the generalized/extended finite element method for threedimensional problems. International Journal for Numerical Methods in Engineering 2009; 78:1220–1257.
 [33] Chin E, Lasserre J, Sukumar N. Modeling crack discontinuities without elementpartitioning in the extended finite element method. International Journal for Numerical Methods in Engineering 2016; doi:10.1002/nme.5436.
 [34] Fries TP, Belytschko T. The intrinsic xfem: a method for arbitrary discontinuities without additional unknowns. International journal for numerical methods in engineering 2006; 68(13):1358–1385.
 [35] Belytschko T, Moës N, Usui S, Parimi C. Arbitrary discontinuities in finite elements. International Journal for Numerical Methods in Engineering 2001; 50(4):993–1013.
 [36] Fries TP. A corrected xfem approximation without problems in blending elements. International Journal for Numerical Methods in Engineering 2008; 75(5):503–532.
 [37] Mohammadi S. Extended finite element method: for fracture analysis of structures. John Wiley & Sons, 2008.
 [38] Pommier S, Gravouil A, Moes N, Combescure A. Extended finite element method for crack propagation. John Wiley & Sons, 2013.
 [39] Béchet É, Minnebo H, Moës N, Burgardt B. Improved implementation and robustness study of the xfem for stress analysis around cracks. International Journal for Numerical Methods in Engineering 2005; 64(8):1033–1056.
 [40] Laborde P, Pommier J, Renard Y, Salaün M. Highorder extended finite element method for cracked domains. International Journal for Numerical Methods in Engineering 2005; 64(3):354–381.
 [41] Fish J. Finite element method for localization analysis. PhD Thesis, Northwestern University 1989.
 [42] Sukumar N, Moës N, Moran B, Belytschko T. Extended finite element method for threedimensional crack modelling. International Journal for Numerical Methods in Engineering 2000; 48(11):1549–1570.
 [43] Natarajan S, Bordas S, Roy Mahapatra D. Numerical integration over arbitrary polygonal domains based on schwarz–christoffel conformal mapping. International Journal for Numerical Methods in Engineering 2009; 80(1):103–134.
 [44] Sukumar N, Tabarraei A. Conforming polygonal finite elements. International Journal for Numerical Methods in Engineering 2004; 61(12):2045–2066.
 [45] Bordas S, Menk A. XFEM  The Extended Finite Element Method. Wiley, 2016. URL https://books.google.co.in/books?id=2ygSywAACAAJ.
 [46] Ewalds H, Wanhill R. Fracture mechanics. Edward Arnold, 1984. URL https://books.google.co.in/books?id=4rNRAAAAMAAJ.
 [47] Saputra AA, Birk C, Song C. Computation of threedimensional fracture parameters at interface cracks and notches by the scaled boundary finite element method. Engineering Fracture Mechanics 2015; 148:213–242.
 [48] Agathos K, Chatzi E, Bordas SP, Talaslidis D. A well conditioned and optimally convergent XFEM for 3D linear elastic fracture. International Journal for Numerical Methods in Engineering 2016; 105(9):643–677.
 [49] Agathos K, Chatzi E, Bordas SP. Stable 3D extended finite elements with higher order enrichment for accurate non planar fracture. Computer Methods in Applied Mechanical and Engineering 2016; 306:19–46.