## 1 Introduction  
In the late 1980s, a new material concept called functionally graded material (FGM) was proposed to resolve the inherent problem of traditional lamination-type composites [1]. The sharp material discontinuity across the layer interface causes the stress concentration, which may trigger the layer delamination. This crucial stress concentration can be significantly minimized by inserting a graded layer between two distinct homogeneous material layers [2; 3]. This is because the material discontinuity completely disappears according to the material composition gradation, where the constituent particles of two base materials are mixed up in a random microstructure within a graded layer to maximize the desired performance [4; 5; 6]. Naturally, a functionally graded material is an inhomogeneous material, with spatially non-uniform material properties characterized by continuity and functionality. In addition to the suppression of stress concentration, the material concept of FGM rapidly and continuously spread throughout engineering and fields [7; 8; 9; 10]. Early research efforts were concentrated on material characterization, fabrication, modeling, and analysis [1; 11; 12]. This was because the mechanical behaviors of FGMs are governed by the geometric dimensions and orientation, microstructure, and volume fractions of constituent particles. Recently, however, the concern toward the crack problems has increased because the structural failure of FGMs is dominated by micro-cracking [7; 13; 14]. In this regard, an accurate numerical prediction of stress intensity factors and the crack propagation has been an essential research subject [15; 16]. For these subjects, an accurate reproduction of the \(\frac{1}{\sqrt{r}}\) singularity in the near-tip stress field in highly heterogeneous media becomes a key task [17; 18; 19; 20]. To evaluate the stress intensity factors of FGMs with cracks, one can consider the use of the well-known \(J-\) or \(M-\)integral methods. However, these conventional indirect integral methods cannot reflect the spatially varying material properties of FGMs. The studies on the fracture mechanics of inhomogeneous bodies were initiated in the 1970-80s by assuming the spatially varying elastic modulus as an exponential function [17; 21]. The standard integral methods were modified and/or refined to reflect the spatial non-uniformity of material properties by subsequent investigators. The most works were made by utilizing the finite element method [15; 22; 23; 24; 25]. However, since late 1990s, the employment of meshless methods to crack problems of inhomogeneous bodies has been actively advanced, particularly for the computation of SIFs by the modified \(M-\)integral method [26; 27]. Here, the extension of the natural element method (NEM) is worth noting, even though it was restricted to 2-D homogenous material [28]. The natural element method was introduced to overcome the demerits of conventional meshless methods [29], the difficulty in enforcing the displacement constraint and the numerical integration. The Laplace interpolation functions in NEM strictly obey the Kronecker delta property so that the imposition of displacement constraint becomes easy. In addition, Delaunay triangles, which were produced during the process for introducing Laplace interpolation functions, also serve as a background mesh for the numerical integration. In particular, PG-NEM can further improve the numerical integration accuracy by maintaining the consistency between the Delaunay triangle and the integration region [30]. Although Laplace interpolation functions provide the high smoothness of \(C^{1}\)-continuity, there is still room for further improvement in capturing the high stress singularity at the crack tip. In this context, this paper introduces an enriched PG-NEM to explore whether and how much the enrichment of interpolation function increases the prediction reliability of stress intensity factors for FGMs.