## 3 Enriched Petrov-Galerkin Natural Element Method
### Numerical Implementation of \(\widetilde{M}^{(1,2)}\)  
Figure 4 schematically represents the identification of domain of integration \(A\) and the definition of weighting function \(q(\mathbf{x})\), which are prerequisite for the numerical implementation of \(\widetilde{M}^{(1,2)}\) in Eq. (6). For these two things, a circle of the radius \(r_{int}\), which is originated at the crack tip, is first imaginarily drawn on the NEM grid. Then, the nodal values of \(q(\mathbf{x})\) are assigned based on this circle, such that unity is specified to all the interior nodes while zero to the outer remaining nodes. Then, referring to Figure 2, the boundary of a square composed of eight darkened Delaunay triangles serves as an internal path \(\Gamma\). On the other hand, the grayed region outside the internal path \(\Gamma\) automatically becomes a domain of integration \(A\). According to the properties of Laplace interpolation functions [29], the defined weighting function \(q(\mathbf{x})\) becomes unity within a set of darkened Delaunay triangles while it has the values between zero and unity within the grayed domain of integration \(A\).  
Both the weighting function \(q(\mathbf{x})\) and the derivative \(\partial q/\partial x_{j}\) become zero outside the domain of integration \(A\). Furthermore, the derivative \(\partial q/\partial x_{j}\) vanishes within the darkened squares. So, only the grayed Delaunay triangles are taken for the modified interaction integral \(\widetilde{M}^{(1,2)}\) such that  
\[\widetilde{M}^{(1,2)}=\sum_{K=1}^{M_{G}}\widetilde{M}^{(1,2)}_{K} \tag{40}\]  
Here, \(M_{G}\) indicates the total number of grayed Delaunay triangles, while \(\widetilde{M}^{(1,2)}_{K}\) stand for the triangle-wise modified interaction integrals, which are defined by  
\[\begin{array}{rl}\widetilde{M}^{(1,2)}_{K}=&\sum_{\ell=1}^{INT}\left(\left[ (\sigma^{(1)}_{ij})\frac{\partial u^{(2)}_{ij}}{\partial x_{1}}+\sigma^{(2)}_ {ij}\frac{\partial u^{(1)}_{i}}{\partial x_{1}}-W^{(1,2)}\delta_{1j}\right] \frac{\partial q}{\partial x_{\ell}}\right|_{x_{\ell}}w_{\ell}|\!|\!|_{x_{\ell}} \\ &\left.+\frac{1}{2}\right[\sigma^{(1)}_{ij}\frac{\partial c^{(2)}_{ij}}{ \partial x_{1}}-\frac{\partial c^{(2)}_{ij}}{\partial x_{1}}\varepsilon^{(1)}_ {ij}+\sigma^{(2)}_{ij}\frac{\partial c^{(1)}_{ij}}{\partial x_{1}}-\frac{ \partial c^{(1)}_{ij}}{\partial x_{1}}\varepsilon^{(2)}_{ij}\right]q(x_{\ell })w_{\ell}|\!|\!|_{x_{\ell}}\end{array} \tag{41}\]  
with \(INT\), \(\mathbf{x}_{\ell}\), and \(w_{\ell}\) being the total number, co-ordinates, and weights of sampling points for Gauss numerical integration for triangles. Meanwhile, \(T_{K}\) in Figure 4 indicates the geometry transformation between the master element \(\Omega\) and the grayed triangle \(\Omega_{K}\) and the master element \(\Omega\), which is defined by  
\[T_{K}:\ x_{\ell}=\sum_{i=1}^{3}x_{i}\psi_{i}(\xi,\eta)_{\ell},\ \ \ y_{\ell}=\sum_{i=1}^{3}y_{i}\psi_{i}(\xi,\eta)_{\ell} \tag{42}\]  
where \(\psi_{i}\) stand for the constant-strain FE basis functions, which are employed for the test displacement field \(\mathbf{v}(\mathbf{x})\).  
Figure 4: A crack tip-originated circle and a graded domain of integration \(A\).